Microsecond Simulation of Electron Transfer in DNA: Bottom-Up

Subscriber access provided by UNIV OF REGINA

Article

Microsecond Simulation of Electron Transfer in DNA: Bottom-Up Parametrization of an Efficient Electron Transfer Model Based on Atomistic Details Mario Wolter, Marcus Elstner, Ulrich Kleinekathöfer, and Tomas Kubar J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b11384 • Publication Date (Web): 03 Jan 2017 Downloaded from http://pubs.acs.org on January 7, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Microsecond Simulation of Electron Transfer in DNA: Bottom-up Parametrization of an Efficient Electron Transfer Model Based on Atomistic Details Mario Wolter,†,§ Marcus Elstner,† Ulrich Kleinekathöfer,‡ and Tomáˇs Kubaˇr∗,¶ Institute of Physical Chemistry, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany, Department of Physics and Earth Sciences, Jacobs University Bremen, 28759 Bremen, Germany, and Center for Functional Nanostructures & Institute of Physical Chemistry, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany E-mail: [email protected]

∗

To whom correspondence should be addressed IPC, Karlsruhe Institute of Technology ‡ Jacobs University Bremen ¶ CFN & IPC, Karlsruhe Institute of Technology § Present address: Institute of Physical and Theoretical Chemistry, Technische Universität Braunschweig, 38106 Braunschweig, Germany †

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 62

Abstract The transfer of electrons over long distances in complex molecular systems is a phenomenon of significance in both biochemistry and technology. In recent years, we have been developing efficient models to study ET in complex systems, including DNA as a prominent example. Ab initio and model approaches have been combined in an ‘on the fly’ calculation of ET parameters, which can be used to propagate nuclear and electronic degrees of freedom simultaneously. These previous efforts have aimed at deriving an efficient non-adiabatic QM/MM simulation scheme for ET, making nanosecond simulations of ET in realistic systems possible. This, however, is still insufficient for the treatment of large donor-bridge-acceptor systems, like the ET in DNA overcoming long adenine bridges. Therefore, we have constructed a theoretical model in a bottom-up way.

All

quantum-chemical as well as force-field calculations are substituted by theoretical models of the involved phenomena on a molecular level, including the polarization and relaxation of the molecular environment, which is often omitted in other recently developed theoretical models of ET. A non-adiabatic simulation scheme is employed and no assumptions regarding the ET mechanism are needed. Thus, the predictive power of the simulations is preserved, while pushing the limits of the accessible time scales beyond microseconds. This model-based simulation scheme is applied to ET in various DNA species. Good agreement with the ‘full’ atomistic non-adiabatic QM/MM scheme is observed for the archetypal DNA ET systems, the polyA sequence as well as the sequences GTn GGG containing adenines as bridge sites. Furthermore, ET in larger, more complex DNA sequences is simulated, and the results are discussed.

1

Introduction

Long-range electron transfer (ET) in DNA has been in the focus of research interest for at least 25 years by now. Various experimental approaches have been applied, and among 2 ACS Paragon Plus Environment

Page 3 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the key experimental investigators, at least the following should be mentioned: Barton et al., 1,2 Schuster et al., 3 Giese et al., 4 Lewis et al. 5 and Porath et al. 6 These experimental findings and potential applications stimulated a large body of experimental and theoretical work in order to elucidate the mechanism of ET. 7–9 Two principal physical mechanisms were proposed for DNA ET, one-step superexchange-mediated tunneling over shorter distances and thermally induced hopping (TIH) for longer distances. 10–14 Several hypotheses were then brought up regarding the molecular mechanisms driving the ET process, taking into account the dynamics of the molecular system – this class of models include the conformational gating, 15,16 the ion-gated phonon-assisted hopping, 3 or the polaron drift. 17 Rather recently, time-resolved spectroscopic experiments have made it possible to follow the kinetics of ET reactions in DNA more directly. 18,19 The considerable attention they have received confirms the unrelenting interest in the topic, and they have motivated continuing theoretical work aimed to obtain a clearer picture of the mechanisms of ET. Ratner et al. 20 set up a stochastic surrogate Hamiltonian and propagated the density matrix of the hole using a Liouville equation. They found a unified mechanism of ET that involves a (partially) delocalized hole state on the (adenine) bridge, rather than a discrete hopping from one bridge site to another. This model was also used, e.g., to interpret the dependence of the ET rates on the sequence of DNA species where CG base pairs were substituted by AT. 21 Beratan et al. 22 proposed a mechanism based on a transient, ‘flickering resonance’ of energy levels of multiple charge-carrying sites facilitating hole transfer, possibly accompanied by a delocalization of the hole. Moreover, Zhang et al. 23 presented a model Hamiltonian with emphasis on spatial and temporal correlations of site energies and couplings and propagated the wave function of the hole using a stochastic Schrödinger equation (SSE). Utilizing this stochastic approach, Tao et al. 24 interpreted the outcome of their experiments on the conductivity of DNA, also involving a delocalized hole. 25 The current state of research has been reviewed, 26 while yet additional proposals for ET mechanisms appear like the ‘quantum mechanical unfurling’. 27 In general, such theoretical models provide a reduced description of the physical sys-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tem. The virtue of this is that irrelevant degrees of freedom are integrated out, and the process is explained by focusing on the important degrees of freedom, which are represented by parameters in the model. Theoretical models are indispensable for the rationalization of experimental results, supplementing details not accessible by experiments, predicting new properties and possibly even supporting the design of new molecules. When developing model Hamiltonians, however, there always is a danger of missing an important component. An issue common to all the above mentioned theoretical investigations is that the polarization and relaxation of the molecular environment, prominently the aqueous solvent, is not considered, and its effect on the energetics of ET is not included. Note that this relaxation manifests itself in the outer-sphere reorganization energy, to use the language of Marcus’ theory of ET, 28,29 which exceeds in magnitude both the differences of site energies and their fluctuations in DNA. In a brief discussion on relaxation in one of the above mentioned studies, 20 it was argued that the relaxation is much slower than the transfer itself, therefore it may be neglected. This statement is debatable, as may be seen from the comparison of the ET kinetics reported in that paper and the relaxation kinetics in DNA, e.g., below in this paper. The major part of relaxation occurs actually on a sub-picosecond time scale and thus on a time scale comparable to that of ET. An unbiased description of ET can be achieved when employing ab initio non-adiabatic dynamics techniques. 30 In an all-atom picture, the time-dependent Schrödinger equation for the electrons is solved, considering a classical description of the nuclei. This approach provides high accuracy, however such simulations are limited to tiny spatial and temporal scales. To pass to larger systems and slower processes, we derived an efficient simulation scheme for ET, using the idea of fragmentation of electronic structure, the QM/MM formalism and applying further well-controlled approximations. 31–34 Nanosecond simulations of ET in realistic systems can be performed routinely, still avoiding any bias, e.g. in the form of assumptions on the mechanism of ET. The advantage of this approach may be illustrated by comparing to the original version of Marcus’ theory, which involves several assumptions


Page 4 of 62

Page 5 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


regarding the localization of wave function, the separation of electronic and ionic time scales, and the degree of adiabaticity of the process. Our simulation approach proved to be a viable alternative, e.g., in the study of ET in the E. coli DNA photolyase; 35 the reaction consists of two consecutive ET steps. Since the time scales of ET and of outer-sphere relaxations overlap, Marcus’ theory underestimates the rate of the second step, while the multi-scale simulation scheme reproduces it fairly well. The idea of the present work is to combine the advantages of the non-adiabatic QM/MM scheme and of the model approaches: the unbiased description of ET of the former, and the unbeatable computational efficiency of the latter. Contrary to previous approaches, 20,22,23 which were based on phenomenological models and performed simulations in such frameworks, we will develop a theoretical simulation scheme in a bottom-up way: Our existing multi-scale non-adiabatic framework will be used as a basis, into which our quantitative knowledge of the ET system will be introduced. In the end, all of the quantum-chemical as well as force-field calculations will be substituted by theoretical models of the involved phenomena on a molecular level. Still, no assumptions will be needed regarding the mechanism of ET itself, so that the predictive power of the simulation will be preserved. Such a development is both natural, as it is necessary to simulate for very long time scales, and unique among the DNA ET research published so far, to our knowledge. Further, while the current work concentrates on ET in dsDNA, it will be possible to extend the model to other DNA systems like, e.g., the recently reported three-way junctions, 36 and to create similar schemes for ET in other molecules or materials. Thus, this work will start with an analysis of the electronic structure of DNA, followed by the parametrization of the fragment-based Hamiltonian; this step is quite similar to previous approaches by other laboratories. 20,22,23 Then, the electronic and nuclear relaxation effects, as well as the time-dependence of the latter will be taken from our multi-scale dynamics scheme partly. This component is new in comparison with previous approaches. Finally, a theoretical simulation scheme based on a non-adiabatic propagation method is developed. It



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

will turn out that three main components are necessary to construct the desired theoretical model: (i) A model Hamiltonian for the electronic structure of a neutral DNA molecule. This point was addressed by many researchers in the past. (ii) A model of the polarization response of the environment on the presence of excess charge (electron). This property was also described previously in some studies; interestingly, not all researchers in the field seem to have recognized its significance. (iii) The dynamics of the polarization response, or the relaxation of environment upon structural changes of the environment. The newly developed model-based simulation scheme will be applied to several archetypal DNA ET systems.

2

Model

2.1

Initial considerations

In this work, a model of solvated DNA is being developed, such that the Hamiltonian, the electron–phonon coupling and the coupling to the bath are obtained from analysis of explicit all-atom QM/MM simulations of solvated DNA. This route of parametrization guarantees all of the important modes of motion of both DNA and its molecular environment (solvent) to be taken into account. These data are being fed into a propagation scheme, similar to the original one based on QM/MM simulation. Thus, the explicit time-dependent character of the simulation is kept even with the theoretical model. Importantly, the model allows to simulate ET in DNA in an unbiased fashion, i.e. the questions of promoting modes, time scale separation, adiabaticity and delocalization are not assumed, rather, they may be answered. The developed model allows for efficient simulations of long DNA oligonucleotides, which are not possible with explicit QM/MM schemes. Notably, simulation of ET in complex molecular systems relies on the sampling of rare events with correct probabilities, and these will obtained as long as the parameters of the model have been determined accurately. Thus, the model simulation scheme will provide correct ET rates in principle, even if the simulated time scales exceed those of the all-atom QM/MM simulations used previously for 6 ACS Paragon Plus Environment

Page 6 of 62

Page 7 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the parametrization, as long as those QM/MM simulations were long enough to provide converged values of the parameters. Our current work makes use of this important point. An additional benefit is the model-based understanding of the ET process generated by the reduction of degrees of freedom. A decomposition into the different, conceptually simple contributions of the model Hamiltonian provides insight into the molecular modes that are the driving forces of ET in DNA. Such knowledge may not be easy to obtain from full ab initio or DFT calculations. Further, the derivation of the model outlines a general scheme for an analysis of ET in complex materials, which can be applied to other systems like ET in organic materials, where also the question of de/localization of electron is paramount. Thus, the concept of this work is as follows: To obtain a time-dependent picture of the ET process, coupled to the modes of motion of DNA or the molecular environment, an explicitly time-dependent Hamiltonian, or total energy expression is to be derived on the basis of explicit all-atom QM/MM simulations of solvated DNA. The analysis determines the most important vibrational modes of DNA and its environment, from which an explicitly time-dependent total energy expression is built then. Finally, the equations of motion for the electronic and ionic system are derived from the total energy expression.

2.2

Electronic structure of DNA

ET in DNA occurs along electronic states supported by the nucleobases. The corresponding electronic structure may be computed in different ways: (i) With a standard linear combination of atomic orbitals (LCAO) approach using any quantum chemical method like the Hartree–Fock (HF) or the density functional theory (DFT), the orbitals of a ‘super-molecule’ consisting of several stacked nucleobases may be calculated with localized atomic orbitals constituting a basis set. The highest occupied molecular orbital (HOMO, ΨHOMO ) resulting from such a calculation is shown in Fig. 1.a (ii) Alternatively, the orbitals ϕm of the individual nucleobases (also ‘sites’) are computed in a first step, for the ‘super-molecular’ orbitals a

Note that such extended (delocalized) orbitals are obtained for static, idealized double-helical DNA structures. Static and dynamic disorder as well as polarizable solvent exhibit localizing effects on the orbitals.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 62

Figure 1: HOMO of a polyA DNA oligonucleotide in an idealized B-DNA geometry, obtained on the DFT level TPSS/def2-SVP. to be obtained as linear combinations of ϕm ,

Ψi =

X

aim ϕm .

(1)

m

This approach is suggested visually by the appearance of this HOMO, which is formed by a combination of HOMOs of the individual nucleobases, ϕm . Importantly, it can be exploited quantitatively by using a linear-scaling ‘fragment molecular orbital’ (FMO) based method. Then, the Hamilton matrix elements Hmn are obtained using ϕm , E ˆ ϕm H ϕm , D E ˆ = ϕm H ϕn ,

Hmm ≡ εm = Hmn ≡ Tmn

D

(2) (3)

