Calculated Reaction Cycle of Cytochrome c Oxidase - The Journal of

Oct 10, 2007 - Department of Theoretical Chemistry, Lund University, Chemical Centre, P. O. Box 124, S-221 00, Lund, Sweden. J. Phys. Chem. B , 2007, ...
0 downloads 0 Views 284KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 12543-12550

12543

Calculated Reaction Cycle of Cytochrome c Oxidase Markus Kaukonen† Department of Theoretical Chemistry, Lund UniVersity, Chemical Centre, P. O. Box 124, S-221 00, Lund, Sweden ReceiVed: January 23, 2007; In Final Form: July 14, 2007

The catalytic cycle of cytochrome c oxidase has been simulated by means of quantum mechanical calculations. The experimental energetics of the catalytic cycle is nearly reproduced. The atomic structures of the intermediates are suggested. In particular, the structures of nonactive “resting” intermediates are proposed.

Introduction Cytochrome c oxidase (CcO) is the fourth enzyme in the aerobic respiratory chain located (in eukaryotic cells) in the inner mitochondrial membrane.1,2 It reduces oxygen to water. The energy gain in this redox reaction, O2(gas) + 4CytC2+(aq) + 4H3O+(aq) f 4CytC3+(aq) + 6H2O(aq), is about 220 kJ/mol (at pH 7). CcO converts this chemical energy to the energy of the electrochemical field across the mitochondrial inner membrane. This energy is later harvested by adenosinetriphosphatase (ATPase) to form adenosine 5′-triphosphate (ATP). The structure of fully oxidized and fully reduced CcO has been known in atomic detail since the mid 1990s.3,4 Recently, some of the water molecules present in the working enzyme have been located by X-ray diffraction techniques.5,6 The structure of the binuclear center at which the water-forming reaction takes place possesses nearly identical X-ray structures independent of the oxidation state and the species. However, the mechanism for the proton-pumping activity of the enzyme remains enigmatic. Protons are invisible with the traditional crystallographic methods making the assignment of the ligands of the important metal atoms difficult. Crystallographic studies are also hindered by the short lifetime (20 µs to 2 ms) and instability of the intermediates in the reaction cycle. These possess specific ligands in the binuclear center (hydroxyl anion, oxygen, water) in combination with particular water (chain) arrangements and are difficult to trap. A further complication arises because this enzyme can enter a resting state when oxidized (and probably also when partially reduced). This/these state(s) does/do not contribute to proton pumping.7 In this article, the states participating in the pumping cycle are labeled as “working” (Oworking, Eworking), while the nonpumping states are labeled “resting” (Oresting, Eresting). In this paper, recent experimental results from the Wiksto¨m group on CcO are used as a starting point. It is assumed that one proton is pumped for each electron (and for each chemistry proton) arriving at the binuclear center.8 The pumping steps are P f F, F f Oworking, Oworking f Eworking, Eworking f R (Figure 1). Thus, two protons are pumped in both reductive and oxidative parts of the cycle. In the past, quantum mechanical (QM) simulations have examined the oxygen binding and the O2 bond-splitting reactions over the reductive part of the cycle and over the whole catalytic cycle yielding locations of water atoms and amino acids in the †

Corresponding author. Tel.: (46)46-2224502. Fax: (46)46-2224543.

vicinity of the binuclear center or in the proton-conducting channels.9-15 Apparently, the current results for the chemistry agree reasonably well with older studies by Moore for the reductive part of the cycle.12 There is also agreement with Siegbahn et al. in that with atom O2A (the oxygen atom not making a hydrogen bond to aspartic acid) at heme a3, R-propionate could serve as a proton-storing site.13 However, the results obtained here are in conflict with the theoretical results by Siegbahn et al. and the ones by Popovic´ and co-workers (His291 is postulated to constitute the proton pump element).13,16,17 Recent theoretical work by Fadda et al. does not support the idea of His291 as a proton pump element.18 Inspired by work of Popovic´ and Stuchebrukhov, the possibility that the QM system would be “lacking” one proton at CuB neighboring His291 was considered. This would mean that His291 would carry formal charge of -1|e| during the cycle.16,17 With the current QM system, this did not bring the even energy distribution to the catalytic cycle. Additionally, the enthalpy of their intermediate steps F f Fp and E f Ep was calculated, but it was found here that deprotonating of His291 was energetically very unfavorable (for notation see Figures 4 and 5 in ref 16). Thus, the current study does not support their model for the proton-pumping mechanism. The discrepancy may arise because of the fully classical treatment of the system in their study with rather arbitrary partial charges on metals and their ligands. Furthermore, here is a comment to another proposed pumping mechanism by Siegbahn and co-workers.13 In their model, no particular “switch” or “gating element” is needed, because depending on the redox (and protonation) state of the binuclear center it is energetically more favorable to add a proton either to R-propionate of heme a3 or to the binuclear center. This study failed to reproduce that result. For example, one tested the feasibility of protonating R-propionate instead of proton arriving to the binuclear center in their proposed O f ET catalytic step (for the notation see Figure 9 in ref 13). Their result is that the step O f ET is exergonic by of the order of 60 kJ/mol. On the contrary, the result of this study for the same step in the bigger QM system is that O f ET is endergonic. One can address this discrepancy to the fact that Siegbahn et al. did not have the hydrogen-bonding partners of R-propionate of heme a3 in their system (Asp364, His368, H2O) and assumed a constant hydrogen bond strength from R-propionate to Asp364. Here, it is found that this assumption may be incorrect, because depending on the protonation and redox state of the binuclear center the

10.1021/jp070578w CCC: $37.00 © 2007 American Chemical Society Published on Web 10/10/2007

12544 J. Phys. Chem. B, Vol. 111, No. 43, 2007

Kaukonen

Figure 1. The schematic presentation of the proposed catalytic cycle of cytochrome c oxidase. Energies are in kJ/mol. The energies are based on QM total energies and experimental free energies of H2O, H3O+, H+, e- and cytochrome c. The electron arriving at the binuclear center is marked as “e-” in the figure. The chemistry and pumped protons are indicated by Hc and Hp, respectively. The calculated unpaired spins are indicated by arrows.

