Ionic Coulomb blockade controls the current in a short narrow carbon nanotube Open Access

Biological ion channels are vital to living creatures and are widely targeted in pharmacology through the development of channel-specific drugs.^1–3 Artificial ion channels, their close analogues, are expected to shape the next generation of nanodevices for DNA sequencing,⁴ nanoscale energy harvesting,^5,6 desalination,⁷ gas⁸ and isotope⁹ separation, and molecular sensors. However, understanding, predicting, and controlling charge transport in ion channels still present formidable challenges to both nanotechnology¹⁰ and biophysics.^1,11 Critical gaps recently identified^12,13 in our understanding of ionic conduction in nano-channels include correlated transport and the conditions needed for the enhanced selectivity of particular ions. Coulomb blockade^14–18 is perhaps the best-known phenomenon believed to underlie the selectivity and conductivity properties of nano-channels.

The electrostatic potential of the system is shown in Fig. 1(a) as global minima of the discrete set of parabolic potentials given by Eq. (1) for different values of n vs Q_ex. The minima of the global potential correspond to stable states with corresponding minima in the permeating current and a fixed number of mobile charges in the channel, resulting in occupancy plateaus as shown in Fig. 1(b). At the peaks of the potential, there are fluctuations between n and (n + 1) in the number of mobile charges, yielding peaks in channel current at the edges of occupancy steps. The lines correspond to strong (dashed) and weak (solid) ICB effects.¹⁶ Full blockade (i.e., zero conduction at minima) is attained when U_c ≫ 1 k_BT, but the effects are much weaker when $Uc∼O(1kBT)$⁠. This approach immediately demonstrates the deep analogy that exists between ICB and electron transport in single electron devices,^20,21,25,26 and it predicts Coulomb blockade and oscillations in ionic transport through narrow channels.

The analytic predictions in Fig. 1(b) show that the current is blocked for small Q_ex, while current oscillations and the corresponding occupancy steps as a function of wall charge occur for larger negative values of Q_ex. A similar behavior was predicted analytically¹⁸ in long channels conducting hydrated ions. As mentioned above, changing U_c affects the shapes and locations of the peaks/steps. This effect was noted in Ref. 15 and described as weak/strong Coulomb blockade; see also Ref. 27 and the supplementary material.

Ionic Coulomb blockade without oscillations was initially predicted in 2D nanopores using a kinetic model.¹⁴ ICB and current oscillations were later predicted in a wide range of biological selectivity filters (including KcsA and NaChBac²⁷) and in carbon nanotubes.^18,27

Feng et al.¹⁷ observed an experimental indication of ICB in 2D nanopores, but the strongly nonlinear current–voltage characteristics that they observed should probably not be considered as an unequivocal demonstration of ICB.²⁸ A systematic sequence of mutation experiments and an elaborated analysis²⁹ showed that ICB oscillations apparently underlie the observed mutation-induced transformations of selectivity and conductivity in biological sodium and calcium filters. However, full ICB oscillations could not be observed, either in the experiments or in the corresponding molecular dynamics (MD) simulations. The latter constitute an important link between theory and experiment. In particular, the experimental observation of ICB¹⁷ was confirmed by MD analysis³⁰ demonstrating the trapping and blocking of conducting ions in graphene-embedded 18-crown-6 ether pores followed by sharp increases in conduction via a knock-on mechanism at large bias voltages. In Refs. 17 and 31, blockade was measured as a function of voltage, preventing observation of ICB oscillations.

A single strong current peak as a function of rim charge Q_ex, and corresponding transitions of the ion translocation mechanism from the diffusion regime to knock-on in these structures, was observed in Ref. 32 at Q_ex = −3.24e. Similarly, a peak in the potassium conductivity through a charged graphene nanopore was observed in Ref. 33. Closely related changes in the selectivity and conductivity mechanisms, depending on the charge and the arrangement of a functionalized graphene nanopore,³⁴ also suggest that ICB can be observed in nanopores. However, observation of ICB current oscillations has remained elusive.^17,28

Here, we seek evidence of ICB oscillations in a charged carbon nanotube (CNT). These structures^{18,27,35–38} are in themselves of special interest because they are widely used in nanofluidics, e.g., to sequence single-stranded DNA³⁶ and to build field-effect transistors,³⁹ ion pumps,⁴⁰ and water⁴¹ pumps. They can be integrated with biological nanopores.⁴² Their structural and functional similarity to biological selectivity filters suggests that CNTs may serve as useful biomimetic devices to reproduce biological functionality in nanotechnology systems.^35–37 

