Water-assisted electronic transport in graphene nanogaps for DNA sequencing

TL;DR

This paper presents a theoretical study of a graphene nanogap device for DNA sequencing, highlighting water's crucial role in enhancing electron transport and enabling nucleotide detection with high sensitivity.

Contribution

It introduces a hybrid quantum-classical methodology to analyze water-assisted electron transport in graphene nanogaps for DNA sequencing, demonstrating high potential for practical applications.

Findings

Water significantly increases conductance in nanogap devices.

DNA nucleotides can be reliably detected and distinguished.

Water's electrostatic and electronic effects are crucial for device performance.

Abstract

Innovative methodologies for reliably and inexpensively sequencing DNA can lead to a new era of personalized medicine. In this work, we performed a theoretical investigation of a nanogap-based all electronic DNA sequencing device. To do so, we used a nitrogen-terminated nanogap on a graphene sheet with the environment fully taken into account. Our investigation is performed using a hybrid methodology combining quantum and classical mechanics coupled to non-equilibrium Green's functions for solving the electron transport across the device. The obtained results show that the DNA nucleotides can be both detected and distinguished in such device, which indicates that it can be used as a DNA sequencing device providing very high sensitivity and selectivity. Furthermore, our results show that water plays a major role in electronic transport in nanoscopic tunneling devices, not only from an electrostatics point of view, but also by providing states that significantly increase the conductance in nanogap-based DNA sequencing devices.

Equations

$$\begin{gathered} G^R\left(E,V\right)=\underset{\eta \to 0^{+}}{\lim }[ϵS_{\mathrm{SR}}-H_{\mathrm{SR}}- \\ {{\Sigma }_{\mathrm{L}}\left(E,V\right)-{\Sigma }_{\mathrm{R}}\left(E,V\right)]}^{-1}, \end{gathered}$$

$$\begin{gathered} I\left(V\right)=\frac{2e}{h}\int T\left(E,V\right)(f\left(E-{\mu }_L\right)- \\ f\left(E-{\mu }_R\right))dE, \end{gathered}$$

$$\begin{gathered} T\left(E\right)=Tr[{\Gamma }_{\mathrm{R}}\left(E,V\right)G^{\mathrm{R}}\left(E,V\right) \\ {\Gamma }_{\mathrm{L}}\left(E,V\right)G^{\mathrm{A}}\left(E,V\right)], \end{gathered}$$

$$G=G_{0}T\left(E_F\right).$$

$$i\left(E\right)_{N\rightarrow M}=\\ 4\frac{e}{h}\sum_{\tiny\begin{array}[]{c}n\in N\\ m\in M\end{array}}\Im\left[\left\{G^{R}\left(E\right)\Gamma_{\mathrm{L}}G^{A}\left(E\right)\right\}_{mn}H_{nm}\right]~{},$$

$$\Sigma (E^{\prime })=\frac{{\overline{T}}_x(E^{\prime })-{\overline{T}}_L(E^{\prime })}{{\overline{T}}_L(E^{\prime })}\times 100\%,$$

$$\sigma (E^{\prime })=\frac{{\overline{T}}_x(E^{\prime })}{{\overline{T}}_G(E^{\prime })}\times 100\%,$$

Full text

††journal: Carbonuspuspfootnotetext: Current address: Institute of Physics, University of São Paulo (USP), São Paulo/SP, Brazil

1 Introduction

Personalized DNA sequencing is the next frontier in health care, as it could be used to identify predisposition towards a number of genetic illnesses, and ultimately provide precision treatments [1, 2]. In order for this to be fully accomplished, however, the speed and cost of the procedure is yet to drop further [3, 4, 5], and it has become increasingly clear that this can only be achieved via novel sequencing methodologies [6, 7].

For this purpose a number of methods have been proposed [8]. One such approach is to use nanopores [9, 10, 11, 12, 13] or nanogaps [14, 15, 16, 17] acting as sieves for DNA segments either measuring electronic or ionic currents as the strands are driven through. The main idea behind an all electronic nanogap sequencing device is to measure conductance differences across the nanogap containing DNA strands. The current is then a signature of the electronic structure of each nucleotide and how its states couple with the nanogap [15, 18, 19].

