Charge Transfer in E. coli DNA Photolyase: Understanding

Article pubs.acs.org/JPCB

Charge Transfer in E. coli DNA Photolyase: Understanding Polarization and Stabilization Effects via QM/MM Simulations Gesa Lüdemann, P. Benjamin Woiczikowski, Tomás ̌ Kubař, Marcus Elstner, and Thomas B. Steinbrecher* Department for Theoretical Chemical Biology, Institute for Physical Chemistry, Karlsruhe Institute for Technology, Kaiserstr. 12, 76131 Karlsruhe, Germany S Supporting Information *

ABSTRACT: We study fast hole transfer events in E. coli DNA photolyase, a key step in the photoactivation process, using a multiscale computational method that combines nonadiabatic propagation schemes and linear-scaling quantum chemical methods with molecular mechanics force fields. This scheme allows us to follow the time-dependent evolution of the electron hole in an unbiased fashion; that is, no assumptions about hole wave function localization, time scale separation, or adiabaticity of the process have to be made beforehand. DNA photolyase facilitates an efficient long-range charge transport between its flavin adenine dinucleotide (FAD) cofactor and the protein surface via a chain of evolutionary conserved Trp residues on the sub-nanosecond time scale despite the existence of multiple potential trap states. By including a large number of aromatic residues along the charge transfer pathway into the quantum description, we are able to identify the main pathway among alternative possible routes. The simulations show that charge transfer, which is extremely fast in this protein, occurs on the same time scale as the protein response to the electrostatic changes; that is, time-scale separation as often presupposed in charge transfer studies seems to be inappropriate for this system. Therefore, coupled equations of motion, which propagate electrons and nuclei simultaneously, appear to be necessary. The applied computational model is shown to capture the essentials of the reaction kinetics and thermodynamics while allowing direct simulations of charge transfer events on their natural time scale.

■

INTRODUCTION Among the countless biochemical charge transfer (CT) systems of high scientific interest, the photosensitive enzymes of the photolyase/cryptochrome family continue to attract considerable attention from both experimentalists and theoreticians. Photolyases (PL) are an ancient class of blue-light activated flavoproteins that repair radiation-damaged DNA in a catalytic CT reaction.1−3 The PL enzymes are divided into the two classes CPD (cyclobutane pyrimidine dimers) and 6-4 (pyrimidine-[6,4]-pyrimidone) PL, depending on the types of UV damage they repair. PL do not play a role in human (or other mammalian) DNA repair, but cryptochromes are found in animals (including humans), as well as plants and bacteria.4 The biochemical function of PL has been known for a long time5,6 and especially over the past decade, the mechanism of PL and related enzymes has been studied in detail.7−12 Beyond the ability to repair DNA, another mechanism of action, the photoactivation process, has been found in which an electron is transferred along a chain of evolutionary conserved tryptophan residues to rereduce the flavin cofactor. It remains unclear if this process is physiologically relevant in PL. Still, it has been suggested that it is functionally important in signaling cryptochromes.13,14 The photoactivation and DNA repair processes are difficult cases for classical theory of CT and have attracted considerable © 2013 American Chemical Society

interest from theoreticians (reviewed in ref 15, and see refs 16−27 for newer work). With CT in photoactivation occurring over nanometer distances and nanosecond times,28−30 direct, atomistic QM/MM simulations are powerful enough to explicitly simulate this reaction in PL.31 Modeling of CT in complex environments is commonly based on the seminal Marcus theory of CT rates32−35 and its extensions,36−40 since the involved time and length scales make it impossible to use time dependent ab initio wave function approaches, such as nonadiabatic molecular dynamics of the entire system, on solvated biomolecules. Many applications use the Marcus formula to compute transfer rates between donor and acceptor sites. Since the molecular environment in biomolecules is crucial to understanding the CT process, theoreticians quickly turned to molecular modeling efforts to better understand the atomistic details of these processes.23,41−51 It is thus possible to estimate the Marcus parameters from MD simulations, although their absolute values are still difficult to determine in sufficient accuracy for quantitative rate calculations. Received: June 26, 2013 Revised: August 20, 2013 Published: August 21, 2013 10769

dx.doi.org/10.1021/jp406319b | J. Phys. Chem. B 2013, 117, 10769−10778

The Journal of Physical Chemistry B

Article

One recent development is a stronger focus on the coupling between atomic conformational fluctuations and electronic structure changes in biomolecules,52−56 which has grown out of earlier work on model systems.57−62 Energetic and structural fluctuations as well as their relative time scales play a vital role for system properties at the transition state and characterize different CT mechanisms like gated CT.63−65 The application of the basic Marcus formula reaches a limit in cases where electron density rearrangements occur on the same time scale as environmental relaxation, as described recently for the initial electron transfer step in the bacterial reaction center.66,67 While solvent relaxation times often seem to be fast compared to hole transfer, as, for example, for the case of DNA,68,69 it has recently been shown that the response of coupled hydration and protein conformational dynamics to sudden electrostatic changes of chromophores inside a protein can range from few picoseconds to several hundreds of picoseconds.70 This is on a comparable time scale to the hole transfer in photolyase, and our recent work indicates that such an effect may be important to understand the fast CT process.31 In this work, we will address this question in more detail and in particular examine effects not considered in our earlier work. We investigate the role of residues near the main CT pathway, which could in principle be intermediate charge carriers, as well as the role of protein electronic polarization, which can have a large impact on CT efficiency. For this purpose, we apply our recently developed CT methodology based on a combination of nonadiabatic propagation schemes, linear-scaling quantum chemical methods, and molecular mechanics force fields.69,71,72 Our wave function based nonadiabatic propagation method offers an alternative way to study nonadiabatic processes in physics, chemistry, and biology. As usual in this approach, we use coupled equations of motion where the electronic degrees of freedom are propagated simultaneously to the classical nuclear ones. A dramatic speedup is gained by introducing several approximations, which make a careful testing before application mandatory. In principle, we achieve an unbiased treatment of the CT process similar to the ab initio techniques from which our methodology is derived, provided that additional approximations are well controlled and that the QM region has an appropriate extension.

treated quantum mechanically, while the remainder of the total energy is represented by force field terms. The total hole (electron) wave function is represented as a superposition of the HOMO orbitals ϕHO m at the fragments m, which are the selected amino acid side chains, namely, tryptophans and tyrosines in this application, as shown in Figure 1. Ψ(t ) =

∑ am(t )ϕmHO m

(1)

Figure 1. Fragment orbital QM scheme. Shown are the HOMO orbitals of the three most important tryptophan fragments in color, embedded in the classical MM surroundings of the PL protein, FAD cofactor, and solvent (gray). Figure was prepared using vmd99,100 and povray.101

