Solvent Fluctuations Drive the Hole Transfer in DNA - American

Sep 2, 2009 - We studied the hole transfer across adenine bridges in double-stranded DNA by means of a multiscale approach, propagating the hole in th...
3 downloads 0 Views 5MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 13107–13117

13107

Solvent Fluctuations Drive the Hole Transfer in DNA: A Mixed Quantum-Classical Study Toma´sˇ Kubarˇ,† Ulrich Kleinekatho¨fer,‡ and Marcus Elstner*,† Institute of Physical Chemistry, Technische UniVersita¨t Braunschweig, 38106 Braunschweig, Germany, and School of Engineering and Science, Jacobs UniVersity Bremen, 28759 Bremen, Germany ReceiVed: July 31, 2009

We studied the hole transfer across adenine bridges in double-stranded DNA by means of a multiscale approach, propagating the hole in the framework of time-dependent DFT coupled to classical molecular dynamics simulation using a QM/MM scheme. The hole transfer in DNA is codetermined by the large fluctuations of site energies on the order of 0.4 eV, induced by the solvent degrees of freedom. These fluctuations lead to charge-transfer active conformations with large transfer efficiency, which are characterized by a favorable alignment of site energies along the DNA strand. This reduces the barrier for the hole transfer dramatically. Consequently, we find that a charge hopping mechanism is operative already for short bridges with fewer than four adenines, in contrast to the charge-transfer models assuming static DNA structures, where only tunneling occurs. The solvent fluctuations introduce a significant correlation between neighboring sites, enhancing the charge-transfer rate, while the fluctuation of electronic couplings has only a minor impact on the charge-transfer characteristics. Our results emphasize the importance of an accurate description of solvent effects as well as proper sampling, and it is suggested that charge transfer in DNA is gated by the dynamics of solvent. Introduction Apart from the primary role as the carrier of genetic information, DNA has been the focus of intensive research in the past decade due to its ability to transport electric charge over large distances. Long-range hole transfer was discovered in natural double-stranded DNA1-4 and also probed in singlemolecule conductivity experiments.5-7 This phenomenon was proposed to be of fundamental significance for such processes like oxidatively generated damage and DNA repair.8-10 Indeed, the electrochemical properties of DNA have been studied since 50 years ago,11 and electrochemical bionsensors based on DNA have been constructed.12,13 Besides that, the suggested application of DNA in the design of nanoelectronic devices relies on sufficient electric conductivity of DNA.14,15 Charge transport in DNA is determined by the energetics of excess charge on the nucleobases. Since the ionization potential (IP) of guanine (G) is the smallest, followed by the 0.4 eV larger IP of adenine (A),16-18 a hole will predominantly occupy the low-lying G sites. The hole transport may then be considered between G’s across a run of A’s, with the pyrimidine bases standing out of the way of transport. Depending on the number of intervening A’s, two distinct mechanisms of such transport were proposed. On the one hand, the hole can tunnel through short A tracts, as investigated in the experiments by Giese et al.19 on GTnGGG sequences. The reported exponential decay of charge-transfer (CT) rates for A bridges with up to three bases was explained in terms of the one-step superexchange-mediated tunneling through the short bridges. In this model, the hole is transferred from one G to another, tunneling through intermediate bridges of A’s (G hopping).19-24 However, this mechanism becomes inaccessible for longer A bridges. Here, the experiments by Giese et al. show weaker distance dependence, * To whom correspondence should be addressed. E-mail: [email protected]. † Technische Universita¨t Braunschweig. ‡ Jacobs University Bremen.

indicating a different transport mechanism. The thermally induced hopping model (TIH) proposes the hole to be thermally activated to the A tract, where it hops from one A to another. The relatively large activation energy of 0.4 eV leads to a quite small Boltzmann factor for the transfer of a hole to A.25 Nevertheless, CT over longer A tracts26,27 can only be explained assuming a thermal excitation of the hole from G to A, followed by a sequential hopping along the A tract (A hopping).28-31 Although thermal energy plays a crucial role in the transport mechanism, all early theoretical studies were based on static structures, mostly idealized B-DNA. For these geometries, the relevant CT parameters could be readily evaluated.32 These are the site energies εi, which are given by the IP of the individual bases, and the electronic couplings Tij between neighboring sites (which are in the range of 0-0.1 eV). They are the key to describe CT in the tunneling as well as in the incoherent hopping regime, and the kinetics of CT in DNA may be described by means of the Marcus theory33,34

k)

[

(∆G0 + λ)2 1 2π |HDA | 2 exp p 4λkBT √4πλkBT

]

(1)

where the CT parameters enter the formula by means of the energy difference ∆G0 and the bridge-mediated electronic coupling HDA, and the reorganization energy λ is the energy required to change the molecular structure along the reaction coordinate. It has two components, corresponding to the change of molecular structure of the nucleobases themselves (internal λ) and to the reorganization of solvent (solvent λ). The CT parameters εi and Tij can also be cast into a tight-binding Hamiltonian of the form

10.1021/jp9073587 CCC: $40.75  2009 American Chemical Society Published on Web 09/02/2009

13108

J. Phys. Chem. B, Vol. 113, No. 39, 2009

ˆ) H

∑ εiaˆi+aˆi + ∑ Tijaˆi+aˆj i

Kubarˇ et al.

(2)

i*j

where aˆ+ i is the creation operator for the excess charge on site i and aˆi is the annihilation operator. The distance dependence of charge transfer across the A tracts has been studied in the framework of Marcus theory35-40 as well as using tight-binding theory29 and kinetic-quantum models.30,31 However, all of these model calculations assume static structures and thus constant values of the CT parameters. In contrast to the static picture, large fluctuations of the CT parameters have been reported in all computational studies, where those have been computed for snapshots along classical MD trajectories.25,41-51 This raises the question, how do the observed fluctuations affect the hole transfer and, in particular, the CT model proposed on the basis of a static description? Recently, electron-transfer reactions in biomolecular systems have been revisited by Balabin et al., emphasizing the role of structural fluctuations on CT and concluding that a few structures with the largest values of HDA can dominate the CT characteristics.52 In the case of proteins, it seems that no particular mode has been identified to promote CT. In the case of DNA, however, several microscopic CT models postulate an active role of specific molecular motions for CT. All of these models focus on various singular aspects of CT and propose quite different molecular modes to drive the transfer. For instance, Barton et al. proposed that a collective motion of DNA (a “local” compression of the DNA structure) could facilitate the hole transfer. This model suggests the existence of certain CT-active conformations, in which the electronic couplings between sites tend to be significantly enlarged (conformational gating model).53,54 In another concept, Schuster et al. highlighted the impact of counterion motion on CT within the ion gating model.55,56 Here, the motion of counterions and the electric field induced by them causes the hole to be either propagated along the strand or pinned to an individual nucleobase. As analyzed so far computationally, the fluctuations of CT parameters are caused by two factors; on the one hand (i), they depend sensitively on conformational changes of the DNA double helix, and on the other hand (ii), they are influenced by the interaction with the solvent, that is, water molecules and counterions. In computational studies, both factors can be investigated independently, and in particular, the interaction with solvent is not included in all theoretical studies. (i) The variability of DNA conformation leads to fluctuations of Tij on the order of 0.05 eV,42,45,47-51,57,58 and those of εi are on the order of 0.1-0.15 eV. Fourier analysis of the autocorrelation function shows the site energies to oscillate with a frequency of 1600 cm-1 (period of 20 fs). This is typical for the double bond vibrations of the molecular skeleton,58 which are not correlated between the individual nucleobases.51 (ii) The interaction of the DNA electronic structure with the solvent has still been neglected in many computational studies.42,45,47-51 If incorporated, the coupling to the solvent and ionic motion introduces remarkable fluctuation of εi on the order of 0.4 eV.25,46,58 These large oscillations with a period of 40 fs were shown to be induced by the fluctuation of solvent,58 corresponding most likely to the vibration of hydrogen bonds, and the fluctuations of neighboring sites are highly correlated. Besides that, a recent study shows the hydrogen bonding with water to control the localization of the hole as well.59 A more detailed analysis showed that the motion of ions, which induces electric field changes in the picosecond regime (compared to the 40 fs modes of water), is grossly shielded by the water motion.58 Thus,