Amongst possible materials that could be used for such a device, graphene [20] could be an excellent choice as a substrate, since its one-atom thickness could provide single-molecule resolution [18]. In fact, we have recently shown that functionalizing the edges of graphene nanogaps with nitrogen could enhance the sensitivity and selectivity of graphene-based devices, boosting their potential for sequencing DNA nucleotides [15]. Furthermore, we also showed that graphene can be used for detecting specific DNA sequences providing good sensitivity and selectivity [21].

A key issue in any realistic simulation of bio-sensing devices is that the molecules are embedded in an aqueous solution - typically under physiological conditions. Usually simulations addressing electronic transport in systems containing molecules neglect the effects of the solvent [14, 15, 22, 23, 24], or introduce it in the form of a classical potential that only influences the Hamiltonian electrostatically [12, 25]. A few works have explicitly included the solvent [26, 27] in theoretical simulations of molecular electronic devices, and in most cases the effect is usually seen as an overall shift of energy levels arising from a chemical gating effect, which are correctly captured by a classical electrostatic effect. Recently we have combined quantum mechanical/classical mechanical methods (QM/MM) with electronic transport to simulate the transport properties of a graphene DNA sequencing device. There we find that the effects of a solvent are important for the correct description of the problem, but explicitly including the electronic states of water does not exert a significant influence compared to a classical external potential added to the quantum mechanical Hamiltonian of the problem. In all cases the devices were conducting.

In this work we show that in the case of tunneling devices, the exact opposite effect occurs. Using a combination of classical and quantum mechanical methods [12, 25] in conjunction with the non-equilibrium Green’s functions (NEGF) formalism, we show that in a tunneling environment where electrons cross a nanogap, explicitly quantum-mechanically described water molecules play a major role in the overall transport properties of a device to the extent that a significant portion of the current passes through electronic states localized in the water molecules [11, 18, 28, 29].

2 Methodology

The device modeled consists of two semi-infinite graphene sheets with nitrogen termination and separated by 17 Å  as shown in Figure 1a. The four naturally occurring nucleotides are shown in Figure 1b: deoxyadenosine monophosphate (A), deoxycytidine monophosphate (C), deoxyguanosine monophosphate (G) and deoxythymidine monophosphate (T). Those nucleotides are placed inside the nanogap-based device, and subsequently filled with water molecules and counter ions with a concentration of 0.2 M. In order to balance the negative charge of the nucleotide, we introduce an imbalance in positive and negative counter ions in such a way that we have 3 Na and 2 Cl atoms in the simulation box. Here we have used isolated nucleotides to simulate the DNA sequencing in line with our previous works [13, 15] that have shown little or no effect in the electronic transport due to adjacent bases.

The first step in our approach is to time-evolve the system to realistically simulate dynamical effects. From the structural point of view, as there are no bonds breaking or forming, classical molecular dynamics (MD) yields good results for biological systems [12]. In our simulations the nanogap is kept fixed and only the direction perpendicular to the nanogap is allowed to change (anisotropic NPT simulation). We also apply a harmonic potential in the whole nucleotide to prevent it from moving outside the gap. The systems are simulated in a box with periodic boundary conditions.

In the employed methodology we first perform a simulation in the isothermal-isobaric ensemble (NPT) using a Berendsen barostat and Nosé-Hoover thermostat to equilibrate the density of the system at a given pressure (1 atm). We then extract one snapshot in which the height of the box in the direction perpendicular to the plane of the nanogap ($y$ direction) coincides with the corresponding average. From this initial configuration we perform an isothermal-isochoric ensemble (NVT) production run (same thermostat and restrictions) to generate the structures. Both NPT and NVT simulations were carried out with a 2 fs time step generating 10 ns of dynamics. We have used the AMBER99SB [30] force field and the SPC [31, 32] model for water, as implemented in GROMACS [33].

For each nucleotide, we then select 50 equally separated uncorrelated frames from the NVT production run. Typical snapshots are shown in Figure 1c. We then use a QM/MM methodology where each frame is divided [34, 35] into quantum mechanical (QM) and classical mechanic (MM) regions depending on the setup (see Figure 1d for one example). In this work we considered two setups, namely (I) the QM partitions consist of the nitrogenated graphene, the nucleotide and one counter-ion to equilibrate the charge; all water molecules are treated classically, as well as the remaining counter-ions. In the second setup (II) one layer of water is included in the QM partition as shown in Figure 1d. A third setup not including solvent effects is presented in the supplementary information for comparison.