partial charges at R-propionate vary and therefore the hydrogen bond strengths also change substantially. A strong motivation for studying the water-forming chemistry of CcO (i.e., the “power supply” of this enzyme) is that this information is an essential prerequisite for understanding the pumping mechanism in CcO. If the underlying chemistry is incorrect, then the value of the proton-pumping study is questionable. This may be the case for the study of the pumping mechanism by Olsson et al., which is based on the work of Siegbahn.13,19 The remainder of this article is organized as follows. The calculation method and approximations are described in the next section, followed by the results. The results section begins with general considerations on the reaction cycle (Figure 1) followed by studies of its particular intermediates (R, A, P, F, Oworking, Oresting, Eworking, Eresting). The atomic geometries, spin densities, and energy differences are the observables that are compared with the experimental data. Method Quantum Mechanics Framework. In this work, a parallel implementation of TURBOMOLE QM program package with density functional theory and B3LYP exchange correlation functional is used.20-22 Thus, the electrons are described by oneelectron wave functions, and they move in an average field of all the electrons in the system. Exchange and correlation effects are embedded into the partially empirical B3LYP exchange correlation potential. The TURBOMOLE computer code has wellcontrolled convergence properties, and it is parallelized efficiently.23 The resolution of the identity approximation is utilized. In this study it reduces the computational cost approximately by a factor of 10. Split valence plus polarization basis sets are used for for C, N, and O atoms; split valence basis for H atoms and a more complete triple-ζ split valence plus polarization basis are used for for Fe and Cu.24 All electrons are taken into account; no pseudo-potentials are used. Unrestricted spin polarization is allowed, hence, high-spin states or antiferromagnetic spin-

coupling can be accounted for (as will be the case with the A-intermediate). The edge atoms of QM system are the Cβ atoms of the amino acids. These are held fixed at their crystallographic positions during the simulations. The broken C-C bonds of these fixed edge atoms are saturated by hydrogen atoms. The positions of all the other atoms are relaxed by conjugate gradient minimization. This minimization is continued until the maximum force in the system is less than 10-3 au (≈5 kJ/mol × Å) and the change in the total energy less than 10-5 au (≈0.026 kJ/mol). The effect of the surrounding protein is taken into account by approximating it as a dielectric continuum with the relative dielectric constant r ) 4. To this end, the conductor-like screening model is utilized.25 The cavity for the conductor-like surface was defined by using a union of spheres around each atom with radius Rsolv + Ri where Rsolv was chosen as 1.3 Å and Ri for each atom type: RC ) 2.000 Å, RCu ) 2.223 Å, RFe ) 2.223 Å, RH ) 1.300 Å, RN ) 1.830 Å, and RO ) 1.720 Å. The dielectric continuum approach for energy calculations is most valid when the total charge of the QM system in question remains constant. When the total charge between adjacent states changes (as it does in pKa and electron affinity calculations, for instance), one should consider taking into account the polarization effects of the environment in a more detailed fashion.26,27 The spin density estimates are based on electrostatic potential (ESP) fitting of the TURBOMOLE 5.9 package. Ten surfaces of overlapping spheres around atoms are generated using the radius R ) Radd + Ri with following radii Ri for the atom types: H, 2.04; C, 2.55; N, 2.55; O H, 2.38; Cu, 3.40; Fe, 3.40 Å. The Radd of the ten shells have values of -0.302, -0.201, -0.101, 0.000, 0.101, 0.201, 0.302, 1.30, 2.30 and 6.30 Å. On each of the surfaces (at 75 000 points), the electrostatic potential arising from the fitted ESP charges and the one originating from the true QM calculation are required to be as close to each other as possible.

Calculated Reaction Cycle of Cytochrome c Oxidase

J. Phys. Chem. B, Vol. 111, No. 43, 2007 12545

Figure 2. The heme a3 part of the QM system is shown in detail (in the R state). The R-propionate of heme a3 is protonated. The residues hydrogen bonded to the propionates are included in the QM system (W126, R438, H368, and D364). The histidine covalently bonded to iron (H376) is also taken into account. Hydrogen bonds are shown with dashed lines. In the figure, PRA and PRD stand for R- and δ-propionates, respectively.

The Selected QM Fraction. A recent X-ray structure with numerous water molecules of bovine CcO (1V55) is selected in its reduced state as a starting structure for the calculations.6 Amino acid numbering in this article follows the bovine scheme. The QM system used in this study is depicted in Figures 2 and 3. The following components are included in the QM system. The heme a3 is included up to the hydroxyl group at atom C12 on its long farnesyl chain (Figure 2). The proximal histidine 376 ligated to heme a3 iron is included (Figure 2). CuB with its three histidine ligands (240, 290, 291) and the crosslinked tyrosine 244 are included (Figure 3). The covalent His240-Tyr244 bond is assumed (Figure 3). The δ-propionate of heme a3 is unprotonated. The R-propionate is protonated at atom O2A (i.e., the oxygen atom not making a hydrogen bond to Asp364). The arginine 438 and tryptophan 126 groups (making hydrogen bonds to the δ-propionate); the neutral histidine 368 and protonated (neutral) aspartic acid 364 groups (making hydrogen bonds to the R-propionate) are included (for discussion see the next section). The edge atoms of the QM system are the Cβ atoms of the amino acids. These are held fixed at their crystallographic positions during the simulations. Two crystallographic water molecules are taken into account. These are the ones that lay between the propionate groups (and make an additional H-bond to His291) and the H2O that is hydrogen bonded to the farnesyl hydroxyl group (believed to belong to the K-channel pathway).6 In some calculations, additional water molecule(s) are put in the space between Cu, Fe, and Tyr244. In the calculations, the position of His368 between the Mg cation and R-propionate of heme a3 changes substantially, and it sometimes even makes a hydrogen bond to Arg438. The movements of this histidine are likely to be an artifact of the calculation due to the limited size of the QM volume (in particular, the Mg cation and its other ligands are not included in the calculation). Therefore these movements of His368 are

