Electronic Couplings and On-site Energies for Hole Transfer in Dna: Systematic Quantum Mechanical/molecular Dynamic Study

The electron hole transfer ͑HT͒ properties of DNA are substantially affected by thermal fluctuations of the ␲ stack structure. Depending on the mutual position of neighboring nucleobases, electronic coupling V may change by several orders of magnitude. In the present paper, we report the results of systematic QM/molecular dynamic ͑MD͒ calculations of the electronic couplings and on-site energies for the hole transfer. Based on 15 ns MD trajectories for several DNA oligomers, we calculate the average coupling squares ͗V 2 ͘ and the energies of basepair triplets XG + Y and XA + Y, where X, Y = G, A, T, and C. For each of the 32 systems, 15 000 conformations separated by 1 ps are considered. The three-state generalized Mulliken–Hush method is used to derive electronic couplings for HT between neighboring basepairs. The adiabatic energies and dipole moment matrix elements are computed within the INDO/S method. We compare the rms values of V with the couplings estimated for the idealized B-DNA structure and show that in several important cases the couplings calculated for the idealized B-DNA structure are considerably underestimated. The rms values for intrastrand couplings G-G, A-A, G-A, and A-G are found to be similar, ϳ0.07 eV, while the interstrand couplings are quite different. The energies of hole states G + and A + in the stack depend on the nature of the neighboring pairs. The XG + Y are by 0.5 eV more stable than XA + Y. The thermal fluctuations of the DNA structure facilitate the HT process from guanine to adenine. The tabulated couplings and on-site energies can be used as reference parameters in theoretical and computational studies of HT processes in DNA.


INTRODUCTION
2][3][4] It has recently been shown that the electrical conduction of single molecules, in particular, DNA, is closely related to electron transfer rate, [5][6][7] and therefore, computational approaches designed to calculate absolute rates for charge transfer in DNA may be employed to estimate the electrical conductance of single DNA sequences.
The efficiency of charge transport is mainly determined by the electronic coupling V of donor and acceptor sites.The coupling values are usually computed for idealized molecular DNA geometries assuming the validity of the Condon approximation ͓V does not change with the molecular geometry in the charge transfer ͑CT͒ reaction͔.Several years ago, it was demonstrated that the electronic coupling for HT between nucleobases in DNA is very sensitive to structural fluctuations and the Condon approximation is rather limited. 8ombining molecular dynamics with an approximate scheme for estimating electronic couplings, Troisi and Orlandi showed that in DNA the standard deviation of the coupling of nucleobases is much larger than its average value. 9They also concluded that HT in DNA occurs in specific conformations, which may considerably deviate from the canonical B-DNA structure.1][12][13][14][15][16][17][18][19][20] The obtained results suggest that the electronic couplings found for idealized B-DNA ͑Refs.21-23͒ may substantially differ from the corresponding values averaged over thermally accessible DNA configurations, and therefore, a combined quantum mechanical/ molecular dynamic ͑QM/MD͒ approach should be used to obtain more accurate estimates of the electronic couplings.
Based on comparison to the CASPT2 and CASSCF results, 24 it has been demonstrated that the INDO/S method provides surprisingly good estimates for the hole transfer ͑HT͒ electronic couplings and energetics in DNA stacks. 25n particular, the INDO/S couplings are in better agreement with the reference data than the Hartree-Fock values.Because the INDO/S method is computationally very efficient, one can use this scheme to treat a huge number of structures extracted from MD trajectories.In the present work, we report the results of a systematic QM/MD study of electronic couplings and on-site energies for hole transfer in DNA.For each of the basepair triplets XGY and XAY, we carry out calculations of 15 000 conformations separated by 1 ps.Then, we compare the rms coupling values with the data estimated for the idealized B-DNA structure.The tabulated QM/MD couplings and on-site energies may be used as reference parameters in theoretical and computational studies of hole transfer phenomena in DNA.

