Optimal transport and colossal ionic mechano-conductance in graphene crown ethers

See allHide authors and affiliations

Science Advances  12 Jul 2019:
Vol. 5, no. 7, eaaw5478
DOI: 10.1126/sciadv.aaw5478


Biological ion channels balance electrostatic and dehydration effects to yield large ion selectivity alongside high transport rates. These macromolecular systems are often interrogated through point mutations of their pore domain, limiting the scope of mechanistic studies. In contrast, we demonstrate that graphene crown ether pores afford a simple platform to directly investigate optimal ion transport conditions, i.e., maximum current densities and selectivity. Crown ethers are known for selective ion adsorption. When embedded in graphene, however, transport rates lie below the drift-diffusion limit. We show that small pore strains (1%) give rise to a colossal (100%) change in conductance. This process is electromechanically tunable, with optimal transport in a primarily diffusive regime, tending toward barrierless transport, as opposed to a knock-on mechanism. These observations suggest a novel setup for nanofluidic devices while giving insight into the physical foundation of evolutionarily optimized ion transport in biological pores.


The search for universal transport mechanisms is an underlying theme in ion channel research (1). Nonetheless, even foundational questions, such as the primary mechanism for selectivity of K+ over Na+ in the potassium ion channel family, remain contentious: Is the smaller dehydration energy of the K+ ion responsible or is it the “snuggle fit” of K+ in the selectivity filter (24)? Solving these puzzles requires a detailed means to control and characterize ion transport. Biological ion channels are notoriously difficult to manipulate, requiring mutagenesis and near-native conditions for patch-clamp studies, thus limiting experimental characterization and tunability. These large macromolecules are also theoretically complex, requiring hefty resources to conduct statistically meaningful simulations. In contrast, biomimetic synthetic pores—or, as addressed here, pores that retain core aspects of some biological channels—may provide a simple platform to explore ion transport mechanisms under a broad range of conditions (5).

In this work, we demonstrate that graphene crown ether pores (6) have competing electrostatic and dehydration contributions to transport, reminiscent of biological ion channel mechanisms and selectivity. The discovery of crown ethers half a century ago triggered a cascade of research that evolved into the fields of host-guest and supramolecular chemistry (7). The isolation of graphene (8) has likewise inspired diverse fields of study, including those that leverage atomically thin pores for biosensing and nanoscale separation (9). Graphene crown ether pores lie at the crossroads, affording a means to address transport and selectivity in a synthetic analog of biological ion channels. Using a graphene-embedded 18-crown-6 ether pore as a model system, we examine the mechanism for ion translocation under variable strain and applied voltage (see Fig. 1). Our results indicate that a few percent changes in pore size can induce a several-fold increase in ionic current, underscoring the subnanometer fine tuning between electrostatics and dehydration seen in biological contexts. We similarly find that small changes in pore radius and charge can optimize transport and even alter the underlying mechanism of ion translocation. These results collectively suggest that graphene crown ether pores are an effective mechanistic probe for understanding transport in both biological and artificial ion channels alike.

Fig. 1 Potential transport mechanisms in graphene crown ether pores.

(A) Using partial charges qO = −0.24e, consistent with the electrostatic potential from DFT, the ion transport mechanism is drift-diffusion. In this case, a K+ (purple sphere) finds an empty pore and translocates through it; the pore then remains empty again for several nanoseconds (see the Supplementary Materials). (B) At larger partial charges (qO = −0.54e or qO = −1.0e), the vicinity of the charge separation at the pore rim results in an energetic well deep enough to trap a K+ (light purple). For a current to be present, an incoming K+ (dark purple) knocks out and replaces the trapped K+. This yields a two-step knock-on mechanism reminiscent of some biological potassium ion channels. However, for qO = −0.54e, the mechanism shifts toward a majority drift-diffusion process at moderate (1%) strain due to a shallowing of the free-energy well. Oxygen, carbon, and hydrogen atoms are small red, gray, and white spheres, respectively.


Partial charge distribution

