Calculation of Redox Properties: Understanding Short- and Long

Basu, G.; Kitao, A.; Kuki, A.; Go, N. J. Phys. Chem. B 1998, 102, 2076−2084. [ACS Full Text ACS Full Text ], [CAS]. (46) . Protein Electron Transfer...
0 downloads 0 Views 381KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 3969-3976

3969

Calculation of Redox Properties: Understanding Short- and Long-Range Effects in Rubredoxin Marialore Sulpizi,*,† Simone Raugei,‡ Joost VandeVondele,§ Paolo Carloni,‡ and Michiel Sprik† Department of Chemistry, UniVersity of Cambridge, Lensfield Road, CB2 1EW Cambridge, United Kingdom, SISSA, Via Beirut 2, 34014 Trieste, Italy, and Institute of Physical Chemistry, UniVersity of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland ReceiVed: NoVember 8, 2006; In Final Form: February 20, 2007

In this computational study we show that for rubredoxin, a small and comparatively simple iron-sulfur protein, it is possible to combine a full ab initio description of the electronic structure of the protein in explicit solvent with sampling of the relevant time scale of the protein dynamics by using a hybrid method based on a force field molecular dynamics/density functional theory scheme. Applying this scheme within the framework of Marcus theory we are able to reproduce the experimental redox potential difference of 60 mV between a mesophilic and thermophilic rubredoxin within an accuracy of 20 mV and explain it in terms of short-range contributions from a few residues close to the metal center. We also compute the reorganization free energy for oxidation of the protein obtaining 720 meV for the mesophilic and 590 meV for thermophilic variant. Decomposition of the reorganization energy by using the classical force field shows that this quantity is largely determined by the solvent, with both short-range (an oxidation induced change of coordination number) and long-range (dielectric) contributions. The 130 meV higher value for the mesophilic form is analyzed in terms of detailed differences in the solvent structure around the metal center and the dielectric response. These results underline the importance of a molecular description of the solvent and of a correct inclusion of the polarization effects.

I. Introduction Bioinorganic oxidation/reduction enzymes and metalloproteins represent more than 40% of IUBMB classified proteins and are not only vital to biological energy conversion in photosynthesis and respiration, but are also critical to a growing number of signaling processes governing gene regulation and expression.1 Understanding the mechanism of electron transfer (ET) between two metal sites or metal site and organic substrate is therefore of both theoretical and practical importance. The key question is how metalloproteins control which processes are thermodynamically feasible (i.e., reduction potentials) and how fast they occur (i.e., rate constants). Some of the issues that have been raised in this context concern the competition between short-range effects, in essence the coordination chemistry of the metal ion, and long-range effects, for example, the reorganization of the protein, the placement of charged and polar residues, and the access to the solvent. Related to this is the question of the relative importance of electronic relaxation effects, such as the difference between hard and soft ligands and electronic polarization of the protein, versus the reorganization due to atomic motion. ET in biological systems has attracted an enormous interest and effort and many of the important questions have been at least partly answered.2 Marcus theory3 has provided a powerful conceptual and theoretical framework for these developments. Computation has also contributed, mainly in the interpretation of experiment and by underpinning specific explanations with * To whom correspondence should be addressed. Phone: +44-1223336423. Fax: +44-1223-336362. E-mail: [email protected]. † University of Cambridge. ‡ SISSA. § University of Zurich.

numerical evidence. Investigation of questions such as competition between short- and long-range reorganization requires large systems and simulation over long periods of time.4 Moreover, for bioorganic systems the detail of electronic structure often matters and more quantitative comparisons can only be performed by using ab initio electronic structure calculation methods. These two computational requirements are usually not compatible, and the preferred solution to this problem in biosimulation is to use QM/MM methods in which only a small part of the system is treated quantum mechanically.5-7 However, for the class of metalloproteins for which the ET activity falls within the Marcus regime, it is possible to combine sampling of the protein motions on the nanosecond time scale with calculation of the electronic structure of the entire system. The idea is to sample from an extended classical molecular dynamics (MD) simulation and compute the electronic structure of the full protein in solution at the Density Functional Theory (DFT) level for a statistically significant number of configurations. This allows us to access the vertical ionization energy, which, within the Marcus approach to ET, will permit us to directly compute ab initio redox properties.4,8 Chemical bonding to the metal center, charge transfer to nearby ligands, and long-range polarization are treated in a homogeneous way, free of artificial boundary effects (which is critical particularly in the case of redox properties). We focus on half reactions, namely the oxidation/reduction of a single protein. Our approach is therefore restricted to computation of the free energies relevant for longrange ET, i.e., the redox potential and reorganization energy, and will not give any information on the electronic coupling. This paper reports on an application of the hybrid force field MD/DFT scheme to Rubredoxin (Rd),1 the simplest of the ironsulfur proteins containing only one FeS4 prosthetic group. Rd

10.1021/jp067387y CCC: $37.00 © 2007 American Chemical Society Published on Web 03/28/2007

