Water-Mediated Electron Transfer between Protein Redox Centers

obtained through the Empirical Pathways (EP) method. ..... In spite of its age, the EP method is still ..... (13) Miyashita, O.; Okamura, M. Y.; Onuch...
0 downloads 0 Views 191KB Size
Download PDF
3774

J. Phys. Chem. B 2007, 111, 3774-3781

Water-Mediated Electron Transfer between Protein Redox Centers Agostino Migliore,*,†,‡ Stefano Corni,*,† Rosa Di Felice,† and Elisa Molinari†,‡ National Center on nanoStructures and bioSystems at Surfaces (S3) of INFM-CNR, Modena, Italy, and Dipartimento di Fisica, UniVersita` di Modena e Reggio Emilia, Via Campi 213/A, Modena, Italy ReceiVed: December 20, 2006; In Final Form: January 30, 2007

Recent experimental and theoretical investigations show that water molecules between or near redox partners can significantly affect their electron-transfer (ET) properties. Here we study the effects of intervening water molecules on the electron self-exchange reaction of azurin (Az), by performing a conformational sampling on the water medium and by using a newly developed ab initio method to calculate transfer integrals between molecular redox sites. We show that the insertion of water molecules at the interface between the copper active sites of Az dimers slightly increases the overall ET rate, while some favorable water conformations can considerably enhance the ET kinetics. These features are traced back to the interplay of two competing factors: the electrostatic interaction between the water and protein subsystems (mainly opposing the ET process for the water arrangements drawn from MD simulations) and the effectiveness of water in mediating ET coupling pathways. Such an interplay provides a physical basis for the found absence of correlation between the electronic couplings derived through ab initio electronic structure calculations and the related quantities obtained through the Empirical Pathways (EP) method. In fact, the latter does not account for electrostatic effects on the transfer integrals. Thus, we conclude that the water-mediated electron tunneling is not controlled by the geometry of a single physical pathway. We discuss the results in terms of the interplay between different ET pathways controlled by the conformational changes of one of the water molecules via its electrostatic influence. Finally, we examine the dynamical effects of the interfacial water and check the validity of the Condon approximation.

Introduction The kinetics of biomolecular ET reactions represents a main issue in the current nanoscale research, not only for its central role in many biological processes,1,2 but also because modern electronics aims at exploiting the intrinsic functions of biomolecules to implement nanoelectronic devices.3-5 The accurate prediction of ET reaction rates is a long-standing challenge, especially for long-distance tunneling, which is the major ET mechanism in proteins.6,7 Water is the nearly ubiquitous molecular environment for biological electron transfer. Water interactions with protein systems can play a crucial role in determining the activation free energies and the most effective ET pathways: for instance, water appears to contribute to the best ET pathway between the copper centers in the enzyme PHM.8 However, the efficiency of water in mediating ET reactions is still largely debated.8-13 Even the existence of specific electron tunneling pathways (not directly observed) is under discussion,6,14 although supported by recent theoretical analyses.15,16 In intermolecular ET, the solvent can be the main bridge or a partial bridge17 between the redox groups, lowering the barrier for electron tunneling with respect to vacuum. Several recent investigations17-20 have been focused on the importance of the bridge conformational changes in determining the ET kinetics. In the present work we solvate the azurin dimer and include * To whom correspondence should be addressed. E-mail: smigliore@ unimore.it (A.M.); [email protected] (S.C.). † National Center on nanoStructures and bioSystems at Surfaces (S3) of INFM-CNR. ‡ Universita ` di Modena e Reggio Emilia.

two water molecules experimentally observed in the crystal phase11 at the interface between the copper active sites (Figure 1). We then perform a conformational sampling on the water medium at fixed protein coordinates, for isolating and estimating the ET role of water. We focus our attention on the effects of the interfacial water molecules, introducing a theoretical framework that can give insight into the conformational dependence of water-mediated ET reactions. In particular, we show how the same water molecules between the redox sites of Az dimers can aid or oppose the interprotein ET, depending on the interplay between the water-protein electrostatic interaction and the water quantum bridging (or pathway mediation) for each given conformation. The main computational efforts concern the ab initio calculation of the electron-transfer matrix elements or transfer integrals21 for the electron self-exchange reaction

[Az(I)-Az(II)]solv f [Az(II)-Az(I)]solv where Az(I) [Az(II)] denotes the reduced (oxidized) azurin and [Az(I)-Az(II)]solv refers to the solvated protein dimer. Transfer integrals provide a compact link between the electronic structure of the relevant donor-acceptor interactions and the overall ET rate.22,23 Within the context of Marcus’ ET theory,24,25 they can be easily combined with estimates of the reorganization energy to evaluate the electron-transfer rates, measured in kinetic experiments. Methods MD Simulations. The structure of the Pseudomonas aeruginosa azurin dimer, including the oxygens of the interfacial water, was drawn from the PDB file 1JVL and reasonably mimics the

10.1021/jp068773i CCC: $37.00 © 2007 American Chemical Society Published on Web 03/21/2007

Water-Mediated ET between Proteins

Figure 1. Atomic model of the Az dimer used in the ab initio calculations. C, N, O, S, and H atoms are shown in green, blue, red, yellow, and white spheres, respectively. Cu ions are shown as pink spheres. D and A indicate the donor group and the acceptor group, respectively. All five Cu-ligand amino acids are included, except for a pruning along the backbone.

