**1. Introduction**

Sensing and analysis of DNA sequences and epigenetic structure modifications in the human genome are essential for the understanding of the mechanism of gene expression in cells and the development of diseases. Efforts to detect and map DNA sequences and modifications with 2D materials such as graphene and hexagonal BN have been made and have utilized a number of methods, such as nanopore-based ionic and transverse current, and bonding to complimentary molecules in nanopores and STM. The development of the surface-enhanced and tip-enhanced Raman spectroscopy (surface-enhanced Raman scattering (SERS) and TERS) [1–5] with nanoscale resolution [6–9] applicable to single biomolecule [10,11] analysis, as well as stimulated Raman scattering and ultrafast twodimensional infrared (2D-IR) spectroscopy of small organic molecules [11–14] opens a way to identify molecular structures in a label-free way and investigate mechanisms of interaction with the environment in biomolecules. Attention to DNA molecules in manifold applications, from sequencing, epigenetic studies, disease relation, etc., to nanobiostructures, makes identification of nucleotides and bases in a variety of states highly desirable in the non-contact label-free way of nanophotonic vibrational spectroscopy [15]. SERS is a promising method for nucleic acid [16] and protein detection [17,18]. In terms of its abundant fingerprint information, anti-interference with water, and extraordinarily high sensitivity, SERS holds many attractive advantages over other techniques. Recent research discusses probing of proteins by SERS in solution using gold nanoparticles in the

**Citation:** Zolotoukhina, T.; Yamada, M.; Iwakura, S. Vibrational Spectra of Nucleotides in the Presence of the Au Cluster Enhancer in MD Simulation of a SERS Sensor . *Biosensors* **2021**, *11*, 37. https://doi.org/10.3390/ bios11020037

Received: 20 December 2020 Accepted: 25 January 2021 Published: 29 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Department of Mechanical Engineering, University of Toyama, Toyama 930-8555, Japan

absence of immobilizing agents so that neither the tertiary structure of the proteins nor the surface properties of the nanoparticles are affected [17]. Methods for the rapid and sensitive detection of extended ( ≥100-base) nucleic acids with reduced preparation are also being developed [16]. They utilize the DNA sequence-specific assembly of silver nanoparticles labelled with a Raman tag to provide molecular recognition of the target DNA with probe orientation and hybridisation procedures found to be critical for the methods.

Another path to the identification of small organic molecules, such as nucleotides, can employ interaction with a nanopore in a 2D solid film, especially graphene [19–27], though experiments and calculations with hexagonal BN were also carried out [28–31]. Detection of the DNA/RNA mononucleotides [32–34] and nucleobases [35] with SERS techniques has also been carried out. Despite the many advantages and simplicity of the 2D nanopore sequencing, a few issues [20] remain to be resolved, including controlling the DNA translocation rate, suppressing stochastic nucleobase motions, and determining the signal overlap between different nucleobases. The graphene-based nanopores suffer from significant signal noises due to graphene hydrophobicity and trapping of DNA bases during translocation [36]. The study of the influence of the interaction force between the nucleobases and graphene nanopore on translocating molecules helps to find a solution to some of the mentioned issues.

Surface-enhanced Raman scattering (SERS) using plasmonic metal nanoparticles has been applied to characterize complex biological samples, ranging from biomacromolecules such as nucleic acids and proteins, to living eukaryotic and prokaryotic cells in whole animals for several decades now. The sensitivity of the SERS signal with respect to interactions at the surface of the plasmonic nanoparticles makes it specifically useful for studies of molecular contacts of nanoparticles in biological systems and does not require chemical or biological functionalization.

Control of the 3D structure of DNA and proteins in SERS spectra acquisition still remains a complex question. In the case of DNA molecules, a translocation of the single or double strand through the nanopore substrate has been well established. For proteins, conformation and unfolding are considered. The discrimination of protein conformations is of critical importance for identifying the unfolding states. Attempts to measure dynamic SERS of proteins [37] were carried out over a long acquisition time. Characteristic protein signatures at different time points were observed and compared to conventional Raman with SERS. TG-SERS can distinguish discrete features of proteins, such as the secondary structures, and is therefore indicative of folding or unfolding of the protein [38]. Other studies [39] mention that at a low gold-nanoparticle: protein-molar ratio, significant unfolding of the protein is observed at the surface of the gold NP. When using unfolded proteins, they can be translocated through the nanostructures [40] and nanopores [41]. Therefore, acquisition of SERS spectra of unfolded proteins at the translocation through a nanopore with an attached enhancer NP should be possible in principle. The process of unfolding and conformation should involve the interaction of the building blocks with system elements, as well as between amino acids in chain molecule and tertiary structures. For secondary and tertiary structures, the amino acids closest to the NP can have slight changes in SERS spectra due to the interactions, and simulation of such changes would help to distinguish and register them in acquired spectra of proteins.