the effects of water and ions are highly anticorrelated, and no evidence for the ion-gated model could be found in the data analyzed so far. The interaction with water and counterions, on the other hand, has no important impact on Tij, and no correlation between Ti-1,i and Ti,i+1 could be detected. Therefore, no evidence for a local compression could be found as proposed in the conformational gating model.58 Until recently, the effect of dynamic disorder was only studied qualitatively, for example, using model Hamiltonians to solve the time-dependent Schro¨dinger equation (TDSE).22,60,61 On the basis of classical MD, two recent computational studies have particularly challenged the static picture of CT in DNA. Both studies used CT parameters were derived along MD trajectories, thereby including the parameter fluctuations due to DNA conformational changes (point (i)) but still neglecting the direct coupling to the solvent (ii). Hatcher et al.51 reported that the fluctuations of εi of 0.1-0.15 eV can introduce CT-active conformations with a vanishing energy barrier between G and A, leading to a significant hole density on the bridge sites. This finding suggests that the A hopping mechanism may become relevant already for short A tracts. Grozema et al.45 studied CT using a tight-binding model based on classical MD trajectories. The referenced experimental data could only be reproduced if dynamical effects had been included, proving the incompleteness of the static bridge model. Further, the strong distance dependence of the CT rate for short A tracts did not result from tunneling but rather from the electrostatic interaction of the hole charge with the negatively charged hole donor in the particular DNA species. Although a lot of new knowledge has been gained due to the increasingly realistic simulations in the last years, several questions still remain open. A key question concerns the impact of fluctuations of εi on the order of 0.4 eV, induced by the solvent, which matches the energy difference between the IP of G and A. Further, the importance of correlation between neighboring sites and the mechanistic role of Tij fluctuations have to be investigated in more detail. To gain a principal insight into the CT mechanism, a computational model has to be applied that involves, on the one hand, the CT parameters εi and Tij evaluated with a good accuracy and, on the other hand, the account for dynamical and environmental effects. The environmental effects can be included by using combined QM/MM methods,62 as has been already done for CT in DNA by several authors.25,46,58 To describe the dynamical effects, the CT parameters have to be calculated along MD trajectories extended for pico- to nanoseconds, which can only be achieved using approximative or semiempirical methods. We have shown recently that the CT parameters obtained with the approximative DFT method SCC-DFTB63 are quite accurate and therefore comparable to those computed with higher-level methods. The dynamical and environmental effects are very well captured with SCC-DFTB/MM, which compares favorably to full DFT QM/MM methods.57 In this work, we apply the SCC-DFTB based coarse-grained DFT Hamiltonian to study the hole transfer in the framework of the time-dependent density functional theory (TD-DFT). The computational efficiency and accuracy of this method allows us to simulate multinanosecond dynamics of CT, a prerequisite to study the hole transfer between G’s across A tracts of varied length. Computational Methodology The subjects of this study were double-stranded DNA species of sequences GTnGGG with n ) 1, 2, 3, 4, 5, 7, 10, and 14, which were previously used in the experimental work of Giese

Solvent Fluctuations Drive the Hole Transfer in DNA

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13109

et al.27 Similar DNA sequences (in the form of hairpins) were used in the experimental work of Lewis et al.64 Following these experiments, we studied the hole transfer from the first G to the terminal GGG triplet, passing across a bridge formed by an intervening run of A’s on the other strand of the duplex. The hole was propagated in a mixed quantum-classical scheme employing a coarse-grained tight-binding (TB) model as developed recently.57 First, a classical MD simulation of the DNA species in question was performed, followed by the calculation of CT parameters (TB Hamiltonian consisting of site energies and electronic couplings) based on the SCC-DFTB method57,63 and the fragment-orbital approach.61,65 Here, the HOMO energies εi represent the site energies. The electronic couplings Tij and the overlap integrals Sij are computed from the HOMO orbitals φi obtained for isolated nucleobases in the atomic orbital basis (φi ) ∑µ cµi ηµ), as follows

ˆ |φj〉 ) Tij ) 〈φi |H

∑∑ µ

ˆ |ην〉 cµi cjν〈ηµ |H

)

ν

∑∑ µ

cµi cjνHµν

ν

(3) Sij ) 〈φi |φj〉 )

∑ ∑ cµi cjν〈ηµ|ην〉 ) ∑ ∑ cµi cjνSµν µ

ν

µ

ν

(4) The use of Tij instead of Hµν to describe the electronic system can be viewed as a coarse graining of the electronic description; Tˆ is thus a coarse-grained Kohn-Sham Hamiltonian. For reasons of applicability in a simplified Schro¨dinger equation, this Hamiltonian needs to be orthogonalized, and this is performed by using the Loewdin transformation in the form66

T' ) S-1/2TS-1/2

(5)

As described in our previous work on this methodology, the diagonal elements of the Hamiltonian matrix, which correspond to the site energies, were taken from the SCC-DFTB calculation of the individual nucleobases. The mentioned orthogonalization procedure thus affects only the off-diagonal elements. ˆ |ην) ) 〈ηµ|H The evaluation of the AO Hamilton matrix HDFTB µν using the DFTB Hamilton matrix, which is precalculated and stored,63 makes the described calculation extremely efficient. Thus, this procedure can be used to calculate CT parameters for snapshots generated by a classical MD simulation even for several nanoseconds (as applied below). For the calculation of the electronic couplings, the LCAO basis used in DFTB has been optimized,57 leading to an excellent agreement with higherlevel DFT61 and CAS-PT2 methods.67 As discussed above, the coupling of electronic parameters to the environmental degrees of freedom is crucial for a reliable description of CT in DNA. This coupling is efficiently introduced using a QM/MM implementation68

