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

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 q_O = −0.21e and q_O = −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 q_O = −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 q_O = −0.54e as a representative example. An upper bound of q_O = −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 q_O should span the low end of the range −0.2e to −0.6e, as evidenced by electrostatic fits of q_O = −0.34e for the 18-crown-6 alkyl ether and q_O = −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 q_O = (−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 q_O (as an indicator of the oxygen-carbon charge separation) and atomic structure. The mechanism is normal drift-diffusion for q_O = −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 q_O = −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 q_O 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 q_O = −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 q_O and V is in Fig. 2C. In this case, the current increases with q_O between q_O = −0.24e and q_O = −0.54e while giving a turnover to decreasing current at a higher charge.

• Download high-res image
• Open in new tab
• Download Powerpoint

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, ΔF_K, is approximately the sum of a dehydration energy barrier ΔE_deh and three electrostatic terms (E_KO, E_KC, and E_KK) that capture the interaction of K⁺ with nearby charged atomsΔFK=ΔEdeh+EKO+EKC+EKK(1)where E_deh = η∑_if_iE_i is the barrier due to dehydration, f_i (E_i) 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πϵ₀ϵ(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. E_KK 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 ΔF_K 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 ΔE_deh [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 q_O 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.

• Download high-res image
• Open in new tab
• Download Powerpoint

Figure 3 shows that when q_O = −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⁺ (E_KO) 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 E_KK term will not contribute in this case. Equation 1 confirms this assignment of the peaks and wells.

This situation is inverted for q_O = −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ρμνda≈qKEρ(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 χ ≈ 100k_BT/nm and χ ≈ 140 k_BT/nm for the q_O = −0.24e and q_O = −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 q_O = −0.24e (the charge q_O = −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.14k_BT 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.7k_BT for q_O = −0.24e and about 1.0k_BT for q_O = −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 q_O = −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 q_O = −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 q_O = −0.24e, this will be the jump from about z = 0.1 to 0.2 nm, and for q_O = −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)I≈Ae−Δ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ΔII≈e(−ΔF⋆(s)+ΔF⋆(s=0))/kBT−1≈eα a s χ/kBT−1(4)where −ΔF^⋆(s) + ΔF^⋆(s = 0) ≈ dΔF/da ⋅ δa ≈ αasχ. The parameter α = 1/a ⋅ da/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 q_O = −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.

• Download high-res image
• Open in new tab
• Download Powerpoint

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 q_O = −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 = 6k_BT/χ), 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 q_O = −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.