Unlike an isolated crown ether, an 18-crown-6 pore in graphene is believed to assume a planar conformation (6). However, while critical to molecular dynamics (MD) simulations, the charge distribution and response of this pore is unknown. Previous simulations on this and related systems assume oxygen partial charges that span between qO = −0.21e and qO = −0.74e (1015), over which transport properties will vary markedly. To provide a more comprehensive approach, we reexamine this assignment using density functional theory (DFT) (see Methods and the Supplementary Materials for details). Determinations from electrostatic potential fitting without an ion present—consistent with additive force fields and the electrostatic environment predicted by DFT—yield a uniform value of qO = −0.24e, decreasing to −0.23e when a K+ is present. Chemically distinct, but oxygen containing, fragments in optimized potentials for liquid simulations (OPLS) and Chemistry at Harvard Macromolecular Mechanics (CHARMM) force fields give substantially higher partial charges, ≥ −0.4e, of which we take qO = −0.54e as a representative example. An upper bound of qO = −1.0e is given by Bader analysis, which tends toward larger values of partial charge and represents the maximum expected from any electronic structure calculation. Physically reasonable values of qO should span the low end of the range −0.2e to −0.6e, as evidenced by electrostatic fits of qO = −0.34e for the 18-crown-6 alkyl ether and qO = −0.54e for a highly dipolar glycine dipeptide backbone carbonyl. The graphene crown ether should be lower than both. We perform MD simulations using three distinct assignments qO = (−0.24, −0.54, −1.0)e to capture the full range of potential behavior, showing both universal features of transport and aspects particular to different assignments. We test the effect of strain on ion transport using different applied biases and deformations, well within the elastic limit of graphene (16, 17).

Ion transport mechanism

Since this pore is tiny and negatively charged, only K+ contributes to the current in a KCl solution (set at a concentration of 1 M in our case). The mechanism of transport depends on the ion’s binding strength in the pore, largely regulated by qO (as an indicator of the oxygen-carbon charge separation) and atomic structure. The mechanism is normal drift-diffusion for qO = −0.24e, shifting to a knock-on type for −0.54e and −1.0e (see the illustration in Fig. 1). For drift-diffusion, the residence time (on the order of 0.1 ns at 0.25 V) is much shorter than the delay to the next ion-crossing event (on the order of 1 ns at 0.25 V). Conversely, for the two-step knock-on mechanism seen at higher charge, the residence time (on the order of 1 ns) is much larger than the delay time (less than 10 ps for most events; see the Supplementary Materials). However, even for qO = −0.54e, drift-diffusion starts to dominate at high strain and/or voltage, both of which weaken ion binding to the pore.

Ionic mechano-conductance

Figure 2A depicts a potential experimental setup, in which a graphene sheet is selectively strained by deforming an underlying substrate. Our simulations mimic these conditions by increasing the cross-sectional area of the membrane to induce a strain. The relative change in the current as a function of strain, for different values of qO and applied voltages, is shown in Fig. 2B. At low voltages (0.25 and 0.5 V), the ionic current increases rapidly with strain, orders of magnitude larger than expected from the increase in pore area alone. Instead, this colossal increase is indicative of a change in the energetics and dynamics of ion transport. In the large voltage regime (1 V), the current remains constant with increasing strain when qO = −0.54e. This is not altogether unexpected, as a large voltage drop across the membrane suppresses the influence of the free-energy profile on translocating ions. The dependence of the current on qO and V is in Fig. 2C. In this case, the current increases with qO between qO = −0.24e and qO = −0.54e while giving a turnover to decreasing current at a higher charge.

Fig. 2 Colossal ionic mechano-conductance.

(A) Schematic of graphene on (or embedded within) a polymeric matrix support (such as ∼50-nm-thick epoxy resin) with a window where the crown ether pore is located. The membrane can be strained by stretching or bending (e.g., via a piezoelectric actuator offset from the window). Alternative experimental setups are possible, such as metallic regions that bind the graphene and are used to apply strain. (B) Relative change in current versus strain in the crown ether pore (for qO = −0.24e and −0.54e) at different voltages. At lower voltages, the current increases substantially with strain, as shown by the fitted dashed lines. (C) Current without strain (I0) for different values of qO and V. Going from qO = −0.24e to qO = −0.54e, the current increases as highly charged oxygen atoms compensate the loss of waters from the hydration layers. However, going from qO = −0.54e to qO = −1.0e, the current decreases because a further increase in the pore charge makes it harder for an ion to escape the potential well in the pore. The typical drift-diffusion limit—set by an access resistance to an uncharged pore with an approximate radius of 0.1 nm—at 0.25 V is about 0.35 nA. Only the intermediate charge case—the one most analogous to biological systems—approaches this limit: It is about a factor of 2 lower at a 2% strain. Further strain allows it to reach this limit. The smaller charge case can approach within a factor of 6. The range of accessible currents is due to the extensive mechanistic leeway permitted by the intermediate charge. The error bars are the standard error (SE) from five parallel runs.

