Linking chemical electron-proton transfer to proton pumping in cytochrome c oxidase: broken-symmetry DFT exploration of intermediates along the catalytic reaction pathway of the iron-copper dinuclear complex.

After a summary of the problem of coupling electron and proton transfer to proton pumping in cytochrome c oxidase, we present the results of our earlier and recent density functional theory calculations for the dinuclear Fe-a3-CuB reaction center in this enzyme. A specific catalytic reaction wheel diagram is constructed from the calculations, based on the structures and relative energies of the intermediate states of the reaction cycle. A larger family of tautomers/protonation states is generated compared to our earlier work, and a new lowest-energy pathway is proposed. The entire reaction cycle is calculated for the new smaller model (about 185-190 atoms), and two selected arcs of the wheel are chosen for calculations using a larger model (about 205 atoms). We compare the structural and redox energetics and protonation calculations with available experimental data. The reaction cycle map that we have built is positioned for further improvement and testing against experiment.

Introduction

The question of how to effectively activate O2 for reaction is of great interest in biochemistry and organic and inorganic chemistry (including materials science) because the products and byproducts of O2 reactions can be variable and also because O2 often encounters kinetic barriers to the reaction, “the O2 overpotential”. Further, coupled electron–proton transfers, both sequential and concerted, are a major theme in the reaction chemistry of redox-active systems, including many biologically important metalloenzymes. Proton-pumping redox enzymes are functionally essential for energy storage but diverse in their mechanisms and structures, comparing cytochrome c oxidases (CcOs, complex IV) to cytochrome bc1 centers (complex III), for example.

Two main classes of mechanisms have been proposed for proton pumping in CcOs: (1) a direct mechanism where the chemistry of O2 reduction is directly coupled to proton pumping within the O2 binding dinuclear complex (DNC) containing an Fe-heme and Cu sites (Fe-a3 and CuB); (2) an indirect mechanism where conformational or electrostatic interactions with the surrounding protein, in particular with the neighboring redox donor centers in the enzyme, are required. Our current and past work has centered on potential direct mechanisms for both fundamental and practical reasons. The driving force for proton pumping resides in the metal redox and O2 binding/reduction chemistry of the DNC, so this is clearly a good place to start. By analyzing what the DNC can do with the other relevant sites as simple inputs, we want to identify the distinctive positive features as well as the potential limitations of our models. Our viewpoint is modular, which simplifies analysis where separability is valid, and reveals its limitations, where predictions prove questionable or false. Our work is an illustration of the theme that while “the map is not the territory”,[1] to understand the territory, it is best to build and improve the map.

Following up on our earlier work,[2] we explore how different metal and tyrosine redox and electronic states interact with molecular O2 binding to promote proton transfer and proton pumping in the DNC of CcO. The theoretical/computational tools are spin-polarized broken-symmetry density functional theory (BS-DFT) in conjunction with continuum solvation dielectric modeling for the surrounding protein–solvent environment.

Biochemical Background

Mitochondrial and bacterial respiration are the major sources of functional energy via adenosine triphosphate (ATP) synthesis in aerobic organisms.[3−7] CcOs are the terminal electron acceptors in the respiratory chain of mitochondria and many bacteria.[3−5] Essentially, O2 linked to the metal center is the final sink for electrons from the earlier electron-transfer complexes I–III and mobile carriers [ubiquinone (Q) and cytochrome c]

There are a number of critical functional requirements for the operation of eq 1 within the dinuclear heme iron–copper complex of CcO. Electrons and protons must be taken up at various steps of the catalytic cycle, O2 must be efficiently bound and activated for reaction without releasing oxygen-radical intermediates like superoxide anion O2•– or other toxic products [for example, hydrogen peroxide (H2O2)], and protons must be transported to and react with metal-bound O2 and/or its subsequent intermediates. Equation 1 summarizes the relevant net chemical reaction but not the mechanism or the essential coupling to proton pumping.

Proton pumps like CcO are enzymes coupled to energy transduction machines.[6,7] Beyond the chemistry of eq 1, the free energy of this reaction in CcO has to be effectively utilized for proton pumping from the inside (in) to the outside (out) of the membrane against the electrochemical gradient (Figure 1). Figure 1 is a modified and annotated version of Figure 1 in ref (2) representing the structure of the bacterial B-type enzyme from Thermus thermophilus (Tt). Important inputs into the dinuclear Fe–Cu reaction center (DNC Cu is CuB and the corresponding Fe is Fe-a3) are shown. The dinuclear Cu–Cu center called CuA is shown by two small blue dots. Heme-b, which is the immediate donor of electrons to the DNC Fe–Cu center, is not shown; heme-b lies behind heme-a3 in this orientation. Many spectroscopic and kinetic studies show that electrons are donated from CuA to heme-b Fe and then to the DNC heme-a3. Further, O2 binding is ordered compared to electron (e–) transfer, so that typically all four metal redox centers are 1e–-reduced before O2 binds to the DNC reaction site between Fe-a3 and CuB.

Figure 1

Schematic (left) of cytochrome ba3 with noted functionalities and cartoon of the same (right) showing the O2 channel inside the lipid bilayer, the reaction chamber between Cu and Fe, and the electron- and proton-transfer “pipes”. Heme-b is not shown.

In standard bioenergetics language, the inside is the “n” or negative side (Figure 1, bottom, cytoplasm) and the outside is the “p” or positive side of the membrane (Figure 1, top, periplasm) corresponding to the polarity of the membrane electrochemical gradient. The detailed mechanism of proton pumping is far from evident, but the overall reaction stoichiometry isover a full reaction cycle. Some language is useful for bookkeeping: of the (4 + n)H+in entering the reaction chamber, 4H+, called “scalar protons”, react with O2 to produce 2H2O, while the other nH+, called “vector protons”, are pumped across the membrane. Measurements of n have yielded values near n = 4 for A-type CcOs, but lower values near n = 2 have been measured for B- and C-type CcOs. (For example, extrapolating the charge-transfer measurements for the oxidative half of the reaction cycle to the entire oxidative plus reductive reaction cycle for B-type Tt ba3 yields n = 2.5 approximately.[8]) Finding accurate experimental values for n is difficult, and significant uncertainties remain.

In common with respiratory chain complexes I and III[5−7] and ordered sequentially by increasing redox potentials, complex IV links e– transfer down the electron-transport chain to proton pumping across the inner membrane of the mitochondria and, analogously, across the plasma membrane in most aerobic bacteria. The resultant proton electrochemical gradient, or proton motive force, is utilized for ATP synthesis from adenosine diphosphatase (ADP) plus inorganic phosphate by complex V (ATP synthase), which also resides within the same membrane. The proton electrochemical gradient has a higher positive potential on the outside (p) of the membrane and a lower negative potential on the inside (n) of the membrane (matrix side in mitochondria). The proton leak “downhill” across the membrane is largely funneled through the ATP synthase active site to generate ATP from ADP+Pi. As expected from energy conservation, continued proton pumping is required to maintain the membrane potential against the leak and to provide the driving force for ATP synthesis. Alternatively, if the electrochemical gradient is depleted, proton pumping is required to restore the gradient. The electrochemical potential is the sum of an electrical potential (voltage) difference, ΔΨ, and a pH difference, ΔpH, across the membranewhere F is the Faraday constant.

While the reaction inputs and outputs of the respiratory complexes are understood in broad stroke, in no case is the actual coupled electron redox/proton pumping mechanism understood in even reasonable detail. Many features are known from spectroscopy and kinetics about the chemical mechanism of O2 reduction and subsequent mechanistic steps of the reaction cycle, but there is no current resolution about whether or how proton pumping is directly coupled to the chemistry of O2 reduction at the Fe–Cu DNC or indirectly coupled through conformational or electrostatic changes in the surrounding protein, particularly involving the nearby redox centers (the CuA homodinuclear and mononuclear Fe-heme donor centers).[9−11] It is possible that the different enzyme classes A1 versus B could have significant differences in their mechanisms (see the later discussion), including whether a direct or indirect mechanism (or some mixture) applies.