geometry of the short-living Az complex in solution. In fact, as shown in ref 11, the dimer was obtained from Asn42Cys mutated azurins, by covalently linking the Cys 42 sulfurs via a bis-maleimidomethylether symmetric linker. Note that the mutation is far from the protein regions essentially involved in intradimer ET, to which our electronic structure calculations are limited. In MD simulations, charges +1e and +2e were assigned to Cu(I) and Cu(II), respectively,26 and the relevant force-field parameters were taken from ref 27. The structure was solvated in bulk TIP3P28 water, by using the tleap program29 of AMBER 7.30 Na+ ions with standard force-field parameters31 were added to neutralize the total charge of the system. The resulting simulation box included 12 028 water molecules and had dimensions of ca. 63 × 70 × 87 Å3, in such a way that the closest distance between its walls and the protein atoms was 10 Å. Simulations were performed through the NAMD program,32 with the AMBER force field,33 periodic boundary conditions, a 12.0 Å cutoff distance for van der Waals interactions, and full electrostatics particle-mesh Ewald calculations34 (the PME grid spacing was less than 1 Å). The system structure (with fixed protein coordinates) underwent energy minimization for 2000 steps (by using the conjugate gradient method; see ref 35), followed by 30 ps equilibration dynamics at a temperature of 298 K without pressure control, a further 50 ps equilibration, and a 0.5 ns production MD at the same temperature and pressure p ) 1 bar, thus simulating an isothermal-isobaric (NpT) ensemble (note that the duration of the production MD is similar to those for systems of comparable size in recent works10). Protein coordinates were held fixed in all MD runs. Indeed, in this work we do not consider the thermal protein motion, because we wish to distinguish and quantify more clearly water effects. The temperature was controlled by means of a Langevin thermostat,36 with a damping coefficient of 1 ps-1. The barostat was simulated by a Langevin piston Nose-Hoover method,37-39 with an oscillation period of 200 fs, a damping time of 100 fs, and a noise temperature of 298 K. The achievement of the desired thermodynamic equilibrium conditions was controlled at each step. Furthermore, for the production MD, the root-mean-square deviations of the interfacial water configurations with respect to the average conformation were checked to level off (Figure 2). The MD time step was 0.5 fs, while coordinates were collected every 2 ps, thus obtaining 250 system configurations with relative water conformations.

J. Phys. Chem. B, Vol. 111, No. 14, 2007 3775

Figure 2. Root-mean-square deviations (rmsd) of the interfacial water configurations with respect to the average configuration.

Transfer Integral Calculations. Many ET reactions in protein systems are characterized by a weak electronic coupling between the donor and acceptor groups; thereby, the reaction rates are approximately given by the high-temperature, nonadiabatic expression17,25

kET )

[

]

2 (∆G0 + λ)2 π 〈UIF 〉 exp λkBT p 4λkBT

x

(1)

where λ is the reorganization energy, ∆G0, is the reaction free energy, kB is Boltzmann’s constant, T is the temperature, and 〈U2IF〉 is the mean-square value of the transfer integral (or electron-transfer matrix element) that measures the coupling between the initial state I and the final state F of the ET reaction. A preliminary estimate of the transfer integral for each sampled configuration was obtained through a best-pathway calculation, by using the HARLEM program.40 Hence, the conformations with suitably similar transfer integrals and structures were grouped, a conformation was selected for each group, and the inherent transfer integral was derived from electronic structure quantities computed through ab initio calculations. Electronic structure calculations were limited to the atomic protein model of Figure 1, with the different conformations of the interfacial water molecules. Indeed, due to the tunneling nature of the electronic coupling between protein partners, many ET reactions involve relatively small protein regions around the redox sites.41 Transfer integrals were derived from the formula42

|UIF| )

|

∆EIFab

a2 - b2

|

(2)

where ∆EIF is the energy difference between the ET diabatic states I and F, and a and b are their respective overlaps with the ground state of the system. The approximations done to derive eq 2 are discussed in detail in ref 42, where transfer integrals obtained from our method were compared to those from an energy splitting method43 and turned out also to be in general agreement with the empirical average-packing-density model.6 Within the theoretical framework of ref 42, eq 2 provides the best estimate of the transfer integral also when the two-state approximation does not hold. All the involved quantities were obtained through ab initio computations, in a density-functional theory (DFT) scheme, which is the best compromise between accuracy and computational feasibility in studying large metal-ion complexes. A discussion on the application of DFT to calculate transfer integrals is presented in ref 42. The energy difference ∆EIF is given by the equation

3776 J. Phys. Chem. B, Vol. 111, No. 14, 2007

∆EIF ) (ED + EA) - (ED+ + EA-) + UD-A - UD+- A-

Migliore et al.

(3)