and the expansion coefficients aim may be calculated with a diagonalization of this Hamiltonian. The results of this two-step procedure can be compared to those of a direct calculations of Ψi from AOs, as was illustrated in Ref. 34 (Fig. 2 in that publication). The energies of Ψi were computed in both ways, and the good agreement of the shifts and splitting of orbital energies due to the hybridization of ϕm confirmed the approximation in Eq. 1 to be valid.b b

This idea was explored in the early work on electron tunneling through organic molecules 37–39 showing that an efficient method for large systems can be based on the calculation of matrix elements for smaller molecules, as long as the chemical environment is comparable. Further, Kurnikov and Beratan 40 elaborated the idea of fragment Hamiltonian, partitioning a large molecule into overlapping fragments, from which the Hamiltonian may be built in a computationally efficient way.


Page 9 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Such a simple H¨ uckel-type Hamiltonian can be used in cases where the support of the wave function results from an assembly of weakly coupled fragments, like in DNA or organic semi-conductors. In these cases, the coupling between the fragments as given by Tmn is weak, typically below 0.1 eV, implying that the HOMOs of the fragments are only weakly disturbed by the presence of the neighboring fragments. Then, the weak coupling of these orbitals results in the formation of the mini-bands. Many various ways to compute εm and Tmn with quantum chemical methods have been developed. 41,42 In our previous work, we have used the semi-empirical DFTB method, 33,34,43 which provides the matrix elements accurately and efficiently. It has been reported recently that matrix elements of excellent quality are obtained with a uniform scaling of DFTB results by a constant factor of 1.54; 42,44 this approach is employed in the current work as well. Last but not least, the electronic structure responds sensitively to the static and dynamic disorders as well as solvation effects. Both factors can be estimated when computing the ET parameters εm and Tmn along classical molecular dynamics (MD) trajectories, including the molecular environment using combined QM/MM approaches.

2.2.1

Site energies

The site energies εn determine the energy landscape of the ET reaction. They can be identified with the ionization potential (IP) in case of hole transfer (hole, missing electron) and with the electron affinity (EA) in case of excess electron transfer based on Koopmans’ theorem.c In the case of DNA hole transfer, the HOMOs of guanine (G) and adenine (A) bases are the only relevant orbitals, as illustrated in Fig. 2, because these nucleobases are the parts of DNA with the lowest IP. This is an entirely static picture of ET, however. As soon as the thermal motion of DNA and the aqueous solvent is accounted for, the site energies c

Note that only Janak’s theorem is available in the framework of DFT, allowing to identify the HOMO energy with the IP, as long as an exact exchange–correlation functional is applied. Since approximate XC functionals are applied practically, εn often deviates from the IP. Fortunately, it is sufficient if the site energy difference is a good estimate of the IP difference. This has to be tested for every system prior to application, 43 or can be corrected using a suitable reference. 45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

start to fluctuate.

Figure 2: A simple idea of hole transfer in DNA, superexchange tunneling (red) and thermally induced hopping (blue).

The QM/MM procedure combined with classical MD simulation makes it possible to separate different sources of site energy fluctuation: 46 (i) when the QM/MM terms are switched off, the environment (solvent and DNA backbones) does not affect the electronic structure of the nucleobase, which is considered in the ‘gas-phase’ effectively; still, the internal motion of the nucleobase leads to fluctuations of εn . (ii) with QM/MM for the environment, the complete fluctuation pattern of εn is obtained, considering both the internal dynamics of the nucleobase and the dynamics of the environment.

Figure 3: IP of an adenine embedded in a DNA oligonucleotide, calculated using the DFTB method along a classical MD simulation. Distribution of IP obtained with (black) as well as without (red) the QM/MM implementation covering the effect of the molecular environment. The obtained distributions of site energies are Gaussian, and the standard deviation is 0.10 eV and 0.30 eV, in the gas-phase and QM/MM calculations, respectively, 46 as shown 10 ACS Paragon Plus Environment

Page 10 of 62

Page 11 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Fig. 3. In the gas-phase calculations, only the internal vibrational modes of the nucleobase are considered, and they lead to rather small fluctuations of site energies. The power spectrum of site energy in Fig. 4 displays the strongest signal around 1800 cm−1 , and the responsible modes of motion can be identified with the stretching modes of the conjugated double bonds. 47 The QM/MM calculations reveal the critical importance of the solvent, with the amplitude of fluctuations increased by a factor of three. The most important modes of motion are in the range below 1000 cm−1 , see Fig. 4. Therefore, the internal degrees of freedom of the nucleobase and the solvent modes are located in quite distinct regions of the power spectra, a finding also reported by other groups, 48–50 and they are statistically uncorrelated (shown in the Supporting Information).

Figure 4: Power spectra (Fourier transform of the autocorrelation function) of the time series of site energies (IP) obtained for the gas-phase as well as QM/MM data. So, the internal modes of nucleobases are mostly above 1000 cm−1 , and they contribute a small part of the fluctuations of site energies. On the other hand, the modes that are attributed to the environment (solvent) are below 1000 cm−1 , and they dominate the fluctuations of site energies. This finding leads to the impression that the solvent modes of motion are the driving forces for the ET, and the internal modes of the nucleobases are of minor importance. Site energy fluctuations of different sites, which are due to solvent effects can, however, be rather strongly correlated. Thus, while there are large fluctuations, the 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 62

difference of neighboring site energies may not fluctuate so strongly, and the fluctuations due to internal modes of motion, which are uncorrelated between the individual sites, may play a non-negligible role. Still, this issue requires further investigations. An effect of conceivable potential to gate the ET is the correlation of IP of close nucleobases. Reported in our earlier work for DNA sequences composed of identical nucleobases, 46 such a correlation exists in dsDNA with non-homogeneous sequence as well, see Tab. 1. Site– Table 1: Pearson correlation coefficients of the time series of the IP in different DNA sequences. 1–n are the ordinal numbers of purine nucleobases in the particular sequence; for example, 1–4 is the correlation of site energies of the first and fourth nucleobase. ‘Dickerson’ denotes the palindromic dsDNA oligonucleotide (CGCGAATTCGCG)2 , PDB ID 1BNA. A8 (GA)4 Dickerson

1–1 1.00 1.00 1.00

1–2 0.57 0.57 0.42

1–3 0.34 0.31 0.31

1–4 0.22 0.21 0.20

1–5 0.14 0.12 0.11

1–6 0.10 0.09 0.05

1–7 0.03 0.05 0.01

1–8 0.00 0.00 0.02

site correlations have been considered in the recent work by Newton, 51 in the framework of Hopfield’s ET model, 52 and their roles were discussed in a different model work by Zhang et al. 23 These correlations mean that the fluctuations of site energies of neighboring sites are quite strongly concerted rather than independent. This phenomenon may have a large impact on the ET efficiency because the site energy difference between neighboring sites determines the energy profile of ET. If neighboring sites fluctuated independently, random barriers would result; however, this does not seem to be the case. Also, the presence of site energy correlations was shown to impact the tunneling through DNA to some degree, however bringing no large-scale changes to the properties of transport. 53

2.2.2

Modeling of site energies

On the basis of the presented extensive analysis, it is possible to derive a model of the time series of site energies. These are represented by a sum of cosine functions with characteristic periods τi and Ti and two different amplitudes for the fluctuations due to the internal modes


Page 13 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of motion (a) as well as the environmental fluctuations (A):

ε∗m (t)

=

ε0m

+

internal X i

environ. X 2π(t + te + m · tc ) 2π(t + tm ) + A · cos a · cos τi Ti i

(4)

The parameter ε0m is chosen to reproduce the desired mean value. The first sum represents the internal modes of the nucleobases, with periods τi . Two distinct sets of periods were created, one for each of A and G nucleobases. For this, the most prominent peaks in the region above 1000 cm−1 of the power spectrum were identified and their periods used as τi ; the number of different internal modes is 14 and 15 for A and G, respectively. The amplitude a was determined so that the standard deviation of the obtained data series fits the result from gas-phase calculations performed along a MD simulation (0.10 eV). A possibility to extend the model of site energies would be to assign every mode i a separate amplitude ai in the summation in Eq. 4, which might be determined from a Fourier transform as in Fig. 4. For the second sum, representing the modes of the environment, a set of 15 periods Ti was chosen to represent the broad band below 1000 cm−1 in the power spectra; note that there are no prominent ‘peaks’, rather a band of largely overlapping environmental fluctuations. The assignment of periods τi and Ti is illustrated in Fig. 5. Since the periods Ti represent the environmental modes and should not depend on the identity of the nucleobase, one common set of periods was developed.d The amplitude A was fitted to reproduce the standard deviation of 0.30 eV yielded by the QM/MM calculations. Further, it has to be made sure that the time series of site energies is not the same in every simulation or even for every nucleobase. To achieve this, the time series are initialized with a randomly generated starting time. The individual nucleobases shall fluctuate independently as their internal dynamics are independent of each other, thus a separate starting time is Clearly, the peak at ca. 600 cm−1 as well as a certain amount of noise below 1000 cm−1 is due to the internal dynamics of the nucleobase. This might explain a fair portion of the difference of the environmental bands of adenine and guanine. As of now, this observation is not implemented in the ET model, and only the cosine functions representing the environment contribute to fluctuations below 1000 cm−1 . d



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5: Assignment of the characteristic frequencies of the ET model (blue bars) on the basis of power spectra belonging to times series of site energies obtained from QM/MM simulations (black). generated for each of them (tm in Eq. 4) randomly.e By contrast, a common starting time (te ) is generated for the fluctuations of the environment randomly. When the environmental fluctuations, which are the same for every nucleobase in the system, are switched on, they affect the fluctuations of the nucleobases and their fluctuations become correlated. This way, the environmental fluctuations would be identical on all nucleobases, and the site energies of all nucleobases would be overly correlated. In order to obtain the right amount of correlation of site energies of neighboring nucleobases from the model, a phase shift of the environmental fluctuations between two neighboring nucleobases (tc ) is introduced. By virtue of the term m·tc in Eq. 4, the difference of phase (time) of each pair of neighbors equals tc . A crucial parameter to monitor here is the distribution of the site energy difference between the neighboring nucleobases, which plays an important role in ET and is closely related to the correlation of the site energies. The standard deviations of site energy differences expresses the width of the respective distributions and are shown in Tab. 2 for a range of values of tc . With the choice of tc = 9.5 fs, the width of these distributions reproduces the data from e

Actually, a separate starting time could be generate for each of the different cosine functions on one site to make them perfectly independent. Any undesired correlations between the different cosine functions are unlikely to appear, however, due to their different frequencies that are not multiples of each other. Therefore, the simplified representation with one common randomly chosen phase is retained.


Page 14 of 62

Page 15 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


QM/MM calculations based on MD simulations. Table 2: Standard deviation of the difference of site energies of neighboring nucleobases as a measure of correlation of these site energies. Data obtained from theoretical simulations of the AAAA oligonucleotide using different choices of the phase shift parameter tc in Eq. 4, as well as from a QM/MM simulation of 1 ns. tc / fs 5.5 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 13.5 QM/MM

std. dev. of εm − εm−1 / eV 0.173 0.192 0.197 0.202 0.207 0.213 0.218 0.223 0.228 0.233 0.251 0.214

The numerical values of all of the parameters in Eq. 4 are presented in the Supporting Information. An example of the time course of site energies as well as the distributions of IP of A and G obtained with the theoretical model in comparison to data from MD simulations are shown in Fig. 6. The obtained distributions of site energies are Gaussian, and the data from the theoretical model agree with those from QM/MM simulations. Note that neither of these data involves any scaling, which would be appropriate to compensate for the missing description of electronic polarization of the environment, as discussed in detail in Sec. 2.4.3. In the applications, a scaling by a factor of 1/1.4 is applied to the cosine functions representing the environmental fluctuations, while those describing the internal dynamics of the nucleobases remain untouched. Let us also look at the spatial and temporal correlation properties of the site energies generated with Eq. 4. (i) The correlation coefficients of site energies of close nucleobases, up to the third neighbor are shown in Fig. 7 (left) for a range of values of the phase shift tc tested in the parametrization. The correlations are too strong in the theoretical model with a small 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6: Left: Visual comparison of site energy of an adenine embedded in a DNA oligonucleotide, obtained from a QM/MM MD simulation (black) and with the theoretical ET model (blue). Right: Distribution of IP of adenine and guanine evaluated with QM/MM along an MD simulation and yielded by the theoretical model. tc , and they get weaker as tc increases. Before the correlation coefficients approach the values from the QM/MM simulations, the correlations between third neighbors start rising again at ca. tc = 9 fs. From this point of view, the final choice of tc = 9.5 fs represents the best compromise between two requirements – the correlation coefficients being as close as possible to the QM/MM values, and the third-neighbor correlation being as weak as possible. (ii) The autocorrelation functions (ACF) of the time series of site energies generated with the theoretical model are compared to those obtained from the QM/MM data in Fig. 7 (right). An important observation is that the ACF decays very fast, with a decay time around 0.3–0.4 ps. The fastest part of autocorrelation from the QM/MM simulations, including the minimum at ca. 20 fs is reproduced quantitatively by the the theoretical model. The course of the ACF on longer time scales is not, however, and the autocorrelation decays too fast and exhibits too large oscillations that do not decay. This is due to the way the time series of site energies are built: the series of cosine functions are constructed from significant frequencies, and the correlations between sites are introduced quite easily. The temporal correlation behavior, however, cannot be reproduced but on the shortest time scales. Now, why should the site energies be modeled with a time series in Eq. 4 and not with a statistical