X-ray crystallographic structures are now available for mitochondrial [bovine heart (Mt)] and several bacterial CcOs, including Paracoccus denitrificans (Pd), Rhodobacter spheroides (Rs; all A1 type), and T. thermophilus (Tt; B type)[3,12−16] The Mt structure is by far the most complex, with 13 subunits (3 encoded by mitochondria and 10 by the nucleus). However, the three mitochondrial subunits (I–III) of Mt closely resemble the corresponding subunits in the bacterial enzymes Pd and Rs, and all are of class aa3. The Tt enzyme is of a different class (namely, ba3) but still has many similarities in sequence, structure, and mechanism to the other CcOs. The bacterial enzymes are easier to work with than the more complex Mt protein and have, therefore, been subjected to studies using systematic mutagenesis of different sites. Our main focus will be on analysis of the Tt enzyme, for which Fee, Stout, and co-workers have obtained high-resolution X-ray structures, including several different mutants, and measured kinetic spectroscopic parameters of many intermediates over the redox and catalytic cycle.[12,13,17−19] Kinetic analysis is enhanced in ba3-type enzymes including Tt because the absorption spectrum of heme-b is well-separated from that of heme-a3 in contrast to aa3 (A type) enzymes, where the heme-a and -a3 bands overlap. Also, importantly, the Tt enzyme has a single input path (KB path) for protons into the dinuclear Fe–Cu center, which simplifies analysis of the proton pathways compared to two paths for A1-type enzymes. While the KB path occupies a structurally analogous region compared to the K path in A1 enzymes, the KB path does not contain lysine, in contrast to the K path in A1.

The biological e–-transfer sequence proceeds[17] in one electron (1e–) steps from cytochrome c to the dinuclear mixed-valence CuA site in subunit II (having a 2Cu1.5+ center in the oxidized state and 2Cu+ in the reduced), to a mononuclear Fe-heme redox center in subunit I (which in Tt is heme-b but in Pd, Rs, and Mt is heme-a), and finally to the dinuclear heme-a3–Fe–CuB reaction center, which is also in subunit I and is where O2 binds. The proton entry pathway is simpler in the Tt ba3 enzyme, having a single entry pathway called the “KB path”, compared to two potential entry pathways, called the “K and D paths”, in the class aa3 enzymes. O2 enters via a well-defined hydrophobic pocket near the dinuclear Fe-a3-CuB site. Figure 1 (left) shows an overall structural model of CcO from Tt with the inputs and outputs superposed, and Figure 1 (right) is a cartoon of the reaction chamber with reactants and products. In Tt ba3, the subunit I residue Ile235 is present (see Figure 1) in place of a glutamate residue present in A1-type enzymes, so that the second proton-transfer pathway (the D or Q proton-transfer pathways seen in A1-type enzymes) is absent in B-type enzymes. This also eliminates the corresponding path to the “proton loading site” as proposed in A-type enzymes. The significant implication is that the path into the DNC is common to both the scalar and vector protons. Any subsequent bifurcation of proton flow between the scalar and vector pathways occurs at or after the entry point into the DNC.

Mitochondrial CcOs play an essential role in human health. Adequate ATP supplies are necessary for all important metabolic cellular functions.[4] Also, disruptions in O2 binding and electron- or proton-transfer reactions at CcO can lead to an excess of reactive oxygen species including hydroxyl and superoxide radicals. These oxygen species can arise directly at CcO or can be induced earlier in the electron-transport chain, for example at complexes I and III, because of a blockage in electron or proton transport at CcO. Finally, disruptions in the normal proton-transport pathways may also have damaging effects on the membrane potential and acidity in the matrix or periplasmic space. Understanding the structures, mechanisms, and functions of mitochondrial CcO is therefore important for better analysis of genetic and metabolic diseases and cancer and is also relevant to pathologies in aging.[3−7,20,21] Mitochondrial DNA is intrinsically subject to elevated mutation frequency possibly because of its proximity to leaky electron-transport chains that produce reactive oxygen species. The close functional and mechanistic relationship between mitochondrial CcO and various bacterial CcOs allows for the development of an overview of the catalytic proton-pumping cycle, while variations in the structures and energetics/kinetics can provide insight into the interactions between different redox centers and proton sources.

The A-type enzymes in Pd and Rs reside in α proteobacteria, some of which were probably ancestors of eukaryotic mitochondria. B-type enzymes, while simpler in some ways, also show significant variations in different species/genera. Further, mutagenesis studies in bacteria provide a good way of examining blockage and inhibition of proton and electron transfer and O2 binding. There have also been very significant accomplishments in making synthetic analogues of the DNC, and testing the reaction, physical, and spectroscopic properties of these systems.[22−24]

Electron Flow within CcO

Electron flow in both A1- and B-type CcOs has been studied extensively. Electrons transfer from cytochrome c to CuA, then on to heme-a (or heme-b in A1 and B types, respectively), and from there to Fe-a3–CuB.[8,25] After O2 binding to the DNC, a number of oxygen-derived intermediates are formed. A good deal is known about the timing of O2 binding and e– transfer by binding CO to heme-a3 and then introducing O2, followed by photolyzing the CO in the fully reduced or half-reduced (mixed-valent) forms of the bovine and related bacterial enzymes. However, CO photolysis experiments have some difficulties, as discussed below.

Reaction Cycle of the DNC in CcO

Because of their considerable structural and functional similarities, A- and B-type enzymes probably have a common mechanism of O2 reduction. Using optical spectroscopy to monitor transient kinetics, Siletsky et al.[8] and Szundi et al.[26] observed four intermediates for reduction and O2 binding and reactions in Tt cytochrome ba3 at the DNC, as shown in Figure 2. To obtain proper timings for kinetics, the reaction of the fully reduced enzyme is started by photolysis of the reduced Fe-a3–CO adduct of the enzyme in the presence of O2, and then changes in a wide array of properties are followed from the nanosecond to second time scale. Here the fully reduced enzyme includes the fully reduced state (R) of the DNC and the reduced states of heme-b and CuA.

Figure 2

Catalytic “wheel” of cytochrome ba3 focusing on the DNC. Concentration of O2 for R → A is 90 μM for τ ≈ 9 μs, from Szundi et al.[26]