where ED (ED+) is the ground-state energy of the isolated group D (see Figure 1) in the reduced (oxidized) state of charge, EA (EA-) is the same for group A in the oxidized (reduced) state, and UD-A and UD+-A- are the energies of (essentially electrostatic) interaction between the donor and acceptor groups in the initial and final diabatic states, respectively. As usual for long-distance through-space ET reactions between chargelocalized states,44 in implementing eq 2 we distinguished a donor group D on which the transferring electron is initially localized (that is the initially reduced azurin active site), and an acceptor group A, receiving the transferred electron. Thus, we defined the initial vector state of the system as |ψI〉 ) |D〉|A〉 and the final vector state as |ψF〉 ) |D+ 〉|A- 〉, where |D〉 (|D+〉) denotes the reduced (oxidized) ground state of the isolated donor redox site, and |A- 〉 (|A〉) the reduced (oxidized) ground state of the isolated acceptor redox site. In defining the D and A groups, each interfacial water molecule was added to the adjacent protein site. The water medium surrounding the azurin molecules was explicitly included in the MD simulations, because it could influence, to some extent, the geometry of the solvent inside the cavity, as a consequence of solvent-solvent interactions essentially mediated by the protein matrix. However, the external water was not considered in the ab initio computation of U2IF, in that it is far from the region between the redox sites and thus it is not directly involved in the relevant ET processes. Protein Atomic Model. Figure 1 shows the atomic model used for the ab initio electronic structure calculations. The Cuto-Cu distance is 14.6 Å. In each Az redox site the copper ion is surrounded by five ligands in a distorted trigonal bipyramidal geometry. Cys112, His46, and His117 are the strong equatorial ligands, while Met121 and Gly45 are weakly bound axial ligands. The two His117 moieties play an important role in the interprotein ET. Indeed, such residues protrude into the cavity between the redox sites, thus enlarging the interprotein throughbond region and significantly lowering the tunneling barrier for electron transfer between azurins.42 Moreover, the two His117 conjugated imidazole rings are expected to efficiently support charge transfer, because of their delocalized electronic character. Computational Details. The pruning that yielded the model in Figure 1 was validated elsewhere.42 H atoms were added to saturate the dangling bonds and the O atoms of the interfacial water molecules. Their relaxation and the computation of nearly all the electronic quantities involved in eq 2 were carried out through the PWSCF code,45 in the repeated supercell approach, by using the DFT generalized gradient approximation (GGA), with a PW91 exchange-correlation functional, a plane-wave basis set with an energy cutoff of 25 Ry, and ultrasoft pseudopotentials46 for representing the atomic cores. The repeating cell had dimensions of 26 × 27 × 33 Å3. This choice ensured that the minimum vacuum thickness between two replicas was above 15 Å and the minimum interdimer Cu-toCu distance was about 23 Å, which is significantly larger than the intradimer value of 14.6 Å. A jellium background balanced the positive charge of the oxidized site, thus establishing the electrical neutrality of the cell. The relevant wave functions were derived as single Slater determinants of the lowest lying occupied KS spin-orbitals. It was discussed elsewhere42 that spin contamination47 is negligible for the system under study. First the ground-state wave function of the whole system, as well as the ground-state energies and wave functions of the isolated D and A groups (in both oxidation states) were obtained by means of the PWSCF code; then the

Figure 3. Time evolution of T2IF along the MD trajectory. Snapshots were taken every 2 ps. Note that not all points corresponding to higher couplings are isolated. Squares denote the snapshots selected for ab initio calculations. The selection criterion is explained in the text.