In the case of PL this constitutes a good approximation since the QM system is built up from weakly interacting fragments; that is, the electronic couplings between the sites (amino acid side chains) are small.71,72 In a next step, we compute the fragment molecular orbitals (FMOs) ϕmHO applying an electronic structure method. For every fragment (tryptophan, tyrosine), an individual SCF calculation is conducted using a QM/MM formalism; that is, the FMOs are polarized due to the inclusion of the MM point charges of the environment (protein, cofactors, and water). In principle, this could be done using DFT or HF methods, but the application of the approximate DFTB2 (SCC-DFTB)80 model allows the computation of these FMOs roughly 2−3 orders of magnitude faster than with standard ab initio methods, without a significant loss in accuracy. For example, a previous benchmarking study showed that DFTB reproduces experimental IP values for nucleobases with a standard deviation of 0.26 eV. Relative IP values for DNA bases were reproduced with a standard deviation of 0.23 eV, compared to 0.09 eV for DFT. Furthermore, DFTB yields energetic fluctuations over time that agree well with reference calculations.79 Similar tests have been performed for tryptophans for the case of PL.31 Now, we can represent the electronic Hamiltonian of the system Ĥ in this FMO-basis, which results in an m × m matrix if m amino acid side chains are included in the CT pathway:

■

MODELS AND METHODS Nonadiabatic CT processes can be captured by electronic structure theory approaches which go beyond the Born− Oppenheimer (BO) approximation. Built onto, for example, the density functional theory (DFT) framework, two approaches are frequently applied,73 the mean-field Ehrenfest scheme74−76 and the surface hopping (SH) algorithm.77,78 Our methodology is an approximation to such nonadiabatic DFT approaches and has been reviewed extensively,69,71,72,79 and therefore, we give only a very brief introduction. In a first step, a QM/MM partitioning is applied where the QM part includes those residues being active in CT while the remainder of the systems is treated by the MM force field. The QM description is very efficient, since (a) a linear scaling fragmentation technique is used and (b) a partitioning of the total QM energy is applied, where only the frontier orbitals are treated quantum mechanically (QM) while the remainder of a fragment is described using molecular mechanics (MM). This is similar in spirit to the Hückel or Pariser−Parr−Pople models, where, for example, only one p-orbital per carbon atom is

0 Hmn = ⟨ϕmHO|Ĥ |ϕnHO⟩

(2)

H0mn

Note that the are affected by the electrostatic potential of the environment due to the inclusion of QM/MM terms in the 10770



Article

Here, EQM stands for the single point energy of the neutral (F) or charged (F+) fragment in brackets, using the optimized geometry of the fragment charge state in subscript. Geometry optimization and energy calculations were conducted using DFT with the B3LYP/cc-pVTZ level of theory. Values for λi of 0.21 and 0.18 eV were obtained for Tyr and Trp, respectively The big advantage of the present formalism is the derivation of an expression for the total energy (eq 5). This allows the computation of nuclear forces, which can be used to drive Born−Oppenheimer molecular dynamics (BO-MD) or coupled electron-ion dynamics, where the electrons move according to the time dependent Schrödinger equation and the nuclei according to classical Newtonian equations of motion. The H0mn in eq 2 are computed for the neutral system, while the Hmn for the N − 1 electron system result from the H0mn by including the influence of the excess charge on the matrix elements in second order perturbation theory.72 During the simulations, the force field charges on the MM subsystem q0A are fixed, while the partial charges on atoms in the QM region qα are updated every MD step to include the contribution of the moving excess charge. The whole framework has been derived for fixed MM point charges, that is, nonpolarizable empirical force fields. As discussed in our previous work, the effects of electronic polarization are quite important and should be considered.69,72 The most simple, although approximate, way to include these effects is to scale the electrostatic potential which is applied in the QM/MM calculation of the H0mn matrix elements. We will use a scaling factor εMM of 1.0, 1.4, and 2.0, respectively, in the calculations described below. The external potential is divided by εMM, leading to a reduced environmental electrostatics effect for larger εMM-values. It has been shown repeatedly that standard force fields using a fixed charge distribution neglecting electronic polarization strongly impact computed CT and other molecular properties.43,49,82−84 Particularly the reorganization energy λ is affected by electronic polarization, while ΔG0 seems to be less sensitive.49 λ can be overestimated by 30−50% when computed with a fixed charge model. Therefore, it was suggested that estimates of λ from MM models need to be scaled by the high-frequency dielectric constant of the medium,47 which can be rationalized by considering the Pekar factor of the Marcus continuum formula.49 The alternative of performing MD simulations utilizing a variable dipole, polarizable force field could be considered a more exact solution but incurs computational costs too high for the large systems and long simulations required here. It is difficult to estimate the effect of electronic polarization on CT reaction rates, since including it in simulations decreases both energy differences between sites arising from environmental stabilization and fluctuations in site energy levels.85 The former effect would tend to increase rates (by reducing static disorder of site energies) while the latter can reduce them (by damping dynamic disorder, an important driver of dynamics). Recently, Leontyev and Stuchebrukhov86 have proposed to scale the charge of charged residues and ions by a factor of roughly 0.7, since the unscaled treatment of charges in MM force fields may be considered to be inconsistent. This procedure is novel and not firmly established yet, but we point out that the introduction of a scaling factor of 1.4, as discussed above, is completely in line with this reasoning. Effectively, the bare hole charge is reduced by a factor of ca. 0.7 due to this scaling. In our context, the scaling factor is

calculation of the FMOs. The electron hole then propagates according to the electronic time dependent Schrödinger equation via the time-dependent orbital expansion coefficients am in eq 1. The “coarse grained” (CG) molecular Hamiltonian, H0mn is time dependent through the motion of the nuclei. This is due first to the motion of the atoms within the fragment, but even more importantly, a large portion of the fluctuations results from the environment through the QM/MM coupling. These fluctuations have been investigated in detail for DNA in, for example, ref 81 and in ref 31 for PL. In the formulation up to now, the Hamiltonian does not depend on the excess charge in the system. To derive coupled equations of motion for the electronic and nuclear systems, we start from the total energy functional of DFT and consider a singly ionized system; that is, we remove one electron from the N electron system. Then, the total energy of the N − 1 electron system with density ρ is expanded up to second order around the charge neutral one with density ρ0, arriving at E N − 1[ρ] = E 0[ρ0 ] + E1[ρ0 , δρ] + E2[ρ0 , (δρ)2 ]

(3)

We assume the orbitals of the charged system Ψi to be identical to those of the neutral system Ψ0i . This approximation may appear severe at first, but it is typically valid (with orbital overlap between Ψi and Ψ0i of 0.99 and higher) for the HOMO orbitals of extended π-systems, as discussed in detail in ref 72, especially Figure 3 therein and accompanying discussion. We may now rearrange the energy expression in second order to obtain E N − 1[ρ] ≈ E N [ρ0 ] − ⟨Ψ 0N /2|H[ρ0 ]|Ψ 0N /2⟩ + E2nd[(δρ)2 ] (4) N