Figure 3. The CuB part of the QM system (in the R state). For schematic orientation, in the figure the direction to which protons are pumped is roughly “upward” (the positive side or the intermembrane space side) and the matrix side (or the negative side) of the mitochondrial membrane is “downward” in the picture. Hydrogen bonds are shown with weak lines.

not commented in the subsequent sections. However, they may reflect the changes in the nature of bonding in the region of R-propionate of heme a3, which in the real system also includes many water molecules that are not accounted for in this study. The Total Charge of the Selected QM Fraction. Adding up formal charges of the selected QM system with deprotonated propionates would yield zero total charge. Blomberg et al. have proposed that the presence of an additional proton in the binuclear center is necessary to obtain a reasonable energy barrier between the states A and Pm.10,11 The results of this study indicate that this assumption is plausible. The experimental midpoint potential for heme a is of the order of 350 mV.28 The calculated midpoint potential (E0) for the Pm f Pr electrontransfer yields +278 mV, which is perhaps slightly less than expected, with the “extra” proton. On the other hand, the calculated E0 for the Pm f Pr is ∼ -1 V without the extra proton, a result that is clearly far off from the experimental value. However, this evidence is weakened by the fact that the model system is in a dielectric continuum, hence, the polarization of protein, lipid bilayer and water is taken into account in a very coarse grained fashion. Additionally, it may be possible that in the Pm state a proton would temporarily move to either of the propionates of heme a3. The source of this proton could be, for instance, the Arg438-Arg439 pair, or the His368 group. In this study, an extra proton is set at the R-propionate of heme a3 with a resulting total charge of +1|e| of the QM system (except that the Pr state has zero total QM charge).

12546 J. Phys. Chem. B, Vol. 111, No. 43, 2007 Finally, it is found with various tested intermediates that without this extra proton the energy drop in the oxidative half of the cycle (R f O) tends to be too large. This observation is in line with the recent QM results by Blomberg.29 Calculation of the Energetics of the Catalytic Cycle. From the QM-calculations (enthalpy part) and experimental results by other groups (entropy contributions and solvation energies), one obtains the Gibbs free-energy change at ∆G0 of a given reaction by subtracting the total formation energy of the reactants from the one of the products. As an example consider the Pm f F step of the reaction cycle of CcO:

Pm + CytC2+ + H3O+ f Pm + CytC2+ + H2O + H+(aq, exp) f Pm + CytC2+ + H2O + H+(vac, calc) f F + CytC3+ + H2O Here Pm and F are the calculated QM total energies of intermediates of the cycle in vacuum. All the cytochrome c data is experimental (in water). Similarly, all the water and hydronium cation data is experimental (in water). The experimental midpoint value for cytochrome c in water yields E0m ) +4.43 V + 0.26 V ) 4.69 V.30,31 Converting this to Gibbs free energy by ∆G0 ) -FE0m, where F is the Faraday 0 3+ ) -452.5 kJ/ constant 96.48 C/mol, gives ∆GcytC 2+fcytC 32 mol. Thus, in the calculations an electron is taken from 452.5 kJ/mol below the vacuum level. Similarly, the experimental Gibbs free energy of hydration for a proton (H+(aq, exp) f H+(vac, exp)) at pH 0 and at the room temperature ∆G0 ) -1105 kJ/mol is used.33 The gasphase free energy of proton -26 kJ/mol is added to get H+(aq, exp) f H+(vac, calc). The thermodynamic equilibrium implies that H3O+(aq) f H2O(aq) + H+(aq) as 0 kJ/mol in 55 M water. Therefore, in the calculations, the proton is taken -1105 - 26 ) -1131 kJ/mol below the calculated energy of a proton in vacuum (0 kJ/mol). It is assumed that the vibrational energy and entropy remain the same for the proton at the initial state (the hydronium cation in water) and in the final state (the hydrogen atom in the protein, for instance, here in the F state). Errors of the order of 10 kJ/mol are expected, because the vibrational spectra of the intermediates are not calculated for the large QM-system. Furthermore, there are uncertainties in estimating the experimental solvation Gibbs free energy of proton as discussed by Liptak and Shields.34 To compare the energetics of a water molecule in bulk water and in protein, the experimental enthalpy of a water molecule in 55 M water H ) -285.8 kJ/mol is used.35,36 Furthermore, it is assumed that the entropy of a H2O molecule is the same in the protein as it is in ice, -TS ≈ -10 kJ/mol.37 The fact that the entropy of a water molecule may be nearly the same as in the protein and water has been suggested by Denisov et al. and by Olano and co-workers.38,39 To estimate the binding energy of oxygen in the R f A step, the experimental entropy contribution of O2 in the gas phase (≈ -60 kJ/mol at the room temperature and 1 bar) is assumed.35,36 This term consists of the kinetic and rotational entropy (the vibrational end electric entropies are virtually zero). The kinetic, rotational, and vibrational energy of Ogas 2 Ekin+rot+vib ) 7/2 kbT is of the order +10 kJ/mol. The experimental change in the Gibbs free energy +16.4 kJ/mol in aq 36 Ogas 2 f O2 is adopted.