The energy landscape

Our results indicate that the energetics of pore permeation lie in a regime where small electromechanical variations modify current magnitudes and give rise to qualitative changes in behavior. The free energy of a potassium ion in the vicinity of the pore, ΔFK, is approximately the sum of a dehydration energy barrier ΔEdeh and three electrostatic terms (EKO, EKC, and EKK) that capture the interaction of K+ with nearby charged atomsΔFK=ΔEdeh+EKO+EKC+EKK(1)where Edeh = η∑ifiEi is the barrier due to dehydration, fi (Ei) is the fractional dehydration (energy) of the ith hydration layer, and the factor η ≈ 1/2 accounts for the nonlinear effects (18). Each electrostatic component is Eμν = qμnνqν/4πϵ0ϵ(zμν)rμν, where ϵ(zμν) is a position-dependent relative permittivity (19), qμ(ν) is the charge of ion species μ(ν), nν is the number of proximate species v (6 for oxygen, 12 for carbon, and 1 or 0 for a nearby K+), zμν (rμν=zμν2+ρμν2) is the axial (total) distance between species μ and v, and ρμν is the respective radial distance. EKK is the energetic coupling of an incoming K+ to a trapped K+. The van der Waals interactions, the coupling to other solvated ions, and the entropy change should have only marginal contributions when ΔFK varies with strain (e.g., the entropy change into the pore is important, but this varies little with moderate strain). To quantify the variable electrostatic environment reflected in ϵ(z), as well as nonlinear effects in ΔEdeh [stronger polarization of bound water molecules as an ion becomes dehydrated (18)], we use all-atom MD simulations to find the free-energy profile of K+ around the pore (Fig. 3 and Methods). Equation 1 explains the key features in the free-energy profiles for the range of qO under consideration. Together with the MD data, this expression dissects the contributions to the free energy. We will use this equation below to characterize the response of ion transport to strain.

Fig. 3 Free-energy landscape.

Free-energy profile, ΔFK, for a K+ translocating along the z axis of an 18-crown-6 graphene pore at different strains. The charge on the oxygen atom, qO, is presented at the top of each plot. The peaks and valleys in ΔFK(z) are due to the balance between dehydration energy penalties and electrostatic interactions with the charged pore atoms. For qO = −0.54, there is an additional contribution from the K+ occupying the deep potential well, giving a barrier at z ≈ 0.2 nm when the incoming ion attempts to go into the already occupied pore (this is confirmed by a free-energy profile with only 1 K+ and 1 Cl present; see the Supplementary Materials). The illustrations show the position of the K+ in the peaks and valleys. In general, the free-energy profile flattens with increasing—but still small—strain, tending toward barrierless transport and making it easier for the ion to translocate. The error bands are the SE from five parallel runs.

Figure 3 shows that when qO = −0.24e, a free-energy barrier appears outside the pore due to the progressive dehydration of K+. The barrier decreases with increasing proximity to the pore opening as the flanking oxygen atoms begin to compensate for water molecules displaced from the ion’s first hydration layer. This results in a potential well on both sides of the membrane. Further dehydration occurs as the ion enters the pore, giving rise to an overall barrier at z = 0 since the electrostatic attraction between the pore oxygen atoms and the K+ (EKO) is insufficient to compensate for the dehydration energy penalty. This central barrier—and only shallow adjacent wells—ensure that an approaching K+ will almost always find an empty pore, and hence, the EKK term will not contribute in this case. Equation 1 confirms this assignment of the peaks and wells.