which has the form of eq 3. E [ρ0] is the total energy of the charge neutral system, at this stage given in the DFT formalism. In the following, we substitute it by the total energy computed using an empirical force field, EN[ρ0] ≈ EMM. For the second term in eq 4, we insert the expansion eq 1 and for the third term we apply a monopole expansion of the charge density fluctuations.71 With these approximations, the total energy (eq 4) reads: E N − 1[ρ] ≈ E MM −

0 − ∑ am*anHmn mn

1 C ∑ Q Q Γmn 2 mn m n (5)

The Qm represent the fraction of the hole charge on site m and Um = Γmm is the chemical hardness of site m (UW = 5.33 eV/e2 for tryptophan and UY = 5.86 eV/e2 for tyrosine), while Γmn = 1/Rmn represents the Coulomb interaction of the excess charge between the sites m and n for the off-diagonal parts. C is a constant which accounts for the DFT delocalization error; that is, it represents an efficient treatment for the DFT selfinteraction error.71 Note that due to the force field treatment the fragment equilibrium geometry does not change upon charging. Therefore, we add the inner sphere contribution to the reorganization energy λi to the second order term as well, as discussed in detail in detail in refs 31, 69, 71, and 72. In brief, λi is incorporated as an effective parameter, determined from ab initio QM calculations of the individual fragments: QM QM + λi = EFQM(F +) + E FQM (F ) + (F ) − E + (F ) − EF F

(6) 10771



Article

motivated by the otherwise dramatic overestimation of the solvent (outer sphere) reorganization energy λS, which is in line with the analysis of ref 86. Molecular Dynamics Simulations and Model Building. The protein structure model has been prepared as described before,31 based on the widely used15,21 X-ray crystal structure of E. coli DNA photolyase by Park et al.87 The improved Amber99 force field was used for the protein.88,89 The FAD cofactor was built from riboflavin and ADP models, using force field parameters from ref 90. for the latter and gaff parameters91 for the former as well as for the light-harvesting cofactor MTHF. Atomic charges for cofactors in all required charge states were determined from RESP calculations.92 Charges for a radicalcation Tyr side chain species were obtained from calculation on a 4-methylphenol (para-cresol) model while maintaining consistency with existing amino acid partial charge parametrization, as described before.31 Additional force field parameters for cationic species remained unchanged compared to the neutral moieties. Note that the internal relaxation of side chain geometries due to their charge state; that is, the internal reorganization energy λi, is accounted for as a separate parameter determined by ab initio quantum chemical calculations (see above). The protein model was solvated in a rectangular periodic boundary box of ca. 10 × 10 × 10 nm3 size filled with TIP3P model water molecules.93 All simulations were carried out using the GROMACS modeling package94 with modifications to carry out coupled QM/MM calculations as described above. The LINCS algorithm95 was employed to keep bonds involving hydrogen at fixed length, and a time step of 1 fs was used to account for the effect of changing charges.

Figure 2. Structure of E. coli DNA photolyase. Left: Complete protein based on the X-ray crystal structure with PDB code 1DNP (cartoon representation with color according to secondary structure elements). The enzyme has a two-domain structure, with an N-terminal α/βdomain and an α-helical FAD binding domain. Right: Details of the CT system, including antenna cofactor MTHF (gray), catalytic cofactor FAD (orange), established CT cascade residues Trp306, Trp359, and Trp382 (red), and potential additional CT active amino acid side chains Trp 316, Trp384, and Tyr464 (light red).

potential fourth hole-carrying residue. Finally, Tyr464 could be suggested as an alternative end point for the CT process if it is sufficiently stabilized by the environment to overcome its unfavorable redox potential. One goal of the presented study is to determine how likely it is that these three additional residues influence the electron transfer in photolyase photoactivation. The average center of mass distance between all six QM residues is summarized in Table 1. We have selected only

■

RESULTS Photolyase Structural Dynamics and QM Region Selection. Starting structures for all calculations presented were taken from previously described31 30 ns long classical MD simulations, in which the pre-equilibrated photolyase enzyme exhibited only little conformational fluctuations, using a total of 29 structural snapshots taken at 1 ns time intervals to run independent QM/MM simulations. For the classical simulation, the charge state of the protein was set to a neutral flavin molecule as well as neutral Trp and Tyr residues, corresponding to the state before the initial charge abstraction that starts the photoactivation process. The 29 snapshots selected represent a typical ensemble of photolyase structural variety, with an average Cα-RMSD value of 0.8 Å compared to their common average structure. Each individual simulation, independent of the QM region selected and the details of the observed CT process, remained structurally stable and only small conformational changes are observed (with Cα-RMSD values for all 29 QM/MM simulations in the range of 0.8 to 1.5 Å). For the QM region, that is, the selection of amino acid side chains the electron hole is able to occupy, six residues were chosen: the three evolutionary conserved Trp residues that make up the established electron transfer chain (Trp382, Trp359, and Trp306) as well as three additional residues close enough to the first three to potentially play a role in the CT process: Trp384, Trp316, and Tyr464. See Figure 2 for a description of the QM region. Among the three additional residues, Trp384 lies close to both the FAD cofactor and Trp382 and could be involved in the initial stages of CT. Trp316 is close enough to both Trp382 and Trp359 to act as a

Table 1. Average Distance between All QM Region Residues From 1 ns Length MD Simulationsa 382 384 316 359 306 a

Trp384

316

359

306

464

1.077

0.735 1.587

0.806 1.294 1.374

1.158 1.634 1.374 0.727

1.574 1.437 1.796 1.282 0.903

All distances are given in nm.

residues containing extended π-systems for the QM region. A larger QM zone, including, for example, peptide bond πsystems, could have been selected, but we have restricted it to Tyr and Trp side chains which are most likely to form radical cations, to limit computational costs. Influence of Additional QM Sites. The ensemble of photolyase structure models described above was used as starting point for 29 unrestrained fully coupled QM/MM simulations of 1 ns length each. The occupation (|am|2, see eq 1) of each of the six potentially charge-carrying QM sites was recorded for each time step of the simulations. Figure 3 shows a summary of the site occupations. A similar picture to that of previous simulations, where the charge was constrained to only three side chains, emerges: Trp382, Trp359, Trp306, and to a small degree Trp316 are involved in the CT process. Neither Trp384, where charge could localize transiently in the initial stages of electron transfer, nor the alternative end point Tyr464 are occupied by the positive electron hole during any of the simulations. 10772



Article

