Abstract
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.
INTRODUCTION
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 (2–4)? 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.
(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.
RESULTS
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 (10–15), 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.
(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
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
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)
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.
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, 22–24), 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.
DISCUSSION
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 (26–28), 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.
METHODS
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)
SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/7/eaaw5478/DC1
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.
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.
REFERENCES AND NOTES
- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).