We approach the problem by implementing all-atom molecular dynamics simulations of ionic conduction in a short, single-wall, carbon nanotube with charged walls.

To study ion conduction through the CNT, we have developed an atomistic model of the system shown in Fig. 2, including graphene sheets and a CNT connecting two bulk solutions with TIP3⁴³ water molecules and Joung–Cheatham⁴⁴ ions. The model was built using the J-OCTA software⁴⁵ and a mixed force field: (i) Amber20⁴⁶ for water, ions, and graphene sheets and (ii) the DREIDING force field⁴⁷ for the CNT. The simulations were performed with the GROMACS software⁴⁸ using the NPT ensemble with a 2 fs time step.

The potential of the mean force (see above and the supplementary material) was calculated using the well-tempered metadynamics approach⁴⁹ implemented in the GROMACS Colvars module.⁵⁰ One or two ions for estimating the 1D or 2D potentials were kept in a cylinder of 10 Å radius and placed on the tube’s axis of symmetry. Other ions were excluded from the cylinder entrance by a repulsive potential.

Theoretical analysis considered ionic conduction through a set of binding sites connected sequentially.²⁷ The conductivity of the system was calculated by analyzing the fluctuation intensity in the number of ions at each site, taking into account correlations of ionic fluctuations at different sites. The susceptibility was calculated using linear response theory, by expanding the mean occupancy of a binding site to the first order approximation with respect to the applied potential. The predictions of Fig. 1 were made for 5 identical binding sites able to accommodate up to 4 ions. The theory was also applied to estimate the model parameters.

The graphene sheets separating the CNT from the bulk solutions have nanopores that allow for non-selective passage of ions from the bulk to the CNT. These nanopores are almost round in shape, with vertical and horizontal openings of 13.3 and 12.85 Å, respectively. Their charged rims play the role of the charge found at the pore entrances of biological channels, thus mimicking the internal and external vestibules of biological channels. Their geometry in the present study remains fixed. The total charge of the rim atoms is Q_np = −1.92e. A constant electric field of either 0.04 or 0.08 V/nm is applied in the z-direction. The charged carbon atoms of the CNT are in the form of rings, mimicking the binding sites of the biological KcsA potassium channel. Note that, although the charge distribution might affect the magnitude of the permeating current, it would not be expected to affect the ICB phenomena that are of interest here.

The model is conceptual. It was developed to seek evidence for ICB oscillations. Because they were predicted^{15,16,18,27,29,52,53} to occur over a wide range of parameters, we allowed for variations in the charge and Lennard-Jones parameters of the atoms. The distribution of the compensating charge and further details of the model are provided in the supplementary material. We now discuss the results obtained.

Simulations of the current through the CNT as a function of wall charge Q_ex for two different applied fields yielded the results shown in Fig. 3. Several features predicted by the theory are immediately evident. A region of complete Coulomb blockade occurs at small negative values of Q_ex. The occupancy of the CNT (shown by blue circles) is an integer and changes from zero to one, and from one to two, in steps as Q_ex becomes more negative. Current peaks corresponding to these occupancy steps are clearly resolved, as indicated by black arrows in the figure. According to theory,^{15,16,18,27,29,52,53} the first peak is related to single-ion transitions as shown in video No. 1 of the supplementary material, while the second peak corresponds to the 2 → 1 knock-on mechanism shown in video No. 3 of the supplementary material. All of these features have long been anticipated on the basis of ICB theory^{15,16,18,27,52} but had not hitherto been demonstrated.

However, a number of unexpected features are also evident in the MD simulations. First, the observed separation of the occupancy steps in terms of the wall charge is significantly larger than the value Q_ex ≈ 1e expected from ICB theory for monovalent ions.^{15,16,18,27,29} Second, additional current peaks (secondary peaks) may appear for constant values of occupancy as indicated by red arrows in the figure. The comparison of the currents for the two different applied fields in Figs. 3(a) and 3(b) shows that the secondary peaks are not always present. For example, the (2 → 1)′ peak is suppressed for F = 0.04 V/nm. We now consider in more detail the ionic dynamics responsible for each set of peaks, considering them in turn from right to left, and we will show that they correspond to two different types of conduction.