3970 J. Phys. Chem. B, Vol. 111, No. 15, 2007 has been one of the most popular model system to investigate ET in metalloproteins. Rd is also the ideal prototype protein for a detailed computational analysis. Indeed, it is a small (≈6 kdalton) protein found in bacteria and archea, and its size permits a full ab initio description of the protein embedded in explicit solvent with periodic boundary conditions. The 3D structure has been widely explored with X-ray, neutron, and NMR techniques to very high resolution.1 Crystal structures for both the oxidized and reduced form obtained under similar conditions are available at 1.5 Å resolution.9,10 The fold belongs to the R + β class, with 2 R helices and 2-3 β strands. The active site is characterized by an iron ion that is coordinated by the sulfurs of four conserved deprotonated cysteine residues in an almost regular tetrahedral configuration. These residues reside on two loops, which are the most conserved regions of the protein. The formal redox couple involves Fe(III) and Fe(II) and the protein primary structure modulation of the redox potential has been thoroughly studied experimentally in several different mutants11 and in metal-substituted proteins.12 Several calculations are also available aiming to explain the molecular origin of the differences in the redox potential among different natural rubredoxin variants or mutations. Some important features have been identified, such as the effect of residue mutations11,13-15 in the first and second coordination shell, the role of the second shell backbone, and its influence on the redox potential.16-18 For what concerns the reorganization free energies only some cluster calculations are available which confirm a rigidity and small reorganization energy of the active site.19,20 However, computation of the reorganization energy of the whole aqueous protein with a realistic model has not been attempted yet (to our knowledge) while the available experimental estimates are indirect. There are several natural varieties of Rd. In this study we will compare the redox properties of two of them, namely the mesophilic Clostridium pasteurianum Rubredoxin (CpRd) and the hyperthermophilic Pyrococcus furiosus Rubredoxin (PfRd). The redox potential of these two forms of Rd differ by 60 mV, with PfRd being marginally the more oxidative.21 Anticipating our results, we can correctly predict the redox free energy differences for the single protein half reactions. We also computed the reorganization free energy for the mesophilic and the thermophilic Rd. A difference of about 130 meV is found between the two species. We show that the relative redox potential between CpRd and PfRd can be understood in terms of short-range contributions and can be reproduced with appropriate cluster models for the protein surrounding the metal site. For a proper account of the reorganization free energy, however, solvent contributions are crucial and an accurate description of the dielectric response of the full model system, protein plus solvent, is required. This study has also a second more technical objective. The first principles calculations on the solvated rubredoxin protein have been performed with the density functional package CP2K/ Quickstep.22 This study is therefore also meant to demonstrate that state-of-the-art first-principles tools and methodologies are ready for production calculations (yielding a number that can be verified experimentally) on protein systems in a realistic (solvated) state. II. Methods Redox properties are calculated according to a MD method based on Marcus theory as originally proposed by Warshel.8 As documented in a series of previous publications, we have developed an implementation of this scheme for use in DFT

Sulpizi et al. based ab initio molecular dynamics simulation.23-26 Here we just recall the equations that will be employed in calculations. Starting from atomistic simulations of the system in the oxidized and in the reduced state the oxidation free energy (∆A) and the reorganization free energy (λ) for the half reaction can be computed from ensemble averages of the vertical ionization energy (∆E), according to the following equations:

1 ∆A ) (〈∆E〉red + 〈∆E〉ox) 2

(1)

1 λ ) (〈∆E〉red - 〈∆E〉ox) 2

(2)

Subscripted angular brackets denote averages over equilibrium trajectories of the system in reduced (red) and oxidized (ox) states. In the hybrid scheme used here ∆E is still computed by using DFT but the atomic configurations are extracted from classical simulations. The aim of the method is to combine the long time scale accessibility of the classical model with the quantum-mechanical methods to calculate ionization energies. For this to work we must assume that the MD and electronic structure model are sufficiently close, so that the discrepancy in total energy is fairly uniform over time.4,27 For Rd, with a, from the point of view of inorganic chemistry, normal unstrained tetrahedral iron coordination, these conditions are relatively easy to fulfill. Reliability of this combined approach already has been demonstrated for small organic molecules in different solvents.26 Initial structures were chosen from available high-resolution crystal structures. For both CpRd and PfRd the experimental oxidized and reduced structures are known. This saves us having to go through the relaxation process associated with the oxidation, which could require simulation over long times. Simulation now only has to provide equilibrium runs. In particular for CpRd the 1IRO9 and 1FMH10 pdb database entries were chosen respectively for the oxidized and the reduced form, while for PfRd we used the 1CAA28 and 1CAD28 structure, respectively, for the oxidized and the reduced forms. All the structure were solvated and placed in an orthorhombic box of edges 31.136, 28.095, and 30.502 Å. In the case of 1IRO and 1FMH 678 waters were added, while in the case of 1CAA and 1CAD 676 waters were added (see Figure 1 for a graphical representation of a configuration of the full CpRd system). Despite the different sequence for CpRd and PfRd the total protein charge is the same, -9, so 9 Na+ counterions were included to ensure charge neutrality in the oxidized form. As a consequence the reduced systems carried a total charge of -1. The total number of atoms for the 1IRO and 1FMH systems is 2825 atoms, while the 1CAA and 1CAD systems contained 2826 atoms. The same box size is chosen for both CpRd and PfRd, in order to cancel finite size effects in the relative free energy of the two Rd variants (see ref 23 for a detailed discussion). Indeed, the absolute value of the redox free energy of a single protein has little physical meaning in our approach. The reason is that ionization energies are relative to a reference potential that is determined by the technique used to compute the electrostatic interaction and by the size of the simulation cell.23,24,26 Since for both Cp and Pf Rd the total charge is the same for a given oxidation state and the employed simulation cell is identical, we expect relative values (differences) to be independent of the calculation technique and to be the useful for a comparison with experiments. The classical molecular dynamics simulations were carried out with the Amber8 program package,29 using the SPC/E model for water30,31 and standard RESP parametrization32,33 for the

Calculation of Redox Properties

J. Phys. Chem. B, Vol. 111, No. 15, 2007 3971

Figure 1. Representative MD configuration of CpRd generated from a 1IRO crystal structure as employed in the DFT calculations. The periodically repeated simulation cell, with edges 31.136, 28.095, and 30.502 Å, contains the protein, 678 water molecules, and 9 Na+ counterions. The orange isosurface represents the spin density for the oxidized state (charge 0, spin 5/2) at 0.005 au.