overlap integrals a and b were derived by the DTI program.48 UD-A and UD+-A- were computed as electrostatic interactions, by using the NWChem package,49 with an Ahlrichs-pVDZ basis set for all atomic species.50 Quadrupole-quadrupole and higher order terms were safely neglected. Results Dependence of the Best-Pathway Couplings on Water Conformation. The EP model has been recently employed10,13 in various electronic coupling calculations, also in the presence of the water medium. In spite of its age, the EP method is still trustfully employed either for describing situations in which the single-pathway picture suffices or as a useful starting point from which to build up multipathways analyses.16,51 Thus, we believe that, despite the simplicity and poor refinement of the EP model with respect to ab initio approaches, the EP results constitute a natural benchmark for our ab initio analysis. According to the EP model, the electronic matrix element is proportional to the product (“pathway product”,16 here named TIF) of the decay factors for each step in the dominant pathway tube connecting the donor and acceptor groups.52 Each decay factor describes the exponential decay of the electronic coupling along a step, with a decay constant depending on the throughbond or through-space character of the step. The decay constant weighs the length of the step, resulting in an associated (unitless) effective distance. In this picture, TIF corresponds to the shortest overall effective distance between donor and acceptor.52,53 Figure 3 shows the evolution of T2IF (the quantity proportional to the ET rate) along our MD trajectory. Once the protein configuration was fixed, we obtained the same value of T2IF for all the “bulk” conformations (or snapshots) in which the water molecules are not involved in the best ET path (named “b-conformations” in the following). In other words, the employed empirical method52 is not very sensitive to geometry (as also stated in ref 16) and does not distinguish between different unfavorable water conformations. Thirteen percent of configurations (named “wconformations”) include a water molecule in the best tunneling path, which always passes through the N atoms of the His117’s. However, the corresponding couplings do not considerably affect the mean-square pathway product 〈T2IF〉. In fact, the value 〈T2IF〉 ) (7.06 ( 0.15) × 10-11, obtained by averaging all the snapshots, resulted in a value only ∼9% larger than that excluding w-conformations, that is 6.49 × 10-11.54 Indeed, according to the EP model, only in a few w-conformation water molecules appear to establish a strongly coupling pathway. As compared to what was recently found for the trypsin-solubilized

Water-Mediated ET between Proteins

J. Phys. Chem. B, Vol. 111, No. 14, 2007 3777

Figure 4. Selection of w-conformations (involving a water molecule in the strongest ET path) for the ab initio calculations. The selected snapshots are indicated as squares in both panels. (a) T2IFvs R (i.e., through-space length of the best ET path, modified to include the relevant water hydrogen). (b) T2IFvs Rbp (i.e., through-space length of the best ET path).

TABLE 1: Square Pathway Products and Transfer Integrals w-c w-c w-c w-c w-c w-c w-c b-c w-c b-c Figure 5. Ab initio square transfer integrals (U2IF) vs square pathway products (T2IF).

bovine liver cytochrome b5,10 this points to a “direct contact” regime, where ET can be directly facilitated by protein atoms, thus reducing the relative importance of water-mediated ET. This is supported by the fact that the distance between His117’s is around 7 Å. Snapshot Selection for ab Initio UIF Computation. The existence of different coupling regimes may be a robust characteristic of solvated protein systems, which can be caught by empirical or semiempirical methods (e.g., the extendedHu¨ckel method in ref 10). On the other hand, the accuracy of semiempirical approaches to electronic coupling calculations depends on the parametrization and the inherent approximations, which is a particularly delicate matter when molecular contacts include through-space regions. Moreover, potentially important electrostatic effects are not considered at an extended-Hu¨ckel (or similar) level. Thus, for obtaining a reliable estimate of 〈U2IF〉, we resorted to complete ab initio calculations of square transfer integrals for suitably selected conformations and then averaged them. In fact, ab initio calculations on all MD configurations are presently unfeasible. The selection method was the following. First, w-conformations were grouped as in Figure 4a, where R is the through-space length of the best ET path and each group is identified with a color. R values were computed by including the relevant water hydrogen, whereas the protein through-bond portion common to all the snapshots was excluded, as well as the length of the O-H bond that turned out to be nearly independent of the configuration. Each group in Figure 4 is characterized by R and T2IF ranges within 25% and 20% of the respective global ranges. Then eight wconformations (one for each group; see squares in Figure 4a) were selected for ab initio calculations, as a compromise between their central positions in the inherent R ranges and the requirement of a maximum spreading of the corresponding T2IF values. As shown in Figure 4b, the chosen configurations are also rather evenly distributed along the Rbp coordinate, repre-

time (ps)

T2IF

U2IF (eV2)

10 30 50 54 82 84 110 180 304 380

7.28 × 1.36 × 10-10 6.63 × 10-11 2.89 × 10-10 1.16 × 10-10 1.77 × 10-10 2.11 × 10-10 6.49 × 10-11 1.02 × 10-10 6.49 × 10-11

8.82 × 10-4 3.03 × 10-5 3.89 × 10-7 4.56 × 10-5 5.59 × 10-5 2.69 × 10-5 4.94 × 10-5 5.30 × 10-5 2.52 × 10-5 2.94 × 10-5

without water

8.82 × 10-5

10-11

senting the through-space length of the best ET pathway directly derived from the HARLEM program (which neglects H atoms). Two further snapshots were chosen among b-conformations (at MD times of 180 and 380 ps) for representing the flat regions on the graph of Figure 3, in such a way as to warrant the maximum geometrical spread of the overall set of selected conformations. Ab Initio UIF Calculations: Performance vs EP Model and Comparison with Experiment. Table 1 shows the values of the squared transfer integrals U2IF from ab initio calculations and the corresponding T2IF for the selected conformations and for the protein system without water. Let us first observe that the two b-conformations give transfer integrals similar to those of most w-conformations. Thus, the relative enhancement of T2IF for w-conformations appears insufficient to meet a similar trend of the corresponding ab initio quantities. This is even more evident from Figure 5, where we plot U2IF values (in a log scale, for clarity) against T2IF values for all the selected conformations. Although a quantitative comparison between the T2IF and U2IF sets is not possible, the general mismatch of their relative values shown in Figure 5 is a clear sign of the existence of ET controlling factors beyond the EP model. We remark that this discrepancy cannot be ascribed to the fact that a too small set of MD snapshots was selected for the evaluation, because we indeed opted for a “smart choice”, based on the following facts. (i) The groups of snapshots identified in Figure 4 represent a wide range of transfer integral values, and we have extracted one configuration for each set, thus covering all the possible transfer integral range, but with a minimal set of structures. (ii) From a temporal point of view, the selected conformations are widely spread along the MD run, as a consequence of the time occurrence of peaks shown in Figure 3. (iii) Bulk and wconformations have similar couplings. Thus, the richer choice among the latter ones is not expected to bias the estimation. Taking jointly into account arguments i-iii, we are confident

3778 J. Phys. Chem. B, Vol. 111, No. 14, 2007

Migliore et al.

that the sampled conformations represent a reliable “smart set” for averages that give the best estimate of the mean-square transfer integral. Furthermore, additional analysis on different w-conformations, taking into account the respective U2IF values, can shed light on the interplay of physical factors that determine the ET kinetics beyond the geometrical and parametrical properties of the best ET path. Our best estimate of the square transfer integral is 〈U2IF 〉 ) (1.2 ( 0.8) × 10-4 eV2, to be compared with the value 0.88 × 10-4 eV2 for the protein system without interposed water.42 The averaging error on 〈U2IF〉 is large, as is generally the case for electronic couplings computed by high-level quantum chemical methods. Indeed, electronic structure calculations certainly reduce the systematic errors of average values, but their computational cost limits the number of data points used in the averaging procedure, thus hardly addressing statistical errors.55 It is important to stress that the presence of the water molecules does not give the dramatic increase of transfer integral suggested in ref 11. Note that the value of 〈U2IF〉 turns out to be larger than the square transfer integral in the absence of water due to the conformation at 10 ps, thus indicating that ET kinetics is controlled by particularly suitable far-from-equilibrium conformations of the medium between donor and acceptor. Moreover, possibly, some (here not sampled) rare configurations could further increase the ET rate in the presence of the water dimer. By inserting 〈U2IF〉 in eq 1, with a temperature T of 298 K, a reorganization energy λ of 0.8 eV,56 and ∆G0 ) 0 (as is the case for a self-exchange reaction), we obtain the ET rate kET ) 9 × 108 s-1. In ref 11 a lower limit of 5 × 104 s-1 is established for kET. Indeed, a more cogent comparison with experiment would require the computational analysis of further factors, such as the effects of protein conformational fluctuations. On the other hand, experiments able to yield less uncertain values are desired for a more stringent test of the theory. Dynamical Characterization of Water Fluctuations. When the charge transfer occurs through a fluctuating medium, eq 1 provides just the leading term k(0) ET in a series expansion of the ET rate.17,18 If the conformational changes of the intervening medium significantly affect the electronic coupling, the Condon approximation generally fails57 and further terms must be considered in the expression of the ET rate. Moreover, the dynamics of conformational changes affects ET rates when such changes happen on a sufficiently fast time scale.17,18 For a selfexchange reaction and “classical” motions of the medium, the first nonzero correction term can be expressed as18

