RT-TDDFT study of hole oscillations in B-DNA monomers and dimers

Abstract

We employ Real-Time Time-Dependent Density Functional Theory to study hole oscillations within a B-DNA monomer (one base pair) or dimer (two base pairs). Placing the hole initially at any of the bases which make up a base pair, results in THz oscillations, albeit of negligible amplitude. Placing the hole initially at any of the base pairs which make up a dimer is more interesting: For dimers made of identical monomers, we predict oscillations with frequencies in the range $f\approx$ 20–40 THz, with a maximum transfer percentage close to 1. For dimers made of different monomers, $f\approx$ 80–400 THz, but with very small or small maximum transfer percentage. We compare our results with those obtained recently via our Tight-Binding approaches and find that they are in good agreement.

Keywords:

• DNA
• charge transfer in biological systems
• electronic transfer in nanoscale materials and structures

AMS Subject Classifications:

• 87.14.gk 82.39.Jn 73.63.-b

Public Interest Statement

Understanding charge transfer (CT) in DNA is important for many biological and technological applications. CT depends on many extrinsic and intrinsic factors. The purpose of this article is to study the effect of one of the intrinsic factors, i.e. the base sequence. This work is a systematic ab initio study of hole (electron void) oscillations in small DNA segments. It discusses in detail many aspects of the phenomenon and compares its results with those of simpler, semi-empirical methods.

1. Introduction

Carriers can be inserted in DNA via electrodes or generated by UV irradiation or by oxidation and reduction. We had better discriminate between the words transport (implying the use of bias between electrodes), transfer (a carrier moves without bias) and migration (a transfer over rather long distances). Charge transport, transfer and migration in DNA are affected by many external factors such as chemical environment, interaction with the substrate or quality of contacts (Maciá, 2009), as well as by many intrinsic factors like the sequence of bases (Cuniberti, Maciá, Rodriguez, & Römer, 2007). Charge transfer in DNA is important for DNA damage and repair (Dandliker, Holmlin, & Barton, 1997; Giese, 2002,   2006; Rajski, Jackson, & Barton, 2000) and it may be used as an indicator for the separation between cancer tumor and healthy tissue (Shih, Cheng, Wells, Hsu, & Römer, 2011). Although (unbiased) charge transfer in DNA nearly vanishes after 10 to 20 nm (Lambropoulos et al., 2015,  2016; Simserides, 2014), DNA is still a promising component for (biased) charge transport in molecular electronics, e.g. as a short molecular wire (Wohlgamuth, McWilliams, & Slinker, 2013). It may also be exploited in nanotechnology for nanosensors (Huang, Li, He, Zhang, & Wang, 2016) and molecular wires (Genereux & Barton, 2010; Kawai & Majima, 2015).

The scope of this article is to study the effect of base sequence on charge transfer, in all possible monomers (i.e. within one base pair) and dimers (i.e. within two base pairs), using Real-Time Time-Dependent Density Functional Theory (RT-TDDFT) (Theilhaber, 1992; Yabana & Bertsch, 1996). The DNA sequence is one of the most important intrinsic factors that affect charge transfer, since it has been indicated that e.g. the emergence of strand breaks preferentially at guanines is connected to hole migration to it from other bases (Burrows & Muller, 1998). Furthermore, there is a possible relationship between information storage and physical properties of DNA, which are directly related to its electronic structure, as determined by the base sequence (Maciá, 2009). In previous articles, we studied small and long B-DNA segments (Lambropoulos, Kaklamanis, Georgiadis, & Simserides, 2014; Lambropoulos et al. 2015,  2016,  2016; Simserides, 2014) with Tight Binding (TB) approaches, specifically with a wire model and an extended ladder model. We determined the spatio-temporal evolution of an extra carrier (hole or electron) along a N base-pair DNA segment and showed that the carrier movement has frequency content mainly in the THz domain, i.e. depending on the case, in the NIR, MIR and FIR range (ISO 20473, 2007). This part of the electromagnetic (EM) spectrum is significant because it can be used to extract complementary to traditional spectroscopic information e.g. about hydrogen bond stretching, bond torsion in liquids and gases, low-frequency bond vibrations and so on, and because it is relatively non-invasive compared to higher-frequency regions of the EM spectrum (Yin, Ng, & Abbott, 2012).