This situation is inverted for qO = −0.54e (see Fig. 3), where a deep potential well forms at the center of the pore due to a strong electrostatic attraction with the K+. The central pore acts as a binding site in this case, ensuring that it is almost always occupied. When K+ approaches the pore, the presence of an interstitial K+ repels it, yielding a potential barrier on either side of the pore. A free-energy calculation with single K+ and single Cl confirms this (see the Supplementary Materials). We can estimate ϵr(z) using the free-energy profile from Fig. 3 and the fractional dehydration from MD. For the physical range of pore charge, ϵr(z) is around 4 to 6 within 0.2 nm of the pore, which is consistent with experimental values for the dielectric constant within 1 nm of an interface (20). A similar value for dielectric constant is given in another recent computational study (19).

Barrierless transport

While the transport mechanism is different for unstrained pores, the colossal mechano-conductance change is present at low voltage, regardless of the partial charge. In both partial charge scenarios, the free-energy profile becomes smoother with increasing strain, tending toward “barrierless transport” (19, 21): Pulling charged oxygens away from the pore center reduces their coupling to translocating ions. The colossal change in conductance indicates that there are optimal structural positions for the oxygens: Picometer changes in their position do not change the dehydration contribution to the central barrier (see the Supplementary Materials), but it vastly changes the transport rate, maintaining exclusion of other ions but enhancing transport rates by several fold. Biological ion channels can further make use of the partial charge, engineering not only structural characteristics but also the electronic environment. As we discussed earlier, the larger partial charges we consider—in line with biological channels (22)—show a shift from knock-on to drift-diffusion mechanisms as transport starts to become barrierless. In this manner, the electromechanical environment may exert a large effect on the current.

Mechanical susceptibility

To understand the mechanics underlying this process, we define the susceptibility, χ, of the free energy to a small change in pore sizeχ=dΔFda=μνΔFρμνdρμνdaqKEρ(2)where a is the effective pore radius, 〈μν〉 indicates the pairs KO and KC (we ignore the role of KK interactions, although this is important for knock-on aspects of transport), and dρμν/da ≈ 1. The latter approximation is reasonable since, as the pore size increases, the oxygen and carbons move outward by the same amount. For simplicity, we took this for an ion at the origin (z = 0, ρ = 0), where the susceptibility is just the local radial electric field, Eρ, multiplied by the charge of the translocating ion. Because of the lack of electrostatic screening in the pore [i.e., ϵ(z = 0) is small] and the proximity of the oxygen atoms, this field is very large, resulting in a large χ: Eq. 1 gives χ ≈ 100kBT/nm and χ ≈ 140 kBT/nm for the qO = −0.24e and qO = −0.54e, respectively (using ϵ = 3.9 and 6.4, which we extract via Eq. 1 and the all-atom MD data). These two susceptibilities are closer than their difference in charge would indicate because of the smaller permittivity for qO = −0.24e (the charge qO = −0.54e more strongly attracts counter charges, specifically K+ and hydrogens on water molecules, which help screen the interaction.)

The large susceptibility means that even a 1-pm change in the pore size gives a 0.1- to 0.14kBT change in the free energy. A 1% strain gives about a 7 pm change in radius. Using this in Eq. 1, free energy changes by about 0.7kBT for qO = −0.24e and about 1.0kBT for qO = −0.54e, in line with the all-atom results in Fig. 3. We note that there are other effects occurring: As we will discuss below, the qO = −0.24e case is already in a limit where dehydration and electrostatic effects are comparable, and thus, the change in the free-energy profile versus z reflects both (the barrier reduction at z ≈ 2 nm is due to a changing dehydration contribution, whereas in between those two outer barriers, electrostatic effects dominate and ΔF increases overall). For qO = −0.54e, the potential well shallows, causing a transition from a knock-on to a diffusive mechanism and removing the ion-ion repulsion that causes the external barrier. Nevertheless, in both cases, the electrostatic contribution to the susceptibility captures the increased free energy for an ion in the pore. While a general susceptibility will need to include the totality of effects, electrostatics and dehydration will typically dominate and are key to understanding amplification.

