1618
J. Phys. Chem. B 2005, 109, 1618-1625
Ab Initio Based Calculations of Electron-Transfer Rates in Metalloproteins Tatiana R. Prytkova,† Igor V. Kurnikov,†,‡ and David N. Beratan*,† Departments of Chemistry and Biochemistry, Duke UniVersity, Durham, North Carolina 27708, and Department of Chemistry, Northwestern UniVersity, 2145 Sheridan Road, EVanston, Illinois 60208 ReceiVed: September 18, 2004
A long-standing challenge in electron-transfer theory is to compute accurate rates of long-distance reactions in proteins. We describe an ab initio Hartree-Fock approach to compute electronic-coupling interactions and electron-transfer rates in proteins that allows the favorable comparison with experiment. The method includes the following key features; each is essential for reliable rate computations: (1) summing contributions over multiple tunneling pathways, (2) averaging couplings over thermally accessible protein conformations, (3) describing donor and acceptor electronic structure explicitly, including solvation effects, and averaging coupling over multiple energy-level crossings of the nearly degenerate donor-acceptor ligand-field states, and (4) eliminating basis set artifacts associated with diffuse basis functions. The strong dependence of coupling on donor-acceptor distance and on pathway interferences causes large variations of the computed electroncoupling values with protein geometry, and the strongest coupled conformers dominate the electron-transfer rate. As such, averaging over thermally accessible conformers of the protein and of the redox cofactors is essential. This approach was tested on six ruthenium-modified azurin derivatives using the high temperature nonadiabatic rate expression and compared with simpler pathways, average barrier, and semiempirical INDO models. Results of ab initio Hartree-Fock calculations with a split-valence basis set are in good agreement with the experimental rates. Predicted rates in the longer-distance derivatives are underestimated by 3-8fold. This analysis indicates that the key ingredients needed for quantitatively reliable protein electron-transfer rate calculations are accessible.
1. Introduction Protein-mediated electron transfer (ET) plays a central role in biochemistry.1-5 Models developed over the past 15 years have aimed to predict the influence of a protein’s folded threedimensional structure on ET rates. Tunneling pathway6-8 and packing density models9,10 describe the qualitative differences between through-bond and through-space donor-acceptor interactions in proteins. Extended-Hu¨ckel11-14 and neglect of differential overlap (NDO)15,16 methods capture multipathway effects. Linking electronic structure schemes with geometry sampling incorporates the influence of protein conformational flexibility.17,18,23,24 Ab initio electronic structure methods improve the description of the redox cofactor electronic structure and of multiple pathway interferences.15,19,20 The goal of this paper is to assess whether a quantitative theoretical description of the extensive ground-state Ru-modified protein ET rates of Gray and Winkler21,22 can be constructed. We find that the answer is yes, but a consistent description requires that we (1) sum contributions from multiple tunneling pathways, (2) average couplings over thermally accessible protein conformations, (3) treat donor and acceptor electronic structure explicitly, including solvation effects, and average coupling over multiple energylevel crossings of the nearly degenerate donor-acceptor ligandfield states, and (4) eliminate basis set artifacts associated with diffuse basis functions. Many protein electron-transfer reactions involve donors and acceptors that are weakly coupled electronically. In this case, †
Duke University. Northwestern University. * Address correspondence to this author. E-mail:
[email protected].
‡
the ET rate is described approximately by the high-temperature nonadiabatic expression:1
kET )
2π p
1
(
|HDA|2exp -
x4πλkBT
)
(∆G0 + λ)2 4λkBT
(1)
HDA is the electronic coupling between donor and acceptor mediated by the protein, λ is the reorganization energy, ∆G° is the reaction free energy, kB is Boltzmann’s constant, and T is the temperature. The challenge of computing electronic coupling is a quantum chemical problem for large systems involving thousands of atoms. It is complicated by the dynamical flexibility of the protein and solvent, and by the smallness of HDA. Progress in computing electronic couplings has been made by many groups.7,8,15,20,23,25-27,30-39 Proteins are too large to permit direct ab initio donor-acceptor couplings using standard methods. Semiempirical approaches do allow analysis of full proteins, but the results are sensitive to parametrization,15 and we present comparisons between ab initio and semiempirical ZINDO results. Our approach here is based on computations of couplings using a protein-fragment approach. The calculations include the essential parts of the protein that mediate donoracceptor coupling. PATHWAY analysis,6 importance analysis,40 protein pruning,41 and tight-binding computations42,69 suggest strategies to identify the key mediating residues. Ab initio Hartree-Fock methods have been applied recently to compute electronic couplings in proteins and protein models.15,23,25,33 While these results are promising, a systematic assessment of the accuracy of these methods is needed. Moreover, it is not clear what errors are introduced by the varied
10.1021/jp0457491 CCC: $30.25 © 2005 American Chemical Society Published on Web 01/08/2005
Electron-Transfer Rates in Metalloproteins
Figure 1. Ru(bpy)2(im)HisX modified azurins.
treatments of structural fluctuations, structure truncation, solvation, and approximate descriptions (or neglect) of the cofactors. Our goal is to develop a systematic, practical, and reliable approach to compute donor-acceptor interactions and ET rates in proteins. We describe a method to compute donor-acceptor interactions and test it on a family of ruthenium-modified azurins of Pseudomonas aeruginosa studied by Winkler, Gray, and coworkers.5 The extensive structural and kinetic data available for these proteins make them an outstanding test bed for theory. Moreover, the Ru-modified azurins were designed so that the rates are essentially coupling limited. That is, the reaction free energy was matched to the reorganization energy. 2. Methods proteins44
The azurin were modified at the surface HisX position (X ) 83, 107, 109, 122, 124, 126) with Ru(bpy)2(im)(HisX)2+ (bpy ) 2,2′-bipyridine, im ) imidazole). ET occurs from Cu+ to Ru3+ in a flash-quench kinetic scheme.45 The reaction driving force was matched to the reorganization energy of approximately 0.8 eV,46 and as such, the rate is nearly activationless.21,47,48 We explore the influence of the ruthenium complex solvation, geometry, positions, and intervening protein medium on the donor-acceptor couplings and ET rates. Modeling Protein Structure and Dynamics. The ET rate calculations described below were performed on the six azurin derivatives (Figure 1). We began the analysis with the His83 azurin crystal structure.49 Hydrogen atoms were added to the protein structure using standard amino acid geometries as templates.50 We mutated the amino acid at the modification sites (as needed) to histidines using the Insight II program and attached Ru complexes at the N imidazole ring positions. Energy minimization was performed with the Amber94 force field for the ruthenated-protein structures to relax atom-atom contacts that arise upon protein modification. The structures were solvated by placing them into a box (55 × 55 × 55 Å) of about 5000 TIP3P water molecules. The Amber94 force field was used to describe the protein. Equilibrium bond lengths and valence angles for Cu and Ru groups were chosen from geometries reported in the X-ray structures.49 The force constants of bonds to the metal atoms were set to 100 kcal/(mol Å2), and the force constants in the valence angles involving the metal atoms were set to 63 kcal/(mol rad2). These values are similar to those used to model similar Ru and Cu complexes.28,29 Atomic charges on metals and ligands were based on ab initio calculations and Mulliken population analysis. Tables of atomic charges can be found in the Supporting Information. We performed 400
J. Phys. Chem. B, Vol. 109, No. 4, 2005 1619 ps MD simulations50-52 for each of the solvated structures. Periodic boundary and constant pressure conditions were used. We selected 10 snapshots from the last 250 ps of each MD trajectory (with an interval of 25 ps) for electronic structure calculations. Computing Donor-Acceptor Couplings. We performed ab initio calculations on pruned-protein structures. Earlier importance analysis40 and pruning methods relied on extended-Hu¨ckel or related methods.25,53,54 We adopted a simpler PATHWAYSbased approach. We first determined the strongest pathway from the donor (Cu) to acceptor (Ru). Then we identified all atoms that belong to coupling pathways with coupling values 50%100%, 5%-100%, and 1%-100% of the strongest pathway coupling value. We then built “subproteins” that include all atoms in these classes. In the fragment extraction procedure, single bonds that are cut are terminated with H atoms. To avoid cutting multiple bonds or fragmenting delocalized (e.g., amide) or aromatic units, we increased the number of atoms in the pruned proteins to include all atoms of these chemical groups. We used energy-splitting methods to compute the donoracceptor interaction, HDA, in the Condon approximation.55-57 We assume that electron transfer occurs in an activated state with degenerate donor and acceptor localized electronic states (φD and φA). In this state, φD and φA are mixed via bridge states (superexchange). The unmixed system is quasi-degenerate, with approximate eigenstates ψ0 ) (1/x2)(φD + φA) and ψ1 ) (1/x2)(φD - φA). |HDA| is one-half of the energy splitting (|E1 - E0|) between these eigenstates. This splitting is computed from the HOMO and HOMO-1 energies of the N + 1 electron system (eq 2). This approach is based on Koopmans’ theorem.55-57 (N+1)e -EHOMO ) IP1 ) E Ne 0 - E0
(2a)
(N+1)e -EHOMO-1 ) IP2 ) E Ne 1 - E0
(2b)
Ne (EHOMO - EHOMO-1) ) E Ne 1 - E 0 ) 2|HDA|
(2c)
Ne In eq 2, IP represents the ionization potential, E Ne 0 and E 1 are the energies of the quasi-degenerate ground and first excited states of the electron-transfer system, and E (N+1)e is the 0 ground-state energy of the system with an extra electron. After constructing PATHWAYS-based pruned proteins, we performed Hartree-Fock calculations with a modified 3-21G basis set using the Gaussian program.58 The modified basis set does not have diffuse s and p basis functions on Cu. These diffuse functions introduce significant coupling artifacts at long distance (observed earlier,59,60 vide infra). We projected the effective Hamiltonian (Fock matrix) to the valence atomic orbital space15 to reduce the computational cost. This procedure gives results in quantitative agreement with full calculations. The ruthenium d-orbitals are closely spaced in energy, and their energy ordering changes as the protein geometry fluctuates. Therefore, we computed the root-meansquare (rms) value of the donor-acceptor couplings to the three upper occupied ruthenium-localized orbitals. The geometries taken from the MD trajectories and truncated based on PATHWAY analysis does not generally place the orbitals into an activated quasi-degenerate geometry. To move them into degeneracy, we applied an electric field along the donoracceptor axis. This field modifies the effective Hamiltonian matrix by
H′(DA) ) Heff(0) + [dxex + dyey + dzez]DA
(3)
1620 J. Phys. Chem. B, Vol. 109, No. 4, 2005
Prytkova et al.
Figure 2. Energy of MOs vs electric field.
where ex, ey, and ez are component of the unit vector along the Cu-Ru axis, DA is the strength of the electric field, and dx, dy, and dz are the electric dipole matrices in the valence-orbital basis. Similar methods were used previously by Larsson, Stuchebrukhov, and Newton.14,61-63 We computed the dependence of the energy eigenvalues on the electric field strength and computed |HDA| as one-half of the minimum energy splitting found as DA was varied. The energy eigenvalue dependence on electric field for Ru(bpy)2(im)(His83) azurin is shown in Figure 2. ET occurs between the Cu- and Ru-localized orbitals. In our electronic structure calculations, we choose a formal charge on Ru equal to +2 (instead of +3) to model the system using Koopmans’ theorem as described above. The electric field moves the Cuand Ru-localized orbitals into degeneracy. The upper positive slope line (Cu) and the highest negative slope line (Ru) cross when the electric field is 4 × 10-4 au/e (au is hartree/bohr). In this crossing region, the orbitals that were localized on Ru and Cu are mixed. The HOMO at zero field is localized on Cu. The difference in energy between the upper occupied orbital on Cu and the HOMO-1 orbital is much larger than kT, so we assumed that ET occurs mostly from the HOMO of Cu. For the Ru complex, the energy of the three orbitals is sensitive to the geometry of the complex. The complex is not rigid, and geometry changes switch the orbital energy ordering. Because this energy reordering is rapid on the ET rate time scale, we compute the mean-square value for the crossing of Cu and Ru orbitals, 〈|HDA|2〉, for the three crossings indicated in Figure 2. Modeling Solvation. The relative energies of the donor and acceptor depend strongly on their interactions with the solvent. We account for the influence of the solvent in our coupling calculations using point charge models. In Ru-modified azurin, the Ru complex acceptor is positively charged (+2/+3) and exposed to solvent. To describe the energetics of the acceptor state even qualitatively, this solvent interaction must be included. Our simplest model places a set of five point charges at the non-protein nitrogen atoms surrounding the Ru (-0.4e on each N) to compensate the +2 charge on Ru. This is an extremely simple solvation model. To build a more realistic model, we compute the solvent-induced charges using continuum electrostatics. To compute the induced charge distribution at the interface between the molecule and the solvent, we use
Find ) -Fext -
1 2 ∇φ 4π
(4)
where Fext is the external charge density and Find is the induced (polarization) charge density. Equation 4 is a form of the Poisson equation64
Ftot ) -
1 2 ∇φ 4π
(5)
Figure 3. Induced charge density localized on the surface of the Ru complex for His83 contour level of charge density 0.001 e/Å3
Figure 4. Cu-Ru distance along the MD trajectory.
where Ftot is the total charge density and φ is the computed electrostatic potential. The induced charge density Find is then included as a set of point charges in the protein electronic structure calculations. Figure 3 shows that the induced charge density is localized mostly on the surface of the Ru complex. 3. Results and Discussion Dependence of the Computed D/A Couplings on the Protein Conformation, Dynamics, and Fragment Size. Classical molecular dynamics (MD) trajectories of the ruthenated azurins reveal a substantial mobility of the attached Ru(bpy)2(im)(HisX)2 + units. The Cu-Ru distance fluctuates by about 2.0 Å on the picosecond time scale (see Figure 4). The Ru complex also rotates around the N(His)-Ru bond. We computed the electronic coupling for 10 snapshots with an interval of 25 ps from the last 250 ps along the MD trajectories of each azurin derivative. We built protein fragments using the PATHWAYS pruning procedure described above (5% cutoff criterion). The protein fragments for this cutoff criterion had 180-280 atoms, including hydrogens. The rms donoracceptor couplings were calculated as described above. Nonadiabatic electron-transfer rates were computed for each fragment using eq 1 and assuming -∆G° ) λ ) 0.8 eV. In these calculations, solvent polarization was modeled using point charges on the five non-protein nitrogen atoms coordinating Ru, as described above. The ET rates computed for 10 snapshots were averaged and are shown in Figure 5. Note that the conformational average is really a double average. That is, we average over the accessible protein geometries as well as the multiple ligand-field state crossings (Figure 2) for each of
Electron-Transfer Rates in Metalloproteins
J. Phys. Chem. B, Vol. 109, No. 4, 2005 1621
Figure 5. Theory vs experiment for 5% cutoff fragments averaged over snapshots and over the ligand-field states.
TABLE 1: Rms |HDA| Values Calculated for Six Derivatives of Azurina modified protein
rms |HDA| (eV)
theor ET rate (s-1)
expt ET rate (s-1)
Ru(bpy)2(im)(His83)Az Ru(bpy)2(im)(His107)Az Ru(bpy)2(im)(His109)Az Ru(bpy)2(im)(His122)Az Ru(bpy)2(im)(His124)Az Ru(bpy)2(im)(His126)Az
7.19 × 10-6 5.74 × 10-8 6.22 × 10-6 1.92 × 10-5 4.37 × 10-7 4.42 × 10-8
8.68 × 105 5.53 × 101 6.50 × 105 6.19 × 106 3.20 × 103 3.28 × 101
1.00 × 106 2.4 × 102 8.5 × 105 7.1 × 106 2.2 × 104 1.3 × 102
The rms |HDA| values are averaged on the MD trajectory. The ET rates are calculated assuming an activationless reaction. a
Figure 6. Computed vs experimental ET rates for the six azurin derivatives at specific MD snapshot geometries.
these geometries. Computed values of 〈|HDA|2〉1/2 appear in Table 1. Agreement between the computed average ET rates and the experimental rates is good. Excellent agreement is observed for the His83, His109, and His122 derivatives, and the computed average rates for the other derivatives with longer Cu-Ru distance are smaller than the experimental values. The largest difference between theory and experiment (about a factor of 7 in the ET rate) is observed for the His124 derivative. Figure 5 summarizes the main results of this paper. The figure indicates that ab initio split-valence Hartree-Fock calculations of the electronic coupling averaged over MD snapshots (including solvation effects modeled with point charges) succeed in making order of magnitude electron-transfer rate estimates. It should be noted that the ET rates computed for different MD snapshots of the same azurin derivative differ by as much as 3 orders of magnitude. This is shown in Figure 6, where ET rates computed for specific snapshots are compared with experimental rates. Averaging over geometries is critical for a quantitative description of experiments, as suggested in earlier studies.17,18,69 Importantly, the ET rate is dominated by a minority of modified protein geometries with the strongest ET
Figure 7. Computed ET rates vs the number of snapshots.
TABLE 2: Calculated |HDA| Values (eV) for Several Snapshot Geometries of the His83 Derivative of Azurina His 83
50% cutoff
5% cutoff
1% cutoff
snapshot 1 snapshot 2 snapshot 3 snapshot 4 snapshot 5 snapshot 6 His83 X-ray structure rms |HDA|
1.25 × 10-5 7.92 × 10-6 6.11 × 10-6 1.07 × 10-6 8.37 × 10-7 3.90 × 10-6 2.40 × 10-6 6.32 × 10-6
1.64 × 10-5 8.16 × 10-6 2.78 × 10-6 1.58 × 10-6 2.85 × 10-6 3.09 × 10-6 3.37 × 10-6 7.32 × 10-6
9.49 × 10-6 1.63 × 10-5 3.05 × 10-6 5.21 × 10-7 5.65 × 10-7 1.52 × 10-6 1.86 × 10-6 7.29 × 10-6
a
Values are calculated for different size protein fragments with 50%, 5%, and 1% pathway cutoff criteria.
couplings. This conclusion is consistent with a recent study of tunneling currents in scanning tunneling experiments on flexible systems70 as well as studies of interprotein ET.71 Electrontransfer rates computed for a randomly chosen thermally accessible geometry of the system most often will “predict” an ET rate substantially smaller than the experimental value. Since tunneling rates vary approximately exponentially with distance, and particularly small coupling values may arise as a result of destructive interference, modest geometry changes can have a large effect on the coupling. The large variations of the coupling with geometry can cause convergence of the rms HDA values to be slow. Figure 7 shows a running average of the ET rate computed using 20 snapshots of Ru-azurin derivative 83. The convergence in this specific case is rapid, and 10 snapshots (the number used in most of the our calculations) should be sufficient for convergence of the computed ET rate within 50%. Recall that 10 snapshots averages represent rms averaging of 30 invidual donor/acceptor coupling values (representing interaction between highest Cu orbital and three Ru orbitals). More extensive calculations will be needed to explore the convergence of ET rates as a function of the number of conformers included in the analysis. To explore the dependence of coupling on the protein flexibility and fragment size further, we calculated ET coupling on fragments of the azurin His83 structure using the PATHWAY pruning cutoff values of 50%, 5%, and 1% (vide supra). The number of atoms in the fragments was approximately 120, 240, and 340, respectively. Seven geometries of the azurin His83 derivative were examined: six MD snapshots and the reported X-ray structure. Table 2 shows the computed electronic coupling elements for these geometries using the different cutoff criteria. We also computed rms values of the electronic couplings for each criterion. Electronic couplings computed for the same geometry of the protein with different fragment sizes can differ substantially, and the couplings can converge slowly with increasing fragment size. Indeed, order of magnitude changes
1622 J. Phys. Chem. B, Vol. 109, No. 4, 2005
Prytkova et al.
Figure 8. Change in the relative energies of donor, bridge, and acceptor molecular orbitals with and without modeled solvent. In the systems without solvent, the molecular orbitals of ruthenium are lower in energy than the bridge-localized orbitals.
TABLE 3: Rms |HDA| Values Calculated for Ru(bpy)2(im)(His83)Az Using the Simple Point Charge and Induced Charge Modelsa solvent model point charges solvent model induced charges solvent model
rms |HDA| (eV)
ET rate theory (s-1)
ET rate expt (s-1)
7.33 × 10-6
9.02 × 105
1.00 × 106
8.60 × 10-6
1.24 × 106
1.00 × 106
Figure 9. Computed vs experimental rates using single snapshots from the end of the 400 ps MD trajectory of the six derivatives using the 3-21G basis set. Note that the His107 and His126 species fall far from the correlation line. Even averaging over multiple snapshots does not move the His126 rate closer to the line.
a The theoretical ET rates assumes an activationless reaction and λ ) 0.8 eV.
in the computed |HDA| values were observed for some snapshots (see Table 2) when the fragment cutoff was changed from 50% to 5% and from 5% to 1%. This behavior indicates “delicate” interferences among the coupling pathways. While slow convergence with fragment size was observed for individual geometries, average (rms) values of the ET couplings for a given value of the fragment cutoff criterion do not vary much with the cutoff criterion. This important observation suggests that it is essential to consider several (5-10) representative thermally accessible geometries. It is then sufficient to limit fragment sizes to moderate PATHWAYS screened values (50% to 5% cutoff criteria). Solvent Effects. In addition to providing superexchange pathways of their own,72,73 water can strongly perturb donor and acceptor orbital energies, especially for exposed ions (i.e., the ruthenium acceptor). If solvent-ion interactions are neglected in Ru-azurin, the molecular orbitals of ruthenium are computed to be lower in energy than the bridge-localized orbitals, and the computed electronic energy spectrum is qualitatively incorrect. It is not clear a priori how accurate an ion-solvent description is required to describe donor-acceptor superexchange interactions. We argue that a simple point charge model will suffice. Figure 8 shows how solvent modeling changes the relative energies of the donor, bridge, and acceptor centered orbitals. Here, the solvent polarization was modeled by placing negative point charges on the five nitrogen atoms coordinating the ruthenium. The solvent polarization drives the energies of the acceptor orbitals up by more than 5 eV, placing them in quasi-resonance with the donor orbitals and above the occupied bridge-centered orbitals, as expected. We also compared an induced-charge model described above to the simple point-charge calculation. The results obtained with the two models are similar (see Table 3). Basis Set Effects on the Couplings. Gaussian split-valence basis sets are known to be adequate56 to describe electronic couplings in covalent organic-bridged systems, although larger basis sets that include diffuse functions are needed to describe coupling across intermolecular gaps of 4 Å or more.15 We
Figure 10. Dependence of the computed ET coupling on the cutoff distance for long-range interactions among atomic orbitals.
TABLE 4: |HDA| Values and ET Rates Calculated for Ru(bpy)2(im)(His126)Az with a Standard 3-21G Basis Set His126 snapshot no.
|HDA| (eV)
theor ET rate (s-1)
expt ET rate (s-1)
1 2 3 4 5 average
8.35 × 10-6 2.11 × 10-6 1.72 × 10-6 8.79 × 10-7 1.50 × 10-6 2.08 × 10-6 (rms)
1.88 × 105 7.48 × 104 4.94 × 104 1.29 × 104 3.77 × 104 7.27 × 104
1.3 × 102 1.3 × 102 1.3 × 102 1.3 × 102 1.3 × 102 1.3 × 102
computed donor-acceptor electronic coupling values for MD snapshots along trajectories of the six azurin derivatives (one snapshot per derivative) using the standard split-valence 3-21G basis set. Results of the ET rate calculations using these couplings and an assumption of -∆G° ) λ ) 0.8 eV appear in Figure 9 and Table 4. The ET rates computed for derivatives with large Cu-Ru distances are substantially overestimated (by 2-3 orders of magnitude). We hypothesized that artificially long-range interactions are introduced by the diffuse orbitals in the basis set.59,60 To test this hypothesis, we zeroed out the one-electron effective Hamiltonian matrix elements between atomic orbitals centered on atoms separated by distances larger than a defined cutoff distance. We investigated the dependence of the computed ET coupling matrix element in the His126 derivative as a function of this cutoff distance (Figure 10). This strategy is equivalent to assuming that all interactions beyond the cutoff distance are mediated by superexchange rather than by direct interactions. We found that the computed ET coupling element decreased 20-fold when the cutoff distance decreased from 25 to 10 Å. Our previous studies15 indicate that electronic coupling is dominated by short-range interactions between atoms separated by no more than 6-7 Å. The unusually long-distance interaction elements in the effective Hamiltonian arise from diffuse s and p basis function in the 3-21G basis set for copper. Since local interactions that contribute to the superexchange are
Electron-Transfer Rates in Metalloproteins
J. Phys. Chem. B, Vol. 109, No. 4, 2005 1623
TABLE 5: Calculated |HDA| Values (eV) for Snapshot Geometries of the Azurin His83 Derivativea His 83
DFT
Hartree-Fock
snapshot 1 snapshot 2 snapshot 3 snapshot 4 snapshot 5 rms |HDA|
1.46 × 10-5 8.19 × 10-6 2.76 × 10-6 1.34 × 10-6 2.87 × 10-6 9.74 × 10-6
1.89 × 10-5 7.78 × 10-6 1.19 × 10-6 1.82 × 10-6 3.51 × 10-6 9.32 × 10-6
a Values are calculated with DFT (B3LYP) and Hartree-Fock methods.
not influenced substantially by the diffuse functions, we excluded them from the 3-21G basis set on copper. The computed ET couplings in the longest distance Ru-modified azurin derivatives were dramatically reduced as a result. The calculations reported throughout this paper were carried out with this modified 3-21G basis set. Hartree-Fock vs DFT Calculations of Donor/Acceptor Couplings. Density functionl theory (DFT) calculations with gradient-corrected functionals are similar to Hartree-Fock calculations in computational cost but give more accurate results for structure and vibrational frequencies. However, the relative advantages of DFT vs HF methods to compute donor/acceptor electronic couplings are not yet clear.78 Indeed, DFT may provide a poor description of electronic coupling between molecules separated at large distance because DFT does not reproduce the asymptotic long-distance decay of one-electron potentials and one-electron functions.77 In typical DFT approximations, the HF exchange potential is replaced by an approximate exchange potential, and, as a result, electron selfinteraction errors introduce a spurious energy shift to the occupied orbitals and corresponding changes to the long-distance asymptotic decay of the electronic states. Concerns regarding the use of DFT to describe charge-transfer excited states and electron-transfer reactions have been been dicussed.74,79 In longrange molecular electron transfer, the long-distance electronic wave function vacuum decay will not come into play because superexchange interactions dominate. One might expect, therefore, that DFT and HF calculations will give similar results.80 Table 5 compares HF and DFT (B3LYP functional) results for several MD snapshots of the His83 Ru-azurin species. Results obtained from DFT and HF calculations are essentially the same. Comparison of Hartree-Fock and Simple-Model Results. Figure 11 shows the experimental rates vs the PATHWAYS couplings and the MD averaged donor-acceptor distance. A linear correlation is observed between the averaged donoracceptor distance and the logarithm of the experimental ET rates, except for the His109 derivative. An average exponential decay constant β of 1.0 Å-1 is fit with a prefactor of 4.3 × 1013 s-1. A comparison of the PATHWAYS model results with experiment indicates that the usual set of PATHWAYS parameters probably overestimates the through-space decay. Improved agreement with experiment is obtained when the through-space (coupling element) decay exponent is reduced from 1.7 to 1.0 Å-1,10 and the smaller value is used in the figure. We also performed semiemprical calculations of donor/acceptor couplings in the azurin derivatives using the ZINDO/1 method75 implemented in the HARLEM program.50 ZINDO/1 calculations were performed using the point charge solvent model, averaging over ligand-field states and over the geometries (10 structures per derivative, the same fragments that were used in the Hartree-Fock analysis). This semiempirical method was chosen because of the availability of parameters for the transition metals.
Figure 11. (a) Exponential distance model averaged over the MD trajectory vs experimental rates. The decay factor fit is to β ) 1.0 Å-1. (b) Square of HDA calculated with the PATHWAY model (throughspace decay exponent of 1.0 Å-1) vs experimental ET rates.
Figure 12. Theory vs experiment for 5% cutoff fragments averaged over snapshots calculated with the ZINDO/1 method.
A comparison of computed and experimental electron-transfer rates based on ZINDO/1 couplings is shown in Figure 12. Computed electron transfer rates for His83, 122, and 126 azurins are in good agreement with experiment while rates for His107, 109, and 124 azurins are substantially overestimated (2-3 orders of magnitude). ZINDO/1 calculations likely do not capture the balance between through-bond and through-space electronic propagation. It is possible that tuning the semiempirical parameters that define the proportionality between off-diagonal matrix elements of the one-electron Hamiltonian and the overlap matrix could lead to better agreement. 4. Conclusions For the first time, systematic ab initio calculations of protein electron-transfer rates were performed for a series of ruthenated
1624 J. Phys. Chem. B, Vol. 109, No. 4, 2005 azurins and compared with high-quality experimental rate data. It was found that electronic coupling values based on HartreeFock calculations with a 3-21G (split-valence) basis set reproduce the coupling-limited experimental rates (within a factor of 7) when the rates are averaged over thermally accessible geometries and ligand-field states. The detailed mechanism of electron transfer in proteins remains a subject of intense interest. Waldeck and co-workers recently argued, for example, that ET rates in ruthenated proteins with small donor-acceptor distances (e.g., His122) can be limited by protein relaxation.76 The excellent agreement found here between computed and measured ET rates for the His122, His109, and His83 derivatives argues against relaxation limited rates. Large fluctuations of the electronic couplings (1-2 orders of magnitude) were computed for thermally accessible geometries. Computed couplings also depend strongly on the amount of protein included in the quantum calculations. The rms couplings, however, vary modestly with the amount of protein included. Strong basis set artifacts enter the computed donoracceptor electronic couplings at very large distances, as understood earlier.59,60 Caution should be exercised when computing donor-acceptor electronic couplings with diffuse functions in the basis set. We suspect that the strategy of computing geometry and ligand-field state averaged couplings will provide theoretical access to a number of important functional biological ET systems. Deviations between theory and experiment for derivatives with large Cu-Ru distances is significant enough to warrant further study in this regime. With this computational framework in hand, the challenge now is to improve the geometry sampling and the description of the donor-acceptor cofactors, including their solvation. Acknowledgment. We thank NIH (GM-048043) and Duke University for support of this research. We also thank H. B. Gray and J. R. Reimers for helpful discussions. Supporting Information Available: Table of atomic partial charges added to the Amber94 parameter set. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Bendall, D. S., Ed.; Protein Electron Transfer; Bios Scientific Publishers: Oxford, 1996. (2) Lippard, S. J.; Berg, J. M. Principles of Biorganic Chemistry; University Science Books: Mill Valley, CA, 1994. (3) Bertini, I.; Gray, H. B.; Lippard, S. J.; Valentine, J. S. Bioinorganic Chemistry; University Science Books: Mill Valley, CA 1994. (4) Balzani, V., Ed.; Electron Transfer in Chemistry; Wiley: Weinheim, Germany, 2001; Vols. 1-5. (5) Gray, H. B.; Winkler, J. R. Q. ReV. Biophys. 2003, 36, 341-372. (6) Beratan, D. N.; Betts, J. N.; Onuchic, J. N. Science 1991, 252, 1285-1288. (7) Regan, J. J.; Onuchic, J. N. AdV. Chem. Phys. 1999, 107, 497553. (8) Skourtis, S. S.; Beratan, D. N. AdV. Chem. Phys. 1999, 106, 377452. (9) Page, C. C.; Moser, C. C.; Chen, X. X.; Dutton, P. L. Nature (London) 1999, 402, 47-52. (10) Jones, M. L.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. A 2002, 106, 2002-2006. (11) Regan, J. J.; Risser, S. M.; Beratan, D. N. J. Chem. Phys. 1993, 97, 13083. (12) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 6111-6114. (13) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 13078-13082. (14) Daizadeh, I.; Gehlen, J. N.; Stuchebrukhov, A. A. J. Chem. Phys. 1997, 106, 5658-5666. (15) Kurnikov, I. V.; Beratan, D. N. J. Chem. Phys. 1996, 105, 95619573. (16) Priyadarshy, S.; Risser, S. M.; Beratan, D. N. J. Phys. Chem. 1996, 100, 17678-17682.
Prytkova et al. (17) Wolfgang, J.; Risser, S. M.; Priyadarshy, S.; Beratan, D. N. J. Phys. Chem. B 1997, 101, 2986-2991. (18) Daizadeh, I.; Medvedev, E. S.; Stuchebrukhov, A. A. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 3703-3708. (19) Babini, E.; Bertini, I.; Borsari, M.; Capozzi, F.; Luchinat, C.; Zhang, X. Y.; Moura, G. L. C.; Kurnikov, I. V.; Beratan, D. N.; Ponce, A.; DiBilio, A. J.; Winkler, J. R.; Gray, H. B. J. Am. Chem. Soc. 2000, 122, 45324533. (20) Stuchebrukhov, A. A. AdV. Chem. Phys. 2001, 118, 1-44. (21) Gray, H. B.; Winkler, J. R. Annu. ReV. Biochem. 1996, 65 537561. (22) Winkler, J. R.; Di Bilio, A. J.; Farrow, N. A.; Richards, J. H.; Gray, H. B. Pure Appl. Chem. 1999, 71, 1753-1764. (23) Zhang, L. Y.; Friesner, R. A.; Murphy, R. B. J. Chem. Phys. 1997, 107, 450-459. (24) Balabin, I. A.; Onuchic, J. N. Science 2000, 290, 114-117. (25) Kobayashi, C.; Baldridge, K.; Onuchic, J. N. J. Chem. Phys. 2003, 119, 3550-3558. (26) Stuchebrukhov, A. A. J. Chem. Phys. 2003, 118, 7898-7906. (27) Newton, M. D. Coord. Chem. ReV. 2003, 238, 167-185. (28) Ungar, L. W.; Scherer, N. F.; Voth, G. A. J. Phys. Chem. B 1999, 103, 7367-7382. (29) Ungar, L. W.; Newton, M. D.; Voth, G. A. Biophys. J. 1997, 72, 5-17. (30) Zheng, X. H.; Stuchebrukhov, A. A. J. Phys. Chem. B 2003, 107, 9579-9584. (31) Newton, M. D. Int. J. Quantum Chem. 2000, 77, 255-263. (32) Cave, R. J.; Newton, M. D. Chem. Phys. Lett. 1996, 249, 15-19. (33) Kim, J.; Stuchebrukhov, A. A. J. Phys. Chem. B 2000, 104, 86068613. (34) Reimers, J. R.; Hush, N. S. J. Phys. Chem. A 1999, 103, 30663072. (35) Broo, A.; Larsson, S. Chem. Phys. 1992, 161, 363-378. (36) Calzado, C. J.; Sanz, J. F. J. Am. Chem. Soc. 1998, 120, 10511061. (37) Hasegawa, J.; Nakatsuji, H. J. Phys. Chem. B 1998, 102, 1042010430. (38) Regan, J. J.; DiBilio, A. J.; Winkler, J. R.; Richards, J. H.; Gray, H. B. Inorg. Chim. Acta 1998, 276, 470-480. (39) Regan, J. J.; DiBilio, A. J.; Langen, R.; Skov, L. K.; Winkler, J. R.; Gray, H. B.; Onuchic, J. N. Chem. Biol. 1995, 2, 489-496. (40) Skourtis, S. S.; Regan, J. J.; Onuchic, J. N. J. Phys. Chem. 1994, 98, 3379-3388. (41) Gehlen, J. N.; Daizadeh, I.; Stuchebrukhov, A. A.; Marcus, R. A. Inorg. Chim. Acta 1996, 243, 271-282. (42) Gruschus, J. M.; Kuki, A. J. Phys. Chem. B 1999, 103, 1140711414. (43) Baik, M. H.; Crystal, J. B.; Friesner, R. A. Inorg. Chem. 2002, 41, 5926-5927. (44) Adman, E. T.; Jensen, L. H. Isr. J. Chem. 1981, 21, 8-12. (45) Chang, I. J.; Gray, H. B.; Winkler, J. R. J. Am. Chem. Soc. 1991, 113, 7058. (46) Langen, R.; Chang, I. J.; Germanas, J. P.; Richards, J. H.; Winkler, J. R.; Gray, H. B. Science 1995, 268, 1733-1735. (47) Crane, B. R.; Di Bilio, A. J.; Winkler, J. R.; Gray, H. B. J. Am. Chem. Soc. 2001, 123, 11623-11631. (48) Regan, J. J.; DiBilio, A. J.; Langen, R.; Skov, L. K.; Winkler, J. R.; Gray, H. B.; Onuchic, J. N. Chem. Biol. 1995, 2, 489-496. (49) Faham, S.; Day, M. W.; Connick, W. B.; Crane, B. R.; Di Bilio, A. J.; Schaefer, W. P.; Rees, D. C.; Gray, H. B. Acta Crystallogr. D: Biol. Crystallogr. 1999, 55, 379-385. (50) HARLEM, Molecular Modeling Package, Kurnikov, I.V. 2002, http://www.kurnikov.org/harlem download. (51) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman 2002, AMBER 7, University of California, San Francisco. (52) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. S.; Cheatham, T. E., III; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Comput. Phys. Commun. 1995, 91, 1-41. (53) Kuki, A.; Wolynes, P. G. Science 1987, 236, 1647. (54) Stuchebrukhov, A.; Marcus, R. A. J. Phys. Chem. 1995, 99, 7581. (55) Onuchic, J. N.; Beratan, D. N.; Hopfield, J. J. J. Phys. Chem. 1986, 90, 3707-3721. (56) Jordan K. D.; Paddon-Row: M. N. J. Phys. Chem. 1992, 96, 11881196. (57) Liang C. X.; Newton, M. D. J. Phys. Chem. 1992, 96, 2855-2866. (58) 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.;
Electron-Transfer Rates in Metalloproteins Nanayakkara, A.; Challa- combe, 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.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. GAUSSIAN 94, Revision C.2; Gaussian, Inc.: Pittsburgh, PA, 1995. (59) Beratan, D. N.; Hopfield, J. J. J. Am. Chem. Soc. 1984, 106, 15841894. (60) Cave, R. J.; Baxter, D. V.; Goddard, W. A., III; Baldeschwieler, J. D. J. Chem. Phys. 1987, 87, 926-935. (61) Newton, M. D. Chem. ReV. 1991, 91, 767-792. (62) Rodriguez-Monge, L.; Larsson, S. J. Phys. Chem. 1996, 100, 62986303. (63) Ivashin, N.; Kallebring, B.; Larsson, S.; Hansson, O. J. Phys. Chem. B 1998, 102, 5017-5022. (64) Portis, A. M. Electromagnetic Fields; Wiley: New York, 1978. (65) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (66) Liu, Y. P.; Newton, M. D. J. Phys. Chem. 1995, 99, 12382-12386. (67) Kurnikov, I. V.; Zusman, L. D.; Kurnikova, M. G.; Farid, R. S.; Beratan, D. N. J. Am. Chem. Soc. 1997, 119, 5690-5700. (68) Sharp, K. A. Biophys. J. 1998, 74, 1241-1250. (69) Kawatsu, T.; Kakitani, T.; Yamato, T. J. Phys. Chem. B 2002, 106, 11356-11366.
J. Phys. Chem. B, Vol. 109, No. 4, 2005 1625 (70) Lin, J.; Beratan, D. N. J. Phys. Chem. B 2004, 108, 5655-5661. (71) Liang, Z. X.; Kurnikov, I. V.; Nocek, J. M.; Mauk, A. G.; Beratan, D. N.; Hoffman, B. M. J. Am. Chem. Soc. 2004, 126, 2785-2798. (72) Miller, N. E.; Wander, M. C.; Cave, R. J. J. Phys. Chem. A 1999, 103, 1084-1093. (73) Ponce, A.; Gray, H. B.; Winkler, J. R. J. Am. Chem. Soc. 2000, 122, 8187-8191. (74) Dreuw, A.; Head-Gordon, M. J. Am. Chem. Soc. 2004, 126, 40074016. (75) Bacon, A. D.; Zerner, M. C. Theor. Chim. Acta 1979, 53, 21-54. (76) Khoshtaria, D. E.; Wei, J.; Liu, H.; Yue H.; Waldeck, D. H. J. Am. Chem. Soc. 2003, 125, 7704-7714. (77) van Leeuwen, R.; Baerends, E. J. Phys. ReV. A 1994, 49, 24212431. (78) Wang, J.; Stuchebrukhov, A. Int. J. Quantum Chem. 2000, 80, 591597. (79) Reimers, J. R.; Cai, Z. L.; Bilic, A.; Nush, N. S. Ann. N.Y. Acad. Sci. 2003, 1006, 235-251. (80) Kim, K.; Jordan, K. D.; Paddon-Row, M. N. J. Phys. Chem. 1994, 98, 11053-11058.