In this article we study B-DNA monomers and dimers with RT-TDDFT, which is one of the few computationally viable techniques to model carrier dynamics in molecular systems (Liang, Chapman, & Li, 2011; Lopata & Govind, 2011; Marques, 2006; Ullrich, 2012). RT-TDDFT can simulate density fluctuations of the order of attoseconds to femtoseconds and can model the effects of an ultrafast, intense, short-pulse laser field, including the full (nonlinear) response of the electron density (Provorse & Isborn, 2016). RT-TDDFT also provides the full frequency dependence of properties of interest via Fourier transform from the time domain to the frequency domain. In addition, the determination of many-particle eigenstates is avoided, and the calculation is reduced to an independent particle problem carried out in real time (Takimoto, Vila, & Rehr, 2007). Recently, RT-TDDFT was applied (Partovi-Azar & Kaghazchi, 2017) to study the charge transfer between the oxidized or reduced a-sulfur monomers and the neighbouring neutral ones. Another recent work used RT-TDDFT to study electron transfer through oligo-p-phenylenevinylene (OPV) and carbon bridged OPV (COPV) (Grånäs, Kolesov, & Kaxiras, 2017).

The rest of this article is organized as follows. In Section 2 we sketch the basic formalism of RT-TDDFT. In Section 3 we give some computational details and explain our notation. In Section 4 we present our results. Finally, in Section 5 we state our conclusions.

2. RT-TDDFT formalism

Density Functional Theory (DFT) (Hohenberg & Kohn, 1964; Kohn & Sham, 1965) is an efficient method for treating ground state properties of many electron systems e.g. molecules or solids. Years after its development, it was extended (Runge & Gross, 1984) to time dependent systems (TDDFT). Specifically, the Time-Dependent Kohn-Sham (TDKS) equations with an effective potential energy ${\mathit{\upsilon }}_{}\text{KS}(\mathbf{r},t)$, uniquely described by the time-dependent charge density, $\mathit{\rho }(\mathbf{r},t)$, are, in atomic units,(1) $$\begin{array}{cc}& i\frac{\mathit{\partial }}{\mathit{\partial }t}{\mathrm{\Psi }}_j(\mathbf{r},t)=\left[-\frac{1}{2}{\mathrm{\nabla }}^2+{\mathit{\upsilon }}_{}\text{KS}(\mathbf{r},t)\right]{\mathrm{\Psi }}_j(\mathbf{r},t)\hfill \\ \hfill & =\left[-\frac{1}{2}{\mathrm{\nabla }}^2+{\mathit{\upsilon }}_{}\text{ext}(\mathbf{r},t)+{\mathit{\upsilon }}_{}\text{H}(\mathbf{r},t)+{\mathit{\upsilon }}_{}\text{xc}[\mathit{\rho }](\mathbf{r},t)\right]{\mathrm{\Psi }}_j(\mathbf{r},t).\hfill \end{array}$$(1)

The charge density is the sum over all occupied orbitals $j=1,2,\dots N_{occ}$, i.e.,(2) $$\begin{array}{c}\hfill \mathit{\rho }(\mathbf{r},t)=\sum _{j=1}^{N_{\text{occ}}}{|{\mathrm{\Psi }}_j(\mathbf{r},t)|}^2,\end{array}$$(2)

where, ${\mathit{\upsilon }}_{}\text{ext}(\mathbf{r},t)$ includes external fields and nuclear potentials and ${\mathit{\upsilon }}_{}\text{H}(\mathbf{r},t)$ is the Hartree potential energy. Exchange and correlation effects are included in ${\mathit{\upsilon }}_{}\text{xc}[\mathit{\rho }](\mathbf{r},t)$.