Kaukonen A preliminary calculation with a very small QM system (Cu, 3NH3, Fe, 5NH3, (O2)) was carried out to estimate the change in the zero-point energy and vibrational entropy upon oxygenbinding reaction (the R + O2(aq) f A). The result that these factors disfavor the oxygen binding by 50 kJ/mol is in agreement with a theoretical QM/MM study of O2 binding in hemerythrin.40 However, in that QM/MM study these effects were compensated by the protein environment (the van der Waals (vdW) interaction between the bound dioxygen and the protein atoms). It is therefore chosen to neglect these cancelling factors in this study. Results and Discussion General Observations on the Catalytic Cycle of CcO. The calculated energetics and corresponding intermediate atomic configurations of the catalytic cycle of cytochrome c oxidase are summarized in Figure 1. The detailed energetics of the calculated intermediates are presented in Table 1. In the results, the changes in entropies are not calculated, but they are partly implicitly taken into account by adopting the experimental values for the free energies of O2, H2O, H3O+, H+ and cytochrome c (for discussion see the previous section). Next, we consider the catalytic cycle of CcO in more detail. In general, the unpaired spin density is localized at iron, copper, and their O or OH- ligands and in some cases also at the His240-Tyr244 moiety. A set of configurations for each state of the cycle and their energy with respect to the lowest energy one of that intermediate are listed in Table 1. Starting from the reduced state, the changes in free energies in the reactions R f A f Pm steps are close to zero. This result assumes that the increase in the vibrational entropy and zero-point energy is compensated by the van der Waals interaction between the bound dioxygen and the protein atoms (see the subsection entitled Calculation of the Energetics of the Catalytic Cycle). In the calculated F state, it is energetically more favorable to return the proton to the Tyr244 group than to the hydroxide ligand of copper. This result agrees with FTIR spectroscopic observations by the Rich group that suggest reprotonation of the Tyr244 residue occurs during the Pm f F step.41,42 For the oxidized states, the calculated results agree with early calculations by Moore et al.12 The higher energy of Oworking/ CuOH- FeOH- TyrOH compared with Oresting/CuOH- H+ FeOH- TyrO- is a consequence of the electrostatic repulsion between the two OH- groups. In the Eworking state, an electron and a proton arrive at copper and form a water molecule that dissociates, leaving Cu+1 without a ligand and Fe with a hydroxide ligand. However, the most favorable E state (Eresting) in terms of energy is such that the iron has a water ligand and the copper has a hydroxide ligand (Figure 1). It can be reached from the Oresting state by protonation of Tyr244 and a slight proton movement in the binuclear center. It is known experimentally that there are two “proton wires” to supply protons into the binuclear center (so-called D- and K-channels). The D-channel leads to the Glu242 group and the K-channel to the Tyr244 group. Mutagenetic experiments point out that the D-channel is active during the oxidative half of the cycle, while the K-channel opens in the reductive part.1,8,43 Under the assumption that the K-channel is only operating when reprotonating Tyr244, the current QM calculations suggest that the K-channel would be active in the Pm f F and Oresting f Eresting steps. Because this conclusion does not agree with the experimental data, more advanced computational reaction path studies including the water chains will be needed.8,44

Calculated Reaction Cycle of Cytochrome c Oxidase

J. Phys. Chem. B, Vol. 111, No. 43, 2007 12547

TABLE 1: The Calculated Energetics of the Spin States of the Intermediates of the Catalytic Cycle of Cytochrome c Oxidasea species

total spin (S) R (the reduced state) 0 css 0 oss 1 triplet 2 2 2

HisTyrOH/Cu/Fe HisTyrOH/Cu/Fe HisTyrOH/Cu/Fe HisTyrOH0.0/Cu0.01u/Fe3.83u HisTyrOH/Cu-H2O/Fe HisTyrOH/Cu/Fe-H2O HisTyrOH/Cu-O-O-Fe HisTyrOH/Cu0.18d-O0.55d-O0.24d-Fe0.96u HisTyrOH/Cu-OsO-Fe HisTyrOH/Cu-OsO-Fe

relative energy [kJ/mol] 50 49 20 0 >100 13

A (dioxygen bound between Fe and Cu) 0 css 0 oss 1 triplet 2

>80 3.6 14 >100

Pm (OsO bond broken, “peroxy”) 0css 0oss 1 triplet 2

>100 >100 -8.0 28

HisTyrO0/Cu-OH/FedO HisTyrO0/Cu-OH/FedO HisTyrO00.95d/Cu0.64u-OH0.17u/Fe1.29udO0.77u HisTyrO0/Cu-OH/FedO

Pr (O-O bond broken, “reduced peroxy”) HisTyrO-/Cu-OH/FedO 0.5 HisTyrO-0.07u/Cu0.61u-OH0.17u/Fe1.32udO0.78 1.5 HisTyrO-/Cu-OH/FedO 2.5

34 0 >100

F (“ferryl”) 0.5 1.5 2.5 0.5 1.5 2.5

39 0 >50 33 33 >50

HisTyrOH/Cu-OH/Fe-OH HisTyrOH/Cu-OH/Fe-OH HisTyrOH/Cu0.50u-OH0.27u/Fe1.03u-OH0.08u HisTyrOH/Cu-OH/Fe-OH

Oworking (oxidized, working, “fast”) 0 css 0 oss 1 triplet 2

>100 >100 60 >100

HisTyrO-/Cu-H2O/Fe-OH HisTyrO-/Cu-H2O/Fe-OH HisTyrO-/Cu-H2O/Fe-OH HisTyrO-/Cu0.63u-OH0.13u/Fe2.88u-H2O0.08u

Oresting (oxidized, resting, “slow”) 0 css 0 oss 1 triplet 2

>150 24 23 0

HisTyrOH/Cu-OH/FedO HisTyrOH0.07u/Cu0.61u-OH0.17u/Fe1.34udO0.77u HisTyrOH/Cu-OH/FedO HisTyrO-/Cu-H2O/FedO HisTyrO-/Cu-H2O/FedO HisTyrO-/Cu-H2O/FedO

HisTyrOH/Cu0.01d/Fe0.93u-OH0.13u HisTyrOH/Cu/Fe-OH HisTyrOH/Cu/Fe-OH HisTyrOH/Cu-H2O/Fe-OH

Eworking (one electron reduced, working, “fast”) 0.5 1.5 2.5 all

Eresting (one electron reduced, resting, “slow”) 0.5 HisTyrOH/Cu0.60u-OH0.16u/Fe0.09u-H2O0.0 HisTyrOH/Cu-OH/Fe-H2O 1.5 HisTyrOH/Cu-OH/Fe-H2O 2.5

24 >100 >100 unstable 0 19 4.8