The free-energy profile and its variation (according to χ) determines the current. In particular, the largest free-energy hurdle, ΔF, in the profiles of Fig. 3 are the main factors influencing the conductance. For qO = −0.24e, this will be the jump from about z = 0.1 to 0.2 nm, and for qO = −0.54e, this will be the jump out of the well (z = 0) to the maxima (at about z = 0.2 nm to z = 0.4 nm depending on strain). The current will take the form (see Methods)IAeΔF/kBT(3)where A is a constant that only weakly depends on pore parameters (such as size), but we will take it as fixed. The total amplification of the current then becomesΔIIe(ΔF(s)+ΔF(s=0))/kBT1eα a s χ/kBT1(4)where −ΔF(s) + ΔF(s = 0) ≈ dΔF/da ⋅ δa ≈ αasχ. The parameter α = 1/ada/ds quantifies the relative response of the pore radius to the application of a strain. Both MD and DFT predict that α ≈ 2, meaning that if the graphene has a 1% strain, then the atoms at the pore rim move by about 2%. This is due to a distortion of the rim structure since the pore is an effective impurity and can relieve strain by relaxing angular coordinates. The straightforward relation in Eq. 4 indicates that the gain in current is the exponential of the work performed by the applied strain on a K+ within the pore. For the values in our case, Eq. 4 gives an amplification of 100 to 200% for a 1% strain, in line with the MD results in Fig. 2B.

Optimum transport

In all parameter regimes, the amplified current is due to an overall flattening of the free-energy profile. Ion transport thus tends toward barrierless transport for small values of strain and displays a colossal increase in conductance along the way. However, for qO = −0.24e, the dominant contribution to ΔF shifts from a pathway connecting a minimum (at z = ±0.1 nm) and the external solution to a pathway that hops over the central barrier (the well minima already lie at about ΔF = 0). After this point, ΔF will continually increase due to a decreasing compensation by the oxygen partial charge. There will be a corresponding decrease in current (see Fig. 4), thus indicating an optimum near this value of strain—a regime that may be experimentally accessible.

Fig. 4 Optimum ion transport and selectivity.

Ionic current versus strain in a pore with qO = −0.24e and an applied bias of 0.25 V. For small strain (blue line), the current increases rapidly, commensurate with a decrease in the outer barriers with increasing strain. This gives an overall flattening of the free-energy profile. Transport (and selectivity; see below) is optimal near a 3% strain (green line) when the outer barrier is sufficiently diminished but the central barrier is not too high. Further increases in strain continue to increase the central barrier (dehydration remains unchanged in this strain regime), decreasing the current (red line). Schematics of the free-energy profiles are above each regime. For every value of strain and voltage, the selectivity of K+ over Cl is perfect for all practical purposes, as the negative partial charge of the oxygens on the pore rim create a strong barrier for Cl transport. Separate simulations (each of duration 1 μs with an equal mixture of NaCl and KCl at 1 M Cl) reveal that Na+ does not cross the pore (except for a 0% strain where a single Na+ crossing event occurs during the 1-μs time frame). Because of the near complete exclusion of Na+ in the time frame of the simulation, we can determine a high-confidence lower bound on the selectivity, which goes from 15 at 0% strain to 52 at 3% strain (see the main text for details). In other words, the Na+ and Cl exclusion is maintained, while K+ current changes with strain, and the lower bound on selectivity exactly traces the K+ current in the plot.

In separate calculations with an equal mixture of NaCl and KCl (with Cl at 1 M), we find that only a single Na+ traversed the pore (compare to 47 K+) during 1 μs of simulation for the unstrained pore at a 0.25-V bias. At higher strains, only K+ translocates through the pore within the 1-μs simulation time. The lack of Na+ transport is due to the larger hydration energy of Na+ compared to K+, giving it a higher penalty to go through the pore. The presence of an overall barrier should result in Na+ crossings following Poisson statistics. Because of the rarity of these crossings, we consider two values for the Na+ current: 0.5 and 0.16 pA. The former would give a 5% probability to detect no Na+ crossings during 1 μs and thus is a high confidence upper bound to the Na+ current. The latter gives equal probabilities for 0 and 1 Na+ crossing and thus represents an estimate of the Na+ current, assuming that we are near the threshold to start seeing Na+ crossings, for which the single crossing at 0% strain suggests is the case. Using the upper bound on the Na+ current, the corresponding high-confidence lower bound on the selectivity ranges from about 15 at 0% strain to 52 at 3% strain (note that the K+ currents in the equal mixtures are 7.5 pA at 0% strain and 26.2 pA at 3% strain, which are approximately half the value in the 1 M KCl case shown in Fig. 4 due to the mixture having half the K+ concentration). The estimate of selectivity, however, ranges from 47 at 0% strain to 164 at 3% strain. Very long simulations (much greater than 10 μs) will be necessary to tightly bound the selectivity, but it will be large nonetheless.

