Reaction pathways, proton transfer, and proton pumping in ba3 class cytochrome c oxidase: perspectives from DFT quantum chemistry and molecular dynamics

After drawing comparisons between the reaction pathways of cytochrome c oxidase (CcO, Complex 4) and the preceding complex cytochrome bc1 (Complex 3), both being proton pumping complexes along the electron transport chain, we provide an analysis of the reaction pathways in bacterial ba3 class CcO, comparing spectroscopic results and kinetics observations with results from DFT calculations. For an important arc of the catalytic cycle in CcO, we can trace the energy pathways for the chemical protons and show how these pathways drive proton pumping of the vectorial protons. We then explore the proton loading network above the Fe heme a3–CuB catalytic center, showing how protons are loaded in and then released by combining DFT-based reaction energies with molecular dynamics simulations over states of that cycle. We also propose some additional reaction pathways for the chemical and vector protons based on our recent work with spectroscopic support.


Introduction
The flow of electrons through an ordered series of proteins embedded in the inner membranes of mitochondria and in aerobic bacteria provides the driving force for proton transfer across these membranes (Nicholls and Ferguson, 2002;Rouault and Tong, 2008;Noodleman et al., 2020).(Figure 1A) In an intact inner membrane containing all the relevant integral membrane protein complexes (Complexes 1-4), with sufficient reducing equivalents from NADH (a two-electron-one-proton donor plus an added proton (2e − , 2H + )), and with access to adequate molecular oxygen O 2 (the terminal electron-proton acceptor for the chemistry) in complex 4, a substantial transmembrane electrochemical potential difference is generated.The negative potential side is inside and the positive potential side outside of the membrane, the proton flow being from inside to outside (Noodleman et al., 2020) (Figure 1B).This potential difference in turn drives proton transfer back through the neighboring protein complex 5 (across the membrane in the opposite direction), turning a rotor within complex 5, and produces ATP from ADP plus inorganic phosphate (Pi) by a dehydration reaction.
2 Comparing cytochrome bc1 and cytochrome c oxidase (complex 3 with complex 4) To have a frame of reference for the mechanism and significance of reaction steps in complex 4 (cytochrome c oxidase), we compare some broad features of complex 4 with those in complex 3 (cytochrome bc 1 ).In the electron transport chain (ETC), protons are pumped across the membrane at complexes 1, 3, and 4 (Figure 1A).Already at complex 1, two electrons and two protons (2e − , 2H + ) are added to the mobile carrier molecule, oxidized ubiquinone (Q), to produce the 2e − reduced form (QH 2 ), a two-electron carrier.Complex 2, which does not pump protons, still adds reducing equivalents to Q, producing additional QH 2 .In its reduced form, QH 2 readily flows through the membrane and is mobile as well within protein complex 3, being able to transit from the upper to the lower part of the protein in both its oxidized form Q and reduced form QH 2 .Both forms Q and QH 2 are uncharged.Furthermore, the long hydrophobic sidechain of the quinone makes it amphiphilic overall, so it can move easily between hydrophilic and hydrophobic regions.
Complex 3 acts as a transducer of electron flow, converting the 2e − transfer by QH 2 into 1e − transfers within complex 3, on to cytochrome c, the mobile protein carrier for 1e − delivery to complex 4. A relevant question for the mechanism of complex 3, and more generally across biological systems, is how this and related electron transduction systems work.The key feature is a bifurcation (splitting) of electron transfer between a high potential pathway (here through the Rieske protein subunit and on to cytochrome c 1 within subunit b) and a low potential pathway (through cytochrome b L and ending with cytochrome b H , from lower to higher redox potential).The initial reaction works in two steps: 1): QH 2 → Q − (radical) + 2H + + 1e − ; then, an oxidation step occurs, 2): Q − (radical) → Q + 1e − , with the acceptors having corresponding redox potentials, high potential for step 1 and low potential for step 2. The high potential pathway of step 1 has as its first acceptor the nearby [2Fe-2S] cluster of the Rieske protein subunit.The one electron (1e − ) transfer via this pathway is required to produce the high-energy Q − (radical) state.
Unlike typical [2Fe-2S] clusters with four cysteines 4(SCys − ) linked to 2Fe, the Rieske FeS center has 2(SCys − ) bound to one Fe and 2(NHis) bound to the alternate Fe (Figure 2) (Ullmann et al., 2002).So how does the relatively high redox potential of the Rieske [2Fe-2S] cluster and protein arise?It was discovered that the Rieske protein has a pH-dependent redox potential in complex 3 systems attributed to additional proton binding on 1e − reduction to the εnitrogen of the histidines, those N atoms not directly bound to Fe.The pH dependence of the redox potential was examined in detail FIGURE 1 (A) Schematic of the four integral membrane protein complexes I-IV in the electron transport chain in mitochondria and aerobic bacteria.The inner (lower) region below the membrane as depicted is the n-side (−) corresponding to the matrix space in the mitochondria and to the cytoplasm in bacteria.The outer (upper) region above the membrane depicted is the p-side (+), intermembrane space in mitochondria, and periplasm in bacteria.The electrons are delivered to cytochrome c oxidase (CcO) by cytochrome c (cyt c), as shown on the p-side of the membrane.Reprint with permission from Figure 1 of Rouault and Tong (2008), copyright 2008 Elsevier.(B) The electrochemical potential gradient across the membrane.
experimentally in the Rieske protein fragment isolated from Thermus thermophilus (Tt) bacteria (Zu et al., 2001), and this property is similar in the mitochondrial (bovine) Rieske protein fragment (Kuila et al., 1992;Link et al., 1992;Baymann et al., 1999).Independently at about the same time, our group demonstrated that one electron reduction of the Rieske cluster involves a coupled single-proton transfer using a combination of density functional theory (DFT) calculations on the Rieske [2Fe-2S] cluster and electrostatic methods from the Poisson-Boltzmann theory for the surrounding protein (Ullmann et al., 2002).Analysis of the pH dependence of the redox potential of the Rieske protein fragment (Tt) gives experimental effective pK a values for two sites in the oxidized form pK a = 7.8 and 9.6 and for the reduced form pK a = 12.5 (approximately) for both sites from the fit.If we consider protonation of an empty site at a typical pH = 7 on one of the two His during 1e − reduction, the calculated ΔG = 1.37 × (7-12.5)kcal/mol = −7.54kcal/mol, equivalent to a positive shift in the redox potential equal to +330 mV.This clearly is a large contributor to the high redox potential of the Rieske protein connecting to the high potential pathway.Both pH-dependent NMR and chemical modification species support the deprotonation of histidine 154 (Tt) (Hsueh et al., 2010;Konkle et al., 2010), the first site deprotonated with increasing pH in the oxidized Rieske protein.
However, we encounter another more subtle issue.In the oxidized form of Tt, the lowest affinity site still has a pK a = 7.8.At pH = 7, the fractional protonation (occupation) of that site is 0.86 from the Henderson-Hasselbalch equation, so protonation of that site from the reduced quinone QH 2 will be strongly inhibited, unless the site is previously emptied by another mechanism.Rieske proteins isolated from other cytochrome bc 1 species show similar oxidized pK a values, including from bovine mitochondria (Link et al., 1992) and Rhodobacter sphaeroides (Rs) cyanobacteria (Zu et al., 2003).We noticed that all the pH-dependent redox measurements were carried out on the isolated Rieske fragment and after a search found that in the intact bc 1 enzyme, there is a critical lysine residue with a positively charged mobile sidechain in the cytochrome b subunit, very near the Rieske iron-sulfur cluster.Assuming this lysine sidechain coordinates close to the histidines, as is feasible in the oxidized redox state, this interaction will favor deprotonation of the lower-affinity histidine site.This lysine is expected to move away from the histidines in the reduced state, where the pK a values are much higher.In a collaboration with two experimental groups, this proposal was successfully tested (Francia et al., 2019).This lysine sidechain was found to be isostructurally conserved across many species, including both prokaryotes and eukaryotes.Furthermore, experimental tests on Rhodobacter capsulatus (Rc) cyanobacteria, acting in photosynthetic mode, showed that both the native lysine residue and an arginine mutant in that position are fully active for both electron transfer and proton pumping, while anionic amino acid residues in that position are completely inactive, and neutral amino acids have diminished activity.Rc cyanobacteria are functionally and structurally similar to Rs bacteria.
Once the Rieske cluster is 1e − reduced and protonated by 1H + , the following electron and proton transfer involves a large-scale movement of the [2Fe-2S] cluster and protein (Rieske head group) by several angstroms to dock to a site near the embedded cytochrome c 1 cofactor.The proton proceeds through the exit path to the intermembrane region, while the single electron is transferred from cytochrome c 1 to cytochrome c, which is the 1e − donor to cytochrome c oxidase, see Figure 1A and the discussion and Figure 2 in Noodleman et al. (2020).Thus, the proton transfer that occurs in complex 3 in the Rieske protein results from a deprotonated histidine accepting an H + along with the electron as the QH 2 is oxidized to semiquinone.The second proton is transferred to a different histidine (not within the Rieske cluster).
The structure and mechanism of complex 4 are quite different.There are no mobile electron and proton carrier molecules like ubiquinone or the Rieske protein head group.Instead, as we will see, water plays a major role, with smaller changes in protein geometry.There are no bifurcations of electron flow in CcO.Instead, there are bifurcations of proton flow, which are connected with oxygen binding and electron transfers.There are also differences in the mechanisms and function of the bacterial ba 3 class of CcO, exemplified by Thermus thermophilus, and aa 3 class CcO found in other aerobic bacteria and in mitochondria.