In order to design and simulate the system that would allow for single molecule selection and detection of its state by high-resolution SERS methods, we tried to predict direction and possible range of spectral changes in the interacting molecule-SERS enhancer systems. The placement or growth of the small Au nanoparticle at the edge of the graphene nanopore is considered as shown in Figure 1. Plasmonic enhancement by the Au NP will make the vibrational signal measurable, while graphene nanopores would allow for control localization and transient conformation of passing DNA nucleotides, strand fragments, or proteins. Interaction at the Van der Waals interaction distances with both components, graphene and metal NP, has high probability of changing the vibrational spectral maps of the passing molecules. The vibrational mode changes in nucleotides with respect to the

interaction with graphene, which were shown in a molecular dynamic (MD) simulation to be present at short distances in our previous calculations. At present, the estimation of the influence of the small Au nanoparticle on the vibrational spectral modes of the cytosine nucleotide has been attempted to be clarified in MD calculations as a model of nucleotide–Au NP interaction. If the spectral modes of both components of interaction are studied relative to localization and molecule conformation, changes in nucleotide spectra due to interaction conditions can be mapped into a kind of library. We performed the initial steps by clarifying the possibility of such mapping. Bond-dependent spectra of nucleotide and NP were tested on their response to interactions to reveal the marker frequencies of the Au20—nucleotide interaction.

**Figure 1.** Schematic of the flow-through setup that allows single nucleotide to flow through in a graphene nanopore at plasmonic resonance upon the laser excitation.

#### **2. Model and Methods**

The present study considers the MD approach to the identification of nucleotides in interactions with the nanopore environment for vibrational spectroscopy. The nucleotide vibrational frequencies we obtained have been attributed to stretching, bending, or ringbreathing modes. By comparison of nucleobase spectra to the ones of nucleotides, the presence of the 2-deoxyribose in the nucleotides can be separated into spectral mapping and mode relaxation. The vibrational density of states is calculated in the transient regime during passing time through the graphene nanopore for each atom of the molecule to resolve the difference in the structure of nucleotides. In dynamics, the vibrational spectra are being evaluated from MD propagation velocities computed in the anharmonic interaction potential, in graphene, and in the molecule, with the Lennard-Jones (LJ) potential between them. Fourier transfer *I(f)* of the velocity autocorrelation function G(*τ*) is as follows:

$$\mathcal{G}(\tau) = \frac{<\upsilon\_i(t\_0) \cdot \upsilon\_i(t\_0 + \tau) >}{<\upsilon\_i(t\_0) \cdot \upsilon\_i(t\_0) >} I(f) \, = \, \int\_{-\infty}^{\infty} \mathcal{G}(\tau) \exp(-2\pi i f t) d\tau \tag{1}$$

where *τ* is the duration of correlation, *vi*(*<sup>t</sup>*0) is the velocity of the atom at time *t*0, and *vi*(*<sup>t</sup>*0 + *τ*) is the velocity of the atom during correlation time. According to the Wiener– Khinchin theorem, *I*(*f*) defines the vibrational density of states (DOS) of the system. The potential used for DNA nucleotides is the MM2/MM3 force field potential [42]. We investigate the transient interaction with graphene modelled by REBO potential [43] and the molecule-graphene interaction for C-X (X = H, O, C, N) by the LJ potential [42]. The interaction causes a shift in some frequencies of vibrations of the nucleotide DOSs. The vibrational spectra in MD exhibit shifts and intensity changes due to interaction with the environment that causes intramolecular vibrational mode dynamics.

The calculation setup is presented in Figure 2. The graphene sheet with the 1.5 nm in diameter pore at its center is shown in two projections, a top and front view, with a location of the nucleotide molecule relative to the pore plane and center. The graphene sheet is oriented in the x-y-plane, the edges along the y-axis are fixed, and the edges along the x-axis are free. The nucleotide can move with a given fixed velocity of the center of mass (vc.o.m.) in the positive z-direction that reproduces the motion of DNA fragments through the pores [44,45] in a driving constant electric field in experimental setups. The vc.o.m. is optimized to enhance the interaction force between the nucleotide

and edge of graphene pore and reduce rotation of the nucleotide in the translocation as much as possible. All atoms of the system except fixed ones are thermally relaxed to the temperatures Tgraphene = 300 K and Tnucleotide = 30 ◦C, as in Figure 2, before sampling. In the translocation process, the single-layer graphene sheet interacts with nucleotides in graphene nanopore. Graphene affects the interaction field of the passing molecule close to the pore edge. The graphene-molecule C-X (X = H, O, N, C) potential is considered a VdW one to avoid bond creation and nucleotide attachment to the pore.

**Figure 2.** (**a**) The structure and initial position of cytosine nucleotide relative to the graphene pore. The spatial orientation of the nucleotide's cyclic plane is at 30◦ to the z-x plane. Atoms are shown: C in grey, N in blue, O in red, and H in light grey. (**b**) Initial configuration of cytosine (right). Hydroxymethylcytosine (HMC) base is shown on the left as an example of corresponding atom numbering for bonds in molecular dynamic (MD) calculations. (**c**) Optimized by density functional theory with generalized gradient approximation (DFT GGA) calculations Au20 nanoparticle. (**d**) Reaction coordinates for stretching and bending in velocity correlation calculation use the →*v* projections along the bonds.

