One-Electron Oxidation of Individual DNA Bases and DNA Base

Jan 6, 2010 - Amitava Adhikary , Anil Kumar , Brian J. Palmer , Andrew D. Todd , and ... Amitava Adhikary , David Becker , Brian J. Palmer , Alicia N...
0 downloads 0 Views 3MB Size
Download PDF
1860

J. Phys. Chem. A 2010, 114, 1860–1867

One-Electron Oxidation of Individual DNA Bases and DNA Base Stacks David M. Close* Department of Physics, Box 70652, East Tennessee State UniVersity, Johnson City, Tennessee 37614 ReceiVed: July 22, 2009; ReVised Manuscript ReceiVed: NoVember 25, 2009

In calculations performed with DFT there is a tendency of the purine cation to be delocalized over several bases in the stack. Attempts have been made to see if methods other than DFT can be used to calculate localized cations in stacks of purines, and to relate the calculated hyperfine couplings with known experimental results. To calculate reliable hyperfine couplings it is necessary to have an adequate description of spin polarization which means that electron correlation must be treated properly. UMP2 theory has been shown to be unreliable in estimating spin densities due to overestimates of the doubles correction. Therefore attempts have been made to use quadratic configuration interaction (UQCISD) methods to treat electron correlation. Calculations on the individual DNA bases are presented to show that with UQCISD methods it is possible to calculate hyperfine couplings in good agreement with the experimental results. However these UQCISD calculations are far more time-consuming than DFT calculations. Calculations are then extended to two stacked guanine bases. Preliminary calculations with UMP2 or UQCISD theory on two stacked guanines lead to a cation localized on a single guanine base. Introduction There is much discussion in the literature on the nature of the one-electron oxidation products in purine bases. Guanine, which has the lowest ionization potential of the DNA nucleobases is therefore the most easily oxidized. An electron loss center created in a DNA duplex by one-electron oxidation ultimately ends up on a guanine base via hole migration through the DNA stack. Interest now has shifted to the study of runs of guanines in duplex DNA as these have been shown to be the more easily oxidized than isolated guanine sites. Experimental studies by Sugiyama and Saito showed that the 5′-G of the 5′-GG-3′ sequence is the most easily oxidized site in the B-form of DNA.1 They followed this work with a theoretical study of stacked DNA bases, which showed that the localization of the guanine cation (hole) depends on the stacking geometries. For small angles of twist between the stacked guanines, the highest occupied molecular orbital (HOMO) of the guanine cation is distributed equally on the two bases. For the normal ranges of twist found in B-DNA (-25° to -45°) 95% of the HOMO is localized on the 5′-G.2 The original series of articles by Sugiyama and Saito appeared in 1995-96. Since then numerous studies on stacked purines have been reported. Among these is an article by Barnett et al. on the four base pair duplex 5′-G1A2 G3 G4-3′ showing the hole being shared between the G3G4 pair.3 This study also considered the influence of Na+ counterions on hole localization. In a study of the AGGA/TCCT duplex, simulation by Bongiorno showed the hole delocalize over the two guanines.4 Parrinello and coworkers have shown that problems can arise in characterizing the electronic structure of the guanine cation which is equally shared between two stacked guanines in their B-DNA conformation in density functional theory (DFT), but is isolated on one of the guanines with ROHF theory, or with a self-interaction correction procedure to be discussed below.5-7 In studies of stacked adenine Adhikary et al. see similar delocalization of * To whom correspondence should be addressed. E-mail: Closed@ etsu.edu.

the adenine cation onto both bases in a stack of two adenines.8 Improta has studied pentamer stacks of 9-methyladenine in the B-DNA configuration and see the HOMO of the adenine cation delocalized over the three central adenines.9 It seems then that there are numerous studies of stacked purine bases that show the purine cation delocalized over several purines. At one time, calculations on two purines in a stack were very time-consuming. With the advent of density functional theory these open-shell calculations are tractable. However, density functional theory has a tendency to yield results that delocalize the purine cation. The construction of most KohnSham DFT functionals lead to a nonzero interaction of a single electron with its own density.10 This self-interaction of the electron with itself is called the self-interaction error. Thus, it is important here to look at the results obtained so far and to ask if some of the results presented in the literature arise from the self-interaction error. To address the problem of possible artifacts, it would be important to compare the calculated results with experimental results. This has been done unambiguously by looking at spin densities and the resulting hyperfine couplings.11,12 The hyperfine couplings for various nucleobase oxidation products have been accurately determined from detailed EPR/ENDOR experiments.13 It should be noted that in most cases these experimental results represent deprotonated cations. This point is considered below in the section on stacked guanines. Now it needs to be determined if theoretical calculations of hyperfine couplings are reliable enough to be of use to distinguish between localized and delocalized holes in stacked purines. For a single nucleic acid base there is no problem with using DFT to compute accurate hyperfine couplings. Results have been reported in a series of papers by Wetmore, Boyd, and Eriksson.14-17 A few results in these original studies contained nonplanar geometries of the bases. These commonly arise for optimizations on a single base without considerations of the H-bonding environment that the free radical is imbedded in. With the proper inclusion of these neighboring interactions, and with the inclusion of a few waters of hydration, very good

10.1021/jp906963f  2010 American Chemical Society Published on Web 01/06/2010

One-Electron Oxidation of DNA Bases and Stacks

J. Phys. Chem. A, Vol. 114, No. 4, 2010 1861

TABLE 1: Calculated Hyperfine Couplings (MHz) in the Allyl Radical atom 13

