Ab Initio Properties of Hybrid Cove-Edged Graphene Nanoribbons as Metallic Electrodes for Peptide Sequencing via Transverse Tunneling Current.

Recently synthesized metallic cove-edged graphene nanoribbons are considered for use as one-dimensional (1D) electrodes for ideal atomistically resolved recognition of amino acids. To this purpose a narrow nanogap device is considered, and the transversal tunneling current flowing across it is calculated during the translocation of a model Gly homopeptide using the nonequilibrium Green function scheme, based on density functional theory. We show that the signal collected from the metallic spin states is characterized by a double peak per residue in analogy with the results obtained with 1D graphene nanoribbon template electrodes. The presented results pave the way for experimentally feasible atomistically resolved tunneling current recognition using metallic edge engineered graphene electrodes obtained by bottom-up fabrication strategies.

Introduction

Future strategies for personalized medicine require an exceptional effort to push ahead proteomics, the protein counterpart of genomics, and face its challenges that are intimately related to the huge population of the human proteome, which is almost 2 orders of magnitude larger than the human genome. Therefore, sequencing the primary structure of proteins and peptides is crucial to evidence mutations or post-translational modifications that affect the protein three-dimensional (3D) conformation and can result in wrong behavior in human cells and, hence, human illnesses. However, such a task is exceptionally challenging these days because the current sequencing methods, which rely on degradation of the protein into small peptides and on the recognition of such fragments on the basis of the current content of the protein data banks, are time-consuming and limited because mutations are often missed.

Therefore, new sequencing methods and devices inspired by nanotechnology strategies that have been successfully employed for DNA sequencing[1−4] are also being considered for proteins.[3−12]

Basically, two methods are considered: the first one exploits the blockage of the ionic current flowing along a nanopore axis, either biological or inorganic, during peptide translocation; here, the recognition benefits from the support of signal processing and machine learning techniques,[5−10,12,13] but the needed resolution of single amino acids (AAs) is still to come because the axial ionic current is affected by a large piece of the protein.[14]

The second method employs the measurement of the transversal tunneling current flowing between two nanoelectrodes during protein translocation across a nanogap[15] and has the advantage that the signal measured (and the AAs’ recognition) comes directly from the quantum mechanical nature of the chemical and physical properties of the piece of molecule occupying the nanogap at a given time.[11,16,17]

Both these schemes require a controlled translocation dynamics of the protein (in its primary structure state), which is still an open technological issue, but the second one allows, in principle, the measurements of AAs-related signals with atomistic resolution provided that two-dimensional (2D) or one-dimensional (1D) electrodes are employed.

Using the transversal tunneling current, one can obtain atomistic resolution by using 2D electrodes, and graphene has been considered as a natural choice for that.[4,18−26]

Recently we have shown that an ideal device conceived as an array of subnanometer gaps in graphene nanoribbons (GNRs) was able to sense single peptide bonds (PB) with atomistic resolution with clear specific features from the atoms involved.[27−30] However, these proof-of-principle results have been obtained using GNR electrodes that, besides being very narrow and thus still out of the present technology, were considered in their metastable metallic phase, not the ground-state one. These results concerning some small model peptides have been obtained in the context of a nonequilibrium Green function approach based on density functional theory (DFT-NEGF) and the well-known Landauer–Büttiker formula[31−34] (see section ). The described phenomena and results are intimately related to the electronic properties of the device and the electrodes involving the pseudo-π and pseudo-π* orbitals of the GNR electrodes and the PB molecular orbitals (MOs), with a crucial role also played by the hybrids formed by the resonant p orbital of the PBs along the GNR direction.

In these articles we chose the unpolarized metallic phase of two semi-infinite zigzag hydrogenated graphene nanoribbons,[35] 2-ZGNR, as a paradigmatic template of metallic or half-metallic 1D graphene-based electrodes. Indeed, while it is known that the electronic ground state of such ribbons is a semiconducting ferromagnetic state,[36,37] the spin-polarized ground state might undergo a transition to a semimetal when the ZGNR is doped with N atoms[38] or when a transverse electric field is applied.[39] Moreover, bottom-up strategies of edge engineering, such as cove-edged and zero-mode superlattices, have been employed recently to obtain metallic zigzag nanoribbons.[40,41] Therefore, in the present article we are showing that, if ground-state cove-edged metallic GNRs are employed as electrodes, the main features of the previous findings obtained with template metallic nanoribbons still hold.