a The lowest energy spin state of each intermediate with the same number of protons has the energy 0 kJ/mol, except for the R, A, and Pm states where the energies are with respect to the R, S ) 2 state. In the Oresting and Oworking, the energies are with respect to the lowest energy one Oresting, S ) 2. Similarly, all the energies of the Eresting and Eworking are with respect to the Eresting S ) 0.5 state. Energies are in kJ/mol. The unpaired spin density of each intermediate with the lowest energy is shown as a superscript and the spin “up” and “down” densities projected to atoms or ligands are marked as “u” and “d”, respectively. In the table, “oss” stands for an open shell singlet and “css” for a closed shell singlet.

In the following subsections the calculated geometries and other properties of the intermediates are described in more detail. Reduced State: R. The catalytic cycle starts with the copper and iron of the binuclear site in their formally reduced forms: Cu+1 and Fe+2. The spin state S ) 2 with four unpaired electrons at the iron atom is found to be energetically the most favorable one. This is in agreement with classical experimental Raman and optical spectroscopy results, which all conclude that heme a3 is a high-spin species.45-47 The S ) 1 triplet state is of the order 20 kJ/mol higher in energy. By the energetics it seems possible that there could be a water molecule between the iron and copper atoms. The lowest structure with this extra water

molecule is of the order of 10 kJ/mol higher in energy, but a water chain from the D-channel may bring this energy down to a reasonably low level. In the QM relaxed structure, the ironcopper distance is 5.12 Å in the S ) 2 state. The iron of heme a3 has moved out of the porphyrin plane in toward the proximal histidine 376, similar to the observations by Rovira et al. for small porphyrin systems.48 The distance from the copper atom to each of the nitrogen atoms belonging to the three neighboring histidine groups are 1.97 Å(H291), 1.95 Å(H240), and 2.13 Å(H290). This agrees well with EXAFS data for the fully reduced enzyme by Ralle and co-workers.49 They find Cu-N distances of 1.92, 1.92, and 2.10 Å for a bo3 type oxidase. Thus one can

12548 J. Phys. Chem. B, Vol. 111, No. 43, 2007

Figure 4. The A intermediate is an S ) 0 open shell singlet biradical state. The dioxygen pair is bound to the iron atom. There is also bonding with it and the copper atom. The unpaired spin density is shown as two isosurfaces. One is located at the iron and the opposite one at the copper and at the two bridging oxygen atoms resulting an antiferromagnetic spin-coupling. They form a π-bond system.

assign the longest bond (2.13 Å in the QM calculation) to the one between the N atom of the His290 group and the copper atom. Dioxygen Bound State: A. The dioxygen bound intermediate (A state) forms in about 10 µs after the R state in the working enzyme. In experiments, it has been studied by preparing a reduced enzyme with a CO molecule bound to the binuclear center, followed by the removal of the CO molecule by a pulse or light from a laser, and monitoring changes in the optical spectrum.50 The A intermediate has a distinct optical spectrum compared with the preceding R and the subsequent Pm (Pr) states. In accordances to the calculations, the lowest energy structure is an open shell singlet. The antiferromagnetic spin coupling arises from one unpaired electron at iron and one at copper and the two bridging oxygen atoms (cf Table 1). The resulting geometry and the unpaired spin density are depicted in Figure 4. The energy difference between the open shell singlet and triplet state (10 kJ/mol) is in qualitative agreement with QM results by Yoshioka et al.14 They estimate that the energy difference between the singlet and triplet state is 30 kJ/mol, perhaps due to fact that their structure for the A intermediate is of the form Fe-O-O‚‚‚HOH rather than Fe-O-O‚‚‚Cu found here. According to present calculations, it is energetically favorable (∼10 kJ/mol) to have a water molecule bridging between the Tyr244 group and the oxygen ligand of the iron atom, in agreement with QM calculations by Yoshioka et al.14 This water molecule has hydrogen bonds to the Tyr244 group, to the oxygen ligand of the iron atom, and to a hydrogen atom on the His290 group. Therefore, in this lowest energy configuration of the A state there is a hydrogen-bonded pathway from the Tyr244 group to the dioxygen bridge. A proton could then

Kaukonen diffuse from the Tyr244 group to form the hydroxyl ligand of the copper atom in the subsequent Pm state. The enthalpy of this state is ∆H ≈ 4 kJ/mol higher compared with the R intermediate (see Figure 1 and Table 1). Zero-point energy effects, the vibrational entropy, and vdW factors are neglected here, as discussed in the section entitled Calculation of the Energetics of the Catalytic Cycle, but they are likely to push the relative free energy of this state upward. Dioxygen Bond Broken State: Pm. The results for the spin density of the Pm intermediate (Figure 1, Table 1) support the suggestion that the “fourth” electron to break the O)O double bond originates from the Tyr244 group that is covalently linked to the His240 group beside the CuB atom. Two other electrons are donated by the heme a3 and one by the CuB atom. Additionally, this tyrosine group donates a proton to the hydroxide ligand of the copper atom, yielding a neutral tyrosyl radical in the Pm state (see Table 1). The Fe)O bond length (1.64 Å) is typical for a double bond. The copper atom adopts nearly a square planar geometry. The CuB-OH bond is 1.89 Å long. The hydroxide ion has a hydrogen bond to the oxygen atom that is attached to the iron atom (the O‚‚‚H distance is 1.86 Å). The mixed nature of the unpaired spin density of the lowest energy spin state S ) 1 may explain why electron paramagnetic resonance (EPR) signals of the Pm state have not been identified.51 However, it remains unclear why also the F state is EPR silent (Figure 1).51 Positioning a water molecule in between the oxygen and OH ligands of the Cu and Fe atoms yields even lower energy structures for Pm. These types of water structures, however, are not possible in nature because a well-conserved residue, Val243, occupies the water position in space. Therefore, it can be postulated that the function of Val243 is to prevent CcO from being too exergonic in the oxidative half of the cycle. This hypothesis is in line the experimental result: mutation of the Val243 residue (in bovine numbering) to Leu, Asn, Ala, or Glu in Paracoccus Denitrificans makes the turnover rate ∼1% of the wild type, while mutations to Ile or Thr, which would have the CH3 group at the same position as in the valine thus preventing the water binding between ligands of Cu and Fe, show unchanged or ∼60% of wild type turnover rates.52,53 “Extra” Electron State: Pm + e- ) Pr? The Pr state is obtained in our QM model by adding one electron to the Pm state (Figure 1). In nature, this electron comes via heme a. Experimentally, heme a has a midpoint potential of the order of 350 mV.28 The calculated value of the midpoint potential of the binuclear center for Pm f Pr is 278 mV. This means that the electron transfer from heme a to heme a3 would NOT be feasible with the structure studied here (with the deprotonated Tyr244 in the Pr state). This result is in agreement with the suggestion by Morgan et al. (cf Figure 1 in ref 51). Thus, in the fully reduced enzyme the Pr state might have a protonated tyrosine. In an initially two-electron-reduced enzyme, however, the conversion Pm f F step could proceed by first reprotonating the tyrosine followed by electron transfer from heme a to the binuclear center (thus Pr with deprotonated Tyr244 would never form). Structurally, Pm and Pm + e- in this study are very similar. The added electron is localized on the Tyr244 group, making it a tyrosinate (Table 1). The calculated result that the Pm + e- intermediate has S ) 3/2 is not in agreement with EPR results by Morgan et al. who find the antiferromagnetic S ) 1/2 state in the EPR experiments.51 On the other hand, one must bear in mind that they