metal center charges. Note that the model for the iron-sulfur complex is different for the reduced and oxidized state with parameters optimized to reproduce the experimental geometry of each oxidation state. The runs were initiated by solvating a protein in the experimental configuration. Counterions were introduced according to the electrostatic energy minima as in the Xleap module from Amber8.29 These structures were equilibrated by using constant pressure-constant temperature (NPT) MD methods for at least 0.5 ns followed by constant volume-constant temperature (NVT) production runs for 5 ns. Selected configurations were extracted every 100 ps and were used to compute the electronic structure properties on the entire solvated protein system. For each snapshot two different DFT calculations were carried out respectively in the oxidized (total charge 0, multiplicity 634) and the reduced forms (total charge -1, multiplicity 534). DFT calculations were performed by using the Becke exchange35 and the Lee-Yang-Parr36 correlation functionals37 with the freely available DFT package CP2K/quickstep,22 which is based on a hybrid Gaussian and plane wave method.38 We employed an efficient wave function optimization scheme that uses an orbital transformation39 technique and analytic Goedecker-Teter-Hutter40,41 pseudopotentials. A Gaussian basis set of QZV2P quality was used for Fe, TZV2P for the four cysteinates bound to the metal, and finally DZVP for all the other atoms. A density cutoff of 280 Ry was used. This setting represents a minimum and reasonable compromise between accuracy and feasibility given the available supercomputing resources. All calculations were carried out on the IDRIS sp4 in France, as part of a DEISA project (http://www.deisa.org/). Cluster calculations were used to isolate the first shell contribution to the redox potential and reorganization free energy. The geometries of these clusters were extracted from protein configurations sampled from the classical trajectory. The smallest cluster (cluster I) consisted of the iron and the four thiolates (from C6, C9, C39, and C42, numbering according to 1IRO pdb entry). For this system the oxidized form has charge

Figure 2. Schematic representation of the three clusters used in the characterization of short-range effects. Cluster I, contains the Fe center, and the four thiolates Cys6, Cys9, Cys39, and Cys42. Cluster II contains the same four thiolates and a portion of the backbone so as to include the NH groups from residues 8, 9, 41, and 42. Cluster III contains in addition to the atoms of cluster II the water molecules whose oxygen is closer then 6 Å from the Fe center. Hydrogen bonds which play a role in the redox properties are indicated by yellow dashed lines.

-1 and multiplicity 6, while the reduced form has charge -2 and multiplicity 5. We also studied an enlarged cluster containing backbone NH groups from nearby residues which form an H-bond to a thiolate sulfur. In particular the NH backbone groups from residues C8, C9, V41, and C42 (sequence and numbering according to 1IRO pdb entry) were included (cluster II). Finally in the largest cluster we considered in this calculation, the cluster II atoms were supplemented by all the water molecules with their oxygen atoms within 6 Å from the iron center (cluster III). See Figure 2 for details. Geometries for all the three clusters were cut out from the sampled configurations along the MD simulation used for the calculations on the full protein. Fragments with severed bonds were terminated with hydrogen atoms placed at a standard bond distance to saturate the valence. For the cluster calculations a QZV2P basis set was used for iron and TZV2P for all the remaining atoms. III. Results and Discussion Redox Properties of the Solvated Protein. The redox properties for CpRd and PfRd, calculated according to eqs 1

3972 J. Phys. Chem. B, Vol. 111, No. 15, 2007

Sulpizi et al.

TABLE 1: Redox Properties of CpRd and PfRd As Calculated for Different Sections of the System and Different Methodsa Cp Pf

cluster I

cluster II

cluste III

full QM

∆A λ ∆A λ

-2.30 0.10 -2.28 0.08

-0.86 0.18 -0.80 0.12

-0.48 0.52 -0.40 0.31

-0.36 0.72 -0.32 0.59

classical ff 6.93 1.27 7.40 0.97

∆∆A ∆λ

-0.02 0.02

-0.06 0.06

-0.08 0.21

-0.04 0.13

-0.47 0.30

a The free energy of oxidation ∆A and reorganization free energy λ (in units of eV) were computed according to eqs 1 and 2. ∆∆A is the difference between ∆A values for CpRd and PfRd, while ∆λ is the corresponding difference between λ values. Clusters are defined in the Methods section (see also Figure 2). The column labeled “classical ff” gives the classical force field values.

Figure 3. Convergence of the oxidation free energies and reorganization free energies as a function of time (ps). Oxidation free energies are indicated by solid lines, reorganization free energies by dashed dashed lines. The curves corresponding to CpRd are in blue, the curves for PfRd are in red. Only relative oxidaton free energies have experimental meaning (for the comparison to experiment recall that reduction potentials and free energies of oxidation have the same sign).

and 2, are summarized in Table 1. The free energy difference between CpRd and PfRd is 0.04 eV, with the higher oxidation free energy for the thermophilic species. This is in good agreement with the experimental measure of 0.06 eV. Usually simulation methods for the calculation of redox properties have to face a host of difficulties connected with the accuracy of the model Hamiltonian, correct sampling of the protein motions, and the effects of periodic boundary conditions and limited system size. The results for our, in this respect, somewhat special system show that the use of DFT, even with a limited number of configurations, which however span a time of 5 ns, can give quantitative results. The accumulated averages of the oxidation free energy and of the reorganization free energies are shown in Figure 3. From this graph it can be seen that our convergence with simulation time is just about sufficient to make the small free energy difference of 40 meV statistically significant. This result also gives us some confidence in the calculation of the reorganization free energy, where no direct experimental measure is available and only indirect estimates have been suggested. The calculated reorganization free energies are 0.72 and 0.59 eV for CpRd and PfRd, respectively. For an impression of the statistical accuracy we refer again to the accumulated averages plotted in Figure 3. Our predictions are in agreement with estimations from the literature, which suggest a value of λ in between 0.4 and 0.7 eV for Rd based on biological ET processes.42,43 Indeed Solomon and co-workers, for example, use a value of 0.7 eV for their estimates of transfer rates19 in self-exchange processes.

Figure 4. Reorganization free energies (plotted as accumulative averages) as calculated from the classical MD trajectories. The solid lines correspond to values for CpRd, while the dashed lines to values for PfRd. Blue curves give the total reorganization free energy, while the contributions from the protein and the electrolitic solution are given in red and green, respectively. The major difference between CpRd and PfRd is in the contribution of the electrolitic solution.

A major advantage of a DFT scheme is that it accounts for electronic relaxation effects in response to oxidation/reduction. Such effects include ligand-metal charge transfer and the adjustment of hydrogen bond strength of coordinated protein residues and solvent. Within a DFT description of the vertical ionization energy for the entire solvated system, as applied here, also the instantaneous equilibration of electronic polarization to the new solute charge distribution is included. This can lead to a significantly lower estimate of the reorganization energy compared to a classical model with fixed charges. We can already see this from the continuum expression of the reorganization energy derived by Marcus, which we give here in the form appropriate for a half reaction

λ)

