Protein dynamics and electron transfer: Electronic decoherence and nonCondon effects
See allHide authors and affiliations

Edited by Harry B. Gray, California Institute of Technology, Pasadena, CA (received for review December 6, 2004)
Abstract
We compute the autocorrelation function of the donoracceptor tunneling matrix element 〈T_{DA} (t)T_{DA} (0)〉 for six Ruazurin derivatives. Comparison of this decay time to the decay time of the timedependent FranckCondon factor {computed by Rossky and coworkers [Lockwood, D. M., Cheng, Y.K. & Rossky, P. J. (2001) Chem. Phys. Lett. 345, 159165]} reveals the extent to which nonCondon effects influence the electrontransfer rate. 〈T_{DA} (t)T_{DA} (0)〉 is studied as a function of donoracceptor distance, tunneling pathway structure, tunneling energy, and temperature to explore the structural and dynamical origins of nonCondon effects. For azurin, the correlation function is remarkably insensitive to tunneling pathway structure. The decay time is only slightly shorter than it is for solventmediated electron transfer in small organic molecules and originates, largely, from fluctuations of valence angles rather than bond lengths.
The interplay among nuclear motion and electronic dynamics is the subject of increasing focus in the field of electron transfer (ET) processes (18). Recent research has focused on the effects of bridge nuclear motion on ET, with chemical, biological, and electronic device applications (see refs. 911 for reviews). Early theoretical analysis indicates that tunneling matrix element modulation by bridge dynamics can alter the free energy dependence of ET reaction rates by causing the BornOppenheimer (12) and FranckCondon approximations to fail (1315). More recently, theoretical studies explored ET kinetics in systems with fluctuating donoracceptor matrix elements (1624). Bridge motion can cause large and rapid donoracceptor matrix element fluctuations, affecting the tunneling pathway structure and the interferences among pathways (2537). The coupling matrix element fluctuations may have large contributions from solventpolarization fluctuations (23), and these fluctuations facilitate electronically forbidden and gated ET (3847). Finally, bridgenuclear relaxation creates inelastic tunneling pathway channels (11, 20, 24, 4749) that can change the mechanism of ET from superexchange to resonant tunneling to sequential hopping (11, 4144, 5062) and can lead to breakdown of the BornOppenheimer approximation (63, 64).
The goal of this work is to characterize tunneling matrix element fluctuations in azurin and, in particular, to examine their influence on the ET rate and on the validity of the FranckCondon approximation. FranckCondon breakdown can reduce the ET rate in the case of activationless ET reactions and enhance the rate for activated ET (18). ET in Rumodified azurin is nearly activationless, and the protein is often approximated as being a rigid medium for tunneling because the tunneling pathways traverse a β sheet. In this work, we compute the effects of tunneling matrix element fluctuations on the rate as a function of distance, temperature, protein structural fluctuations, and intervening pathway structure. Further, we identify the types of motion that cause the coupling to fluctuate.
The general derivation of the nonadiabatic rate expression for fluctuating donoracceptor matrix elements cannot assume the validity of the FranckCondon separation. As explained below, the FranckCondon approximation is only applicable when matrix element fluctuations are slower than the decay time of the thermally weighted nuclear FranckCondon overlap factor. This decay time, τ _{FC} , measures the escape time of the system from regions of FranckCondon allowed transitions. τ _{FC} is roughly the time it takes for two nuclear wavepackets to lose overlap when they are initially placed at the crossing point between the donor and acceptor diabatic energy surfaces and each is allowed to move on a separate surface. , where λ is the total reorganization energy and T is the temperature (65, 66).
ET Rates and Correlation Functions
The argument above about the validity of the FranckCondon approximation can be derived from a timedependent formulation of the ET rate. In the vibronic state representation, the nonadiabatic rate is a thermally weighted sum of energyconserving transitions between initial D; v, a〉 and final A; w, b〉 BornOppenheimer vibronic levels, i.e., (17, 20, 24, 67), where D〉 and A〉 are the donor and acceptor electronic states, and v〉 and w〉 are vibrational levels associated with the donor and acceptor diabatic energy surfaces. These vibrations modulate the energies of D〉 and A〉 but not the tunneling barrier between them. a〉 and b〉 are vibrational levels associated with vibrations that modulate only the tunneling barrier (dominantly orientational modes for longdistance ET). The energies of D; v, a〉 and A; w, b〉 are E_{D} + ε _{v} + ε _{a} and E_{A} + ε _{w} + ε _{b} , respectively, where E_{D} and E_{A} denote the minimum electronic state energies and ε _{i} denotes vibrational state energies. The thermal probability of the initial vibronic state D; v, a〉 is taken to be the product P_{Dv}P_{a} of the thermal probabilities of the vibrational levels v〉 and a〉 (consistent with the assumption that motions affecting the donor and acceptor energies are independent of those affecting the tunneling barrier). The scattering operator T̂ in Eq. 1 is given by where V̂ is the donorbridge/acceptorbridge electronic coupling operator and Ĝ ^{br} is the bridge electronic/vibrational Green's function. T̂ describes tunneling transitions between the electronic D〉 and A〉 states that are accompanied by transitions between bridge vibrational states, i.e., 〈D; v, aT̂A; w, b〉 = 〈vw〉〈aT̂_{DA} b〉 where T̂_{DA} ≡ 〈DT̂A〉 is an operator in the space of bridge vibrational levels. This characteristic means that the tunneling electron can interact with the bridge vibrations (e.g., refs. 14, 20, 47, and 49).
By using the Fourier representations of the δfunction in Eq. 1 (68), the nonadiabatic rate is where is the thermal correlation function of the tunneling operator in Eq. 2 . T̂_{DA} (0) = 〈DT̂A〉 and , where is the bridge vibrational Hamiltonian. The term is the thermally weighted timedependent FranckCondon factor ( and are the vibrational Hamiltonians for the donor and acceptor energy surfaces).
Eq. 3 is a timedomain expression for the rate, equivalent to the energydomain expression in Eq. 1 (its Fourier transform). It provides a timedependent view of the Condon approximation breakdown. The time integral in Eq. 3 is the Fourier transform of the product of two decaying functions, and , evaluated at the DA energy gap frequency. To recover the standard expression for the FranckCondon rate, we replace C_{TDA} (t) in Eq. 3 by and factor it out of the integral, i.e., A necessary condition for the validity of this approximation is that C_{FC} (t) decays much more rapidly than C_{TDA} (t). The time integral in Eq. 3 is equal to the convolution integral of the Fourier transforms of C_{FC} (t) and C_{TDA} (t). If C_{FC} (t) decays much faster than C_{TDA} (t), then the width of C̃ _{FC}(ω) is much greater than the width of (assuming roughly Lorentzian or Gaussian envelopes for both functions). Therefore, in the frequency region around , where is nonzero, C̃_{FC} (ω) is a slowly varying function. Substituting in the convolution integral, and keeping only C̃_{FC}(ω_{DA}), gives (where ω′ = ω _{DA}  ω). This zerothorder expression is equivalent to Eq. 6 , because , and is the time integral of Eq. 6. Corrections to the FranckCondon rate can be derived by retaining higherorder terms in the Taylor expansion of C̃_{FC} (ω). NonCondon corrections generally depend on E_{D}  E_{A} , in addition to the decay times of C_{FC} (t) and C_{TDA} (t) (1719, 22, 24), because the energy gap determines a third time scale in Eq. 3 , the oscillation period of the exponential.
The advantage of the timedomain description (Eq. 3 ) is that the decay times of C_{FC} (t) and C_{TDA} (t) can be associated readily with structural fluctuations of the molecule and its solvent, by using molecular dynamics (MD) simulations and semiclassical approximations of the correlation functions. In the classical limit for bridge motion, the T̂_{DA} (t) matrix elements averaged over bridge vibrational states (Eq. 4 ) are replaced by a superexchange matrix element that is modulated by the classical bridge trajectories R(t) {i.e., T_{DA} (t) = T_{DA} [R(t)]}. In this case, the averages in Eqs. 4 and 5 become averages over classical trajectories.
Correlation Functions in Azurin
To our knowledge, the only computation of the shorttime decay of C_{FC} (t) for a biological ET system was carried out for azurin by Lockwood et al. (69). The decay was computed by using the frozen Gaussian approximation for the nuclear wavepackets (70, 71). The multidimensional nuclear wave functions were approximated as products of atomcentered Gaussian wavepackets of constant widths moving along classical trajectories (obtained by MD simulations). In this case, C_{FC} (t) is the thermal average of the timedependent overlap between an initial nuclear wavepacket propagated with the electron in the donor state and the same wavepacket propagated with the electron in the acceptor state. The time it takes to lose overlap is τ _{FC} . In this semiclassical approach, the rate in Eq. 3 is where T_{DA} (t) is the superexchange matrix element evaluated along classical nuclear trajectories, and 〈 〉 _{T} denotes the classical thermal average. Lockwood et al. (69) reported τ _{FC} = 34 fsec with contributions from both solvent and protein matrix relaxation. Their value is consistent with the estimate fsec, with λ ≃ 1 eV (1 eV = 1.602 × 10^{19} J) at room temperature. Although τ _{FC} is not the only time scale associated with the thermally weighted FranckCondon factor (65, 66), it seems that, at least in the case of azurin, it is the shortest time scale.
T_{DA} (t) fluctuations have been studied for many ET systems using MD simulations and quantum chemical semiempirical calculations (20, 25, 2737). In most cases, large fluctuations in T_{DA} were observed, sometimes on time scales as short as tens of femtoseconds. A measure of the fluctuations of T_{DA} (t) is the coherence parameter, defined as (32). For a flexible bridge with multiple interfering tunneling pathways R _{coh} approaches zero. The coherence parameter is the infinite time limit of the normalized autocorrelation function, R _{coh} = lim _{t→∞}C_{TDA}(t)/C_{TDA} (0) (in the classical limit for the nuclear motion). The coherence parameter does not contain information about how fast coherence is lost as a result of the fluctuations in the matrix element. This information is contained in the decay time of C_{TDA} (t), τ_{coh}. Analysis of the spinboson model with an exponentially decaying C_{TDA}(t) (C_{TDA}(t) ∝ exp[t/τ_{coh}]) has shown a complex dependence of the nonCondon rate on τ_{coh} (18). Matrix element fluctuations may enhance or reduce the ET rate (with respect to its static, infinite τ_{coh} value), depending on the relative magnitudes of the reorganization energy and the driving force (18).
Recently, a perturbative theory was developed (24) to calculate nonCondon corrections to the ET rate. The theory is based on the Taylor expansion of C̃_{FC} (ω) described in the previous section. It is equivalent to an expansion of the rate in terms of R _{coh} and τ_{coh}, and it gives a lowestorder correction to the classical Marcus rate equal to where k ^{(0)} is the Marcus rate, λ the reorganization energy, and ΔE ^{0} is the driving force of the reaction. k ^{(2)} also depends on (τ _{FC} )^{2}, as seen by substituting in Eq. 8 . R _{coh} and τ_{coh} were computed from MD/semiempirical simulations of a Cclamp donoracceptor molecule in different solvents (37). C_{TDA} (t) showed a decay time of τ_{coh} ≈ 0.1 psec and a coherence parameter close to zero. These values gave a negligible nonCondon correction to the rate (0.04%) based on Eq. 8 .
To our knowledge, this Cclamp molecule analysis is the only computation of C_{TDA} (t) for an ET system. Because FranckCondon breakdown depends on the relative decay times of C_{TDA} (t) and C_{FC} (t), it is important to explore C_{TDA} (t) for a biological ET system like azurin, where there is a wealth of experimental information on the rates (7276), as well as prior analysis of C_{FC} (t) (69). Here, we compute C_{TDA} (t) for the Rumodified azurin with Protein Data Bank ID code 1BEX (77) and some of its mutants. We use MD simulations and extendedHückel electronic structure calculations to compute C_{TDA} (t). Given the femtosecondtimescale decay of C_{FC} (t) (69), we focus on the short time decay of C_{TDA} (t) and interpret this decay in terms of protein structure and dynamics. Fig. 1 shows the structure of azurin and positions of the histidine residues linked to the Rubipyridine complex in the wild type (WT) and the derivatives as follows: H83G Q107H, H83G M109H, H83G K122H, H83G T124H, and H83G T126H. The donoracceptor pairs were chosen to study C_{TDA} (t) as a function of the following: (i) average (MDbased) donoracceptor distance (metaltometal), T122H (14 Å), WT (17 Å), M109H (19 Å), T124H (22 Å), Q107H (25 Å), and T126H (25 Å); (ii) distance along the same covalent pathway, (T122H vs. T124H and T126H); and (iii) intervening medium structure (covalent, hydrogenbonded, or mixed pathways) for approximately the same average donoracceptor distance (Q107H and T126H). Further, for all donoracceptor pairs, C_{TDA} (t) is computed for a range of tunneling energies to assess the sensitivity of the correlation function to D/A energetics.
By choosing donoracceptor pairs with similar immediate environments, we aim to isolate the effects of bridge structure and dynamics on C_{TDA} (t). Further, because of the approximations contained in MD simulations and extendedHückel calculations, we seek generic features of the correlation functions that are robust with respect to parameter choice (e.g., tunneling energy) and geometric constraints. Only the shorttime dynamics of C_{TDA} (t) are explored because the initial decay of the FranckCondon factor is very rapid.
MD Simulations. The simulated Pseudomonas aeruginosa azurin structures were based on the Protein Data Bank structure 1BEX (77). Five other structures with different locations of the Rubipyridine complex were modeled by performing a H83G mutation and Q107H, M109H, K122H, T124H, or T126H mutations using the sybyl molecular modeling package (Version 6.5, Tripos Associates, St. Louis). The Rubipyridine complex was added to each derivative by using the program vmd (www.ks.uiuc.edu/Research/vmd) and the 1BEX structure as a template. All structures were solvated in bulk TIP3 water, and Na^{+} and Cl^{} counterions were added to an ionic strength of 200 mM. The resulting structures included 18,26719,420 atoms and had dimensions of ≈55 × 57 × 53 Å^{3}.
The structures underwent energy minimization for 500 steps with harmonic constraints imposed on the protein backbone and Rubipyridine complex, followed by equilibration for 200 psec with the same constraints and 1 nsec without constraints at 310 K. Simulations were performed on up to 300 processors at the Duke University Center for Computational Science, Engineering and Medicine cluster (www.csem.duke.edu) and on the Biophysics Cluster at the University of Cyprus, by using the programs charmm (78) and namd (www.ks.uiuc.edu/Research/namd) with the charmm force field (79), isothermalisobaric (NpT) ensemble, periodic boundary conditions, a Langevin thermostat, and full electrostatics particlemesh Ewald calculations (80). The forcefield parameters for the copper center were taken from ref. 81, and the parameters for the Rubipyridine complex were taken from ref. 27. Protein conformation samples for the effective coupling calculations were obtained by performing three series of MD simulations (5,000 conformations saved every 1, 10, and 100 fsec, respectively) at 310 K for each derivative.
To explore how protein structural motion affects the correlation function decay time, additional MD simulations were performed for the WT 1BEX as follows: (i) Simulations with no constraints at 100, 30, 10, and 3 K (preceded by a 1nseclong equilibration at these temperatures); (ii) simulations at 310 K with SHAKE constraints on all bonds to suppress bondlength fluctuations [the program charmm (78) was used for this simulation, because implementation of the SHAKE algorithm in namd2 is limited to bonds with hydrogen atoms); and (iii) simulations at 310 K with SHAKE constraints on all bonds with hydrogen atoms and all valence angle force constants (including UreyBradley terms) scaled up by a factor of 1,000 to suppress fluctuations of both bond lengths and valence angles (these simulations were performed with integration time step of 0.1 fsec).
T_{DA} Calculations. In the superexchange regime, the effective tunneling coupling is given by refs. 82 and 83 where S̃ and H̃ are the orbital overlap and Hamiltonian matrices of the donor [index (d)], bridge [index (b)], and acceptor [index (a)], respectively. The calculations were performed by using Slatertype valence orbitals and an extendedHückel Hamiltonian. Calculations at this level are rapid and provide reasonable qualitative estimates of the effective tunneling matrix elements in proteins and small molecules (32, 35, 8390)
The calculations were performed by using code based on the program forticon8 (91), itpack2c routines (92), and parameters for the extendedHückel Hamiltonian and the Slater orbitals based on density functional calculations (93). To describe the extent of mixing the initial with final states onto the bridge, the effective electronic coupling was renormalized with the method described in refs. 35, 90, and 94 as follows: Ab initio calculations of Solomon and coworkers (95, 96) indicate that the azurin donor state is dominantly composed of a d_{xy} orbital of copper () and a p_{y} orbital of sulfur in Cys112 (). In our calculations, the molecules were oriented so that the CuS metal bond is aligned along the x axis, and the nitrogens in His46 and His117 are positioned in the xy plane [more precisely, close to the xy plane, because the strongest copper ligands form a slightly nonplanar structure (95, 96)]. The donor orbital was defined as . Because the t_{2g} states fluctuate in energy as the protein fluctuates, the effective electronic coupling was computed to each of the five dorbitals of ruthenium. The values of R _{coh} and τ_{coh} were computed for each dorbital and then averaged.
Discussion
Fig. 2 shows the autocorrelation function of T_{DA} , normalized by its initial value, the initial ensemble average of , . The averaging was performed over a trajectory of 5,000 × 1fsec MD snapshots. The coherence parameter is the long time limit value of , and τ_{coh} is approximately the amount of time it takes for to drop to 1/e of its initial value. in Fig. 2 is not exact; it contains both systematic and statistical error. The systematic error is related to the approximations of classical MD, time step, and extendedHückel electronic structure analysis. The statistical error is related to the use of finite duration trajectories for the averaging. This error increases rapidly with t because fewer points are used for the computation of for large t. Therefore, only the short time dynamics of is reliable, and τ_{coh} should be derived from these dynamics. The method used to compute τ_{coh} and its error is described in ref. 97. Given for a trajectory of length N (where ), the coherence time (in femtosecond units) is defined as , where the sum is terminated when (97). The corresponding error is , and it reflects the statistical error for the shorttime component of the autocorrelation function. An alternative definition for the coherence time is , where Δ is restricted to short times (98). In this case, τ_{coh} can be computed by plotting τ_{coh}(Δ) as a function of Δ until the first plateau is reached. We found that the two methods give similar results for τ_{coh}. Because of the considerable error in the autocorrelation function for large t, the coherence parameters were not derived by taking the large t limit of normalized autocorrelation functions. Rather, they were computed directly from the averages 〈T_{DA} 〉 and by setting . Error analysis for 〈T_{DA} 〉 and (and thus R _{coh}) is based on the block renormalization method described in ref. 99.
Fig. 3 shows the dependence of 〈T_{DA} 〉 (Upper) and (Lower) on the number N of MD conformations used to compute the averages (WT with acceptor on the Ru d_{x2y2} orbital). All conformations were chosen from the same MD trajectory (0.5nsec total length) at equal intervals. The first point in both graphs corresponds to averaging over 10 MD conformations separated by 50 psec, the final point to 5,000 conformations separated by 100 fsec. The averages and their errors stabilize as the number of points increases. At N = 5,000, the values are similar to those obtained from 5,000 conformations separated by 1 fsec. This result is not surprising because all of the computed coherence times are <100 fsec, i.e., a 5,000 × 1fsec trajectory contains at least 50 independent blocks. It should be noted that this kind of error analysis is often neglected in the computation of 〈T_{DA} 〉 and for proteins. Rather, emphasis is usually placed on the level of electronic structure theory. Even though a highlevel quantum chemical method may reduce the systematic errors of average values, improving electronic structure methods does not address statistical errors. Statistical errors can be large because of the limited number of data points used in the averaging as the computational cost of the quantum chemical calculations grows.
Tables 1, 2, 3 show the values of R _{coh} and τ_{coh} (and their statistical errors ε _{Rcoh} and ε _{τcoh} , respectively) averaged over the five donor Ru dorbitals. For each Ru dorbital, both R _{coh} and τ_{coh} were computed from a 5,000fsec trajectory (1fsec time step) at E _{tun} = 10.8 eV [the value used in previous studies of azurin (88, 90)]. Table 1 shows R _{coh} and τ_{coh} for the WT and the mutants at 310 K. The entry order in the table tracks donoracceptor centertocenter distance (the last column shows the average distances for >5,000 fsec). Table 2 shows R _{coh} and τ_{coh} for the WT as a function of temperature. Table 3 shows these values for the WT at 310 K and different constraints imposed on the protein. The first row of the table reports values obtained for unconstrained MD. The second row data were computed with constraints on all bond lengths (vide supra). The third row shows the values computed with constraints on lengths of bonds with hydrogen atoms and all valence angle force constants scaled up by a factor of 1,000.
The tables are representative of the results we obtained for different 5,000fsec trajectories. The mutants T122H, T124H, and T126H lie along the same covalent pathway, whereas Q107H and T126H have approximately the same donoracceptor distance but different intervening pathways (Fig. 1). Table 1 shows that the decay of 〈T_{DA}(t)T_{DA} (0)〉 for all structures is of the order of a few tens of femtoseconds, and the coherence parameter is of the order of 0.1 or less. Therefore, both R _{coh} and τ_{coh} are robust, changing little with donoracceptor distance and intervening medium. Table 2 shows that R _{coh} and τ_{coh} change little with temperature. An increase in R _{coh} is observed only at 30 K and lower temperatures (not probed in experiments), whereas τ_{coh} does not depend on temperature within the statistical error.
Table 3 summarizes how constraints on protein structural flexibility influence R _{coh} and τ_{coh}. Use of different constraints, along with the temperature dependence, provides a strategy to identify specific kinds of protein motion that control these two parameters. The T_{DA} (t) time series show the major fluctuations on the time scale of tens of femtoseconds. Much smaller amplitude fluctuations of the femtosecond time scale are superimposed on the major fluctuations. These two time scales correspond to valence angle vibrations and covalent bond vibrations, respectively. MD does not accurately describe the time dependence of bond vibrations, but it does give the correct order of magnitude for the rms deviations of bond lengths. Therefore, we expect that the magnitudes of matrix element fluctuations caused by bond vibrations, as predicted by MD, are of the right order of magnitude. The first kind of constraints used in our simulations was the SHAKE constraints on all bond lengths (including the donor and acceptor units). In the T_{DA} (t) time dependence, these constraints eliminate the weak (femtosecond time scale) fluctuations but have a negligible effect on the strong (tens of femtoseconds) fluctuations. Accordingly, as shown in Table 3, neither R _{coh} nor τ_{coh} is significantly affected by the SHAKE constraints, indicating that these parameters are not sensitive to covalent bond vibrations. The other kind of constraints included the SHAKE constraints on lengths of all bonds with hydrogen atoms, combined with scaling up all force constants for valence angles. Suppressing both the fastest bond vibrations and all valence angle vibrations, these constraints visibly increase the time scale of the major T_{DA} (t) fluctuations, consistently increasing τ_{coh} by almost 80% (Table 3). Importantly, the 1,000fold increase in the valence angle force constants has a similar effect on angle fluctuations as a 1,000fold decrease in temperature, i.e., cooling the system down to ≈0.3 K. This result explains why no significant increase in τ_{coh} was observed in the simulations at low temperatures: even at 3 K, angle fluctuations were large enough to cause dephasing. In addition to the direct effect on dephasing, valence angle fluctuations also drive dihedral angle fluctuations that, in turn, modulate the lengths of hydrogen bonds, thereby affecting the electronic interference among the dominant ET pathways and T_{DA} (32). Indeed, we observed a decrease of >30% in dihedral angle fluctuations caused by imposing the above constraints. However, these constraints have no significant effect on R _{coh}, suggesting that the latter is controlled not only by valence angle fluctuations but also by dihedral angle fluctuations, which dominate the slower motion of the protein structure.
Finally, we do not find a significant dependence of R _{coh} and τ_{coh} on the electron tunneling energy for relevant ET tunneling energies in azurin. A large effect of the tunneling energy on τ_{coh} was observed only when the tunneling energy was artificially brought into resonance with bridgecentered eigenstates. In this case, τ_{coh} dropped, and coherence was lost within 5 fsec. This situation is not expected to be relevant to the Rumodified derivatives examined here, but recent experiments with highly oxidizing ET reagents in proteins may access this regime (100, 101). Moreover, nearresonant ET may be relevant to DNA ET (8), as well as to unsaturated organic ET systems with low donor(acceptor)bridge energy gaps (41, 42).
Fast coherence loss in azurin arises from T_{DA} fluctuations on the time scale of tens of femtoseconds. The time scale for coherence loss is long compared with the FranckCondon time of 34 fsec for azurin computed by Lockwood et al. (69) (the latter is in agreement with the estimate fsec for λ ≈ 1 eV at room temperature). Therefore, we expect the FranckCondon approximation to Eq. 6 to be accurate for singlestep tunneling in azurin.
Interestingly, the values of R _{coh} and τ_{coh} in azurin are not very different from those computed (37) for ET in a much smaller Cclampshaped donoracceptor molecule. In the Cclamps, fluctuations in the superexchange interaction through solvent presumably lead to the decoherence. The time scale of this decoherence is within an order of magnitude of the electronic decoherence time computed here for the Rumodified azurins. By using Eq. 8 for azurin with an average τ_{coh} ≃ 30 fsec, λ ≃ 1 eV, and Δ E ^{0} ≃ λ, gives a correction of 0.8% of the Marcus rate. As in the case of ET in Cclamped molecules, for ET in azurin nonCondon effects are not significant.
The type of nonCondon effects considered here are relevant to ET reactions through avoided crossings. If ET involves an electronic transition at a conical intersection (102) (as may be the case in excitedstate ET), nonCondon effects are significant. In this situation, the BornOppenheimer approximation is not valid at the avoided crossing, and corrections to the superexchange expression for the tunneling matrix element must be considered (63).
Conclusions
Simulations of the tunneling matrix element autocorrelation function for Rumodified WT azurin and several of its mutants show that the decay of 〈T_{DA} (t)T_{DA} (0)〉 is rapid, with a time constant of a few tens of femtoseconds. The coherence parameter is of the order of 0.1 or less. This finding means that azurin is a fluctuating medium for tunneling with multiple interfering tunneling pathways. Within the statistical errors of our analysis, we cannot ascertain any significant dependence of R _{coh} and τ_{coh} on the donoracceptor distance or intervening medium structure or any sensitivity to the presence or absence of chemical bond vibrations. Importantly, decoherence occurs on the time scale of valence angle fluctuations, indicating the central role of valence angle fluctuations in T_{DA} dephasing. All computed τ_{coh} values are an order of magnitude longer than the decay time of the thermally weighted FranckCondon factor τ _{FC} , computed by Lockwood et al. (69) to be 34 fsec. Therefore, nonCondon effects are insignificant for the case of azurin in the deep tunneling regime. As tunneling energies approach the charge injection limit, the Condon approximation will likely fail.
The methodology described in this work is suited for the study of dynamicsfunction relationships for ET reactions. An interesting question that should be addressed, by using this methodology, is whether nonCondon effects and dephasing characteristics of proteins vary with their secondary (and tertiary) structures.
Acknowledgments
S.S.S. was supported by University of Cyprus Research Grant “Repair of UVdamaged DNA by DNA Photolyase: Insights from Molecular Dynamics and Electron Transfer Calculations.” D.N.B. was supported by National Institutes of Health Grant GM048043 and by the Duke Center for Computational Science, Engineering and Medicine.
Footnotes