Page 16 of 62

Page 17 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


approach like the one used for electronic couplings (see below), which would reproduce the desired autocorrelation behavior? The reasoning is that the cosine series is a handy way to introduce the strongest effects observed, which are the dominant frequencies in the signal as well as the correlations between sites. The autocorrelation of site energies pertaining merely for few hundreds of femtoseconds is considered less important and not introduced explicitly.

Figure 7: Spatial and temporal correlations of site energies with the theoretical model compared to data from QM/MM simulations. Left: Correlation coefficients of site energies of nucleobases up to third neighbors; data from the theoretical model using different choices of the parameter tc (crosses) and data from QM/MM simulations (dashed lines). Right: Autocorrelation function of site energies from the theoretical model (green) and from QM/MM simulations (other colors); these decays on a longer interval up to 20 ps (inset, logarithmic time scale).

2.2.3

Electronic couplings

An example of time series of Tmn is shown in Fig. 8. Unlike the fluctuations of site energies εm , which are driven by the environment (i.e. the QM/MM contribution, as discussed above), the variation of Tmn depends exclusively on the relative orientation of the nucleobases, and the environment has little effect. 46 Interestingly, the electronic coupling (EC) passes through a minimum value of nearly zero at the DNA equilibrium geometry (with a twist angle of 36◦ ). 43,54 While the specific form of the dependence of Tmn on the orientation of nucleobases is 2 very complex, it is absolutely clear that even rather high average values hTmn i (as obtained

along MD trajectories) are determined by structural fluctuations out of this equilibrium 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 62

structure. Therefore, it shall be emphasized that the dynamics are of great importance and any static picture is misleading. This topic was discussed in detail in our earlier work, 46 and a similar finding was reported for certain classes of proteins. 55–58

Figure 8: Left: Time series of EC calculated for 5 adenine-adenine base-pair steps for 600 fs. Right: Distributions of EC in homogeneous DNA sequences polyA and polyG calculated along MD trajectories of 100 ps. In detail, the distributions of EC are Gaussian and centered around non-zero values, see Fig. 8. They are distinctly different for the different base-pair steps, with the mean values spanning an order of magnitude; 46 the mean values and variances of EC obtained are presented in the Supporting Information. The correlations between EC of neighboring base-pair steps are negligible, see the Supporting Information. Similarly indistinctive are the power spectra, which feature a broad band that rises continually with decreasing frequency, while any specific structures at higher frequencies are weak.

2.2.4

Modeling of electronic couplings

The dominance of the low-frequency band makes it possible to describe the time dependence of EC with the results of an autocorrelation analysis. To this end, the ACF was calculated from the recorded time series of Tmn as

ACF(τ ) = hTmn (t) · Tmn (τ + t)i , 18 ACS Paragon Plus Environment

(5)

Page 19 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and Fig. 9 shows the obtained ACF for the couplings in a polyA DNA sequence. The obtained ACFs feature a tri-exponential decay. Fitting a tri-exponential function of the P form ACF = 3i=1 Ai · exp [−t/τi ] leads to the following amplitudes Ai and decay times τi : 41 % at 17.3 ps, 13 % at 0.43 ps, 46 % at 0.057 ps.

Figure 9: Autocorrelation functions of the EC of three pairs of neighboring adenines in a polyA sequence (black). Autocorrelation function of an EC time series created by the currently developed method (orange). Detailed view of fast decays up to 12 ps (inset). Now, the time series of EC within the theoretical model are created by drawing random values from a normal distribution with the desired parameters and the desired autocorrelation behavior. This approach seems to be more sensible than setting up a series of cosine functions because the broad autocorrelation trace actually dominated the power spectra over any specific frequencies, as shown above. The current implementation makes use of the Ornstein–Uhlenbeck process (OUP), 59 which generates exponentially time-correlated normally distributed colored noise. An OUP produces data series with an ACF of the form

ACF(τ ) = σ 2 · exp [−λτ ] ,

(6)

with σ 2 being the variance of the normal distribution and λ standing for the inverse correlation time. The details of the whole procedure are given in Ref. 60. Since a tri-exponential auto-correlation function is needed, the time series of an EC within the theoretical model is



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 62

set up as a sum of three OUP-generated series,

EC(t) = OUP1 (A1 · σ 2 , τ1 ) + OUP2 (A2 · σ 2 , τ2 ) + OUP3 (A3 · σ 2 , τ3 ) + µ.

(7)

The amplitudes Ai and correlation times τi are obtained from the analysis of time series of EC in MD simulations, see above. The mean value µ and variance σ 2 are those found for the various base-pair steps in MD simulations.

Figure 10: Left: Time series of EC taken from an MD simulation (black) and created by the implemented algorithm (blue; running average over 20 time steps: red). Right: Distributions of EC yielded by the theoretical model (black) compared to data from QM/MM simulation (orange). Fig. 10 shows an example of the generated time series (blue line). A certain deficiency is apparent in the time series generated with the combination of OUPs, namely the fact that the continuity of the time series is not enforced. Such a behavior may be problematic because continuously varying Hamiltonian is desirable when non-adiabatic propagation schemes are to be used (in particular, a surface hopping framework). To overcome this issue, a pragmatic approach has been employed consisting in the calculation of running averages over the last 20 time steps, which makes the generated time series smoother. Finally, the time series of EC obtained with the theoretical model resemble closely those obtained in MD simulations, as illustrated in Fig. 10 for an A|A base-pair step in a dsDNA oligonucleotide. Also, while the specific peaks present in the power spectrum of EC cannot 20 ACS Paragon Plus Environment

Page 21 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


be reproduced, the ACF of the modelled time series of EC has the correct tri-exponential form, see Fig. 9.

2.3

Electronic relaxation from second-order corrections: a Hubbard model

The site energies are interpreted as IP or EA on the basis of Koopmans’ theorem in the HF theory, and Koopmans’ theorem neglects relaxation effects, i.e. the response of electronic structure to the presence of net charge (electron or hole). In DFT, the single-particle energies can be derived from the total energy Etot as the first derivatives with respect to occupation numbers, or the charge on site n, 61 εn =

∂Etot . ∂qn

(8)

The HOMO eigenvalue is related to the chemical potential of the fragment. Then, the electronic relaxation effects are described by second-order contributions, leading to the chemical hardness. This quantity describes the response of the chemical potential to the occupation number, and is known in physics as the Hubbard parameter. In general, one can write

Γmn =

∂ 2 Etot , ∂qm ∂qn

(9)

where qm is the net charge on molecular fragment m. The diagonal term Um = Γmm corresponds to the double of the chemical hardness (Hubbard parameter) of the fragment, and the off-diagonal term Γmn represents basically the Coulomb repulsion between the net charges on the fragments m and n, and is inversely proportional to their distance, Γmn ∝ 1/Rmn . The energy of HOMO is lowered (raised) by the amount Um · ∆qm with ∆qm = ±1 when the hole (excess electron) is localized on the site m completely. Further, the relaxation may lead to a moderate rotation of lower-lying orbitals, while the higher occupied orbitals are rather unchanged in shape, as was shown in our earlier work. 32



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 62

The total electronic energy is written as follows: 32–34

Etot = E 0 +

X

|am |2 εm +

m

X

a∗m an Tmn + C ·

m6=n

X

Γmn ∆qm ∆qn

(10)

m,n

Here, E 0 is the zeroth-order term, or the energy of the system in a charge-neutral state. The second-order term would only contain orbital relaxation effects in the HF theory. In DFT, however, this term is contaminated by the notorious self-interaction error. Therefore, it is in order to apply a self-interaction correction, and this is accomplished by scaling down with a constant factor C as detailed in Ref. 32. No additional parametrization is needed for the diagonal terms because these are constant. As for the off-diagonal terms, the center-of-mass distances of the respective purine nucleobases were obtained from MD simulations of dsDNA oligonucleotides. A constant value of distance is considered in theoretical simulations, corresponding to the mean value observed in MD simulation; see the Supporting Information for the numerical values.

2.4

Nuclear relaxation effects: Inner-sphere and solvent contributions

Nuclear relaxation effects are accounted for within Marcus’ theory, being described by the reorganization energy λ. They can be divided into the so-called inner-sphere reorganization, quantified by the parameter λi , and the outer-sphere contribution with a corresponding λo .

2.4.1

Inner-sphere reorganization

The change of the internal geometry of nucleobases participating in the ET reaction is challenging to describe with classical MD simulation, for two reasons: principally, classical (Newtonian) mechanics may not be a good approximation to describe accurately the dynamics of high-frequency intramolecular vibrations, which may possess a fair amount of quantum character. Secondly, it is not clear if sufficient accuracy can be expected from 22 ACS Paragon Plus Environment

Page 23 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a parametrization of bonding interactions depending on the charge state of the molecule. Therefore, even our previously developed direct atomistic simulation of ET cannot actually involve any dynamic effects that would translate into λi when ensemble averaged. The simulation protocol contains a phenomenological description of the inner-sphere reorganization, 62 which is identical in form to that in Marcus’ theory. The energy difference due to the change of geometry is obtained as

0 0 (neutral geom.) (charged geom.) − Em λim = Em

(11)

where the energies are computed for a charge-neutral molecule in different geometries. Thus, only the energy difference due to change of geometry is considered, 32 and the following contribution can be added to the total energy:f i

Eλ =

X

λim (∆qm )2

(12)

m

Within this approximation, the inner-sphere relaxation occurs instantaneously, at the same time that the change of electronic structure occurs. This is not a necessary choice, and a suitable relaxation function might be applied instead, as detailed for the outer-sphere contribution below. However, in the case of ET in solvated DNA, this term is small and the relaxation is fast compared to the ET event, therefore the assumption of instantaneous relaxation has probably only a minor effect.

2.4.2

Solvent reorganization

An excess charge (hole, electron) located on one of the nucleobases induces a strong polarization of the molecular environment – the solvent as well as the remaining parts of the DNA (backbones, other nucleobases) – and affects several close nucleobases considerably. The enf In the implementation, this term is added to the second-order contribution in Eq. 10 because it has the same functional form.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 62

vironment is polarized by translating and reorienting polar solvent molecules, ions and parts of DNA, and in return, it induces an electric field on the nucleobases. The site energies εn of nucleobases in the vicinity are decreased by as much as 3 eV (or 2 eV when corrected for the missing electronic polarizability, see further below). Consequently, the excess charge is stabilized on the site where it is located, thus the polarization represents a charge-localizing force. The decrease of site energies along the DNA strand for various nucleotide sequences is presented in the Supporting Information. The decrease of IP in response to the presence of excess charge (∆qm on site m) is modelled as follows: ∆εn = c ·

sites X

∆qm · fdrop (|n − m|)

(13)

m

Here, the function fdrop captures the site energy decrease on neighboring nucleobases. The function is based on the decrease of site energies in the polyA oligonucleotide, and its shape is shown in Fig. 11; the numerical values of fdrop may be referred to in the Supporting Information. fdrop has the largest negative value for the site where the excess charge resides (m = n), and it converges towards zero for distant sites (large |m−n|). c is a factor correcting for the missing electronic polarization of the environment described with a classical force field, as discussed in Sec. 2.4.3; it shall be set to 1/1.4.

Figure 11: The effect of environmental polarization due to the presence of excess charge on the DNA: The function fdrop , which models the decrease of IP of nucleobases on an oligonucleotide with a hole fully localized on nucleobase marked (+).


Page 25 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As soon as the charge is transferred to another site (nucleobase), the environment repolarizes, and the newly occupied site experiences the lowest site energy and thus the strongest stabilization of the excess charge. This is the molecular nature of the outer-sphere reorganization, which poses a barrier to the ET reaction and is quantified within Marcus’ theory by the outer-sphere reorganization energy λo . The reorganization of solvent does not happen instantaneously, and a rather slow decay is observed following an ET event instead. To monitor the relaxation, 1000 non-equilibrium classical MD simulations of 20 ps each were performed, starting from structures accumulated along an equilibrium simulation of DNA with no excess charge. At the start of each nonequilibrium simulation, the MM atomic charges were modified manually to simulate a ‘forced’ elementary ET reaction. In each of these 1000 simulations, a hole was introduced onto a nucleobase, and the time series of site energy was recorded. A mean value of site energy was obtained for every time step by averaging over the 1000 simulations; the resulting relaxation constitutes the basis for the function frelax (see below) and is shown in Fig. 12. The averaging over such a large number of simulations is meant to cancel any stochastic fluctuations, and the retained fluctuations of IP are considered a net effect of the ET event.

Figure 12: Decay of the site energy of an adenine in the polyA dsDNA oligonucleotide following the transfer of a fully localized hole from a neighboring nucleobase to the one whose site energy is monitored. Comparison of the relaxation function implemented in the theoretical model (black) and the data obtained by averaging over 1000 non-equilibrium MD simulations of polyA DNA of 20 ps each (red).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