C1 C2 C1-H1 C1-H2 C2-H 13

a

UHF

CIS

CID

CISD

QCISD

MP2

MP4

B3LYP

exp.a

175.80 -158.04 -103.29 -99.37 66.11

38.30 187.13 -5.00 -5.45 -126.94

95.64 -67.99 -60.88 -57.66 18.94

101.01 -74.43 -63.83 -60.55 22.19

89.55 -60.69 -55.50 -52.51 13.83

71.37 -35.82 -50.58 -47.27 2.33

78.81 -48.82 -50.57 -47.40 6.74

74.26 -45.84 -45.03 -42.78 13.15

61.56 -48.23 -41.56 -39.01 11.77

Experimental results in ref 28.

hyperfine couplings are easily calculated. Results have been summarized in a review article.18 For the present study, it is important to know how to calculate the hyperfine couplings of purines in a stack. If DFT calculations on stacks tend to result in delocalization of the purine cation, then it would be interesting to see if post Hartree-Fock methods give the same description of the wave function. Molecular orbital calculations on molecules with open shell electronic configurations may be performed by using restricted Hartree-Fock (RHF) or unrestricted Hartree-Fock (UHF) procedures. In RHF theory, the spatial components of the R and β orbitals are constrained to be the same. As such, the resulting wave function yield pure electron spin states that are eigenfunctions of the S2 operator. However, the resulting wave functions cannot lead to negative spin densities and cannot therefore describe spin polarization which is important in calculating hyperfine couplings. In the UHF procedure, the R and β exchange interactions are evaluated separately and can yield a more realistic description of the spin distribution in open-shell systems. However, the UHF wave functions are not eigenfunctions of the S2 operator, since they are contaminated by states of higher spin multiplicity. Attempts to improve the HF wave function involve adding terms that represent promotion of electrons from occupied to virtual molecular orbitals. The configuration interaction (CI) treatment of electron correlation may involve single excitations (CIS), double excitations (CID), a combination of single and double excitations (CISD), and the quadratic configuration interaction (QCISD).19 Another approach to include electron correlation involves Møller-Plesset perturbation theory (MPn). The second order MP2 method requires the inclusion of double excitations with respect to the HF reference function in the calculations. At fourth order MP4, single, triple, and quadratic excitations are considered, in addition to double excitations. A cluster expansion of the UHF wave function restricted to the double replacement operator leads to the coupled-cluster UCCD method. The additional inclusion of the single replacement operator in the cluster expansion gives the UCCSD method. It has been shown that the UCCSD method may be satisfactorily approximated by the UQCISD method.19 Carmichael20 has shown that (QCISD) calculations are particularly good at calculating reliable hyperfine couplings. Computational Methods. Molecules were first optimized using second order Møller-Plesset perturbation theory (MP2) and hybrid Hartree-Fock density functional theory (B3LYP) in conjunction with the 6-31+G(d) basis set. For the test case of the allyl radical optimizations were also carried out with quadratic configuration interaction (QCISD) theory. For the calculations on the nucleic acid bases only single point QCISD calculations were performed on the optimized structures. All calculations were performed on the Gaussian 98/03 suite of programs.21,22 To recover the proper spin densities it is necessary to use the keywords Density)MP2, or Density)QCI.

To visualize the spin densities it is necessary to use the keywords spin)MP2 or spin)cc in the GaussView program.23 All spin density results presented here are shown with the isodensity value of 0.0004 e(Å)-3. In performing calculations on purine stacks, it became clear that small changes in starting geometries made big changes in the outcomes. To explore different starting geometries, tests were performed using PQS software24 that has the Guess)Atom option.25 This option involves an initial guess based on the sum of atomic densities rather than a Hu¨ckel molecular orbit guess. Results and Discussion Test Calculations on the Allyl Radical. To test the various techniques available, a number of calculations were performed on the allyl radical. The allyl radical (C3H5) is a prototype of conjugated π-radicals.26 It is the simplest example of an odd alternate radical and as such shows sign alternations of the π spin density, with negative spin density at the central carbon. For an excellent review, see Bally and Borden’s book chapter entitled, “Calculations on Open-Shell Molecules”, which uses calculations on the allyl radical in numerous examples.27 For consistency, all calculations were done with the same basis set (6-31+G(d)). Some will argue that this basis set is too small for calculating accurate hyperfine couplings (especially for MPn theory). While this is true, the purpose here is to explore methods that can be scaled up to calculations on stacked purines. Doing so with large basis sets would be very time-consuming. The results of hyperfine coupling constant calculations for the allyl radical are shown in Table 1. The results in Table 1 show the MPn results are in fair agreement with the experimental hyperfine couplings,28 except for the small proton coupling at the central carbon (C2-H). The results indicate that the (QCISD) method has the ability to predict hyperfine couplings that are near to the experimental values for all three proton couplings. For comparison, DFT calculations are also included that are seen to give the best agreement with the experimental results. For the present discussion, it is important to see how visualization programs display spin densities. From Table 1, spin densities are shown (Figure 1) which were derived from the outputs of the MP2, QCISD, and DFT calculations. On the left of Figure 1 the structure of the allyl radical is shown. For the spin densities, the molecule is rotated 90° (shown edgewise) so one can view the spin densities perpendicular to the radical plane. One sees in the visualizations of the spin densities that the MP2 and QCISD calculations are very different. This then nicely shows the differences in the hyperfine couplings at C2-H (Table 1). The visualizations of the spin densities at the QCISD and DFT levels of theory are very similar, also reflected in the computed hyperfine couplings in Table 1. Calculations of the DNA Bases. The goal here is to find a computational approach that gives hyperfine couplings that agree well with the experimental results, and that would be useful