↵ † To whom correspondence may be addressed. Email: skourtis{at}ucy.ac.cy or david.beratan{at}duke.edu.

Author contributions: S.S.S., I.A.B., T.K., and D.N.B. designed research; S.S.S., I.A.B., T.K., and D.N.B. performed research; S.S.S., I.A.B., T.K., and D.N.B. analyzed data; and S.S.S., I.A.B., T.K., and D.N.B. wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: ET, electron transfer; MD, molecular dynamics.
 Copyright © 2005, The National Academy of Sciences
References
 ↵

Kuznetsov, A. M. (1995) Charge Transfer in Physics, Chemistry and Biology (Gordon & Breach, Amsterdam).

Jortner, J. & Ratner, M., eds. (1997) Molecular Electronics (Blackwell Scientific, Oxford).

Kuznetsov, A. M. & Ulstrup, J. (1999) Electron Transfer in Chemistry and Biology (Wiley, Chichester, U.K.).

Jortner, J. & Bixon, M., eds. (1999) Electron Transfer: From Isolated Molecules to Biomolecules, Advances in Chemical Physics (Wiley, New York), Vols. 106107.

May, V. & Kühn, O. (2000) Charge and Energy Transfer Dynamics in Molecular Systems (WileyVCH, Berlin).

↵
Balzani, V., Piotrowiak, P., Rodgers, M. A. J., Mattay, J., Astruc, D., Gray, H. B., Fukuzumi, S., Mallouk, T. E., Haas, Y., de Silva, A. P. & Gould, I. R., eds. (2001) Electron Transfer in Chemistry (WileyVCH, Weinheim, Germany), Vols. IV.