Calculated Reaction Cycle of Cytochrome c Oxidase suggest that the tyrosine is protonated in the Pr state, which is not the case here. Finally, experiments do not establish whether the Pr state is a distinct intermediate of the cycle in real physiological conditions.51 Reprotonated Tyrosine: The F State. In the F state with S ) 1.5, it is clearly energetically favorable to return the proton to the tyrosine 244. This agrees with FTIR spectroscopy observations by the Rich group, suggesting that reprotonation of the Tyr244 group occurs during the Pm f F step.41,42 The resulting geometry and spin density for F(Cu-OH- Fe)O TyrOH) near the copper and iron atoms are essentially the same as in the calculated Pr intermediate. Therefore, it is difficult to understand why experimentally the Pr state has an EPR fingerprint, whereas the Pm and F states are EPR silent.51 Taken the current calculations and the EPR data by Morgan into account would suggest that the Pr state studied here (with unprotonated Tyr244) is not the one measured experimentally.51 Oxidized Working State: Oworking. The calculated energetics suggest that the oxidized, working cytochrome c oxidase Oworking has a hydroxide ligand on both the iron and copper atoms, and that it is a S ) 1 triplet (Figure 1). It is of the order 60 kJ/mol higher in energy compared with Oresting state. The higher energy of Oworking can be explained by the mutual repulsion of the hydroxide ligands. Probably, cytochrome c oxidase is trapped in the Oworking state because an energy barrier for breaking the O-H bond in the tyrosine. Additionally a bridging water molecule with a correct orientation may be required for the Oworking f Oresting process. The fact that in the A f Pm the proton is probably taken from Tyr244 would suggest that it is the lack of the bridging water molecule between Tyr244 and the OHligands of the metals that prevents the fast transition Oworking f Oresting. The copper atom has short bonds to two of the three neighboring histidine groups (1.99 Å His290, 2.00 Å His240). The distance to His291 is longer (2.14 Å) with a geometry that is similar to the R state. The Cu-O distance in the copperhydroxide bond is 1.87 Å. This hydroxide group makes a hydrogen bond to the hydroxide ligand of the iron atom. The iron atom remains on the porphyrin plane and makes a single bond that is 1.83 Å long to its hydroxide ligand. The hydrogen bond between the hydroxide ions is a relatively good hydrogen bond with its O-H distance being 1.75 Å. Two unpaired electrons with the same spin are located on a nonbonding orbital at the iron atom and on a σ-bond between the copper and its ligands (Figure 1), resulting in a ferromagnetic coupling. By comparing the electronic density of states for Oworking and Oresting, it can be deduced that Oresting possesses a gap state associated with the Tyr244 group. This raises the total energy of the system. The calculated highest occupied molecular orbital-lowest unoccupied molecular orbital (HOMO-LUMO) gap is only 1.3 eV for Oresting. However, there is no such gap state in the Oworking intermediate. The forbidden gap is now ∼2.3 eV, but the occupied electronic states associated with the iron and copper atoms are elevated in energy due to the repulsion from the two OH- ligands. Thus, the energy cost of the single gap state in Oresting is less than the many elevated states in Oworking. The elevated states in Oworking are localized on the iron and copper atoms and their ligands. Oxidized Resting State: Oresting. Deprotonation of the Tyr244 group and transferring this proton to the hydroxide ligand of iron making it a water ligand liberates ≈60 kJ/mol ofenergy. This is the reaction Oworking f Oresting (Figure 1). The spin density at Tyr244 remains constant in this process, thus