(

(∆q)2 1 1 2R ∞ 0

)

(3)

R is the radius of the cavity in the dielectric continuum (effective Born radius of the solute) and ∆q the change in charge (∆q ) 1 for our system). ∞ is the dielectric constant at optical (high) frequency and 0 the static dielectric constant for slow (zero frequency) response. The implication of eq 3 for solvents with a high 0 such as water is that λ is approximately inversely proportional to the optical dielectric constant. Standard nonpolarizable classical force fields model electrostatic interactions by fixed point charges which amounts to setting ∞ ) 1. This explains why the reorganization energy when evaluated with use of classical force fields comes out too large. This effect is clearly demonstrated by the estimate of λ for Rd labeled classical in Table 1. These values were obtained by substituting in eq 2 the vertical ionization energy as computed from the force field used to drive the dynamics. The 1.18 and 0.97 eV obtained for CpRd and PfRd, respectively, are about 0.4 eV too high compared to the DFT results. The discrepancy of about a factor of 1.5 is consistent with ∞ ≈ 1.5, which is a typical value for the optical dielectric constant (the squared refractive index) of aqueous systems. The classical analysis, however, can be a useful instrument to decompose the reorganization free energy and it can help in isolating the single contributions from different portions of the system as shown by a number of studies of electron transfer in cytochromes.44-49 The result of such an analysis for rubredoxin is given in Figure 4 showing the various contributions to the (classical) reorganization energy displayed as accumulated averages over time. Water and protein contributions constitute

Calculation of Redox Properties

J. Phys. Chem. B, Vol. 111, No. 15, 2007 3973 TABLE 2: Structure of the Water Solvent Close to the Metal Centera

CpRd reduced oxidized red - ox PfRd reduced oxidized red - ox

Figure 5. Iron-water oxygen radial distribution functions g(r). Distributions relative to the reduced (blue) and oxidized (red) forms of both CpRd (solid line) and PfRd (dashed line) are reported. In the inset the comparison between the g(r) relative to the reduced (blue) and oxidized (red) forms of CpRd calculated for the small production simulation cell (solid line) and a larger cell containing 12 000 water molecules (dashed line) is reported.

75% and 25% of the total reorganization energy for CpRd, and 66% and 34% for PfRd.50 Most of the reorganization energy is therefore due to the solvent. To interpret this result we recall that reorganization energy in the Marcus picture is determined by the inertial (orientational) polarization of the “solvent’’ only, i.e., protein + water (this is what is approximated by eq 3). Because of the constraints on the orientational freedom of polar groups in proteins, the solvent relaxation can be expected to dominate the reorganization. This is why metalloproteins screen the metal center from the water, reducing the reorganization and therefore the activation energy for electron transfer. The relatively high solvent component of ≈70% we find for the λ of Rd must therefore be related to the fact that the metal center is not completely buried within the protein, as happens for other cofactors, like porphyrins, but is quite exposed to the solvent. For cytochrome c, for example, previous calculations have shown a similar contribution to the reorganization energy from the protein and from the solvent.44,48,49 In particular in the reduced form (for both CpRd and PfRd) water can approach the metal center as close as 4.5 Å, as shown by the radial distribution functions depicted in Figure 5. The importance of this observation for Rd has been pointed out already in the discussion of previous classical MD simulations51 and confirmed by high-resolution crystal structure data.10,52 The protein, evidently, provides only partial screening from the dielectric response of the water. The reason is most likely optimization of the balance between reorganization (activation) energy and electronic coupling. Partial exposure of the metal center, while increasing the activation energy, also shortens the ET pathways to the surface. This strengthens the electronic coupling as demonstrated by the calculation in ref 19, leading to the fast rates of electron self-exchange observed experimentally. What is the origin of the 130 meV difference between the values of λ for the mesophilic and thermophilic species? According to Figure 4 the answer to this question must be sought in the solvent contribution. We have therefore subjected the solvent around the iron center to a more detailed structural analysis using the proven statistical tools of classical MD. Figure 5 shows the (water) oxygen-iron radial distribution functions (g(r)) for the two forms of Rd differentiated between oxidation states. Both reduced species show a pronounced short distance peak marking a well-defined first coordination shell extending to the minimum in the rdf at 6.0 Å. This structure has almost vanished in the oxidized state indicating a substantial loss of

N (R ) 6 Å)

p (R ) [0 Å, 6 Å])

N (R ) 8.2 Å)

p (R ) [6.0 Å, 8.2 Å])

4.1 (4.0) 1.6 (1.8) 2.5 (2.2) 3.8 (3.7) 1.9 (1.6) 1.9 (2.1)