Structure, mechanism, and function of cytochrome c oxidases
In Figure 3A, we present a structural global representation of the ba 3 -type cytochrome c oxidase from Thermus thermophilus (Tt) (Chang et al., 2009).Electrons enter by transfer from the oneelectron (1e − ) protein carrier molecule cytochrome c and are passed to Cu A , which is a dinuclear Cu cluster.From Cu A , electrons transfer via tunneling through the protein to the low-spin mononuclear heme Fe b (in this Figure hidden   histidine ligand, a nearby Cu site called Cu B bonded to three histidines, and the surrounding protein and water molecules (see Figure 3B).One of these histidines below Cu(B) has a special covalent linkage to a tyrosine (N-C) via their sidechains (His233 to Tyr237 in Tt ba 3 ).This linkage is universal in both ba 3 -and aa 3 -type CcOs and has not been found so far in proteins other than CcO.Protons enter the CcO protein from the cytoplasm connected to the membrane inner surface, as schematically shown in Figure 3. CcO proteins are embedded in the inner membrane both in the mitochondria and in bacteria.The transfer of protons from the inside (n-side, negative) to the outside (p-side, positive) of the membrane and the opposite transfer of electrons from outside to inside produce an electrochemical potential difference across the membrane (Figure 1B).Importantly, this potential difference arises not only from proton pumping but also from the catalytic reaction chemistry cycle, which removes positive charge (protons) from the inside and negative charge (electrons) from the outside and combines these with O 2 to produce water.The overall reaction is as follows: Here, 4e − and 4H + are used in the chemical reaction cycle, and nH + is pumped from the inside to outside.O 2 enters via a hydrophobic channel from the membrane's interior, passing near Cu B and binding to Fe a3 starting from its reduced Fe 2+ state.(see Figure 3A) In the equation above, the total amount of charge transferred is Q = 4 + n electron charge equivalents per catalytic cycle, where n is between 2 and 4. To have a sufficient number of electrons to produce a full reaction cycle from the oxidized enzyme, 4e − must be added to CcO in 1e − steps originating sequentially from four Cyt c molecules.At that fully reduced (4e − ) stage, Cu A has redox state Cu + -Cu + for the copper cluster starting from the mixedvalence oxidized state Cu  Noodleman et al. (2014).
Figure 4 shows the comparison of structures and proton entry paths for A-and B-family CcOs, with R. sphaeroides (A-family, aa 3type) and T. thermophilus (B-family, ba 3 -type) as examples (Chang et al., 2009;Noodleman et al., 2020).There are two different entry paths (D and K paths) to the DNC in aa 3 enzymes and only one (K path analog) in ba 3 enzymes (Wikström et al., 2018).For the D-path in R. sphaeroides, the proton entry involves an Asp132 (D132), which takes a proton from the cytoplasmic region (Figure 1B) (analogous to the matrix space of the mitochondrial inner membrane in animals).This proton then passes through a series of sites, including several water molecules (green spheres as depicted) and protein sidechains, until it reaches the Glu286 (E286).This residue sidechain is at or near a branch point, where protons may go either into the reactive region where O 2 binds and is transformed or protons may be diverted for proton pumping, depending on the state of the reaction cycle and other conditions.For the other proton input path, in the aa 3 family, a Glu (E101 of subunit II in R. sphaeroides) picks up the proton from the cytoplasm (n-side), and then this proton is transported from site to site, eventually reaching the special tyrosine (Y288).A lysine residue (K362) plays an important role in the transport through the K path in aa 3 enzymes.In comparison, in the B-family (ba 3 as in T. thermophilus), there is no D-path for protons into the reactive complex region or for proton pumping, but only a single-proton input path starting at the Glu (E15 of subunit II in T. thermophilus) and then leading into the reactive region via the special tyrosine (Y237).
Table 1 summarizes the number of vectorial and chemical protons that are transported via the different proton channels in aa 3 -and ba 3 -type enzymes.From proton pumping experiments on cytochrome ba 3 in lipid vesicles, Rauhamaki and Wikstrom found that two or more vector protons (three or four) are pumped per catalytic cycle, compared to four in cytochrome aa 3 , but with considerable uncertainties (Rauhamäki and Wikström, 2014).
They attribute this variability either to the greater sensitivity of ba 3 enzymes to higher membrane potentials or to increased proton backflow compared to aa 3 enzymes plus measurement difficulties.Table 2 shows the different modes of transport feasible for electrons and for protons as well as their respective ranges for single tunneling events.

Catalytic reaction cycle in the DNC
We will focus on the major arc of the DNC catalytic reaction cycle from the reduced state R binding to O 2 (state A) through several intermediates to the reactive oxidized state, labeled O H ,    (Noodleman et al., 2020).Over a full cycle from State R and returning to State R, four electrons (4e − ) are taken up, as shown.Four chemical protons are also taken up from input channels.In Figure 5, the sequence of states and transitions shows a close similarity between the ba 3 class and aa 3 class enzymes from optical difference spectroscopy, but with some differences in O 2 binding affinity and kinetics through the cycle.First, we have a few comments about the sequence of intermediates, kinetics, and other physical parameters portrayed in Figure 5.The kinetic rate constants for various transitions depicted in Figure 5 are taken from the work of Einarsdottir et al. (Szundi et al., 2010) using optical difference spectroscopy on ba 3 CcO intermediates from Tt in the detergent to track the Proton Grotthuss process--fast proton transfer along chains of water molecules, also in combination with amino acid sidechains (bucket brigade mechanism) Modes of proton and electron mobility in proteins and membranes.electron transfer from one metal cofactor (reduced versus oxidized states) to the next.(Their group and others have done comparable work on aa 3 class enzymes, where there is some loss in accuracy for kinetics parameters due to the overlapping optical difference absorption bands for heme a and heme a 3 ).For our purpose, the most important aspects of the reaction cycle kinetics are the very fast rates from R → F and the much slower rates from F → O H (Bloch et al., 2004).State O H is an activated state which promotes the following electron transfer to generate State E H and, therefore, activates the next catalytic cycle (Bloch et al., 2004).If the electron transfer occurs too slowly, O H decays to the "resting" State O, which is more difficult to reduce.While these differences in kinetics are now evident, the structural differences between States O H and O are not, either in aa 3 -type or in ba 3 -type CcOs.We used detailed DFT modeling (te Velde et al., 2001;van Lenthe and Baerends, 2003;Baerends et al., 2017) et al., 2020;Han Du et al., 2022).Important work involving spectroscopy and kinetics continues on this problem.For CcO in detergents [in Einarsdottir's and related work (Siletsky et al., 2007;Szundi et al., 2010)], there is no net proton pumping, no electrostatic potential is developed because there is no membrane, and there is no net sidedness (orientation) in physical space for the CcO protein molecules.In the presence of a membrane, as, for example, in experiments where CcO proteins are embedded in lipid vesicles with preferred directionality, the transfer rates from state to state are slower because a membrane potential is developed producing added energy barriers for proton pumping and barriers involving the chemical reaction cycle steps as well.We also follow the chemical reaction cycle with resonance Raman vibrational spectroscopy from the work of Rousseau et al. (Egawa et al., 2012;Ishigami et al., 2015).The observed Fe-O stretch frequencies for the series of States are shown in Figure 5 for an aa 3 class enzyme; ba 3 intermediates and Fe-O frequencies are expected to be similar.We have calculated DFT Fe-O frequencies for our ba 3 quantum cluster models (Han Du et al., 2019) showing good correspondence with the measured frequencies in aa 3 (Figure 6B), and the sequence of intermediates from DFT tracks well with the observed sequence of different frequencies from resonance Raman spectroscopy.Figure 7 shows an excellent linear correlation between the DFT-calculated Fe-O distances over the reaction pathway with a function of the DFT-calculated Fe-O frequencies (Han Du et al., 2019) [an empirical inverse power law correlation called Badger's rule (Badger, 1935)].A simplified map of the different states through the cycle is also provided (Figure 6A) (Han Du et al., 2019).In the bar graph comparing the DFT-calculated versus experimental resonance Raman frequencies, we want to emphasize some points.There are no experimentally observed resonance Raman frequencies for peroxo or hydroperoxo states, probably due to the very fast kinetic steps going through these intermediates starting from State A to states P M and P R .In related synthetic systems with peroxo or hydroperoxo bridges, the measured resonance Raman frequencies do fit those calculated here by DFT for CcO well.The experimental absence of these two states does raise the question of which is, in fact, present and in what form?
State A, where O 2 binds to ferrous Fe (see the prior State R) to form the "ferrous-adduct," has an important electronic structure and illustrates some of the difficulties in obtaining accurate electronic structure descriptions.Figure 3B shows the geometric structure of State A. State A also sets the stage for the subsequent rapid sequence of reactions.During the formation of State A from the reduced State R, where the high-spin ferrous (heme a 3 ) metal ion Fe 2+ reacts with the high-spin triplet ground state of O 2 , effectively, there is a 1e − transfer from Fe 2+ to O 2 , with partial covalent-ionic bonding between low-spin ferric Fe 3+ and low-spin superoxide radical anion O − 2 .Substantial negative charge transfer occurs from the Fe 2+ to the O 2 , but less than one full electron.The calculated Mayer bond order (Mayer, 1983) is 0.48 for the Fe-O stretch in state A, which slowly increases to 0.54 and 0.59 for the following peroxo and hydroperoxo states, while the O-O Mayer bond order successively decreases from 1.28 for State A to 1.06 for peroxo and 0.911 for hydroperoxo.As expected, the Fe-O bond strengthens and shortens through this sequence of states, while the O-O bond weakens and lengthens in preparation for breaking.The calculated Fe-O and O-O stretch frequencies follow this same pattern, see Figure 7.The DFT-predicted Fe-O stretch frequency for state A of 518 cm −1 deviates by approximately 10% from the experimental value of 571 cm −1 , while the DFT-calculated O-O stretch for superoxide is accurate to approximately 5% compared with experiments in comparable systems, like oxymyoglobin.While the calculated Fe-O stretch frequency error is the largest percentage for the entire observed reaction cycle, this error is still fairly small.The O 2 binding process is intrinsically complex, involving 1e − oxidation of the Fe 2+ site, reduction of O 2 , and weakening of the π bond in forming the superoxide radical anion O − 2 , conversion of both fragments into low-spin S = 1/2 forms, and covalent bonding with partial back charge transfer of electron density from O − 2 to Fe 3+ in the porphyrin site.We have shown (Noodleman et al., 2020) that the overall reaction here involves some compensating errors, particularly involving the oxidation of the Fe 2+ complex and the partial bond breaking of the O 2 π bond.Proceeding to the next steps of the catalytic cycle, after the assumed 1e − transfer to O − 2 from Cu + → Cu 2+ in Cu B to form the Fe 3+ -(O-O) 2− -Cu 2+ , peroxo bridge state, the following proton transfer from the special tyrosine (Tyr-OH) shows a significant positive energy difference (uphill) to form the hydroperoxo bridge state, calculated by DFT, approximately 12 kcal/ mol higher than that of the peroxo state.By contrast, when we followed the actual O-OH bond cleavage of the hydroperoxo state along the O-O reaction coordinate, this energy barrier is quite small, approximately 1-2 kcal/mol (Han Du et al., 2016).The DFTcalculated (approximate) transition state has an O-O bond length near 1.7 Å, which then generates the final state P M , with the final fourth electron to break the O-O bond derived by 1e − transfer from the special tyrosine (Tyr237 − anion in T. thermophilus) through the covalently linked His233 to Cu 2+ and on to O 2 .
Let us now examine in detail some of the structural features of the DNC from T. thermophilus using our DFT calculated structural model for the Fe-superoxide-bound State A (see Figure 3B).Near the top of the DNC is the sidechain ring of His376, a small water pool, the sidechain of Asp372, and the heme a 3 proprionate sidechain (Prop-A).These residues form the start of the proton-loading network (PLN) for the vector proton transit.His376 is central to this network and forms the proton-loading site for vector protons, as our group and others have proposed.In the hydroperoxo state, the Asp372 is thought to be protonated, as discussed later, and His376 is doubly protonated, forming a histidine cation with a high pK a , approximately pK a = 11.0 from our DFT calculations (OLYP potential with D3 dispersion) (Han Du et al., 2016).After O-O bond cleavage in the hydroperoxo state to form State P M , the His376 + pK a remains similar, and DFT-calculated pK a = 11.8.The next steps are formation of State P R by "external" 1e − transfer from heme b Fe 2+ and then proton transfer to Tyr237.The proton transfer source is a glutamate/glutamic acid at the starting point of the K pathway into the DNC along a hydrogenbonded chain of residues and water molecules to the Tyr237.At this stage, after net (1e − ,1H + ) transferred to the Tyr237 neutral radical, we have calculated the His376 + pK a and again found it to be high, pK a = 9.2 (Han Du et al., 2016).These three states, hydroperoxo, P M , and P R after Tyr237 protonation (neutral) form convenient reference states to examine proton transfers in the protonloading network and exit paths.
First, however, we explore the state-to-state analysis from P R → F → O H in more detail, which will illuminate the driving force for proton pumping at this stage of the cycle, see Figure 8 (Han Du et al., 2016).In Figure 8, the number and types of intermediates amplify on  those summarized in Figure 5.The number of intermediates from P R → F → O H has been expanded with more detailed typing.These states are labeled from a → b → c → d → e → f → g → h.We also distinguish three types of hydration models near the DNC central reaction chamber, M1, M2, and M3.These models contain explicit quantum-described water molecules.The first hydration model (M1) is built based on the observed crystallographic water molecules in the ba 3 T. thermophilus X-ray structure of Tiefenbrunn et al. (Hunsicker-Wang et al., 2005;Tiefenbrunn et al., 2011) in the resting oxidized state (state O) and then geometry optimized using DFT in correspondence with the various forms of State P R (see models a and b in Figure 8).Models c, d, e, and f are of state F-type, but with 2H 2 O and protons shifting positions from state to state.We put in the observed water molecules that directly H-bond to the various oxygen species, here the Cu 2+ B -OH − or Cu 2+ B -H 2 O in the P R → F → O H sequence, the other nearby hydrogen bonding (hydrating water molecules) within the DNC active site, as well as the small water pool that is found at the top of the DNC, interacting with His376, Asp372, and Prop-A.Hydration models M2 and M3 have the positions of some water molecules significantly shifted, as we will describe later.All of these water molecules lie within the quantum cluster model.Later on, when we address the classical molecular dynamics (MD) studies, many additional water molecules will be included.
The DFT computational and structural model uses the OLYP-D3(BJ) exchange correlation potential with D3(BJ) for dispersion effects and an overall quantum cluster size of approximately 205 atoms for the DNC with small variations depending on the overall protonation state.Since, in contrast with the MD studies, the surrounding protein, including other redox centers, hydrating water molecules, and the membrane, are not included in the quantum chemical model, we employed a continuum dielectric with dielectric constant ε = 18.5 to represent the overall polar environment outside the quantum cluster model.When protons are added to the quantum cluster (labeled H + uptake), the approximate ΔG calculations assume proton addition at pH = 7 from the donating dielectric medium.In the actual biological system, for the ba 3 enzyme, the protons added to the DNC come from the K path and originate from a glutamic acid/glutamate, as mentioned previously.For proton transfers within the DNC, the ΔG values are calculated from the ΔpK a values between the different groups, including the calculated zero point energy differences between residues or between residues and bound H 2 O or OH − .
Returning to the reaction pathway in Figure 8, state a, P R (M1), corresponds to the state labeled P R within the larger global reaction cycle in Figure 5. Proton uptake gives a DFT-calculated approximate ΔG = −8.8kcal/mol for protonation of Tyr237 − to form neutral Tyr in state b.State b is also the same one where the pK a (His376 + ) = 9.2, as mentioned above.The proton transfer from Tyr237 (neutral, protonated) to generate Cu 2+ -H 2 O from Cu 2+ -OH − has a calculated ΔG = −2.7 kcal/mol.The following proton uptake at Tyr237 has a ΔG = −5.5 kcal/mol, see the angled arrow directed from state c → d.The water model is still in State M1 at this stage.For all the states described with hydration M1, we now focus on the water dimer, labeled (2H 2 O), which is closest to the OH − or H 2 O, which is the Cu 2+ -bound ligand.In this dimer, the first H 2 O is directly H-bonded to the OH − or H 2 O Cu 2+ ligand, while the second H 2 O is H-bonded to the other water of the dimer [these two water molecules correspond to the crystallographic water molecules HOH604 and HOH608 depicted in Figure 3B].We treat the whole water dimer (2H 2 O) as moving together to the upper water pool around His376 after the Tyr237 − is protonated.The calculated free energy cost is small, ΔG wds = +4.2kcal/mol for the process model d (M1) → model d (M2) (wds = water dimer shift).There is a major associated change in the histidine cation (His376 + ) pK a from pK a (i) = 12.2 for (M1) before to pK a (f) = 6.5 for (M2) after the water dimer shift (i = initial, f = final), so ΔG i→f −1.37(pK a (f ) − pK a (i)) +7.8 kcal/mol for the relative ΔG change in proton binding.Deprotonation of the His376 cation has a far more favorable free energy change in the final state (water model M2) compared to the initial state (water model M1).The corresponding ΔG deprot (His376 + ) = 1.37 (pK a (f) − pH) = −0.7 kcal/mol at pH = 7.Therefore, the total ΔG for state d (M1) → state e (M2) with a proton disassociated from His376 + is +3.5 kcal/mol, which is positive and small, followed by a small but negative proton transfer energy (−1.4 kcal/mol) from Tyr237 to restore the full protonation of the His376 + cation in state f.We attribute the large pK a shift in His376 + to the descreening of the Cu 2+ -H 2 O charge distribution due to the water dimer relocation from the lower pool near Cu 2+ -H 2 O to the upper pool near His376 + .We also separately tested the effect of adding two more water molecules to the upper pool while keeping the other water molecules in the M1 arrangement [(M1+2H 2 O) Model].The calculated pK a = 10.7 is just 1.5 pH units below the pK a found in State d (M1).State e (M2) (pK a = 6.5) is much different from either of these M1 states, supporting our descreening proposal.In State g, M3 has a similar hydration model to state M2, with the water dimer in the upper pool near His376 + , but the H 2 O ligand to Cu 2+ has moved to a hydrogen bonding orientation to the ferryl-oxo (FeO)(2+).In the next redox and proton uptake step, the active oxidized complex State O H is formed with a bridging hydroxyl (OH − ) derived from the oxo ligand to Fe. Considering the full catalytic reaction cycle, each full cycle generates 2H 2 O from [O 2 + 4e − + 4H + ], as seen in Equation ( 1), and in Figure 5.These 2H 2 O products per cycle must flow out at some stage.We propose that these two water molecules per cycle become the H-bonding water dimer to the Cu 2+ -OH − or H 2 O ligands, and subsequently the water dimer moves to the upper water pool before being expelled via the water exit pathways.A similar mechanism applies if the two water product molecules exchange with other water molecules in the reaction pathway before exit.