1862

J. Phys. Chem. A, Vol. 114, No. 4, 2010

Close

Figure 1. Spin densities in the allyl radical. Positive spin densities are in blue and negative spin densities are in green. All three calculations were performed with the same 6-31+G(d) basis set.

TABLE 2: Calculated Hyperfine Couplings (MHz) for the 1-MeThymine Cation atom

MP2

QCISD

B3LYP

exp.a

N1 N3 N1-CH3 N3-H C5-CH3 C6-H

-8.48 2.96 10.8 -8.64 63.88 -44.63

24.80 -3.54 24.06 0.16 57.41 -3.40

19.25 -2.57 23.13 -0.11 61.1 -1.39

17.6 2.0 23.0 neglig. 57.9 -2.0

a Experimental results from ref 30. Neglig. indicates a coupling too small to measure.

and practical to use with stacks of bases. Options including RHF, UHF, MP2, hybrid DFT, and various CIS, CISD, QCISD methods are explored. As a first step, it was decided to attempt MP2 and QCISD calculations on single DNA bases. While the QCISD method has been available for some time, it is presently not used very much. When it was being developed (1987), QCISD calculations on molecules with 3-4 heavy atoms were very time-consuming.29 Since most calculations are now done with DFT, there is not much available in the literature on using the QCISD method on molecules the size of the nucleic acid bases. Here MP2 and QCISD calculations are performed on the nucleic acid bases with the goal of computing accurate hyperfine couplings of the one-electron oxidization products. The results can then be compared with experimentally measured hyperfine couplings from detailed EPR/ENDOR experiments.13 1-Me-Thymine Cation. The EPR data on the 1-MeThymine cation indicates a π-radical structure with the spin density shared between N1 and C5 (the numbering scheme used here is shown in the Supporting Information as S1). The dominant hyperfine couplings determined experimentally for the native thymine cation are N1(17.6 MHz), N1-CH3(23.0 MHz) and C5CH3(57.9 MHz).30 Problems are often encountered in doing open-shell calculations with different methods on the pyrimidines. In some cases the optimizations at the MP2 level of theory converge on a structure with an elongated N1-C2 bond.31 The tendency to converge to this ring-opened (σ-radical) structures can be mitigated by studying the N1 methylated pyrimidine (>N1-CH3). Table 2 shows the results of MP2, QCISD, and DFT hyperfine calculations on the 1-MeThymine cation. Here one sees problems with the calculations at the MP2 level of theory. The hyperfine coupling at N1 is of the wrong sign, and the hyperfine coupling to N3-H and C6-H is overestimated. In Table 2, the experimental N3-H hyperfine coupling is listed as negligible indicating it is too small to measure (less than 2 MHz) and only contributes to the overall EPR line width. In comparing the hyperfine couplings from the MP2 calculation with calculations of the hyperfine couplings at the QCISD level of theory, it seems as if the QCISD calculations are in better