1 QM/MM DFTB Hµν ) Hµν + Sµν 2

∑ (Qδ/RRδ + Qδ/Rβδ)

(6)

δ

where Qδ are MM charges in the environment, Sµν is the AO DFTB is the AO Hamilton matrix overlap matrix from DFTB, Hµν elements calculated in DFTB, and RRδ is the distance between the QM atom where the AO orbital ηµ is located and the MM DFTB atom with charge Qδ. In our application, the term Hµν describes the interaction between the individual nucleobases,

while the DNA backbone, the water molecules, and the counterions are represented by their respective partial charges Qδ in the QM/MM term. Therefore, inserting eq 6 into eq 3 accounts for a full coupling to the environment. To study the effect of solvent and counterions independently, we will discard DFTB will be the QM/MM term in some simulations, that is, Hµν ; this will be denoted as in vacuo. used in eq 3 instead of HQM/MM µν Note that this simplified approach is equivalent to that used in some of the previous studies by other authors.45,51 This time-dependent coarse-grained Kohn-Sham Hamiltonian Tˆ′ (eq 5) is then cast into the TDSE for the excess charge (hole)

˙ ) Tˆ′Ψ ipΨ

(7)

Here, Ψ is the wave function of the hole, Tˆ′ is the TB Hamiltonian, Tii ) εi is the site energy of nucleobase i, and Tij is the electronic coupling of sites i and j. Time-dependent equations like this one can be derived in the framework of timedependent DFT using the action principle; the formal justification has been given by Runge and Gross.69 Recently, a TDDFT version for SCC-DFTB has been developed,70 and eq 7 can be seen as a coarse-grained non-self-consistent version of this TD-SCC-DFTB method. The wave function of the hole was expanded in the basis of HOMO orbitals φi on the considered purine nucleobases, Ψ ) ∑i aiφi. The initial state of the hole was localized on the first G, representing a charge injection at this site. An extra negative imaginary potential was applied on the diagonal elements of the Hamiltonian that correspond to the G’s of the terminal triplet (in the same way as Grozema et al.).45 This extra term in the Hamiltonian (Hii ) εi - iτ for i ∈ GGG, τ ) 1000 ns-1) models the trapping of the hole by means of a reaction with water, with rate constant 2τ (q˙i ) ... - 2τqi). The magnitude of τ is discussed in the Supporting Information. Finally, this TDSE was integrated using the Runge-Kutta predictor-corrector integrator implemented in RKsuite,71 yielding the time evolution of the hole. The time-dependent LCAO coefficients ai then determined the charges on the individual nucleobases qi ) |ai|2, and the norm of the wave function

P(t) ) |Ψi(t)| 2 )

∑ |ai(t)|2 ) ∑ qi(t) i

(8)

i

gave the total magnitude of charge remaining in the system, in other words, the hole survival probability. For every DNA species, we proceeded as follows: First, a 50 ns classical MD simulation was performed with the Amber parm99+parmBSC0 force field.72,73 The AMBER 9 simulation package74 was used for the simulation setup, while the simulations were run with the GROMACS 4 package75 (cf. the Supporting Information for details on the simulation setup and the equilibration procedure). One-hundred frames from the original run (separated by 0.5 ns) were used as starting points for 20 ps simulations with the time step of 1 fs. Within these windows of 20 ps, the CT parameters were evaluated, which were subsequently used to propagate the hole according to eq 7. Therefore, 100 independent quantum-classical simulations of 20 ps in length were performed for every DNA species. Such a large number of uncorrelated simulations was chosen to guarantee sufficient sampling and results free of statistical errors. For every DNA sequence considered, it was necessary to

13110

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Figure 1. (a) The survival of hole P(t) versus time in GTGGG. Data from 100 simulations, 20 ps each. (b) Survival with the static model. Note the 50× longer time scale.

evaluate the Hamilton matrix Tˆ′ two million times, showing the necessity of an efficient computational scheme.

Kubarˇ et al.

Figure 2. (a) The occupation of the bridge A by the hole in GTGGG; the averaged time dependence from dynamical simulations and the result with the completely static model. (b) The occupation of the A bridge in all GTnGGG sequences; the averaged time dependence from dynamical simulations, labeled with the value of n.

Results Dynamics of the Hole; The Rate of Transfer. In a first set of simulations, the hole was propagated through the sequences GTnGGG, with n ) 1-14, using eq 7 with the Hamiltonian defined in eqs 3 and 6. Figure 1a shows the hole survival probability P(t) (cf. eq 8) in GTGGG for the 100 trajectories of 20 ps, as described above. A large variance in the rate of survival decay is observed within a single 20 ps simulation as well as between the (100) individual simulations. The curves in Figure 1a feature distinct sections. There are intervals where the survival remains moreor-less constant, that is, no charge is being annihilated, resulting from a negligible rate of transfer across the bridge. On the other hand, other intervals are observed where the survival is dropping steeply as a consequence of the hole transfer at a fast rate. For comparison, Figure 1b shows the survival probability for a static model, where the CT parameters were fixed at their mean values. Note that the decay is about 2 orders of magnitude slower in this case. Grozema et al.45 discussed this case as well, and they concluded that this type of static model is in disagreement with the experimental data available. Figure 2a shows the transient occupation of the bridge site for the dynamic and static model. While the static model shows a nearly vanishing occupation, indicating tunneling through the bridge, the dynamic model exhibits a more substantial occupation of the bridge site, which suggests that the hopping mechanism is operative already for this bridge consisting of a single A. A similar picture arises for longer sequences, as shown in Figure 3 for GT5GGG and GT10GGG. Here, the results with the static model are not shown as no decay could be detected

Figure 3. Survival of the hole P(t) versus time in GT5GGG (a) and GT10GGG (b).

within the available simulation time. Apart from the features discussed in the case of GTGGG, the dynamical simulations

Solvent Fluctuations Drive the Hole Transfer in DNA

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13111

Figure 4. Simulation of hole transfer across A bridges. Presented are the survival probability after 20 ps and the rate constant of hole transfer, as mean values with standard deviations.

Figure 5. Rate constant of hole transfer in GTnGGG. Data from full MD-based calculations as well as those based on constant site energies and on constant electronic couplings.