For qO = −0.54e, the free energy at z = 0 will not turn from a well to a barrier until about s = 6% (i.e., about δa = 6kBT/χ), with the optimum appearing at potentially larger strains and outside of measurable regimes. Up to this point, the current will increase with strain as the free-energy profile becomes truly barrierless. In either case, the dehydration contribution itself (at z = 0) will not substantially change until the pore becomes sizable, as indicated by the continuing decrease of the current in Fig. 4. The selectivity filters of biological channels have partial charges near qO = −0.54e by virtue of the dipolar peptide backbone (4, 2224), lying within this more strongly coupled regime and affording them extensive mechanistic leeway. These systems generically exploit the spatial distribution of charges or functional groups via the protein structure to achieve optimality. Strain in the graphene crown ether pore modulates a single collective variable, acting as a model for this optimization process.


The investigation of 18-crown-6 and other crown ether pores, as well as pores generally, with respect to strain (the experimental extraction of χ and determination of the local interaction environment), voltage, and ion concentration will shed further light into the behavior of selectivity and optimal transport. These characteristics underlie the mechanism of physiological ion channels and are of substantial interest in biology and medicine (1, 25). This is paralleled by a broad applicability to technologies for the sensing and sequencing of biomolecules (2628), molecular sieving and desalination (29, 30), as well as batteries, fuel cells, and other energy harvesting devices (31, 32). In all of these applications, the atomic thickness and physical strength of graphene provide a unique advantage (27, 28), facilitating the detailed examination of molecular processes such as ion dehydration (18, 33). Mechanical modulation of ion transport will thus lend itself to filtration and other technologies.

Our results indicate the ideal pore characteristics for optimal transport: Minute changes in pore size or charge give rise to a colossal change in ionic conductance as transport becomes barrierless, but further change will decrease the conductance. These systems also compose a simple and transparent platform to understand the evolution of biological ion channels. Working by analogy, it is possible to identify structural features and charge distributions that maximize selectivity while maintaining high permeation rates. These considerations, in turn, constrain which conformational transitions may be leveraged to regulate transport (34). Layered crown ether graphene pore heterostructures—mimicking aspects of the selectivity filter in the potassium channel from Streptomyces lividans (KcsA)—are likewise promising for workable models of biomolecular channel dynamics. These systems should be quite sensitive to ionic conditions, potentially gating or doping graphene to access broad regimes of behavior. The susceptibility to transverse strain is also generally interesting for artificial porous membranes as it gives otherwise inaccessible information regarding the local electrostatic and hydration conditions. This platform will thus open new and general opportunities for the study, quantification, and tuning of ion transport at the nanoscale.


We determined point charges using both localized basis and extended plane-wave DFT schemes. The localized basis calculations (B3LYP/6-31G*) use a truncated model of the graphene pore alongside the Gaussian09 code (35), affording point charges through CHELPG electrostatic potential fitting (36). In contrast, for plane-wave DFT, we used Vienna Ab initio Simulation Package (VASP) (37) and a dispersion-corrected scheme (PBE-D2/PAW) for a graphene supercell (1.71 nm by 1.73 nm) containing an 18-crown-6 pore, with periodic images separated by a 2.0-nm vacuum layer (5 by 5 by 1 Monkhorst-Pack mesh, 500-eV cutoff) (38). In this case, point charges were determined using Bader analysis (39). A complementary set of plane wave calculations used a similar pore within a cell (2.96 nm by 2.99 nm), described using CP2K (40) and a mixed Gaussian plane-wave basis scheme (PBE-D3/GPW/DZVP/GTH pseudopotentials, 5442-eV real-space grid cutoff) with Γ-point sampling and RESP charge fitting (41). Point charges were derived from optimized geometries, defined where self-consistent energy convergence were below 1.0 × 10−4 eV and forces were below 0.1 eV nm−1. CHELPG fitting yielded qO = −0.24e without a K+ present in the pore and −0.23e when a K+ is present. The complementary plane-wave calculations confirmed these assignments. The Bader analysis yielded qO = −1e, an upper bound for local partial charges because of its strong bias toward chemically intuitive, but not electrostatically representative, distributions. Both MD (below) and DFT gave an α of about 2 (DFT yields a slightly higher value), indicating that the pore rim has an enhanced response to the strain. Reference fragments of qO = −0.54e included calculated CHELPG carbonyl charges for the Ac-Gly-Gly-CH3 peptide (B3LYP/6-31G*), as well as standard OPLS (qO = −0.50e) and CHARMM (qO = −0.51e) backbone carbonyl parameters. Graphene should have oxygen partial charges lower than any of these fragments’ values, as well as the isolated crown ether’s oxygen partial charges, because of the delocalized nature of the electrons in the carbon sheet (see the Supplementary Materials for details and extended methods).