↵
Schuster, G. B., ed. (2004) Topics in Current Chemistry: LongRange Charge Transfer in DNA (Springer, Berlin), Vols. 236237.

↵
Newton, M. D. (2001) in Electron Transfer in Chemistry, ed. Balzani, V. (WileyVCH, Weinheim, Germany), Vol. 1, pp. 363.

Skourtis, S. S. & Beratan, D. N. (2001) in Electron Transfer in Chemistry, ed. Balzani, V. (WileyVCH, Weinheim, Germany), Vol. I, pp. 109125.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Bixon, M. & Jortner, J. (2003) Russian J. Electrochem. 39 , 510.
 ↵
 ↵
 ↵

Daizadeh, I., Medvedev, E. S. & Stuchebrukhov, A. A. (1997) Proc. Natl. Acad. Sci. USA 94 , 37033708. pmid:9108041
 ↵

↵
Balabin, I. A. & Onuchic, J. N. (2000) Science 290 , 114117. pmid:11021791
 ↵
 ↵
 ↵
 ↵
 ↵

Schlag, E. W., Sheu, S.Y., Yang, D.Y., Selzle, H. L. & Lin, S. H. (2000) Proc. Natl. Acad. Sci. USA 97 , 10681072. pmid:10655485

↵
Schlag, E. W., Yang, D.Y., Sheu, S.Y., Selzle, H. L., Lin, S. H. & Rentzepis, P. M. (2000) Proc. Natl. Acad. Sci. USA 97 , 98499854. pmid:10954730
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Bialek, W., Bruno, W. J., Joseph, J. & Onuchic, J. N. (1989) Photosynth. Res. 22 , 1527.