For the QM partition we performed density functional theory calculations within the Generalized Gradient Approximation (GGA) as parametrized by Perdew, Burke and Ernzerhof (PBE) [36] with double zeta (DZ) basis set for the carbon atoms of the nanogap and double zeta polarized (DZP) for all remaining atoms. For the atomic core electrons we used Troullier-Martins norm-conserving pseudo-potentials [37]. The MM potential is calculated using the same force field, is smoothed out in the electrodes region and is introduced in the DFT calculation of the QM partition. The external potential and DFT calculations were done using the SIESTA package [38]. Then, we finally perform the electronic transport calculation [38, 39, 40] in these snapshots using a non-equilibrium Green’s functions formalism as implemented in Smeagol [40, 41]. The electronic transport properties were investigated via the NEGF method [39, 40, 41, 42, 43]. Within that approach, the system is divided into three parts: a scattering region (SR) enclosed by two semi-infinite electrodes on either side [44]. The main quantity that can be used to describe the system is the retarded Green’s function

[TABLE]

where $\epsilon=E+i\eta$, $S_{\mathrm{SR}}$ and $H_{\mathrm{SR}}$ are the overlap and Hamiltonian for the scattering region, and $\Sigma\left(E,V\right)_{\mathrm{L/R}}$ are the so-called self-energies, which include the effects of the electrodes onto the SR. Here the Hamiltonian for the scattering region is taken as the Kohn-Sham Hamiltonian including solvent effects via the external potential. From that, the current can be calculated

[TABLE]

where the quantity

[TABLE]

is the probability that an incoming electron with energy E, will be transmitted across the scattering region. The operators $\Gamma_{L,(R)}\left(E,V\right)=i\left[\Sigma_{L,(R)}-\Sigma_{L,(R)}^{\dagger}\right]$ define the coupling with the left and right electrodes.

In the limit of zero bias (linear regime), the conductance can be related to the transmission $T\left(E,0\right)$ via the Fisher-Lee relation [45]

[TABLE]

Finally, by projecting the transmission between any pair of adjacent sites $N$ and $M$, one arrives at the following formula:

[TABLE]

where the above sum runs over all localized atomic orbitals $n$ and $m$ of the basis set belonging to atoms in either site, and for a specific energy $E$. As we are in the limit of zero bias, we considered only right moving components for the current. For further details, we refer to the works of Okabayashi et al. [46] and of Paulsson and Brandbyge [47] and Lima et al. [48].

Finally, given the tunneling nature of the electronic transport here, and consequently its exponential dependence on barrier lengths, we calculate geometric averages of the transmission over snapshots, smoothing out fluctuations resulting in a self-averaging quantity [49].

3 Results

Figure 2 shows the average transmission spectra as a function of energy, for all nucleotides in both setups (I and II), as defined in the Methodology section. Interestingly, for setup I, the magnitude of the transmission at Fermi level is comparable to the system under dry conditions (see supplementary information), but the shape of the curve is significantly different as we are both considering an average behavior as well as including the correct electrostatics from a realistic environment. On the other hand, comparing systems I and II, we note that when the water molecules are considered explicitly (setup II), the average transmission increases for all nucleotides by several orders of magnitude. This result contrasts with previous ones [26, 12, 25] where the inclusion of water molecules only incurs in a gating effect as states coming from water lie far away from the Fermi level [25].

One way of determining the contribution of water states on these devices is through the analysis of local currents. For this purpose we selected one particular frame with transmission close to the average and calculated the bond current close to the Fermi Level (Equation 5). Figure 3d shows the total local current projected in the direction of electron flow averaged along the transverse plane. The general profile shown in Figure 3d is a significant decrease as a function of position along the device in the region that comprises the gap, i.e. $10<z<25$ Å. This is expected as we have a tunneling mechanism. We notice that the contribution due to the water molecules is of the same order of magnitude of the local current passing through the nucleotides.

One can also look at the contribution due to different system constituents. Figures 3e and 3f shows this analysis for local currents from one sheet to the other, passing through the nucleotide, and the current that bridges the gap via the water molecules that are treated quantum mechanically.

Figure 4 shows a 3D analysis of the local currents on the example of nucleotide T (see supplementary information for the corresponding graphs for all systems) using different colors to identify the currents of the different parts of the system. The pink arrows represent the current flowing between all the atoms of the sheet, blue arrows represent the current from the sheet to the water molecules and between the water molecules, orange ones represent the current flowing between water molecules and the nucleotide in either direction, and the green ones represent the current flowing between the sheet and the nucleotide in either direction and between the atoms constituting the nucleotide. Since the local current decays exponentially, for clarity we divide the system into 6 regions in the z direction - two inside the gap and two on either electrode - (see supplementary information). The arrow’s thickness in Figure 4 is normalized by the largest bond current within each region.