The potentials of the mean force (PMFs) for a single ion transition (1 → 0) are shown in Fig. 4(a) for three different values of the wall charge Q_ex. Details of the calculations are provided in the supplementary material, Sec. III. A key feature of the equilibrium PMF shown in Fig. 4(a) is the existence of high dehydration barriers at the channel entrance and exit. An increase in Q_ex reduces the entrance barrier but, at the same time, deepens the potential well within the channel, thus making it difficult for the ion to traverse the CNT. An external field [Fig. 4(b)] biases the potential, helping ions to slide through the CNT. It also reduces both dehydration barriers, making the ionic transition nearly barrierless for the optimal values of Q_ex = −1.6e and F = 0.04 V/nm. Sub-optimal values of the wall charge increase either the entrance (Q_ex = −1.4e) or exit (Q_ex = −1.8e) barrier, thus reducing the current. The optimal applied field corresponds to an equal probability of the CNT being occupied by either one or zero ions, thereby maximizing the current peak. These features are in agreement with the predictions of ICB theory.^{15,16,18,27,29} An illustration of this type of transition is given in video No. 1 of the supplementary material.

For comparison, the 1D PMF for the transition peak (2 → 1)′ is shown in Fig. 4(b) for three values of applied field. First, we notice significant deepening of the potential well in the CNT, of up to −20 k_BT so that, even in the presence of the F = 0.04 V/nm applied field, the height of the barrier at the CNT exit is $∼10kBT$⁠. As a consequence, the ion cannot leave the CNT without the assistance of a second ion. At the same time, the occupancies of the CNT and of the channel between two graphene sheets practically do not change (see the blue circles and red triangles in Fig. 3), indicating that the second ion enters the CNT only very briefly to knock out the first ion. Additional studies discussed in the supplementary material, Sec. IV, confirm that such a behavior represents a remote knock-on conduction mechanism^54,55 when the second ion knocks the first ion out of the CNT while still remaining at the channel mouth. This process becomes possible due to the enhanced ion–ion interaction within the CNT.⁵⁶ We note that such a mechanism is likely to be responsible for pushing ions from the S1–S3–S5 configuration to the S0–S2–S4 configuration of the KcsA filter,⁵⁷ where S0 and S5 are the external and internal vestibules of the KcsA filter, which correspond to the channel mouths of our model; see the supplementary material, Sec. IV, for further discussion.

Note in Fig. 3(b) that the exit barrier, of $∼10kBT$ for F = 0.04 V/nm, is too high even in the presence of the remote knock-on mechanism so that the (2 → 1)′ conduction peak is correspondingly suppressed. When the applied field is increased to F = 0.08 V/nm, however, the exit barrier decreases to 3 k_BT and the (2 → 1)′ peak is then clearly resolved. Thus, the remote knock-on mechanism is able to explain both the origin of this secondary current peak and the fact that it does not always appear.

A very similar physical picture underlies the dynamics of the (2 → 1) and (3 → 2)′ transition peaks located at Q_ex = −3.6e and= −4.4e, respectively. The (2 → 1) peak belongs to the set of classical ICB oscillations predicted by the theory, similar to the (1 → 0) peak discussed above. It is located at the population step from 1 to 2 corresponding to the intersection of the second and third branches of the electrostatic energy (1) shown in Fig. 1(a). In addition, as expected during this transition, the CNT contains either one or two ions with equal probability; see video No. 3 of the supplementary material.

Finally, the peak (3 → 2)′ located at Q_ex = −4.4e occurs at nearly constant occupancy of the CNT and is therefore one of the secondary peaks. During this transition, the third ion enters the channel mouth and remotely knocks out the rightmost ion from the CNT. Similarly to the case (2 → 1)′, such intrusions are very brief in time and the CNT predominantly contains only two ions. However, the channel mouth occupancy may increase significantly as shown by the cyan diamonds in Fig. 3, supporting the idea of this remote knock-on mechanism; see also video No. 4 of the supplementary material. The intensity of these peaks is determined by a nontrivial interplay of the applied field, dehydration barrier, and wall charge. They can be partially or completely suppressed for some model parameters, including changes of the applied field and/or scaling of the ionic charge.⁵⁸ 