( )(

p (0) 1 k(2) ET ) kET 8 τckBT

2

1-

)

2kBT (1 - Rc) λ

(4)

where τc is the correlation time of the electronic coupling and Rc ) 〈UIF〉2/〈U2IF〉 is the coherence parameter, which measures the amplitude of the fluctuations of the transfer integral. For the system under study, τc measures the rapidity of the conformational changes of the bridging water molecules that are significant for the ET efficiency. Changes in the geometrical parameters of the water molecules due to stretching or bending vibrations were estimated to be within 2.5%, thus having negligible effects on the relevant electrostatic interactions and effective distances between donor and acceptor. On the other hand, the fastest displacements of the water molecules are oscillations around temporary equilibrium positions with an average period of ∼2 × 10-13 s,58 which can be assumed as a reasonable estimate of the effective correlation time. Therefore, (0) -3 it is k(2) ET/kET ≈ 10 (1 - Rc). From the smallness of this ratio

we conclude that the Condon approximation works well for our system and the correction term k(2) ET is negligible independently of the value of the coherence parameter, which turned out to be Rc ) 0.53. This relatively high value of Rc means that the water molecules are far from extreme flexibility or free-motion condition: in such conditions, in fact, a value of Rc close to zero is generally expected.18 They are also far from a tightbinding condition, necessary for a significant intervention of high-frequency modes. The above results, together with the validation of the Condon approximation for the single (not fixed) azurin molecule given elsewhere,59 exclude significant nonCondon effects for the fluctuating azurin dimer far from the charge injection limit.60 Interplay of Bridging and Electrostatic Effects. It is clear from Figure 3 that the EP model does not allow for the possibility that interfacial water slows down the ET kinetics. In fact, the water dimer can only increase the transfer integral whenever its conformation reduces the effective distance between donor and acceptor. On the other hand, electronic structure calculations show (see Table 1) that various water configurations can negatively affect the ET rate. This behavior can be ascribed to the simultaneous action of two factors: the bridging role of water and the electrostatic interaction between the water and protein subsystems.61 The two factors are obviously tangled, in that the charge distribution on each water molecule is influenced by its quantum interaction with the remainder. However, the lack of a tight binding between the water and protein subsystems, along with the strong dipolar character of water molecules, allow a first-order separation of the two effects. In particular, we verified, through the RESP algorithm,49 that the partial atomic charges on each water molecule changed within a few percent (j5% ) for all the different conformations. Similar changes did not significantly affect the electrostatic interactions between different atomic groups, although they may have a considerable influence on the interprotein ET. At first, we closely analyze and compare the crystallographic conformation and the ET-dominant conformation, the latter reached after 10 ps of MD (see Figure 6). The water molecule near the acceptor site (aw) is involved in the shortest ET pathway, which passes through the N(His117) atoms. This mol ecule is slightly better oriented in the ET-dominant conformation, so that the best path is shortened by 0.1 Å. On the other hand, the water molecule adjacent to the donor site (dw) moved considerably away, so that the distance between the two water oxygens changed from 2.91 Å to 3.39 Å. In both conformations, the dw molecule is out of any favorable ET pathway, because of its distance from the acceptor His117 and of its orientation. In fact, although the dw oxygen can form a hydrogen bond with , this the hydrogen attached to the nearby donor N(His117)  H-bond is not part of an otherwise favorable network of H-bonds connecting the two histidines. On the other hand, the dw molecule can indirectly affect the ET process in a significant manner, via electrostatic interactions. The dipole moment of the dw unit is directed toward the acceptor group. Also taking into account its orientation with respect to the Cu-to-Cu direction, mainly involved in the electron-transfer process, the latter circumstance is expected to electrostatically oppose the electron charge transfer from the rear donor protein.62 After the displacement of the dw molecule, the change of electrostatic effects considerably promotes the ET reaction in the conformation at 10 ps. A measure of these effects is given by the consequent rearrangement of the electronic charge (quantified as a variation of the Mulliken

Water-Mediated ET between Proteins

J. Phys. Chem. B, Vol. 111, No. 14, 2007 3779

Figure 6. (a) Change of the Mulliken population on the donor (D) site and on the water molecule near the acceptor (A) site under the electrostatic field of the water molecule adjacent to the D site, when the last moves from the crystallographic position (gray color) to the one after 10 ps of MD (purple). The color scale is in units of the absolute electron charge. (b) Enlargement around the gap between the proteins, preserving only the His117 residues. The color code refers only to the donor His117. The yellow dashed line represents the best ET pathway. The role of pathways 1 and 2 is discussed in the text. (c) Two additional water conformations are represented for comparison: at 50 ps (yellow) and 138 ps (pink) along the MD trajectory.

population63 in Figure 6) on aw, important for the quantum bridging, and on the donor protein, where the transferring electron is initially localized. First of all, as shown in Figure 6a, the increased distance between the two water oxygens allows a larger electron population on the aw hydrogen (and, slightly, on the aw oxygen) directly involved in the best ET pathway (yellow dashed line in Figure 6b, including the relevant hydrogens). The resulting effect along this path is an electron charge transfer from the donor His117 to the halfway water molecule. Another effect promoting ET comes from the charge redistribution on the aromatic ring of the donor His117 (see Figure 6b). On the whole, there is a transfer of electron charge to the donor Cδ-Hδ group, which is far from the dw oxygen and protrudes toward the acceptor group. On the other hand, in the conformation at 10 ps, the minimum distance between the donor Hδ atom and the acceptor His117 is only 0.13 Å longer than the through-space length of the best ET pathway. Consequently, the increased electron population on the Cδ-Hδ group can appreciably contribute to the interprotein ET. We remark that the rearrangement of electronic charge is relatively small, but significant for a long-range ET reaction with a considerable through-space character. In a multiple ET pathways view,10,64 with reference to Figure 6b, we can say that the change in the water conformation reduced the unfavorable interplay (from an electrostatic, not a quantum interference, point of view) between the best pathway and pathway 1 and enforced the electron transfer through pathway 2. In Figure 6c we extend the comparison to other water conformations. After 50 ps (yellow color), the dw molecule is involved in the shortest ET path (with a through-space length Rbp ) 7.0 Å, to be compared with Rbp ) 6.5 Å after 10 ps), but with an unfavorable orientation. In fact, its dipole moment tends to repel the transferring electron back toward the donor group. Moreover, aw is far apart from the His117-to-His117 direction. This concomitance of unfavorable circumstances agrees with the particularly small transfer integral obtained by ab initio calculations. A further validation of the correspondence between the above quantum-electrostatic arguments and the computed transfer