show much longer hole survival times. As a consequence, the occupation of bridge sites is significantly increased (cf. Figure 2b), indicating that hole hopping becomes the dominant transport mechanism. For any sequence, the CT dynamics in the full QM/ MM framework differs dramatically from the static model. Assuming first-order kinetics for the hole transfer, the rate constant may be calculated from the survival P(t) as k ) -log(P(t))/t, t ) 20 ps, for every individual simulation, and Figure 4 shows the mean and standard deviation of the 100 values obtained for every DNA species. The survival of the hole is considerably depleted already within 20 ps of simulation, even in long A bridges (cf. Figure 3b), and the simulation time of 20 ps is thus considered sufficient. The rate constant is also averaged over smaller numbers of individual simulations (25 and 50) to check the convergence behavior of these values with the increasing number of values; see Supporting Information. This simple test proved satisfactory convergence of the rate constants, taking into account also the relatively large standard deviation of the data. The obtained CT rate exhibits a significant decrease with the length of the A bridge, as expected. However, the distance dependence is algebraic already for short bridges and not exponential as reported by Giese et al.,27 who observed the algebraic dependence only for longer A tracts (with more than three A’s). Alas, a comparison of absolute values of the CT rate is not possible as Giese et al. reported their relative magnitude only. Our results look similar to those of Lewis et al.64 Our rate of transfer across the shortest bridge agrees with the rate of hole arrival reported by Lewis et al. (around 400 ns-1), while the distance dependence of the rate is too weak as well. (It should be noted here that Lewis et al. performed a slightly different reaction than what we simulated; their hole donor kept a negative charge, which may have then affected the transfer of the hole; there was no such negative charge in the system that we dealt with. Grozema et al.45 took this negative charge into account and obtained a correct shape of the distance dependence but with a 2 orders of magnitude faster transfer for the long bridges.) This inaccuracy of our simulations is caused most likely by the missing response of solvent to the hole charge (polarization of the solvent). The contribution of distance-dependent solvent reorganization energy λs to the exponential distance dependence of the CT rate was already established,76 and the inclusion of a correct distance dependence of λs brought about the desired distance dependence of the CT rate in a recent study using a more derived model of CT.40 We have recently estimated λs using classical MD simulations,77 and the observed steep

increase of λs for short A bridges (up to three A’s) would lead to a steeper, exponential-like decrease of the CT rate with distance. Obviously, it is necessary to include the response of the solvent in the computational model in order to obtain reliable absolute rates of CT, and to solve this rather difficult task is an ongoing effort in our research group. The Effect of Site Energies and Electronic Couplings. The impact of fluctuations of εi and/or Tij on CT has been discussed in many recent studies. Balabin et al.52 and Rˇeha et al.47 emphasized the effect of fluctuations of Tij on CT, while Hatcher et al.51 proposed the skeletal motions of nucleobases to “gate” CT via the fluctuations of εi. In the following, we will investigate the relative importance of these two factors in detail. To address this issue, we performed three sets of quantum simulations, (i) with the QM/MM CT parameters evaluated in every MD time step as discussed above, (ii) with the site energies εi from these calculations while keeping the values of electronic couplings Tij constant at their mean values, and (iii) taking the electronic couplings Tij from the QM/MM MD and keeping the εi constant at their mean values. The estimated rate constants obtained within these models are summarized in Figure 5. (For the data in tabular format, see the Supporting Information.) Interestingly, the simulations (ii) with constant Tij give quite similar results as the regular dynamics (i). Therefore, at first sight, the fluctuations of couplings have only minor influence on the CT at least for short A bridges, and they cannot be regarded as the major driving force. However, the observed rate is overestimated by ∼60% for the longest bridge if the couplings are held constant, indicating the necessity of appropriate description of the couplings. This is probably the consequence of their very broad distribution reported, for example, in our previous work,58 which cannot be approximated by a single (mean) value realistically. On the contrary, the computational model may be required to capture the transient moments in the simulation where the couplings are minimized or maximized, far from their mean values. Consistent with that is the finding of Balabin et al.,52 who observed CT in model protein systems being controlled by dynamical disorder, expressed in terms of a low coherence parameter of the overall donor-acceptor electronic coupling. Anyway, the effect of fluctuations of Tij is small compared to those of εi. On the other hand, constant site energies εi in simulations (iii) nearly suppress CT, with the rate converging to zero at machine precision already for GT4GGG. Evidently, the variance of the energy profile (determined by the site energies varying in time) is vitally important for the hole transfer to occur. The

13112

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Kubarˇ et al.

Figure 6. A “slow” (left) and “fast” (right) example of a 20 ps simulation of hole transfer in GTGGG. Electronic coupling of G1 and the bridge A (top), difference of site energies of G1 and A (center), as well as the occupation of these bases and the hole survival (bottom). The vanishing barrier leads to abrupt decay of survival if the electronic coupling is sufficient (vertical arrows); no decay is observed with modest coupling (slanted arrows).

idea of transfer passing over a barrier of constant height would be incorrect, leading to too small transfer rates. Rather, the large fluctuations undergone by the site energies will make the barriers fluctuate as well, which in turn will speed up the transfer considerably. The findings reported so far suggest the existence of certain CT-active conformations on a CT-silent background as proposed by Barton et al. yet in a more generalized way that shall involve not only the DNA conformation but also the configuration of environment (water and counterions). These CT-active conformations are then characterized by favorably adjusted site energies. To make this point more transparent, we can analyze selected trajectories in more detail. Figure 6 shows a couple of 20 ps simulations of hole transfer in GTGGG. In the left panel is shown a slow trajectory with modest mean electronic coupling between the initial G1 and the bridge A (of 4 meV compared to the average of 13 meV), so that the individual CT events are clearly separated. There is virtually nothing happening for extended periods of time, while a significant decay of survival occurs in several quite short intervals. These phases of observable decay always follow a transfer of a part of the hole from G1 to A. Concentrating on these events, we find that the difference of site energies of G1 and A vanishes in these moments (cf. vertical arrows), making a barrierless transfer between these sites possible. This observation provides support to the idea of CT-active conformations consisting of the favorably adjusted site energies. Interestingly, some intervals with the difference of site energies near zero exhibit no measurable transfer, which is caused by insufficient electronic coupling at those points (cf. slanted arrows). It has to be noted again that this particular simulation exhibits electronic coupling below the average, which makes the last-mentioned events more frequent and more visible. On the other hand, the simulation shown in the right panel is “fast”, meaning that the entire hole comes through in a sequence of barely separated transfer events,