COMPUTATIONAL DETAILS
Our quantum mechanical calculations are based on the 15 ns benchmark MD trajectories for several doublestranded DNA sequences ͑Table I͒ obtained within ABC project. 26The trajectory and parameter files were downloaded from the website of the Beveridge group. 27Each of considered 15-basepair oligomers contains three tetranucleotide repeats 5Ј-G D ͑A B C D͒ 3 G-3Ј.The basepair triplets XGY and XAY were extracted from the sequences as listed in Table I.To avoid potential artifacts from end effects, the triplets were taken from a middle part of the oligomers ͑the considered triplets are placed at least five basepairs from the ends of oligomers͒.
MD simulations.It was performed by ABC consortium 26 using well established protocol ͑T = 300 K, P = 1 atm, 2 fs integration step, parm94 force field, 28 TIP3P water molecules, 29 periodic boundary conditions, cutoff of 9 Å for nonbonded interactions, particle-mesh Ewald algorithm 30 for treatment of electrostatic interactions͒.Details of the simulation protocol are described in Ref. 26.
Quantum chemical calculations.As noted above, the INDO/S ͑Ref.31͒ method provides accurate values of electronic coupling and on-site energies. 25The scheme is computationally very efficient; it allows us to treat all 1 ps snapshots for XGY and XAY triplets extracted from 15 ns MD trajectories ͑15 000 different conformations were calculated for each of the 32 basepair triplets͒.The adiabatic energies and dipole moment matrix elements were computed with the INDO/S approach.We used the generalized Mulliken-Hush ͑GMH͒ method introduced by Cave and Newton 32,33 to derive electronic couplings and on-site energies.Within this scheme the diabatic states are chosen diagonal with respect to the component of the dipole moment operator along the main DNA axis.Once the unitary transformation T is defined, the electronic couplings ͑off-diagonal matrix elements of the diabatic Hamiltonian͒ can be calculated as V da = ͚T id E i T ia .The three-state formulation of the method was applied ͑details of calculations are considered elsewhere, Ref. 20͒.The use of the three-state GMH scheme for triplet XBY makes it possible to obtain at once the electronic couplings V͑X-B͒, V͑B-Y͒, and V͑X-Y͒ as well as the on-site energies E͑X + ͒, E͑B + ͒, and E͑Y + ͒ ͑diagonal elements of the diabatic Hamiltonian͒ corresponding to the states with the hole on basepairs X, B, and Y, respectively.Because the onsite energies are sensitive to the nature of neighboring basepairs, 34 below we consider the dependence of the B + state energies on the nature of X and Y neighbors.

Electronic couplings
As a measure of the coupling we use its root mean square V rms , V rms = ͱ͗V 2 ͘ = ͱ ͑1 / n͚͒ i=1 n V i 2 , where ͗¯͘ denotes the arithmetic mean, the sum is taken over n configurations.By definition, V rms 2 = ͗V͘ 2 + 2 , where is the standard deviation of V.The use of V rms rather than ͗V͘ is in line with the expression for the nonadiabatic CT. 35,36 For structurally flexible donor-acceptor systems the rate constant can be expressed as k =2/ បV rms 2 FC , if fluctuations in the coupling are slower than the evolution of the Franck-Condon factor FC . 37he intra-and interstrand couplings derived for different 32 basepair triplets are listed in Tables II and III.
Intrastrand couplings.Each of the intrastrand couplings V͑G-G͒, V͑A-A͒, V͑G-A͒, and V͑A-G͒ was extracted from the calculation of several triplets ͑Table II͒.Let us consider, for instance, the V͑G-G͒ values.In the oligomer 5Ј-GG ͑GGGG͒ 3 G-3Ј, the G-G couplings are expected to be ͑al-most͒ identical.However, the calculation of G 7 G 8 G 9 ͑indices i, j, and k in triplet X i B j Y k point at the basepair position in oligomers͒ gives V͑G 7 -G 8 ͒ and V͑G 8 -G 9 ͒ to be 0.0707 eV and 0.0627 eV, respectively.The difference ͑by ϳ10%͒ encountered in the calculation is due to the fact that G 7 and G 9 are nonequivalent even in the ideal B-DNA stacks because of the helical structure ͑they become equivalent in a completely unwound triplet structure͒.From calculations of the triplet G 8 G 9 G 10 we get V͑G 8 -G 9 ͒ = 0.0701 eV and V͑G 9 -G 10 ͒ It means that the difference of G-G couplings in the GGG triplets is solely caused by the end effects.In principle, these effects can be reduced by treating more extended QM models consisting, for instance, of four or five basepairs.However, such an approach will not only result in considerable computational cost but may also lead to serious difficulties by estimating the couplings using the multistate GMH scheme.As recommended by Cave and Newton, the number of states n treated simultaneously should be kept as small as possible. 32Taking into account that many thousands of calculations must be performed, a consistent multistate analysis with n Ͼ 3 is rather difficult.
Comparison of V͑G-G͒ values extracted from the calculation of GGY and XGG ͑Table II͒ suggests that the coupling depends on the nature of neighboring pairs X and Y.The coupling increases when X or Y changes as follows G Ͻ A Ͻ C Ͻ T. This effect may be associated with the two factors.First, alteration of neighboring pair leads to slightly different conformations of the G-G basepair step.However, this factor appears to be rather small. 26Second, the energy of hole state XG + G and GG + Y becomes higher when purine nucleobases ͑G and A͒ are replaced by pyrimidine nucleobases ͑C and T͒.
Similar behavior is also found for V͑A-A͒, V͑G-A͒, and V͑A-G͒ ͑Table II͒.In particular, the V rms value of A-A coupling ranges from 0.05 to 0.07 eV ͑Table II͒; in AAY and XAA the coupling increases by variation of X and Y in the following manner A Ͻ G Ͻ C Ͻ T.
Overall, V rms values for G-G, A-A, G-A, and A-G are found to be similar and the average value of 0.07 eV can be used to describe the intrastrand coupling of purine nucleobases in DNA -stack.Also, for other combinations of nucleobases ͑Table III͒, the coupling of B / BЈ is considerably stronger than that of B \ BЈ.For A/A and A\A, the average V rms values ͑0.064 and 0.017 eV, respectively͒ differ by a factor of 4, while V͑G / A͒ = 0.037 eV is three times stronger than V͑G\A͒ = 0.012 eV.Note that configurations G/A and A/G ͑as well as G\A and A\G͒ are equivalent.Because of a significant difference of V͑B \ BЈ͒ and V͑B / BЈ͒, the efficiency of hole transfer through a DNA stack with alternating purine-pyrimidine nucleobases will be determined by the product of the corresponding couplings; for instance, in poly͑GT͒-poly͑AC͒, the effective coupling squared is equal to ͉V͑G \ A͒ • V͑G / A͉͒.