integrals came from the search for a w-conformation similar to the one after 10 ps. We found that, in the conformation at 138 ps (pink color in Figure 6c), the dw unit is halfway between the crystallographic position and the conformation at 10 ps. The reciprocal arrangement of the dw and aw molecules is similar to that after 10 ps. Indeed, electronic structure calculations yielded U2IF ) 1.16 × 10-4 eV2, which is the second largest U2IF value, after that at 10 ps. Obviously, a similarly detailed analysis on each snapshot would be required to find otherwise favorable dispositions, also among the b-conformations. Conclusions In this paper we combined several theoretical and computational tools to disentangle and understand the different involved aspects of the ET mediation role of water in Az complexes. After a conformational sampling in a suitably large water box, the analysis was focused on the effects of a water dimer that can fit the spatial gap between two Az redox sites. It was previously argued11,65 that this water dimer, by virtue of its crucial position, could largely increase the interprotein ET rate. Here, instead, both the EP method and the ab initio multielectron calculations point to a slightly enhanced ET kinetics, dominated by quite a few strongly favorable water arrangements. Moreover, contrary to the expectations of the EP method, our quantumelectrostatic calculations illustrate how the ET rate also can be instantaneously decreased. This happens whenever the water configuration does not provide an efficient quantum bridging, while concurring to determine an electron distribution on the protein matrix that is unfavorable to the electron-transfer process. The loose contact between the water and protein subsystems allows very different water configurations, so that the conformational sampling turned out to be essential for a reliable estimate of the water effects. On the other hand, the relevant water dynamics did not appear fast enough to appreciably affect the ET rate. In other terms, the system is in the limit of slow conformational modulation, where dynamical effects can be ignored,18 and it is sufficient to replace the square transfer integral with its mean value (with respect to water configurations) in eq 1.