RT-TDDFT is based on a direct numerical integration of Equation (1). This differs from the traditional linear-response approach, which is not actually a time-resolved method but instead solves Equation (1) in the frequency domain for the excitation energies of a system subject to a small perturbation (Lopata & Govind, 2011). Within RT-TDDFT, we solve the TDKS equations and obtain the electron density at each time step. The electron density is then used for the calculation of the Hamiltonian in the next cycle of the self-consistent process.

3. Computational details and notations

With the notation YX we mean that the bases Y and X of two successive base pairs Y-Y${}_{}\text{compl}$ and X-X${}_{}\text{compl}$ are located at the same strand in the direction $5^{\prime }$–$3^{\prime }$. X${}_{}\text{compl}$ (Y${}_{}\text{compl}$) is the complementary base of X (Y). In other words, the notation YX means that Y-Y${}_{}\text{compl}$ is the one base pair and X-X${}_{}\text{compl}$ is the other base pair, separated by 3.4 Å and twisted clockwise by 36${}^{}\text{o}$ relatively to the first base pair, along the growth axis of the nucleotide chain. For example, the notation AC denotes that one strand contains A and C in the direction $5^{\prime }$–$3^{\prime }$ and the complementary strand contains T and G in the direction $3^{\prime }$–$5^{\prime }$. The notation (10) means that the hole is initially placed at the 1st base pair of the dimer and the notation (01) means that it is initially placed at the 2nd base pair of the dimer.