We may put the observed relaxation behavior in relation to some recent simulation studies. So, Furse and Corcelli studied the relaxation following the de-excitation of a fluorescent dye molecule bound to DNA; 63 this was motivated by time-dependent Stokes shift experiments. Reported was multi-exponential dependence with decay times of 1.4 and 20 ps, possibly not mentioning a fastest decay time. A tri-exponential fit to the relaxation of polarization in our work yields decay times of 8.0 fs, 0.22 ps and 2.4 ps (the respective relative magnitudes are 53 %, 34 % and 14 %). Thus, the decay is faster than that reported previously. There may be various reasons for the differences between our current results and the mentioned previous work: the different reaction (de-excitation vs. electron transfer), different composition of the system (ligand bound to DNA vs. pristine DNA) and possibly different force fields, above all regarding the water model. Further research along these lines is underway in our laboratory currently. To model the time dependence of this relaxation, one of different analytical functional forms could be considered. A straightforward choice would be a single exponential or a sum of multiple exponentials, possibly modulated by cosine functions to represent the most significant vibrational modes being responsible for the relaxation process.g In the current work, a tabulated, non-analytical relaxation function frelax is used in order to increase the flexibility of the approach and to avoid the search for a suitable functional form. The final form of the relaxation function is obtained by normalizing the range of the function to the interval (0,1); the length of the function considered currently is 12,639 fs. To reduce the numerical noise due to a finite number of simulations used to generate the time series, the relaxation function was smoothed with a running-average procedure, starting at the mark of 0.6 ps; the length of the running average increased gradually from 3 fs at the beginning (0.6 ps) to 2.4 ps at the end (12.6 ps). It was made sure that the relaxation function reaches zero at the end. With all of this information available, the contribution of time-dependent environmental g

Further, to prevent a sudden change of site energy at the onset of the exponential function, a Gaussianshaped function may be added.


Page 26 of 62

Page 27 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


polarization to the site energies is modeled with a function constructed as

∆εn (t) = c ·

sites X 0 X

(∆qm (t) − ∆qm (t − ∆t)) · frelax (t) · fdrop (|n − m|),

(14)

m t=−T

where ∆qm (t) is the ‘charge history’ of site m, or the portion of excess charge located on site m at time t before the present time step 0; ∆t is the ‘nuclear’ time step, and T corresponds to the length of the used relaxation function; refer to the Supporting Information for the numerical values of frelax . We note that the evaluation of this convolution is the most timeconsuming part of the simulation with the theoretical model. The computational complexity is linear with the length of the relaxation function implemented, and it is even quadratic with the number of fragments participating in ET. With the employed description of environmental polarization and of its dynamics, the magnitude of polarization is assumed to be directly proportional to the magnitude of excess charge present on the respective molecular fragment (nucleobase). Additionally, this effect is assumed to be additive, to describe the polarization due to a hole or an excess electron that is delocalized over several fragments. Data in support of this approach may be referred to in the Supporting Information.

2.4.3

Effective treatment of electronic polarization

The MM environment has been represented by fixed atomic charges in all of the QM/MM calculations throughout this work, and thus it does not possess any electronic polarizability. This approximation represents a likely source of errors whenever the electronic polarization is changing, and ET reactions are such a case. It has been established that the reorganization energy is overestimated in such a situation. 64–68 Of the results presented in the current work, the affected quantities are the decrease of site energies in the vicinity of the excess charges (which translates into the reorganization energy directly), and further the fluctuations of site energies due to the environment. It is also known that such energy-related quantities



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 62

are overestimated by a common factor equal to the square root of the Pekar factor of the medium, which coincides with its optical dielectric constant, εel , for a medium with a high static dielectric constant. 66,69 Thus, an easy solution to this problem is to scale down the affected quantities by a corresponding constant factor. 70,71 Truly, correct solvation energetics were observed with a uniform scaling of atomic charges of ions. 72–74 It is thus justified to expect the energetics of ET to improve upon such a uniform scaling. In this work, the drop function and the fluctuation of site energies are scaled down with the √ same factor, which shall correspond to c = 1/ εel (e.g. in Eq. 14), and a value of c = 1/1.4 is considered.

2.5

Total energy expression

In our previous work, we have derived a general total energy expression from DFT, which is suitable for an efficient description of ET problems. 32–34 The total DFT energy functional for the charged system containing a hole or an excess electron with total charge density ρ was expanded around the charge neutral system with charge density ρ0 up to the second order, E[ρ] = E 0 [ρ0 ] + E 1 [ρo , ∆ρ] + E 2 [ρ0 , (∆ρ)2 ],

(15)

with ∆ρ = ρ − ρ0 . The zeroth-order contribution depends on the neutral charge density only, i.e. it corresponds to the energy of the neutral DNA oligonucleotide in solvent, and is set to zero in the theoretical model in the current work. The first-order contribution E 1 involves the Hamilton matrix Hmn , and the second-order terms E 2 represent the electronic relaxation effects described in Sec. 2.3. Therefore we can write

Etot =

X m

|am |2 εm +

X

a∗m an Tmn + C ·

X

(Umn + δmn λim )∆qm ∆qn ,

(16)

mn

m6=n

with the two components of site energies εm = ε∗m +∆εm modeled as given by Eqns. 4 and 14, respectively. Taking derivatives with respect to am , the electronic Hamilton matrix elements 28 ACS Paragon Plus Environment

Page 29 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are obtained: Hmn =

2.6

  P  εm + k (C · Ukm + δkm λim )∆qk

for m = n,

  Tmn

otherwise.

(17)

Quantum chemical propagation

Such a total energy expression constitutes a convenient basis of a semi-classical propagation scheme that describes the dynamics of the system in a combined way, for both the atoms (nuclei, ions) and electrons. In our previous work employing all-atom molecular dynamics, 33,34 the surface hopping scheme (SH) 75 was implemented in its fewest-switches variant. 76 While an application of accurate quantum chemical methods yielding the necessary non-adiabatic couplings was possible in an earlier work, 30 our approach has involved the development targeted to semi-empirical methods. 77 The SH implementation described in Ref. 34 is used in the current work as well. Obviously, there is a significant difference between the all-atom simulations and the new theoretical model of electron transfer: There are no explicit atom coordinates in the theoretical model, therefore no Newton equations of motion for the atoms are being solved. The dynamics of the atomic system is represented by the temporal variance of the site energies and EC, which is an integral component of the model and does not require any additional treatment. What is to be treated in the same way as in the all-atom simulation, however, is the dynamics of the electronic system, or more precisely speaking, of the excess charge – the transferring hole. The semi-classical character of the simulation, requiring a coupling of the electronic system to the ‘classical’ system, is accomplished with an additional contribution to site energies, ∆εm in Eqn. 14, see Sec. 2.4.2. This captures the response of site energies to the polarization of the aqueous solvent, which itself has responded to the distribution of excess charge on the DNA molecule. A problem that needs to be resolved in an application of SH is the lack of decoherence, which was shown to make the fewest-switches algorithm fail in regions of weak non-adiabatic



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

coupling. 78 Among the solutions proposed, the ‘smooth decay of coherence’ scheme introduced by Truhlar et al. 79 is constructed in such a way that it can be applied within our theoretical model approach. Recall that within SH, the time-dependent Schrödinger equation is integrated to yield the expansion coefficients of the individual adiabatic states, bi , and after every step of integration, the state/surface is identified on which the consequent step of classical propagation shall be performed, i.e., the so-called populated state. In our theoretical model, this means that the charges of the individual sites, ∆qm , are calculated from the eigenvector of the populated state. The idea of decay of coherence is to let the expansion coefficients bi of all of the adiabatic states but the one currently populated, decay exponentially. To this end, a decay time is determined first for every non-populated state from the energies of the adiabatic states, τi = c · h ¯ / |Ei − Epop |. The multiplicative factor c is modulated by the kinetic energies in the original scheme, while it is considered constant, c = 3, in our application lacking explicit kinetic energies. Then, the expansion coefficients bi of the all of the adiabatic states, but the one currently populated, are scaled down, b0i = bi ·exp [−∆t/τi ]. Finally, the expansion coefficient of the populated state is scaled up by a corresponding factor, to maintain a normalized wave function. The current implementation of SH involves an additional minor modification, which regards the issue of the so-called frustrated hops in SH. 76 These are the occasions where a transition to a higher-energy state shall be performed but there is not enough kinetic energy in the vibrational mode of motion parallel to the non-adiabatic coupling vector. Within the theoretical model of electron transfer, there are no explicit atom coordinates and velocities, as well as no explicit non-adiabatic couplings, and thus it is impossible to even analyze whether a frustrated-hop-like situation occurs. To resolve this issue in an approximative way, a simple energy-based criterion is introduced: A predicted transition to a higher-energy state will be performed as long as the difference of energies of the involved adiabatic states is smaller than a certain threshold. This threshold shall correspond to some average amount of energy available in the system, and we employ a value of 0.2 eV.


Page 30 of 62

Page 31 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Despite looking rather high, corresponding to ca. 8 kB T at 300 K, this value is physically motivated. For one thing, the strongest skeletal vibrations of nucleobases, which are the modes of motion probably playing a role in the relevant non-adiabatic couplings, possess typical frequencies of ν˜ = 1800 cm−1 . The energy e = hc˜ ν of such a mode of motion is ca. 0.2 eV. Another point is related to the magnitude of the fluctuations of the site energies, which are the quantities that control the energetics of ET. The previously reported standard deviation of site energies of nucleobases of 0.3 eV divided by the scaling factor of 1.4 as described in Sec. 2.4.3 gives ca. 0.2 eV as well. Finally, let us discuss the issue of the thermal Boltzmann population, which should be reached at long times, and of the detailed balance between the surface hops up and down in energy. In the usual framework of semi-classical simulation, the adjustment of atomic velocities and the appropriate treatment of frustrated hops are essential to achieve detailed balance. 80 Since our model does not involve any atomic coordinates and thus velocities, such procedures cannot be applied. It has been discussed recently that direct fewest-switches SH dynamics provides detailed balance approximately, 81 while development on this issue is underway for a rather broad family of SH schemes. 80 Thus, it is well plausible that the on-going improvements of SH schemes will make it possible to enforce detailed balance even in the dynamics of such an abstract model like that one currently introduced. Moreover, we note that the other principal approach to non-adiabatic dynamics, the (Ehrenfest) mean-field method, does not provide detailed balance. 82

2.7

The simulation workflow

The algorithm of ET simulations within the theoretical model is derived from the previously developed atomistic multi-scale non-adiabatic scheme. Therefore, the sequence of calculations of the quantities occurring in the model during one step of the dynamics simulation mirrors the sequence of calculations in the atomistic simulations. How the model represents the ingredients that determine the ET process has been explained above in detail, and the 31 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

underlying idea of the theoretical simulation is illustrated in Fig. 13.

Figure 13: Descent of the simulation workflow of the theoretical model (right) from the original atomistic multi-scale scheme (left). In the following, the calculation scheme is described in detail. It is still a variant of a molecular dynamics simulation realized by a numerical integration of the time-dependent Schrödinger equation. Thus, a time step is defined, and one iteration of the loop in Fig. 13 is performed per step. The following computations are involved in one iteration step: 1. The coarse-grained Hamilton matrix is calculated. The diagonal elements, i.e., the site energies, are obtained using the cosine series in Eq. 4, and the contribution from the polarization of the environment from the previous iteration step is added. The off-diagonal elements, i.e., the electronic couplings, are modeled using the Ornstein– Uhlenbeck processes in Eq. 7, and running averages of these time series are employed. The Hamilton matrix is set up as given in Eq. 17. 2. With the Hamiltonian at hand, a step of the quantum mechanical propagation of the excess electron/hole is performed. Different approaches may be applied here. These schemes include, for instance, a simple diagonalization of the Hamiltonian, employing the Born–Oppenheimer approximation, or, in the technically simplest non-adiabatic approach, integration of the time-dependent Schrödinger equation for one step. Following our previous work, a SH scheme is applied and complemented by a simple decoherence correction. Whatever propagation scheme is used, the result is always a linear combination of basis functions represented, actually, by the individual charge carrying sites, which are the nucleobases in the case of DNA ET. Importantly, the 32 ACS Paragon Plus Environment

Page 32 of 62

Page 33 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


expansion coefficients of the wave function aim provide directly the occupation of the site m by the excess electron/hole, or the ‘charge’ of the site ∆qm = |aim |2 . Within SH, aim are the expansion coefficients of the currently populated adiabatic state. 3. The new charges on the sites are used to compute the contributions to site energies induced by the polarizable molecular environment, ∆εm . This calculation is based on the above considerations of how the polarization of environment relaxes upon an elementary ET event, and is performed by evaluating the convolution in Eq. 14. This scheme can be thought of as substituting the integration of Newton equations of motion in an ordinary MD simulation – a substitute that concentrates on the electric field induced by the molecular environment. Note that this calculation represents the largest part of the computational cost; it scales with the square of the number of sites in principle, and depends linearly on the length of the relaxation function frelax .

3 3.1

Results and discussion Homogeneous DNA sequences – AAAA and GGGG

Hole transfer along a homogeneous DNA sequence was examined for the oligonucleotides AAAA and GGGG in 50 model simulations of 10 ns each, for each of the DNA species. The purpose of these simulations is, firstly, to compare the performance of the new theoretical model with the results from previous atomistic multi-scale QM/MM simulations. 33 To make such a comparison possible, ET events were counted in the theoretical simulations in a way that is closely related to the counting employed in Ref. 33: The lowest site energy was identified and the hole was considered as located on this site provided one of the following conditions was met: (i) the hole was located on this site in an immediately preceding period of the simulation already; or (ii) the hole was located on another site though, but the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