3780 J. Phys. Chem. B, Vol. 111, No. 14, 2007 Our analysis indicates a breakdown of the single-pathway picture51 for the ET system under study,66 while it is consistent with the multiple ET pathways picture, recently proposed10,64 for other systems, and is able to read the interplay between ET pathways in terms of quantum and electrostatic effects. The EP method allowed an approximate quantification of exclusively quantum effects strictly related to the water arrangement between the copper sites, while the Mulliken analysis provided evidence of the ET effects of water-protein electrostatic interactions. Both kinds of effects are consistently taken into account by our method to calculate transfer integrals (eq 2). The effects of the water-protein electrostatic interactions on the transfer integrals, described in this paper, are to be distinguished from the electrostatic effects on the reorganization energy. Indeed, the former affect the probability of an ET reaction once the transition state is reached, while the latter affect the probability to reach such a transition state. Note also that our method (as is generally the case for electronic-structurebased methods16) includes, by construction, the effects of the quantum interference between multiple tunneling pathways. We believe that the methodology used in this work, combining pathways and quantum-chemical analyses, can be a fruitful approach for the predictive investigation of ET reactions in biochemical environments. Acknowledgment. We thank Manuela Cavallari, Andrea Ferretti, Carlo Cavazzoni, and Arrigo Calzolari for fruitful discussions and Anna Garbesi for carefully reading the manuscript. Funding was provided by INFM-CNR through the allocation of computer time at the CINECA facilities (Bologna, Italy), by MIUR-IT through project FIRB-NOMADE, by the EC through contract IST-2001-38951. Supporting Information Available: Analysis of the transfer integrals computed with a slightly different implementation of the EP method, according to ref 53: two figures and discussion. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Berg, J. M.; Stryer, L.; Tymoczko, J. L. Biochemistry, 5th ed.; Freeman: New York, 2002. (2) Osyczka, A.; Moser, C. C.; Daldal, F.; Dutton, P. L. Nature 2004, 427, 607-612. (3) Nitzan, A. Annu. ReV. Phys. Chem. 2001, 52, 681-750. (4) Aviram, A.; Ratner, M.; Mujica, V. Molecular Electronics II. Annals of the New York Academy of Sciences; The New York Academy of Sciences: New York, 2003; Vol. 960. (5) Maruccio, G.; Biasco, A.; Visconti, P.; Bramanti, A.; Pompa, P. P.; Calabi, F.; Cingolani, R.; Rinaldi, R.; Corni, S.; Di Felice, R.; Molinari, E. AdV. Mater. (Weinheim, Ger.) 2005, 17, 816-822. (6) Page, C. C.; Moser, C. C.; Chen, X.; Dutton, P. L. Nature 1999, 402, 47-52. (7) Gray, H. B.; Winkler, J. R. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 3534-3539. (8) Francisco, W. A.; Wille, G.; Smith, A. J.; Merkler, D. J.; Klinman, J. P. J. Am. Chem. Soc. 2004, 126, 13168-13169. (9) Ohta, K.; Okada, O.; Yoshida, T. J. Phys. Chem. 1989, 93, 932936. (10) Lin, J.; Balabin, I. A.; Beratan, D. N. Science 2005, 310, 13111313. (11) van Amsterdam, I. M. C.; Ubbink, M.; Einsle, O.; Messerschmidt, A.; Merli, A.; Cavazzini, D.; Rossi, G. L.; Canters, G. W. Nat. Struct. Biol. 2002, 9, 48-52. (12) Wenger, O. S.; Leigh, B. S.; Villahermosa, R. M.; Gray, H. B.; Winkler, J. R. Science 2005, 307, 99-102. (13) Miyashita, O.; Okamura, M. Y.; Onuchic, J. N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 3558-3563. (14) Langen, R.; Chang, I.; Germanas, J. P.; Richards, J. H.; Winkler, J. R.; Gray, H. B. Science 1995, 268, 1733-1735.