For our DFT and RT-TDDFT calculations we used the NWChem open-source computational package (Valiev et al., 2010). Additionally, we performed Fourier analysis, implemented with MATLAB, to determine the frequency content of the charge oscillations. The geometry of monomers and all dimers was produced via the BIOVIA Discovery Studio (http://accelrys.com/resource-center/downloads/freeware/index.html). The backbone was removed and H atoms were added at the corresponding atomic sites.

The range-separated functional CAM-B3LYP (Yanai, Tew, & Handy, 2004), which is appropriate for the correct estimation of the exchange energy, both at short and long ranges, was used for all monomers and dimers. The calculations were performed using the 6-31++G** (Hehre, Ditchfield, & Pople, 1972), 3-21++G (Binkley, Pople, & Hehre, 1980) and aug-cc-pVDZ (Dunning, 1989) basis sets, which include diffuse functions, for all systems. Within a small number of time steps, we noticed that our results were similar (cf. Appendix 1). Hence, we chose to perform our longer simulations with the 6-31++G** basis set (Hehre et al., 1972), i.e. the least computationally expensive between the converging sets.

In a Gaussian basis set, it is most natural to use the single particle reduced density matrix, whose time evolution is governed by the von Neumann equation. The Magnus propagator is used in NWChem’s RT-TDDFT implementation, which is both stable and conserves the density matrix idempotency (Lopata & Govind, 2011). At the end of each time step the fragments’ charge is calculated with an appropriate population analysis method, along with the total dipole moment.

For monomers, the initial state (i.e. hole localized at a specific base) was produced by CDFT, where the charge constraint was calculated with Löwdin population analysis that is also used in the subsequent RT-TDDFT simulation. For dimers, the initial state (i.e. charge localization) was produced by the following procedure: first, a ground state DFT calculation of the charged and neutral isolated monomers was performed, yielding the corresponding monomers’ eigenstates. Then, the dimer’s eigenstates were approximated in a non self-consistent manner (NOSCF) by the orthogonalized monomers’ eigenstates. Since the monomers comprising each dimer are 3.4 Å apart, this is a fairly good approximation that helps circumvent CDFT convergence problems.

Löwdin (1950) population analysis was integrated into RT-TDDFT module of NWChem for the calculation of each fragment’s charge at each time step. Löwdin population analysis is known to be much less basis-set dependent and also it does not suffer from the ultra-fast charge oscillations that Mulliken analysis (which is the default scheme in NWChem’s RT-TDDFT) artificially introduces in RT-TDDFT charge simulations. As a result, Löwdin population analysis gives a more clear picture of charge oscillations convergence towards the basis-set limit, cf. Figures 7 and 8. The main frequencies of charge oscillations are extracted from the results via Fourier analysis. In order to increase the resolution, zero-padding with appropriate signal attenuation was used where necessary.

4. Results

Before discussing our RT-TDDFT results, we give a brief summary of the TB picture. In TB, charge transfer between the two bases of a monomer or between the two base pairs of a dimer can be easily analytically solved (Lambropoulos et al., 2014,   2016; Simserides, 2014). In both cases, it is a movement between the two sites $i,j$ of a two-site system and it is mathematically similar to Rabi oscillations in Quantum Optics (Simserides, 2016). The maximum transfer percentage p, i.e. the maximum probability to find the carrier at the site where it was not initially placed is(3) $$\begin{array}{c}\hfill p=\frac{{(2t_{ij})}^2}{{(2t_{ij})}^2+{({\mathrm{\Delta }}_{ij})}^2},\end{array}$$(3)

where ${\mathrm{\Delta }}_{ij}$ is the difference between the two sites’ on-site energies and $t_{ij}$ is the hopping parameter between the two sites. Since between the two bases of a monomer, ${\mathrm{\Delta }}_{ij}$ is much larger than $t_{ij}$, p is negligible. For movements between the two monomers of a dimer things are different: For dimers made of identical monomers, since ${\mathrm{\Delta }}_{ij}=0$, it follows that $p=1$. For dimers made of different monomers $p<1$; in fact, it is usually negligible with one exception, a hole moving between the two monomers of a GA dimer, where p is of the order of 0.25. In TB, the frequency and period of charge oscillations is given by(4) $$\begin{array}{c}\hfill f=\frac{1}{T}=\frac{\sqrt{{(2t_{ij})}^2+{\mathrm{\Delta }}_{ij}^2}}{h},\end{array}$$(4)

hence, increasing $t_{ij}$ and / or ${\mathrm{\Delta }}_{ij}$ results in higher frequencies. In what follows, this simple picture is qualitatively obeyed by our RT-TDDFT results.

Hole oscillations between the two bases of a base pair are in the THz regime, but have negligible maximum transfer percentage, in agreement with the TB results (Lambropoulos et al., 2016).

Hole oscillations between the two base pairs of a dimer, with the 6-31++G** basis set, are shown in Figure 1. We place the hole initially at one of the monomers which make up the dimer and call p the maximum transfer percentage to the other monomer. For dimers made of identical monomers (AT, TA, GC, CG), the hole is transferred almost completely to the other monomer ($p\approx 0.9$). For dimers made of different monomers (AC, CA, GA, AG), p is close to zero with the exception of the GA dimer, where $p\approx 0.3$. Our calculations for the cases AT(01), TA(01), GC(01), and CG(01) are not included in Figure 1, since they are identical with the cases AT(10), TA(10), GC(10), CG(10), respectively, as expected from base-pair symmetry and found also in Refs. Simserides (2014), Lambropoulos et al. (2014), Lambropoulos et al. (2016).

Figure 1. Net charge versus time for the monomer where the hole was initially placed: (upper panel) dimers made of identical monomers, (lower panel) dimers made of different monomers. The notation (10) means that the hole is initially placed at the 1st base pair of the dimer, while the notation (01) means that the hole is initially placed at the 2nd base pair of the dimer. The 6-31++G** basis set was used.

In Figure 2 we present the mean (over time) probabilities to find the extra hole at each monomer, having placed it initially either at the 1st monomer (10) or at the 2nd monomer (01). Let us call ${\mathcal{P}}_i$ the mean probability at the monomer where the hole is initially placed, and ${\mathcal{P}}_f$ the mean probability at the other monomer. We observe that for dimers made of identical monomers, where $p\approx 1$, we have ${\mathcal{P}}_i\approx {\mathcal{P}}_f\approx 0.5$, while for dimers made of different monomers ${\mathcal{P}}_i\approx$ 1 and ${\mathcal{P}}_f\approx 0$, with the exception of the GA dimer. This, again, agrees qualitatively with the TB picture (Simserides, 2014; Lambropoulos et al., 2014; Lambropoulos et al., 2016).

Figure 2. The mean (over time) probabilities to find the hole at each base pair of a dimer, having placed it initially at the 1st (10) or at the 2nd (01) base pair.

In Figure 3 we depict the maximum transfer percentage p as well as the electronic coupling energy (also called the electron transfer matrix element) $|V_{\text{RP}}|$ between the reactant and product states R and P (Bolton, Mataga, & McLendon, 1991; Farazdel, Dupuis, & Clementi, 1990; Newton & Sutin, 1984; Marcus & Sutin, 1985), which reflects the magnitude of the interaction between the two monomers. In our case, the reactant state corresponds to the hole at the initial monomer and the product state corresponds to the hole at the other monomer. We also show the energy difference $|\mathit{\delta }E|$ between reactant and product states. In analogy with a Rabi oscillation, we expect charge transfer to increase by increasing $|V_{\text{RP}}|$ and/or by decreasing $|\mathit{\delta }E|$. Additionally, for $|\mathit{\delta }E|\approx 0$, we expect p to be close to 1. For the AT, TA, GC and CG dimers, $|\mathit{\delta }E|\approx 0$, hence p is indeed close to 1. For dimers made of different monomers, we observe that AC, CA and AG have large $|\mathit{\delta }E|\approx$ 0.45–0.75 eV and small $|V_{\text{RP}}|\approx 0-1$ eV, therefore their p is insignificant. An exception is the dimer GA which not only has a large $|\mathit{\delta }E|\approx 0.65$ eV, but also the largest $|V_{\text{RP}}|$ of all dimers $\approx$ 5 eV which makes charge transfer significant, with p $\approx 0.3$.

Figure 3. Maximum transfer percentage, i.e. oscillation amplitude p, electron transfer coupling energy $|V_{\text{RP}}|$ and energy difference $|\mathit{\delta }E|$ between reactant and product states. Initially, the hole is placed at the 1st monomer (10) or at the 2nd monomer (01). Note: Dotted lines serve as guides to the eye.

Figure 4. Frequencies of hole oscillations in dimers made of identical monomers (upper panel) and different monomers (lower panel), obtained via Fourier analysis. Note: Only frequencies with non negligible amplitudes are shown.

In Figure 4 we present the frequencies of hole oscillations obtained by Fourier analysis. For dimers made of identical monomers (upper panel) the frequencies $f\approx$ 20-40 THz and their amplitudes $\approx$ 0.2–0.5. For dimers made of different monomers (lower panel) $f\approx$ 80–400 THz and the amplitudes $\approx$ 0.003–0.02, except for the two GA cases which have amplitudes $\approx$ 0.1–0.2. Obviously, larger amplitudes reflect larger transfer. As expected by TB, for dimers made of different monomers we generally get higher frequencies than for dimers made of identical monomers. Additionally, in Figure 5 we observe that for dimers made of identical monomers the frequencies follow $|V_{\text{RP}}|$, which does not hold for dimers made of different monomers. This, again, agrees with the TB prediction (Equation (4)).

Figure 5. Hole oscillation frequency f vs. $|V_{\text{RP}}|$ for all dimers. Note: Dotted lines serve as guides to the eye.

Finally, in Figure 6 we compare our RT-TDDFT results with those of TB (Simserides, 2014; Lambropoulos et al., 2016) for two different sets of TB parameters: (a) used in Ref. Lambropoulos et al. (2016) (calculated in Ref. Hawke, Kalosakas, & Simserides, 2010) and (b) used in Ref. Simserides (2014). We observe that the maximum transfer percentages obtained by RT-TDDFT are in good agreement with those obtained by TB. For the periods, the results, both for RT-TDDFT and TB are quite close especially when we compare them with parametrization (b) used in Ref. Simserides (2014).

Figure 6. Comparison of RT-TDDFT (this work) and TB from (a) Ref. Lambropoulos et al. (2016) and (b) Ref. Simserides (2014): maximum transfer percentage (left panel) and oscillation period (right panel). Note: Dotted lines serve as guides to the eye.

5. Conclusion

Understanding charge transfer in DNA is important for many biological and technological applications and depends on many extrinsic and intrinsic factors. The purpose of this article was to study the effect of the base sequence which is one of the intrinsic factors. Thus, we studied how base sequence affects hole transfer in small DNA segments (between the two bases of a monomer and between the two base pairs of a dimer) with the RT-TDDFT method.

Specifically, we calculated the maximum transfer percentage, i.e. the oscillation amplitude, p, the mean (over time) probability at the site where the hole is initially placed, ${\mathcal{P}}_i$, and the mean probability at the other site, ${\mathcal{P}}_f$, as well as the frequency content of the oscillations. For our calculations we used the CAM-B3LYP functional with three different basis sets (6-31++G**, 3-21++G and aug-cc-pVDZ) and we observed that our results were similar.

Hole oscillations between the two bases of a base pair are in the THz regime, but have negligible maximum transfer percentage.

Next, we studied hole oscillations between the two monomers (base pairs) of a dimer. For dimers made of identical monomers, the maximum transfer percentage is almost 1. For dimers made of different monomers, we observed negligible hole transfer except for the GA dimer where $p\approx 0.3$. We also calculated the electronic coupling energy, $|V_{\text{RP}}|$, between the two states and the relevant energy difference, $|\mathit{\delta }E|$. We saw that when $|V_{\text{RP}}|$ is large and $|\mathit{\delta }E|$ is small, the maximum transfer percentage is large, while for small $|V_{\text{RP}}|$ and large $|\mathit{\delta }E|$, the maximum transfer percentage is insignificant.

For all dimers, we found frequency content in the THz domain and noticed that Fourier amplitudes reflect p and ${\mathcal{P}}_f$. Specifically, for dimers made of identical monomers we found $f\approx$ 20–40 THz and Fourier amplitude $\approx$ 0.2–0.5, while for dimers made of different monomers $f\approx$ 80–400 THz and the amplitudes $\approx$ 0.003–0.02. An exception is GA with amplitude $\approx$ 0.1–0.2. We, also, observed that for dimers made of different monomers we generally get higher frequencies than for dimers made of identical monomers.

We compared the RT-TDDFT results with those obtained by TB in Refs. Simserides (2014), Lambropoulos et al. (2014), Lambropoulos et al. (2016) and we found that they are in good agreement.

Figure 7. Time evolution of hole transfer in dimers made of identical monomers, for the three different basis sets.

Figure 8. Time evolution of hole transfer in dimers made of different monomers, for the three different basis sets.

6. Related work at

http://users.uoa.gr/~csimseri/physics_of_nanostructures_and_biomaterials.html.

7. Cover image

Source: Authors.

Acknowledgements

This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility - ARIS - under [project ID pr002008 - CODNA].

Additional information

Funding

A. Morphis thanks the State Scholarships Foundation-IKY for a PhD research scholarship via ’IKY Fellowships of Excellence’, Hellenic Republic-Siemens Settlement Agreement.

Notes on contributors

M. Tassi

M. Tassi is currently a postdoctoral researcher at the National and Kapodistrian University of Athens (NKUA), Physics Department. Her main research interest is in the area of ab initio calculations studying electronic properties such as charge transfer in molecular systems, e.g. in DNA segments.

A. Morphis

A. Morphis is a PhD student at NKUA, Physics Department. His main research interest is in the study of electronic and optical properties of organic molecules and nanomaterials.

K. Lambropoulos

K. Lambropoulos is a PhD student at NKUA, Physics Department. Currently, he is studying the energy structure and the charge transfer properties of periodic crystalline, quasi-crystalline, fractal, amorphous, random and natural DNA segments.

C. Simserides

C. Simserides is an assistant professor at NKUA, Physics Department. He is mainly interested in charge transfer and transport along DNA, in optical properties of quantum dots and in ferromagnetic properties of diluted magnetic semiconductors.

Appendix 1

In Figures 7 and 8 we show charge in dimers oscillations calculated with different basis sets: 6-31++G** (Hehre et al., 1972), 3-21++G (Binkley et al., 1980) and aug-cc-pVDZ (Dunning, 1989). We observe that the results are similar (cf. Section 3).