Compounds A, P (where P is likely composed of both P and P), and O in Tt ba3 are similar to their analogues in A1-type enzymes. More broadly, the optical absorption difference spectra of heme-a3 do not clearly report on the oxidation state (or related chemical/electronic features) of either CuB or the unusual tyrosine called Y′OH. Recently, Szundi et al.[26] found that the early oxidative steps of the reaction become unreliable with the CO photodissociation method because CO can then bind to CuB, and this inhibits O2 binding to Fe-a3 and CuB. Instead of using CO bound to the reduced state R, they instead utilized a photolabile O2 carrier complex to deliver O2 to the DNC. The kinetics are then clarified by examining different concentrations of O2 delivered to the DNC. Lifetimes (τ = inverse of rate constants r) from this study are assembled in Figure 2. These lifetimes differ from those we previously presented based on Siletsky’s kinetic studies.[2] Starting from state R, O2 binding occurs very rapidly (even at the low effective [O2] = 90 μM) and so does the following formation of ferryloxo states, spectroscopic state P. State P is a combination of electronic/protonation states P and P, and because the heme-b Fe is also reduced, the very fast e– transfer to reduce the tyrosine radical is not detectable in these experiments. (The state P is detected in separate experiments, where the O2 reaction is started from the “mixed-valence state”, with the DNC being in state R, while the heme-b and CuA sites are oxidized. Then the cycle stops at P because of a lack of reducing electrons.) The cycle depicted in Figure 2 omits the next spectroscopically observed intermediate, P, which corresponds to the timing for e– transfer from reduced CuA to oxidized heme-b Fe, which was the redox center oxidized in the P → P transition. After the heme-b Fe is reduced, it is ready to deliver the final electron in the transition P → O, which is comparatively slow. The reaction cycle experiments of Szundi et al. end with state O.[26] Recent pulse radiolysis experiments[27] have shown a comparatively slow 1e–-transfer rate from state O (reducing agent heme-b Fe) to heme-a3 Fe (III → II), lifetime τ = 4800 μs, compared to faster e– transfer from CuA to heme-b, τ = 260 μs in Tt ba3. Therefore, both P → O and O → R have comparatively slow e– transfers. We have argued in our published work[2] that e– transfer from heme-b (reduced) to the DNC may be determined by the nature of the oxidant center present on heme-a3, which is likely to involve both reorganization energy terms and the driving force (redox potentials) for different reaction steps.