virtually all at once. This is a consequence of the vanishing energy barrier between G1 and A, accompanied by a sufficient value of the electronic coupling (with a mean value of 11 meV). The Role of Solvent and Its Dynamics. To assess the effect of the electric potential induced by the solvent on the rate of hole transfer rigorously, we performed a similar set of simulations as that in the first section, neglecting the QM/MM term (representing the environment) in eq 6. Consequently, we obtained the time series of site energies with much smaller fluctuations on the order of 0.1-0.15 eV, similar to the results of other studies, where the QM/MM term was neglected.45,47,51 The distribution of site energies of the initial G1 and the bridge A is shown in Figure 7a, for both cases. As can be seen, the overlap of these distributions increases significantly upon the inclusion of the QM/MM term, in agreement with other studies.25,46 It is important to note that there is clear methodological convergence. The computational investigations mentioned above use different Hamiltonians; however, very similar values of electronic couplings as well as site energies are obtained. The main difference represents the inclusion of the QM/MM term, which brings on the fluctuations of site energies of 0.4 eV consistently, compared to 0.1-0.15 eV obtained in vacuo. The strong correlation of site energies with the electrostatic potential induced by the environment may be inferred from the time series of both quantities presented in Figure 7b. The correspoding correlation coefficients take values larger than 0.9, much the same as that reported previously.58 The resulting rate constants from both sets of simulations (with and without the QM/MM term) are presented in Figure 8 for all sequences. We can see that the CT is slower if the environment is excluded from calculations. However, the CT rate remains far from negligible, which was not the case if the site energies were held constant. The fluctuations are small though, yet they are still sufficient to facilitate hole transfer. The difference of CT rates in Figure 8 is more prominent in the short A bridges. This observation may seem surprising, yet

Solvent Fluctuations Drive the Hole Transfer in DNA

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13113

Figure 9. Difference of site energy of the initial G1 and the bridge A (the barrier), calculated for GTGGG with the inclusion of an environment (QM/MM) and without that (in vacuo).

Figure 7. (a) Distribution of site energies of the initial G and of the bridge A in GTGGG, calculated in vacuo (top) as well as using the full QM/MM model (bottom). (b) Example of the time course of site energies of these nucleobases (ε, red) as well as the electrostatic potential induced on them by the environment (ESP, black).

along the A bridge. These short-lived barriers may then obstruct the further transport of the hole that has already been excited to the bridge. Therefore, the large fluctuations of site energies constitute a double-edged factor; whereas they facilitate the excitation of the hole from the initial G to the A bridge, they may hinder the transfer of the hole along the A bridge itself. Still, the error brought about by the neglect of the environment is significant. The inclusion of the environment within the model of external point charges does not compromise the computational efficiency, and thus, it should be applied in studies of this kind. Correlations along the Bridge. In our previous work, we found that the onsite energies of neighboring nucleobases in the double strand are correlated.58 In a simple picture, this means that the fluctuations of neighboring sites are in phase. Such a pattern would highly enhance the CT compared to the situation with independent fluctuations of neighboring sites. We evaluated the correlation coefficients Fij for the site energies of individual nucleobases εi, εj

Fij )

Figure 8. Rate constant of hole transfer in GTnGGG with parameters calculated with the inclusion of an environment (QM/MM) and without that (in vacuo).

a simple explanation is possible. In shorter bridges, it is probably the excitation of the hole onto the A bridge that determines the rate of the whole process. If the parameters are calculated in vacuo, the intrinsic difference of IP of A and G is retained (due to the smaller fluctuations of IP), and the initial barrier is more difficult to overcome than if the parameters are calculated within the QM/MM formalism. Indeed, this can be clearly seen in our simulations (see Figure 9). Within QM/MM, the barrier reaches values close to zero quite often, while these events are relatively rare without the inclusion of environment, making the CT-active conformations appear much less frequently. In longer bridges, on the other hand, the situation may be reversed. The QM/MM treatment leads to larger fluctuations of site energies, thereby introducing larger instantaneous barriers

〈(εi - 〈εi〉) · (εj - 〈εj〉)〉

√〈(εi - 〈εi〉)2〉 · 〈(εj - 〈εj〉)2〉

(9)

Much the same as in the case of poly(G) DNA in our previous work,58 we found a strong correlation of neighboring sites, with values of Fij around 0.7 between neighbors and 0.4 between second neighbors. (see the Supporting Information for the complete data.) A very interesting and surprising feature is that whereas the A’s are highly correlated among themselves (i.e., neighboring site energies oscillate in a correlated fashion), the A’s are much less correlated with the neighboring G’s. This may be another important factor affecting the height of the initial barrier. If the initial G1 and the first A were highly correlated as well, then the energy difference of 0.4 eV between G and A (which the hole has to overcome) would be less likely to vanish. Since the large fluctuations are caused by the fluctuating electric field generated by the solvent, it is not unplausible to presume that the A’s and the G’s respond differently to the field, probably due to their very different polarity. To estimate the effect of correlation, a similar set of quantum mechanical simulations as those in the previous sections was performed. Prior to that, the distributions of both site energies and electronic couplings were inspected along all of the trajectories analyzed above. The normality of all distributions was checked visually and confirmed in all cases; see Supporting

13114

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Kubarˇ et al.

Figure 10. Comparison of the results for the dynamical and statistical models.

Information for more details. Then, artificial data sets of site energies and electronic couplings were created by generating values using a random-number generator, in the way proposed by Grozema et al.45 A value was drawn from the respective normal distribution, and the time of life s of this value was drawn from an exponential distribution P(s) ) 1/σ exp[-s/σ] (σ ) 100 fs). The particular value of the CT parameter was kept in the data set for the time interval s, after which a new value of the parameter was obtained. This sequence of steps was applied repeatedly until the length of the data set reached 20 ps. This procedure ensured the decorrelation of CT parameters. Finally, 100 simulations of 20 ps were performed for every DNA species using these data sets. The resulting rate constants are summarized in Figure 10 in comparison with the results of QM/MM simulations. The decorrelation of parameters leads to the decrease of the CT rate for all bridge lengths. This decrease is relatively small in short bridges (∼30%) but grows as the bridge becomes longer. The hole transfer across a bridge consisting of 14 A’s is then rendered considerably slower at ∼40% of the rate obtained in simulations based on true time series of CT parameters, which may be clearly seen in the plots of hole survival; see Figure 11. This is probably a consequence of the discarded temporal course of site energies, in particular, the missing correlation of neighbors. The statistical models make it also possible to decouple the effect of the correct time course of site energies and electronic couplings by generating the parameter sets in the following way; the site energies may be drawn randomly from the respective distributions while taking the full time series of electronic couplings constant, or vice versa, or both εi and Tij sets may be generated from the distributions. The quantum simulations performed with these (however artificial) time-dependent Hamiltonians for the sequences GTGGG and GT14GGG yielded the rate constants presented in Table 1. If we discard the correlation of neighboring site energies by drawing these from the distribution (model dist, MD in Table 1), then we observe a decreased CT rate. The decrease is more distinct for the long bridge (∼70%) than that for the short one (∼35%) but is still considerable in both cases. Obviously, it is not quite possible to disregard the actual time course of εi. On the other hand, the use of a statistically generated series of electronic couplings brings on a smaller change of CT rates, making the transfer slightly faster. On the basis of this observation, it seems possible to build a model of hole transfer in DNA which involves the electronic couplings being drawn from predetermined distributions, making the evaluation of Tij unnecessary at the price of slightly overestimated CT rates. Also, these results demonstrate