Computational Methods

The ideal device and the system employed is schematically drawn in Figure . Here, the two electrodes are two semi-infinite hydrogenated zigzag antisymmetric-even cove-edged GNRs with eight zigzag carbon rows (8CEZGNR)[40] (the width is dw ≈ 1.6 nm including the hydrogens). The central device region contains pieces of the cove-edged ZGNR that are connected to two pieces of zigzag GNR with six carbon rows (6ZGNR). The gap between the two pieces of ZGNRs is dG = 5 Å, consistent with the previous literature,[27,29,30] and the Gly homopeptide translocates between them.

Figure 1

Ideal device made of an array of nanogaps and GNRs onto a nanopore. During the peptide translocation across the nanopore, a tunneling current trigger signal, related to the peptide backbone, is collected from the central nanogap.

Initially, the gap device is fully relaxed in the context of DFT using the Quantum ESPRESSO package:[42] electronic structures are calculated in the framework of a plane-wave basis set expansion, using a plane-wave energy cutoff of 70 Ry and an artificial smearing of 0.005 Ry. The PBE (Perdew–Burke–Ernzerhof)[43] functional is used for the electron exchange-correlation potential V[n(r)] and norm-conserving pseudopotentials built with the Troullier–Martins scheme[44] are employed. Self-consistent calculations are performed until the convergence threshold of 0.001 au for the total force is achieved.

The initial stage is the peptide translocation across the gap simulated by nonequilibrium steered classical molecular dynamics (SMD)[45] in water in order to collect the translocation configurations that will be further processed as detailed in the following discussion. In this stage the system contains a total number of atoms N = 27 322 that includes the device, the peptide, and the water at ambient conditions. In order to perform the SMD, we first need to equilibrate the system using a sequence of NVT (T = 300 K) and NPT (T = 300 K, P = 1 atm) runs, with a total equilibration time of teq = 11 ns, using the Langevin thermostat (with a damping coefficient of 1 ps–1) and the Langevin piston method[46] as implemented in NAMD. Then, the SMD simulations are performed at T = 300 K using the same thermostat and a constant velocity steering protocol (v = 0.001 Å/fs) for a total time of tSMD = 1.5 ns. The employed SMD velocity is much higher than the one currently used in translocation experiments due to computational reasons.

The completely linearized peptides have been translocated across the gap several times in sequence employing the periodic boundary conditions perpendicularly to the GNR plane. More details on the classical molecular dynamics (MD) protocol adopted to equilibrate the system and to perform the SMD runs can be found in the previous literature.[27−30]

After having produced the SMD translocation configurations, we have removed the waters[27,47] and calculated the tunneling current for selected atomic configurations with meaningful groups located in the middle of the gap. More precisely, we have considered the following groups: the carboxyl CO and the amino NH groups (that are bonded together in the peptide bond), the side chain (SC), and the middle-bond configurations between these three groups, namely, CBSN (center of the bond between the carbon bonded to the side chain and the nitrogen of the ith residue), CBSC (center of the bond between the carbon bonded to the side chain and the Cα atom of the ith residue), and PB (center of the peptide bond between the ith and the jth residues), with i(j) = 1–5 indicating the amino acid in the peptide central subsequence in order of translocation.

The quantum transport calculations have been carried out in dry ambient conditions since the presence of water does not affect the transport properties of the ZGNR nanogap because of the hydrophobic character of graphene.[27,29,47] According to the relevant literature, the experimental evidence of modified graphene conductivity due to edge doping by water is most probably due to the different status of the edges, to the author’s knowledge.[48] Lastly, the gap sizes considered here prevent the passage, at the same time, of the peptide chain and of water molecules: the last ones, however, can be easily filtered because they give a tunneling current signal 1 order of magnitude smaller than the ones for the peptide.[27,29]

The selected configurations have been further relaxed at T = 0 K, in the context of spin-polarized local density approximation density functional theory (LSDA-DFT),[49,50] using the SIESTA package[51] with a large threshold for the atomic forces, that was set to 0.1 eV/Å,  to conserve the thermal disorder of the SMD configurations. The root mean square deviation (RMSD) per atom measured between the quantum-relaxed and the SMD configurations is approximately 2.9 × 10–2 Å/atom.