↵
Levich, V. G. & Dogonadze, R. R. (1959) Dokl. Acad. Nauk SSSR 124 , 123126.
 ↵
 ↵
 ↵
 ↵
 ↵

Langen, R., Chang, I.J., Germanas, J. P., Richards, J. H., Winkler, J. R. & Gray, H. B. (1995) Science 268 , 17331735. pmid:7792598

Onuchic, J. N., Beratan, D. N., Winkler, J. R. & Gray, H. B. (1992) Annu. Rev. Biophys. Biomol. Struct. 21 , 249377.
 ↵
 ↵

↵
Brooks, B., Bruccoleri, R., Olafson, B., States, D., Swaminathan, S. & Karplus, M. (1983) J. Comput. Chem. 4 , 187217.
 ↵
 ↵
 ↵
 ↵

↵
Balabin, I. A. & Onuchic, J. N. (1998) J. Phys. Chem. 102 , 74977505.

Tan, M. L., Balabin, I. A. & Onuchic, J. N. (2003) Biophys. J. 86 , 18131819.
 ↵
 ↵

↵
Howell, J., Rossi, A., Wallace, D., Haraki, K. & Hoffmann, R. (1977) forticon 8 (ExtendedHückel Program) (Quantum Chemistry Program Exchange, Indiana Univ., Bloomington, IN), QCPE 11, No. 344.

↵
Kincaid, D. R., Respess, J. R., Young, D. M. & Grimes, R.G. (1999) itpack 2c: A Fortran Package for Solving Large Sparse Linear Systems by Adaptive Accelerated Iterative Methods (Univ. of Texas, Austin, TX).
 ↵
 ↵
 ↵

↵
Gerwith, A. & Solomon, E. I. (1988) J. Am. Chem. Soc. 107 , 38113819.

↵
Janke, W. (2002) in Quantum Simulations of Complex ManyBody Systems: From Theory to Algorithms, Lecture Notes, John von Neumann Institute for Computing Series, eds. Grotendorst, J., Marx, D. & Muramatsu, A. (John von Neumann Institute for Computing, Jülich, Germany), Vol. 10, pp. 423445.

↵
Frenkel, D. & Smit, B. (2002) Understanding Molecular Simulation: From Algorithms to Applications (Academic, San Diego).
 ↵
 ↵
 ↵

↵
Klessinger, M. & Michl, J. (1995) Excited States and PhotoChemistry of Organic Molecules (Wiley, New York).
Citation Manager Formats
Article Classifications
 Physical Sciences
 Chemistry
 Biological Sciences
 Biophysics