Figure 11. Survival of the hole in GT14GGG with the parameter set generated in dynamical simulation (a) and the statistical model (b).

TABLE 1: Rate Constant of Hole Transfer (mean value ( standard deviation, ns-1) in GTGGG and GT14GGG with the Various Statistical Modelsa model

sequence

εi

Tij

GTGGG

MD MD dist dist

MD dist MD dist

361 ( 184 417 ( 177 228 ( 116 248 ( 77

GT14GGG 29 ( 19 34 ( 17 8(5 12 ( 5

a Site energies εi and electronic couplings Tij were each either taken from simulations (MD) or drawn randomly from the respective distributions (dist).

clearly that there are no collective motions in DNA which would drive CT via the fluctuations of electronic couplings Tij, that is, there are no regions of local DNA compression with enlarged couplings which would enhance CT. Discussion and Conclusions In this work, we have studied the hole transfer across simple A bridges in double-stranded DNA by means of a coarse-grained TD-DFT approach coupled to classical MD simulations of DNA in a solvent environment. Important ingredients of such simulations are (i) an accurate estimate of the CT parameters εi and Tij, (ii) inclusion of dynamical effects, and (iii) proper treatment of solvation. A study focusing on one of these factors yet neglecting the other ones will miss some of the crucial aspects of the hole transfer mechanism. In the last years, model Hamiltonian approaches have been mostly replaced by bottom-up methods, where the CT parameters are calculated using quantum chemical methods along classical MD trajectories, thereby taking care of the points (i) and (ii). Concerning (i), there is an obvious

Solvent Fluctuations Drive the Hole Transfer in DNA methodological convergence in the field of simulations. Thanks to the ongoing efforts of Voityuk et al., there are plenty of reliable benchmark data for the CT parameters, first at the HF level and later on using the CAS-PT2 method.32,78 In recent simulations, the CT parameters are computed with methods like HF, DFT, SCC-DFTB, INDO/S, or other semiempirical Hamiltonians,45-47,51,57,58 which all reproduce the benchmark data reasonably well. Further, many studies make use of classical MD trajectories, along which the parameters are evaluated. Since the common force fields describe the structure and dynamics of proteins and DNA quite well, no major disagreement is to be expected here. All authors report variations of site energies due to structural fluctuations along the MD trajectories on the order of 0.1-0.15 eV and those of the couplings on the order of 0.05 eV.45-47,51,57,58 Proper inclusion of the solvent effects (iii) is more intricate. Standard QM/MM frameworks are able to take these effects into account though, yet they have been used by only a few authors up to now. The coupling of CT parameters to the solvent degrees of freedom brings on fluctuations of site energies on the order of 0.4 eV.25,46,57,58 Our results show that these large fluctuations are indeed the driving force for the hole transfer, while the effects of structural variance are merely of secondary importance. However, it has to be noted that the QM/MM coupling used in this work still neglects the polarization of solvent by the hole charge. This may be the reason for the weak distance dependence of the CT rate for short A bridges found in our simulations. Therefore, the development of a bottom-up approach including the solvent polarization will be a major task for future work. It is Steinbrecher et al.46 who already made a first step in this direction; yet, they propagated the hole in terms of adiabatic Born-Oppenheimer dynamics. Such an approach may lead to reasonable statistical distributions though but cannot account for the true dynamics of the electronic system whatsoever. Our simulations give further insight into the role of dynamics of both DNA and solvent. Electron-transfer reactions in proteins and solution may be dominated by few fluctuations, as reported by Balabin et al.52 However, this seems to be different in DNA. As we have shown, the fluctuations of Tij are of minor importance compared to those of the site energies. Thus, the variation of DNA structure (or the DNA collective modes) does not drive the CT in the first place. This shall not be interpreted so that the electronic couplings are irrelevant for the transport; as we have demonstrated, there are instances of insufficient coupling of the nucleobases, where CT is effectively suppressed. However, the results of simulations employing a stochastic model of CT parameters indicate that the time course of fluctuations plays no unique CT-promoting role. Therefore, there are no CT-promoting DNA modes, which would enhance CT via the rise of favorable couplings along the DNA, and the original conformational gating model of Barton et al. is not supported by our calculations. Nevertheless, the basic intuition of CT-active conformations is confirmed; it is just necessary to focus on the collective motions that drive the dynamics of site energies. A similar idea was proposed by Hatcher et al.,51 suggesting that the nucleobase skeletal motion drives the fluctuations of site energies. The observation of the internal dynamics of nucleobases affecting the site energies is correct in principle, but this effect is completely overwritten by the large fluctuations introduced by the solvent. The solvent introduces rapid fluctuations with a period of 40 fs, and the CT seems to be determined by largescale changes in the electrostatic potential in the picosecond domain. We find periods on this time scale where the effective

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13115 barrier is dramatically reduced, which explains the very different CT for different trajectories. Since the reduction of the barrier is controlled by the electrostatic potential induced by the solvent, we can conclude that it is the solvent motion that controls the CT in DNA, introducing CT-active conformations. The importance of the factors (i)-(iii) can be discussed conveniently in the framework of Marcus theory (eq 1), as was already done for the exponential distance dependence of the rate of hole transfer in short A bridges (with up to four A’s) in the experiments by Giese27 and Lewis et al.64 The exponential dependence arises due to both terms in eq 1, on the one hand, HDA as the wave function overlap decays exponentially and on the other hand, from the solvent reorganization term owing to the distance dependence of solvent reorganization energy λs, as shown by LeBard et al.76 (∆G0 vanishes in this case as both the donor and the acceptor are guanines.) Therefore, the exponent in the distance dependence of the CT rate

k ∝ exp[-βRDA]

(10)

can be written as a sum

β ) β0 + βs

(11)

where β0 arises from the HDA part and βs originates from the distance dependence of λs.76 β0 was evaluated in several studies using quantum chemical simulations, assuming a static barrier of A’s through which the hole has to tunnel, and the results were in the range of 0.7-1.0 Å-1.32,37,79 Here, the exponential decay rate is a trivial consequence of the tunneling mechanism. βs is more difficult to evaluate; recent studies estimated this exponent at around 1 Å-1.76,80 The combination of these two values, however, leads to a huge overestimation of the decay parameter when compared to experiment, as discussed in detail in ref 76. In the current work, we have focused on the former factor, that is, the effect of the HDA part manifesting itself in β0. Up to now, the hole transfer through short A bridges was described in terms of superexchange tunneling (with decay characterized by the exponent β0). This picture was supported by model calculations based on the CT parameters computed for static structures. Indeed, the G-A energy barrier of about 0.4 eV and the electronic couplings on the order of 0.1 eV result in a decay parameter of around 1 Å-1,32,37,79 in reasonable accordance with the experimental observations.19,81,82 However, this agreement was achieved neglecting the effect of solvent response, represented by βs. Therefore, the agreement with experiment may have resulted from a fortunate error cancellation as both dynamical and solvent effects were neglected. In our work as well as in the mentioned studies based on MD simulations,45,51 a significant hole population was found on the bridge sites, indicating the hopping mechanism to be operative already in short bridges with up to four A’s. Since the interactions with solvent (QM/MM) were disregarded in both previous studies,45,51 the hole population on the bridge was probably even underestimated, as our results suggest. Therefore, it seems very likely that the role of superexchange tunneling is less important than currently assumed. However, the hopping mechanism may lead to a weaker distance dependence of rate, and as a matter of fact, our calculations cannot reproduce the (experimentally observed) exponential distance dependence. An analysis the CT-relevant time scales may help to understand the issue of distance dependence. As discussed above

13116

J. Phys. Chem. B, Vol. 113, No. 39, 2009

and in our previous work,58 the phenomena affecting the CT parameters occur on several different time scales. First, there are fast nucleobase skeletal modes with a period of 20 fs and the water modes discussed above with periods of about 40 fs. Further, we have shown58 that the counterions introduce fluctuations in the picosecond regime but are shielded by water molecules strongly. The mentioned calculations suggested that the electric fields induced by the ions and the water do not cancel completely. Rather, a net total electric potential results, biasing the ionization potentials of nucleobases (site energies) vigorously. This leads to very different electronic configurations (G-A energy barrier), which explains the very different CT efficiency along the individual trajectories. Therefore, the CT is strongly dependent on the conformation, that is, the solvent conformation gates the CT. However, this means that the time scales of ion motions and hole transfer possess comparable magnitudes. In such a case, as discussed in detail by Berlin et al.,83 there are several time scales necessary to characterize the transfer instead of a single one, and the distance dependence of the rate becomes nonexponential (dispersive kinetics). Interestingly, the fluctuations of couplings Tij are very fast compared to the time scale of ion fluctuations. This leads to self-averaging of the couplings,76 meaning that the time series of Tij can be substituted by the average values of the couplings without compromising the accuracy of resulting rates, as demonstrated above as well. Therefore, another main result of this work is that the dynamics of the system and the electrostatic effects of solvent lead to nonexpontial kinetics of hole transfer, and no β0 parameter can be determined for this process. The hole transfer governed by the HDA part is then too fast and nonexponential, that is, it cannot account for the kinetics observed experimentally. This means that the solvent reorganization must be the other key to understand the CT kinetics. According to the analysis in ref 76, the distance dependence of the solvent reorganization energy becomes the determinant for the exponential distance dependence. As a conclusion, our simulations describe the features of hole transfer dynamics arising from the HDA part; however, to render correct transfer rates (and their distance dependence), it is necessary to account for the effect of solvent reorganization on the microscopic scale as well. Acknowledgment. This work was supported by Deutsche Forschungsgemeinschaft (SPP 1243 “Quantum Transport at the Molecular Scale”, Project DFG-EL 206/5-1). Supporting Information Available: Selected data presented as graphics in tabular format. Simulation details. Discussion of the trapping constant τ. Discussion of convergence of chargetransfer rates. Survival of charge in GTGGG if either site energies or electronic couplings are taken constant. Survival of charge in GTGGG without the QM/MM treatment (in vacuo). Correlation coefficients of site energies in all sequences studied. Distributions of CT parameters along MD trajectories. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Hall, D. B.; Holmlin, E.; Barton, J. K. Nature 1996, 382, 731. (2) Hall, D. B.; Barton, J. K. J. Am. Chem. Soc. 1997, 119, 5045. (3) Gasper, S. M.; Schuster, G. B. J. Am. Chem. Soc. 1997, 119, 12762. (4) Giese, B.; Wessely, S.; Spormann, M.; Lindemann, U.; Meggers, E.; Michel-Beyerle, M. E. Angew. Chem., Int. Ed. 1999, 38, 996. (5) Porath, D.; Bezryadin, A.; De Vries, S.; Dekker, C. Nature 2000, 403, 635.

Kubarˇ et al. (6) Xu, B.; Zhang, P.; Li, X.; Tao, N. Nano Lett. 2004, 4, 1105. (7) van Zalinge, H.; Schiffrin, D. J.; Bates, A. D.; Haiss, W.; Ulstrup, J.; Nichols, R. J. ChemPhysChem 2006, 7, 94. (8) Boon, E. M.; Livingston, A. L.; Chmiel, N. H.; David, S. S.; Barton, J. K. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 12543. (9) Holman, M. R.; Ito, T.; Rokita, S. E. J. Am. Chem. Soc. 2007, 129, 6. (10) Long-Range Charge Transfer in DNA I-II: Topics in Current Chemistry; Schuster, G. B., Ed.; Springer: Heidelberg, Germany, 2004. (11) Palecˇek, E. Electroanalysis 1996, 8, 7. (12) Kelley, S. O.; Barton, J. K.; Jackson, N. M.; Hill, M. G. Bioconjugate Chem. 1997, 8, 31. (13) Wong, E. L. S.; Gooding, J. J. Anal. Chem. 2006, 78, 2138. (14) Introducing Molecular Electronics: Lecture Notes in Physics; Cuniberti, G., Fagas, G., Richter, K., Eds.; Springer: Berlin, Germany, 2005. (15) Joachim, C.; Ratner, M. A. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 8801. (16) Hush, N. S.; Cheung, A. S. Chem. Phys. Lett. 1975, 34, 11. (17) Orlov, V. M.; Smirnov, A. N.; Varshavsky, Y. M. Tetrahedron Lett. 1976, 17, 4377. (18) Steenken, S.; Jovanovic, S. V. J. Am. Chem. Soc. 1997, 119, 617. (19) Meggers, E.; Michel-Beyerle, M. E.; Giese, B. J. Am. Chem. Soc. 1998, 120, 12950. (20) Ly, D.; Kan, Y.; Armitage, B.; Schuster, G. B. J. Am. Chem. Soc. 1996, 118, 8747. (21) Jortner, J.; Bixon, M.; Langenbacher, T.; Michel-Beyerle, M. E. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 12759. (22) Grozema, F. C.; Berlin, Y. A.; Siebbeles, L. D. A. J. Am. Chem. Soc. 2000, 122, 10903. (23) Bixon, M.; Jortner, J. J. Phys. Chem. B 2000, 104, 3906. (24) Berlin, Y. A.; Burin, A. L.; Ratner, M. A. J. Am. Chem. Soc. 2001, 123, 260. (25) Voityuk, A.; Siriwong, K.; Roesch, N. Angew. Chem., Int. Ed. 2004, 43, 624. (26) Giese, B.; Amaudrut, J.; Ko¨hler, A. K.; Spormann, M.; Wessely, S. Nature 2002, 412, 318. (27) Giese, B. Top. Curr. Chem. 2004, 236, 27. (28) Bixon, M.; Jortner, J. J. Am. Chem. Soc. 2001, 123, 12556. (29) Berlin, Y. A.; Burin, A. L.; Ratner, M. A. Chem. Phys. 2002, 275, 61. (30) Bixon, M.; Jortner, J. Chem. Phys. 2002, 281, 393. (31) Bixon, M.; Jortner, J. Chem. Phys. 2006, 326, 252. (32) Roesch, N.; Voityuk, A. Top. Curr. Chem. 2004, 237, 37. (33) Marcus, R. A. J. Chem. Phys. 1956, 24, 966. (34) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265. (35) Priyadarshy, S.; Beratan, D. N.; Risser, S. M. Int. J. Quantum Chem. 1996, 60, 1789. (36) Tong, G. S. M.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. B 2002, 106, 2381. (37) Olofsson, J.; Larsson, S. J. Phys. Chem. B 2001, 105, 10398. (38) Cramer, T.; Krapf, S.; Koslowski, T. J. Phys. Chem. B 2004, 108, 11812. (39) Renger, T.; Marcus, R. A. J. Phys. Chem. A 2003, 107, 8404. (40) Jakobsson, M.; Stafstro¨m, S. J. Chem. Phys. 2008, 129, 125102. (41) Voityuk, A.; Siriwong, K.; Roesch, N. Phys. Chem. Chem. Phys. 2001, 3, 5421. (42) Troisi, A.; Orlandi, G. J. Phys. Chem. A 2002, 106, 2093. (43) Cramer, T.; Steinbrecher, T.; Labahn, A.; Koslowski, T. Phys. Chem. Chem. Phys. 2004, 7, 4039. (44) Voityuk, A. Chem. Phys. Lett. 2006, 439, 162. (45) Grozema, F. C.; Tonzani, S.; Berlin, Y. A.; Schatz, G. C.; Siebbeles, L. D. A.; Ratner, M. A. J. Am. Chem. Soc. 2008, 130, 5157. (46) Steinbrecher, T.; Koslowski, T.; Case, D. A. J. Phys. Chem. B 2008, 112, 16935. (47) Rˇeha, D.; Barford, W.; Harris, S. Phys. Chem. Chem. Phys. 2008, 10, 5436. (48) Siriwong, K.; Voityuk, A. A. J. Phys. Chem. B 2008, 112, 8181. (49) Voityuk, A. A. J. Chem. Phys. 2008, 128, 045104. (50) Voityuk, A. A. J. Chem. Phys. 2008, 128, 115101. (51) Hatcher, E.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Beratan, D. N. J. Am. Chem. Soc. 2008, 130, 11752. (52) Balabin, I. A.; Beratan, D. N.; Skourtis, S. S. Phys. ReV. Lett. 2008, 101, 158102. (53) O’Neill, M. A.; Barton, J. K. J. Am. Chem. Soc. 2004, 126, 11471. (54) O’Neill, M. A.; Barton, J. K. Top. Curr. Chem. 2004, 236, 67. (55) Schuster, G. B.; Landmann, U. Top. Curr. Chem. 2004, 236, 139. (56) Barnett, R. N.; Cleveland, C. L.; Joy, A.; Landman, U.; Schuster, G. B. Science 2001, 294, 567. (57) Kubarˇ, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. J. Phys. Chem. B 2008, 112, 7937. (58) Kubarˇ, T.; Elstner, M. J. Phys. Chem. B 2008, 112, 8788. (59) Mantz, Y. A.; Gervasio, F. L.; Laino, T.; Parrinello, M. Phys. ReV. Lett. 2007, 99, 058104.

Solvent Fluctuations Drive the Hole Transfer in DNA (60) Grozema, F. C.; Siebbeles, L. D. A.; Berlin, Y. A.; Ratner, M. A. ChemPhysChem 2002, 3, 536. (61) Senthilkumar, K.; Grozema, F. C.; Guerra, C. F.; Bickelhaupt, F. M.; Lewis, F. D.; Berlin, Y. A.; Ratner, M. A.; Siebbeles, L. D. A. J. Am. Chem. Soc. 2005, 127, 14894. (62) Senn, H. M.; Thiel, W. Curr. Opin. Chem. Biol. 2007, 11, 182. (63) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. ReV. B 1998, 58, 7260. (64) Lewis, F. D.; Zhu, H.; Daublain, P.; Fiebig, T.; Raytchev, M.; Wang, Q.; Shafirovich, V. J. Am. Chem. Soc. 2006, 128, 791. (65) Endres, R. G.; Cox, D. L.; Singh, R. R. P. ReV. Mod. Phys. 2004, 76, 195. (66) Loewdin, P. O. J. Chem. Phys. 1950, 18, 365. (67) Rak, J.; Makowska, J.; Voityuk, A. Chem. Phys. 2006, 325, 567. (68) Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. J. Phys. Chem. B 2001, 105, 569. (69) Runge, E.; Gross, E. K. U. Phys. ReV. Lett. 1984, 52, 997. (70) Niehaus, T. A.; Heringer, D.; Torralva, B.; Frauenheim, T. Eur. Phys. J. D 2005, 35, 467. (71) Brankin, R. W.; Gladwell, I.; Shampine, L. F. RKSUITE: A suite of Runge-Kutta codes for the initial value problem for ODEs, Softreport 92-S1; Department of Mathematics, Southern Methodist University: Dallas, TX, 1992. (72) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049.

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13117 (73) Pe´rez, A.; Marcha´n, I.; Svozil, D.; Sˇponer, J.; Cheatham, T. E., III; Laughton, C. A.; Orozco, M. Biophys. J. 2007, 92, 3817. (74) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668. (75) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (76) LeBard, D. N.; Lilichenko, M.; Matyushov, D. V.; Berlin, Y. A.; Ratner, M. A. J. Phys. Chem. B 2003, 107, 14509. (77) Kubarˇ, T.; Elstner, M. J. Phys. Chem. B 2009, 113, 5653. (78) Blancafort, L.; Voityuk, A. J. Phys. Chem. A 2005, 110, 6426. (79) Jortner, J.; Bixon, M.; Voityuk, A. A.; Roesch, N. J. Phys. Chem. A 2002, 106, 7599. (80) Siriwong, K.; Voityuk, A. A.; Newton, M. D.; Ro¨sch, N. J. Phys. Chem. B 2003, 107, 2595. (81) Lewis, F. D.; Kalgutkar, R. S.; Wu, Y.; Liu, X.; Liu, J.; Hayes, R. T.; Miller, S. E.; Wasielewski, M. R. J. Am. Chem. Soc. 2000, 122, 12346. (82) Lewis, F. D.; Liu, J.; Weigel, W.; Rettig, W.; Kurnikov, I. V.; Beratan, D. N. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12536. (83) Berlin, Y. A.; Grozema, F. C.; Siebbeles, L. D. A.; Ratner, M. A. J. Phys. Chem. C 2008, 112, 10988.

JP9073587