This relaxing stage is necessary for two reasons: first the interactions between the ZGNR leads and the peptide during the SMD are basically incorrect due to absence of the electronic (or band structure) part of the interaction between the atoms; these may be important especially for nonlocal wave functions such as the ones of the ZGNR orbitals. Second, although the peptide configurations selected are the ones with minimum strain, we still have some residual strain, due to the large value of the SMD velocity, that needs to be relaxed.

After the relaxation, we have calculated the transmission function and the tunneling current according to the DFT-NEGF scheme[32] as implemented in the TRANSIESTA code;[33] although it is a ground-state theory and not a steady-state one, DFT-NEGF is the most popular approach for steady-state transport in nanostructures and has been successfully applied in many cases with results similar to the ones obtained from formally correct steady-state methods.[52,53] According to DFT-NEGF, the spin-resolved transmission function iswhere Gσ(ε) = is the Green’s function of the system, ΓL(R)σ(ε) = is the left(right) spin-resolved coupling function, and ∑L(R)σ is the left(right) electrode self-energy for the σ spin.[32] The related tunneling current is obtained through the Landauer–Büttiker formula[31] for an external bias voltage of V = 1 V applied along the z directionwhere f(ε) is the Fermi–Dirac distribution function, μL(R) is the electrochemical potential of the left(right) electrode, and T(ε, V) = T↑(ε, V) + T↓(ε, V).

In the last two stages, split-double-ζ basis sets augmented with polarization orbitals (DZP) are employed for H, O, C, and N atoms, using a mesh cutoff of 250 Ry, the LSDA-PZ (Perdew–Zunger)[49,50] functional for the electron exchange-correlation potential, and norm-conserving Troullier–Martins pseudopotentials.[44]

The spin-dependent tunneling currents have been calculated for all the relaxed configurations and the intermediate positions between the previous ones obtained by rigidly shifting the closest relaxed configurations. Thus, the tunneling current signal during translocation is sampled with five samples per AA. Because the typical experimental translocation time is in the range of 0.1–5.5 ms per amino acid[3] in real experiments, this corresponds to a required sampling rate of at least 50 kHz, which is definitely achievable with the currently available bandwidths of amplifiers and fast picoammeters.

Lastly, we have further analyzed the signals, looking at the atomic current branches between the atoms that are calculated aswithand the projection operator on the n site P = ∑γ|ϕ,γ⟩⟨ϕ,γ| expressed in terms of a complete localized basis set {|ϕ⟩} in a subspace of the full device region, where |ϕ⟩ is an atomic orbital of type γ, located at the atomic site n.[34,54]

Results and Discussion

We first analyze the behavior of the 8CEZGNR electrodes. Their structure is reported in Figure , parts a and b, where the alternate up/down edge bending is evidenced. In the spin-resolved band structure reported in Figure c we evidence that the cove-edged ribbon is nearly metallic (the two spins have identical band structures), the π and π* orbitals being nearly degenerate at the Γ-point of the 8CEZGNR Brillouin zone (BZ), where we have measured a gap of Eg ≈ 3.5 meV (Figure d). This is the result of graphene BZ folding that superimposes the graphene K-point onto the 8CEZGNR Γ-point. Such a negligible gap is easily overcome at 300 K, which is the electronic temperature employed here.

Figure 2

Antisymmetric-even cove-edged GNR with eight zigzag carbon rows (8CEZGNR): side (a) and top (b) views. The spin-resolved band structure of the 8CEZGNR (c). The nearly metallic behavior can be appreciated from the magnified band structure at the Γ-point (d).

We have also found (not reported) that this small gap increases, even though at a small extent, for the narrower 6CEZGNR (a cove-edged GNR with six carbon rows) as expected. Interestingly, the spin-polarized PBE scheme as implemented in Siesta predicts an even more marked, opposite behavior that is unphysical.

The nearly degenerate orbitals at the Γ-point are reported in Figure for the two spin states: the valence band spin states have pseudo-π character with higher densities at the opposite edges similar to the “conduction band” pseudo-π* spin states but with opposite distribution of the up and down spins.

Figure 3

Spin-resolved orbitals of the lower and higher nearly degenerate states at the Γ-point of the 8CEZGNR band structure in Figure .