agreement with the experimental values. For comparison, the hyperfine couplings computed with DFT are included in the Table 2. The results with QCISD and DFT theory have the proper spin densities at N1 and at C5 and are in much better agreement with the experimental results in the smaller hyperfine couplings. The features described in Table 2 are illustrated in Figure 2. In particular, spin densities calculated at the MP2 level of theory show negative spin density at N3-H and at C6-H as shown in green. The negative spin density at N1 is obscured by the positive π-orbital (see Supporting Information S5). The spin densities are much different at the QCISD level of theory. Now there is no negative spin density at N1 or at C6-H. Furthermore, the general shape of the spin densities is quite similar to the spin densities calculated with DFT theory. The calculations at the QCISD and DFT levels of theory are seen to be in good agreement with the spin densities determined experimentally. 1-Me-Cytosine Cation. The EPR data on the 1-MeCytosine cation indicates a structure with the spin density shared between N1 and C5 (the numbering scheme used here is shown in the Supporting Information as S2). The dominant hyperfine couplings determined experimentally for the native cytosine cation are N1 (17.6 MHz) and C5-H (-41.2 MHz).32 In a single crystal study of 5′dCMP an additional isotropic hyperfine coupling to C1′-H was measured to be 41.9 MHz.32 In the 1-MeCytosine model used here this additional hyperfine coupling is equivalent to one large N1-CH3 hyperfine coupling. Table 3 shows the results of MP2, QCISD, and DFT calculations on the 1-MeCytosine cation. One sees that the hyperfine couplings calculated at the MP2 level of theory are in poor agreement with the experimental results. The N1, N3, and C6-H hyperfine couplings are of the wrong sign, and the N3 and C6-H couplings are too large. Alternatively, the hyperfine couplings calculated at the QCISD level of theory are in rather good agreement with the experimental hyperfine couplings. For comparisons the hyperfine couplings computed with DFT are included in the Table 3. The results with both QCISD and DFT theory have the proper spin densities at N1 and at C5, and are in much better agreement with the experimental results in the smaller hyperfine couplings. The features described in Table 3 are illustrated in Figure 3. Here one sees the negative spin density on N1 computed with MP2 theory, and the correct large positive spin density at N1 predicted by QCISD and DFT theory. Likewise, the large negative spin density at C6-H predicted by MP2 theory can be compared with the much smaller values predicted by QCISD and DFT theory. The negative spin density at C5-H in these figures can only be seen by rotating these figures to observe the radical plane (see Supporting Information, S6). Adenine Cation. The EPR data on the native adenine cation indicates a structure with the spin density shared between C8 and the exocyclic nitrogen N6 (the numbering scheme used here

One-Electron Oxidation of DNA Bases and Stacks

J. Phys. Chem. A, Vol. 114, No. 4, 2010 1863

Figure 2. Computed spin densities of the 1-methyl-thymine cation. All three calculations were performed with the same 6-31+G(d) basis set.

TABLE 3: Calculated Hyperfine Couplings (MHz) for the 1-MeCytosine Cation a

atom

MP2

QCISD

B3LYP

exp.

N1 N3 N4 N1-CH3 N4-H1 N4-H2 C5-H C6-H

-7.78 -36.49 7.75 17.20 -14.96 -13.27 -42.80 -25.69

25.20 7.49 -1.30 22.57 -0.16 1.08 -49.15 6.05

17.20 5.52 -0.64 26.64 -0.59 -0.36 -37.06 4.36

17.6 2.0 9.2 41.9 neglig. neglig. -41.2 2.0

a Experimental results from ref 32. Neglig. indicates a coupling too small to measure.

is shown in the Supporting Information as S3). The dominant hyperfine couplings are C8-H (-13.7 MHz), N6-H1(19.7 MHz), and N6-H2(-20.5 MHz).33 Table 4 shows the results of MP2, QCISD, and DFT calculations on the adenine cation. One sees that the hyperfine couplings calculated at the MP2 level of theory are in poor agreement with the experimental results. The C8-H hyperfine coupling is of the wrong sign and is too small. Also the C2-H coupling is too large. However, the hyperfine couplings calculated at the QCISD level of theory are in rather good agreement with the experimental hyperfine couplings. For comparisons the hyperfine couplings computed with DFT are included in the Table 4. The results with both QCISD and DFT theory have the proper spin densities at N3 and at C8, and also seem to agree well with the sites having small spin density such as N7. The calculated spin densities for adenine are illustrated in Figure 4. The spin density at N1 and at C2-H calculated from MP2 is large and negative, and decreases in moving to the QCISD and DFT calculations. Also the spin density calculated with MP2 theory is positive on C8-H, and becomes negative in moving to the QCISD and DFT calculations. The negative spin densities in the QCISD model are best seen by rotating the structure perpendicular to the projection shown here (see Supporting Information, S7). Guanine Cation. Experimentally it has been determined that the spin density of the guanine cation resides on N3, on the exocyclic N2, and on C8 (the numbering scheme used here is shown in the Supporting Information as S4). The dominant hyperfine couplings are at N3(16.8 MHz), N2-H1,2 (-12.1 MHz) and C8-H (-14.6 MHz).34 The experimental hyperfine couplings are given in Table 5. Calculations of these hyperfine couplings at the MP2 level of theory seem somewhat amiss. The N7 hyperfine coupling is too big, and the C8-H hyperfine coupling is too small. One sees that calculation of the hyperfine coupling using the QCISD and DFT level of theory give reasonable agreements with the experimental values especially for the nitrogen hyper-

fine couplings. For the C8-H coupling, the calculations at the QCISD and B3LYP level of theory overestimate the experimental value. This is a problem of doing calculations on an isolated molecule and comparing the results with experimental values determined in a single crystal (which involves an intricate network of hydrogen bonds) When the hydrogen bonding network to nearest neighbors is included the C8-H hyperfine coupling is computed to be -15.35 MHz, in excellent agreement with the experimental value.35 Spin densities computed with each method are shown in Figure 5. One sees the biggest differences in the N7-C8-H region. All of the hyperfine couplings computed at the MP2 level are reflected in the spin densities shown in Figure 5. One sees the negative spin densities at N1 and on the exocyclic N2-H2, and almost no spin density on the C8-H. The C8-H spin density however shifts negative with QCISD or DFT theory. The MP2 calculations show positive spin density at N7, which becomes negative with the QCISD or DFT calculations. Also the figures show that the spin densities are quite similar for the QCISD and DFT calculations. This is expected from the similarities of the hyperfine couplings listed in Table 5. Calculations Involving Stacked Guanines. There are many aspects to the study of stacked guanines. A summary in the Introduction indicates that there is a tendency for some calculations of the cation to spread over several bases in a stack. This is easy to show by doing a calculation on a G:G stack. It is important, however, to first discuss the methods used as it is possible to get different answers from slightly different starting geometries. The first problem faced is having the proper constraints in place for an optimization. By trial-and-error, one sees that the bases tend to move apart, so it is necessary to constrain the bases at a fixed distance to preserve stacking. A good description of the procedure has been given by Kumar and Sevilla.36 For a calculation on a G:G stack with maximum overlap (structure on left panel of Figure 6), one guanine base sits directly above another guanine base (no twisting angle) and the G:G separation is 3.5 Å. For one-electron oxidation an openshell optimization is done on the entire structure. No changes in symmetry occurred with this optimization procedure. With this procedure one will obtain a delocalized cation with a calculation at the B3LYP/6-31+G(d) level of theory (Figure 6, DFT), and a localized cation with a calculation at the UHF/631+G(d) level of theory (Figure 6, UHF). The hyperfine couplings presented in Table 6 show that the DFT calculation has the guanine cation delocalized onto both bases. However, the UHF calculation has the guanine cation localized on a single guanine. Unfortunately, the UHF results do not agree well with the hyperfine couplings determined experimentally (especially for the C8-H) as discussed below.34 It should be noted that while the C8-H hyperfine coupling

1864

J. Phys. Chem. A, Vol. 114, No. 4, 2010