J. Phys. Chem. B, Vol. 111, No. 43, 2007 12549 the electron is not following the proton. In experiments when no further electrons are available for the binuclear center, this decay Oworking f Oresting takes several minutes at room temperature and pH 7. In EXAFS experiments, it is found out that the three Cu-N bonds at the binuclear center of an oxidized enzyme are all 1.98 Å long.54 The fourth bond between the CuB atom and a H2O would be also 1.98 Å long; when an OH ligand was present the CuB-O bond length would be shortened to 1.90 Å.54 In this study, the high spin S ) 2 structure is the lowest energy one for the Oresting intermediate. Three of the unpaired electrons are on iron and its water ligand, and one is localized at copper, its hydroxide ligand, and three nitrogen neighbors yielding a ferromagnetic coupling. The S ) 1 triplet and S ) 0 open shell singlet states are (≈20 kJ/mol) higher in energy compared to the S ) 2 one. Experimentally, it is not clear whether this state is in low- or high-spin state as demonstrated by contradicting EPR results.2,55 The calculated CuB-N bond lengths are 2.00 (His291), 1.99 (His240), and 2.06 Å, the longest bond being again to His290. The Cu-O distance is 1.92 Å. The copper atom has an OH ligand, and the iron atom possesses a H2O ligand, but the bond lengths deviate much from typical O-H bond lengths for OH and H2O. This geometry can also be described as two hydroxide ions sharing a proton. The distances from the oxygen atoms to the shared proton are 1.41 and 1.07 Å neighboring the Cu and Fe atoms, respectively. The bond length between the copper atom and the His290 group is too long by 0.07 Å compared to the experiment. This may stem from the fact that in the calculations a hydrogen bond between the His290 and Thr309 groups is not taken into the account. It may also be a consequence of keeping the Cβ atoms in fixed positions. One Electron Reduced Binuclear Center: Eworking. The calculated atomic geometry of the Eworking state is such that the CuB atom is threefold coordinated, and without an oxygenous ligand. The calculated CuB-N bond lengths are 1.95, 1.96, and 2.14 Å, the longest bond being again to His290. The iron atom of heme a3 has a hydroxide ligand, and the Tyr244 group is protonated. According to the calculations, a H2O molecule initially positioned as a ligand of copper undergoes barrierless abstraction from the CuB atom. The spin state of this structure is S ) 1/2. The unpaired spin density is located on the iron atom and on its OH- ligand (Table 1). There is no unpaired spin density at copper, which is indicative for a formal charge of +1 at Cu. The calculated Eworking state is ≈20 kJ/mol higher in energy compared to the Eresting state with OH- ligand at copper and water on iron (Figure 1). The experimental EPR evidence by Bra¨nden et al. favors the OH ligand at iron in the Eworking state.56 This is in agreement with the current calculation. Additionally, they propose that heme a3 has low spin in state Eworking, which is also supported by the calculations (see Figure 1 and Table 1). One Electron Reduced Binuclear Center: Eresting. In the calculated Eresting state, copper has a hydroxide ligand and iron a water ligand. As with the Eworking state, the lowest energy structure has one unpaired electron (S ) 1/2). On the contrary to the Eworking state, it is now localized at copper and its OH ligand (Table 1). This suggests that the formal charge of Cu remains +2. The atomic structure is similar to the Oresting one with the difference being protonation of Tyr 244, and there are slight changes of the bonding of ligands of Cu and Fe.

12550 J. Phys. Chem. B, Vol. 111, No. 43, 2007 Experimentally, the midpoint potentials of CuB are slightly higher than for heme a3, and thus the electron should go to copper during the Oresting f Eresting step.28 This is not the case here; the electron goes to iron in the Oresting f Eresting step. On the other hand, the EPR data by Blair et al. suggests that the Eresting state has low spin ferrous heme a3 and EPR active copper (CuII, S ) 1/2), which agrees with the current calculations.57 Conclusion Consensus among the scientific community on the energy distribution over the catalytic cycle cytochrome c oxidase has not yet been reached; hence, no unique model for the function of the system exists. In this work, it is shown that it is possible, using a method based on quantum mechanics, to construct a model that fits well with the experimental energetic data of Wiksto¨m’s group.8 Atomic structures and their electronic configurations are suggested based on the QM calculations. In particular, it is found that de/reprotonation of Tyr244 is energetically plausible. Interestingly, it is found that two pairs of intermediates (Pm, F, and Oresting, Eresting) differ in geometry approximately only by the protonation state of Tyr244. Acknowledgment. The generous computer resources of the Center of the Scientific Computing, Espoo, Finland are acknowledged. The author is grateful to Michael Verkhovsky, Marten Wikstro¨m, Anne Tuukkanen, Liisa Laakkonen, Camilla Ribacka, Audrius Jasaitis, Mika Molin, Dage Sundholm, Mikael Johansson, and Christopher Latham for valuable discussions. Supporting Information Available: The coordinates of the QM structures labeled in Table 1 given in xyz and TURBOMOLE formats. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Brezezinski, P.; G., L. Biochim. Biophys. Acta 2003, 1605, 1-13. (2) Malmstro¨m, B. Chem. ReV. 1990, 90, 1247-1260. (3) Iwata, S.; Ostermeier, C.; Ludwig, B.; Michel, H. Nature 1995, 376, 660-669. (4) Tsukihara, T.; Aoyama, H.; Yamashita, E.; Tomizaki, T.; Yamaguchi, H.; Shinzawa-Itoh, K.; Nakashima, R.; Yaono, R. Science 1996, 272, 1136-1144. (5) Svensson-Ek, M.; Abramson, J.; Larsson, G.; To¨rnroth, S.; Brezezinski, P.; Iwata, S. J. Mol. Biol. 2002, 321, 329-339. (6) Tsukihara, T.; Shimokata, K.; Katayama, Y.; Shimada, H.; Muramoto, K.; Aoyama, H.; Mochizuki, M.; Shinzawa-Itoh, K.; Yamashita, E.; Yao, M.; Ishimura, Y.; Yoshikawa, S. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 15304-15309. (7) Moody, A. Biochim. Biophys. Acta 1996, 1276, 6-20. (8) Bloch, D.; Belevich, I.; Jasaitis, A.; Ribacka, C.; Puustinen, A.; Verkhovsky, M.; Wikstro¨m, M. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 529-533. (9) Blomberg, M. R. A.; Siegbahn, P. E. M.; Babcock, G. T.; Wikstro¨m, M. J. Inorg. Biochem. 2000, 80, 261-269. (10) Blomberg, M. R. A.; Siegbahn, P. E. M.; Babcock, G. T.; Wikstro¨m, M. J. Am. Chem. Soc. 2000, 122, 12848-12858. (11) Blomberg, M. R. A.; Siegbahn, P. E. M. Inorg. Chem. 2003, 42, 5231-5243. (12) Moore, D. B.; Martı´nez. J. Phys. Chem. A 2000, 104, 2367-2374. (13) Siegbahn, P. E. M.; Blomberg, M. R. A.; Blomberg, M. L. J. Phys. Chem. B 2003, 107, 400-417. (14) Yoshioka, Y.; Kawai, H.; Yamaguchi, K. Chem. Phys. Lett. 2003, 374, 45-52. (15) Cukier, R. I. Biochim. Biophys. Acta 2004, 1656, 189-202. (16) Popovic´, D. M.; Stuchebrukhov, A. A. J. Am. Chem. Soc. 2004, 126, 1858-1871. (17) Popovic´, D. M.; Quenneville, J.; Stuchebrukhov, A. A. J. Phys. Chem. B 2005, 109, 3616-3626. (18) Fadda, E.; Chakrabarti, N.; Pome´s. J. Phys. Chem. B 2005, 109, 22629- 22640. (19) Olsson, M. H. M.; Sharma, P. K.; Warshel, A. FEBS Lett. 2005, 579, 2026-2034. (20) Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346-354.