Comparison to ideal B-DNA
2][23] In Table IV, we compare the values derived from the QM/MD computations with the couplings calculated for the idealized -stack geometries.As can be seen, in several cases the couplings calculated for ideal B-DNA structure significantly differ from the QM/MD data.More specifically, the intrastrand G-A and A-G cou-  plings as well as the interstrand G/G and A/A couplings in the basepair stacks of the ideal structures are in reasonable agreement with the QM/MD results; however, in all other cases the couplings derived for ideal B-DNA do not agree with the V rms values.Especially large deviations are found for V͑A-A͒, and V͑G / A͒.As already noted, HT couplings are extremely sensitive to conformational fluctuations.As an example, let us consider the data for coupling in the AAG trimer.As illustrated in Fig. 1, the coupling oscillates around zero frequently changing the sign.Its mean value of ϳ−0.015 eV is by a factor of 4 smaller than the standard deviation = 0.058 eV.Thus, the magnitude of V rms ͑A-A͒, which is very close to the , is almost completely determined by structural fluctuations.Such a situation is also observed for other couplings.Taking into account the extreme sensitivity of V to conformational changes, the acceptable estimates obtained for the idealized B-DNA structure appear to be really surprising.,39 On-site energies The energetics of hole states was already studied for DNA -stacks of the ideal geometry. 34,40It was shown that the energetic stabilization of a hole state B + in 5Ј-XBY-3Ј duplexes is considerably influenced by the subsequent base Y while the effect of the preceding base X is rather small. 34nlike electronic couplings, the energies of G + and A + hole states are found to be not very responsive to conformational fluctuations.Table V lists the mean values of on-site energies for hole states G + and A + in the triplets calculated for 15 ns trajectories for DNA oligomers ͑Table I͒.As expected, the most effective trap for cation radical states is guanine in triplets GGG and AGG.Replacing 5Ј-X by T or C leads to destabilization of XG + G by 0.1 eV.As noted above, the Y-3Ј nucleotide has an essentially larger effect on the G + energetics.As compared to 5Ј-XG + G-3Ј, the G + state in 5Ј-XG + A-3Ј are less stable by about ϳ0.10 eV ͑Table V͒, followed by 5Ј-XG + T-3Ј ͑0.27 eV͒ and 5Ј-XG + C-3Ј ͑0.30 eV͒.Thus, the stabilizing effect of 3Ј-Y neighbors decreases in the order G Ͼ A Ͼ T ജ C.
The calculated energies of XA + Y states are also listed in Table V.The mean difference between ionization energies of A and G in the triplets with identical neighbors is about 0.5 eV.The smallest difference of 0.34 eV is found for GA + A and GG + A. There is no overlap between the ionization potential of G and A, independent of neighboring nucleotides in XG + Y; the G + states of the highest energy ͑e.g., in CG + C͒ are found to be more stable by 0.1 eV than the lowest-energy states AA + G and GA + G ͑see Table V͒.
The QM/MD calculations show that the main results on the energetics of G + and A + hole states in DNA, previously established for ideal -stacks, 34 do not change significantly.However, structural fluctuations of the -stack lead to qualitatively new picture, facilitating the thermal hole transfer G + → A + within DNA.The thermal fluctuations cause variations of the diabatic energies of G + and A + .The corresponding standard deviations are almost independent of the neighboring basepairs and found to be of ϳ0.15 eV.If one takes into account also environmental fluctuations ͑movement of