second-lowest site energy is at least 0.64 eV higher than the currently lowest.h Note that it is impossible to count ET events in exactly the same way as in Ref. 33 because the values of electric potential (ESP) used in Ref. 33 are not explicitly available within the theoretical model. Table 3: Delocalization of the hole (Srec , mean ± std. dev.) and the rate of hole transfer (in ns−1 ) in the DNA species AAAA and GGGG obtained with the newly developed theoretical simulation as compared to the previous results from full non-adiabatic QM/MM simulations. 33

QM/MM – Ref. 33 ET model – this work

AAAA Srec rate / ns−1 1.07 ± 0.15 104 1.05 ± 0.16 86

GGGG Srec rate / ns−1 1.01 ± 0.06 77 1.02 ± 0.08 67

The results are presented in Tab. 3 in terms of the rate of ET events per nanosecond and the delocalization of charge (expressed as the inverse sum of charges of sites squared, P Srec = 1/ i Q2i ). The excess hole is localized nearly perfectly on a single nucleobase, on average, while the extremely tiny delocalization is slightly larger in AAAA than in GGGG. The rate of transfer is on the order of ten transfers per nanosecond, and is slightly larger in AAAA than in GGGG. The rates of transfer are only slightly smaller with the theoretical model as compared to the atomistic QM/MM simulations, and the overall agreement of the results is very good. The energy criterion on SH events introduced in the current work (see Sec. 2.6) represents an improvement over the non-adiabatic schemes employed in the previous work. It eliminates some of the possibly occurring highly unlikely state transitions that may have been present in the previous work; these would lead to artificially increased rate of electron transfer conceivably. The results obtained from simulations of the oligonucleotide AAAA that employed the energy criterion are presented in Tab. 4. Very slight delocalization of charge is found in all of the simulations, and the rate of elementary transfer is found in the h

The value of 0.64 eV was motivated by the following consideration: The difference of ESP between the site occupied by a full hole and the neighboring site amounts to ca. 1.2 eV. With the scaling of ESP by the factor of 1/1.4 to correct for the missing electronic polarization, 0.86 eV results. For the criterion to consider ET events, we choose 3/4 of this value arbitrarily, which gives 0.64 eV.


Page 34 of 62

Page 35 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interval of 10–20 ns−1 , which is clearly smaller than obtained without any SH criteria, above. Importantly, these values are in good accordance with the experimental elementary rates of transfer, namely 20 ns−1 by Majima et al., 83 and 1.2 ns−1 by Lewis et al. 18 Table 4: Delocalization of the hole (Srec ) and the rate of hole transfer (in ns−1 ) in the DNA species AAAA obtained with two different energy criteria for surface hopping; mean ± std. dev. presented. In “0.2 eV / 0.2 eV”, an energy threshold was applied for state transitions to higher and lower-lying states; in “0.2 eV / ∞”, a transition to a lower-lying state was always allowed. SH hop criteria 0.2 eV / 0.2 eV 0.2 eV / ∞

AAAA Srec rate / ns−1 1.06 ± 0.13 16.2 ± 7.4 1.06 ± 0.12 12.6 ± 6.8

The favorable comparison with both our previous work, from which the current scheme has descended, and with previous experimental work indicates that the construction of the model is meaningful. The previously developed all-atom simulation scheme is able to describe the microscopic mechanism of ET. The current theoretical model has been developed and parametrized on the basis of the atomistic scheme, and it seems that the quality of the results is not compromised by the abstraction process.

3.2

Transfer over an adenine bridge – GTn GGG

Further, we apply the scheme for the hole transfer over an adenine bridge of varied length. Representing an archetype of distance-dependent electron transfer in polarizable medium, there has been an enormously rich discussion about the mechanism of transfer in such DNA species. 7 The performance of the new theoretical model may be benchmarked by means of simulations of hole transfer in the simple A-bridge containing sequence GAG, which was previously simulated with the atomistic non-adiabatic QM/MM scheme. 33 Here, 100 model simulations of 1 ns each were performed. The rate of hole transfer was obtained by an analysis of the occupation of the sites, i.e. the charges ∆qm = e · |am |2 . A completed transfer event was 35 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 62

counted whenever the initial G was occupied for at least 20 ps and then the terminal G was occupied for at least 20 ps, or vice versa. Thus, no ET event is counted when there is an immediate backward transfer within 20 ps. Table 5: Delocalization of the hole (Srec ), occupation of the bridge A, and the rate of hole transfer (in ns−1 ) in the DNA species GAG obtained with the newly developed theoretical simulation as compared to the previous results from full non-adiabatic QM/MM simulations; 33 mean ± std. dev. presented.

QM/MM – Ref. 33 ET model – this work

Srec 1.00 ± 0.02 1.00 ± 0.01

GAG Q(A) / e 0.005 ± 0.065 0.001 ± 0.016

rate / ns−1 4.6 ± 2.0 3.0 ± 1.3

The results are presented in Tab. 5 in terms of the averaged delocalization of charge, Srec , the averaged occupation of the A-bridge, Q(A), and the rate of ET events per nanosecond. The occupation of the bridge adenine of 0.001 on average is somewhat smaller with the theoretical model than in the atomistic simulations, and this may be why the overall rate of transfer is 35 % smaller. Importantly, the applications of this theoretical model will aim at explaining the mechanism of ET and finding the differences in ET rate between different DNA sequences. In spite of the not quite perfect agreement with the atomistic non-adiabatic QM/MM simulations, we thus consider the performance of the theoretical model acceptable. Similarly to the previous work of ourselves as well as others, we turn our attention to dsDNA of the sequences GTn GGG now, in which the hole is placed on the left-terminal G initially. Sets of simulations of hole transfer extending to 100 ns were performed for each of the DNA species: the number of simulations in the set was 256, 256, 256, 384, and 512 for n = 1, 2, 3, 4, 5, respectively, providing a total sampling of 25.6–51.2 µs. The larger number of simulations for the longer DNA species were needed to improve the statistics, because of the less frequent hole transfer in these molecules. The energy criteria for SH of 0.2 eV/∞ were applied. The site energies and the occupations of sites were recorded in every 100th step only, to reduce the required amount of disk space and time needed for analysis. The rate of hole transfer was obtained by an analysis of the occupation of the sites, i.e. the 36 ACS Paragon Plus Environment

Page 37 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


charges ∆qm = e · |am |2 . A completed transfer event was counted whenever the initial G was occupied for at least 20 ps and then the terminal GGG was occupied for at least 20 ps, or vice versa. Thus, no ET event was counted when there was an immediate backward transfer within 20 ps.

Figure 14: Hole transfer rates in simulations performed with the theoretical model on the DNA species GTn GGG. Data are based on 25.6–51.2 µs of simulation per DNA species. The results are presented in Fig. 14. There is a transfer from G to GGG and back with a rate that is on the order of units to tens per microsecond. In other words, the hole transfer occurs once per 20–1000 ns on average, over an adenine bridge constituting of one to five nucleobases. The dependence of rate on the number of adenines (or, equivalently, on distance) appears exponential over the entire range n = 1, . . . , 5; the exponent β in the distancedependence exp[−βr] takes the value of 0.8 per nucleobase, or 0.2–0.3 ˚ A−1 . This may look surprising when compared to the crossover from exponential to algebraic dependence at ca. n = 3, which has been reported previously and discussed richly; 11,14 there is no such crossover in our simulations. Still, the transfer is ca. 1000× slower than the transfer within a polyadenine sequence in Sec. 3.1, thus the ‘excitation’ from a guanine onto the adenine bridge is the bottle neck within the theoretical model, correctly. Although we are unable to resolve this issue definitely, possible causes of discrepancy are obvious: Not quite equivalent quantities are compared, namely the rate of transfer and the relative efficiency of terminal oxidation, and possibly several different processes are also included in the experimentally reported rates, like 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 62

the charge injections or a terminal charge recombination. Also, importantly, the tunneling model of ET over A-bridges was questioned previously, and the non-zero occupation of the bridge, which was also found previously by others, may have consequences for the distance dependence. The occupation of the bridge during a transfer event is a feature that will require dedicated examination with possibly improved non-adiabatic propagation schemes. To study the microscopic mechanism of transfer in the simulations in more detail, the individual elementary transfer events shall be analyzed now. To this end, additional sets of 400 simulations of 10 ns each (for a total of 4 µs) were performed for every DNA species, and the hole location was recorded in every step of the simulation (i.e., every 1 fs). The limited length of these simulations is due to the very large amount of disk space necessary to save the data and the time needed for analysis. The following elementary transfer events were distinguished and counted: from the initial G onto An and back to G (G→An →G, no net transfer); from the terminal GGG onto An and back to GGG (GGG→An →GGG, no net transfer); from G onto An and further to GGG (G→An →GGG); from GGG on An and further to G (GGG→An →G); from GGG directly to G without An being populated transiently (GGG→G); from G directly to GGG (G→GGG). Only those transfer events are counted where the newly formed state (G+ or GGG+ ) remains stable for 0.1 ps. Table 6: Rate (in µs−1 ) of elementary hole transfer events in surface hopping simulations of hole transfer in DNA sequences GTn GGG. Data are based on 4 µs of simulation per DNA species. n 1 2 3 4 5

G→An →G GGG→An →GGG inactive excitation onto bridge 15 44 105 40 103 24 97 18 112 6

G→An →GGG GGG→An →G transfer – with bridge occupied 18 12 10 9 4 4 2 2 1 1

GGG→G G→GGG – bridge not occupied 7 1 1 0 1 1 0 0 0 0

sum of all 97 165 137 120 120

The results expressed as number of events per microsecond are presented in Tab. 6. The total rate of elementary events is roughly similar for all of the cases n = 1, . . . , 5, in the interval 100–200 µs−1 , and we may rationalize this by considering the rate of excitation onto


Page 39 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the bridge being independent of the bridge length. Notably, a hole that has transferred onto the adenine bridge may return to the guanine where it resided initially, and such an event does not produce any net electron transfer; it will be called ET-inactive in the following. It is of interest to analyze and compare the rate of ET-active and ET-inactive events, depending also on the number of adenines. So, the rate of ET-inactive events G→An →G is rather independent of the adenine bridge length, and it is quite high, so that it represents the majority of events (except for n = 1). Here, the hole hops onto the first, neighboring adenine and then returns to the initially occupied guanine. No neighboring nucleobases play any role, thus it does not matter how many adenines there are. By contrast, the rate of ET-inactive events taking place at the terminal guanine triplet, GGG→An →GGG, decreases rather sharply with the increasing number of adenines in the bridge. This may be explained by the rather limited length of the individual simulations, which was 10 ns each. In order to observe an event of this kind, the hole has to transfer from the initial guanine to the guanine triplet first. This takes a certain amount of time, and only then an event GGG→An →GGG may be observed; in an extreme case, the hole will never make it to GGG within the simulation, and GGG→An →GGG is impossible altogether. The events G→An →GGG and GGG→An →G constitute the larger part of ET-active elementary events. Here, the hole occupies the adenines in the bridge transiently, as opposed to the events in the simulations where this is not observed, G→GGG and GGG→G. Correlating with the overall rate of transfer, the rate of these events seems to decay exponentially with the number of adenines on the bridge. In our opinion, these events fulfill the definition of TIH mechanism. Still, the dependence of rate on the length of the adenine bridge is observed to be exponential. While this is at odds with the initial proposals that identified algebraic distance dependences with the TIH mechanism, it was shown later that exponential dependence may be due to other reasons. 22,84,85 Although we cannot provide a clear explanation of this situation, our results indicate that the relation between the mechanism and the distance



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dependence of the rate is perhaps more complex. Also, a certain number of ET-active events were observed in the GTGGG species (n = 1) during which the bridge was not occupied at all, designated here as G→GGG or GGG→G. Although one may be tempted to consider this a signature of transfer via a tunneling mechanism through the bridge, as opposed to TIH, such a conclusion would be far-fetched, in our opinion. Instead, we shall like to point out a recent finding by Prezhdo et al. 86 of possibly underestimated rate of tunneling with a fewest-switches SH propagator. An alternative ‘global-flux SH’ propagator was proposed; application to electron transfer simulation is underway in our laboratory. Lastly, let us investigate in more detail the energy relations in the electron transfer system in terms of site energies. The site energies determine the energetics and co-determine the kinetics of any electron transfer, therefore they are in the spotlight of any discussion of ET mechanism, including the recently proposed complex mechanisms, e.g., the flickering resonance. 22 We recall that ‘instantaneous values’ or vertical values of IP play the role of site energies as long as hole transfer is considered. We analyzed the hole transfer in the DNA species GTGGG and GT2 GGG during the ET events G+ →A(A)→GGG involving the bridge being populated by the hole. Monitored were the site energies of the initially charged G, of the bridge A or As, and of the terminal GGG (the site energy of the central (G)G(G) was considered, and we note that three values were always very similar, data not shown). The values of differences ε(A) − ε(G+ ) as well as, optionally, ε(GGG) − ε(A) or ε(A2 ) − ε(A1 ) were obtained in two ways: (red) from the intervals of the simulations during which the ET events were taking place; and (blue) from simulations with a hole placed on the initial guanine with no ET allowed in the simulation. The distributions of the obtained data series are shown in Fig. 15. It is apparent from the distributions of ε(A) − ε(G+ ) in GTGGG that ET events occur when the site energies of the initially charged G+ and the bridge A are much closer (red) than they are on average (blue in Fig. 15). This is rather understandable: the hole transfer is facilitated whenever