Figure 4. Averages of IP of QM sites and site−site electronic couplings over 29 simulations. Only electronic couplings larger than 0.5 meV are shown. Values of IP (color as in Figure 3) are influenced significantly by the dynamics of the moving charge, indicating charge stabilization via polarization of the environment. In contrast, electronic couplings mainly fluctuate rapidly around unchanging average values. The clear energetic separation between the charge carrying and non charge carrying sites is evident, as well as the nature of Trp306 as a stable end state of the CT process.

Figure 3. Dynamic electron hole occupation and average occupation over 29 simulations for four of the selected QM sites. Data for Trp384 and Tyr 464 is not plotted, since for these residues no significant charge occupancy was observed in any simulation. The plots for all individual simulations (top, different colors correspond to different QM/MM MD trajectories) show the large variability of CT in this system. The average site occupation plot (bottom) suggests an overall two-step kinetics process, with charge flowing from Trp382 to transiently occupy Trp359 and finally to Trp306 as the end point of CT.

the location of the electron hole (Figure 4). Initially, Trp382, the CT starting point, has the lowest IP, which quickly rises by ca. 1 eV during the initial 100 ps of the simulations. The IP values for Trp359 and Trp306 drop significantly during this time, with the IP of Trp306 decreasing continuously until it has dropped by ca. 2 eV after 500 ps. This is consistent with the picture of fast initial CT from Trp382 to Trp359 and further to Trp306, with environmental polarization by the moving charge leading to a significant drop in IP and dynamic charge stabilization. Interestingly, Trp359 was never found to have the lowest IP of the QM sites in the average plots, even though it does so during many individual simulations when the positive charge localizes on it. The large error bars indicate the large variety of the underlying MD simulation results. The non charge carrying sites Trp384 and Trp316, which are located in the protein center, experience a similar increase in their IP values as Trp382, as the charge moves toward the protein surface, while the IP of Tyr464 drops concurrently with that of Trp306. Still, the energetic separation between Tyr464 and Trp306 remains large, since environmental stabilization of the harder to oxidize Tyr side chain is not sufficient to allow charge localization on it. The FMO approach in eq 1 relies on electronic couplings between fragment orbitals to be small, compared to their energetic separation. The data in Figure 4 shows this approximation to be valid, with couplings in the 10 meV range, compared to IP differences of several hunded meV or more. To connect the calculated occupation trajectories to experimental data, we have conducted a numerical fit to kinetic models (see the Supporting Information for details). Fitting was performed both with and without allowing backward reactions and the model including backward reactions yielded significantly better fits. The rate constants obtained are given in Table 2. The first fit (top block in Table 2) disregards the

The plots for all individual simulations (Figure 3 top, different colors correspond to different QM/MM MD trajectories) show the large variability of CT in this system. The average site occupation plot (Figure 3 bottom) suggests an overall two-step kinetics process, with charge flowing from Trp382 to transiently occupy Trp359 and finally to Trp306 as end point of the CT. The combined plot of all individual CT simulations shows the large variability from simulation to simulation underlying the average occupation plots. The positive charge partially delocalizes to Trp316, mainly in the initial 250 ps of a few of the simulations, but for most of them Trp316 is rarely occupied to more than 10%. In contrast, Trp382, Trp359, and Trp306 are regularly occupied by the majority (>90%) of the charge. Localization on either of the three sites is often only very brief, in the tens of picoseconds range, indicating frequent forward and backtransfers, especially between Trp382 and Trp359. During the final stages of the electron transfer, after ca. 800 ps, the charge predominantly localizes on Trp306 with only occasional, short back-transfers to Trp359. The average characteristic CT parameters show electronic couplings between sites as only weakly affected by the nuclear and electron hole dynamics (Figure 4). The coupling between Trp382 and Trp359 steadily rises from ca. 7 meV at simulation start to ca. 10 meV after 1 ns. The electronic couplings between Trp382 and Trp 316, as well as between Trp359 and Trp306 are fluctuating around their respective average values of 4.98 ± 0.68 and 3.52 ± 0.51 meV. The only other nonzero coupling is that between Trp306 and Tyr464, at 0.26 ± 0.17 meV. In contrast, the average IP values vary significantly, depending on 10773



Article

Table 2. Charget Transfer Rates Using Different Kinetic Models and Polarization Parametersa donor Trp382 Trp359 Trp382 Trp359 Trp382 Trp316 fit Trp382 Trp359 fit Trp382 Trp359

acceptor

kF

kR

fit including Trp382, Trp359, and Trp306 ⇌ Trp359 34.8 18.3 ⇌ Trp306 43.3 7.4 fit including Trp382, Trp359, Trp316, and Trp306 ⇌ Trp359 32.4 5.0 ± 3.3 ⇌ Trp306 47.6 7.7 ⇌ Trp316 8.9 ± 1.8 72.4 ± 21.1 ⇌ Trp359 0.0 9.8 ± 3.4 including Trp382, Trp359, and Trp306 with εMM = 1.4 ⇌ Trp359 88.3 86.7 ⇌ Trp306 64.3 6.3 including Trp382, Trp359, and Trp306 with εMM = 2.0 ⇌ Trp359 179.7 125.4 ⇌ Trp306 162.5 29.6

Figure 5. Average site occupation over 22 simulations for three QM sites under varying polarization scaling factors εMM. Scaling factors are 1.0, 1.4, and 2.0 from top to bottom. Increase in εMM results in faster transfer and increased fluctuations of the average occupation values.

2.0 gives a 6-fold increase in kTrp382→Trp359 and a 4-fold one in kTrp359→Trp306. The backward reaction rates also rise significantly. As shown in ref 31 (Figure 6 right), by removing the solvent contribution from the QM/MM term, the reaction photoactivation CT reaction becomes an endothermic, uphill process when only the interaction with the protein is considered. This can be explained by the fact that the charge is better stabilized inside the protein, i.e. at the first tryptophan, than at the outside when the solvent is neglected. Therefore, the exergonicity of the CT reaction is mainly established by the solvation contribution to ΔG0. As the polarization scaling reduces the values of ΔG0 as well as the reorganization energy λS, the energy difference between the sites is reduced (see Figure 6). With larger scaling factors, the average energetic separation between QM sites drops. The stabilization of the sites due to the solvent is reduced, leading to faster transfers. However, the reduced energetic separation also enables back-transfer more easily. This makes the CT more facile and results in a slightly larger electron hole delocalization. We define as a charge localization parameter the square of the largest orbital expansion coefficient in eq 1 (a different definition based on the root of the sum of all squared fractional charges96 yields very similar results; see Figure 6, bottom). For εMM = 1.0, the electron hole remains localized to more than 85% during most of the simulation trajectory. For εMM = 1.4, the localization parameters drop slightly, and for εMM = 2.0, typically only ca. 75% of the positive electron hole is localized on one site. Notably, in all three cases, the localization quickly drops during the initial stages of the simulations and increases once the end state of the charge on Trp306 is reached. Cause for Efficient Transfer across Trp359. To test the degree to which the electron hole would be capable of localizing on each residue and how much environmental stabilization each site provides, we have conducted six classical MD simulations using a fixed charge of +1 e on each of the six sites, respectively. To study how much fast environmental relaxation contributes to charge stabilization, the average ESP during the first 20 ps of simulation time was computed for each residue, depending on the location of the charge (see Figure 7). We find that even after such comparably short equilibration times, every side chain stabilizes charges well enough to form a stable localized state, with an energetic gap in IPs of more than 1 eV (except for the case of Trp316) compared to that of any other potential site. Nevertheless, we observed fast CT from Trp359 to Trp306 in the simulations above, which seems inconsistent with the existence of a stable charge trap state on