As mentioned, the device region also contains two pieces of 6ZGNR that are semiconducting for both spins with a band gap of Eg ≈ 0.35 eV. The antiferromagnetic spin-polarized ground state of narrow ZGNRs has been demonstrated to be stable up to room temperature (RT) for ZGNRs produced by nanolithography from graphene sheets.[55]

Of course the transmission function of the 8CEZGNR–6ZGNR nanogap is nearly zero in the −0.5 to 0.5 V range for a bias of 1 V between the electrodes (not shown), the subnanometer gap being sufficiently wide to avoid any meaningful tunneling of the electrons when the gap is empty: indeed, we have measured a dark current on the order of the femtoamps, 5 orders of magnitude lower than the signal (see below and the Supporting Information). The properties of this structure, however, evidence interesting features at equilibrium (no bias applied). The spin-polarized model of the 6ZGNR–8CEZGNR junction barrier between the semiconducting 6ZGNR and the nearly metallic semi-infinite 8CEZGNR should predict a Schottky barrier on the basis of the respective band structures.

For the sake of completeness we have also considered a periodic cell structure where we have four 8CEZGNR unit cells and 12 6ZGNR unit cells. The results, reported in the Supporting Information, show the existence of an energy gap in the periodic structure that, however, disappears in the present gap structure because in this case the 6ZGNR pieces are just terminals of metallic ribbons (see below).

In the DFT context (no NEGF), the gap structure, indeed, appears as a periodic repetition of finite asymmetric-even cove-edged nanoribbons with four unit cells (each unit cell contains eight carbon rows) that are terminated with two pieces of 6ZGNR, each containing seven carbon rows, through two identical junctions. The 5 Å gap is wide enough to avoid a strong superposition between the adjacent replicas of the hybrid GNRs. This hybrid nanoribbon retains the metallic behavior of the central cove-edged structure, and the Fermi level is pinned in the middle of a small density of states (DOS) peak (see the Supporting Information): the local density of states integrated in the energy range of this band is reported in Figure and shows that the metallic band extends over the entire structure.

Figure 4

Local density of states in the metallic band energy range of the 8CEZGNR–6ZGNR gap structure.

Having clarified the behavior of the 8CEZGNR electrodes and 8CEZGNR–6ZGNR gap structure, we now consider the tunneling current flowing across the gap during the translocation of a Gly homopeptide, which is reported in Figure for the central part of the peptide.

Figure 5

Tunneling current calculated during the translocation of Gly homopeptides across the nanogap. The configurations where the tunneling signals peak are reported in the insets.

In all the cases examined at various translocation configurations and spin states, we have found just one transmission channel. The signal collected is characterized by a structured peak feature per residue, basically the same behavior found when using template metallic 2-ZGNR narrow electrodes in an unpolarized regime;[27,28] indeed, each residue is characterized by a double peak around each Gly peptide bond with maxima corresponding to the center of the PBi = 2, 3, j = 3, 4 peptide bonds and around the NHi = 2, 3 configurations, with the NH groups lying in the GNR plane and is in the middle of the gap. The signal shape of the current tunneling signal is quite similar to the one obtained previously using template metallic 2-ZGNR electrodes,[27,29] even though the signal intensity is approximately 1 order of magnitude lower. The reason for this lower tunneling current is strictly related to the lower DOS in the −0.5 to 0.5 eV energy range for the cove-edged electrode than the one in the metastable unpolarized metallic narrow zigzag GNR employed previously, as shown in Figure where we compare the DOS of the spin-polarized ground and the spin-unpolarized metastable phases of the 2-ZGNR to the spin-polarized ground state of the 8CEZGNR employed here.

Figure 6

Normalized DOS of spin-polarized 8CEZGNR compared to spin-polarized and unpolarized 2-ZGNR.

The calculated signal intensity, however, is still well above the detection sensitivity limits of the currently available picoammeters and could be easily detected with standard instruments. Moreover, because the tunneling current signal during translocation is presently sampled with five samples per AA and the typical experimental translocation time is in the range of 0.1–5.5 ms per amino acid,[3] in real experiments the sampling of the tunneling current would require a sampling rate of at least 50 kHz, which is definitely achievable with the currently available bandwidths of amplifiers and fast picoammeters. However, the technological issue of a controlled translocation of peptides and proteins remains still unsolved to date.