Close

Figure 3. Computed spin densities of the 1-methyl-cytosine cation. All three calculations were performed with the same 6-31+G(d) basis set.

TABLE 4: Calculated Hyperfine Couplings (MHz) for the Adenine Cation a

atom

MP2

QCISD

B3LYP

exp.

N1 N3 N7 N9 N6 N6-H1 N6-H2 C2-H C8-H N9-H

-12.76 -10.93 42.48 -3.21 3.29 -8.57 -8.21 -29.99 2.46 -1.95

-7.14 6.66 4.05 -3.85 11.90 -17.28 -18.00 -17.5 -19.13 -0.48

0.15 9.11 2.80 -1.51 14.90 -20.82 -21.19 -6.94 -16.73 -1.14

neglig. 6.0 neglig. 3.6 10.00 -19.7 -20.5 neglig. -13.7 neglig.

a Experimental results from ref 33. Neglig. indicates a coupling too small to measure.

computed with DFT seems to be in agreement with the experimental results, this is misleading. In the section above on the guanine cation, it was noted that the DFT overestimates the C8-H hyperfine coupling. In Table 5, the C8-H hyperfine coupling was calculated to be ca. -25.00 MHz. So in Table 6, if one has the C8-H spin density spread over both bases, then the C8-H hyperfine coupling would be ca. -12.5 MHz. The calculations of the type presented in Figure 6 are categorized as procedures designed to maintain symmetry. However, one could break symmetry by optimizing a single one-electron oxidized guanine and then placing this base 3.5 Å above a neutral guanine base. This could lead to a localized cation with a DFT calculation. While this is more or less what is expected, there are counter examples. Parrinello and coworkers have calculations on a G:G stack where one guanine is preoptimized as a cation. A calculation at the ROBLYP level of theory produces a hole delocalized over both bases.6 This then points out that one has to be careful in doing these calculations. Another reason for considering a symmetry break in a G:G stack involves ideas from radiation chemistry. In principal, it is difficult to understand how a native cation can be stably trapped. Ionizing radiation produces both an electron and a hole. These charges are free to move about and if they encounter one another they will simply recombine. To stably trap a defect such as a guanine cation it is usually necessary to separate charge and spin. The DNA base cations are good acids, so there is a strong driving force for the native cation to deprotonate. The net effect of deprotonation is that the spin now resides on a neutral base, making it less susceptible to recombination. Experimentally, this process is observed at low temperatures. In order to observe the primary radiation induced defects, the samples are irradiated at 10 K. Even at 10 K, one normally observes the deprotonated cation. The tendency of a cation to deprotonate has important implications in what can be observed experimentally. In the

discussion above, experimental results were given for the spin distributions of various radiation induced free radicals. For example, the study the guanine cation with EPR/ENDOR spectroscopy involved the use of single crystals. Guanine is not very soluble, so it is necessary to grow crystals in acid. This results in a crystal with the guanine base protonated at N7. Irradiation at low temperatures (10 K) ejects an electron leaving the guanine base with a positive charge which then rapidly deprotonates (at N7). The result is called the “guanine cation” and is what would be created in normal one-electron oxidized guanine. Of course the same chemistry holds for one-electron oxidized guanine. There is a strong driving force for this positively charge molecule to deprotonate. A recent paper by Sevilla and co-workers uses 8-deuteroguanine in DNA-oligomers to determine the location of holes and their protonation states at 77 K.37 The results are that in ssDNA the one-electron oxidized guanine exists as an equal mixture of the native guanine cation (G · +) and the guanine deprotonated cation (G-H) · . In dsDNA, the base pair (G · +:C) exists entirely as the proton transferred radical (C+H+):(G-H) · at 77 K. At ambient temperatures there would be an equilibrium between (G · +:C) and (C+H+):(G-H) · . This study further reports that there is no evidence for delocalization of the hole in G:G:G stacks. This is as expected since the deprotonated guanine cation (G-H) · breaks the symmetry in the stacking arrangement. The calculations reviewed above in the Introduction indicate a delocalized guanine cation. In some cases there is a high degree of symmetry, and in other cases there is a deliberate break in the stacking symmetry. For the DFT calculations, it is not easy to determine if the delocalization is a result of the selfinteraction error. The calculations presented in Table 6 were performed at the B3LYP/6-31+G(d) level of theory. There are newer DFT functionals, some of which are said to be useful in overcoming the self-interaction error. For the present study, several methods based on generalized gradient approximations (GGA) were used; BYLP38 (a combination of Becke’s GGA and the correlation functionals of Lee, Yang, and Parr), BH and HLYP39 (Becke’s half + half nonlocalized exchange in combination with the LYP correlation functionals). Also several new density functionals of the M06-class were used; M06-L40 (with no HF exchange), M06-2X41 (with high percentage of HF exchange) and M06HF42 (with full HF exchange). The results are presented in Table 7. To save space, the results only involve the site of primary spin density (C8-H). None of the DFT calculations performed here on the G:G stack yielded a localized guanine cation. The next step is to reexamine the results in Table 6. The DFT calculations show a delocalized cation, while the UHF results show a localized cation. The problem here is that the UHF results cannot be expected to produce good hyperfine

One-Electron Oxidation of DNA Bases and Stacks

J. Phys. Chem. A, Vol. 114, No. 4, 2010 1865

Figure 4. Computed spin densities of the adenine cation. All three calculations were performed with the same 6-31+G(d) basis set.

TABLE 5: Calculated Hyperfine Couplings (MHz) for the Guanine Cation atom

MP2