Page 40 of 62

Page 41 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 15: Left: Energetics of hole transfer in the dsDNA species GTGGG. Distributions of differences of IP, (top) ε(A) − ε(G+ ) and (bottom) ε(GGG) − ε(A). Data obtained from portions of simulations during which transfer events occurred (red), and from resting-state simulations with the initial guanine charged (blue). Right: ditto in the dsDNA species GT2 GGG, distributions of (top) ε(A1 ) − ε(G+ ) and (bottom) ε(A2 ) − ε(A1 ). the structural fluctuations (that the theoretical model was fitted to reproduce) induce such an electric field on the nucleobases that the energy barrier for transfer is decreased or even vanishes. However, the site energies of the bridge A and of the terminal GGG are by no means closer than they are on average (with no ET taking place). Their difference is actually even slightly higher than in the resting-state situation (with a hole stationary on the initial guanine). The site energy difference of the initial G and the terminal GGG during ET events is similarly large, as well. An observation similar to that described in this paragraph is also made for the hole transfer event GGG+ →A→G, see the Supporting Information. Thus, while the site energies of the initial state and of the bridge are in resonance during ET events, there is no such resonance for the remainder of the system. Still, let us investigate the transfer of the hole onto the A-bridge alone, in the dsDNA species GT2 GGG. What Fig. 15 (right) shows in addition is the difference of site energies within the bridge, ε(A2 ) − ε(A1 ). Evidently, the distribution of this difference during the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ET events differs little from that in the resting state, and we do not observe any resonance over the entire bridge, at the moment that the hole is hopping onto the bridge. There is no energetic basis for a delocalized state of the hole, and thus our simulations do not seem to support the mechanisms of tunneling or flickering resonance.

3.3

Sensitivity analysis

Whenever a parametrized theoretical model is being built, it is of interest to examine how sensitively the results respond to any change of parameter values in the model. In a broader sense, such a manipulation of the model provides an idea of the way that individual physical phenomena on the molecular level affect the ET process. This kind of knowledge may make it possible to compare the properties, such as the rate of ET in different molecular systems as long as they are characterized in terms of the parameters contained in the model. To this end, we have constructed several alternative models of ET in DNA by modifying one of the involved parameters at a time. Subsequently, each of the models was used to simulate hole transfer in simple homogeneous and heterogeneous DNA sequences: for AAAA, eight simulations of 10 ns each were performed for a total of 80 ns, while 80 ET trajectories of 10 ns each were generated for GAG for a total of 0.8 µs. Each of the trajectories within this ensemble was generated with a different initial seed for the random number generator to guarantee the independence of ET trajectories. The different modifications of the model are presented in Tab. 7, together with the resulting values of hole delocalization (Srec ) and ET rate. While the number of simulations performed is rather small, qualitative conclusions may still be drawn from the results obtained. The IP difference between guanine and adenine controls the height of the energy barrier opposing the hole transfer between guanines across adenine bridges, and its unperturbed value is 0.36 eV. In the simulation series 1A and 1B, this difference was increased and decreased by 0.05 eV, respectively, and the rate of transfer decreases to 14 and increases to 35 µs−1 , respectively, from the value of 15 µs−1 with the unmodified model. The roughly 42 ACS Paragon Plus Environment

Page 42 of 62

Page 43 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 7: Sensitivity of the ET rate to the change of parameter values. Srec and ET rate from simulations of hole transfer in the DNA sequences AAAA and GAG, obtained by averaging over the ensembles of eight and 80 simulations of 10 ns, respectively. #

modification of parameters

ref.

0 1A 1B 1C 1D 1E 1F 2A 2B 2C 3A 3B 4A 4B 4C 4D 4E

no modification site E (ε0 ) of adenines increased by 0.05 eV site E (ε0 ) of adenines decreased by 0.05 eV site E – larger fluctuation, std. dev. increased by 0.05 eV site E – smaller fluctuation, std. dev. decreased by 0.05 eV site E – higher correlation of neighbors, tc decreased by 2 fs site E – lower correlation of neighbors, tc increased by 2 fs el. couplings – smaller fluctuation, std. dev. halved el. couplings – larger fluctuation, std. dev. doubled el. couplings – slower decay, decay time doubled off-diagonal Hubbard Γmn decreased by 20 % off-diagonal Hubbard Γmn increased by 20 % polariz. of environment stronger – fdrop increased by 10 % polariz. of environment weaker – fdrop decreased by 10 % polariz. of environment faster – frelax shortened to 6.3 ps polariz. of environment slower – frelax extended to 25.3 ps polariz. of environment slower – frelax extended to 50.6 ps

AAAA Srec ns−1 1.066 5.2

Eq. 4 Eq. 4

Eq. 4 Eq. 4 Eq. 7 Eq. 7 Eq. 7 Sec. 2.3 Sec. 2.3 Eq. 13 Eq. 13 Fig. 12 Fig. 12 Fig. 12

1.089 1.049 1.053 1.079 1.054 1.121 1.064 1.131 1.037 1.047 1.096 1.072 1.058 1.061

17.9 0.7 1.4 12.9 3.0 9.7 5.1 7.6 6.6 1.4 15.6 20.6 1.3 1.0

GAG Srec µs−1 1.0031 15 1.0028 14 1.0036 35 1.0037 181 1.0030 0 1.0030 4 1.0034 49 1.0008 8 1.0145 138 1.0033 24 1.0045 36 1.0023 13 1.0027 1 1.0040 139 1.0031 9 1.0032 14 1.0032 21

doubled rate of transfer with decreased energy barrier meets the expectations, while the fact that the rate hardly changes upon the increase of the site energy difference is surprising. This result may be explained by insufficient sampling due to a rather small number of individual simulations performed, but physical reasons cannot be excluded. One possibility is a change of the ET mechanism upon lowering the site energy difference. This issue deserves a more detailed analysis in a future study. Another effect crucial for ET is the strong temporal fluctuations of site energies, leading to a standard deviation of 0.31 eV. Upon correction for the missing electronic polarizability of the environment, it reduces to a still substantial value of 0.23 eV. The fluctuations of site energies were manipulated in simulation series 1C and 1D, where the standard deviation of site energies was increased to 0.36 and decreased to 0.26 eV, respectively (the applied scaling reduces these values to 0.27 and 0.19 eV, respectively). The effect on ET in both AAAA 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and GAG was huge: the increased fluctuations led to an acceleration of transfer by a factor of 3 and 12 in AAAA and GAG, respectively. At the same time, the rate of transfer was nearly vanishing as soon as the fluctuation of site energies was decreased. The outstanding significance of fluctuations for DNA ET is thus confirmed, and the large sensitivity of the results to the variance of the site energies generated by the model is illustrated. Besides the adenine–guanine IP difference and the magnitude of site energy fluctuations, the correlations of site energies of neighboring nucleobases were postulated to play a role in some ET mechanisms. In simulation series 1E and 1F, these correlations were made stronger and weaker, respectively, by appropriate choices of the parameter tc . It is important to realize how this correlation affects the site energy difference between neighboring nucleobases: while the difference remains the same on average, the magnitude of fluctuations decreases as the site energies become more strongly correlated, and vice versa. This is illustrated in Fig. 16 where distributions of site energy differences are shown, taken from model simulations of the DNA sequences AAAA and GAG with the parameter tc modified. Evidently, a smaller value of tc , which leads to a stronger correlation of site energies, makes the distribution of site energy differences narrower. In other words, values of site energy difference that are far from the mean are reached less often; in the opposite case of a larger value of tc , they are reached more often. We may contemplate that if a hole is present in the system and thus there is a large instantaneous difference of site energies due to the reorganization energy, the resonance state with equal site energies will be reached more often for weaker site energy correlations. The reason would be the broader distribution of the site energy difference. Notably, the effect is thus the same as that of simply broadened distributions in the simulation series 1C. In the case of model ET simulations in the series 1E and 1F, this is truly the case: The weaker site energy correlations due to an increased value of tc led to faster transfer, while stronger correlations slowed down the ET. Another determinant of ET is the electronic coupling of electron carriers. In DNA, the EC between nucleobases exhibits extremely strong fluctuations due to structural variances


Page 44 of 62

Page 45 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 16: Distributions of site energy differences in model simulations performed in absence of excess charges for DNA sequences AAAA and GAG. The different choices of the parameter tc that were employed are indicated. These values correspond to the standard set of parameters ‘0’ as well as the modified sets ‘1E’ and ‘1F’. Note that the widths of the distributions change while the mean values do not. of the DNA double strand. To study the role of these fluctuations, the standard deviation of EC was halved and doubled in simulation series 2A and 2B, respectively. Smaller variance of the couplings led to slower transfer, and larger variance accelerated the transfer largely. The difference was particularly pronounced in the GAG system, for which the rate of transfer increased by an order of magnitude. It seems that the fluctuation of EC may be an important factor for hole transfer over an adenine bridge. It must be noted, however, that doubling the width of the distribution of EC is a rather harsh and unrealistic intervention. Thus, the significance of these observations is not to be overestimated. A related point is a possible role of the auto-correlation behavior of the EC. The simulation series 2C was performed using a time series of the EC constructed from an OUP with a slower temporal decay. The modified temporal correlation behavior of EC hardly affected the rate of ET, however. The current model of ET in DNA involves constant off-diagonal elements of the ‘Hubbard’ matrix Γmn (|m − n| = 1). Note that these elements correspond to the inverse of the centerof-mass distance of neighboring nucleobases. Therefore, it is approximated implicitly that the distances of nucleobases do not change. Modifications to these off-diagonal Hubbard



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

elements reported below are also considered as a simple test of how sensitive the results are to this approximation. To this end, the simulation series 3A and 3B were performed with the off-diagonal Hubbard elements decreased and increased by 20 %, respectively, which is a rather large change taking into account that the distances of the nucleobases are obtained with high accuracy from the MD simulations. In the DNA species AAAA, this change affected the delocalization of hole markedly, while the rate of transfer remained untouched. In GAG, however, the decreased Hubbard value led to hole transfer occurring twice as often. While this is a large change of course, we have to note that the change of off-diagonal Hubbard was actually huge: the 20% decrease of Hubbard corresponds to the distance of nucleobase increased from 3.8 to 4.5 ˚ A for the G|A base-pair step, and from 4.1 to 5.0 ˚ A for A|G. Such large changes of the double-stranded DNA structures can hardly occur under normal conditions. Therefore, ET rates do not appear to be overly sensitive to the value of the off-diagonal Hubbard elements. Finally, the polarization of the molecular environment and its dynamics control the ET as well. The polarization of the environment such as, e.g., an aqueous solvent, is a major localizing force on the excess charge, opposing the transfer of electron or hole, and as such it is quantified by the reorganization energy within Marcus’ theory. The present theoretical model describes this phenomenon using the function fdrop (Eq. 13). In the simulation series 4A and 4B, the values of fdrop are increased and decreased by 10 %, respectively. This makes the environment more and less strongly polarizable, respectively, and the energy barrier opposing ET is modulated directly. The differences in the barrier height between the standard parameter set ‘0’ and the sets 4A/4B are ca. 0.1 eV. As a result, the rate of transfer changes largely, by a factor of roughly 3 and 10 in AAAA and GAG, respectively. Note that for GAG, the effect of this perturbation is similar to that of the variation of site energy of adenine in 1A/1B, which was however smaller, only 0.05 eV. The energy barriers opposing ET are even more strongly modulated in the simulations 4A/4B, and the huge changes of the transfer rates are reasonable.


Page 46 of 62

Page 47 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The dynamics of the environmental polarization is described by the function frelax (Fig. 12) in the current theoretical model. This function determines how quickly the polarization of environment adjusts to the new charge distribution on the DNA after an ET event has taken place. In the simulation series 4C, 4D and 4E, the time scale of the polarization decay determined by frelax is varied – the relaxation is 2× faster, 2× slower, and 4× slower, respectively. The rate of ET in AAAA is affected strongly by the variation of rate of the environmental relaxation, with faster relaxation leading to increased rate of ET, and the other way around. Note that after a completed ET event, it takes a finite amount of time for the energy barrier opposing a backward ET to recover. If the relaxation of environment is accelerated, this time gets shorter, thus the probability of immediate backward transfer is smaller, resulting in an increased occurrence of a completed ET event. On the other hand, a slower relaxation of the environment provides a potential backward transfer with an extended interval of lower energy barriers, leading to a more efficient backward transfer and a decreased rate of completed ET. In GAG, on the other hand, the dependence of the ET rate on the rate of environmental relaxation is much weaker and possibly inconclusive due to the limited sampling in the simulations.

3.4

Showcase – multiple An bridges