The relative size and weight of the nucleotides as compared with nucleobases reduce rotational motion inside the pore [46–49] that is further controlled by translocation velocity. For the given graphene-nucleotide setup, evaluation of the vibrational spectral maps of the DNA nucleotides, cytosine, thymine, adenine, and guanine, has been initially carried out in reactive coordinates by accumulating spectra of individual atom of bonds obtained by autocorrelation functions of Equation (1). The resulting spectral tables (Tables S1–S8, Figure S1) of each base and nucleotide are included in Supplementary Materials. It was confirmed that the nucleotides–graphene pore interaction does not generate recognizable changes in the spectral map of graphene pore atoms in our system. However, the conformation of the molecules responds to the interaction flexibly, exhibiting spectral changes. The center of mass (c.o.m.) of nucleotides was adjusted relative to the center of the pore along the y-axis to maximize force at the pore edge. The orientation of the cyclic planes of each type of nucleotide was selected to have the same tilt angle for comparison of spectra for various nucleotide interactions. An initial position of the nucleotides was chosen to be outside the interaction region in the z-direction with the tilt of the nucleobase plane taken to be 30 degrees relative to z-x-plane. The single-layer graphene was thermally equilibrated at room temperature prior to sampling. The present simulation was performed in the absence of hydration in the system and hydrophobicity of graphene; therefore, the electrostatic interactions term is not included. For correlation evaluation of nucleotide spectra, the interaction time with graphene was taken at the interaction interval of the passing time through the pore. The high-resolution spectral calculations were carried out at the step size

of Δ*t* = 0.05 fs and *τ* = 32,768 steps of correlation delay for the shown states of nucleotides; in the present calculations, cytosine was used to test the system. To attribute obtained frequencies to a particular type of vibrations, such as stretching, bending, or torsion, the autocorrelation function of Equation (1) was calculated over the reaction coordinates that were taken by the projection of velocities on the bond vector, and the respective angle between these vectors as can be seen in Figure 2d. The single frequency obtained in such coordinates can be pinpointed as corresponding to the atom of the vibrating molecule, as well as to the type of vibrational mode that can be pinpointed as corresponding to the particular bond or angle. The ring-breathing modes can be identified by the same frequency present in the spectra of all bonds forming the ring. However, low-level noise from the nearest bonds to the bond studied was also present due to the correlation's calculation method. Therefore, only relatively high intensity spectral modes were considered for evaluation.

The stability of the obtained frequencies relative to the trajectory randomization has also been tested for cytosine molecules. The 10% to 15% randomization in the initial position and velocity distributions of the cytosine nucleotide at five different calculations have shown a preliminary absence of the changes in spectra above a few percent for stretching of the X-Y (X, Y = C, N, O) bonds. Our calculation system does not include solution molecules yet.

To create a calculation model that would be close to an experimental set up in a SERS type of measurement, we decided to introduce a metal nanoparticle in the nucleotide– graphene system of Figure 1c. The Au20 nanoparticle was pyramid-shaped and optimized in ab initio DFT GGA (density functional theory with generalized gradient approximation) calculations with the b3pw91 basis set. The size and shape of the NP was taken to be a stable configuration [50]. The obtained configuration was further relaxed by the MD calculation with the EAM potential [51] and free boundary conditions during 3000 steps with Δ*t* = 0.1 fs at *T* = 300 K to obtain stable nanoparticle samples that can be used separately or on the top of the graphene sheet. The potential parameters for LJ interaction with nucleotide and graphene were defined using the Lorentz–Berthelot mixing rule and based on the LJ parameters ε (Au −Au) = 0.039 kcal/mol, σ (Au −Au) = 2.9 Å [52]. The mixing was done using Van der Waals parameters of the MM2/MM3 force field potential for nucleotides. For the interaction with graphene, mixing used the sp3 C atom parameters from the force field potential, the same values that were applied for the nucleotide–graphene Van der Waals parameters.

The study of the vibrational spectra of the nucleotides in the dynamic interaction with metallic nanoparticles (NP) [53,54] located close to the graphene nanopore can combine translocation localization and nucleotide interaction enhancement or modification due to (1) the graphene LJ interaction and (2) the Au LJ interaction. For the initial interaction of the Au20 nanoparticle with cytosine nucleotide, the Au20 NP was localized close to the cytosine initial position as it is shown in Figure 3. As the first step, both NP and cytosine had zero c.o.m. velocity. Initial orientation and distance of NP relative to the nucleotide was initially controlled to estimate vibrational spectra of both interacting parts of the system.

**Figure 3.** Bonds in cytosine nucleotide (cytidine) for which vibrational spectra are shown below. Numbers of atoms are in internal numbering order.