Migliore et al. (15) Regan, J. J.; Onuchic, J. N. AdV. Chem. Phys. 1999, 107, 497553. (16) Jones, M. L.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. A 2002, 106, 2002-2006. (17) Troisi, A.; Ratner, M. A.; Zimmt, M. B. J. Am. Chem. Soc. 2004, 126, 2215-2224. (18) Troisi, A.; Nitzan, A.; Ratner, M. A. J. Chem. Phys. 2003, 119, 5782-5788. (19) Barnett, R. N.; Cleveland, C. L.; Joy, A.; Landman, U.; Schuster, G. B. Science 2001, 294, 567-571. (20) Ratner, M. A. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 387-389. (21) Newton, M. D. Chem. ReV. 1991, 91, 767-792. (22) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265322. (23) Farazdel, A.; Dupuis, M.; Clementi, E.; Aviram, A. J. Am. Chem. Soc. 1990, 112, 4206-4214. (24) Bixon, M.; Jortner, J. AdV. Chem. Phys. 1999, 106, 35-202. (25) Kuznetsov, A. M.; Ulstrup, J. Electron Transfer in Chemistry and Biology; John Wiley & Sons: New York, 1999. (26) Wang, C. X.; Bizzarri, A. R.; Xu, Y. W.; Cannistraro, S. Chem. Phys. 1994, 183, 155-166. (27) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024-10035. (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (29) Schafmeister, C. E. A. F.; Ross, W. S.; Romanovski, V. LEAP; University of California: San Francisco, CA, 1995. (30) 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.; et al. AMBER 7; University of California: San Francisco, CA, 2002. (31) Aqvist, J. J. Phys. Chem. 1990, 94, 8021-8024. (32) Kale´, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. J. Comput. Phys. 1999, 151, 283-312. (33) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (34) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577-8593. (35) Fletcher, R.; Reeves, C. M. Comput. J. 1964, 7, 149-154. (36) Kolb, A.; Du¨nweg, B. J. Chem. Phys. 1999, 111, 4453-4459. (37) Bhandarkar, M.; et al.: http://www.ks.uiuc.edu/Research/namd/2.5/ ug/. (38) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177-4189. (39) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613-4621. (40) Kurnikov, I. V. HARLEM: http://www.kurnikov.org/. (41) Stuchebrukhov, A. A. Theor. Chem. Acc. 2003, 110, 291-306. (42) Migliore, A.; Corni, S.; Di Felice, R.; Molinari, E. J. Chem. Phys. 2006, 124, 64501. (43) Prytkova, T. R.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. B 2005, 109, 1618-1625. (44) Newton, M. D. Theor. Chem. Acc. 2003, 110, 307-321. (45) Baroni, S.; Dal Corso, A.; de Gironcoli, S.; Giannozzi, P.; Cavazzoni, C.; Ballabio, G.; Scandolo, S.; Chiarotti, G.; Focher, P.; Pasquarello, A.; Laasonen, K.; Trave, A.; Car, R.; Marzari, N.; Kokalj, A.: http://www.pwscf.org/. (46) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892-7895. (47) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory; Wiley: New York, 2000. (48) (a) Ferretti, A.; Ruini, A.; Bussi, G.; Molinari, E.; Caldas, M. J. Phys. ReV. B 2004, 69, 205205. (b) Ferretti, A. DTI program, 2005: http:// www.s3.infm.it/dti. (49) (a) Apra`, E.; Windus, T. L.; Straatsma, T. P.; Bylaska, E. J.; de Jong, W.; Hirata, S.; Valiev, M.; Hackler, M.; Pollack, L.; Kowalski, K.; et al. NWChem, A Computational Chemistry Package for Parallel Computers, Version 4.7; Pacific Northwest National Laboratori: Richland, WA 99352-0999, 2005. (b) Kendall, R. A.; Apra`, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J.; Nichols, J. A.; Nieplocha, J.; et al. Comput. Phys. Commun. 2000, 128, 260-283. (50) The basis set convergence has been tested. In particular, we found that the Ahlrichs-pVDZ, aug-cc-pVDZ, and aug-cc-pVTZ basis sets agreed within 5% in computing the dipole moments of the water molecules. (51) A detailed and general discussion on the breakdown of the singlepathway picture in cases where the multiple pathways become important can be found in: Regan, J. J.; Risser, S. M.; Beratan, D. N.; Onuchic, J. N. J. Phys. Chem. 1993, 97, 13083-13088. (52) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488-4498. (53) Betts, J. N.; Beratan, D. N.; Onuchic, J. N. J. Am. Chem. Soc. 1992, 114, 4043-4046.

Water-Mediated ET between Proteins (54) Compatibly with the relatively small number of w-conformations, we verified that the T2IF values for the w-conformations approximately follow a Gaussian, one tail distribution, with 〈T2IF〉w ) (1.2 ( 0.4) × 10-10. Combining this result with the EP coupling for b-conformations (equal to that in the absence of water, since such conformations do not comprise any water molecule in the best ET pathway), we obtained 〈T2IF〉 ) (7.2 ( 0.5) × 10-11. Therefore, also this more sophisticated statistical approach points to an increase of the EP coupling in the presence of the interfacial water. (55) Let us note that the U2IF evaluated for different conformations, because of their spreading on the MD time scale, can be considered as “independent estimates” of 〈U2IF〉. Hence, the error associated with 〈U2IF〉 can be read as an estimate of the standard deviation of the corresponding normal distribution. This gives us an ∼66% probability that the presence of interfacial water molecules increases the square electronic coupling between donor and acceptor, and thus the ET kinetics. (56) Di Bilio, A. J.; Hill, M. G.; Bonander, N.; Karlsson, B. G.; Villahermosa, R. M.; Malmstro¨m, B. G.; Winkler, J. R.; Gray, H. B. J. Am. Chem. Soc. 1997, 119, 9921-9922. (57) Davis, W. B.; Ratner, M. A.; Wasielewski, M. R. J. Am. Chem. Soc. 2001, 123, 7877-7886.

J. Phys. Chem. B, Vol. 111, No. 14, 2007 3781 (58) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: Oxford, UK, 1969. (59) Skourtis, S. S.; Balabin, I. A.; Kawatsu, T.; Beratan, D. N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 3552-3557. (60) Skourtis, S. S. Chem. Phys. Lett. 2003, 372, 224-231. (61) Migliore, A.; Corni, S.; Di Felice, R.; Molinari, R. J. Phys. Chem. B 2006, 110, 23796-23800. (62) On the other hand, the analysis in terms of the shape of the overall electrostatic barrier is (also theoretically) complicated by the fact that we are considering a three-dimensional situation. (63) The electrostatic field of the water molecule near the donor Az was reproduced by point RESP charges. These were derived through the NWChem Electrostatic Potential (ESP) module. (64) Balabin, I. A.; Onuchic, J. N. Science 2000, 290, 114-117. (65) Mikkelsen, K. V.; Skov, L. K.; Nar, H.; Farver, O. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 5443-5445. (66) The data presented in the Supporting Information, derived from an EP method that takes explicitly into account the hydrogen bonds,53 show the robustness of such a conclusion with respect to the details of the employed EP method.