Proton loading network and transfers by DFT
At this point, let us return to the problems of the proton uptake and proton transfer within the proton loading network (PLN) by quantum chemistry (DFT) and by molecular dynamics (MD) methods (Yang et al., 2016a).A schematic diagram for the various feasible protonation states involving the PLN is shown in Figure 9.All states of this proton uptake and proton transfer cycle have either two or three active protons, as found from our DFT calculations with OLYP-D3 (exchange-correlation potential with D3 dispersion correction).The PW91-D3 form of DFT gives generally similar relative energies.The DNC quantum cluster model we used in Yang et al. (2016a) corresponds to the hydroperoxo state in a model containing 204 atoms (205 atoms when there are three protons in the PLN).[We note that this hydroperoxo state has a free energy of protonation at pH = 7 for His376 → His376 + (−5.5 kcal/ mol) similar to P R state d (M1) (−7.1 kcal/mol) in Figure 8 and also has a hydration model like M1.]In the protonation and proton transfer cycle, state D, with His376 + , is closely analogous to state d (M1); see below.)As shown in Figure 9, the main protonation sites in the PLN are on the Asp372 carboxylate sidechain, the proprionate A sidechain (PRA) of the heme a 3 , and the His376 ring nitrogen (both ε and δ).All of these groups are within the quantum chemical cluster model, and so the relative energies of these states can be calculated using DFT.We consider the sequence of states starting at state A (two protons in the PLN, on Asp372 and His376), and then a shift of one proton from Asp372 → PRA occurs with a histidine ring flip to form state B → proton uptake to form State C′, and another proton transfer from PRA → His376 to form His376 + (doubly protonated), State D. Figure 10 shows that the proton transfer from Asp372 → PRA is substantially endergonic (product State B) based on our DFT calculations both for OLYP-D3 and PW91-D3, and that state C′ lies well above state D with the same number of active protons (Noodleman et al., 2020) in this network.Asp372 is the most geometrically accessible protonation site from residues and water molecules below.As depicted in the figure, from state B → state B2 corresponds to a proton transfer from PRA → His376, which becomes (His376 + , doubly protonated), and State B2 has a similar or lower energy than state B. State B2 also allows the proton uptake at Asp372 to result in the lower energy State D, while only the higher energy state C′ is directly attained by proton uptake from state B. Based on our earlier discussion, P R state d (M1) → M2 will lead to much weaker proton binding in His376 + and readily allow the removal of a proton by Glu126 (subunit II) sidechain anion, as depicted in Figure 9, States E and then F (equivalent to State A′, His376 neutral, protonated at N δ).Furthermore, State C′ has a much higher energy than State D by 15 kcal/mol (OLYP-D3) or 20 kcal/mol (PW91-D3), while the pK a of His376 + is near 11 or 12 in State D. This result implies that even after a strong lowering of this pK a by approximately 6 pH units, backflow of a proton via State C (or C') to the lower H-bonded network (vectorially toward the n-side) will be strongly diminished compared to the forward motion through state D → E → F and out through the residue sidechainlined water pathways.

Proton loading and transfer by combining DFT and MD methods
By shifting the emphasis from the quantum DFT calculations to our classical molecular dynamics (MD) simulations [using the AMBER software package (Götz et al., 2012;Salomon-Ferrer et al., 2013a;Salomon-Ferrer et al., 2013b;Case et al., 2014) and AMBER ff12SB and Lipid14 force fields (Dickson et al., 2014;Maier et al., 2015)], we can see their implications.The details of the methods and the models for these MD calculations are presented in Yang et al. (2016a); Yang et al. (2016b), but a few main points be emphasized.In states A, B, B2, C, and D, all the proton transfers lie within the overall quantum cluster.In States E, F, and A, however, proton transfers occur outside the quantum cluster region.Therefore, to analyze all of these states, MD models include the whole protein along with the metal cofactors (Cu A cluster, heme b Fe, and the dinuclear heme a 3 Fe-Cu B ), as well as a representative membrane, explicit solvent and many mobile water molecules in the protein interior, and ions for overall neutrality.The partial charges of the metal cofactor regions are calculated from DFT (OLYP-D3) to obtain electrostatic potential (ESP)-based charge distributions which are combined with standard charges from an AMBER force field model and other AMBER force field parameters.The protonation states of the various protein sidechains and for proprionate A are given in Table 1 of Yang et al. (2016a), and the fractional hydrogen bond occupancies for the different states from the MD simulations in Table 2 and  Changing the perspective from proton transfer energies, we now consider changes in geometries with related energies, water mobility, and occupancies as described by MD methods.Over a selected finite simulation time (140 ns), we have calculated the number of H 2 O molecules that exchange from the upper water pool to bulk water or vice versa via two pathways (small water channels), P1 and P2.The upper water pool is bounded by His376, Asp372, PRA carboxylate, and other adjacent residues, sidechains of Glu126 II , Gln151 II , Arg225, Asp287, Arg449, Arg450, and Asn377.See Figures 11,  12.We do not count the water molecules within these two channels.The total number of H 2 O molecules exchanged (total for P1 + P2) for States A, B, C, D, E, and F (A') are 15, 2, 5, 0, 5, and 16.Clearly, state D shows the least exchange and states A and F the most exchange.The proton released from state D (His376 + ) to Glu126 II and then onward opens these water exchange channels.We can expect that when the His376 + pK a value is greatly lowered, the proton transfer from His376 + to the Glu126 II anion will eliminate that salt bridge, which then leads to H 2 O exchanges.Another intact salt bridge, Arg225-Asp287, is observed nearby.We have calculated the potential of mean force (PMF using MD) by displacement of the reaction coordinate between C(ζ) of the guanidinium from Arg225 and carbon C(γ) of Asp287 carboxylate, Figure 13.This salt bridge does not break, but the distortion to larger distances costs far more free energy in State D than in States E or A. The resistance to stretching the Arg225-Asp287 salt bridge is also consistent with the very low water mobility we find between the upper water pool and bulk water in State D. Overall, the change in the electrostatic potential within the DNC, particularly near the Cu 2+ B site with a water shift, promotes the water flow and proton exit from the upper water pool through the two exit channels.

Splitting the proton entry path into inner and outer branches
Another large-scale problem arises, i.e., where does the proton entry path split into two, with one path (the inner path) entering the Cu B -Fe a3 "reaction chamber" (language from JA Fee), and the other path (the outer path) providing vector protons for binding to the proton loading network and exit?How are these two paths sequenced or synchronized?It is clear that Tyr237 is involved at the end of the proton uptake K path, but what occurs after that point? Figure 14 depicts that structural region in a different way from that in earlier figures, emphasizing the surroundings of Cu B .For the inner path, let us concentrate on the pathway for forming the Cu 2+ B bound hydroperoxide starting from the reduced state R, which, upon binding O 2 , becomes a ferric-superoxide complex (Fe 3+ -O − 2 ) (state R → state A → hydroperoxide).We propose this negative charge transfer from Fe to O 2 triggers strong geometric changes and a proton transfer involving Cu B , which is in the cuprous Cu + state.Cu + -histidine(N) bonds can be labile, and then the broken bond can protonate at the histidine nitrogen.Starting from Figure 14, we propose that His283, during Cu + -His(N) cleavage, picks up a proton by a ring flip from His282 protonating to His283 + , and this proton loss on His282 − is then restored by Thr302, whose proton is resupplied by Tyr237 and other amino acid carriers and water molecules.Both the local Cu ligands and Thr302 (or similarly a serine) are well-conserved isostructurally across ba 3 and aa 3 CcO′s.Trp229 is also structurally well-conserved and well-placed to form a stabilizing π-cation interaction with His283 + .From this point, a 1e − transfer from Cu + to the O − 2 will be linked to a proton transfer from His283 to form the hydroperoxide state, which will also allow for reformation of the Cu 2+ -His(N) covalent bond.This proton transfer is likely concerted with the 1e − transfer, although these may be asynchronous.This pathway avoids two problems: (Rouault and Tong, 2008) direct electron transfer from Cu + to O − 2 and bond formation without a rapid proton transfer could yield the very-low-energy peroxo state Fe 3+ -(O-O) 2− -Cu 2+ (DFT OLYP-D3BJ energy calculated as 12 kcal/ mol below the following hydroperoxide state, protonated near Cu 2+ ), which would strongly inhibit proton transfer from Tyr237 according to our DFT OLYP-D3-BJ calculations (Han Du et al., 2019).In this same work, we noted that the peroxo bridge structure is considerably distorted geometrically compared to the superoxide state (state A) with a peroxo calculated Fe-Cu distance of 3.88 Å and Cu-O distance of 1.94 Å, compared to Fe-Cu 4.55 Å, and Cu-O distances between 2.95 and 3.2 Å in state A. From a kinetics viewpoint, coordinating large charge transfers and large geometry changes likely requires a high reaction barrier.By contrast, the hydroperoxo state has a bridge structure, calculated Fe-Cu distance 4.31 Å, and Cu-O distance 2.11 Å, closer to these distances in state A (Nicholls and Ferguson, 2002).The alternative proton pathway from Tyr237 to a nearby water molecule and then into the peroxide-bound state is long (approximately 4 Å) and likely results in a proton bound to the peroxide, generating a hydroperoxide at the higher-energy Fe 3+ side instead of placing the OH near Cu 2+ (OLYP-D3 energy difference 11.8 kcal/mol) (Han Du et al., 2016).By contrast, the mechanism we propose is consistent with the mobility and linkage of proton transfers seen in other systems (Grotthuss process) (de Grotthuss, 1806) and the lability of Cu + -histidine bonds (See Table 2).In earlier work, we proposed a similar mechanism from DFT calculations and found strong supporting evidence from the X-ray structure in the reaction pathway of Cu-Zn superoxide dismutase (CuZnSOD) (Konecny et al., 1999;Shin et al., 2009;Perry et al., 2010).
Here, the first 1e − transfer from the initial superoxide bound to Cu 2+ forms Cu + , and O 2 leaves.Then, a following 1H + transfer breaks a Cu-histidine bond (to the histidine bridging Cu to Zn covalently).When a second superoxide enters and reacts with Cu + , the 1e − transfer to O − 2 is accompanied by deprotonation of that histidine, formation of hydroperoxide (OOH) − , and reformation of the Cu 2+ -His bridge bond to Zn 2+ .X-ray structures of CuZnSOD in the reduced Cu + form have a broken Cu-His bond, while in other structures of oxidized proteins, this bond is intact.
In answer to the initial question raised, it is likely that Thr302 (or similar isostructural Thr or Ser in other CcO′s) is a significant branch point: a negative local electrostatic potential in the (Fe 3+ -O − 2 )-Cu + state A drives protonation at His283 near Cu + in the inner path, while in the outer path, a positive electrostatic potential near Cu 2+ -H 2 O and Tyr-OH drives deprotonation of His376 + and proton recovery (from State B2) via proton transit from Thr302 to Asp372 using intervening water molecules and amino acid sidechains.As discussed before, State C' lies well above State D by 15-20 kcal/mol.Using a simple Marcus-type model for the reorganization energy and driving force associated with the proton transfer from state A' (reactant) to state C' (product) versus state A' to state D (product) to lowest order (linear), the ΔΔG (activation) = ΔG (C' → D)/2, so the activation barrier is calculated to be lower by 7.5-10 kcal/mol in State D compared to State C' (see chapters 5 and 7 in Connors (1990).
After O-O bond cleavage, the formation of the radical transient state P M is expected with a Tyr237 radical (ba 3 ) until 1e − is transferred from the mononuclear low-spin heme b (in ba 3 ) or  Frontiers in Chemistry frontiersin.org14 Noodleman et al. 10.3389/fchem.2023.1186022analogously mononuclear low-spin heme a (in aa 3 enzymes).During normal cycling, the electron transfer from heme b or heme a is rapid.
In contrast, Yu et al. (2012) trapped two different Tyr radicals with their spectra determined by X-band and D-band EPR and further analyzed with QM-MM methods to compare calculated spin distributions with ring hyperfine parameters.These radicals were trapped at a low concentration in bovine CcO (aa 3 -type).A mixed valence state that stops at P M was obtained by starting with the oxidized bovine CcO and reacting this with hydrogen peroxide (HOOH).(Hydrogen peroxide does not have access to the active site in ba 3 CcO as in Thermus thermophilus in the oxidized state, so the analogous experiment has not been done.)These two Tyr radicals in bovine aa 3 CcO are Tyr244, covalently linked to His240, and Tyr129, with a nearby Trp236 which are analogous to Tyr237 linked to His233, and Tyr136 near to Trp229 in ba 3 CcO (Thermus thermophilus).Figure 14 shows the Cu B near surroundings based on the Tt ba 3 X-ray crystal structure.The covalent linkage of Tyr237 to His233 is shown.The His233 main-chain amide is hydrogen bonded to the Trp229 main-chain carbonyl.It is also evident from the residue numbering that there is a covalent linkage as well between His233 and Trp229 by a longer pathway.Yu et al.
have proposed (see their Figure 7) that at P M , a radical transfer mechanism can move a hole from the special Tyr244 radical to the higher positioned Tyr129, forming a cation radical with a loss of a proton later from the same tyrosine (Yu et al., 2012).Hole motion from Tyr244 to Tyr129 is equivalent to 1e − transfer from Tyr129 to Tyr244, forming a Tyr244 − anion.The proton on Tyr129 + cation radical will exit via the D-path available in bovine aa 3 CcO, and after Tyr129 is 1e − reduced by the mononuclear heme a, another proton will bind to Tyr129 from another position of H-bonded groups lower along the D-path.Something similar appears likely to us in ba 3 -type CcO.We emphasize more the role of the intervening Trp229 in Tt ba 3 (see Figure 14) (Noodleman et al., 2020).After the formation of the Trp229 + radical cation state, the nearby Tyr136 sidechain will be deprotonated via electrostatic repulsion.The Tyr136 − anion formed may persist or may exchange an electron with the Trp229 + cation radical.In either case, the Trp229 + radical or the Tyr136 neutral radical then receives an electron from the Fe 2+ on heme b.The Tyr136 − anion formed has its proton restored from His376 and adjacent residue sidechains, which are nearby.Toward the end of the oxidative cycle, we see the transition from the state labeled (f,M2) → (g,M3) → (h) (State O H ) (Figure 8).The first two states named are different forms of state F, where the H 2 O ligand on Cu 2+ shifts to a hydrogen bonding position to the ferryl-oxo oxygen (Fe 4+ = O), which is then followed by a high redox potential coupled electron-proton transfer (e − ,H + ), calculated by DFT to have a ΔG = −9.5 kcal/mol, calculated with respect to a typical cytochrome c redox potential (approximately 0.22 V versus SHE).This is one of the slower transitions involving protonation within the reaction chamber, here of an oxo to hydroxo state, along with the electron transfer.We propose a feasible mechanism for this reaction.In this mechanism, the H-bonded water transfers one proton to the oxo after or concerted with the one-electron transfer to the Fe 4+ → Fe 3+ .The newly formed hydroxyl (the one not bound to Fe 3+ ) will then first activate protonation of the Tyr237 − anion and then proton transfer to that hydroxyl to re-form an H-bonded water.This looks to us like the most probable mechanism.In state (g,M3), the DFT-calculated pK a is quite low, pK a = 1.8 to 4.0, but a nearby nearly free OH − should increase the effective pK a of Tyr237 via electrostatic stabilization for subsequent proton transfer to the same OH − .Over the full reaction cycle, the mechanisms for state A → hydroperoxo and state F → state O H that we have proposed account for the entry of two (out of four) chemical protons into the reaction chamber to bond to various oxygen species (see Table 1).The coupled (1e − ,1H + ) transition from State O H → E H reduces Cu 2+ to Cu + , so a mechanism similar to the one proposed for state A → hydroperoxo is likely feasible.The mechanism for the entry of the chemical proton to form Cu 2+ -H 2 O from Cu 2+ -OH − (P R → F) in Figure 5, or in more detail, state (b,M1) → state (c,M1) in Figure 8 is not clear, although the DFT-calculated ΔG = −3 kcal/mol approximately.We have also shown results from DFT and MD calculations supporting a mechanism for pumping of one vector proton in the state P R → state F transition.For the preceding State P M → state P R transition, we think a mechanism similar to the one proposed by Rousseau et al. (Yu et al., 2012) for aa 3 -type CcO based on EPR experiments applies also to ba 3 CcO in Thermus thermophilus.

Aspects of DFT parameters employed in our work
In terms of DFT exchange-correlation potentials, generalized gradient approximations (GGA) (including OLYP and PW91) are much more computationally efficient than hybrid DFT-Hartree-Fock Frontiers in Chemistry frontiersin.orgmethods (like B3LYP and B3LYP*), particularly when optimizing geometries.We have calculated geometries for some families of Fe complexes using OLYP and PW91, particularly iron-sulfur, iron-sulfur-nitrosyl, and iron-carboxylate-nitrosyl systems, and found that the predicted geometries agree better with the experiment than B3LYP or B3LYP*.
For relative spin state energetics, Vancoillie et al. (2010) have found that OLYP yields similar good-quality results with respect to experiment when compared to B3LYP and B3LYP*.Currently, all practical methods for predicting relative spin state energetics for Fe complexes and other transition metal complexes are not completely reliable, but some perform well on many systems.Radoń  D3 force field-type corrections for dispersion effects became available (from Grimme et al.), and we utilized these to further improve geometries in our work combining DFT and molecular dynamics studies (OLYP-D3, PW91-D3) (Yang et al., 2016a) and in our prior Fe-Cu hydroperoxide cleavage work (Han Du et al., 2016).Further improvements to dispersion force field terms yielded the D3(BJ) term, which was incorporated to obtain results with OLYP-D3(BJ) and PW91-D3(BJ) in our papers in 2018 (Han Du et al., 2018) and later (Han Du et al., 2019).Some B3LYP-D3(BJ) singlepoint calculations using PW91-D3(BJ) optimized geometries were also done.
DFT basis sets are Slater-type using the ADF code (te Velde et al., 2001;van Lenthe and Baerends, 2003;Baerends et al., 2017), triple zeta with polarization functions (TZP) for Fe and Cu, and double-zeta with polarization functions (DZP) for all other atoms, with inner core orbitals frozen.Basis sets do vary based on the requirements of the properties being calculated, for example, Mossbauer properties (no frozen core on Fe, Cu) or vibrational frequencies, where we have sometimes expanded the basis set.

Additional connections to the literature on cytochrome ba 3 and aa 3
For a very recent review from a broader viewpoint, we recommend the paper of Wikstrom, Gennis, and Rich, which combines an analysis of structures, kinetics, and spectroscopy with some theory/computation as well (Wikström et al., 2023).For DFT calculations and reaction energies for the aa 3 chemical cycle, the recent paper of Blomberg ( 2021) is valuable as it shows a different approach from ours.
A recent study of the reaction kinetics of ba 3 -type CcO in Tt at a low temperature (10 C) clearly shows the presence of an additional intermediate state just after the hydroperoxo state (see Figure 5), but before the O-O cleavage product state P R (Poiana et al., 2017).In the fully reduced enzyme, after O 2 binding, optical difference spectroscopy shows that beginning at State A, one electron leaves the heme b Fe 2+ well before the O-O cleaved P R state is observed.The relevant time constant ratio is 1/10 for the oxidation of heme b compared to the formation of P R at 10 C, but near 1/1 for this ratio at 45 C, which is closer to the typical physiological temperature range for Thermus thermophilus species (50 C-80 C in hot springs) (Dyer, 2003).These results implicate a newly observed low-lying intermediate state at the dinuclear reaction site, presumably with Cu B again reduced to Cu + .We have proposed that in forming the hydroperoxo state, there is initially a broken bond between His283 + and Cu + B before (1e − ,1H + ) transfer to the Fe 3+ -superoxide state.If after these transfers, the re-formation of the His283 to Cu 2+ B bond is slow, the reduction to Cu + B may occur with a small reorganization energy.In addition, this last 1e − transfer from heme b Fe 2+ reducing Cu B to Cu + can potentially occur either as a single-electron tunneling event or alternatively with a hole-type (amino acid radical) intermediate along the pathway.
There are very significant differences between the proposed aa 3 CcO pathways and what we think are feasible and probable in ba 3 class CcOs.In particular, aa 3 CcOs have two proton input paths K and D. Most proton transfers occur via the D-path, see Table 1.Along the D-path, there are water molecules and H-bonded residues where protons, one at a time, can be diverted by a glutamic acid in some proposed mechanisms toward a proposed proton loading site, suggested to be His291 and/or the nearby proprionate D, which can, in principle, be part of the chemical proton transfer path (His291) or part of the proton pumping path (His291 or Prop-D) (Popović and Stuchebrukhov, 2004;Popović, 2013).Geometrically, the His291 is above the Cu B and near a small water pool.The combined DFT/ electrostatic analysis (Popović et al., 2005;Popović and Stuchebrukhov, 2005;Makhov et al., 2006;Blomberg and Siegbahn, 2012) supports the presence of His291 as an anion (His291 − ), which can then be one site for proton loading, neutralizing the His291 charge, to be further transferred to oxygen species, like OH − , or pumped out.A number of other proton loading and exit pathways have been proposed in aa 3 CcOs, including sites near Prop-D, near Prop-A, or near the lowspin heme a, triggered by electron transfer from heme a Fe 2+ to heme a 3 Fe 3+ or Fe 4+ .In this context, early papers by Blomberg et al. (2000); Siegbahn et al. (2003) used DFT calculations on aa 3 from bovine heart mitochondria to calculate the O-O bond cleavage pathway, redox potentials, and pK a values for metal-bound groups and other groups, including tyrosine.A later paper on aa 3 class CcO in bacterial Rhodobacter sphaeroides used a larger quantum model (approximately 250 atoms) to focus on proton transfer near Prop-A and Prop-D of heme a (Blomberg and Siegbahn, 2010).
In ba 3 class CcOs, there is only one proton input path, the K path, which is located in a similar position to the K path in aa 3 CcOs (see Figure 4) but has some differences in the amino acid and water molecule placements.In any event, the absence of a D-path in ba 3 class CcO (Tt) and the resulting use of the K path for both chemical and vector (pumping) proton pathways leads to different mechanisms.The Cu B ligands in mitochondrial CcO are His240, crosslinked to Tyr244, His290, H-bonded to Thr309, and His291.The corresponding Cu B ligands in Tt CcO are His233 (crosslinked to Tyr237), His282, H-bonded to Thr302, and His283 (in a similar structural position to His291).We note that this Thr or Ser is broadly conserved in CcOs, both aa 3 and ba 3 .With protons coming from the K path in ba 3 CcO, our proposed mechanism for protonation of superoxide, O − 2 radical (state A) involves a His283 + cation which dissociates from Cu B in the Cu + state, a transient His282 − anion, and reprotonation from Thr302 via Tyr237 and the K path.We think that Thr302 is also a branch point for proton transfer via water molecules to Asp372.Therefore, these mechanisms are quite different, including which histidines are charged and when, but so are the proton input paths in aa 3 versus ba 3 enzymes.
For state P R and state O H in ba 3 , we have calculated pK a values and the redox potentials for Tyr237.These results are also relevant for vibrational frequency and pKa measurements for other tyrosines Tyr136 and Tyr133.For the P R state, the pK a (H 2 O) for H 2 O/OH − bound to Cu 2+ B with TyrO − equals 15.4.The pK a (Tyr) for TyrOH/ TyrO − with OH − bound to Cu 2+ B equals 13.4.The ΔG for proton transfer from TyrOH to Cu 2+ B OH − is thus −2.7 kcal/mol.The calculated redox potential for TyrO − /TyrO • with respect to SHE is +0.05V.The calculated redox potential for TyrOH/TyrO • with respect to SHE is +0.43 V, which is considerably higher, as expected, and also higher than that of cytochrome c. pH = 7 is assumed for the protonation state of the surrounding medium acting as the proton source.The cytochrome c redox potential is set empirically to +0.22 V.The Tyr is Tyr237 covalently linked to His233.Broadly, our results and proposed models agree with spectroscopic studies by Koutsoupakis et al. (2011) since the Tyr237 often has a high pK a , while they observe moderate pK a near 7.8 for other Tyr, presumably Tyr136 and/or Tyr133, through different reaction cycle states, showing shifts in the intensities of vibrational frequencies reflecting protonated versus deprotonated tyrosines.Tyr136 is along one exit path in our MD studies (Yang et al., 2016a).Tyr136 also plays a role in our proposed radical transfer mechanism in ba 3 , which follows the original experimental EPR work of Rousseau et al. on aa 3 .These data come from our recent papers (Han Du et al., 2018;Han Du et al., 2019).
For state O H , by comparing the OH − symmetric bridge to the H 2 O bridge in different geometries, we calculated the bridge H 2 O pK a between 11.7 and 12.2 (22).We did not calculate the corresponding TyrOH to TyrO − pK a .Following the trend of prior states in the reaction cycle, we assumed that Tyr237 is deprotonated here (TyrO − ) and then would pick up protons during the reductive part of the reaction cycle.
Cai et al. use electrostatic and molecular dynamics methods to characterize the protonation states and their statistical occupancies including relative energies during proton exit when stepping through the various redox, protonation, and oxygen species binding in the chemical cycle (Cai et al., 2018;Cai et al., 2020).The chemical cycle states interact via electrostatics with the protonation states along the exit path.There is no quantum chemical treatment of the catalytic reaction cycle, so that part of the energy is omitted.For the ba 3 enzyme in CcO, the authors focus particularly on the sites Asp372, Prop-A, His376, and Glu126(II).They find, on average, 1-2 "active protons" on these sites.This agrees with our count of 2-3 protons in this region before proton release with water outflow.There is a difference in language here since they count His + as having one "active" proton, while we count two protons (total) on His + .More significantly, in our work, we find that Prop-A when protonated is in a high-energy state, which helps block proton backflow.Cai et al. find some states where Prop-A is protonated at low energy, which disagrees with our results.
For the combined water-proton exit paths, there is a body of work dealing with the coordination geometry and energies of protonated water clusters, particularly the H 5 O + 2 Zundel-type ion bonded to anionic amino acid sidechains, as found in the electrostatics modeling of bacteriorhodopsin (Spassov et al., 2001).This is one example of a class of related possibilities for proton-water exit mechanisms (see Table 2, last three lines).Next, we discuss aspects of the sequence and structure dependence of proton-water exit pathways in ba 3 versus aa 3 .
In a publication by Sugitani and Stuchebrukhov (2009), the proposed exit pathway of aa 3 and ba 3 class enzymes is compared, using the different ba 3 enzyme structure to propose a different exit pathway than in the aa 3 enzyme.Their proposed ba 3 pathway agrees very well with our proposed pathway, called P1, see Figure 12.See also Supplementary Figure S5 in the Supplemental Material of our work, where we find two pathways P1 and P2, as shown in the figures in Yang et al. (2016a).
Gennis et al. studied reaction cycling and proton pumping via related exit pathways for various site-directed mutants in the exit pathways (Chang et al., 2012).There is no clear explanation for the specific effects of these different site-directed mutations on pumping and on catalytic cycling.We can make a few general comments based on our work, but we have analyzed only the native state of ba 3 CcO from Thermus thermophilus, and not from the mutants.Their experimental observations are that some point mutations of His376 and Asp372 keep proton pumping intact (H376Y, H376F, D372N, and D372V), while a few others disrupt proton pumping (H376N and D372I).In addition, other mutations produce low electron turnover in the catalytic cycle (H376A, H376D, and D372E), so residues in the exit path can disrupt the chemical catalytic cycle as well.Based on these experiments, our view is that there are multiple solutions to these exit paths.The proton pumping mechanism we have proposed for the transition of state P R → state F is connected to a large shift in the pK a of His376 by longrange electrostatic effects.Similar large shifts may be feasible for other residues substituting for His376 or nearby residues.Furthermore, the "typical" pK a values of residue sidechains can vary substantially depending on the surroundings, and even singlepoint mutations can alter local geometries and nearby water clusters.In addition, we find that Prop-A facilitates proton exit by blocking proton backflow, an effect which could well-persist after mutations near Prop-A.In addition, our proposed radical transfer mechanism in state P M leads to deprotonation of Tyr136 at a different location along the water and proton exit path.However, these comments do not answer the detailed issues raised by the proton pumping experiments in Chang et al. (2012).
In a very recent work, Yang et al. (2023) have analyzed the effects of some of these site-directed mutants of ba 3 CcO on the proton exit pathways using electrostatics and MD methods closely related to those we used in our collaborative work on the native enzyme (Yang et al., 2016a).The authors find that often site-directed mutations change structures of adjacent residue sidechains and that the water pathways are shifted often, consistent with our speculations mentioned above.They have been able to rationalize some experimentally observed effects of these mutants on proton exit, but many open questions remain.

Conclusion
We briefly summarize some of the main conclusions from this work.under an existing membrane potential?f.What are the effects of the regulatory subunits in mitochondria?Where are the three core subunits similar or different in aa 3 CcO enzymes from other organisms (various bacterial aa 3 ) or in ba 3 CcO enzymes?g.What are more general themes or insights that we can apply to other biological or chemical systems?
behind the DNC Fe a3 heme) and then on to the dinuclear center (DNC) [also called the binuclear center (BNC)].The DNC consists of heme Fe a3 with an axial
FIGURE 3 (A) Schematic of cytochrome ba 3 with noted functionalities.SU is a subunit.Reproduced from Figure 1 of Noodleman et al. (2014), copyright 2014 American Chemical Society (further permissions for reusing this figure should be directed to the ACS).(B) Dinuclear center (DNC) of cytochrome c oxidase consisting of heme Fe a3 , Cu B , and relevant residues and water molecules in the vicinity.Shown is a structural model for the adduct state A of the DNC in which molecular oxygen O 2 coordinates to Fe a3 .Some apolar hydrogen atoms are omitted for clarity.Reproduced from Figure 3 of ref.Han Du et al. (2019), copyright 2019 American Chemical Society (further permissions for reusing this figure should be directed to the ACS).

FIGURE 4
FIGURE 4Comparison of the proton delivery pathways for the A-and B-types of CcO. Green spheres are water molecules.Reproduced with permission from Figure3ofChang et al. (2009).
prepared for the next redox, proton pumping, and oxygen binding cycle.(In the "fully reduced enzyme," all metal sites are reduced, Fe 2+ or Cu + .)From the viewpoint of the DNC, additional electrons, where needed, are supplied directly by Fe b (in the ba 3 enzyme), and the Cu A dimer then reduces Fe b for the next 1e − donation to Fe 3+ a3 (or to Fe 4+ a3 ) or Cu 2+ B .The Cu A cluster and Fe b are resupplied by electrons from Cyt c as needed, but always in single-electron transfer steps.The DNC catalytic reaction cycle is depicted in Figure5

FIGURE 5
FIGURE 5Intermediate states of the DNC in the catalytic cycle of CcO.States A, P M /P R , F, and O H were identified by resonance Raman experiments on A-family oxidases.Although the peroxo, hydroperoxo, and F H states (red frames) were not observed experimentally, they may exist for a short time.Tyr is the special Y288 in aa 3 and Y237 in ba 3 enzymes.Reproduced with permission from Figure6ofNoodleman et al. (2020).
FIGURE 6 (A) Simplified map of the intermediate states of the DNC in the catalytic cycle of CcO shown in Figure 5. Reproduced from Figure 1 in Han Du et al. (2019), copyright 2019 American Chemical Society (further permissions for reusing this figure should be directed to the ACS).(B) Comparison of Fe a3 -O/ O-O vibrational frequencies for states along the reaction cycle of CcO obtained from DFT calculations on DNC cluster models and from resonance Raman experiments.

FIGURE 7
FIGURE 7 Correlation between the calculated Fe a3 -O bond lengths (r Fe-O , Å) and the corresponding Fe a3 -O stretch frequencies (1/] 2/3 Fe−O , based on Badger's rule) for eight DNC state structures.Reproduced from Figure 7 in Han Du et al. (2019), copyright 2019 American Chemical Society (further permissions for reusing this figure should be directed to the ACS).

FIGURE 8
FIGURE 8 Proton uptake/pumping processes during the P R →F →O H transition.The two water molecules are HOH604 and HOH608 from Figure 3B.The ΔG (kcal/mol) change at pH = 7 (ΔG =1.37 (pK a −7) between two different protonation states and the energy difference corrected by zero point energy between two tautomers are shown.The indices (a through h) label the various redox, protonation, and structural states as discussed in the text.Reproduced from Figure 4 in Han Du et al. -2018), copyright 2018 American Chemical Society (further permissions for reusing this figure should be directed to the ACS).

FIGURE 9
FIGURE 9 Feasible proton transfer pathways in the proton loading network (PLN) based on molecular dynamics simulations of the hydroperoxo state.Protonation state D is lowest in energy.Transferred protons and the ε-nitrogen on the imidazole ring of His376 are colored to guide the eye.(A) A strong hydrogen bond exists between protonated Asp372 and PRAa 3 , facilitating proton transfer.(B) Rotation of His376.(C) Proton uptake by Asp372 results in the formation of a strong hydrogen bond between protonated PRAa 3 and His376.(D) Proton transfer from PRAa 3 to His376 results in a strong hydrogen bond between His376 and Glu126 II .(E) Proton transfer from doubly protonated His376 to Glu126 II and subsequent release into the water pool.(F) Structure and protonation state equivalent to (A) with the exception of the orientation of His376 and protonation site δ-nitrogen.Reproduced from Figure 1 in Yang et al. (2016a), copyright 2016 Elsevier.
the Supporting Material there (states A′, B′, and C′ are alternate tautomers to States A, B, and C (protonated at N ε) on His376 (neutral).These alternate tautomers are expected to have very similar energies.During these dynamics simulations, in addition to fluctuating hydrogen bond linkages and water mobility, ring flips of His376 and movements of nearby sidechains are frequent for simulations from any of the starting States A-F (or States A', B', and C').
FIGURE 10 (A): DFT energies for proton transfer and proton uptake from the N-side of the membrane (A'/B/B2 = cluster model + H + at pH = 7) to His376 (D).The initial step is a proton transfer from Asp372 to PRAa 3 (B).Protonation of His376 can occur either before proton uptake (B2) or after proton uptake (C′).The labels correspond to the nomenclature in Figure 9.The prime indicates δ-protonation in place of ε-protonation of His376.(B): complete DNC cluster model for State A'.Reproduced from Figure 3 in Yang et al. (2016a), copyright 2016 Elsevier.

FIGURE 11
FIGURE 11 Water pool above the DNC in the crystal structure of ba 3 -type CcO from Thermus thermophilus (PDB ID 3S8F).(A): Location of the water pool within the enzyme with surrounding residues highlighted.Purple spheres are metals in the protein (Cu A , Cu B , and Fe a3 ).Subunits I, II, and III are shown in yellow, green, and purple, respectively.(B): Details of the water pool and surrounding residues.The orange region is the water pool consisting of eight crystal water molecules.Reproduced from Figure 5 in Yang et al. (2016a), copyright 2016 Elsevier.

FIGURE 12
FIGURE 12General view of the two water exit pathways that connect the water pool above the DNC to the P-side of the membrane.Residues in blue color line pathway P1, and residues in orange color line pathway P2.The positions of water molecules along one branch of each of the pathways are depicted with red spheres.Reproduced from Figure6inYang et al. (2016a), copyright 2016 Elsevier.
FIGURE 13 (A): Free energy profiles for the salt bridge Arg225-Asp287 in different protonation states of CcO.The distance is measured between the carbon atom C ζ of the guanidinium group in Arg225 and the carbon atom C γ of the carboxylate group in Asp287.(Dotted line) (B): Snapshot from a simulation of state E with protonated Glu126 II , highlighting relevant residues.Reproduced from Figure 4 in ref.Yang et al. (2016a), copyright 2016 Elsevier.

FIGURE 14
FIGURE 14Major residues and water molecules in the DNC of the Tt ba 3 X-ray crystal structure 3S8G.For clarity, only Prop-A and Prop-D of the heme ring are presented here; the whole heme ring is not shown.Reproduced with permission from Figure15ofNoodleman et al. (2020).

1 )
DFT calculations find structures and vibrational frequencies along a reaction pathway that closely correlate with Fe-O frequencies observed by resonance Raman spectroscopy.Calculated Fe-O and O-O Mayer bond orders provide further insights into bonding and associated changes in geometries.2) We develop the water shift mechanism for proton pumping by combining DFT and molecular dynamics work.Two water molecules per cycle are reaction products, but these two water molecules (or equivalently other water molecules that exchange with them) shift their positions when both Cu(II)-OH − and Tyr(O) − within the dinuclear center (DNC) are protonated, descreening and affecting histidine deprotonation within the proton loading network.Small pools of water lead to water exit channels, which in combination with protein sidechains, act as proton exit carriers.Both long-range electrostatic effects (Coulombic interactions) and short-range chemical bonding with charge transfer (for example, Tyr → His → Cu) and proton/electron transfers are critical for all these mechanisms.3) Cu(I)-histidine bond lability is probably critical to histidine protonation, activating the rapid electron transfer/proton transfer reaction (eT/pT):O − 2 + 1e − + 1H + → OOH − Fe 3+ -Superoxide → Fe 3+ -Hydroperoxide, with 1e − from Cu + )4) We propose additional mechanisms for transfers of chemical protons into the reaction chamber of the dinuclear center for the proton-coupled electron transfer processes state F → O H and state O H → E H .We propose a pumping mechanism in ba 3 CcO for a single-vector proton via a radical transfer mechanism during the State P M + 1e − → P R transition similar to the one proposed byYu et al. (2012) in aa 3 -type CcO enzymes.5) In Complex 3, the splitting (bifurcation) of electron flow between the high potential upper pathway and the low potential lower pathway is controlled by a coupled proton transfer from reduced ubiquinone to the Rieske [2Fe-2S] cluster.In Complex 4, the split pathways (bifurcations) of proton flow depend both on electron transfers, for example, from Cu 2+ → Cu + (or conversely), on other electron transfers (one electron transfer to superoxide among others), and on the interactions of chemical and vector protons.These comparisons illustrate common themes between Complex 3 and Complex 4 mechanisms, while substantial differences exist in the structures of the metal complexes and their functions, the types of mobile carriers of electrons and protons, and the proton pumping mechanisms.6) There are many important open questions.a.Chemical reaction paths and associated intermediate energies and reaction barriers are incomplete either from an experimental or a quantum chemical/molecular dynamics perspective.b.Proton-loading network (PLN) mechanisms and water protein pathways for proton exits are partly developed.c.What additional supporting experimental and quantum chemical information can be developed for radical transfer pathways of proton pumping?d.Where are reaction mechanisms and energies similar or different in aa 3 -type versus ba 3 -type cytochrome c oxidases?e.How efficient is the reaction cycle in practice while pumping 1.5+-Cu 1.5+ , Fe b is in the Fe 2+ state, 1e − reduced from the oxidized Fe 3+ state, and the DNC is also reduced,

TABLE 1
Proton and electron transport pathways and mechanisms in CcO.Reproduced from Noodleman et al. (2020).Differences in proton channels and number and types of transported protons in bacterial and mitochondrial CcO.
NO, and O 2 to two model heme systems.OLYP and BP86 best reproduced the binding energies of CASSCF/CASPT2.Potentials BP86 and PW91 often give good binding energies, but over-stabilize low-spin states for Fe and other transition metal complexes, while OLYP and OPBE are more balanced in providing accurate relative spin state energies.Since the reaction cycle for cytochrome c oxidase involves both oxygen species binding, Fe and Cu spin state changes, redox changes, and protonation/deprotonation, we have emphasized OLYP in our more recent work published in 2014 and later.