QCISD

B3LYP

exp.a

N1 N3 N7 N9 N1-H N2 N2-H1 N2-H2 C8-H N9-H

-4.25 8.29 40.42 -7.06 0.28 10.14 -15.12 -16.19 2.06 3.29

-3.61 14.02 -2.47 -7.00 0.76 6.02 -10.10 -8.79 -30.00 2.33

-2.28 10.60 -0.40 -4.16 -0.16 11.78 -8.00 -9.04 -21.71 1.46

neglig. 16.8 neglig. neglig. neglig. 10.00 -12.1 -12.1 -14.6 neglig.

a Experimental results from ref 34. Neglig. indicates a coupling too small to measure.

couplings. So an effort was made to find good post-HF methods that can be used to obtain good hyperfine couplings. Knowing that any calculation of a G:G stack is going to be very demanding, it is necessary to find a trade-off between speed and accuracy. It was shown above that for the single DNA bases QCISD/6-31+G(d) calculations can calculate hyperfine couplings that agree reasonably well with the experimental results. It is possible to perform MP2 calculations on a G:G stack if one has enough scratch disk space. At the 6-31+G(d) level of theory, the storage requirements are ∼9 Gbs of scratch space and require 8 days on a quad processor to complete the calculation. At the MP2/6-31+G(d) level of theory, the majority of the spin density resides on one base. The results are shown in Table 8. The hyperfine couplings shown in Table 8 compare almost exactly with the results shown in Table 5 for the cation in a single guanine base. In this calculation, there is only a very small amount of spin density on the second guanine. For example the hyperfine coupling at C8-H (base 2) is only 1.98 MHz. To perform a QCISD calculation at the 6-31+G(d) on a G:G stack requires >60 Gbs of memory. This can only be done on a parallel machine with lots of memory per node. A preliminary single point calculation using NWChem43 on a parallel machine at ORNL executed, but did not return hyperfine couplings. Unfortunately, the TCE module that performs the QCISD calculations was not developed to evaluate the properties of the wave functions. It is possible to complete a single point QCISD calculation on a G:G stack using a smaller basis set (3-21G). The default convergence criteria on the energy of 10-7 was lowered to 10-6. This calculation requires ∼14 Gbs of scratch space, and 1600 mb of memory for 2 64-bit processors. No optimization was attempted since a single point calculation required 162 h of CPU time. The results of the QCISD/3-21G calculation are shown in Table 8. Again the hyperfine couplings a localized to a single base, and are very close to the values given in Table 5 for one-

electron oxidization of a single guanine base. As a precaution test calculations on a single guanine base were run at the QCISD/6-31+G(d) and 3-21G levels of theory. There were only small differences in the hyperfine couplings between these two calculations. In summary, the G:G model shown in Figure 6 (on the left) was used in a series of calculations. A calculation at the B3LYP/ 6-31+G(d) level of theory produced a delocalized cation, while a UHF calculation produced a localized cation. Further calculations on the same model at the MP2 and QCISD level also produced a localized cation. This is not surprising since MP2 calculations are based on UHF reference functions. However, from the present discussion of localization/delocalization of the guanine cation, it is evident that the results depend on small changes in the starting geometries of the models used. From the new results on a G:G stack, it is possible to comment on the previous results summarized in the Introduction. The results of the Saito et al. study2 involved a single point HF/6-31G(d) on a G:G stack as shown in Figure 6. For this configuration, they report a delocalized HOMO. Since they did not do an open shell calculation, there is no way of determining the localization or delocalization of the cation. However, one can see from the discussion above that for a UHF/6-31+G(d) calculation the cation is localized on a single guanine. The calculations by Parrinello and co-workers show delocalization in a G:G stack with a separation of 2.5 Å and localization with a separation of 6.0 Å using ROHF theory.6 It is difficult to compare these results with the present work since the emphasis here is on calculating hyperfine couplings. ROHF theory is not suitable for calculating hyperfine couplings since it does not properly treat spin polarization in the doubly occupied orbitals. The studies by Barnett et al.,3 and by Bongiorno4 on four base pair duplexes do show delocalization of the cation over neighboring guanines. While these studies are based on DFT, it is not possible to extrapolate from the results presented here on a simple G:G stack to a DNA duplex. It would seem that the randomness in the guanine coordinates in a DNA model would be enough to induce a localized cation. However, one never knows if this is actually the case without doing the actual calculation. Conclusions This study has looked at density functional theory calculations that have a tendency to yield results that delocalize the purine cation in calculations involving stacked bases. UHF calculations on the same models result in purine cations localized on a single base. The emphasis in the present study therefore has been to find convenient methods for performing calculations on stacked bases independent of DFT. While UHF theory does produce localized cations in some calculations on stacked bases, it is unsuitable for the calculation

1866

J. Phys. Chem. A, Vol. 114, No. 4, 2010

Close

Figure 5. Computed spin densities of the guanine cation. All three calculations were performed with the same 6-31+G(d) basis set.

Figure 6. Spin densities calculated for a cation in the guanine:guanine stack. All calculations were performed with the same 6-31+G(d) basis set.

TABLE 6: Calculated Hyperfine Couplings (MHz) for a Guanine:Guanine Stack B3LYP

TABLE 8: Calculated Hyperfine Couplings (MHz) for the Cation in a Guanine:Guanine Stack

UHF a

atom

guanine 1

guanine 2

guanine 1

guanine 2

exp.

N1 N3 N7 N9 N2 C8-H NH2a NH2b

-1.20 4.98 -0.61 -2.22 3.22 -12.95 -3.14 -3.21