a

kF and kR indicate the forward and backward transfer rates, respectively. All rates are given ns−1.

involvement of Trp316 and is equivalent to the model presented in ref 31 (Supporting Information, Figure S14), but it is based on a different set of starting structures and a different QM zone selection during simulations. A very similar result, with kTrp382→Trp359 smaller than kTrp359→Trp306 and small but significant rates for backward transfers, is found. Extending the model to include transfer via Trp316 (second block in Table 2) changes kTrp382→Trp359 and kTrp359→Trp306 only slightly. Since kTrp382→Trp359 is ca. 4 times larger than k Trp382→Trp316 , and k Trp316→Trp382 is much larger than kTrp382→Trp316, the Trp382→Trp316 transfer only plays a small role, and Trp316 never accumulates a large fraction of charge. The uncertainties given for several rates in Table 2 indicate that different numerical fits did not converge to the same value. This involves cases where the elementary reaction in question does not influence the overall kinetics much, for example, transfer originating from the rarely populated Trp316. The charge occupation trajectories are too noisy to properly fit to these cases. No Trp316→Trp359 transfer is found at all (kTrp316→Trp359 is approximately zero). A comparison of these two kinetic models shows that extending the QM zone to include additional residues does not influence the overall kinetics significantly. Trp316, even though it occasionally has a nonzero charge occupation in some simulations, can be safely disregarded in the kinetic model. Electrostatic Potential Scaling. To explore the effect of electronic polarization on CT, we have conducted a series of 22 QM/MM simulations for scaling factors εMM of 1.0, 1.4, and 2.0, respectively. Here, the QM zone was restricted to Trp382, Trp359, and Trp306 to increase computational efficiency, since the above analysis showed that alternative charge carrying sites have little effect on the electron transfer process. The average occupation of the three QM sites is markedly affected by changes in εMM. Figure 5 and Table 2 show that rates increase markedly when εMM is changed from 1.0 to 1.4 and increase slightly further when εMM is set to 2.0. Additionally, fluctuations of the average occupation values increase with higher εMM, most notably for εMM = 2.0. For the forward transfer steps, much larger rate constants are found in both cases (compare third and fourth block of rates in Table 2 to the first block). εMM = 1.4 leads to approximately tripling kTrp382→Trp359 and doubling kTrp359→Trp306, while εMM = 10774



Article

Figure 8. Logarithmic residence time histogram of charge on Trp359 (top) and typical charge residence trajectories (bottom). When measuring the time for which the positive charge localizes on Trp359 (defined as >0.6 e on the site), residence times of less than 10 ps before CT back to Trp382 or forward to Trp306 is observed in the majority of cases. Histogram bars show the occurrence of residence times 100 ps. Bars are located at the average residence time of each bin. Plots of typical CT trajectories (bottom, colors as in Figure 6), indicating on which site (if any) the charge is localized during 1 ns of simulation times, reinforce this finding and show that localization of the positive charge on Trp359 is brief in most cases.

Figure 6. Effect of the polarization scaling factor εMM on average IP of QM sites (top) and electron hole localization (bottom). Scaling factors are 1.0, 1.4, and 2.0 from top to bottom. Increases in εMM reduce the energetic separation between sites and dampen energetic fluctuations. This results in significantly reduced hole localization properties. Especially at εMM = 2.0, a significant drop in localization is observed. The localization parameter, again averaged over all simulations, is defined as the fraction of the charge located on the highest occupied site (red) or as the square root of the sum of all fractional charges (black). Both definitions give very similar results.

histogram generation. Longer residence times of 100 ps or more are found regularly for the other two main charge carrying sites instead.

■

DISCUSSION Finding Electron Transfer Pathways. The calculations with extended QM zones, including six potentially charge carrying sites, still yielded an essentially unchanged picture of CT in PL: A positive electron hole moves from Trp382 to Trp359 and from there to Trp306, with forward rate constants in the 30−40 ns−1 range. Backward transfers play a small role, Trp316 does not accumulate significant amounts of positive charge, and the CT reaction is complete after less than 1 ns. The accumulation of a small amount of partial charge on Trp316 may be purely a simulation artifact, due to the tendency of Ehrenfest simulations to overly delocalize charge distributions. These finding are unsurprising, since the dominant role of the three Trp side chains on CT has been well established experimentally. Nevertheless, a structural analysis and the energetics shown in Figure 7 suggest at least the possibility of additional side chains involved in electron transfer. Extending the QM zone to include such additional amino acid residues and following the evolution of the moving charge distribution allowed us to demonstrate that the effect of the additional residues on CT is negligible. This shows our coarse-grained QM/MM approach to be a valuable predictive tool to establish transfer pathways in proteins where conclusively experimental validation is lacking, such as other members of the photolyase/ cryptochrome protein family. This findings make direct

Figure 7. ESP on each potentially charge carrying residue averaged from short 20 ps length MD simulations, depending on the location of the electron hole. As soon as the charge localizes on any site, fast environmental relaxation stabilizes it significantly, pushing the chargecarrying site to the lowest energy in each case. Neighboring sites experience additional, smaller stabilization. The sites close to the solvated protein surface (Trp306 and Tyr464) experience stronger stabilization.

Trp359. This can be explained by the fact that the positive charge never localizes on Trp359 long enough to trigger a strong environmental response. Figure 8 shows that the charge residence time on Trp359 is below 1 ps in most cases. The very large number of extremely short occupations in the 10 fs range is caused by rapid fluctuations of the noisy charge distribution trajectories. These very short “site occupations” occur clustered in rare parts of the trajectories, disappear when smoothed curves are used, and can be considered artifacts from the 10775



Article