Kaukonen (21) Payne, M. C.; Teter, M. P.; Allan, D. C.; Arias, T. A. ReV. Mod. Phys. 1992, 64, 1045-1097. (22) Curtiss, L. A.; Raghavachari, K.; Redfern, R. C.; Pople, J. A. J. Chem. Phys. 2000, 112, 7374-7383. (23) von Arnim, M.; Ahlrichs, R. J. Comp. Chem. 1998, 19, 17461757. (24) Scha¨fer, A.; Huber, C.; Ahlrichs, R. J. Chem. Phys. 1994, 100, 5829-5835. (25) Klamt, A.; Schu¨u¨rmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 5, 799-805. (26) Kuhn, B.; Kollman, P. A.; Stahl, M. J. Comput. Chem. 2004, 25, 1865-1872. (27) Schutz, C. N.; Warshel, A. Proteins: Struct., Funct., Getet. 2001, 44, 400-417. (28) Gorbikova, E. A.; Vuorilehto, K.; Wikstro¨m, M.; Verkhovsky, M. I. Biochemistry 2006, 45, 5641-5649. (29) Blomberg, M. R. A.; Siegbahn, P. E. M. J. Comput. Chem. 2006, 27, 1373-1384. (30) Schroedl, N. A.; Hartzell, C. R. Biochemistry 1977, 16, 49664971. (31) Reiss, H.; Heller, A. J. Phys. Chem. 1985, 89, 4207-4213. (32) Chang, R. Physical Chemistry for Chemical and Biological Sciences; University Science Books: Sausalito, CA, 2000. (33) Palascak, M. W.; Shields, G. C. J. Phys. Chem. A 2004, 108, 36923694. (34) Liptak, M. D.; Shields, G. C. J. Am. Chem. Soc. 2001, 123, 73147319. (35) NIST Chemistry WebBook. http://webbook.nist.gov/chemistry/ (date accessed May, 2007). (36) CRC Handbook of Chemistry and Physics, 85th ed.; Lide, D. R., Ed.; CRC Press: Boca Raton, FL, 2004. (37) Fischer, S.; Smith, J. C.; Verma, C. S. J. Phys. Chem. B 2001, 105, 8050-8055. (38) Denisov, V.; Venu, K.; Peters, J.; Ho¨rlein, H. D.; Halle, B. J. Phys. Chem. B 1997, 101, 9380-9389. (39) Olano, L. R.; Rick, S. W. J. Am. Chem. Soc. 2004, 126, 79918000. (40) Wirstam, M.; Lippard, S. J.; Friesner, R. A. J. Am. Chem. Soc. 2003, 125, 3980-3987. (41) Iwaki, M.; Puustinen, A.; Wikstro¨m, M.; Rich, P. R. Biochemistry 2003, 42, 8809-8817. (42) Iwaki, M.; Puustinen, A.; Wikstro¨m, M.; Rich, P. R. Biochemistry 2004, 43, 14370-14378. (43) Belevich, I.; Tuukkanen, A.; Wikstro¨m, M.; Verkhovsky, M. Biochemistry 2006, 45, 4000-4006. (44) Wikstro¨m, M.; Jasaitis, A.; Backgren, C.; Puustinen, A.; Verkhovsky, M. Biochim. Biophys. Acta 2000, 1459, 514-520. (45) Babcock, G.; Callahan, P.; Ondrias, M.; Salmeen, I. Biochemistry 1981, 20, 959-966. (46) Carter, K.; Palmer, G. J. Biol. Chem. 1982, 257, 13507-13514. (47) Vanneste, W. Biochemistry 1966, 5, 838-848. (48) Rovira, C.; Kunc, C.; Hutter, J.; Ballone, P.; Parrinello, M. J. Phys. Chem. A 1997, 101, 8914-8925. (49) Ralle, M.; Verkhovskaja, M. L.; Morgan, J. E.; Verkhovsky, M. I.; Wikstro¨m, M.; Blackburn, N. J. Biochemistry 1999, 38, 7185-7194. (50) Verkhovsky, M. I.; Morgan, J. E.; Wikstro¨m, M. Biochemistry 1994, 33, 3079-3086. (51) Morgan, J.; Verkhovsky, M.; Palmer, G.; Wikstro¨m, M. Biochemistry 2001, 40, 6882-6892. (52) Riistama, S.; Puustinen, A.; Verkhovsky, M.; Morgan, J. E.; Wikstro¨m, M. Biochemistry 2000, 39, 6365-6372. (53) Riistama, R. Structural and functional studies of bacterial hemecopper oxidases. http://ethesis.helsinki.fi (date acessed May, 2007). (54) Fann, Y. C.; Ahmed, I.; Blackburn, N. J.; Boswell, J. S.; Verkhovskaja, M. L.; Hoffman, B. M.; Wikstro¨m, M. Biochemistry 1995, 34, 10245-10255. (55) Moody, A.; Cooper, C.; Gennis, R.; Rumbley, J.; Rich, P. Biochemistry 1995, 34, 6838-6846. (56) Bra¨nde´n, M.; Namslauer, A.; Hansson, O ¨ .; Aasa, R.; Brzezinski, P. Biochemistry 2003, 42, 13178-13184. (57) Blair, D.; Witt, S.; Chan, S. J. Am. Chem. Soc. 1985, 107, 73897399.