0.70 (0.73) 0.32 (0.36) 0.38 (0.37) 0.64 (0.65) 0.39 (0.35) 0.25 (0.30)

24.1 (23.0) 18.5 (18.8) 5.6 (4.2) 23.8 (22.9) 19.5 (19.8) 4.3 (3.1)

0.26 (0.28) 0.15 (0.16) 0.11 (0.12) 0.25 (0.23) 0.10 (0.13) 0.15 (0.10)

a Given are the coordination number (N) and orientational order Parameter (p) for two values of cutoff radius R. p is defined as the product between the water dipole (symmetry axis) and the water charge center-iron center displacement vector averaged over all the waters within a distance R from the metal atom. The two cutoff radii of 6 and 8.2 Å correspond to the first and second shell, respectively (see Figure 5). The N and p values are calculated averaging over 5000 snapshots sampled every ps from the 5 ns trajectory. In parentheses are the values obtained when averaging is restricted to the subset of 50 snapshots used for the computation of vertical ionization energies (Table 1).

hydration upon oxidation. Indeed, determining coordination numbers by spherical integration of the rdf up to 6 Å we find that two water molecules have been expelled from the first coordination shell (see Table 2), which is responsible for a significant fraction of the reorganization energy as will be quantified later on by the cluster analysis. What matters for a comparison between CpRd and PfRd, however, are the subtle differences in the rdf values. On closer inspection of Figure 5 we note that the breakdown of the first solvation shell of oxidized PfRd is somewhat less drastic compared to that of CpRd, namely a decrease of coordination number of ∆N ) 1.9 for PfRd versus ∆N ) 2.5 for CpRd (see Table 2). This difference is consistent with the lower reorganization energy for PfRd. The effect is, however, not very large to the point that it disappears in the noise when the active site solvation numbers are calculated over the limited number of snapshots used in the calculation of λ (number in parentheses in Table 2). A similar trend is observed for the orientations of water in the vicinity of the metal site. The orientational effect is, however, much enhanced compared to the response of coordination number. Indeed the orientation of water molecules can be strongly influenced by the proximity to a protein. For example, it has been shown by experimental studies of chymotrypsin that associated water exhibits an anomalously strong dielectric reorganization, significantly exceeding that of bulk water.52 To characterize the degree of water ordering around the metal center we have introduced the following order parameter. The orientational order of water molecules within a given spherical shell was expressed as the average scalar product p(r) between the water dipole and the vector defined by the iron position and the center of charge of the water molecule. A positive value of +1 corresponds to a water with a dipole moment completely aligned toward the iron center, as for example in the configuration of a solvated anion. Plots of this order parameter are shown in Figure 6. Average values for molecules within the first and second active site solvation shells are again reported in Table 2. As can be seen from Figure 6, water molecules as far as the second solvation shell exhibit an orientational ordering characteristic for solvation of an anion. The maxima in the p(r) also clearly parallel the peaks in the g(r) profiles. The ordering is stronger for the reduced state consistent with the bigger negative charge of the complex. Comparing the two Rd forms we see that the orientational order in CpRd is generally larger than that in PfRd. Most important, as can be deduced from Table

3974 J. Phys. Chem. B, Vol. 111, No. 15, 2007

Figure 6. Water orientation order parameter p(r) (see text) as a function of the oxygen-iron center distance. In the main panel the calculation relative to the reduced (blue) and oxidized (red) forms of both CpRd (solid line) and PfRd (dashed line) is reported. In the inset the comparison between p(r) relative to the reduced (blue) and oxidized (red) forms of CpRd calculated for the production simulation cell (solid line) and a larger cell containing 12 000 water molecules (dashed line) is reported.

2, the water in the first solvation shell of CpRd undergoes a larger orientational change than that in PfRd. This difference remains appreciable also when the calculation is performed considering only the snapshots employed in the QM calculation supporting the hypothesis of a correlation with the larger reorganization free energy of CpRd. The discussion above focused on the interaction of the solvent with the negatively charged metal redox center. However, also the electric charge distribution on the protein surface can strongly affect the water polarization. The metalloproteins in the Rd family can have a variety of charged groups on the protein surface which may compete with the charged metal center in orienting water layers around the protein. In particular in PfRd there are two additional charged residues, (K6 and G47) on the surface located at about 9 and 11 Å from iron respectively replacing the polar T7 and Q48 of CpRd. These residues, characteristic for PfRd, add to the interaction of surface charges with the water around the protein interfering with the ordering by the metal center. The net effect could be a lowering of the dielectric constant for solvent around iron in PfRd, explaining the smaller reorganization free energy compared to CpRd. This effect is, however, more difficult to quantify and must remain, therefore, somewhat speculative. We conclude this section with a technical comment. Errors related to the limited time scale and system size of MD simulations are always a concern. On the basis of the accumulative averages of Figures 3 and 4, the duration of the simulation is probably adequate to ensure convergence. The calculations were performed, however, in a relatively small simulation cell. The number of water molecules accounts just for the first and second solvation shell of the protein. Further shells are shared with the periodic images. It is therefore important to verify to what extend the solvation structure of the protein and in particular of the redox site is influenced by the boundary conditions. For this purpose we carried out additional simulations of an enlarged CpPf model solvated by more than 12 000 solvent molecules. The results for the iron-water radial distribution functions and water orientational order are shown in the insets of Figures 5 and 6. For the rdf there is little or no difference between the smaller and larger systems. The orientational order shows some size effects, which, however, do not alter the picture of the reorientational reorganization sketched above.