For future reference, we comment briefly on the issue of the electrochemical potential (load versus no or low load) in kinetic studies and proton-pumping experiments.[7] When experiments are conducted in detergents and/or buffers, there is no net directionality (vectorial character) to the orientation of the CcO proteins in the media, and so there is no net electrochemical potential. For kinetic studies, these types of conditions are often used,[17,26,27] and they lead to faster kinetics than those in the presence of an electrochemical potential (experiments under load). Even X-ray crystallographic experiments done in the lipidic cubic phase have no net directionality (see Tiefenbrunn et al.[19]) because the individual CcO proteins are aligned antiparallel in the lipidic cubic phase. By contrast, experiments designed to measure the proton-pumping stoichiometry or net charge (q; e–,H+) transfer across the membrane (usually vesicle experiments) are set up so that the CcO proteins have specific directional alignment with respect to the inside versus outside of the vesicles. Proton-pumping stoichiometry (vector protons) is typically determined by running a full or partial cycle and collapsing the electrical potential difference across the vesicle using an ionophore, which does not conduct protons but does conduct some other charged particles. (For example, valinomycin allows countertransport of K+ ions when protons are pumped in the CcO reaction cycle.[7,28]) In that case, ΔpH is measured in the corresponding partial cycle. To measure full charge transfer, corresponding to electrons, scalar protons, and vector protons driven against the existing gradient, measurements of the electrical potential change over a partial cycle are obtained. Both for proton pumping and net charge transfer, the standing electrochemical potential (and its electrical component) is typically much less than that under actual physiological conditions. (Usually, the electrical term dominates eq 3 under physiological conditions in mitochodria and analogously in aerobic bacteria.[7] So, in practice, most proton-pumping and charge-stoichiometry experiments are under low load.

Our DFT calculations for the reaction cycle are performed under idealized pH gradient conditions, where we assume that each pumping of a vector proton across the DNC corresponds to transfer from a neutral pH = 7 cytoplasm (n side) to an acidic pH = 3 periplasm (p side), and there is no electrical potential gradient, as in the vesicle experiments with an ionophore described above. This makes bookkeeping for the cost of proton pumping easy because only the vector protons have a cost, which is accounted for as the proton exits from the DNC. By contrast, a standing electrical potential across the membrane will lead to a pump cost for the vector protons and to a cost for charge transfer for the electrons and scalar protons because these combine in the reaction chamber of the DNC. The energy accounting here is more complex, and we will defer it to future work. (For a simple analysis, see Siletsky et al.[8]) However, there is one general observation that is very useful. If the same energy cost over an entire cycle (say for pumping four vector protons from pH = 7 to 3) is instead distributed over both the vector proton transfers and the scalar proton plus electron transfers as occurs for an electrical potential difference, rather than a pH difference, the charge-transfer cost across the DNC over individual steps is considerably less (roughly 1/2 or less because there are q = 8 charges transported across the membrane versus n = 4 vector protons).

Objectives of the Theoretical/Computational Cycle

In our published work (2008),[2] we have introduced the first complete catalytic reaction cycle “wheel” for the cytochrome ba3 DNC, based on a detailed chemical model supported by large-scale DFT calculations. Our work is anchored in several known features of the experimental reaction scheme but goes well beyond what is presently established. Very recently, we have further refined the DFT model by examining alternative tautomeric and protonation states.

We emphasize some important general points starting from Figure 2.[8,29](1) Electron and Proton Delivery. The points of electron entry into the cycle as well as O2 binding are well established, but only four protons are accounted for within this scheme. These are the scalar protons as in eq 1, and the vector protons are not included (see eq 2). On the basis of eq 2, the optimal proton-pumping ratio is 4H+/4e–, equivalent to 1H+/1e–. Typically, experimental proton-pumping ratios are near 1, but early published measurements with ba3 yielded values of about 0.4–0.5. In any event, a pumping ratio of 1 requires a proton uptake ratio for H+in/1e– of (4 + n)/4 = 2. Any uptake ratio greater than 1 requires special chemistry. Moreover, under physiological conditions (and in some experimental setups), the pumped protons must be moved against the electrochemical gradient at a cost in ΔG°, using the chemical driving force available from the redox/protonation reactions with O2. (2) Important Chemical Intermediates. There are critical and distinctive chemical questions raised by the catalytic reaction cycle and the structure/composition of the DNC. (2a) At least one Fe–Cu bridging peroxo intermediate is likely to form subsequent to the initial oxy complex (A), which has been identified, but the peroxo intermediate(s) has not been spectroscopically identified. The peroxo complex is likely to be very short-lived[26] but is clearly essential to the energetics of O–O bond cleavage and for avoidance of a substantial overpotential. Both the energetics and properties of these intermediates are of great interest and are connected to binding modes and protonation states. (2b) There is a unique tyrosine (Y′OH in Figure 2) that is formed from a covalent bond between the tyrosine phenol and a histidine imidazole ring (residues Y237 and H233 in Tt ba3). Further, H233 is coordinated to CuB. What are the roles of this unusual linkage in the catalytic cycle? (3) Proton Pumping by the DNC. What is the physical mechanism connecting the redox and chemical bonding cycle to the proton pumping? The lack of mobile proton carriers within CcO precludes a “redox loop” mechanism, but CcO pumps protons, nonetheless. We have already proposed a proton exit channel into the watery surface between subunits I and II (see Figures 1 and 3). Specifically, this involves deprotonation of H376 (from the positive imidazolium form). As an alternative, deprotonation may instead involve a nearby bound water (probably a hydronium ion), which also looks reasonable. Our current working hypothesis is that the DNC complex itself and the KB path conducting protons into the DNC are the major factors driving proton pumping. We also think that the redox potentials and timing of e– transfer from heme-b Fe to the DNC play a role, as will be described later.

Figure 3

Polar representation of the new 14-step catalytic mechanism. The central polar plot of the running sum, ∑ΔG°, shows the energetic course of catalysis as computed at the level of PW91.

A combination of electrostatic and large-scale DFT calculations with associated energy analysis has been used to examine the proton-pumping mechanism in class A1 (aa3-type) enzymes.[11,30,31] In a number of papers, Blomberg and Siegbahn[32−34] have proposed that 1e– transfer from the CuA center to the mononuclear heme-a drives protonation at a pump-loading site (PLS) close to a propionate in heme-a3, as well as protonation of the glutamate (Glu286 in bacterial Rs and a homologous glutamate in mitochondria). Both protons originate from the D proton channel. The Glu286 proton then activates proton exit from the PLS after 1e– transfer from heme-a to the DNC Fe site in heme-a3, with a coupled “chemical” proton transfer close to the reduced Fe-a32+ site (the 1e–-transfer step modeled is O → E, the first 1e– transfer in the O → R 2e– reduction). The chemical proton transfer is also modeled as originating from the D channel. Presumably, at other steps of the catalytic reduction O2 binding reaction cycle, some “chemical” protons originate from the K channel. In any event, in class B (ba3-type) enzymes, there is no D-type proton-transfer pathway, and in Tt ba3, there is an isoleucine (Ile235) in a position homologous to Glu286 (see Figure 3 in ref (35)). Nonetheless, ba3 enzymes do pump protons, although probably at lower efficiency than aa3.[8](4) Proton Entry and Flow. Therefore, we will focus here on the catalytic reaction cycle, including proton flow into the DNC from the KB path and particularly within the DNC itself. The vector and scalar protons almost certainly take different paths through the DNC. The currently available high-resolution X-ray structures of Tt cytochrome ba3[35,36] from Fee, Stout and collaborators[19] show a much richer array of waters bound to the DNC than that seen in earlier X-ray structures at lower resolution.

Models and Methods: Catalytic Wheel

(a) In our earlier paper (2008),[2] we constructed the entire catalytic reaction “wheel” for the cytochrome ba3 DNC, employing large-scale DFT calculations. (b) On the basis of subsequent work, we have found an improved energy path by modifying the sequence of proton uptake and proton-transfer events and examining alternative proton tautomeric states. These two sets of models each contain 183–189 atoms and are called S1 and S2, for the original and for the newer model encompassing different tautomers for several states. Both models are based on the Hunsicker-Wang and Fee X-ray structure of ba3.[12] Both of the smaller models are computed with a Perdew–Wang (PW91) exchange-correlation potential for geometry optimization in a dielectric medium with dielectric constant ε = 18.5 for the surrounding environment. Then single-point PW91 and B3LYP* calculations are compared over the reaction path. (c) For selected parts of the catalytic reaction path, we have now introduced a larger model (L), incorporating additional waters and residues and with greater flexibility in the electronic states described. The larger model is built based on the high-resolution (1.8 Å) X-ray structure of Tiefenbrunn et al.[19] in the lipidic cubic phase, which has a higher density of water in the protein interior. As discussed in a recent paper,[37] the observed X-ray structure contains a bound dioxygen species that is probably hydroperoxide. As in the Hunsicker-Wang and Fee (2.3 Å resolution) X-ray structure,[12] the DNC in the protein is probably radiolytically reduced, but within DFT, we can control the oxidation state and the types of bound ligands that are introduced in the calculations. The model L calculations utilized an OLYP exchange-correlation potential. Because of the increased cluster size, we found that the calculations could be performed with fewer geometric constraints. In particular, only the link atoms at the outside boundary of the quantum cluster were constrained to their positions in the Tiefenbrunn X-ray structure, as specified in ref. (37). The dielectric model is the same as models S1 and S2, and all use the COSMO solvation model.[38−41] All calculations use the ADF codes.[42,43]

Catalytic Reaction Cycle Results and Discussion

Small Models

The catalytic wheel diagram for the DNC from our new model S2 is presented in Figure 3. From geometry optimization and energy evaluation of 30 different tautomeric/net protonation states (defined as model Sall) for the 14-step cycle, we have selected 14 states in model S2. Compared to model S1,[2] if a different tautomer is chosen in S2, this is indicated by a subscript on the state number, for example, 7 versus 7. Where a proton has been added, this is indicated by an H, so; for example, 2H compared to 2 and 8H compared to 8 indicate both a tautomeric shift and an additional proton added.

Figure 3 shows a detailed map of the catalytic reaction cycle for model S2. See also the related cycle in ref (2). We have augmented the figure by labeling intermediates O, R, A, P, and P, which have all been observed spectroscopically. All of these intermediates have approximate counterparts in Figure 2, and the theoretical scheme is supported by detailed DFT calculations. Most of the spectroscopic intermediates correspond to more than one detailed intermediate in the 14-state reaction cycle. For example, states 13, 14, and 1 correspond to experimental intermediate O; all of these states are FeIII–CuII species, but they differ with respect to protonation state(s) or H2O binding. In addition, our catalytic wheel includes “true peroxo” states not seen spectroscopically, which are testable by energetic calculations and potentially by spectroscopic measurements when these (or related states) are isolated. Structures 3, 6, 10, and 13 are detailed 2D maps of the dinuclear Fe-a3–CuB center, while the other states are represented by simple chemical formulas. A full 3D picture of the Fe3+–(O–O)2––Cu2+ bridging peroxo state, intermediate 5, is shown in Figure 4. By comparison with Figure 3, the coordination, covalent, and hydrogen bonds are quite clear in the 2D maps, but the 3D links are distorted; for example, the heme proprionate side chains are attached to the Fe-heme edge.

Figure 4

Active site structure of compound 5 of ba3. Left: smaller model used in the 2008 calculations.[2] Right: larger model used in recent calculations,[37] which is taken from the DNC of the high-resolution (1.8 Å) X-ray crystal structures (PDB entry: 3S8G) of ba3 CcO from Tt.[19] For ease of visualization, these two models are not on the same scale.

The new 14-step cycle differs from the previous one[2] in (1) having protons enter the reaction path earlier at two steps, 1 → 2H and 7 → 8H, (2) retaining water in state 11 compared to state 11, and (3) utilizing different proton tautomers, for all of the states where there are subscripts in Figure 3. (The original tautomers from model S1 are not subscripted.) By convention, tautomers from model S1 are each tautomer 1; for example, 6 = 6. In constructing the intermediates in Figure 3, we mainly stayed with the lowest-energy tautomers. The one exception is state 6, which is critical to the peroxide cleavage reaction. In our model S2, this is 6 kcal mol–1 higher (in PW91) than the lowest tautomer, 6, where the special tyrosine (Y237) is protonated in place of the peroxide. As we will show later, both tautomers 6 and 6 are energetically feasible along the reaction path. The outcome of replacing our original reaction path (model S1) with one constructed from the lowest-energy tautomers (also alternate protonation states) is shown in Figure 5, comparing parts A and B. The effect is striking. These changes to the previous mechanism result in a substantially “smoother” catalytic cycle. The largest positive free-energy difference ΔG°, accounting for successive uphill steps, is lowered to about 20 kcal mol–1 for PW91 in model S2 (31 kcal mol–1 in B3LYP*), while in the reaction pathway of model S1 for both PW91 and B3LYP*, this value is close to 30 kcal mol–1. It is easy to see that both exchange-correlation potentials run in a similar manner for energies over the reaction cycle and that Figure 5B shows a more reasonable cycle for both PW91 and B3LYP*. The ΔG° values for individual steps are considerably less in all cases (about 15 kcal mol–1 maximum). Accounting for transition states using standard transition state theory will also give barriers larger than those from ΔG° steps between intermediates, while concerted processes and/or proton tunneling may well lower the effective barriers.

Figure 5

Plots of ∑ΔG° (kcal mol–1) along the trajectory from 1 to 14. The top panel (A) is from the 2008 mechanism (from Table 4 for energies);[2] the bottom panel (B) represents new work presented here. Solid black circles represent energies obtained from PW91 calculations, while reddish diamonds represent energies obtained from B3LYP* calculations. The experimental energies from 1 → 2 and 2 → 3 are obtained from refs (26) and (44)–[46].

One protonation step that promotes the catalytic cycle is to add one additional proton to Y′O– in the transition 1 → 2 to make the state 2H, so that the added proton transfer from the KB path to Y237 is effectively concerted with 2e– reduction of Fe3+–Cu2+ to give Fe2+–Cu+. Our computed redox potentials for 1 → 2H are 0.24 V for PW91 and −0.02 V for B3LYP*, versus 0.275 V (experimental average[44,46] for CuII → CuI, +0.34 V, and FeIII → FeII, +0.21 V, for mitochondria) or 0.314 V using the experimental redox potential for FeIII → FeII for Tt ba3.[45] The corresponding error in ΔG° is about 1.6–3.1 kcal mol–1 for PW91 and about 13.5–15.5 kcal mol–1 for B3LYP*. In contrast to state 1 → 2H (O → R), the calculated O2 binding affinity (for standard state [O2] (aqueous) = 1 M) is much poorer for PW91 than B3LYP*, ΔG° = −17 and −10 kcal mol–1, respectively, versus −8.3 kcal mol–1 experimental for Tt ba3 from the reaction kinetics study of Szundi et al.[26] (see also Figure 2; state 2H → 3, R → A). These few available experimental values are plotted in Figure 5B. Experimental kinetics work of Hallen and Nilsson[47] shows that there is no proton uptake when O2 binds (noting that these are experiments in mitochondria), which is consistent with state 2H → 3 (Figure 5B) rather than 2 → 3 (Figure 5A) and is correctly represented in PW91 but not B3LYP*.

The KB path into the DNC based on Tiefenbrunn, Fee, and Stout’s high-resolution X-ray structure[19] is depicted in Figure 6. All residue side chains are constructed based on the X-ray density map, except for one added water (center purple), which completes the hydrogen-bonding network from Glu15(II) (from subunit II) to Y237 of the DNC and involves a3-ROH from the lower heme edge. The orientation of Figure 6 is 180° rotated about the vertical axis with respect to Figures 3 and 4 so that the Cu site is toward the left (above Y237) and Fe is toward the right. Figure 6 shows the structure of the KB path for protons into the DNC in Tt ba3. The Glu15(II) anion is positioned to accept a proton from the solution (n side) and return to a hydrogen-bond alignment. Then a proton can be readily transferred via a concerted Grotthuss-type process[48] (in common language, “a bucket brigade mechanism”) to Y237–. Finally, the Glu15(II) anion picks up another proton. Usually we expect that the destination for the proton is Y237–, but the geranylgeranyl a3-ROH group provides an additional entry point into the complex where an alternative pathway is needed.

Figure 6

Structure of the K path in cytochrome ba3 as deduced from a combination of structural and mutational analyses. The top two residues are actually part of the active site structure per se. Both Glu15(II) and Y237 are protonated. This figure is taken from Figure 2 of ref (35). Reprinted with permission from ref (35). Copyright 2009 National Academy of Sciences.

To examine the physical/chemical forces that drive proton pumping across the DNC, we now look at the path from state 4 to 5 and include an additional tautomer 5. Compared to state 4, state 5 has a proton shifted from H282 to H376 to form a state with H376H+, H282–, and YOH. The proton-shift distance is much smaller than that directly from state 4 to 5, and the predicted energy cost is small (7 kcal mol–1 for state 4 → 5, followed by −5 kcal mol–1 for 5 → 5, in place of 2 kcal mol–1 for state 4 → 5 in PW91; in B3LYP*, both energy differences are near 2 kcal mol–1). State 5 is then generated by shifting a proton in 5 from neutral Y237 (YOH) to H282– to form YO–, H282H. In Figure 7, the structures of these states are presented with the dipole moment vectors superimposed in red and the magnitudes labeled [μ in Debye (D) units]. Because dipole moment vectors are oriented from negative to positive, it is clearly seen that the proton shifts strongly reorient the very large dipole moment vector from 4 → 5 → 5, so that state 4 is almost oppositely aligned to state 5. The corresponding electrostatic potential becomes more positive in the vicinity of H376H+, or alternatively in the nearby water cluster (top of Figures 3 and 4). The following proton entering via the KB path ends in state 6 (see Figure 3), but the proton entry state is tautomer 6 (FeIII–O–O2––CuII, YOH). State 6 is 6 kcal mol–1 below state 6 in PW91, so the path 5 → 6 → 6 has energy differences (−8.7 and +6 kcal mol–1). After proton entry in state 6, the proton shift to form the hydroperoxide in state 6 produces an increase in the dipole moment magnitude of about 13 D, in a direction similar to that in state 5, again enhancing the electrostatic potential near H376. As shown in Figures 3 and 5B, the energy difference 6 → 7 is only mildly endergonic (11 kcal mol–1 in PW91; pumping against an assumed proton gradient from the pH = 7 n side to the pH = 3 p side, equivalent to 5.5 kcal mol–1).

Figure 7

Dipole moments of intermediates 4, 5, and 5. These three structures are overall neutral.

As expected on physical grounds, reducing CuII → CuI in the preceding step 1 → 2 increases the proton affinity of the covalently linked Y327–H233 pair bonded to Cu (ΔG° = −11.4 kcal mol–1 for PW91 for 2 → 2H, so ΔG° is lowered for 14 → 2H to 11.1 kcal mol–1 compared to that for 14 → 2, +22.5 kcal mol–1). This is also consistent with the experimental observation that O2 binding is not pH-dependent because that proton enters earlier at state 2H.

After O2 binding and formation of a peroxo-type bridge and higher Fe and Cu oxidation states, (Fe3+, Fe4+) and Cu2+, we expect that Cu2+ would acidify Y237–H233 (i.e., the Y′OH site) via through-bond interactions with Cu2+. This is supported by relative ΔG° calculations for the process Y′OH(H282–) → Y′O–(H282H), where ΔG° = −8.7 ± 2.7 and −6.5 ± 3.7 kcal mol–1 for PW91 and B3LYP*, respectively (average, range). This driving force strongly facilitates proton transfer from Y237 to H282. This average and range is based on our calculations for the tautomers of seven different states (5, 6, 9, 10, 11, 13, and 14; Figure 3), which are all Cu2+ species.

From Figures 5B and 3, it is evident that the entire pathway from states 3 → 9 has only low-energy barriers, at least judging from the energy differences between these intermediates, and the overall pathway is highly exergonic, which is consistent with the rapid rate observed from states A → P (observed as P). The e– transfer from heme-b (FeII) to the tyrosine radical P → P is clearly very fast, faster than 5 μs from experiment (Figure 2), and the earlier “true peroxo “ intermediates are not observed (Figure 3). The calculated ΔG° value for 8H → 9 is −23 kcal mol–1 (corresponding to a redox potential of about +1.22 V; we assume that the redox potential of the donor is near that of cytochrome c, about 0.22 V). If the proton coming in must go to Y237 first, then ΔG° for tautomer 9 is more positive, −10.4 kcal mol–1 (redox potential of about +0.67 versus SHE). However, this is still quite exergonic. For the next redox step, from state 11 to 12, ΔG° is again very negative, −21 kcal mol–1 (redox potential of +1.12 V versus SHE). There are two coupled protonations assumed here, and it is assumed that e– transfer from heme-b Fe is to CuII → CuI. If state 12 is skipped in favor of state 13, so that e– transfer goes from heme-b FeII to heme-a3 Fe FeIV → FeIII directly, the predicted ΔG0 value is almost the same. The corresponding spectral state changes of the DNC are P → F′ and P → O for reduction at Cu and Fe, respectively. As is seen in Figure 2, the transition from P → O is observed to go fairly slowly, and P → F′ is not directly observed. CuI is, however, not a good spectroscopic reporter using optical difference spectroscopy, and spectral analysis is focused on heme-a3 Fe. Typically, the state F′ is not seen in the normal reaction cycle of CcO but can be generated by binding H2O2 to the reduced state R (FeII and CuI) in the DNC. This shortcircuits the normal catalytic cycle, and the intermediate state F′ slowly changes back to state O (resting oxidized).[49,50] In any event, the comparative slowness of the observed P → O transition may well be a consequence of the coupled proton transfers and associated proton rearrangement. (All calculations above are with PW91, unless noted otherwise.)

Large Model: Structures and Energies

In Tables 1 and 2 (see also Figure 4, right), we present structures and energies for a number of states of the large model. We focus on two arcs with OLYP potential, one encompassing the “true peroxo” states (4 → 5 → 6), Table 1, and the second covering the last 1e– reduction ending in the resting state O (11 → 12 → 13 → 14 → 1; P → F′ → O), Table 2. The size of the quantum cluster is increased to 205–207 atoms by adding structural waters observed in the Tiefenbrunn X-ray structure and a glycine fragment that hydrogen bonds to internal waters near the reaction site. The larger size of the cluster allows us to model the reaction geometries with less constraints. The Fe–Cu distance is no longer constrained to X-ray geometries as in work by Fee et al.[2] Convergence of the electronic structure self-consistent field and geometries is improved. Also, the OLYP potential used typically has a better balance in describing the relative energies of different site spin states for Fe and typically yields more accurate spin coupling energies between transition-metal sites compared to PW91. In work by Fee et al. (2008),[2] we employed parallel spin on (ferromagnetic coupling, F) Fe and Cu sites for ease of computation and kept the observed experimental site spins over the reaction cycle. By contrast, the newer calculations take account of both F and AF (antiferromagnetic) alignment of the site spins, and different possible site spins as well (see Table 1), and our recent paper.[37]

Table 1

OLYP-Calculated Energies (E, offset by −27900 kcal mol–1), S2 Expectation Value, Net Charge of the Cluster (Q), and Key Geometric Data (Distances in Angstroms and Angles in Degrees) for Both F- and AF-Coupled Low-Spin Fe3+–Cu2+ Compounds 4, 5, 5, 6, 6H, and 6 with a Larger Modela

		Eb	O1–O2	Fe–O1	Fe–N(H384)	Cu–O2	Cu–N(H233)	Cu–N(H282)	Cu–N(H283)	Fe···Cu	O···O(Y237)	∠Fe–O1–O2	∠Fe–O1–O2–Cu	S2	Q
4	F	–33.5	1.33	1.87	2.19	2.08	2.14	2.12	2.18	4.60	3.28	125	–137	2.04	0
	AF	–32.8 (+0.3)	1.34	1.82	2.29	2.07	2.19	2.13	2.12	4.40	3.21	121	–123	0.54	0
5	F	–35.1	1.36	1.87	2.20	2.03	2.10	2.12	2.17	4.42	2.77	115	–135	2.04	0
	AF	–35.2 (+0.0)	1.34	1.84	2.30	2.10	2.13	2.14	2.14	4.44	2.83	122	–123	0.61	0
52	F	–24.7	1.35	1.87	2.19	2.02	2.16	2.02	2.21	4.49	3.35	123	–135	2.04	0
	AF	–27.3 (−1.1)	1.35	1.82	2.27	2.07	2.23	2.03	2.19	4.39	3.40	122	–123	0.59	0
6	F	–26.4	1.53	1.77	2.18	2.35	2.08	2.06	2.11	4.78	3.10	117	–128	2.02	1
	AF	–28.5 (−2.2)	1.49	1.79	2.15	2.66	2.03	2.06	2.07	5.07	3.08	116	–134	1.02	1
6H	F	–19.0	1.46	1.82	2.13	2.64	2.03	2.06	2.03	5.07	3.41	115	–137	2.04	2
	AF	–15.0 (+1.7)	1.48	1.78	2.11	2.68	2.03	2.06	2.03	5.04	3.47	118	–125	0.6	2
62	F	–34.4	1.35	1.88	2.21	2.05	2.14	2.13	2.16	4.46	3.44	113	–136	2.04	1
	AF	–36.3 (−0.9)	1.34	1.83	2.24	2.10	2.18	2.14	2.14	4.45	3.47	121	–127	0.63	1

Note that structure 4 of the larger model here does not have the water molecule above His376. O1 is the atom that is closer to the Fe site.

F stands for ferromagnetic coupling and AF for antiferromagnetic coupling. The energies given in the AF-coupled states are broken-symmetry state energies. The spin projection corrections (kcal mol–1) are given in parentheses.

Table 2

OLYP Calculations on Larger DNC Clusters in the Later Cycle with the Scalar Water on Cu+/2+a

	geometry									net spin
stateb	Fe–N(H384)	Fe–O1	Cu–O2	O1···O2	Fe···Cu	O···O(Y237)	E	Q	S2	Fe3+	O1	O2	Cu2+
112_2Ht1 (F)c	2.07	1.83	2.32	2.72	5.68	3.35	–49.26	2	3.79	1.29	0.15	0.02	0.37
112_2Ht2 (F)d	2.18	1.66	2.06	2.65	5.24	3.32	–44.91	2	3.78	1.31	0.77	0.09	0.54
12_H2O (F)e	2.20	1.65	3.83	2.95	5.74	3.41	–162.48	1	2.02	1.23	0.84	0.01	0.00
13_H2O (F)f	2.07	1.84	2.39	2.67	5.55	3.12	–158.30	1	2.03	0.95	0.14	0.02	0.33
14_H2O (F)g	2.07	1.85	2.20	2.60	5.39	3.15	–155.05	0	2.02	1.01	0.13	0.02	0.34
14_H2O (AF)g	2.08	1.85	2.28	2.64	5.46	3.14	–156.36	0	1.02	–0.94	–0.13	0.02	0.34
1_H2O (F)h	2.07	1.84	2.44	2.68	5.57	3.10	–152.10	0	2.02	0.94	0.14	0.02	0.37
1_H2O (AF)h	2.07	1.84	2.43	2.68	5.57	3.10	–152.40	0	1.02	0.94	0.14	–0.02	–0.37

OLYP calculations have been performed on the larger model. The calculated properties given here include energies (E, offset by −28000 kcal mol–1), S2 expectation values, net charges of the cluster (Q), the Mulliken net spin polarizations for the Fe, Cu, O1, and O2 atoms, and the key geometric data (distances in angstroms).

F stands for ferromagnetic coupling and AF for antiferromagnetic coupling. In the following notes, LS represents low spin and IS represents intermediate spin.

Formally adding two protons to state 11, one on the H282– side chain and another one on Fe4+=O2–. This is tautomer 1 (t1) of 11_2H with lower energy (see footnote d for tautomer 2). The optimized structure is a mixed state of Fe3+,LS–OH–···H2O–Cu2+(H376H+,H282H,Y237•), Fe4+,IS–OH–···H2O–Cu2+(H376H+,H282H,Y237–), and Fe4+,IS–OH–···H2O–Cu+(H376H+,H282H,Y237•). The calculated net spin on the Y237• side chain is 0.93.

State Fe4+,IS=O···H2O–Cu2+(H376H+,H282H,Y237). This is the second tautomer (t2) of state 11_2H with higher energy.

State Fe4+,IS=O···H2O···Cu+(H376H+,H282H,Y237). This is formal state 12 + H2O, where H2O is scalar water on Cu+.

State Fe3+,LS–OH–···H2O–Cu2+(Y237–). This is formal state 13 + H2O.

State Fe3+,LS–OH–···H2O–Cu2+(H376,H282H,Y237–). Note that in this larger model no water molecule is added on the top of H376.

State Fe3+,LS–OH–···H2O–Cu2+(H376H+,H282–,Y237–). Spin projection corrections are: −1.3, −0.3 kcal mol−1 for 14_H (AF) and 1_H (AF).

In Tables 1 and 2, simple redox potentials without coupled protonation events can be calculated with the equationwhere the absolute standard hydrogen electrode energy value is taken as ΔE(SHE) = −4.34 eV = −100.1 kcal mol–1 and Eox and Ered are the energies of the oxidized and reduced states. ΔE(SHE) corresponds to the best measured solvation free energy of a proton ΔGsolv = −264.0 kcal mol–1, where the ionization threshold energy of a free electron is taken as 0, with no account of electron entropy (see the Appendix and references therein). Similarly, ΔGdeprot can be computed from ADF energies; we take a typical ΔE(ZPE) (ZPE = zero-point energy) for deprotonation of hydroxyl oxygen as about −8.7 kcal mol–1, and this gives the protonation energy correction for single proton addition as −4.5 kcal mol–1 at pH = 7 (see the Appendix, Table A3). This is a DFT correction, −2.9 kcal mol−1 plus an empirical correction (−1.6 kcal mol−1), to the ADF energy differences E(HA) – E(A–), from Table 1. For different tautomers in Tables 1 and 2, energy differences are directly meaningful without correction, assuming only that ΔE(ZPE) cancels for deprotonation/protonation. In Table 1, we include the small spin projection energy corrections from the broken symmetry (AF) to the proper pure spin singlet states, in parentheses. From Table 1, state 4, and ref (37), 1e– reduction of state 4 is very unfavorable, redox potential E° = −0.73 V vs SHE, so cytochrome c (or heme-b) cannot further reduce the FeIII–(O–O)2––CuII complex. Protonation of state 6 at the tyrosine anion (Y237–) is very unfavorable (ΔGprot = 7.2 kcal mol–1, about 5.3 pH units more acid than pH = 7), which is important because the tyrosine anion is well adapted to become a radical after 1e– transfer to FeIII–(OOH–)–CuII, which would be suppressed if the neutral tyrosine (YOH) were very stable. The reaction pathway 4 → 5 → 5 → 6 → 6 has ΔG = +5.1, −6.8, −4.1, and +6.5 kcal mol–1 and is similar to the energy profile for the smaller model S2. The protonation free energy for 5 → 6 −4.1 kcal mol–1 includes an empirical correction (see the Appendix, Table A3). Structurally, it is clear that the hydroperoxide in state 6 is very important for breaking the O–O bond because the activated AF state has an O–O distance of 1.49 Å. The calculated F-coupled state 6 has an even longer O–O bond, 1.53 Å, and is only 4.3 kcal mol–1 above AF. Even more importantly, the Cu–O2 (bond to hydroperoxo O atom) bond is significantly shorter in the F state, 2.35 Å compared to 2.66 Å in the AF state, and is on the path to a CuII–OH structure after O–O bond breaking. The Fe–O1 bond is also close to that expected for ferryl oxo in both F and AF. Also, the charged a3 ROH···OY– hydrogen bond is very short, O···O = 2.83 Å in state 5, and still short in state 6 but much longer for all states with neutral Tyr–OH (Y237). This is entirely consistent with what we saw in ref (37), where the X-ray structure of Tiefenbrunn[19] was proposed to have a charged a3 ROH bond. Turning now to the reaction pathway for the final e– transfer P → O, it is clear that obtaining FeIII–OH– would be facilitated by having a nearby proton, but in the earlier small model (S2; Figure 3), we proposed that H2O exited the reaction pocket directly between Fe and Cu at this stage. Instead, we now retain this scalar H2O in states 12_HO, 13_HO, 14_HO, and 1_HO (F′ and O). In place of state 11 in Figure 3, we now add 2H+ to form two tautomeric states (the first two lines in Table 2). We find that the addition of these two protons is very favorable and so should precede 1e– reduction. These two states are related by a proton shift from Y237 (YOH) to replace the proton shifted from H2O–CuII to FeIV=O. Our calculations indicate that to some extent an electron may accompany this proton transfer onto the Fe, but this state is a resonance mixture of valence isomers. The proton-shifted tautomer 11_2H (tautomer 1; Table 2, line 1) is lower in energy than 11_2H (tautomer 2) by about 5 kcal mol–1, while upon 1e– reduction, the CuI state is predicted to be lower, by about 4 kcal mol–1 (0.17 eV). The corresponding redox potential is about 0.57 V (reduction of Cu) or 0.40 V (reduction of Fe). The actual reduction process may not follow the lowest energies because it is favorable to preserve the proton shift to form Fe–OH– once this is formed. Further, kinetic measurements of Farver, Fee, and co-workers[27] implicate total reorganization energies in the range of 1.0–1.5 eV for e– transfer between heme-b and heme-a3 Fe. Both the driving force and the reorganization energies then are important for the timing of e– transfer.

Table A3

ΔGprot,shift Calculations (kcal mol–1)

	ΔGprot,shift
potential	DFT	emp	final
OLYPa	–2.9	–1.6	–4.5
OLYPb	–5.2	–1.6	–6.8
PW91c	–2.0	–3.0	–5.0
OPBEc	–5.0	0.0	–5.0
B3LYPc	+7.2	–3.0	+4.2
B3LYP*c	+5.6	–4.4	+1.2

ΔGprot,shift(DFT) depends on −ΔZPE. Here −ΔZPE = +8.7 kcal mol–1 for use in the catalytic reaction cycle for OLYP.

–ΔZPE = +6.5 kcal mol–1 from the average OH bond ZPE for H2O in OLYP. For use in Table A2.

–ΔZPE = +6.6 kcal mol–1 for all other potentials for both the catalytic reaction cycle and thermochemistry in Table A2.

The FeIII–OH–···H2O–CuII state (13_HO) will provide a focal point for the last part of the reaction cycle (second arc) extending from 11_2H → 13_HO → 14_HO → 1_HO. In Figure 8, we present the energy profile from the last 1e– reduction to the last proton pumping as assumed within the model cycle to the final reset to the initial state 1. We also compare the energetic pathway with one where the scalar H2O on CuII has dissociated (Figure 9). Typically, the energy profile is similar to that without H2O. However, the presence of H2O on CuII facilitates proton transfer between CuII–H2O and FeIV=O2– or FeIII–OH–. Further, there is recent X-ray structural evidence for two waters or OH––H2O between FeIII and CuII in the resting state O.[19,37] We do find that ΔG(H2O dissociation) = −9 kcal mol–1 (OLYP), favoring dissociation, and less, −0.6 kcal mol–1 (PW91) and +3.2 kcal mol–1 (B3LYP*), but based on the active site structure, kinetic trapping for the CuII-bound H2O is likely. Two waters are probably released after 2e– reduction of state 1_HO → 2H (O → R) with a lowering of the free energy. This conclusion is supported by X-ray structural evidence.[18,19,37]

Figure 8

Plots of ∑ΔG° (kcal mol–1) relative to state 13_HO along the trajectory from 11_2H (lower energy tautomer 1) to 1_HO in the cycle with the scalar water on CuII. The PW91 and B3LYP* energies were obtained for the smaller model (see Table S2 in the Supporting Information, SI). OLYP calculations were performed on the larger models (Tables 2 and A3).

Figure 9

Plots of ∑ΔG° (kcal mol–1) relative to state 13 along the trajectory from 11_2H to 1 in the cycle without the scalar water on CuII. The PW91 and B3LYP* energies were obtained for the smaller model (see Tables S2 and S3 in the SI). OLYP calculations were performed on the larger models (Tables 3 and A3).

Table 3

OLYP Calculations on Larger DNC Clusters in the Later Cycle without the Scalar Water on Cu+/2+a

				net spin
stateb	E	Q	S2	Fe2+/3+/4+	Cu2+/+
112 (F)	–28052.20	0	3.78	1.28	0.42
12	–27838.93	1	2.02	1.20	0.00
13 (F)	–27836.04	1	2.03	1.12	0.18
13 (AF)	–27835.24	1	1.04	0.94	–0.22
14 (F)c	–27836.30	0	2.03	1.14	0.19
14 (AF)c	–27836.23	0	1.03	0.94	–0.21
1 (F)	–27830.82	0	2.03	1.09	0.25
1 (AF)	–27829.53	0	1.03	0.93	–0.25
2H	–27708.87	1	6.22	3.84	0.00

OLYP calculations have been performed on the larger model. The calculated properties given here include energies (E in kcal mol–1), S2 expectation values, net charges of the cluster (Q), and the Mulliken net spin polarizations for Fe and Cu sites.

F stands for ferromagnetic coupling and AF for antiferromagnetic coupling.

Note that in this larger model no water molecule is added on the top of H376.

Returning to the energy profiles in Figure 8, the calculated ∑ΔG° values for 1_HO relative to state 13_HO are +14.0 and +15.6 kcal mol–1 for the large OLYP model (without and with the empirical deprotonation correction) and +19.0 and +18.3 kcal mol–1 for PW91 and B3LYP*. In Figure 9, the corresponding energies are +13.6, +15.2, +21.2, and +14.8 kcal mol–1 for OLYP, OLYPcorr, PW91, and B3LYP*. These are not high energy differences for proton pumping considering that a +5.5 kcal mol–1 pumping cost (for pH = 7 → pH = 3) is included.

Bifurcation Points, Chemical and Vector Protons, and the Exit Pathway

In light of the experimental structural, kinetics, and mutagenesis data for ba3 from Tt, the idea that there is a largely common pathway for scalar and vector protons and that this derives from the KB pathway should not be surprising. We can see two reasonable bifurcation points. The first, mentioned earlier (see Figure 6), is at a3-ROH and the adjacent Tyr237. This bifurcation point is mainly utilized as needed to bypass Tyr237 when it is in a radical state, and a proton needs to be passed on, especially a chemical proton entering the reaction chamber between Fe and Cu at that stage. The other bifurcation point, which we think is much more common in the reaction path, occurs at His282. From this location, a proton on His282 (neutral) can move onto a CuII-bound hydroxo or peroxo (scalar proton), with the different tautomers being close in energy from our DFT calculations. Alternatively, a vector proton can move onto the exit position (PLS) and proceed on the exit pathway. Specifically, starting with His282–, a proton enters typically from Tyr237 to generate His282 neutral. This will weaken the hydrogen bond to threonine (T302), allowing its side chain to rotate and deposit a proton on His376 or more probably in the nearby water pool (see Figure 4, large model, right). (This location is near propionate A, as previously proposed for the PLS location.) Above His376 in Figures 3 and 4 is an Asp287–Arg225 salt bridge. Catalytic turnover, proton pumping, and mutagenesis studies of Chang et al.[28] indicate that Asp287 is probably important for efficient proton pumping; these results are still incomplete, and only one mutation has been examined for Arg225, which is so far inconclusive. We propose that the Asp287–Arg225 salt bridge is potentially a simple electrostatic switch or gate, which could be opened by a strongly positive electrostatic potential nearby, for example, by having a proton in the adjacent water pool. The energy barrier to rupturing the Asp–Arg salt bridge could be paid when the proton enters the water pool or during transit. Also, the local electric field of this proton (and from the surrounding charge distribution) possibly acts in opposition to the long-range electric field from the membrane potential. This model is currently under study.

Conclusions

On the basis of our earlier work on the catalytic cycle of the cytochrome ba3 DNC, we have found an improved reaction path by modifying the sequence of proton uptake and proton-transfer events and examining alternative tautomeric states. The new 14-step reaction cycle results in a much “smoother” energy pathway for the catalytic reaction cycle, comparing parts A and B of Figure 5 with lower barriers between intermediates. The entire pathway from state 3 → 9 has low-energy barriers and proceeds downhill. We have also tracked the DFT-calculated dipole moments of peroxo and hydroperoxo intermediates after O2 binding (Figure 7). Here we find that there are very substantial changes in the magnitude and direction of the dipole moment of the DNC with tautomeric shifts, and a progressive increase in the positive electrostatic potential near His376 (in addition to that due to protonation of His376 itself). We then constructed large model structures and evaluated geometries and energies (OLYP potential) for two arcs of the cycle, one for “true peroxo” states 4 → 5 → 6 and the other for the last 1e– reduction from state 11 → 1 (spectroscopic states P → O). Important results from these DFT calculations are presented in Figures 8 and 9 and Tables 2 and 3. We have tested the effect of keeping the scalar water bound to CuII in state 11 through the remainder of the catalytic cycle both for the large model (OLYP potential) and for the smaller model (PW91 and B3LYP*) potentials. The energy profiles in both Figures 8 and 9 are similar, but the water on CuII facilitates protonation of Fe-oxo. The energy barriers for the final proton-pumping step 13 → 1 are about 15–16 kcal mol–1 for OLYP, large model, and somewhat higher for B3LYP* (15–18 kcal mol–1 for B3LYP*, small model) and for PW91 (19–21 kcal mol–1). These results include an assumed 5.5 kcal mol–1 cost for proton pumping at one step. Further, we discuss a potential proton exit pathway involving breaking of an Asp–Arg salt bridge. The exiting proton would be carried or transferred via a water pool from inside to outside. Earlier on the vector proton pathway, we have proposed that His282 and Thr302 play an important role.

There are many important issues still to be understood in the reaction cycle and proton-pumping mechanism of CcO both in ba3 class enzymes and more broadly including also aa3 enzymes (bacterial and mitochondrial). We have not yet addressed the kinetics in detail, including transition states, effects of sequence differences, and mutational effects. The proton exit path from the DNC is probably complex.[28] More needs to be done to integrate the KB path and the other redox centers into the models presented here. The scalar water exit pathway is also essential for proper analysis. In general, however, we expect that the structural and energetic maps of the reaction cycle in the DNC that we have generated will contribute to a better representation of the coupling of electron and proton transfer to proton pumping in CcO.

Table A1

Empirically (emp) Corrected E(H+) (eV)

potential	E(H+)	ΔEcorr(H+)	final E(H+)	ΔGprot,shift(emp) (kcal mol–1)
OLYP	12.64	+0.07	12.71	–1.6
PW91	12.51	+0.13	12.64	–3.0
OPBE	12.64	0.00	12.64	0.0
B3LYP	12.11	+0.13	12.24	–3.0
B3LYP*	12.17	+0.19	12.36	–4.4

Table A2

Reaction Free Energies for Five Different Exchange Correlation Potentials

reaction	PW91	OLYP	OPBE	B3LYP	B3LYP*	Exp
O2(gas) + 2H2(gas) → 2H2O(gas)	–99.9	–91.1	–93.6	–97.2	–96.1	–109.3
O2(aq) + 4e– (from cyt c) + 8H+(aq,pH=7) → 2H2O(liquid) + 4H+(aq,pH=3)a	–27.3 (−15.3)	–19.2 (−12.8)	–21.4 (−21.4)	–26.2 (−14.2)	–24.7 (−7.1)	–36.9
4ΔGprot,shift(emp)b	–12.0	–6.4	0.0	–12.0	–17.6

Corrected energies are given without parentheses, and uncorrected energies are in parentheses.

4ΔGprot,shift(emp) is the total protonation free-energy empirical correction for 4H+ added to O2, forming 2H2O.