The detailed electron coupling and transmission across the peptide can be studied by calculating the cumulative atomic bond currents injected from the left lead (electron source) into the various peptide chemical groups, namely, the CO and the NH groups involved in the PB and finally the CH2 groups containing the α-carbon and the H elementary side chain, according towhere L is the left lead and G is one of the previous groups we can employ to ideally decompose the peptide. We report in Figure the injected electron currents from the source electrode for meaningful configurations where the current peak occurs.

Figure 7

Bond currents injected from the left lead into the peptide for selected configurations corresponding to the maxima of the tunneling current. The (0), (+1), and (−1) superscripts indicate the shift from the closest relaxed configuration; for instance, PB23(+1) means that the peptide has been shifted from the relaxed PB23(0) configuration one step forward along the translocation direction.

Similar to the case study reported previously using metallic template 2-ZGNR electrodes, the major contributions to the signal come from the coupling of the CH2 groups that include the CαH of the peptide backbone and the elementary H side chain of Gly. This coupling is enhanced for those translocation configurations where the CH2 groups are located above and below the terminal 6ZGNR in the device region (the graphene plane) so as to favor the hybridization between the CH2 orbitals and the left/right pseudo-π/pseudo-π* orbitals. The Gly homopeptide, indeed, behaves as a special case due to the elementary side chain with the reduced size perpendicular to the peptide backbone. Because of this, the translocation configurations of the Gly amino acids close to the GNR plane are such that the CH2 group is placed almost symmetrically to form a bridge between the two 6ZGNR terminals of the electrodes. This behavior differs from the one of peptides with larger amino acids where a very important role is played by the CαH groups of the peptide backbone that is aligned along the ribbon direction and oriented either toward the source or the drain terminals.[29] It is important to evidence that this behavior is strictly related to the constraint applied by the narrow subnanometer gap on the peptide during the translocation. It should be emphasized that the required manufacturing of such a device requires an exceptional control at the atomistic level that is not far from the present nanotechnological capabilities but, maybe, still to come. We have already discussed in the previous literature and for the template electrodes how the tunneling signal could be affected by the gap size, i.e., by the precision of the gap manufacturing.[28] In the present case, with “real” metallic electrodes, we have considered the same problem by looking at the PB23 configuration as a test case showing that the signal drops with the gap size by approximately 1 order of magnitude per angstrom (see the Supporting Information). Therefore, we think that for gap sizes above 8 Å the bias should be increased in order to have reliable and stable current signals. Because the electrodes are nearly 5 nm apart in the present device, this circumstance should not represent a serious drawback. Contrarily, one should engineer and study new electrodes in order to increase the DOS at the Fermi level.

Conclusions

We conclude by remarking that we have studied a newly conceived nanoribbon junction as part of a hybrid cove-edged nanoribbon gap device for peptide sequencing via transverse tunneling current across the gap during peptide translocation. The nanodevice electrode is made of an asymmetric-even cove-edged nanoribbon that can be experimentally manufactured using current bottom-up fabrication nanotechnology. The device contains a junction between the cove-edged nanoribbon and a zigzag nanoribbon with six rows. First, we have studied the nanodevice structure in the context of spin-polarized ground-state DFT, evidencing the metallic nature of this hybrid structure that is employed in the sequencer at the two sides of the subnanometer gap. Then, using the DFT-NEGF method in the elastic regime, we have calculated the tunneling current flowing across a nanogap in narrow graphene nanoribbons during the translocation of a Gly model peptide taken as reference with the previous literature where metastable spin-unpolarized metallic zigzag GNR template electrodes were employed. We have shown that the signal obtained in this new sequencer is characterized by a structured double peak per residue, where the major contributions come from the tunneling across the Gly CH2 groups that include the Cα group and the H side chain, with a minor contribution from the peptide bond groups. The signal level calculated is lower than the one obtained using the template metallic zigzag GNR electrodes but, nevertheless, still well within the measurable range of currently available picoammeters. Therefore, realistic GNR-based devices can be definitely employed as peptide sequencers because recent bottom-up strategies have allowed the synthesis of cove-edged zigzag graphene nanoribbons as narrow as the ones considered here. Of course, the signal level depends on the gap size, and therefore, the atomistic control of this parameter is critical. The present results pave the way toward the fabrication of realistic GNR metallic electrodes and devices for atomistically resolved recognition of the amino acids in peptides and proteins.