6894
J. Phys. Chem. B 2000, 104, 6894-6902
Calculation of Quantum Parameters for Nonadiabatic Redox Reactions. Application to Photoreduction of Flavin in DNA Photolyase EunJoo Lee, Emile S. Medvedev,† and Alexei A. Stuchebrukhov* Department of Chemistry, UniVersity of California, DaVis, California 95616, and Institute of Chemical Physics Problems, Russian Academy of Sciences, 142432 ChernogoloVka, Moscow, Russia ReceiVed: January 12, 2000; In Final Form: May 10, 2000
A simple practical way to account for the effect of quantum modes on electron transfer (ET) is to use the Jortner expression for ET rate. The expression includes two quantum parameters which describe the properties of high-frequency modes in the system. In our recent paper (J. Chem. Phys. 2000, 112, 9015), we developed a method to calculate these parameters for redox cofactors from ab initio data on their potential energy surfaces. In this paper, we extend our method to include the solvent, which is treated as a continuous dielectric medium. As an example, two reactions describing photoreduction of the flavin cofactor in DNA photolyase (an enzyme that repairs thymine dimers in DNA) are investigated. We calculated the quantum parameters for each cofactor (flavin radical, tryptophan, tryptophan cation radical, and tyrosine), as well as for water, which was used as a model medium, and then obtained the dependence of ET rates kET on the driving force ∆G0 in a wide range of ∆G0. As expected, the quantum modes tend to flatten the log kET vs ∆G0 dependence in the inverted region, and the calculation provides a quantitative estimate of this effect. A similar effect has been predicted to originate from inelastic tunneling (Medvedev, E. S.; Stuchebrukhov, A. A. J. Chem. Phys. 1997, 107, 3821). This effect makes kET relatively insensitive to ∆G0 variations, compared with the classical Marcus theory, and might be of importance for biological ET systems to ensure their stable performance.
1. Introduction A useful expression for nonadiabatic electron transfer (ET) rate which takes into account contributions of both classical and quantum modes was given by Jortner1
Sqn kET ) e-Sq kclET(∆G0 + npωq) n! n)0 ∞
∑
(1.1)
where
kclET(∆G0)
|TDA|2 ) p
x
[
]
(λ + ∆G0)2 π exp λkBT 4λkBT
(1.2)
is the classical Marcus rate. This expression depends on five parameters: the driving force of the reaction, ∆G0; the electrontransfer matrix element, TDA; the reorganization energy of the classical modes, λ ) λcl; the average frequency of the quantum modes, ωq; and Sq which represents a dimensionless coupling between the tunneling electron and the quantum modes. Equation 1.1 assumes that the quantum effect can be described by a single “effective” quantum mode which is characterized by only two parameters, ωq and Sq. These quantum mode parameters are usually considered as adjustable in fitting experimental data.2 Kumar et al.3 have recently analyzed various methods to estimate all five parameters, combining theoretical calculations with experimental data. There also exists an exact quantum expression for the nonadiabatic transition rate, which takes into account all quantum modes of the system with their shifts of equilibrium positions, frequency changes, and Duchinsky rotation.4 The †
Russian Academy of Sciences.
advantage of the exact expression is that all relevant parameters can be calculated for a given system using available methods of quantum chemistry. However, it is rarely used in practice because of its complexity. Recently, we have presented5 a method to calculate the quantum parameters ωq and Sq from harmonic potential energy surfaces (PES) of the cofactors participating in a reaction. We have derived simple equations that relate ωq and Sq to the full set of parameters characterizing the PES of the cofactors, which include shifts of equilibrium positions, frequency changes, and Duchinsky rotation matrix. The resulting Jortner approximation for the rate, eq 1.1, with our definition of ωq and Sq was shown to be in excellent agreement with the exact rate over a wide range of ∆G0 where the rate varied by several orders of magnitude. In the traditional approach3 an estimate of ωq is used, assuming that the frequency dispersion of the quantum modes is small. In contrast, we show that for any dispersion there is an optimal way of choosing ωq which results in the best approximation of the exact rate. In the reactions under study in this paper, the internal reorganization energy is distributed in the entire frequency range, including both low ( ωcl) regions, so that R’s are represented as
R1 ) R1cl + R1q = Rcl + R1q
(2.10)
R2 ) R2cl + R2q = Rcl + R2q
(2.11)
where the classical parts of R’s are approximately the same and Rcl is defined as arithmetic average of R1cl and R2cl. Being multiplied with Λ, they give the medium reorganization parameters
λscl ) ΛRcl, λs1q ) ΛR1q, λs2q ) ΛR2q
(2.12)
The medium quantum parameters ωsq and Ssq are defined in terms of λs1q and λs2q by the same eqs 2.5 and 2.6 as those for donor and acceptor. From these equations it is seen that ωsq can be defined solely in terms of the medium R’s, whereas Ssq requires knowledge of Λ too. It is essential to emphasize that λcl, λ1q, and λ2q are additive in the sense that their partial values found for each cofactor and medium individually sum up for the whole reaction. This is because it is assumed that in a long-distance biological electron transfer the donor and acceptor complexes are well separated from each other. Then, the quantum parameters for the whole reaction are expressed in terms of those for donor, acceptor, and medium by these simple relations5
ωq )
Sq )
(ωdq)2 Sdq + (ωaq)2 Saq + (ωsq)2Ssq ωdq Sdq + ωaq Saq + ωsq Ssq (ωdq Sdq + ωaq Saq + ωsq Ssq)2 (ωdq)2Sdq + (ωaq)2Saq + (ωsq)2Ssq
(2.13)
photoinduced electron transfer from FADH- to a thymine dimer, which splits upon the one-electron reduction14-17 (for theoretical work see ref 18 and references therein). Under certain conditions, FADH- can lose one electron, to become FADH•, and the enzyme becomes catalytically inactive. The catalytically inert cofactor can regain its activity, however, if the enzyme is exposed to near-UV-visible light. The reaction responsible for the flavin reactivation is believed to be an electron transfer from a tryptophan residue of the protein to the photoexcited oxidized flavin.19-25 The oxidized tryptophan is then reduced by an electron from a tyrosine residue, presumably to prevent a back reaction between the reduced flavin and oxidized tryptophan. The oxidized tyrosine, in turn, later is reduced by an external agent from the solvent. Such a chain of ET reactions was observed in A. nidulans photolyase and was proposed to be the mechanism of FADH reduction in the enzymes from other organisms.26 The unusual redox reactions in which actual ionization of the protein side chains are involved are gaining much interest in the current literature.27,28 The two ET reactions that lead to reactivation of the flavin cofactor and involve ionized states of tryptophan and tyrosine are the subject of the discussion in the following sections of this paper. The reactions of tryptophan and tyrosine in fact involve both electon and proton transfer since the neutral radicals of these side chains are observed in the experiments. Although little is known at present about the proton transfer in these systems,26 it is believed that the primary reaction is electron transfer, which is followed by the proton transfer as the system relaxes in the product electronic state. According to the literature data,23,29,30 the absorption of a photon by the ground-state FADH• yields the excited doublet state of the radical, which within 100 ps undergoes intersystem crossing to the quartet state 4FADH•. It is this long-lived 4FADH• that abstracts an electron from the tryptophan residue and eventually forms catalytically active two-electron reduced singlet FADH-. The final state of the flavin reduction reaction, however, is the triplet 3FADH-. Therefore, the following reactions will be considered:
FADH• + TrpH f 3FADH- + TrpH•+
4
(2.14)
where superscripts d, a, and s mark the donor, acceptor, and medium parameters, respectively. Hence, the quantum parameters of each cofactor and medium need to be calculated only once, and then they can be used in calculations of the ET rate as a function of driving force for any other reaction where this cofactor participates. It is worthwhile to emphasize that the direct methods of calculating the exact ET rate (i.e., numerical integration and saddle-point calculation) require full-length calculations for each value of the driving force (see, e.g., ref 7). The rest of the paper describes an example of implementation of the above scheme to two particular reactions between redox cofactors of DNA photolyase. 3. Redox Cofactors in DNA Photolyase Photolyase is an enzyme that catalyses the photorepair of thymine dimers in UV-damaged DNA by an electron transfer reaction. The enzyme’s catalytic cofactor is doubly reduced flavin adenine dinucleotide, FADH-. The repair reaction is a
TrpH•+ + Tyr f TrpH + Tyr•+ Details of these reactions and the models of the reaction cofactors are shown in Figure 1. As a model for the flavin cofactor, acceptor in our ET reaction, we take lumiflavin, 7,8-dimethylisoalloxazine with a methyl group at nitrogen 10. Tryptophan is modeled by 3-methylindole, and tyrosine by p-cresol. Our first scheme in Figure 1 illustrates the first electron transfer mechanism between flavin radical and tryptophan. Scheme 2 in Figure 1 shows the second electron-transfer reaction between tryptophan cation radical and tyrosine. 4. Results and Discussion 4.1. Geometries and Spin Densities. The spin-restricted Hartree-Fock method and the hybrid density functional theory (Beck3-Lee-Yang-Parr functional) were used for geometry optimization and frequency calculations in closed-shell systems (tryptophan and tyrosine cofactors). The open-shell systems (flavin radical, tryptophan cation radical, and tyrosine cation radical) were treated with the spin-unrestricted open-shell
Quantum Parameters for Nonadiabatic Redox Reactions
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6897
Figure 1. Model for ET reactions in DNA photolyase: (scheme 1) the ET reaction between flavin radical and tryptophan; (scheme 2) the ET reaction between tryptophan cation radical and tyrosine.
Hartree-Fock (UHF) method, spin-restricted open-shell HF (ROHF) method, and spin-unrestricted open-shell B3LYP functional, in order to compare the values of reorganization energy and quantum parameters obtained by different methods. Single-point energy calculations at the HF/6-31G* geometries were performed at the second-order Mo¨ller-Plesset theory (MP2). All calculations were performed employing the 6-31G* basis set and GAUSSIAN 94 code.31 To obtain redox quantum parameters by our method, the two electronic states are analyzed at the same geometry. In this calculation the equilibrium geometry of the initial state of the redox cofactor (see Figure 1) is used for such an analysis. 4.1.1. FADH. Frequency calculation for flavin neutral radical 4FADH• by UHF and ROHF methods gave two imaginary frequencies at the planar optimized geometry, suggesting the structure to be a second-order saddle point rather than a true minimum. In order to obtain the true minimum structure we distorted the molecule along the normal mode which corresponds to the highest imaginary frequency and rerun the optimization followed by frequency calculation. The modified structure yielded no imaginary frequencies, and its total energy was lowered by 0.04 eV for UHF method and 0.1 eV for ROHF method. This new minimum strucutre is slightly bent along the N5-N10 axis by about 10°. However, no imaginary frequencies was obtained for the planar structure of 4FADH• using the DFT method. Therefore, depending on the method, we obtained two different geometrical structures of the flavin radical. Geometry optimization of the two-electron-reduced flavin (3FADH-) using HF method also yielded a nonplanar structure whereas the DFT method produced a planar structure. This nonplanar structure of flavin was also obtained by other authors. In their NMR studies of reduced flavin in solution, Moonen et al.32 suggested that the reduced flavin is slightly bent in nonpolar solvents and is planar in polar solvents. Zheng and Ornstein33 observed that energy of the puckered flavin semiquinone radical is lower than that of the planar structure by 0.05 eV at the Hartree-Fock level, which is in good agreement with our estimated energy
difference (0.04 eV) between the nonplanar and planar structures obtained by the Hartree-Fock method. We calculated by different methods the total electronic energy of neutral quartet flavin radical (4FADH•). Since there are no published data of ab initio or semiempirical calculations on the quartet flavin radical, we had to use the doublet flavin radical as a reference system for comparison. Zheng and Ornstein33 theoretically studied the flavin structure in different oxidation and protonation states. As a model compound, lumiflavin was used by these authors. The redox active site of this cofactor is the isoalloxaine system, which can evolve in three reversible oxidation steps, i.e., via quinoid (oxidized), semiquinoid (radical), and hydroquinoid (reduced) states. These flavin radicals can exist in different protonation states depending on the pH of the solution, viz., in cation, neutral, and anion forms. On the basis of Zheng and Ornstein’s results, we selected the flavin semiquinone radical (compound No.5 in ref 33), which has two protons at N3 and N5, since its structure most closely resembles our neutral one-electron reduced flavin 4FADH•. The energy obtained by us excellently matches that from ref 33. The agreement between our total electronic energy and geometrical parameters (discussed later in this section) of the doublet flavin radical and those by Zheng and Ornstein33 and Meyer et al.34 served as a solid starting point for our ab initio calculation of the quartet state of the neutral flavin radical (4FADH•). The single point MP2 corrections to the HF energies (both UHF and ROHF) resulted in a significant decrease of the energy of quartet flavin radical. We now examine the geometrical parameters of the quartet flavin radical obtained by different methods and compare them with the doublet state geometries obtained in refs 33 and 34. Table 1 lists the nonplanar geometrical parameters for 4FADH• calculated by the HF method and the planar structure parameters by DFT method. The UHF geometrical parameters were compared with Zheng and Ornstein’s compound No. 5 (2FADH•) structure, and the ROHF parameters were compared with Meyer et al.’s model (10-methylisoalloxazine) geometry. From Table
6898 J. Phys. Chem. B, Vol. 104, No. 29, 2000
Lee et al.
TABLE 1: Structural Parameters of the One-Electron-Reduced Neutral Flavin Radicala parameter
TABLE 3: Spin Density of the Tryptophan Cation Radical (Reaction II)
QT-NPSUHF DBUHFb QT-NPSROHF DBROHFc QTB3LYP
N1-C2 C2-N3 N3-C4 C4-C4a C4a-C10a C10a-N1 C4a-N5 N5-C5a C5a-C9a C9a-N10 N10-C10a C5a-C6 C6-C7 C7-C8 C8-C9 C9-C9a C2-O2 C4-O4
1.381 1.393 1.363 1.456 1.437 1.295 1.361 1.403 1.505 1.395 1.365 1.380 1.392 1.529 1.388 1.394 1.191 1.201
1.374 1.401 1.366 1.447 1.432 1.290 1.344 1.374 1.406 1.411 1.365 1.398 1.391 1.411 1.400 1.397 1.192 1.205
1.376 1.400 1.360 1.454 1.440 1.287 1.344 1.389 1.456 1.401 1.361 1.331 1.489 1.496 1.331 1.459 1.191 1.202
1.373 1.402 1.359 1.451 1.433 1.289 1.343 1.387 1.399 1.407 1.362 1.384 1.382 1.383 1.385 1.399 1.191 1.202
1.375 1.417 1.382 1.448 1.430 1.338 1.358 1.387 1.490 1.372 1.395 1.373 1.412 1.484 1.367 1.432 1.227 1.234
a Distance in Å. The 6-31G* basis set was used for all methods. DB and QT stand for the doublet and quartet states, respectively. NPS stands for the nonplanar structure. b Zheng and Ornstein.33 c Meyer et al.34
TABLE 2: Spin Density of the Neutral 10-Methylisoalloxazine (4FADH•) Radical atom
UHF
N1 0.4 C2 -0.14 O2 0.18 N3 -0.006 C4 -0.17 O4 0.18 C4a 0.82 N5 0.14
ROHF UB3LYP atom 0.08 0.01 0.02 0.01 0.04 0.04 0.5 0.2
0.35 -0.04 0.21 0.02 -0.02 0.11 0.44 0.22
C5a C6 C7 C8 C9 C9a N10 C10a
UHF 0.83 -0.75 0.92 0.89 -0.76 0.81 -0.02 -0.28
ROHF UB3LYP 0.04 0.15 0.9 0.09 0.03 0.7 0.08 0.06
0.31 -0.12 0.58 0.23 -0.05 0.46 0.15 0.13
1 it is seen that the quartet and doublet state geometries differ insignificantly except for a few bonds. In particular, an elongation of C7-C8 bond in the quartet state by approximately 8% as compared to the doublet state was obtained by all methods. The UHF method gave a 7% elongation of the C5aC9a bond with respect to the doublet state, and a 7.7% C6-C7 bond elongation was obtained with the ROHF method. The bond elongations in the quartet state can be explained by two extra unpaired spins. This effect can be elucidated further by studying spin density distributions in both states obtained with different methods. Spin density of flavin neutral radical (4FADH•) is presented in Table 2. The DFT and ROHF methods provide a high spin density on C4a, C7, and C9a. The UHF method bears a high spin density on C4a, C5a, C7, C8, and C9a, but also has a large negative spin density in positions at C6 and C9, which will introduce some spin contamination from higher multiplets. Despite these differencies, all three methods showed that a large spin density was on the quinoxaline ring, which has been also observed by the ENDOR study of a flavosemiquinone model system.35,36 Recent theoretical studies of the doublet neutral flavin radical33,34,37 showed that the highest spin density is at the C4a site, which was also observed by us in the quartet state. In the quartet state we observe two more sites of high spin density, mainly on C7 and C9a, which contribute to the bond elongation. The high spin density at the C7 site results in a C7-C8 bond elongation using all methods, and the high density on the C5a site obtained with the UHF method contributes to the C5a-C9a bond elongation. 4.1.2. TrpH and Tyr. We now turn to the tryptophan cation radical, which serves as an electron acceptor in the second
atom
ROHF
B3LYP
ROHFa
B3LYPb
exptc
N1 C2 C3 C4 C5 C6 C7 C8 C9 C10
0.100 0.266 0.415 0.062 0.002 0.054 0.026 0.042 0.008 0.004
0.138 0.186 0.36 0.259 -0.088 0.178 0.07 0.021 -0.113 -0.02
0.1 0.24 0.4 0.09 0.00 0.08 0.04 0.04 0.01
0.136 0.189 0.36 0.259 -0.089 0.179 0.07 0.022 -0.113 -0.020
0.14 0.35 0.41
a
Krauss et al.45
b
Jensen et al.46
c
-0.07
Huyett et al.42
electron transfer reaction (see Figure 1). It has attracted much attention in the literature because it exists in many important biological systems.19,20,23,38-44 We calculated the geometry and the Mulliken spin densities of tryptophan cationic radical using the ROHF method and the B3LYP functional with the 6-31G* basis set. Table 3 lists three sets of spin densities from papers42,45,46 for the purpose of comparison with our data. Our B3LYP and ROHF results are in good agreement with Jensen’s spin densities on 3-methylindole (B3LYP),46 and with Krauss et al.’s results (ROHF).45 All methods have shown a significant spin density on the C3 site, which was also observed in the ENDOR experiments on measuring spin densities for tryptophan-191 of cytochrome-c-peroxidase.42 Our results on the geometrical parameters and the total energy for the tryptophan cation radical agree with those by Jensen et al.46 The tyrosine cofactor in its initial state is in a neutral form, and its geometry analysis was quite standard. We have used the DFT method using the B3LYP functional with the 6-31G* basis set for the analysis. We compared our calculated geometry parameters with Qin and Wheeler’s ab initio calculation47 on tyrosine phenoxyl radical (TyrO•) using the SVWN functional with the 6-31G* basis set. Except for the C-O bond, the differences between two calculated bond distances are under 4%. Since we have an O-H hydrogen bonding, our C-O bond was larger by 9%. 4.2. Harmonic PES Parameters of Redox Cofactors. Once the equilibrium configuration of the initial state of a cofactor is established, we then determine the normal-mode frequencies ωk and eigenvectors in this state. One of our basic assumptions is that the vibrational modes can be divided into the classical (ωk < ωcl ≡ kBT/p) and quantum (ωk > ωcl) modes, and that the classical and quantum modes play fundamentally different roles in the ET dynamics. The quantum modes contribute to forming two sets of quantum levels between which the transitions occur, whereas the classical modes cause Gaussian fluctuations of these levels; due to these fluctuations the resonance condition can be fulfilled for electron transfer to occur. As a consequence of this difference, the classical modes are described collectively by a single parameter, the root mean square of the energy fluctuation, whereas the quantum modes are characterized by their individual parameters. Therefore, in the first step of calculation, we eliminate the classical modes (along with the translational and rotational modes) by projecting the complete Hesssian of the initial state onto the subspace of the quantum modes. Our second assumption is that the quantum modes can be described in the harmonic approximation. Usually, in the systems under study, the energy dumped into vibrations is small, on the order of a few typical vibrational quanta. Since this
Quantum Parameters for Nonadiabatic Redox Reactions TABLE 4: Calculated Half-Reaction λ’s, eVa cofactor
method
λ1q
λ2q
ΛT
4FADH•
UHF UHF/MP2 ROHF ROHF/MP2 UB3LYP UHF UHF/MP2 ROHF ROHF/MP2 UB3LYP UHF UHF/MP2 ROHF ROHF/MP2 UB3LYP UHF UHF/MP2 ROHF ROHF/MP2 UB3LYP
0.43
1.35
0.38 0.25 0.55 0.16 0.12 0.41 0.23 0.31 0.19 0.18 0.41 0.17 0.32 0.33 0.15 0.35 0.32 0.34 0.32 0.21
TrpH
TrpH•+
Tyr
a
0.45
1.81
0.13 0.30
0.38 1.13
0.36
1.29
0.21 0.49
0.57 1.44
0.33
1.13
0.17 0.34
0.60 1.25
0.35
1.17
0.21
0.64
For all methods the 6-31G* basis set was used.
energy is shared among several modes, the average excitation at a single mode does not exceed one vibrational quantum. Therefore, the harmonic approximation for the quantum modes should work well. This latter consideration should be taken into account when performing the next-step calculation related to the final state of the given cofactor. The frequency job for the final state is run at the same configuration as for the initial state, i.e., at the equilibrium configuration of the initial state.48 As a result, we obtain not only a Hessian matrix but also a nonvanishing force vector since the configuration is not the equilibrium of the final state. Since we already have got the eigenvectors for the classical modes, as well as for the translational and rotational modes, we can project the Hessian matrix and the force vector onto the subspace of the quantum modes, and then we diagonalize the projected Hessian to obtain the normal modes in the final state. It may occur, however, that a few modes have large displacements of equilibrium positions between the initial and final electronic states due to strong anharmonicity. (For instance, nearly freely rotating CH3 groups may turn over large angles.) According to our basic assumption these modes should be treated as classical, and therefore they have to be projected out from the Hessian before diagonalization. The above procedure have been described in more detail in ref 5. After it is done, we obtain shifts of equilibrium positions, frequency changes, and the Duchinsky rotation matrix for the quantum modes, and the next step is to calculate the quantum parameters governing the ET rate. 4.3. Quantum Parameters for Redox Cofactors. The values of λ1q and λ2q calculated for individual cofactors by different ab initio methods are given in columns 3 and 4 of Table 4. Column 5 contains the total (classical plus quantum) reorganization energy λT obtained by single-point calculation as the energy difference between the initial and final states, both at the equilibrium configuration of the final state. In principle, if all modes were exactly harmonic, we would obtain λT ) λcl + λ1q. Therefore, λT approximately represents an upper boundary on the values of λ1q and will be used below for comparison with λ1q. We observe from Table 4 that the HF methods give significantly higher reorganization energies than DFT. It is known that the neglect of electron correlations tends to overestimating the vibrational frequencies by about 10%, and
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6899 TABLE 5: Quantum Parameters Calculated by Eqs 2.5 and 2.6 Using the DFT Data of Table 4 half-reaction 4FADH•
f 3FADHTrpH f TrpH•+ TrpH•+ f TrpH Tyr f Tyr•+
ωq, cm-1
Sq
1220 1130 1470 1270
0.9 1.5 0.9 1.3
this is indeed observed in our calculations. However, this effect is too small to explain the large differences in λ1q. Therefore, the main source of discrepancy is in the overestimated forces in the final state. On average, the DFT method gives lower forces than HF, thus resulting in lower λ1q’s. This effect is more pronounced in the open-shell cofactors (the radicals 4FADH• and TrpH•+) than in the closed-shell ones (TrpH and Tyr) since the former are more sensitive to the electron correlations. If we compare the λT values obtained by different methods, we will see the same trend. That is, the HF methods overestimate λT, as compared to DFT, by a factor of 3 for the open-shell cofactors, and by a factor of 1.5 for the closed-shell ones. In order to elucidate this point further, we calculated λT using the second-order Mo¨ller-Plesset (MP2) perturbation theory. The results are also shown in Table 4. It is seen that inclusion of the electron correlations significantly diminishes λT, making it closer to the DFT result. This is in line with the results of Kumar et al.3 (see their Table 1). For the closed-shell cofactors the effect of electron correlations is smaller, as expected. Thus, we conclude that the electron correlations significantly contribute to the values of quantum parameters. An additional indication of the failure of the HF methods for the open-shell cofactors comes from a comparison of λT and λ1q values. For 4FADH• and TrpH•+, λ 1q exceeds λT by 0.05 and 0.08 eV, respectively, which contradicts to the definition of these quantities (λ1q e λT). The DFT parameters also feature this inconsistency, but to a lesser extent, i.e., 0.03 eV at most. Therefore, the DFT method seems to be more appropriate for such calculations. The DFT values are used to calculate the quantum parameters ωq and Sq which are shown in Table 5. As was pointed out at the beginning of section 4.2, the full vibrational space was reduced by removing the classical modes, and additionally some modes featuring large displacements between the initial and final states. In practice, there were only a small number of modes that needed to be removed (11 out of 90 in the case of flavin). Therefore, the fact that the λT and λ1q DFT values are very close to each other just implies that the reorganization energy of the inner-sphere classical modes is too small to be separated out of λT by subtracting λ1q. Therefore, we will neglect the inner-sphere contribution to λcl, which is in line with common practice to derive λcl from the outer-sphere degrees of freedom. The above argument implies that a complicated calculation of λ1q could be replaced with a simpler calculation of λT. Still, there is the second parameter, λ2q, for which no such a simple prescription can be given. This parameter governs ωq according to eq 2.5. The values of ωq shown in Table 5 fall into the interval from 1100 to 1500 cm-1 for the cofactors under study, which is similar to estimates of Kumar et al.3 for other class of systems (1400 ( 200 cm-1). The variation of ωq by a factor of 1.4 is significant because λ2q (hence, ωq) governs the width and asymmetry of the log kET vs ∆G0 curve. Therefore, an accurate calculation of λ2q is necessary in order to obtain a valuable estimate for the rate. The accuracy of the rate calculation mostly depends on the accuracy of the ab initio method employed. The above-mentioned discrepancy of 0.01-0.03 eV between λT and λ1q can be both due to anharmonicity and inaccuracy of the
6900 J. Phys. Chem. B, Vol. 104, No. 29, 2000
Lee et al.
TABLE 6: Reorganization Parameters (r’s) and the Quantum Mode Parameters (ωsq and Ssq) for Water ωcla
Rclb
R1q
R2q
R1c
R2c
ωsq, cm-1
Ssq d
1 1.5 2
0.33 0.38 0.41
0.25 0.20 0.17
0.54 0.49 0.45
0.58
0.87
900 1020 1100
2.2 1.6 1.2
a In units of kBT/p ) 3.9 × 1013 rad/s ) 207 cm-1 (T ) 300 K). The arithmetic average of R1cl and R2cl. c R1 and R2 are independent of ωcl. d Calculated for Λ ) 1 eV.
b
Figure 2. Calculated rate for a hypothetic reaction between spherical donor and acceptor with no vibrational degrees of freedom in water. Thin solid, dotted, and dashed lines (nearly indistinguishable from each other): eq 1.1 with parameters from Table 6 for three possible values of ωcl. Thick solid line: the exact rate by saddle-point calculation. Thick dashed line: eq 1.2 with λ ) λscl. The ordinate’s zero is arbitrary. Λ ) 1 eV.
DFT method. If we consider this mismatch as an estimate for the accuracy of the DFT method, then we obtain that our calculated values of λ1q and λ2q bear uncertainties within 1015%, and hence, ωq and Sq are uncertain within 20-30%. 4.4. Quantum Parameters for the Medium. Example of Water. The calculation of the quantum parameters for the medium by eqs 2.10-2.12 deserve special comments since the separation of the medium reorganization parameters into their classical and quantum components contains an uncertainty depending on the choice of ωcl. In this connection we emphasize that no such a problem arose with the inner-sphere parameters (see at the end of section 4.3) since there the contribution of the classical modes is small and not discernible with the ab initio methods used. This is not the case, however, with the outer-sphere parameters. In Figure 5 of paper7 the integrand in eq 2.7 as a function of log ω for water is shown. From that figure it is seen that ωcl is located in a region where the integrand is essentially nonvanishing. Hence, a change in ωcl will affect the parameters in the right-hand sides of eqs 2.10-2.12, as well as the quantum parameters ωq and Sq. Explicitly, this is shown in Table 6, where the calculated values of R’s and ωq are given for three values of ωcl. The values for Sq are also shown for Λ ) 1 eV. Only the total integrals, R1 and R2, are independent of this choice. However, both the reorganization parameters Rcl, R1q, and R2q and the effective quantum mode parameters ωsq and Ssq are affected by the choice of ωcl. Therefore, we have to find out how accurately our approximation represents the exact rate and how ωcl influences the calculated rate. In Figure 2 we show the calculated ET rate for a hypothetic reaction between structureless (with no vibrational degrees of freedom) donor and acceptor in water. The thick solid line is exact calculation using the first saddle-point approximation.6,49
TABLE 7: Parameters for the Whole Reactions Calculated by Eqs 2.13 and 2.14 Using the Data of Tables 5 and 6 (Λ ) 1 eV) full reaction
λcl, eV
ωq, cm-1
Sq
+ TrpH f 3FADH- + TrpH•+ TrpH•+ + Tyr f TrpH + Tyr•+
0.33 0.33
1050 1180
4.5 4.2
4FADH•
The thick dashed line is calculated by eq 1.2 with λ ) λscl. This is to show that the quantum effects are large in the inverted region. Three thin lines are calculated by eq 1.1 for three values of ωcl given in Table 6. The abscissa is intentionally expanded deep into the inverted region in order to show that the divergence of the thin curves is small. It is seen, first, that the Jortner formula represents the exact rate rather accurately in a wide region of ∆G0 where the rate varies by 3 orders of magnitude. Therefore, we can state that the Jortner formula captures most of the quantum effect correctly provided that the quantum parameters are chosen according to our prescription. Thus, we have demonstrated that the quantum effect of a dielectric medium can be described by the same eq 1.1 as that of the inner-sphere vibrational modes. Second, the variation of the Jortner curve with ωcl is insignificant (thin lines are nearly indistinguishable from each other in the figure) even though the quantum parameters themselves, ωsq and Ssq, vary essentially (see Table 6). We remind that the primary parameters, λs1 and λs2, were kept constant in these calculations, and only their dividing into the classical and quantum parts was varied by different choice of ωcl. Therefore, we conclude that a reasonable variation of ωcl around kBT/p does not affect the rate calculated by eq 1.1. Now we are in position to calculate the quantum parameters for the whole reactions between the selected cofactors in a given medium. Yet, these are hypothetical reactions since we do not know (ω) for the biological medium which includes the protein itself, the bound water, and the free water where the protein is dissolved. The dielectric measurements for such media are difficult to do, and only limited data are available on powedered proteins (lysozyme) with a few layers of bound water.50,51 These measurements are restricted by the low-frequency region (ω < 10-3ωcl), so that the quantum effects of the medium cannot be considered. Therefore, we select water as a medium, and for this same reason we do not calculate precise values of Λ for the particular reactions under study. Instead, we use the Marcus formula to obtain a rough estimate. Thus, for instance, taking ion radii 7 and 10 Å and their separation 15 Å we obtain Λ ) 0.8 eV. Shown in Table 7 are the quantum parameters for the whole reactions calculated for Λ ) 1 eV. λcl is equal to λscl since donor and acceptor do not contribute, and λscl is calculated by eq 2.12. 4.5. Rate Calculations. The dependence of log kET vs ∆G0 calculated by eqs 1.1 and 1.2 with parameters from Table 7 for the reactions of flavin with tryptophan and tryptophan with tyrosine are shown in Figure 3, where the Marcus curve is also shown for the 4FADH• + TrpH reaction (the Marcus curve for the TrpH•+ + Tyr is nearly the same and is not shown). It is seen from Figure 3 that the quantum effect for the two reactions is especially large in the inverted region, as expected. The rates of the two reactions are very similar since the cofactors have very close values of total λ’s, which results in similar ωq and Sq (see Table 7). It is interesting to compare the quantum effects due to innersphere and outer-sphere modes. We remind that each λ in eqs 2.1 and 2.2 is a sum of inner-sphere and outer-sphere terms. Then, it is easy to manipulate them in order to obtain the desired
Quantum Parameters for Nonadiabatic Redox Reactions
Figure 3. Rates of the reactions 4FADH• + TrpH (thick solid line) and TrpH•+ + Tyr (thin solid line) in water calculated by eq 1.1 with parameters from Table 7. The Marcus rate for the 4FADH• + TrpH reaction is also depicted (thick dashed line). The Marcus rate for the TrpH•+ + Tyr reaction is nearly the same (not shown).
parameters. Consider, for instance, the medium as purely classical. Then for the parameters of the whole reaction we obtain the following equations:
λcl ) λscl + λs1q λ1q ) λd1q + λa1q λ2q ) λd2q + λa2q In other words, the medium’s quantum mode contributions to λ1q and λ2q are omitted and their reorganization energy is included into λcl. Similarly, if inner-sphere modes are considered as classical, then
λcl ) λscl + λd1q + λa1q λ1q ) λs1q λ2q ) λs2q As soon as λcl, λ1q, and λ2q are found, parameters ωq and Sq are calculated by eqs 2.5 and 2.6, and the rate is obtained from eq 1.1. The calculated rates of the flavin-tryptophan reaction are shown in Figure 4. The thick solid line includes the quantum effects from donor, acceptor, and medium. The thick dashed line is the rate when all inner-sphere and outer-sphere vibrations are classical. These lines are the same as in Figure 3. The thin dashed line is calculated for a classical medium whereas the inner-sphere quantum vibrations are properly described. The thin dotted line is for classical donor and acceptor vibrations whereas water’s quantum contribution is included. It is seen that the quantum effect from water is of minor significance as compared to that from inner-sphere vibrations. The results for the tryptophan-tyrosine reaction are very similar (not shown). 5. Conclusions In this paper, we have shown that the quantum effect of a continuous dielectric medium can be reasonably well approximated by the same Jortner formula, eq 1.1, as that of the cofactors provided that the quantum parameters are chosen correctly. The method to calculate these parameters was
J. Phys. Chem. B, Vol. 104, No. 29, 2000 6901
Figure 4. Calculated rates for the 4FADH• + TrpH reaction. Thick (solid and dashed) lines are the same as in Figure 3. Thin dashed line: classical medium. Thin dotted line: classical inner-sphere vibrations (see text).
discussed in this paper. We compared our approximation with the exact rate and found that the most of the quantum effect is well reproduced in a wide range of the driving force. As an example, we have calculated the quantum parameters for three cofactors of DNA photolyase and water as a representative medium with a relatively large quantum effects as compared to other media. The quantum effect from water is found to be insignificant, on the background of that from the cofactors, unless the reaction is in the deep inverted region. We expect that the quantum effect of the actual protein medium will be smaller than that of water, so that our conclusion is that the medium mostly plays its specific classical role whereas the cofactors mainly contribute to the quantum effect, giving nearly nothing to the classical reorganization energy. It is worthwhile to outline the similarities and differences between our approach and the traditional method used, e.g., in the recent paper by Kumar et al.3 Both methods are based on the Jortner formula, and both of them calculate the internal reorganization energy by ab initio methods. The difference is in the way the second parameter, the quantum frequency ωq, is found. Kumar et al. use only an estimate of ωq. They indentify the modes contributing mostly to geometry changes (carboncarbon lengths in their case), and then the average frequency of these modes, 1400 ( 200 cm-1, is assigned to ωq. In our systems we found that the dispersion of quantum frequencies is much larger than 400 cm-1,53 and in this case it would be difficult to chose a single frequency using only qualitative arguments, as it was done in ref 3. In contrast to the qualitative estimate, which can be used only when the dispersion is small, our method provides a rigorous procedure which for any dispersion allows one to calculate an optimal frequency ωq that results in the best approximation of the exact rate with the Jortner expression in a wide range of ∆G0. Shifts of frequencies and Duchinsky rotation are taken into account in our approach. Comparison of different quantum chemical methods showed that electronic correlations strongly affect the calculated quantum parameters for open-shell cofactors and to a lesser extent for closed-shell ones. Therefore, the DFT methods were found to be more appropriate for open-shell systems. Our method enabled us to make an estimate of the quantum effects in a model which describes two redox reactions in DNA photolyase. For the purpose of having a complete picture, we investigated these reactions at different driving forces and found that the quantum effect in the inverted region is sufficiently
6902 J. Phys. Chem. B, Vol. 104, No. 29, 2000 strong to flatten the rate vs driving force dependence. We have not yet estimated a similar effect caused by inelastic tunneling,10,11 which works in the same direction. The effect of the rate flattening in the inverted region may be utilized in biology to ensure stable performance of ET systems, i.e., to make the ET rate insensitive to small driving force variations. In the example discussed in this paper we find that indeed the quantum effects from redox cofactors are strong enough to stabilize the rate in a vicinity of its maximum. Recent experiments by Gray et al.52 showed no inverted region for Ru-modified azurins, which may result from both the quantum effects and inelastic tunneling. Acknowledgment. We are grateful to X. Song for providing us with data files for dielectric constant of water. This work has been supported by the research grants from the National Institutes of Health (GM54052-02) and NIH Fogarty International Center (1 R03 TW00954-01). A.A.S. acknowledges the fellowships from the Sloan and Beckman Foundations, and E.S.M. acknowledges the support of the Russian Foundation for Basic Research (98-03-33155a). References and Notes (1) Jortner, J. J. Chem. Phys. 1976, 64, 4860. (2) Closs, G. L.; Calcaterra, L. T.; Green, N. J.; Penfield, K. W.; Miller, J. R. J. Phys. Chem. 1986, 90, 3673. (3) Kumar, K.; Kurnikov, I. V.; Beratan, D. N.; Waldeck, D. H.; Zimmt, M. B. J. Phys. Chem. A 1998, 102, 5529. (4) Kubo, R.; Toyozawa, Y. Prog. Theor. Phys. 1955, 13, 160. See also: Perlin, Yu. E. SoV. Phys. Usp. 1964, 6, 542. (5) Lee, E.; Medvedev, E. S.; Stuchebrukhov, A. A. J. Chem. Phys. 2000, 112, 9015. (6) Ovchinnikov, A. A.; Ovchinnikova, M. Ya. SoV. Phys. JETP 1969, 29, 688. (7) Song, X.; Marcus, R. A. J. Chem. Phys. 1993, 99, 7768. (8) Siders, P.; Marcus, R. A. J. Am. Chem. Soc. 1981, 103, 741. (9) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265. (10) Daizadeh, I.; Medvedev, E. S.; Stuchebrukhov, A. A. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 3703. (11) Medvedev, E. S.; Stuchebrukhov, A. A. J. Chem. Phys. 1997, 107, 3821. (12) Ungar, L. W.; Newton, M. D.; Voth, G. A. J. Phys. Chem. B 1999, 103, 7367. (13) Vath, P.; Zimmt, M. B.; Matyushov, D. V.; Voth, G. A. J. Phys. Chem. B 1999, 103, 9130. (14) Sancar, A. Biochemistry 1994, 33, 2. (15) Heelis, P. F.; Hartman, R. F.; Rose, S. D. Chem. Soc. ReV. 1995, 24, 289. (16) Heelis, P. F.; Kim, S. T.; Okamura, T.; Sancar, A. J. Photochem. Photobiol. B: Biol. 1993, 17, 219. (17) Okamura, T. A.; Sancar, A.; Heelis, P. F.; Begley, T. P.; Hirata, Y.; Mataga, N. J. Am. Chem. Soc. 1991, 113, 3143. (18) Antony, J.; Medvedev, D. M.; Stuchebrukhov, A. A. J. Am. Chem. Soc. 2000, 122, 1057. (19) Essenmacher, C.; Kim, S. T.; Atamian, M.; Babcock, G. T.; Sancar, A. J. Am. Chem. Soc. 1993, 115, 1602. (20) Kim, S. T.; Sancar, A.; Essenmacher, C.; Babcock, G. T. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 8023. (21) Kim, S. T.; Heelis, P. F.; Sancar, A. Methods Enzymol. 1995, 258, 319. (22) Cheung, M. S.; Daizadeh, I.; Stuchebrukhov, A. A.; Heelis, P. F. Biophys. J. 1999, 76, 1241.
Lee et al. (23) Heelis, P. F.; Okamura, T. A.; Sancar, A. Biochemistry 1990, 29, 5694. (24) Kim, S. T.; Sancar, A.; Essenmacher, C.; Babcock, G. T. J. Am. Chem. Soc. 1992, 114, 4442. (25) Sanders, D. B.; Wiest, O. J. Am. Chem. Soc. 1999, 121, 5127. (26) Aubert, C.; Mathis, P.; Eker, A. P. M.; Brettel, K. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 5423. (27) Sigel, H. D.; Sigel, A. Metalloenzymes InVolVing Amino Acid Residue and Related Radicals, Metals Ions in Biological Systems; Dekker: New York, 1994; Vol. 30. (28) Stubbe, J.; van der Donk, W. A. Chem. ReV. 1998, 98, 705. (29) Okamura, T. A.; Sancar, A.; Heelis, P. F.; Hirata, Y.; Mataga, N. J. Am. Chem. Soc. 1989, 111, 5967. (30) Heelis, P. F.; Sancar, A.; Okamura, T. A. J. Photochem. Photobiol. B: Biol. 1992, 16, 387. (31) Gaussian 94, Revision E.2; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian, Inc.: Pittsburgh, PA, 1995. (32) Moonen, C. T.; Vervoot, J.; Mu¨ller, F. Biochemistry 1984, 23, 4859, 4868. (33) Zheng, Y. J.; Ornstein, R. L. J. Am. Chem. Soc. 1996, 118, 9402. (34) Meyer, M.; Hartwig, H.; Schomburg, D. J. Mol. Struct. (THEOCHEM) 1996, 364, 139. (35) Kurreck, H.; Bock, M.; Elsner, B. M.; Kraus, H.; Lubitz, W.; Mu¨ller, F.; Geissler, J.; Kroneck, P. N. H. J. Am. Chem. Soc. 1984, 106, 737. (36) Hausen, H.-D.; Schulz, A.; Kaim, W. J. Chem. Soc., Perkin Trans. 1993, 2, 343. (37) Platenkamp, R. J.; Palmer, M. H.; Visser, A. J. W. G. J. Mol. Struct. 1980, 67, 45. (38) Li, Y. F.; Heelis, P. F.; Sancar, A. Biochemistry 1991, 30, 6332. (39) Hoffman, B. M.; Roberts, J. E.; Kang, C. H.; Margoliash, E. J. Biol. Chem. 1981, 256, 6556. (40) Sivaraja, M.; Goodin, D. B.; Smith, M.; Hoffman, B. M. Science 1989, 245, 738. (41) Houseman, A. L. P.; Doan, P. E.; Goodin, D. B.; Hoffman, B. M. Biochemistry 1993, 32, 4430. (42) Huyett, J. E.; Doan, P. E.; Burbiel, R.; Houseman, A. S. P.; Sivaraja, M.; Goodin, D. B.; Hoffman, B. M. J. Am. Chem. Soc. 1995, 117, 9033. (43) Sahlin, M.; Lassmann, G.; Po¨tsch, S.; Slaby, A.; Sjo¨berg, B.-M.; Gra¨slund, A. J. Biol. Chem. 1994, 269, 11699. (44) Lendzian, F.; Sahlin, M.; MacMillan, F.; Bittl, R.; Fiege, R.; Po¨tsch, S.; Sjo¨berg, B.-M.; Gra¨slund, A.; Lubitz, W.; Lassmann, G. J. Am. Chem. Soc. 1996, 118, 8111. (45) Krauss, M.; Garmer, D. R. J. Phys. Chem. 1993, 97, 831. (46) Jensen, G. M.; Goodin, D. B.; Bunte, S. W. J. Phys. Chem. 1996, 100, 954. (47) Qin, Y.; Wheeler, R. A. J. Am. Chem. Soc. 1995, 117, 6083. (48) It is important to run both jobs with the NoSymm keyword in order to prevent GAUSSIAN from reorienting the coordinate system between the initial and final states. (49) Medvedev, E. S.; Osherov, V. I. Radiationless Transitions in Polyatomic Molecules; Springer: Berlin, 1995; Chapter 6. (50) Harvey, S. C.; Hoekstra, P. J. Phys. Chem. 1972, 76, 2987. (51) Jordanides, X. J.; Lang, M. L.; Song, X.; Fleming, G. R. J. Phys. Chem. B 1999, 103, 7995. (52) Mines, G. A.; Bjerrum, M. J.; Hill, M. G.; Casimiro, D. R.; Chang, I.-J.; Winkler, J. R.; Gray, H. B. J. Am. Chem. Soc. 1996, 118, 1961. (53) We found that for flavin 50% of the quantum reorganization energy comes from modes with frequencies 300-700 cm-1 and 50% from 1200 to 1800 cm-1. For tryptophan 30% come from 400 to 500 cm-1 and 70% from 1000 to 1700 cm-1. For tryptophan cation radical 15% come from 600 to 1000 cm-1 and 85% from 1200 to 1600 cm-1. For tyrosine, 25% come from 400 to 700 cm-1 and 75% from 1000 to 1700 cm-1.