We simulated ion transport through the pore using all-atom MD via the NAMD2 code (42) with an integration time step of 1 fs. We took a simulation cell aspect ratio of 1.2—the golden aspect ratio—where the simulation captures both access and pore resistances (43, 44). While access resistance should be less than a 10% correction for some parameter regimes, it must be taken into account when understanding how the drift-diffusion limit is approached. Its exclusion by choosing inappropriate (either too small or a poor aspect ratio) simulation cells hinders the comparison with the drift-diffusion limit (also see the Supplementary Materials).

We used the adaptive biasing force method (45) to calculate the equilibrium free-energy ΔF profile for a K+ going through the pore. We computed ΔF versus z in a cylindrical region with a radius of 0.1 nm and a length of 3 nm centered around the pore. As a K+ approaches the pore, it experiences either a free-energy barrier or a potential well. Assuming a small applied voltage V, the current over a barrier ΔF > 0 is (46)IK=2qKkineΔF/kBTsinh(qKV/2kBT)(5)and the current across a potential well ΔF < 0 isIK=2qKkinkout c eΔF/kBTsinh(qKV/2kBT)kin c+kout eΔF/kBT cosh(qKV/2kBT)(6)where kin and kout are the rate constant for incoming and outgoing ions, respectively, and c is the ion concentration in bulk solvent. Both of these expressions give Eq. 3 for large ΔF (see the Supplementary Materials for an additional discussion and results regarding the influence of the free-energy barrier, including a multiple-barrier rate model).


Supplementary material for this article is available at

Section S1. Simulation methods

Section S2. Equilibrium free energy and many-body effects

Section S3. Ion transport mechanism: Knock-on versus drift-diffusion

Section S4. Ion transport through multiple barriers

Fig. S1. Models used for charge calculations.

Fig. S2. Free-energy profile when only one K+ and one Cl are present in the solution.

Fig. S3. Free-energy dependence on strain for qO = −0.24e near the pore.

Fig. S4. Dehydration effects and the dielectric constant in confined geometries.

Fig. S5. Residence time for a K+ translocating through the crown ether pore in 1 M KCl solution with various values of qO, strain, and applied voltage.

Fig. S6. Delay time between one K+ leaving the pore and another K+ replacing it for various values of qO, strain, and applied voltage.

Table S1. Geometric parameters and Bader charges from plane-wave DFT.

Table S2. Oxygen point charges within an 18-crown-6 graphene pore using electrostatic potential fitting (CHELPG and RESP) and Bader analysis.

Table S3. Force-field parameters to calculate the total bonded energy.

Table S4. Pore strain, supercell edge lengths (ℓx, ℓy), nominal pore radii (rn), and geometric pore radii (rp) as a function of supercell strain.

References (4767)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: S. S., J. E., and C. R. acknowledge support under the Cooperative Research Agreement between the University of Maryland and the National Institute of Standards and Technology Center for Nanoscale Science and Technology, award 70NANB14H209, through the University of Maryland. Author contributions: S.S. analyzed the ion transport behavior and performed the MD simulations. J.E. and C.R. performed the DFT calculations. S.S., J.E., and M.Z. formulated and derived the mathematical expressions. All authors developed the ideas and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Certain commercial products are identified in this paper in order to specify the theoretical procedure adequately. Such identification is not intended to imply recommendation or endorsement by the National Institute of Standards and Technology nor is it intended to imply that the products identified are necessarily the best available for the purpose. Data and materials availability: All data necessary to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article