Furthermore, the solvent significantly affects the kinetics of hole transfer. For the example, at the starting point of the CT reaction Trp382→Trp359, a positive electron hole is located on Trp382. A quantity of interest may then be the difference of IP for Trp359 and Trp382 under these conditions. In terms of Marcus theory, this difference is linked to λS and contributes an important component of the activation energy. It follows that the solvent treatment, especially the consideration of electronic polarization via scaling as discussed here, markedly affects the kinetics of the hole transfer, via changes in the outer-sphere reorganization energy. Again, these solvation effects increase when proceeding toward the protein surface. Now we can discuss the effect of the MM scaling scheme. The site energies of the three tryptophans (H0mm in eq 2) for different scaling factors are shown in Figure 6. Stronger scaling lowers the differences in the effective site energies, thus reducing both ΔG0 and λS at the same time. Therefore, the scaling of the MM electric field (a) makes the reaction less exergonic and (b) lowers the reaction barrier. These two competing effects combined lead to the accelerated rates shown in Table 2. The smaller ΔG0 in the scaled simulations also lowers the barrier for back-transfers. Therefore, the wave packet will experience increased rates for back-transfers, and, as a result, a stronger delocalization of the wave packet is found as shown in Figure 6. Efficient Transfer due to Incomplete Environmental Relaxation. Classical calculations show that each of the five selected tryptophan potentially charge-carrying sites could maintain a localized positive electron hole, indicated by the lower IP in comparison to all other sites. Therefore, each Trp side chain can in principle form a stable radical cation (on short time scales) from a thermodynamics perspective. However, in our coupled QM/MM calculations we only find long-term stable charge localization on the first and the final sites Trp382 and Trp306. This can be explained by the CT kinetics: Low electronic couplings prevent Trp384 from becoming populated at all, while for Trp316, a combination of a higher IP (as compared to Trp359) and a small electronic coupling to Trp382 prevent it from becoming an alternative intermediate site in the transfer pathway. This leaves the case of Trp359 on which the charge localizes, but only briefly, even though thermodynamics suggest it could serve as a long-term trap state. This discrepancy can be understood when taking the delayed environmental response into account. As soon as a positive charge localizes on a given side chain, it will polarize its environment, resulting in charge stabilization. Some of these stabilization processes like solvent reorientation occur fast compared to the transfer reaction time scale or extremely fast in the case of electronic polarization. However, in larger biomolecules, much slower relaxation time scales of up to several nanoseconds play important roles as well.68,70 Figure 7 shows that a large amount of charge stabilization is completed after 20 ps, but even longer simulations with charges fixed at a given site would allow additional slow protein conformational equilibrations to take place and increase the energy gap significantly (as shown for the case of three QM sites in Figure 12d in ref 31). This slower contribution to relaxation does not influence the dynamics for the case of Trp359, where charges typically localize for only a few picoseconds (Figure 8). The CT event from Trp359 to Trp306 therefore typically occurs out of nonequilibrium conformations and needs to overcome a significantly smaller

simulations of CT a more versatile tool than methods where charge donor and acceptor need to be defined explicitly. Electronic Polarization of the Environment. In our previous work, we have conducted calculations applying a fixed point charge model,31 and the estimated rate especially for the first step was underestimated by roughly 1 order of magnitude. One obvious reason for this disagreement could be the missing effect of electronic polarization, which was a main motivation for the calculations in the present study. As discussed above, nonpolarizable force fields tend to overestimate the solvent reorganization energy λS by roughly a factor 1.4−1. This motivated us to investigate effective ways to include polarization effects onto CT problems.69,72 It has been proposed that the effect of a bare charge in molecular mechanics context should be scaled down by roughly a factor of 1.4−1 to cover polarization implicitly.86 Here, this is done by scaling the QM/MM interaction, that is, the interaction of the excess charge with the MM environment. Including polarization this way scales down the interaction of the hole charge with the MM environment and leads, in a Marcus picture, to a lower solvent reorganization energy λS. This in turn speeds up the hole transfer rates by a factor of 3−6, depending on εMM, in our simulations. When comparing our computed rates to spectroscopic data, we find that experiments show the CT between Trp382 and Trp359 to occur within 33 ns−1.10,29,30,97,98 Comparing these experimental rates directly to the results in Table 2 shows that with a polarization scaling factor of εMM = 1.4 rates that are within a factor of 2 of experimental results are obtained for the forward CT rate constants. In our previous work where QM/MM simulations were conducted without polarization scaling, rates almost an order of magnitude below experimental numbers were found. Additionally, using εMM = 1.4 instead of εMM = 1.0 yields a higher rate for the initial Trp382→Trp359 CT than for the following Trp359→Trp306 one, in good agreement to experiment. Note that a direct comparison between theory and experiment is made more difficult since our kinetic modeling scheme required the inclusion of significant very fast rates for immediate backward CT reactions for a good fit, which are not accessible experimentally. Nevertheless, directly computing CT rates from unbiased QM/MM simulations in qualitative and even semiquantitative agreement to experimental data for a system as large and complex as PL confirms the power and wide applicability of our simulation approach. As discussed above, the energy profile for CT results mainly from interaction with solvent in the case of photoactivation in PL. The three tryptophan sites are chemically identical, and the protein environment alone does not modify the energetics drastically (indeed, when only the protein influence is considered in the QM/MM term, no CT away from Trp382 takes place at all). Trp306 is located on the surface of the protein, thus the solvent molecules rearrange (polarize) markedly when the hole is transferred to Trp306, providing a strong contribution to the stabilization of an electron hole on this site. This stabilizing effect of solvent becomes progressively weaker for charge-carrying sites located further away from the surface of the protein, for example, Trp359 and Trp382 in this case. Therefore, the polarizability of solvent drives the energetics of hole transfer in the DNA photolyase, making the hole transfer reaction exergonic overall. 10776



■

ACKNOWLEDGMENTS T.B.S. is indebted to the Baden-Württemberg Stiftung for the financial support of this research project by the Eliteprogramme for Postdocs.

energy gap than equilibrium simulations would suggest, resulting in much faster transfer rates. Presumably, the structure and dynamics of the photolyase enzyme have evolved to prevent Trp359 from becoming a charge trap state, thereby leading to a fast and efficient photoactivation process.

■

■