Sulpizi et al. Cluster Analysis of Short- and Long-Range Redox Free Energies. For further analysis of the various contributions to the redox free energies we have performed a series of calculations on clusters of increasing size lifted out of the snapshots used for the DFT calculations on the whole solvated protein. The idea is to understand and quantify the role of first and second shell ligands and of approaching water molecules which can form hydrogen bonds to the sulfur atoms. The Fe(SCH3)41-/2- complex, which can be stabilized in free form in certain solvents, has often been used as a minimal cluster model for the redox center in Rd. This is our cluster I of Figure 2. The electronic properties of the iron-sulfur core have been probed both experimentally (electron detachment spectroscopy, photoelectron spectroscopy19,54-57) and theoretically.19,20,54,55,57-60 This already large body of published information has been recently extended by the embedded cluster studies of Hillier and co-workers.61 The oxidation free energies ∆A we compute for cluster configurations sampled from the CpRd and PfRd trajectories are very similar (see Table 1). Differences in the iron first shell of both proteins appear to be minimal as is also suggested by comparing Fe-S distance measurements in crystal structures.9,10,28 Our values for ∆A are negative and in very good agreement with adiabatic detachment energies (ADE) as calculated in previous work.57-61 The reorganization energies involved must be small as is confirmed by the data in Table 1. The apparent rigidity of the tetrahedral iron sulfur complex is also consistent with the small Fe-S bond change upon reduction. Indeed from the crystal structural data only a 0.1 Å shortening of the Fe-S bonds is measured upon reduction. Our results are in agreement with previous DFT calculations on the same models.19,20 In the cluster model one step up in size (cluster II) we have added backbone portions which have NH groups with direct hydrogen bonds to the sulfur atoms (see Figure 2). These hydrogen bonds are supposed to play a role in the redox potential modulation. The evidence for this is mainly derived from 15N NMR spectroscopy17 measurements which show that there is a correlation between HN-Sγ bond lengths and the variation of the reduction potential over a series of 10 CpRd mutants. Also the HN-Sγ bond length was observed to be shorter in the variants with higher reduction potential. Cluster II has been designed to capture this effect. There are in total six HN-Sγ bonds to the ligand cysteines that could play a role. The two more buried cysteines, namely C6 and C39 for CpRd, have two hydrogen bonds each, coming from NH respectively from V8 and C9 and from V41 and C42. These are strong hydrogen bonds. The surface exposed cysteines C9 and C42 also form hydrogen bonds, namely with the NH of Y11 and V44 NH, respectively, but these bonds are comparatively weaker. Only the four strong bonds with the buried cysteines C6 and C39 have been included in cluster II. The ∆A and λ listed in Table 1 have again been computed for cluster configurations extracted from the same snapshots used for the DFT calculations on the whole solvated system. The calculations therefore indirectly include the effect of the protein environment in constraining relative distances and geometries. This is how we can understand that the HN-Sγ distances involving C6 and C39 are shorter in the reduced form (see Table 3 for the average values over the MD simulation). Similarly the hydrogen bonds involving C6 are stronger for the PfRd than for CpRd. From this point of view PfRd behaves like the CpRd mutants with reduction potential higher than in the wild type. Comparing to the results for the isolated Fe(SCH3)41-/2-, we see that the protein backbone has a stabilizing effect on the reduced species leading to less

Calculation of Redox Properties

J. Phys. Chem. B, Vol. 111, No. 15, 2007 3975

TABLE 3: HN-Sγ Averages Distances (Å) and Standard Deviations (Å) in the Oxidized and Reduced Form of CpRd and PfRda CpRd

PfRd

a

C39(Sγ)-C42(HN) C39(Sγ)-L41(HN) C6(Sγ)-C9(HN) C6(Sγ)-V8(HN) C9(Sγ)-Y11(HN) C42(Sγ)-V44(HN) C39(Sγ)-C41(HN) C39(Sγ)-I40(HN) C5(Sγ)-C8(HN) C5(Sγ)-I7(HN) C8(Sγ)-Y10(HN) C41(Sγ)-A43(HN)

oxizidized

reduced

∆