As can be observed, the largest currents comes from the graphene sheet (pink arrows) and are initially transmitted from the sheet to the water molecules (blue arrows). Then this current goes from the water molecules to the nucleotide (orange arrows). After that the current flows inside the nucleotide (green arrows) and is transmitted into the right electrode.

The results presented in Figures 3 and 4 show why the inclusion of water molecules in the QM part of the system has the effect of increasing the total transmission. As one can observe, the current goes through the water molecules, especially the ones close to the nucleotide and is of particular importance decreasing the barrier between the edges and the core. This is the case because the transmissions at the $E_{F}$ are not dominated by resonances coming from the nucleotide molecular levels. The alignment of such states at the Fermi level have appeared in a number of simulations (including our own) and can be ascribed to the absence of the solvent and the treatment of the nucleotide as a neutral molecule. Thus, in our case the tail of states from the base and water molecule play similar roles on transport. It is worth mentioning that for the empty gap case the transmittance is null by construction, since there will be no overlap between the states from the left and right electrodes.

Finally, we also performed an analysis of the device’s potential sensitivity and selectivity, defined as:

[TABLE]

and

[TABLE]

respectively. Here $\overline{T}_{x}(E^{\prime})$ is the average transmission for each nucleotide   (x = A, T, C and G), $\overline{T}_{G}(E^{\prime})$ is the average transmission for Guanine (used as our reference) and $\overline{T}_{L}(E^{\prime})$ is the average transmission for the nanogap in liquid without the nucleotide (the reference system), at energy $E^{\prime}=E-E_{F}$. Figure 5 shows the sensitivity and selectivity for each nucleotide in three different energies: $E^{\prime}=-0.05\ eV$, $0.0\ eV$ and $0.05\ eV$.

By means of sensitivity calculations, we estimate how the presence of any nucleotide compares to the basal conductance. As we can see in Figure 5a, the four nucleotides have extremely high sensitivities for all energies. This means that the presence of one of the bases inside the gap significantly increases the transmission.

At the same time our selectivity calculations, shown in Figure 5b for energies above and below $E_{F}$ all nucleotides are indeed distinguishable. For $E^{\prime}=-0.05\ eV$, for example, Adenine has a transmission around $150\ \%$ larger than Guanine and the other nucleotides have transmissions between $100\ \%$ and $20\ \%$ larger than Guanine at the same energy. For $E^{\prime}=E_{F}=0.0\ eV$, the nucleotides are also distinguishable, except Thymine, which has virtually the same transmission as the reference. Finally, for $E^{\prime}=0.05\ eV$ we observe that the nucleotides are still distinguishable and Thymine presents the larger selectivity compared to the reference (Guanine). These results for sensitivity and selectivity show that this device for DNA sequencing is highly sensitive and selective for the chosen gate voltages, which are close to Fermi Energy.

4 Conclusions

In conclusion we presented a full dynamical treatment to include solvent effects in the investigation of an all-electronic device for DNA sequencing using a graphene sheet containing a nitrogen-terminated nanogap. We used a hybrid methodology that combines quantum mechanics and molecular mechanics (QM/MM) coupled to non-equilibrium Green’s functions formalism (NEGF). Our results showed that water plays a fundamental role in the electronic transport properties of this system and that the inclusion of explicit water molecular states increases the total conductance by several orders of magnitude. This result is explained in terms of electron flow along states belonging to water molecules. As the system is away from resonant tunneling the tail of the wave functions of water and base give similar contributions. Finally we also showed that the sensitivity and selectivity of this device is very high, opening up possibilities for further investigations to use it as a DNA sequencing device.

5 Acknowledgments

The authors acknowledge financial support from FAPESP (Grant # FAPESP 2017/02317-2 and 2016/01343-7). A.R.R. acknowledges support from the ICTP-Simons Foundation Associate Scheme. R.G.A. acknowledges financial support from CNPq (2535/2017-1 and 437182/2018-5) and R.H.S. thanks the Swedish Research Council. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001, and Brazilian Institute of Science and Technology (INCT) in Carbon Nanomaterials, and FAPEMIG and CNPq. The authors also thank the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, which have contributed to the results reported here.