Here, we adopt the construction of DNA ET systems by Giese and Schuster, and study the oligonucleotides GTn GTn GTn GGG (n = 2, 3, denoted as ‘Long2’ and ‘Long3’ in what follows). Simulation of ET in these DNA species is already quite challenging because of the rather large distance and the long time scales that need to be simulated. While the latter can be seen as the consequence of the former, the problem of extremely long simulation time necessary is clearly amplified. Hole transfer was simulated in Long2 for 100 ns, and in Long3 for 200 ns. For each molecule, 128 of such simulations were performed, leading to an accumulated simulation time of 12.8 µs for Long2, and 25.6 µs for Long3. Notably, a 100ns simulation takes 11 days 47 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for Long2, and 15 days for Long3, on one core of an Intel Xeon E5-2630v3 CPU. Still, and most importantly, simulation of ET in such systems on such time scales with our all-atom fragment-orbital QM/MM scheme developed previously would be prohibitively expensive, requiring many months to complete, in spite of the high efficiency of the method in comparison with usual QM/MM approaches. The dynamics of the hole observed in these simulations is presented in Fig. 17. Here, each trace represents the location of the hole along the DNA strand as a function of time, and every graphic contains traces from all (128) of the simulations for the respective DNA species. Apparently, the simulated hundred-nanosecond time scale is sufficient to observe hole transfer along these DNA species. The hole transfer proceeds in individual transfer steps from one guanine to another, a.k.a. G-hopping. (Note that the temporal resolution of the data analyzed here is insufficient to clarify whether the adenine bridge is being occupied during the transfer events or not.) The terminal GGG triplet is reached in several of the 128 simulations, in both Long2 and Long3. Still, in nearly 50 % of the simulations performed, the hole remains located on the initial guanine for the entire simulation time.

Figure 17: Traces from the hole transfer simulations in DNA sequences Long2 (left) and Long3 (right). Every of the 128 lines represent the location of the hole during a single model simulation of 100 ns (Long2) or 200 ns (Long3). The hole was initially located on the outermost G (bottom). If the rate of reaction of guanine radical cation with water was known, it would be possible 48 ACS Paragon Plus Environment

Page 48 of 62

Page 49 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to estimate the ratios of yields of these reactions on the individual guanines. This, however, does not seem to be the case: Giese et al. wrote in 1998 4 that “the reaction rate of G+· with H2 O is not known yet;” Beratan et al. in 2014 22 considered “the trapping rate kr to be ca. 2 ns−1 , consistent with the earlier estimate of Jortner et al. [from 1999]. 10 Giese and Spichty [in 2000] 87 estimated the water-trapping rate to be as low as 0.00006 ns−1 , [from fitting] based on a tunneling rate through two [bridge] sites of 0.0025 ns−1 .” On a qualitative level, we could argue that the (averaged) state of the system after a certain simulation time may correspond to the outcome of a real experiment provided the reaction with water takes a corresponding amount of time. Thus, the observation in a 100ns simulation may be considered a prediction of the outcome of an experiment assuming for the reaction of water a time scale of 100 ns, or, equivalently, a rate of 0.01 ns−1 . Under these conditions, the largest yield of final reaction would be observed on the initial guanine, and gradually decreasing amounts would be present on G2, G3 as well as the terminal GGG.

in 100 ns G2 G3 GGG

Long2 55 % 15 % 5%

Long3 41 % 10 % 2%

Figure 18: Top: Hole transfer in the DNA species Long2 (solid line) and Long3 (broken line) expressed as hole population on the initial guanine G1, relay guanines G2 and G3 as well as the terminal GGG. Bottom: The probability of reaching G2, G3 and GGG in Long2 and Long3 within 100 ns, data obtained as the respective percentages of the performed simulations (128 of Long2 and Long3 each). An alternative way to analyze these data is to consider a kinetic scheme in which the 49 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 50 of 62

population of hole on the initial G decays with first-order kinetics. Then, assuming an exponential decay of population on the initial G, a fit to the data in Fig. 18 provides decay times of 210 and 490 ns for Long2 and Long3, respectively. The number of simulations performed is not large, therefore these figures shall be understood as merely qualitative. Briefly, simulations employing the newly developed theoretical model confirm a G-hopping mechanism and a multi-hundred-nanosecond time scale of the hole transfer reaction in these DNA species.

4

Conclusion

The phenomenon of electron transfer in DNA has attracted a huge amount of experimental, simulation as well as theoretical work. As for the latter, many theoretical models of ET have been proposed, starting actually with the seminal Marcus’ theory, and continuing with more elaborate recent efforts. 88,89 What seems to be common to most, if not all of this previous work is the top-down character, actually starting with an abstract idea of the ET process. A disadvantage of such approaches is that some features of the ET process may be missed. What makes our current approach unique is the completely different underlying idea: The abstraction process starts with a full, atomistic description. Then, chemical knowledge of the molecular system is collected and used to parametrize selected features of the atomistic model. This bottom-up way provides a clear idea of what exactly has been parametrized or possibly neglected. Therefore, the current work represents the best sense of multi-scale modeling. The benefits of the current theoretical model are obvious: Microsecond sampling is possible for the simulation of ET processes in molecular systems containing hundreds of atoms in the ‘quantum’ region, while little accuracy of the original atomistic representation is compromised. The multi-scale character of the scheme is useful in particular because rare events are studied, occurring much less frequently than once per nanosecond, while the underlying


Page 51 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


elementary processes fluctuate on a femto-to-picosecond scale. The current approach can also be used to create parametrized theoretical models of ET in other systems of interest, like proteins or organic semi-conductors. Further, the approach may be applied for simulation of electron transport through molecular junctions. 90 Also, other semi-empirical implementations of surface hopping might benefit from such an extension. 91–93

Acknowledgement This work has been supported by the Deutsche Forschungsgemeinschaft (DFG), joint grant EL 206/13-1 & KL 1299/14-1, as well as by the bwHPC initiative and the bwHPC-C5 project, funded by the state of Baden-W¨ urttemberg and the DFG, through the services of the JUSTUS high performance computing facility located at the University of Ulm.

Supporting Information Available Statistical independence of internal and environmental degrees of freedom; parameters of the site energy function; mean values and variance, as well as correlation coefficients of electronic couplings; distances of neighboring nucleobases; decrease of site energies induced by environmental polarization and the numerical values of the drop function; numerical values of the temporal relaxation function; additional data on the environmental polarization; distribution of site energies for the elementary transfer reaction GGG+ →A→G in the DNA species GTGGG. This material is available free of charge via the Internet at http://pubs. acs.org/.

References (1) Murphy, C. J.; Arkin, M. R.; Jenkins, Y.; Ghatlia, N. D.; Bossmann, S. H.; Turro, N. J.; Barton, J. K. Long-range photoinduced electron transfer through a DNA helix. Science



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1993, 262, 1025–1029. (2) Hall, D. B.; Holmlin, E.; Barton, J. K. Oxidative DNA damage through long-range electron transfer. Nature 1996, 382, 731–735. (3) Henderson, P. T.; Jones, D.; Hampikian, G.; Kan, Y.; Schuster, G. B. Long-distance charge transport in duplex DNA: The phonon-assisted polaron-like hopping mechanism. Proc. Natl. Acad. Sci. USA 1999, 96, 8353–8358. (4) Meggers, E.; Michel-Beyerle, M. E.; Giese, B. Sequence dependent long range hole transport in DNA. J. Am. Chem. Soc. 1998, 120, 12950–12955. (5) Lewis, F. D.; Kalgutkar, R. S.; Wu, Y.; Liu, X.; Liu, J.; Hayes, R. T.; Miller, S. E.; Wasielewski, M. R. Driving force dependence of electron transfer dynamics in synthetic DNA hairpins. J. Am. Chem. Soc. 2000, 122, 12346–12351. (6) Porath, D.; Bezryadin, A.; De Vries, S.; Dekker, C. Direct measurement of electrical transport through DNA molecules. Nature 2000, 403, 635–638. (7) Schuster, G. B., Ed. Long-Range Charge Transfer in DNA I & II ; Topics in Current Chemistry; Springer: Heidelberg, 2004; Vol. 236 & 237. (8) Wagenknecht, H.-A., Ed. Charge Transfer in DNA: From Mechanism to Application; Wiley: Weinheim, 2005. (9) Genereux, J. C.; Barton, J. K. Mechanisms for DNA charge transport. Chem. Rev. 2010, 110, 1642–1662. (10) Bixon, M.; Giese, B.; Wessely, S.; Langenbacher, T.; Michel-Beyerle, M. E.; Jortner, J. Long-range charge hopping in DNA. Proc. Natl. Acad. Sci. USA 1999, 96, 11713–11716. (11) Giese, B.; Amaudrut, J.; Köhler, A. K.; Spormann, M.; Wessely, S. Direct observation of hole transfer through DNA by hopping between adenine bases and by tunnelling. Nature 2001, 412, 318–320. 52 ACS Paragon Plus Environment

Page 52 of 62

Page 53 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12) Berlin, Y. A.; Burin, A. L.; Ratner, M. A. Charge hopping in DNA. J. Am. Chem. Soc. 2001, 123, 260–268. (13) Bixon, M.; Jortner, J. Long-range and very long-range charge transport in DNA. Chem. Phys. 2002, 281, 393–408. (14) Lewis, F. D.; Zhu, H.; Daublain, P.; Fiebig, T.; Raytchev, M.; Wang, Q.; Shafirovich, V. Crossover from superexchange to hopping as the mechanism for photoinduced charge transfer in DNA hairpin conjugates. J. Am. Chem. Soc. 2006, 128, 791–800. (15) O’Neill, M. A.; Becker, H.; Wan, C.; Barton, J. K.; Zewail, A. H. Ultrafast dynamics in DNA-mediated electron transfer: Base gating and the role of temperature. Angew. Chem. Int. Ed. 2003, 42, 5896–5900. (16) O’Neill, M. A.; Barton, J. K. DNA charge transport: Conformationally gated hopping through stacked domains. J. Am. Chem. Soc. 2004, 126, 11471–11483. (17) Conwell, E. M. Charge transport in DNA in solution: The role of polarons. Proc. Natl. Acad. Sci. USA 2005, 102, 8795–8799. (18) Conron, S. M. M.; Thazhathveetil, A. K.; Wasielewski, M. R.; Burin, A. L.; Lewis, F. D. Direct measurement of the dynamics of hole hopping in extended DNA G-tracts. An unbiased random walk. J. Am. Chem. Soc. 2010, 132, 14388–14390. (19) Harris, M. A.; Mishra, A. K.; Young, R. M.; Brown, K. E.; Wasielewski, M. R.; Lewis, F. D. Direct observation of the hole carriers in DNA photoinduced charge transport. J. Am. Chem. Soc. 2016, 138, 5491–5494. (20) Renaud, N.; Berlin, Y. A.; Lewis, F. D.; Ratner, M. A. Between superexchange and hopping: An intermediate charge-transfer mechanism in poly(A)-poly(T) DNA hairpins. J. Am. Chem. Soc. 2013, 135, 3953–3963.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(21) Renaud, N.; Berlin, Y. A.; Ratner, M. A. Impact of a single base pair substitution on the charge transfer rate along short DNA hairpins. Proc. Natl. Acad. Sci. USA 2013, 110, 14867–14871. (22) Zhang, Y.; Liu, C.; Balaeff, A.; Skourtis, S. S.; Beratan, D. N. Biological charge transfer via flickering resonance. Proc. Natl. Acad. Sci. USA 2014, 111, 10049–10054. (23) Liu, C.; Beratan, D. N.; Zhang, P. Coarse-grained theory of biological charge transfer with spatially and temporally correlated noise. J. Phys. Chem. B 2016, 120, 3624–3633. (24) Liu, C.; Xiang, L.; Zhang, Y.; Zhang, P.; Beratan, D. N.; Li, Y.; Tao, N. Engineering nanometre-scale coherence in soft matter. Nat. Chem. 2016, 8, 941–945. (25) Xiang, L.; Palma, J. L.; Bruot, C.; Mujica, V.; Ratner, M. A.; Tao, N. Intermediate tunnelling–hopping regime in DNA charge transport. Nat. Chem. 2015, 7, 221–226. (26) Blumberger, J. Recent advances in the theory and molecular simulation of biological electron transfer reactions. Chem. Rev. 2015, 115, 11191–11238. (27) Levine, A. D.; Iv, M.; Peskin, U. Length-independent transport rates in biomolecules by quantum mechanical unfurling. Chem. Sci. 2016, 7, 1535–1542. (28) Marcus, R. A. On the theory of oxidation-reduction reactions involving electron transfer. I. J. Chem. Phys. 1956, 24, 966–978. (29) Marcus, R. A.; Sutin, N. Electron transfers in chemistry and biology. Biochim. Biophys. Acta 1985, 811, 265–322. (30) Plasser, F.; Lischka, H. Semiclassical dynamics simulations of charge transport in stacked π-systems. J. Chem. Phys. 2011, 134, 034309. (31) Kubaˇr, T.; Kleinekathöfer, U.; Elstner, M. Solvent fluctuations drive the hole transfer in DNA: a mixed quantum–classical study. J. Phys. Chem. B 2009, 113, 13107–13117. 54 ACS Paragon Plus Environment

Page 54 of 62