We note from the observation of Li et al.⁵⁹ that a polarizable force field is needed for ions to enter the channel. In the present work, we have found, however, that the charge on the CNT reduces the dehydration barrier sufficiently for the charge to enter, without use of polarizable force field. The origin of this apparent discrepancy is currently unclear, but it is evident that there is some interesting physics here that we intend to explore in the near future.

In conclusion, our all-atom MD simulations of the ionic current through a short, single-walled, CNT with charged walls have revealed ICB oscillations and an occupancy staircase. Although both effects had been predicted by ICB theory,^{15,16,18,27,52} neither had previously been demonstrated in a CNT.

We have shown, however, that the current peaks as a function of Q_ex are attributable to two distinct conduction mechanisms. The first, corresponding to the set of transitions denoted by {n + 1} → {n}, occurs at occupancy steps, as predicted by ICB theory. During this type of transition, the CNT is occupied by either n or (n + 1) ions, with equal probability, i.e., P_n ≈ P_n+1. The second mechanism corresponds to the transitions denoted by ({n + 1}→{n})′, which occur at constant occupancy of the CNT. This latter type of transition contradicts the expectations based on ICB theory. We showed that the CNT is then predominantly occupied by n ions (P_n+1 ≪ P_n), and knock-on is accomplished remotely by an ion arriving transiently at the channel mouth. This type of transition varies in probability and can be attenuated or eliminated for certain choices of model parameters.

These discoveries confirm the conjecture that ICB plays a fundamental role in controlling ionic conduction in artificial, as well as biological, nano-channels. The new insights pave the way for systematic studies of a wide range of important physical phenomena through their individual effects on ICB oscillations. These include dehydration, wall fluctuations, local dielectric permittivity, wall functionalization, polarization, channel mouth geometry, and other physical properties. Systems of the type described by Pang et al.³⁹ and Rabinowitz et al.⁴⁰ appear to be suitable for experimental tests. We comment also that the system considered here is an excellent candidate for mimicking the celebrated selectivity and conductivity properties of KcsA-like biological channels.^60–63 We aim to extend our studies of ICB oscillations to encompass longer CNTs, the role of the charged vestibules at the CNT entrances, and DNA sequencing in short CNTs.

The results reported promise to facilitate the design and optimization of the controllable nanoscale ionic devices needed for, e.g., blue energy harvesting,⁶⁴ DNA sequencing,³⁶ nanofluidic neuromorphic computing,⁶⁵ ionotronics,⁶⁶ medical dialysis,⁶⁷ biosensing,⁶⁸ and biointerfacing.⁶⁹ 

The supplementary material PDF document provided with this paper includes (i) a summary of the statistical theory of ionic Coulomb blockade; (ii) a description of the charge distribution on the CNT and on the adjacent graphene sheets; (iii) a discussion of the use of the potential of the mean force in molecular dynamics computations; and (iv) a detailed discussion of the remote knock-on permeation process. In addition to the PDF, there are the 4 videos alluded to in the main paper showing how one or more K+ ions permeate the CNT. Each of them is for rim atom charge Qnp = −1.92e.1.Permeation dynamics for a single-ion transition (1 → 0): Qex = −1.60e and voltage F = 0.08 V/nm.2.Permeation dynamics for a knock-on transition (2 → 1)′: Qex = −2.20e and voltage F = 0.08 V/nm.3.Permeation dynamics for a knock-on transition (2 → 1): Qex = −3.60e and voltage F = 0.08 V/nm.4.Permeation dynamics for a knock-on transition (3 → 2)′: Qex = −4.40e and voltage F = 0.08 V/nm.

We are grateful to the J-OCTA team for providing the J-OCTA software, for support, and for valuable discussions. This work was supported by the Leverhulme Trust (UK) (Grant No. RPG-2017-134) and by the Engineering and Physical Sciences Research Council (UK) (Grant No. EP/M015831/1).

The authors have no conflicts to disclose.

William A. T. Gibby: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Miraslau L. Barabash: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Igor A. Khovanov: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Writing – review & editing (equal). Dmitry G. Luchinsky: Conceptualization (lead); Formal analysis (equal); Funding acquisition (equal); Investigation (lead); Methodology (equal); Supervision (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Peter V. E. McClintock: Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.