-1.06 4.78 -0.55 -2.02 3.60 -11.75 2.93 -2.43

-3.40 17.82 -32.58 -13.06 4.86 -76.14 -6.12 -7.65

0.00 0.11 -0.21 -0.12 0.05 -0.39 -0.40 -0.01

neglig. 16.8 neglig. neglig. neglig. -14.5 -12.1 -12.1

atom

MP2/6-31+G(d)

QCISD/3-21G

N1 N3 N7 N9 N1-H N2 N2-H1 N2-H2 C8-H N9-H

-4.40 7.94 43.97 -7.08 0.24 10.16 -14.45 -15.25 6.88 3.94

-3.76 18.87 -8.74 -8.02 0.43 5.86 -7.10 -7.89 -30.07 2.98

a

Experimental results from ref 34. Neglig. indicates a coupling too small to measure.

TABLE 7: Hyperfine Couplings (MHz) Calculated for a Cation in a Guanine:Guanine Stack HFCC (MHz)

B3LYP

BLYP

BH and HLYP

M06L

M062X

M06HF

C8-H(base 1) C8-H(base 2)

-12.95 -11.75

-10.18 -10.28

-16.82 -17.59

-14.99 -15.16

-10.31 -10.70

-11.31 -9.43

of spin densities. This led to an investigation of post Hartree-Fock methods that properly treat the configuration interaction. For all four calculations on the individual bases presented here, hyperfine couplings calculated at the MP2 level of theory do not agree well with the experimental results. The reasons for this were discussed some time ago by Carmichael; “While UMP2 spin densities are usually closer to observation than those obtained at the UHF level, the degree of accord with experiment is critically dependent on the participation of odd excitations in the correlated wave function and the extent to which the UMP2 procedure approximates the effect of amplitudes due to double replacements”.20 However, UQCISD computed hyperfine couplings agree much better with the experimental results because of the

inclusion of the effect of amplitudes due to single replacements in the UHF determinant, and the concomitant relaxation of the doubles term. Results presented here show that reliable hyperfine couplings can be computed for all four nucleic acid bases at the UQCISD/6-31+G(d) level of theory. However, these UQCISD calculations are far more time-consuming than DFT calculations. In fact, the results presented here show that DFT calculations of hyperfine couplings on the individual oneelectron oxidized DNA bases yield results that are in excellent agreement with the experimental results. The discussion then moved to calculations on G:G stacks. It was shown that if one began with a geometry that produced a localized cation at the UHF level of theory, then post Hartree-Fock calculations on the same geometry also produced a localized cation. While UQCISD theory is practical to use on the isolated nucleobases, it is very time-consuming for a G:G calculation. The preliminary calculations presented here were on a symmetrical G:G stack with maximum overlap, a situation that does not occur in DNA. Future work will continue by applying these results to calculation on stacked bases in geometries more relevant to DNA with the goal of studying the extent of charge delocalization in high symmetry and in symmetry-broken configurations.

One-Electron Oxidation of DNA Bases and Stacks Acknowledgment. Thanks to Ian Carmichael (Rad. Lab., Notre Dame) for help with using the QCISD method for calculating hyperfine couplings. Thanks to Olexandr Isayev and Leonid Gorb (Chem. Dept., Jackson State Univ.) and Tage Øhman (Physics Dept., Univ. Oslo) for help with performing the QCISD calculations presented here. Preliminary QCISD calculations on a G:G stack were performed by Scott Kirkby (Chem. Department, ETSU) and by Bobby Sumpter (Computer Science and Mathematics Div., ORNL). Many thanks to Einar Sagstuen and Tage Øhman (Physics Dept., Univ. Oslo), Ewald Pauwels (Center for Molecular Modeling, Ghent University) and Jerzy Leszczynski (Chem. Dept., Jackson State University) for helpful discussions. This work was completed during the Fall 2009 semester when I was awarded a Non-Instructional Assignment by ETSU. Supporting Information Available: Additional figures displaying spin densities are included. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Saito, I.; Takayama, M.; Sugiyama, H.; Nakatani, K. J. Am. Chem. Soc. 1995, 117, 6406–6407. (2) Sugiyama, H.; Saito, I. J. Am. Chem. Soc. 1996, 118, 7063–7068. (3) Barnett, R. N.; Cleveland, C. L.; Joy, A.; Landman, U.; Schuster, G. B. Science 2001, 294, 567–571. (4) Bongiorno, A. J. Phys. Chem. B 2009, 112, 13945–13950. (5) Gervasio, F. L.; Boero, M.; Parrinello, M. Angew. Chem., Int. Ed. 2006, 45, 5606–5609. (6) Mantz, Y. A.; Gervasio, F. L.; Laino, T.; Parrinello, M. J. Phys. Chem. A 2007, 111, 105–112. (7) In the Parrinello and co-workers papers, only two separations of the stacked bases were used (3.0 and 7.0 Å). In the follow-up to the present study, an investigation will be made of the geometries found in DNA (i.e., 3.5 Å, and various twists as described in ref 1). (8) Adhikary, A.; Kumar, A.; Khanduri, D.; Sevilla, M. D. J. Am. Chem. Soc. 2008, 130, 10282–10293. (9) Improta, R. Phys. Chem. Chem. Phys. 2008, 10, 2656–2664. (10) Bally, T.; Sastry, G. N. J. Phys. Chem. A 1997, 7923–7925. (11) Lægsgaard, J.; Stokbro, K. Phys. ReV. Lett. 2001, 86, 2834–2837. (12) Pacchioni, G.; Frigoli, F.; Ricci, D.; Weil, J. Phys. ReV. B 2000, 63, 054102-1–8. (13) Close, D. M. Radiat. Res. 1993, 135, 1–15. (14) Wetmore, S. D.; Boyd, R. J.; Eriksson, L. A. J. Phys. Chem. B 1998, 102, 9332–9343. (15) Wetmore, S. D.; Boyd, R. J.; Eriksson, L. A. J. Phys. Chem. B 1998, 102, 5369–5377. (16) Wetmore, S. D.; Himo, F.; Boyd, R. J.; Eriksson, L. A. J. Phys. Chem. B 1998, 102, 7484–7491. (17) Wetmore, S. D.; Boyd, R. J.; Eriksson, L. A. J. Phys. Chem. B 1998, 102, 10602–10614. (18) Close, D. M. In Computational Chemistry: ReViews of Current Trends; Leszczynski, J., Ed. ; World Scientific: Singapore, 2003; Vol. 8. (19) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Phys. Chem. 1987, 87, 5968–5975. (20) Carmichael, I. J. Phys. Chem. 1991, 95, 6198–6201. (21) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, Jr., R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andre´s, J. L.; Gonza´lez, C.; Head-