Page 55 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(32) Kubaˇr, T.; Elstner, M. Coarse-grained time-dependent density functional simulation of charge transfer in complex systems: Application to hole transfer in DNA. J. Phys. Chem. B 2010, 114, 11221–11240. (33) Kubaˇr, T.; Elstner, M. Efficient algorithms for the simulation of non-adiabatic electron transfer in complex molecular systems: application to DNA. Phys. Chem. Chem. Phys. 2013, 15, 5794–5813. (34) Kubaˇr, T.; Elstner, M. A hybrid approach to simulation of electron transfer in complex molecular systems. J. Royal Soc. Interface 2013, 10, 20130415. (35) Woiczikowski, P. B.; Steinbrecher, T.; Kubaˇr, T.; Elstner, M. Nonadiabatic QM/MM simulations of fast charge transfer in Escherichia coli DNA photolyase. J. Phys. Chem. B 2011, 115, 9846–9863. (36) Young, R. M.; Singh, A. P. N.; Thazhathveetil, A. K.; Cho, V. Y.; Zhang, Y.; Renaud, N.; Grozema, F. C.; Beratan, D. N.; Ratner, M. A.; Schatz, G. C. et al. Charge transport across DNA-based three-way junctions. J. Am. Chem. Soc. 2015, 137, 5113– 5122. (37) Newton, M. D.; Boer, F. P.; Lipscomb, W. N. Molecular orbital theory for large molecules. Approximation of the SCF LCAO Hamiltonian matrix. J. Am. Chem. Soc. 1966, 88, 2353–2360. (38) Closs, G. L.; Miller, J. R. Intramolecular long-distance electron transfer in organic molecules. Science 1988, 240, 440–447. (39) Shephard, M. J.; Paddon-Row, M. N.; Jordan, K. D. Electronic coupling through saturated hydrocarbon bridges. Chem. Phys. 1993, 176, 289–304. (40) Kurnikov, I. V.; Beratan, D. N. Ab initio based effective Hamiltonians for long-range electron transfer: Hartree-Fock analysis. J. Chem. Phys. 1996, 105, 9561–9573. 55 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 56 of 62

(41) Newton, M. D. Quantum chemical probes of electron-transfer kinetics: The nature of donor-acceptor interactions. Chem. Rev. 1991, 91, 767–792. (42) Kubas, A.; Hoffmann, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic couplings for molecular charge transfer: Benchmarking CDFT, FODFT, and FODFTB against high-level ab initio calculations. J. Chem. Phys. 2014, 140, 104105. (43) Kubaˇr, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. Efficient calculation of charge-transfer matrix elements for hole transfer in DNA. J. Phys. Chem. B 2008, 112, 7937–7947. (44) Kubas, A.; Gajdos, F.; Heck, A.; Oberhofer, H.; Elstner, M.; Blumberger, J. Electronic couplings for molecular charge transfer: benchmarking CDFT, FODFT and FODFTB against high-level ab initio calculations. II. Phys. Chem. Chem. Phys. 2015, 17, 14342– 14354. (45) Heck, A.; Woiczikowski, P. B.; Kubaˇr, T.; Welke, K.; Niehaus, T.; Giese, B.; Skourtis, S.; Elstner, M.; Steinbrecher, T. B. Fragment orbital based description of charge transfer in peptides including backbone orbitals. J. Phys. Chem. B 2014, 118, 4261–4272. (46) Kubaˇr, T.; Elstner, M. What governs the charge transfer in DNA? The role of DNA conformation and environment. J. Phys. Chem. B 2008, 112, 8788–8798. (47) https://www2.chemistry.msu.edu/faculty/reusch/virttxtjml/Spectrpy/InfraRed/infrared.htm (last accessed on 22 December 2016). (48) Voityuk, A. A.; Siriwong, K.; Rösch, N. Environmental fluctuations facilitate electronhole transfer from guanine to adenine in DNA π stacks. Angew. Chem. Int. Ed. 2004, 43, 624–627. (49) Steinbrecher, T.; Koslowski, T.; Case, D. A. Direct simulation of electron transfer reactions in DNA radical cations. J. Phys. Chem. B 2008, 112, 16935–16944. 56 ACS Paragon Plus Environment

Page 57 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50) Hatcher, E.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Beratan, D. N. PNA versus DNA: Effects of structural fluctuations on electronic structure and hole-transport mechanisms. J. Am. Chem. Soc. 2008, 130, 11752–11761. (51) Newton, M. D. Extension of Hopfield’s electron transfer model to accommodate site–site correlation. J. Phys. Chem. B 2015, 119, 14728–14737. (52) Hopfield, J. J. Electron-transfer between biological molecules by thermally activated tunneling. Proc. Natl. Acad. Sci. USA 1974, 71, 3640–3644. (53) Woiczikowski, P. B.; Kubaˇr, T.; Gutiérrez, R.; Caetano, R. A.; Cuniberti, G.; Elstner, M. Combined density functional theory and Landauer approach for hole transfer in DNA along classical molecular dynamics trajectories. J. Chem. Phys. 2009, 130, 215104. (54) Senthilkumar, K.; Grozema, F. C.; Guerra, C. F.; Bickelhaupt, F. M.; Lewis, F. D.; Berlin, Y. A.; Ratner, M. A.; Siebbeles, L. D. A. Absolute rates of hole transfer in DNA. J. Am. Chem. Soc. 2005, 127, 14894–14903. (55) Balabin, I. A.; Beratan, D. N.; Skourtis, S. S. Persistence of structure over fluctuations in biological electron-transfer reactions. Phys. Rev. Lett. 2008, 101, 158102. (56) Beratan, D. N.; Skourtis, S. S.; Balabin, I. A.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Xiao, D. Steering electrons on moving pathways. Acc. Chem. Res. 2009, 42, 1669–1678. (57) Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Fluctuations in biological and bioinspired electron-transfer reactions. Annu. Rev. Phys. Chem. 2010, 61, 461–485. (58) Beratan, D. N.; Liu, C.; Migliore, A.; Polizzi, N. F.; Skourtis, S. S.; Zhang, P.; Zhang, Y. Charge transfer in dynamical biosystems, or the treachery of (static) images. Acc. Chem. Res. 2015, 48, 474–481.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(59) Uhlenbeck, G. E.; Ornstein, L. S. On the theory of the Brownian motion. Phys. Rev. 1930, 36, 823–841. (60) Fox, R. F.; Gatland, I. R.; Roy, R.; Vemuri, G. Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise. Phys. Rev. A 1988, 38, 5938–5940. (61) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (62) Olofsson, J.; Larsson, S. Electron hole transport in DNA. J. Phys. Chem. B 2001, 105, 10398–10406. (63) Furse, K. E.; Corcelli, S. A. Molecular dynamics simulations of DNA solvation dynamics. J. Phys. Chem. Lett. 2010, 1, 1813–1820. (64) Kuharski, R. A.; Bader, J. S.; Chandler, D.; Sprik, M.; Klein, M. L.; Impey, R. W. Molecular model for aqueous ferrous-ferric electron transfer. J. Chem. Phys. 1988, 89, 3248–3257. (65) King, G.; Warshel, A. Investigation of the free energy functions for electron transfer reactions. J. Chem. Phys. 1990, 93, 8682–8692. (66) Schulten, K.; Tesch, M. Coupling of protein motion to electron transfer: Molecular dynamics and stochastic quantum mechanics study of photosynthetic reaction centers. Chem. Phys. 1991, 158, 421–446. (67) Alden, R. G.; Parson, W. W.; Chu, Z. T.; Warshel, A. Calculations of electrostatic energies in photosynthetic reaction centers. J. Am. Chem. Soc. 1995, 117, 12284–12298. (68) Blumberger, J.; Lamoureux, G. Reorganization free energies and quantum corrections for a model electron self-exchange reaction: Comparison of polarizable and nonpolarizable solvent models. Mol. Phys. 2008, 106, 1597–1611.


Page 58 of 62

Page 59 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(69) Marchi, M.; Gehlen, J. N.; Chandler, D.; Newton, M. Diabatic surfaces and the pathway for primary electron transfer in a photosynthetic reaction center. J. Am. Chem. Soc. 1993, 115, 4178–4190. (70) Warshel, A.; Hwang, J.-K. Simulation of the dynamics of electron transfer reactions in polar solvents: Semiclassical trajectories and dispersed polaron approaches. J. Chem. Phys. 1986, 84, 4938–4957. (71) Sulpizi, M.; Raugei, S.; VandeVondele, J.; Carloni, P.; Sprik, M. Calculation of redox properties: Understanding short- and long-range effects in rubredoxin. J. Phys. Chem. B 2007, 111, 3969–3976. (72) Leontyev, I. V.; Stuchebrukhov, A. A. Electronic continuum model for molecular dynamics simulations of biological molecules. J. Chem. Theory Comput. 2010, 6, 1498– 1508. (73) Leontyev, I. V.; Stuchebrukhov, A. A. Accounting for electronic polarization in nonpolarizable force fields. Phys. Chem. Chem. Phys. 2011, 13, 2613–2626. (74) Vazdar, M.; Pluhaˇrová, E.; Mason, P. E.; Vácha, R.; Jungwirth, P. Ions at hydrophobic aqueous interfaces: Molecular dynamics with effective polarization. J. Phys. Chem. Lett. 2012, 3, 2087–2091. (75) Tully, J. C.; Preston, R. K. Trajectory surface hopping approach to nonadiabatic molecular collisions: The reaction of H+ with D2 . J. Chem. Phys. 1971, 55, 562–572. (76) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061–1071. (77) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. J. Chem. Phys. 2001, 114, 10608–10615.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(78) Granucci, G.; Persico, M. Critical appraisal of the fewest switches algorithm for surface hopping. J. Chem. Phys. 2007, 126, 134114. (79) Zhu, C.; Jasper, A. W.; Truhlar, D. G. Non-Born-Oppenheimer Liouville-von Neumann dynamics. evolution of a subsystem controlled by linear and population-driven decay of mixing with decoherent and coherent switching. J. Chem. Theory Comput. 2005, 1, 527–540. (80) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent progress in surface hopping: 2011–2015. J. Phys. Chem. Lett. 2016, 7, 2100–2112. (81) Jain, A.; Subotnik, J. E. Surface hopping, transition state theory, and decoherence. II. thermal rate constants and detailed balance. J. Chem. Phys. 2015, 143, 134107. (82) Doltsinis, N. L. In Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms; Grotendorst, J., Marx, D., Muramatsu, A., Eds.; NIC Series; John von Neumann Institute for Computing: J¨ ulich, 2002; Vol. 10; Chapter Nonadiabatic Dynamics: Mean-Field and Surface Hopping, pp 377–397. (83) Takada, T.; Kawai, K.; Cai, X.; Sugimoto, A.; Fujitsuka, M.; Majima, T. Charge separation in DNA via consecutive adenine hopping. J. Am. Chem. Soc. 2004, 126, 1125–1129. (84) LeBard, D. N.; Lilichenko, M.; Matyushov, D. V.; Berlin, Y. A.; Ratner, M. A. Solvent reorganization energy of charge transfer in DNA hairpins. J. Phys. Chem. B 2003, 107, 14509–14520. (85) Renger, T.; Marcus, R. A. Variable-range hopping electron transfer through disordered bridge states: Application to DNA. J. Phys. Chem. A 2003, 107, 8404–8419. (86) Wang, L.; Trivedi, D.; Prezhdo, O. V. Global flux surface hopping approach for mixed quantum-classical dynamics. J. Chem. Theory Comput. 2014, 10, 3598–3605. 60 ACS Paragon Plus Environment

Page 60 of 62

Page 61 of 62

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(87) Giese, B.; Spichty, M. Long distance charge transport through DNA: Quantification and extension of the hopping model. ChemPhysChem 2000, 1, 195–198. (88) Cuniberti, G.; Maciá, E.; Rodriguez, A.; Römer, R. A. In Charge migration in DNA: perspectives from Physics, Chemistry and Biology; Chakraborty, T., Ed.; Springer, Berlin, 2007; p 1. (89) Gutiérrez, R.; Cuniberti, G. In NanoBioTechnology: BioInspired device and materials of the future; Shoseyov, O., Levy, I., Eds.; Humana Press, 2007; p 107. (90) Kubaˇr, T.; Gutiérrez, R.; Kleinekathöfer, U.; Cuniberti, G.; Elstner, M. Modeling charge transport in DNA using multi-scale methods. Phys. Stat. Sol. B 2013, 250, 2277–2287. (91) Ren, J.; Vukmirović, N.; Wang, L.-W. Nonadiabatic molecular dynamics simulation for carrier transport in a pentathiophene butyric acid monolayer. Phys. Rev. B 2013, 87, 205117. (92) Pal, S.; Trivedi, D. J.; Akimov, A. V.; Aradi, B.; Frauenheim, T.; Prezhdo, O. V. Nonadiabatic molecular dynamics for thousand atom systems: A tight-binding approach toward PYXAID. J. Chem. Theory Comput. 2016, 12, 1436–1448. (93) Spencer, J.; Gajdos, F.; Blumberger, J. FOB-SH: Fragment orbital-based surface hopping for charge carrier transport in organic and biological molecules and materials. J. Chem. Phys. 2016, 145, 064102.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic


Page 62 of 62

Microsecond Simulation of Electron Transfer in DNA: Bottom-Up

Recommend Documents