2.56 ( 0.24 2.69 ( 0.23 3.00 ( 0.46 3.02 ( 0.31 3.89 ( 0.20 3.17 ( 0.20 2.68 ( 0.24 2.77 ( 0.21 2.62 ( 0.26 2.82 ( 0.25 2.69 ( 0.19 2.77 ( 0.23

2.38 ( 0.22 2.50 ( 0.17 2.27 ( 0.16 2.44 ( 0.17 3.98 ( 0.24 3.25 ( 0.24 2.37 ( 0.22 2.48 ( 0.15 2.33 ( 0.20 2.49 ( 0.17 2.83 ( 0.33 2.90 ( 0.31

0.18 0.19 0.73 0.58 -0.09 -0.08 0.31 0.29 0.29 0.33 -0.14 -0.13

∆ is the difference between the oxidized and the reduced values.

negative ∆A. From the data in Table 1 it is clear that the modulation of the redox potential is mostly determined by shortrange contributions, coming from the protein backbone. Indeed the cluster II model exactlty reproduces the measured shift in the redox potential between CpRd and PfRd (∆A ) 0.06 eV). On the other hand, from our analysis it also emerges that the protein backbone contribution to the reorganization free energy λ is limited to only 0.08 eV for CpRd with a smaller value for the PfRd. For a correct description of λ, the solvent contribution has to be added in to the calculations. Analysis of the radial distribution functions (Figure 5) revealed a pronounced response of the first hydration shell of the metal center. It was found that as much as two water molecules are added in the reduced state. The implications for the reorganization λ have been investigated by extending the cluster calculations by including the strongly interacting first shell water molecules. Thus, cluster III contains all the waters within 6 Å radius from the iron center (see also Table 2). As expected, we find that the added waters are responsible for a sizable increase in reorganization free energies, when compared to the “dry’’ cluster I and cluster II results. However, the values of λ that we obtain are still about 0.2-0.3 eV too low with respect to the full quantum calculation (Table 1). Note, however, the effect on the difference ∆λ in reorganization of CpRd relative to PfRD, which is in the opposite direction. ∆λ is maximal for the cluster III and the long-range interactions evidently decrease ∆λ. Again we see the importance of longrange contributions requiring us to go beyond the simple cluster description to get quantitative results. IV. Conclusions We have employed a combination of DFT and classsical molecular dynamics to provide quantititative estimates of key redox properties for two variants of rubredoxin. We were able to reproduce the small difference in the reduction potential between CpRd and PfRd giving a validation of the comptational methodology. The hydrogen bond interactions of the ligand cysteines with the NH group of residues V8, C9, V41, and C42 were found to be the biggest factor in the redox potential modulation, with a stronger network of hydrogen bonds leading to more positive reduction potentials. We computed the reorganization free energy for two rubredoxin variants, CpRd and PfRd, and were able to rationalize the origin of the different contributions. It emerges that the protein controls the reorganization free energy in a balance between shielding of the high dielectric solvent medium and preserving some water accesibility to the metal center for a fast electron transfer to the surface.

Discussion of the redox properties of clusters of increasing size showed that a full DFT calculation on the whole solvated system is necessary to correctly include all the polarization effects which are associated with the solvated protein reorganization. Our analysis of the reorganization free energy of the CpRd and PfRd suggests that the sizable difference (130 meV) between these two variants of rubredoxin can be attributed to the different dielectric response of the solvent around the active site. These results underline the importance of a molecular description of the solvent and of a correct inclusion of polarization effects. Our work confirms that in the case of rubredoxin the redox potential modulation is determined by the short-range contributions, while the reorganization free energy is dominated by longrange effects. We also hope to have shown that with modern computational methods, calculations on a full small size protein are within reach and can offer a powerful predictive instrument to quantify properties that are not easily measured by experiment. Acknowledgment. We thank Professor U. Rothlisberger for useful discussions. M.S. is grateful to the EPSRC for financial support. We also thank the DEISA Consortium (co-funded by the EU, FP6 project 508830), for support within the DEISA Extreme Computing Initiative (www.deisa.org), which made this study possible. References and Notes (1) Meyer J.; Moulis J. M. Rubredoxin. In Handbook of Metalloproteins; Messerschmidt, A,. Huber, R., Wieghardt, K., Poulos, T., Eds.; Wiley: New York,2001; pp 505-517 and references therein. (2) Gray, H. B.; Winkler, J. R. Annu. ReV. Biochem. 1996, 65, 537561. Gray, H. B.; Winkler, J. R. Q. ReV. Biophys. 2003, 36, 341-372. Gray, H. B.; Winkler, J. R. Proc. Natl. Acad. Sci. U.S.A. 2005, 102 (10), 3534-3539. (3) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265322. (4) Olsson, M. H. M.; Hong, G.; Warshel, A. J. Am. Chem. Soc. 2003, 125, 5025-5039. (5) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227-249. (6) Colombo, M. C.; Guidoni, L.; Laio, A.; Magistrato, A.; Maurer, P.; Piana, S.; Ro¨hrig, U.; Spiegel, K.; Sulpizi, M.; VandeVondele, J.; Zumstein, M.; Ro¨thlisberger, U. Chimia 2002, 56, 13-19. (7) Gao, J.; Truhlar, D. G. Annu. ReV. Phys. Chem. 2002, 53, 467505. (8) Warshel, A. J. Phys. Chem. 1982, 86, 2218-2224. (9) Dauter, Z.; Wilson, K. S.; Sieker, L. C.; Moulis, J. M.; Meyer, J. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8836-8840. (10) Min, T.; Ergenekan, C. E.; Eidsness, M. K.; Ichiye, T.; Kang, C. Protein Sci. 2001, 10, 613-621. (11) Xiao, Z.; Lavery, M. J.; Ayhan, M.; Scrofani, S. D. B.; Wilce, M. C. J.; Guss, J. M.; Tregloan, P. A.; George, G. N.; Wedd, A. G. J. Am. Chem. Soc. 1998, 120, 4135-4150. (12) Maher, M.; Cross, M.; Wilce, M. C. J.; Guss, J. M.; Wedd, A. G. Acta Crystallogr. Sect. D 2004, 60, 298-303. (13) Meyer, J.; Gaillard, J.; Lutz, M. Biochem. Biophys. Res. Commun. 1995, 212, 827-833. (14) Kummerle, R.; Zhuang-Jackson, H.; Gaillard, J.; Moulis, J. M. Biochemistry 1997, 36, 15983-16991. (15) Park, I. Y.; Eidsness, M. K.; Lin, I. J.; Gebel, E. B.; Youn, B.; Harley, J. L.; Machonkin, T. E.; Frederick, R. O.; Markley, J. L.; Smith, E. T.; Ichiye, T.; Kang, C. Proteins 2004, 57, 618-624. (16) Eidsness, M. K.; Burden, A. E.; Richie, K. A.; Kurtz, D. M.; Scott, R. A.; Smith, E. T.; Ichiye, T.; Beard, B.; Min, T.; Kang, C. Biochemistry 1999, 38, 14803-14809. (17) Lin, I. J.; Gebel, E. B.; Machonkin, T. E.; Westler, W. M.; Markely, J. L. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 14581-14586. (18) Tan, M. L.; Kang, C.; Ichiye, T. Proteins: Struct., Funct., Bioinf. 2006, 62, 708-714. (19) Kennepohl, P.; Solomon, E. I. Inorg. Chem. 2003, 42, 696-708. (20) Sigfridsson, E.; Olsson, M. H. M.; Ryde, U. Inorg. Chem. 2001, 40, 2509-2519. (21) Stephens, P. J.; Jollie, D. R.; Warshel, A. Chem. ReV. 1996, 96, 2491-2513.

3976 J. Phys. Chem. B, Vol. 111, No. 15, 2007 (22) The CP2K developers group, 2005: http://cp2k.berlios.de/. VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Comput. Phys. Commun. 2005, 167, 103-128. (23) Tateyama, Y.; Blumberger, J.; Sprik, M.; Tavernelli, I. J. Chem. Phys. 2005, 122, 234505 (1-17). (24) Blumberger, J.; Sprik, M. Theor. Chem. Acc. 2006, 115, 113126. (25) Blumberger, J.; Tavernelli, I.; Klein, M. L.; Sprik, M. J. Chem. Phys. 2006, 124, 064507 (1-12). (26) VandeVondele, J.; Sulpizi, M.; Sprik, M. Angew. Chem., Int. Ed. 2006, 45, 1936-1938. (27) The bias created by averaging over configurations generated by a Hamiltonion deviating from the true Hamiltonian can be eliminated by including a weight factor consisting of the Boltzmann exponent of the difference in energy of the target Hamiltonian and sampling Hamiltonian.4 However, for our system, these energy differences are dominated by errors in the DFT description of the intermolecular interactions stabilizing the protein structure. The classical force field model, which has been designed for this purpose, is in fact much more reliable. Correcting for the difference between the DFT and classical model would therefore deteriorate rather than improve the accuracy of our results for the redox properties, which is why we have not applied such a procedure. (28) Day, M. W.; Hsu, B. T.; Joshua-Tor, L.; Park, J. B.; Zhou, Z. H.; Adams, M. W.; Rees, D. C. Protein Sci. 1992, 1, 1494-1507. (29) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. S.; Kollman, P. A. AMBER 8, University of California, San Francisco, 2004. (30) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269-6271. (31) Straatsma, T. P.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 5876-5886. (32) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269-10280. (33) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993, 115, 9620-9631. (34) Phillips, W. D.; Poe, M.; Weiher, J. F.; McDonald, C. C.; Lovenberg, W. Nature 1970, 227, 574-577. (35) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100. (36) Lee, T. C.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785789. (37) The high spin multiplicity for the system discussed here might raise questions about the validity of our approximation to the exchange and correlation functional (BLYP). However, in our case both oxidized and reduced states have high multiplicity. Vacuum calculations for similar systems indicate that DFT can give quite good results in contrast to, for example, haemes, where high and low spin states are compared.

Sulpizi et al. (38) Lippert, G.; Hutter, J.; Parrinello, M. Mol. Phys. 1997, 92, 477488. (39) VandeVondele, J.; Hutter, J. J. Chem. Phys. 2003, 118, 43654369. (40) Goedecker, S.; Teter, M.; Hutter, J. Phys. ReV. B 1996, 54, 17031710. (41) Hartwigsen, C.; Goedecker, S.; Hutter, J. Phys. ReV. B 1998, 58, 3641-3662. (42) Blankman, J. I.; Shanazad, N.; Miller, C. J.; Guiles, R. D. Biochemistry 2000, 39, 14806-14812. (43) Fraga, E.; Webb, M. A.; Lopponow, G. R. J. Phys. Chem. 1996, 100, 3278-3287. (44) Simonson, T.; Perahia, D. J. Am. Chem. Soc. 1995, 117, 79878000. (45) Miyashita, O.; Okamura, M. Y.; Onuchic, J. N. J. Phys. Chem. B 2003, 107, 1230-1241. (46) Basu, G.; Kitao, A.; Kuki, A.; Go, N. J. Phys. Chem. B 1998, 102, 2076-2084. (47) Basu, G.; Kitao, A.; Kuki, A.; Go, N. J. Phys. Chem. B 1998, 102, 2085-2094. (48) Miyashita, O.; Go, N. J. Phys. Chem. B 2000, 104, 7516-7521. (49) Muegge, I.; Qi, P. X.; Wand, A. J.; Chu, Z. T.; Warshel, A. J. Phys. Chem. B 1997, 101, 825-836. (50) Counterions contribute very little to the total reorganization energy (10% for CpRd and 6% for PfRd). Over the sampled time scale of 5 ns, none of them was seen to approach the metal center closer than 7 Å, which could perturb the metal coordination properties. (51) Yelle, R. B.; Park, N. S.; Ichiye, T. Proteins: Struct., Funct., Genet. 1995, 22, 154-167. (52) Park, I. Y.; Youn, B.; Harley, J. L.; Eidsness, M. K.; Smith, E.; Ichiye, T.; Kang, C J. Biol. Inorg. Chem. 2004, 9, 423-428. (53) Mertz, E. L.; Krishtalik, L. I. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 2081-2086. (54) Kennepohl, P.; Solomon, E. I. Inorg. Chem. 2003, 42, 679688. (55) Kennepohl, P.; Solomon, E. I. Inorg. Chem. 2003, 42, 689-695. (56) Yang, X.; Wang, X. B.; Fuand, Y. J.; Wang, L. S. J. Phys. Chem. 2003, 107, 1703-1709. (57) Niu, S.; Wang, X. B.; Nichols, J. A.; Wang, L. S.; Ichiye, T. J. Phys. Chem. A 2003, 107, 2898-2907. (58) Bair, R. A.; Goddard, W. A., III J. Am. Chem. Soc 1978, 100, 5669-5676. (59) Mouesca, J. M.; Chen, J. L.; Noodelman, L.; Bashford, D.; Case, D. A. J. Am. Chem. Soc 1994, 116, 11898-11914. (60) Koerner, J. B.; Ichiye, T. J. Phys. Chem. B 1997, 101, 3633-3643. (61) Sundararajan, M.; Hillier, I. H.; Burton, N. A. J. Phys. Chem. 2006, 110, 785-790.