CONCLUSIONS
In many theoretical studies on hole transport in DNA, the electronic couplings have been used which were derived for basepair stacks of idealized structure.However, these parameters are known to be very sensitive to conformational changes in the stack.To justify this situation, we carried out the systematic study of electronic coupling for hole transfer between nucleobases in DNA using a QM/MD approach.Based on 15 ns MD trajectories for several DNA oligomers, we calculate average coupling squares ͗V 2 ͘ and mean on-site energies in the basepair triplets XGY and XAY, where X, Y = G, A, T, and C.
We found that in several important cases the couplings estimated for the idealized B-DNA structure differ considerably from the QM/MD results, especially large deviations are found for the V͑A-A͒ and V͑G / A͒.Because those couplings are much smaller than the corresponding QM/MD values, the HT efficiency calculated ignoring thermal fluctuations can be underestimated by several orders of magnitude.
Unlike electronic couplings, the energies of the hole states G + and A + are not very sensitive to conformational changes.The effects of neighboring nucleobases on the HT energetics and numerical estimates, which were established for ideal -stacks, remain almost unchanged.While the states XG + Y are found by 0.5 eV more stable than the XA + Y states, structural fluctuations of DNA allow the hole to be transferred from G to A.
The tabulated QM/MD results for the electronic couplings and on-site energies can be used as reference parameters in theoretical and computational studies of electrical conductance and other phenomena related to hole transport through DNA.
Interstrand couplings.In poly-͑GC͒ double-strand sequence the hole transfer rate depends on the interstrand couplings of guanines ͑see Scheme 1͒.Even in the ideal -stack, the mutual position of G 1 -G 2 and G 2 -G 3 are quite different.The structural dissimilarity of G 1 G 2 and G 2 G 3 stacks is reflected in the coupling values G\G and G/G, respectively ͑hereafter to designate the mutual position of the nucleotides in opposite strands we use B \ BЈ for the configuration like G 1 G 2 or G 3 G 4 , and B / BЈ for the position G 3 G 2 or G 5 G 4 , see Scheme 1͒.Our calculations of V rms for G/G and G\G show that these couplings are quite different ͑by a factor of 4͒; the average values of V͑G / G͒ and V͑G\G͒ are found to be 0.0076 and 0.031 eV ͑Table III͒.

TABLE I .
15-mer oligomers used to extract geometries of trimers XGY and XAY.

TABLE III .
QM/MD interstrand couplings of neighboring pairs ͑eV͒ derived for different triplets.The calculated triplets are in parentheses.The bold G and A mark the basepairs for which the coupling is computed.

TABLE IV .
Comparison of the QM/MD couplings, with estimates derived for ideal B-DNA.
b Unsigned values are given ͑Ref.23͒.FIG. 1. Fluctuations of electronic coupling A-A calculated in the basepair trimer AAG.

TABLE V .
On-site energies for hole states XB + Y localized on G and A ͑eV͒.Then, the standard deviation of the free energy ⌬G for hole transfer G + → A + can be estimated as 2 Ϸ 2 ͑E͑G + ͒͒ + 2 ͑E͑A + ͒͒ Ϸ 2 ϫ 0.22 2 Ϸ 0.1 eV, 2 leading to ͑⌬G͒Ϸ0.3 eV.Because the estimated value of ͑⌬G͒ is comparable to the energy gap E͑A + ͒-E͑G + ͒ of ϳ0.5 eV, the thermally induced G + → A + hole transfer in DNA appears to be quite possible.