CONCLUSION We have applied a newly developed method to study the fast CT in DNA photolyase. This method does not presuppose mechanistic details, like charge localization, time scale separation, or adiabaticity of the transfer. In contrast, a fragment based Ansatz for the wave function is made and the whole system is propagated according to the derived coupled equations for hole wave function and classic nuclei. The methodology therefore follows the spirit of the electronic structure based nonadiabatic propagation schemes. We have shown that this is an important factor for the description of systems with fast CT like in PL. Here, the protein response to the moving charge is on the same time scale as the charge movement itself. Therefore, a straightforward application of the Marcus formula for CT fails, because standard ways for the calculation of the thermodynamic parameters ΔG0 and λS involve a complete relaxation of the environment. Interestingly, the protein electrostatic field itself does not contribute much to the driving force of this process. As shown recently, it is the solvent stabilization which leads to an exergonic process.31 That means that the solvent reorganization energy is larger for the more solvent exposed sites. Additional possible CT sites, not included in previous studies, do not contribute significantly to the CT pathway, since they are higher in energy and only barely populated during our simulations. Our previous work used a QM/MM methodology with fixed point charges: that is, electronic polarization effects, which are known to be imporant for CT, were not included. This led to CT rates about 1 order of magnitude too low. Electronic polarization has been shown to lower the solvent reorganization energy by roughly 30%. Therefore, in a fixed point charge QM/ MM simulation, the hole charge interacts too strongly with the solvent, leading to too large a solvent polarization energy. This matches the analysis of Leontyev and Stuchebrukhov,86 who argue that a scaling of net charges by 0.7 leads to an inclusion of electronic polarization in an effective way. It is similar to what is done here when scaling the ESP in the QM/MM calculations,69,72 which had not been applied to the study of CT in DNA PL before. Thereby, CT rates are corrected in the right direction, leading to a very good agreement with experimental rates.

■

REFERENCES

(1) Sancar, A. Chem. Rev. 2003, 103, 2203−2237. (2) Byrdin, M.; Sartor, V.; Eker, A. P. M.; Vos, M. H.; Auberta, C.; Brettel, K.; Mathisa, P. Biochim. Biophys. Acta, Bioenerg. 2004, 1655, 64−70. (3) Carell, T.; Burgdorf, L. T.; Kundu, L. M.; Cichon, M. Curr. Opin. Chem. Biol. 2001, 5, 491−498. (4) Lin, C.; Todo, T. Genome Biol. 2005, 6, 220. (5) Kelner, A. Proc. Natl. Acad. Sci. U.S.A. 1949, 35, 73−79. (6) Sancar, A. J. Biol. Chem. 2008, 283, 32153−32157. (7) Weber, S. Biochim. Biophys. Acta, Bioenerg. 2005, 1707, 1−23. (8) Moldt, J.; Pokorny, R.; Orth, C.; Linne, U.; Geisselbrecht, Y.; Marahiel, M. A.; Essen, L. O.; Batschauer, A. J. Biol. Chem. 2009, 284, 21670−21683. (9) Essen, L. O.; Klar, T. Cell. Mol. Life Sci. 2006, 63, 1266−1277. (10) Brettel, K.; Byrdin, M. Curr. Opin. Struct. Biol. 2010, 20, 693− 701. (11) Kodali, G.; Siddiqui, S. U.; Stanley, R. J. J. Am. Chem. Soc. 2009, 131, 4795−4807. (12) Beukers, R.; Eker, A. P. M.; Lohman, P. H. M. DNA Repair 2008, 7, 530−543. (13) Giovani, B.; Byrdin, M.; Ahmad, M.; Brettel, K. Nat. Struct. Biol. 2003, 10, 489−490. (14) Biskup, T.; Schleicher, E.; Okafuji, A.; Link, G.; Hitomi, K.; Getzoff, E.; Weber, S. Angew. Chem., Int. Ed. 2009, 48, 404−407. (15) Harrison, C. B.; ONeil, L. L.; Wiest, O. J. Phys. Chem. A 2005, 109, 7001−7012. (16) Chatgilialoglu, C.; Guerra, M.; Kaloudis, P.; Houee-Levin, C.; Marignier, J. L.; Swaminathan, V. N.; Carell, T. Chemistry 2007, 13, 8979−8984. (17) Tachikawa, H.; Kawabata, H. J. Phys. Chem. B 2008, 112, 7315− 7319. (18) Prytkova, T. R.; Beratan, D. N.; Skourtis, S. S. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 802−807. (19) Acocella, A.; Jones, G. A.; Zerbetto, F. J. Phys. Chem. B 2010, 114, 4101−4106. (20) Sadeghian, K.; Bocola, M.; Merz, T.; Schütz, M. J. Am. Chem. Soc. 2010, 132, 16285−16295. (21) Krapf, S.; Koslowski, T.; Steinbrecher, T. Phys. Chem. Chem. Phys. 2010, 12, 9516−9525. (22) Masson, F.; Laino, T.; Rö t hlisberger, U.; Hutter, J. ChemPhysChem 2009, 10, 400−410. (23) Koslowski, T.; Burggraf, F.; Krapf, S.; Steinbrecher, T.; Wittekindt, C. Biochim. Biophys. Acta, Bioenerg. 2012, 1817, 1955− 1957. (24) Zheng, X. H.; Garcia, J.; Stuchebrukhov, A. A. J. Phys. Chem. B 2008, 112, 8724−8729. (25) Zheng, X.; Ly, N. M.; Stuchebrukhov, A. A. Int. J. Quantum Chem. 2007, 107, 3126−3131. (26) Medvedev, D.; Stuchebrukhov, A. A. J. Theor. Biol. 2001, 210, 237−248. (27) Antony, J.; Medvedev, D. M.; Stuchebrukhov, A. A. J. Am. Chem. Soc. 2000, 122, 1057−1065. (28) Byrdin, M.; Eker, A. P. M.; Vos, M. H.; Brettel, K. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 8676−8681. (29) Aubert, C.; Vos, M. H.; Mathis, P.; Eker, A. P. M.; Brettel, K. Nature 2000, 405, 586−590. (30) Lukacs, A.; Eker, A. P. M.; Byrdin, M.; Brettel, K.; Vos, M. H. J. Am. Chem. Soc. 2008, 130, 14394−14395. (31) Woiczikowski, P. B.; Steinbrecher, T.; Kubar, T.; Elstner, M. J. Phys. Chem. B 2011, 115, 9846−9863. (32) Marcus, R. A. J. Chem. Phys. 1956, 24, 966−978. (33) Marcus, R. A. J. Chem. Phys. 1956, 24, 979−989.

ASSOCIATED CONTENT

S Supporting Information *

Details of the kinetic model for rate constant calculations and the procedure adopted to obtain numerical solutions for the rate equations. This material is available free of charge via the Internet at http://pubs.acs.org.

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 10777



Article