J. Phys. Chem. A, Vol. 114, No. 4, 2010 1867 Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, ReVision A.11; Gaussian, Inc.: Pittsburgh, PA, 1998. (22) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, ReVision E.01; Gaussian, Inc.: Wallingford, CT, 2004. (23) GaussView, ReV. 2.08, Gaussian, Inc, Pittsburgh, Pa. 1999. (24) Parallel Quantum Solutions, Fayetteville, AK. (25) van Lenthe, J. H.; Zwaans, R.; van Dam, H. J. J.; Guest, M. F. J. Comput. Chem. 2006, 27, 927–932. (26) In organic radicals, the unpaired electron is found in π-orbitals, such as in >C · -H. The spin of the unpaired electron partially polarizes the spin of the one of the electrons in the C-H bond in the parallel sense. This results in an opposite or negative spin polarization in the 1 s-orbital of the proton. This gives rise to the isotropic proton hyperfine splitting (typically, 25G or -70 MHz) assigned to any unpaired spin density at the nuclei. (27) Bally, T.; Borden, W. T. Calculation on Open-Shell Molecules: A Beginner’s Guide. In ReViews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; Wiley-VCH: New York, 1999; Vol. 13, pp 197. (28) Kochi, J. K.; Krusic, P. J. J. Am. Chem. Soc. 1968, 90, 7157– 7159. (29) One problem being computing time. MP2 calculations scale as N5 where N is the number of orbitals. QCISD calculations scale as N6. (30) Geimer, J.; Beckert, D. J. Phys. Chem. A 1999, 103, 3991–3998. (31) Close, D.; Forde, G.; Gorb, L.; Leszczynski, J. Structural Chemistry 2003, 14, 451–454. (32) Close, D. M.; Hole, E. O.; Sagstuen, E.; Nelson, W. H. J. Phys. Chem. 1998, 102, 6737–6744. (33) Nelson, W. H.; Sagstuen, E.; Hole, E. O.; Close, D. M. Radiat. Res. 1992, 131, 272–284. (34) Close, D. M.; Sagstuen, E.; Nelson, W. H. J. Chem. Phys. 1985, 82, 4386–4388. (35) Unpublished results (DMC). (36) Kumar, A.; Sevilla, M. D. J. Phys. Chem. B 2006, 110, 24181– 24188. (37) Adhikary, A.; Khanduri, D.; Sevilla, M. D. J. Am. Chem. Soc. 2009, 131, 8614–8619. (38) Becke, A. D. Phys. ReV. A 1988, 38, 3098–3100. (39) Becke, A. D. J. Chem. Phys. 1993, 98, 1372–1377. (40) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101-1-18. (41) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241. (42) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 13126–13130. (43) Bylaska, E. J.; de Jong W. A.; Govind, N.; Kowalski K.; Straatsma, T. P.; Valiev M.; Wang, D.; Apra, E.; Windus, T. L.; Hammond J.; Nichols, P.; Hirata S.; Hackler, M. T.; Zhao Y.; Fan, P.-D; Harrison, R. J.; Dupuis, M.; Smith, D. M. A.; Nieplocha, J.; Tipparaju, V.; Krishnan, M.; Wu Q.; Van Voorhis, T.; Auer, A. A.; Nooijen, M.; Brown, E.; Cisneros, G.; Fann, G. I.; Fruchtl, H.; Garza, J.; Hirao, K.; Kendall, R.; Nichols, J. A.; Tsemekhman, K.; Wolinski, K.; Anchell, J.; Bernholdt, D.; Borowski, P.; Clark, T.; Clerc, D.; Dachsel, H.; Deegan, M.; Dyall, K.; Elwood, D.; Glendening, E.; Gutowski, M.; Hess, A.; Jaffe, J.; Johnson, B.; Ju, J.; Kobayashi, R.; Kutteh, R.; Lin, Z.; Littlefield, R.; Long, X.; Meng, B.; Nakajima, T.; Niu, S.; Pollack, L.; Rosing, M.; Sandrone, G.; Stave, M.; Taylor, H.; Thomas, G.; van Lenthe, J.; Wong, A.; Zhang, Z. NWChem, A Computational Chemistry Package for Parallel Computers, Version 5.1; Pacific Northwest National Laboratory: Richland, Washington 99352-0999, USA, 2007.

JP906963F