(34) Marcus, R. A. Annu. Rev. Phys. Chem. 1964, 15, 155−196. (35) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265− 322. (36) Levich, V. G.; Dogonadze, R. R. Dokl. Akad. Nauk SSSR 1959, 124, 123−126. (37) Hush, N. S. Trans. Faraday Soc. 1961, 57, 557−580. (38) Hopfield, J. J. Proc. Natl. Acad. Sci. U.S.A. 1974, 71, 3640−3644. (39) Jortner, J. J. Chem. Phys. 1976, 64, 4860−4867. (40) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 355, 796−802. (41) Warshel, A. J. Phys. Chem. 1982, 86, 2218−2224. (42) Hwang, J. K.; Warshel, A. J. Am. Chem. Soc. 1987, 109, 715−720. (43) Algen, R. G.; Parson, W. W.; Chu, C. T.; Warshel, A. J. Am. Chem. Soc. 1995, 117, 12284−12298. (44) Parson, W. W.; Chu, Z. T.; Warshel, A. Biophys. J. 1998, 74, 182−191. (45) Warshel, A.; Chu, C. T.; Parson, W. W. J. Photochem. Photobiol. A 1994, 82, 123−128. (46) Nonellam, M.; Schulten, K. J. Phys. Chem. 1991, 95, 2059−2067. (47) Schulten, K.; Tesch, M. Chem. Phys. 1991, 158, 421−446. (48) VandeVondele, J.; Sulpizi, M.; Sprik, M. Chimia 2007, 61, 155− 158. (49) Blumberger, J. Phys. Chem. Chem. Phys. 2008, 10, 5651−5667. (50) Kubař, T.; Elstner, M. J. Phys. Chem. B 2009, 113, 5653−5656. (51) Koslowski, T.; Steinbrecher, T. Z. Phys. Chem. 2009, 223, 739− 752. (52) Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Annu. Rev. Phys. Chem. 2010, 61, 461−485. (53) Beratan, D. N.; Skourtis, S. S.; Balabin, I. A.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Xiao, D. Acc. Chem. Res. 2009, 42, 1669−1678. (54) Balabin, I. A.; Beratan, D. N.; Skourtis, S. S. Phys. Rev. Lett. 2008, 101, 158102. (55) Grozema, F. C.; Siebbeles, L. D. A.; Berlin, Y. A.; Ratner, M. A. ChemPhysChem 2002, 3, 536−539. (56) Grozema, F. C.; Berlin, Y. A.; Siebbeles, L. D. A. J. Am. Chem. Soc. 2000, 122, 10903−10909. (57) Vekhter, B. G.; Ratner, M. A. J. Chem. Phys. 1994, 101, 9710− 9715. (58) Vekhter, B. G.; Ratner, M. A. Phys. Rev. B 1995, 51, 3469−3475. (59) Berlin, Y. A.; Siebbeles, L. D. A. Chem. Phys. Lett. 1998, 291, 85−93. (60) Siebbeles, L. D. A.; Berlin, Y. A. Chem. Phys. 1998, 238, 97−107. (61) Olson, M.; Mao, Y.; Windus, T.; Kemp, M.; Ratner, M. A.; Leon, N.; Mujica, V. J. Phys. Chem. B 1998, 102, 941−947. (62) Grozema, F. C.; vanDuijnen, P. T.; Berlin, Y. A.; Ratner, M. A.; Siebbeles, L. D. A. J. Phys. Chem. B 2002, 106, 7791−7795. (63) Hatcher, E.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Beratan, D. N. J. Am. Chem. Soc. 2008, 130, 11752−11761. (64) Beratan, D. N.; Skourtis, S. S.; Balabin, I. A.; Balaeff, A.; Keinan, S.; Venkatramani, R.; Xiao, D. Acc. Chem. Res. 2009, 42, 1669−1678. (65) Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Annu. Rev. Phys. Chem. 2010, 61, 461−485. (66) Wang, H.; Lin, S.; Allen, J. P.; Williams, J. C.; Blankert, S.; Laser, C.; Woodbury, N. W. Science 2007, 316, 747−750. (67) Skourtis, S. S.; Beratan, D. N. Science 2007, 316, 703−704. (68) Furse, K. E.; Corcelli, S. A. Phys. Chem. Lett. 2010, 1, 1813− 1820. (69) Kubar, T.; Elstner, M. Phys. Chem. Chem. Phys. 2013, 15, 5794− 5813. (70) Golosov, A. A.; Karplus, M. J. Phys. Chem. B 2007, 111, 1482− 1490. (71) Kubar, T.; Elstner, M. J. Phys. Chem. B 2010, 114, 11221− 11240. (72) Kubar, T.; Elstner, M. J. R. Soc., Interface 2013, 10, 20130415. (73) Doltsinis, N. L. Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms; John von Neumann Institute for Computing, Julich, Germany: 2002; Vol. 10; Chapter Nonadiabatic Dynamics: Mean-Field and Surface Hopping, pp 377−397. (74) Ehrenfest, P. Z. Phys. 1927, 45, 455−457.

(75) Meyer, H. D.; Miller, W. H. J. Chem. Phys. 1979, 70, 3214− 3223. (76) Micha, D. A. J. Chem. Phys. 1983, 78, 7138−7145. (77) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562−572. (78) Tully, J. C. J. Chem. Phys. 1990, 93, 1061−1071. (79) Kubař, T.; Woiczikowski, P. B.; Cuniberti, G.; Elstner, M. J. Phys. Chem. B 2008, 112, 7937−7947. (80) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260. (81) Kubař, T.; Elstner, M. J. Phys. Chem. B 2008, 112, 8788−8798. (82) King, G.; Warshel, A. J. Chem. Phys. 1990, 93, 8682−8692. (83) Ceccarelli, M.; Marchi, M. J. Phys. Chem. B 2003, 107, 5630− 5641. (84) Blumberger, J.; Lamoureux, G. Mol. Phys. 2008, 106, 1597− 1611. (85) Tipmanee, V.; Oberhofer, H.; Park, M.; Kim, K. S.; Blumberger, J. J. Am. Chem. Soc. 2010, 132, 17032−17040. (86) Leontyev, I.; Stuchebrukhov, A. Phys. Chem. Chem. Phys. 2011, 13, 2613−2626. (87) Park, H. W.; Kim, S. T.; Sancar, A.; Deisenhofer, J. Science 1995, 268, 1866−1872. (88) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (89) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (90) Meagher, K. L.; Redman, L. T.; Carlson, H. A. J. Comput. Chem. 2003, 24, 1016−1025. (91) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (92) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (93) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (94) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (95) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463−1472. (96) Pipek, J. Int. J. Quantum Chem. 1989, 36, 487−501. (97) Lukacs, A.; Eker, A. P. M.; Byrdin, M.; Villette, S.; Pan, J.; Brettel, K.; Vos, M. H. J. Phys. Chem. B 2006, 110, 15654−15658. (98) Byrdin, M.; Lukacs, A.; Thiagarajan, V.; Eker, A. P. M.; Brettel, K.; Vos, M. H. J Phys. Chem. A 2010, 114, 3207−3214. (99) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (100) Frishman, D.; Argos, P. Proteins: Struct., Funct., Genet. 1995, 23, 566−579. (101) Persistence of Vision Raytracer; Persistence of Vision Pty. Ltd.: Williamstown, Victoria, Australia, 2004.

10778


Charge Transfer in E. coli DNA Photolyase: Understanding

Recommend Documents