Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biochim Biophys Acta. Author manuscript; available in PMC 2012 April 1.
Published in final edited form as:
PMCID: PMC3055932

Exploration of the Cytochrome c Oxidase Pathway Puzzle and Examination of the Origin of Elusive Mutational Effects


Gaining detailed understanding of the energetics of the proton-pumping process in cytochrome c oxidase (CcO) is a problem of great current interest. Despite promising mechanistic proposals, so far, a physically consistent model that would reproduce all the relevant barriers needed to create a working pump has not been presented. In addition, there are major problems in elucidating the origin of key mutational effects and in understanding the nature of the apparent pKa values associated with the pH dependencies of specific proton transfer (PT) reactions in CcO. This work takes a key step in resolving the above problems, by considering mutations, such as the Asn139Asp replacement, that blocks proton pumping without affecting PT to the catalytic site. We first introduce a formulation that makes it possible to relate the apparent pKa of Glu286 to different conformational states of this residue. We then use the new formulation along with the calculated pKa values of Glu286 at these different conformations to reproduce the experimentally observed apparent pKa of the residue. Next, we take the X-ray structures of the native and Asn139Asp mutant of the Paracoccus denitrificans CcO (N131D in this system) and reproduce for the first time the change in the primary PT pathways (and other key features) based on simulations that start with the observed structural changes. We also consider the competition between proton transport to the catalytic site and the pump site, as a function of the bulk pH, as well as the H/D isotope effect, and use this information to explore the relative height of the two barriers. The paper emphasizes the crucial role of energy-based considerations that include the PT process, and the delicate control of PT in CcO.

Keywords: Cytochrome c oxidase, proton pump, apparent pKa, mutational effects, energetics of proton transfer, isotope effect, electrostatic energy

1. Introduction

Cytochrome c oxidase (CcO) couples the four electron reduction of O2 to water and trans-membrane proton pumping from the negative (n) to the positive (p) side of the membrane (e.g.), which results in an electrochemical proton gradient that drives, for example, ATP synthesis. The elucidation of the structure of CcO combined with experimental and theoretical investigations of the wild-type and mutant forms of the enzyme (e.g., for reviews see refs.) has offered an opportunity to analyze, on a molecular level, this intriguing molecular machine. Unfortunately, the exact details of the action of CcO continue to present an extremely challenging problem and we still do not have a clear idea of the nature and the factors that control the pumping process and block the leaking processes. This is further complicated by the fact that for each reaction that is observed experimentally many events occur with the same rate. For example, in the last step of the reaction of the four-electron reduced CcO with O2 (ferryl (F)-oxidized (O) reaction) there is electron transfer to the catalytic site, proton transfer (PT) to the catalytic site through the so-called D proton-transfer pathway, uptake of a pumped proton through the same pathway to an internal acceptor and release of the pumped proton. All these reactions occur with an observed time constant of ~1 ms at pH 7.5. Consequently, it is difficult to investigate specific single events, which would be required to be able to address the overall problem of proton pumping using theoretical tools.

A significant progress has been made in defining the conditions that would allow CcO to pump protons against a pH gradient, in assessing the electrostatic energy of different possible intermediates, in examining the energetics of the key water chains and of some non-rate limiting PT steps. Furthermore, examinations of the energetics of the overall pumping process have been performed using a semi-macroscopic model. However, the crucial relationship between the protein structure, the PT energetics and directionality has not been established and we are not aware of any study so far that could generate a proton pump purely based on the structure and well-defined energy considerations. More specifically, theoretical evaluations of the pump features should determine the barriers for the PT steps when starting from the observed structure (e.g.). Despite major progress we still cannot establish the exact sequence of events or the relevant energetics of the states involved. For example, the identity of the primary acceptor for pumped protons is not known, even though the D propionate of heme a3 (Prd in Fig. 1) is one of the most likely candidate for being at least a “transient” acceptor. There has been some progress to use structures in functional analyses, but there are still problems and so far it has not been possible to reconstruct the whole reaction cycle from structure-based energy calculations. In this respect we would like to emphasize that although phenomenological attempts to assign the energetics of key states can provide very useful insights, we do not consider such approaches as structure-function correlation studies.

Figure 1
The key proton pathways in CcO at the end of the D pathway. The notations used in this figure are taken from ref.. Here B designates the Fe-bound OH and the Wi (i=1,2,3,4) designate the water molecules that may be involved in the PT ...

In resolving many of the above-discussed problems, one may greatly benefit from focusing on a simpler problem and studying mutant forms of the CcO, which do not pump protons (referred to as uncoupled mutants), but in which the rates of intra-molecular proton-transfer reactions are unaffected. In these cases the PT parameters are almost known. On one hand, studies of the uncoupled mutants are blind to the PT pathway to the acceptor site for pumped protons (Propionate D or A or some other site, see Figure 1); yet they provide much clearer and unique information related to PT through the Dpathway and to the Fe bound OH at the binuclear catalytic site (this specific OH will be referred to as B in the subsequent text. See Fig. 1 for clarification).

It must be mentioned that in other systems such as bacteriorhodopsin (bR) it has been possible to rationalize the energetics of the pumping process and its coupling to the pKa changes quite early. Furthermore, the initial ideas of using the light induced charge separation to change pKas and then for inducing irreversible conformational changes turned out to be consistent with detailed structure-based calculations. Similarly it is clear how light energy is used to drive an irreversible charge separation in photosynthesis. However, the situation is much less clear in CcO.

The present paper focuses on structure based analysis of the nature of the PT in the N139D mutant (N131D in the P. denitrificans CcO) and relates it to the apparent pKa of Glu286 (Glu278 in P. denitrificans) in different conformations.

2. Background

The action of CcO (under steady-state turnover conditions) involves steps where electrons are transferred from cytochrome c to CuA and then consecutively to heme a and the catalytic site composed of heme a3 and CuB. Each such electron transfer is accompanied by proton uptake to the catalytic site (B in Fig. 1) as well as the uptake (to Prd) and release of a proton that is pumped across the membrane. Thus, each step of the reaction involves both electron transfer through the CcO as well as the uptake of two protons from the n-side and release of a proton to the p-side. An alternative approach that allows studies of single specific reaction steps in the time-resolved domain involves preparation of the “fully reduced” CcO in which each of the redox centers, CuA, heme a, CuB and heme a3 are reduced. A blocking carbon monoxide ligand is then bound to heme a3 at the catalytic site. After addition of O2 the CO ligand is removed by means of a short laser flash, which allows O2 to bind (with a time constant of ~10 μs at 1 mM O2) and react at the catalytic site. In the first step an electron is transferred from heme a to the catalytic site (which now carries 3 electrons) with a time constant of ~30μs forming a state that is called PR and in which the O-O bond is broken. This state is unique in that its formation is not linked to any proton uptake from solution and the reaction only involves local proton movement within the catalytic site. Consequently, the PR state carries one negative charge in excess to the other states. In the next step the excess negative charge at the catalytic site is neutralized by proton uptake from solution to the catalytic site forming state F. The PRF reaction does not involve any electron transfer within the CcO. However, it involves the uptake of a substrate proton to the catalytic site (B−) as well as the uptake (to Prd) and release of a proton that is pumped across the membrane, all displaying a time constant of ~100μs at neutral pH. In the last step of the reaction the fourth electron is transferred to the catalytic site to form the oxidized CcO, accompanied by proton uptake to the catalytic site and the uptake and release of a proton that is pumped across the membrane, all displaying a time constant of ~1ms at neutral pH.

A generic model showing electron and proton-transfer reactions that occur in a typical cycle of CcO (see above) is described in Figs. 1 and and2.2. The primary event involves an electron transfer to heme a and subsequent PT from E286 to Prd, followed by a reprotonation of E286 and then a PT to the catalytic site (binuclear center, B). Here, a particularly interesting issue is the role of E286. Possible conformational as well as ionization states of this residue have been considered by several workers. Several of these studies assumed that this residue can serve as a valve (e.g.), but this may violate microscopic reversibility as a true valve must involve changing the barriers between alternative paths (see discussion and the results presented in, as well as the concluding remarks of this work). Studies of the ionization state of E286 involve more consensus but also some misunderstandings. For example, it has been suggested that the work of Olsson et al. concluded that E286 cannot be deprotonated while results from this research group have reproduced the observed E286 pKa since the first publication. The possible confusion reflects the fact that ref. considered a mechanism where the first proton was assumed to be transferred from the bulk solution. This assumed mechanism gave reasonable energy only with a concerted PT. However, a subsequent detailed consideration of the energetics of the different options could also provide a physically based rationalization for a mechanism where the primary proton donor is E286, as was proposed on the basis of early experimental studies. Of course, none of these studies argued that E286 cannot be ionized at high pH.

Figure 2
(a) A schematic depiction of the reaction steps in CcO. The numbers indicate the order of reactions. Due to space limitation, we do not draw each step separately and thus the charges correspond to the first step in each figure. Here, we designate a reduced ...

At any rate, as outlined in great detail in it is extremely difficult to consistently account for proton pumping and for the blockage of proton leaks while using energy-based calculations. Even the recent study presented in did not provide the energetics for the entire pumping process and left open questions related to the nature of the entire mechanism, and probably even the identity of Prd as the primary acceptor. In fact, there are no experimental data directly showing that Prd is involved in proton pumping (see e.g. ref.), although the Prd is one likely candidate. Thus, as we stated in the Introduction section the focus on mutations where PT to the Prd is blocked may provide a major simplification that should allow one to have much more confidence in the detailed analysis.

More specifically, we simplify the problem in two ways: First, we focus on a mutant CcO in which the catalytic reaction is uncoupled from proton pumping so that step 2 in Fig. 2 does not occur. Second, we focus on a simple PT through the D pathway to the catalytic site, where deprotonation and reprotonation of E286 is studied without an accompanying ET, which would complicate the conclusions.

3. Experimental studies of the proton-transfer event

Detailed studies of the pH dependence of the PRF reaction with the wild-type CcO showed that the rate titrates according to a Henderson-Hasselbalch equation with a pKa of 9.4, which is attributed to the titration of E286:

equation M1

where α EH is the fraction of protonated E286, pKE286 is its pKa and kH is the proton-transfer rate from E286 to the OH in the catalytic site (B) estimated to be 1.1×104 s−1. Thus the model is considered with the assumption that the internal PT from E286 to the catalytic site displays a rate constant of 1.1×104 s−1 and that the proton equilibrium between solution and E286 is rapid (>104 s−1).

Despite the possibility to model the PRF reaction using a simple titration model, as already mentioned above, the reaction involves transfer of both the substrate proton to the catalytic site and the transfer of a proton that is pumped across the membrane. In other words, two protons are transferred simultaneously through the Dpathway. As already mentioned above, to further simplify the process, we consider mutant CcOs in which the O2-reduction reaction is uncoupled from proton pumping, i.e. the PRF reaction only involves transfer of a substrate proton to the catalytic site (B−).

Studies of the pH dependence of the PRF reaction in a number of such mutant CcOs showed that the apparent E286 pKa is typically shifted to lower or higher values (Table 1) (reviewed in ref.). This is independent whether or not charged residues are introduced as a result of the mutation. The present work is focused on the N139D mutant CcO because it is one of the rare cases where we have structural information. Even though the structure was determined for the CcO from P. denitrificans, that enzyme displays essentially the same functional characteristics as the R. sphaeroides CcO counterpart. Furthermore, in this mutant CcO of the R. sphaeroides CcO the PRF reaction displays the same PT rate as with the wild-type CcO, but it titrates with a higher pKa of ≥11.

Table 1
Kinetic parameters for the wild-type (WT), N139D and N139T mutant CcOs.

We start our analysis by considering the fact that the apparent pKa is more complex than what is frequently assumed as it does not reflect the true pKa of the given residue at the specified state. To deal with this issue we start with the tentative model of (which will be refined below). In this model we consider two positions of E286 denoted by E1 and E2 where in each of these states E286 can be either deprotonated (E1, E2) or protonated (E1H, E2H). The pKas of E286 are different in the two states, namely pKa1 and pKa2. PT to the catalytic site can only take place when the E286 is protonated and we assume that the PT takes place only from state E2, i.e. E2H. The scenario is summarized in Fig. 3. Representative structures for the two positions of E286 are shown in Fig. 7 as found in the X-ray structure of the N139D mutant.

Figure 3
A schematic diagram for the states that are involved in the determination pKapp in mutants where the transfer to Prd is blocked. Here, PT designates proton transfer to B. E1H denotes the protonated E286 in the downward configuration and E2H denotes ...
Figure 7
The two configurations of E286 as found in the crystal structure of the N139D mutant (PDB code: 3EHB) in P. denitrificans CcO. They are named as follows. Configuration E1: downward E286 found in the crystal structures of both native and mutant enzymes; ...

Using the diagram of Fig. 3 we define the relevant equilibrium constants by:

equation M2

We also assume a rapid equilibrium between the different states so that the observed proton-transfer rate from E286 to the catalytic site, kH,CS, is proportional to the fraction CcO in state E2H, [E2H]:

equation M3

where, kH,max is a constant, which is equal to the proton-transfer rate when [E2H] equals unity.

For simplicity, we assume first that the above-discussed mutations within the Dpathway have no effect on the pKa values pKa1 and pKa2 and only lead to changes in the equilibrium constants between states E1 and E2. In reality the mutations are likely to affect also the pKa values, for example in the case of the N139D mutant the electrostatic interaction between the ionized D139 would shift the pKa of E286 by a small amount. Furthermore, the mutations can change the distribution of the water molecules in the D pathway and thus change their ability to solvate to ionized E286.

The apparent pKa value, pKapp (defined as [E2H]max/2) is given by:

equation M4

Using this expression allows us to explore the effect of the local environment on the apparent pKa and to try (with the help of computational analysis) to assess the energetics of different configurations of E286. The above relationships provide major restrictions on the possible values of K1 and K2 for a given Kapp (see Figure 4). In addition, the model can provide an effective experimentally-based free energy as is illustrated in Figure 5 for a set of assumed pKa1 and pKa2. However, as is demonstrated below, a more unique analysis must involve the use of computational approaches that can relate the pKa values to the specific environment and also to estimate K1 and K2.

Figure 4
Regions of allowed pKa1 and pKa2 for feasible ranges of K1 and experimentally known pKapp, based on Eq. 6.
Figure 5
The energetics of the states depicted in Fig. 3 based on a representative set of computed values of K1, K2 and pKapp as presented in Table 4. The computed values are for the native CcO: K1=0.2, K2=3.0×10−5, pKa1 = 9.4, pKa2 = 13.2, pK ...

It is also useful to introduce at this step the qualitative considerations of the relationship between the pKas of the two structures and the mechanism of CcO. This is done in Figure 6 for the tentatively assumed values (Figure 5) and then considered for the case when the transfer to Prd is blocked. These figures are only given as an introduction to guide the reader and a more quantitative description is given below using structural information.

Figure 6
A tentative description of the energetics of the first few reaction steps in the native and the N139D CcOs. A more quantitative description is given in Fig. 9 (where the state notations are given in a different and more general way). Here E, P, B ...

Finally, we note that the apparent pKa may reflect the competition between different PT paths and this issue will also be considered in the present work.

4. Theoretical Methods

In order to explore the ideas presented in the previous section it is crucial to find a relation between structures and, pKas and free energies. This was done in our work by using models that were verified in previous studies of PT reactions in proteins. These include the empirical valence bond (EVB) and the free energy perturbation/adiabatic charging (FEP/AC) microscopic models (both with a polarizable force field) and the protein dipole Langevin dipole semi-macroscopic version (the PDLD/S-LRA model). All of these models are described extensively elsewhere and thus they are not discussed in detail here.

The simulation system was constructed by starting from the X-ray structure of CcO from P. denitrificans [PDB entry: 3HB3 for native and 3EHB for mutant]. The protein was surrounded by a simplified polarizable membrane, as was done in.

The PDLD/S-LRA calculations involved two steps (e.g., references in): first, running molecular-dynamics simulations (MD) to generate protein configurations for the charged and uncharged states, and then averaging the PDLD/S results for the generated configurations. The MD runs were performed with the polarizable ENZYMIX force field with a 0.5 fs time step. All the PDLD/S-LRA calculations were performed by the automated procedure of the MOLARIS program where we generated 30 configurations for the charged and uncharged states using MD simulations at every 2 ps of the run. The microscopic FEP and EVB/US calculations were performed using the ENZYMIX force field, with the solute parameters described in Refs. . The simulation included the use of 22 Å of the SCAAS spherical constraints and the local reaction field (LRF) long-range treatment (see). The simulation system represented the membrane by a grid of induced dipoles (e.g., see ref.) which were treated explicitly in our polarizable model.

In view of the difficulties of obtaining converged results in the CcO interior (see) we found it useful to apply a hybrid of the EVB and the PDLD/S-LRA methods. In this approach we performed the EVB calculations first, with 31 frames with each of them being 20 ps long, for the transfer between each step in the given PT process. Next, we adjusted the EVB minima (by changing the gas phase shift) forcing them to reproduce the PDLD/S-LRA energetics. This type of strategy has been found to be very effective in our previous studies of PTR in proteins (e.g.).

The free energy of rotating E286 was evaluated by free energy perturbation (FEP) method. In addition to the original configuration (Configuration E1 in Fig. 7) of the side chain, a replica of the target geometry (Configurations E2/E3/E4 in Figs. 7 and and8)8) was constructed and constrained to respective configurations. Atoms in the original configuration were mutated during the FEP run to dummy atoms (no non-bonded interaction with the environment) starting from real atoms (with full interactions), whereas atoms in the target geometry were mutated from dummy atoms to real atoms. This whole process corresponds to moving the side chain from the original configuration to the target geometry. These calculations were done within MOLARIS using FEP method with the whole run extending over 31 frames in the forward direction, with each frame being 20 ps long.

Figure 8
The rotated E286 configurations generated by us. Configuration E3: within 6Å from B; Configuration E4: within 6Å from Prd.

In all of the calculations we determined the ionization states of the cofactors according to the given step in the cycle. The charges of the protein ionizable groups in the neighborhood of the relevant groups (those considered in Fig. 2) were determined by PDLD/S-LRA calculations and those of distant groups were set to the corresponding values in water at pH=7. The interaction between the ionizable groups and the distant groups were evaluated with εeff=30. The interaction with closer residues were evaluated with εeff=20. The considerations in using εeff are discussed in detail for example in Refs. .

As will be pointed out in the subsequent sections, we are dealing with a very complex system and must consider different kinetic options. Thus, we have found it essential to formulate the kinetics of the system in terms of the following master equations (see Ref. for a related study) based on the competing pathways as described in Fig. 9. The actual simulations involved the detailed kinetic diagram of Fig. 10, where the relevant rates are given by:

equation M5
Figure 9
A schematic diagram for the energetics of the early competing steps in CcO (based on ref. and subsequent estimates) this diagram does not provide definitive conclusive results but rather it is given to clarify the nature of the competing paths. The paths ...
Figure 10
The subset of kinetic pathways shown in Fig. 9 that has been used for our simulated kinetics study. Here all species names correspond to the different states as described in Fig. 9. The species I can transfer a proton to a site in the D pathway, which ...

Here, all species names correspond to the different states described in Figs. 9 and and10.10. [A] indicates the concentration of the species “A” at a given time. We have allowed species I and VI′ to remain in proton exchange equilibrium with the bulk (D channel). Species I′ (not shown in Fig. 10) designates the deprotonated form of I after it releases a proton to the D channel as shown in Figure 12. This is assumed to be at equilibrium (before the initial electron transfer to a) and the initial population was determined by the relation equation M6 at all pH values, where equation M7 is the pKa of the state I. We have used equation M8 and equation M9. The rate expressions for PT from/to the bulk solution have been developed earlier in ref. and shown in Fig. 10. The proton-transfer rate for moving to the D channel is given by: equation M10, where the estimated barrier for PT from the bulk to the lowest energy point (called equation M11 here) in the D channel is 2 kcal/mol. Similarly, the proton-transfer rates from the D channel to the bulk is given by: R−b = [A2H]k−b = [A2H]×6×1012 exp[− {2+1.36(pKa(A2H)−(−1))}/RT], where pKa of a bulk H3O+ is considered to be −1. Here we take pKa(A2H) as 4, but any number in the range of 0 to 5 would give similar result. The rate constant for moving from A2H to A1H (our E286) is given by kD = 6×1012 exp{−6/RT} and the for the reverse process it is k−D = 6×1012 exp[−{6+1.36(pka(A1H)−pka(A2H))}/RT]. Here pka (A1H) is taken as 11 and 8 for the species I and VI′, respectively while pKa (A2H) is assumed to be larger than zero. Under these conditions we can assume that the kinetics of moving the proton from the bulk to the equation M12 follows the Michaelis-Menten kinetics as was assumed in Ref., however this assumption was not used here.

Figure 12
Simulated kinetics for the pathways and rate constants given in Fig. 10. The simulations show the concentration of relevant species as a function of time for (a) native: k1=2×104, k−1=105, k2=103, k−2=102, k3=2×104, k−3 ...

5. Results and Discussion

5.1 Exploring the relationship between the apparent pKa and the orientations of E286

Following the discussion in Section 3, we attempted to move to a structure-energy description that integrates the model outlined in Fig. 3 with the energetics of CcO. In doing so we explored the structural changes of E286, keeping in mind that one of the structures should have a high barrier for PT to B. The obvious candidates are the structures of the native and the N139D mutant CcO, which are depicted in Fig. 7 (note that this mutant has two positions for E286). However, these two structures involve only a small rotation around the Cα-Cβ and Cβ-Cγ bonds, and the resulting change in distance to B is unlikely to be enough to change the PT reaction barrier significantly. Thus, we also generated another structure with much more significant change in position by about 120° rotation around the Cα-Cβ bond such that the COO group is within about 6Å distance from B (See Fig. 8). The structural parameters of the two models are given in Table 2. We also evaluated the potential of mean force (PMF) for moving between the different configurations and the corresponding results are given in Table 3. Now, as seen from Table 4 the PMF-based energetics of the “small rotation” (configuration E2) of the mutant is close to zero and thus consistent with the observation of two configurations in the X-ray study.

Table 2
The structural parameters of the four configurations of E286 (shown graphically in Figs. 7 and and88)a
Table 3
The energetics of rotating E286 to different configurationsa
Table 4
Calculated pKa values for different configurations and the corresponding pKapp, K1 and K2 valuesa

With the above structural models we evaluated the pKa of E286 for both the native and mutant forms in the two configurations of E286: Config. E1 and Config. E2 (for structural parameters see Table 2 and Figs. 78). The corresponding results for pKa calculations are given in Table 3. Using the calculated pKas and the observed pKapp we can also examine the range of reasonable K1 and K2 and this is done in Fig. 4, where we describe the regions of possible solutions of Eq. 6. As seen from Table 4 we were able to obtain a reasonable approximation for the observed pKapp with reasonable K1 and K2 values. Note in this respect that the range of the dielectric for the self-energy (εp) used here is within the range allowed with the PDLD/S-LRA (which includes a full relaxation of the protein).

The calculations of the pKas of the native and mutant were done for different configurations. That is, the pKa has been calculated for state VI′ in the native (where Prd is protonated), and for state XI for the mutant (see Fig. 10 for clarification of the nature of the different states). Note that the treatment is completely consistent with the physics of the system, since the pKapp reflects different path in each case.

It is important to note at this point that the small value of εp does not mean that the effective dielectric for charge-charge interaction, εeff is small. The relevant considerations of εp and εeff (which is typically around 20 and very rarely or never can be smaller than 8) are given elsewhere (e.g. see refs. and). This is important in view of the fact that some workers (e.g. ref.) have assumed εeff values smaller than 5, which is hard to justify by models that were used previously to reproduce experimental observations in other proteins (e.g. ref.). In this respect it would be important to explore the meaning of the low apparent dielectric constants deduced in ref. .

Of course, the real question is whether the tentative analysis of Fig. 4 and Table 4 is consistent with the actual energetics of the system. That is, a consistent analysis requires that the K1 value obtained from the PMF of the real system would result in a pKapp that is close to the corresponding observed value. Indeed the results given in Table 4 show a very reasonable consistency between the pKapp obtained with the calculated pKa and calculated K1 values. For example, we find that for εp=3.5 the calculated values for the native and mutant systems are 9.4 and 10.8, whereas the experimentally observed values are 9.4 and >11, respectively.

5.2 Exploring the competition between the two proton-transfer paths

Another interesting issue that is explored here is the rate of the PT from E286 to B. As mentioned above, in the case of the N139D mutant CcO we know with confidence that the experiments reflect this specific PT rate and thus we are able to explore the reliability of our modeling. The calculated activation barriers obtained by combining the PMF and the PDLD/S-EVB energies are summarized in Table 5 and Fig 11(a), for cases with one or two bridging water molecules. We find that the observed low barrier of about 13 kcal/mol is reproduced with the bridging two water molecules, whereas the one bridging water molecule has a higher barrier because of the penalty associated with the large rotation of E286.

Figure 11
a) The barriers for PT from E286 to B in the native (solid line) and mutant (dashed line) enzymes. b) The barriers for PT from E286 to Prd in the native (solid line) and mutant (dashed line) enzymes. Two alternatives are shown in each panel; ...
Table 5
The energetics of two competing paths for PT to B−:a

Although our initial focus was placed on the study of the path for PT from E286 to B, we realized that the present study could provide interesting insights into the nature of the illusive mutational effect. That is, while we know from the experimental observations that the barrier for PT to B is not modified significantly in the N139D mutant (the PT rate to B− is the same as in the wild-type CcO), the effect of the mutation has not been elucidated. Here, it is possible that it is due to an increase in the barrier for PT to Prd. However, it is not entirely clear that the primary event is indeed a PT to Prd, and it is clearly not obvious why proton pumping is impaired in the N139D mutant CcO. Thus, we attempted to explore the change in the PT barrier as a result of the N139D mutation. In doing so we combined the PMF and PDLD/S-LRA calculations and the results are summarized in Tables 5 and and66 as well as Fig. 11. The energy terms considered in these tables are similar to those considered in ref., where the reference energy for the given PT in water ( equation M13) is given by:

equation M14

where, equation M15 and equation M16 are the pKa values of respective species in water, and ΔGqq is the charge-charge interaction of A and H3O+ at the distance where PT occurs. Here we have used Δ Gqq =3.9 at a distance of 3Å and εeff=30.

Table 6
The energetics of competing paths for PT to Prd−:a

The results of our analysis (which are given for the cases with a bridge of one water molecule and two water molecule) provide interesting insights into the nature of the effects of mutations in the D pathway and in particular in the N139D mutant CcO. That is, as seen from Fig. 11(b) our calculations for the PT to Prd show almost no change in the energy profile for the single water bridge while the mutation results in an appreciable increase in free energy in the case of the two water bridge. If the two-water molecules path provides a lower barrier than the single water molecule, we might have a consistent explanation for the blockage of PT to Prd in the mutant. However, we also found that the barrier for PT to B undergoes a decrease upon mutation. Thus, with the current results we have two options: (i) The barrier for PT from E286 to Prd increases and the barrier for PT to B decreases (or stays unchanged), and (ii) the overall barrier for PT to Prd does not change significantly (if the barrier with two water molecules is not lower than that with one water molecule) and the barier for protonating B is reduced. Note that in both cases we have to satisfy the experimntal observation that the overall rate for protonating B does not change with the mutation. The first possibility provides a simple explanation for the mutational effect, where effectively the path to Prd is blocked (outcompeted). The second option is more challenging (as it seems to lead to faster rate of forming B in the mutant) and we need to explore the kinetics of the system in a more complete way, considering the possibility that a fraction of protons that goes to Prd returns to E286 and joins the fraction that went directly from E286 to B, as well as other kinetic effects which are hard to deduce by just looking at the rate constants. Thus we felt that it is essential to have a quantitative understanding of the kinetic complexity and applied the master equation treatment as described in Section 4. This was done by integrating the equations 7–16, while using the rate constants described in Fig. 10 and a time step of 10−12. We have used 3 sets of rate constants (described in Fig. 12) to reproduce the native (Fig. 12(a)) and mutant (Fig. 12(b–c)) systems. Note that the Fig. 12(b) does reproduce the observed trend of not pumping, still having a reasonably small change in the time it takes to reach 50% population of BH (which determines the apparent rate constant). The overall formation rate of BH (combination of XV′, XIV′ and X) was comparable for both the native and mutant CcOs above t=10−4 (~100μs), although it was impossible to satisfy this condition at very short times. In order to test for option (ii) mentioned above, we have considered a case where the barrier to Prd (and the rate constant k1) was kept the same in both cases, whereas the barrier to B (and the rate constant k2) was lowered for the mutant (Fig. 2(c)). In this case, it was hard to prevent pumping without a large increase in k2, and it was also required to suppress the final step by reducing k4 (X→X′) to in the range of 2ms, which seems to be too slow. Obviously option (ii) is more exotic but it should also be considered. It should also be pointed out that our calculated rate constant is not exactly the rate constant determined experimentally by fitting an exponential function, A = A0exp (k*t), to the observed time dependence, since a significant part of the observed apparent rate constant reflects the change in A0 (e.g. reduction in the amount of protonated E286 at high pH). Furthermore, the current approach seems to give the most stable results. Thus, we prefer to leave the current “effective rate constant” until we reach more unique conclusions in analyzing the kinetics of the reprotonation step in the mutant (where we are sure about the available pathways).

Another issue that could be explored by our approach is the effect of having two PT paths on the apparent pKa. This complex issue can be examined by changing the barriers and looking for the rate of populating BH. Typical simulation results are presented in Fig. 13. As seen from the simulations the apparent pKa may also reflect the competition between the two paths. That is, the case when the Prd path is blocked the apparent pKa reflects the protonation of E286 in state I (which is discussed above). However, when the Prd path is open the apparent pKa can also reflect the reprotonation of E286 in state XIII″. Here we have a lower pKa in state VI′ due to the change in electrostatic interaction between E286 and Prd (where Prd is either protonated or its proton is in another close by pump site). That is, now we have a case where the pKa of E286 in state XIII″ of Fig. 9 (that corresponds to PT to Prd) is ~8 and since this state is accessible at the lower pH we have a reduced apparent pKa. Here again we have used the master equations emphasizing the pH dependence and the fact that the population of AH is determined at equilibrium before the experimental initiation of the cycle. The actual results depend on the rate constants used and some typical results are shown in Fig. 13, where we considered the total rate for protonating B from both paths, with different ratios between the competing rates. Considering the simulation results it is tempting to assume that this trend is reflected in the difference between the native (where the path to Prd is open) and the mutant where this path is closed. However, here we have a complication since in the case of the native CcO the reprotonation of Glu 286 should depend on the bulk pH (a feature that is reflected in the simulations) and no such dependence is observed experimentally. This complication may also reflect the fact that there are 4 different PT reactions during Pr--F, where one is the rate-limiting step (100 micros) probably from Glu to B. However, at present our estimate of the rate for proton transfer from the bulk to A2 is slower than 10 s−1 at pH 8. This may force one to invoke additional complications with regard to the nature of the “bulk” protons, an issue that will be explored in subsequent studies. Despite this problem it seems that the pH profile does contains crucial information about the ratio between the barriers of PTR along the two paths.

Figure 13
Simulated kinetics for the rate of proton transfer to B as a function of pH for both the native and N139D mutant CcOs. Here the rate has been computed as [BH](t)/t = ([XV ′](t) + [XIV′](t) + [X](t))/t because XV′, XIV′ ...

A related piece of information has been obtained from studies of the deuterium isotope effect on the pumping efficiency of the N139D mutant CcO. Results from earlier studies have shown that in D2O the mutant CcO pumps protons with the average stoichiometry of 30% of that of the wild-type CcO. Here, we studied the pD dependence of the PRF reaction rate (Figure 14). As seen in the Figure, the maximum rate of the transition was a factor of ~1.6 slower in D2O than in H2O, which is approximately the same as that found for the wild-type CcO. The data also show that qualitatively the pKa associated with the reaction is lower, i.e. it moves towards that of the wild-type CytcO in D2O (9.6, ref.). It should be noted that the PRF in Figure 14 were obtained from fitting the kinetic data to a single exponential (see ref.). However, considering the fact that the N139D mutant CcO displays a 30 % pumping efficiency, it is unlikely that pKa change is fully manifested when simply fitting the data to a single exponential component (a more realistic evaluation should probably include a 30/70% distribution fractions displaying different kinetics or a 30/70% scaling of the pKa). Unfortunately, the signal-to-noise of the data did not allow for a meaningful distinction between different models. Nevertheless, qualitatively the data show a decrease in the effective pKa value.

Figure 14
The PRF rate as a function of pH/pD for the wild type and N139D mutant CcO in H2O and in D2O. The N139D data in H2O are from ref.. The wild-type data in D2O are from ref.. Experimental conditions were as in ref..

The observation that in deuterium N139D mutant CcO pumps with an efficiency 30% of that of the wild-type CcO may mean that in D2O we have almost equal barriers for PT to Prd and B. The posibility that the deuteration increases the activation barrier to B− more than the barrier to Prd indicates that the two barriers have different shapes (perhaps the difference is between having two water molecules or one water molecule). Another instructive point is the fact that the apparent pKa of the mutant was lower in D2O. Since the typical pKa changes upon deuteration are relatively small (ref.). This effect is likely to reflect the change in the paths, rather than the pKas in the different E286 configurations or the change in the intrinsic pKa due to deuteration.

5.3 Some Thoughts on Analyzing The Effect of Mutations

Mutational studies are extremely instructive in many cases, but it is important to keep in mind that the interpretation of functional effects of mutations is sometimes problematic. This issue has been previously addressed in papers of some of us (e.g.), but it is discussed again here. For example, one can claim that since N139D and N139T (see Table 1) have similar effect of blocking the E286 to Prd path, it is unlikely that N139D contributes to this effect by its electrostatic charge–charge interaction. The problem with this argument is that it is not deduced by any energy-based analysis of the PT process, but rather from the assumption that a comparison of N139D and N139T is based on the change in electrostatic charge-charge interaction. However, since CcO is very well tuned, many factors can lead to changes in the balance between the two competing paths. The effect of the electrostatic interaction between D139 and E286 is perhaps the most solid fact in the analysis of this system, since long range electrostatic effects are well understood and do not require the knowledge of the exact structural changes upon charging residue D139 (usually this results in compensating changes that reduce the effect of charge formation at site D139 on the charge at E286).

On the other hand, mutations can lead to structural changes around E286 and such effects can also lead to changing the balance between the PT paths. This effect is hard to quantify without structures of the native and mutant proteins (or at least calculated structural changes), and clearly cannot be interpreted without a model of the competing PT processes. Thus, we strongly argue that while mutation experiments provide a remarkable benchmark for modeling, their interpretation requires studies of the type reported here. Otherwise, the interpretation may not provide mechanistic information. Of course, if removing a given acidic residue results in changing the PT rate by several orders of magnitude, this group is most probably involved in the PT, but smaller effects from groups that are not directly involved are not so informative without very careful modeling studies.

6. Concluding Remarks

This work explored the molecular meaning of pKapp of E286 and provided a structure-based analysis of the effect of the N139D mutation on this pKa. The analysis of the apparent pKa involved a novel cycle that considers the conformational dependence of the pKa of E286 and relates it to the corresponding free energy of the conformational transition. Our analysis provides results that are simultaneously consistent with structure-based pKa calculations and with the PMF of the E286 conformational changes, both for the native and N139D mutant CcOs. Interestingly, the phenomenologically derived equilibrium constant for the conformational changes (needed to obtain the pKapp value), were found to be in a reasonable agreement with the corresponding results obtained using the PMF calculations. We are not aware of any previous study that has exploited the observed structural changes upon mutations in this way.

It is useful to note that a related pKa cycle has been considered in, however, that study did not consider the difference in PT rates for each structure, nor the relationship between the pKas of the two conformations and the pKapp. Furthermore, the treatment in Ref. involved basically hypothetical pKas rather than the pKas calculated using the structure of the given system.

The present work provides new insights into the nature of the effects of mutations in the D pathway and in particular the N139D mutant CcO. That is, our study reproduced an appreciable increase in the barrier of PT to Prd (but still not much larger than our error range) and a decrease in the barrier for PT to B. This finding presents an interesting dilemma: one feasible, but unlikely possibility is that Prd is not the primary proton acceptor. In case the rate for PT to Prd is unchanged then the reduction in the barrier for PT to B is the most likely reason for the inability of N139D to pump. In this case, we must also account for the fact that the overall rate for protonating B does not increase with the mutation. In order to account for this observation and still accepting the calculated reduction in the barrier for PT from Glu286 to B, we need to consider the possibility that the barrier for a direct transfer from Glu286 to B− is reduced in the mutant, but it required us to considerably suppress the last step in our kinetic scheme (X→X′). Still, the idea of a barrier to B being reduced in the mutant seems to be supported by one observed fact, namely the lower free energy cost for rotating E286 in the mutant (which is established by the availability of the two observed conformations, which are absent in the wt CcO). The finding that the availability of a rotated configuration of E286 appears to help in reducing the barrier is fully consistent with our previous study, where it was found that having a water chain is not so advantageous and that (in contrast to the belief of some) the PT rate is faster when we have as little water molecules as possible between the donor and acceptor. This, of course, depends on the probability that the donor and acceptor are at a close distance.

The present analysis of the effect of the N139D mutant CcO provides a more consistent picture than that obtained in our previous studies. That is, the study of was based on the assumption that the electron transfer step is on the same time range as that of the PT, which seems now unlikely. Our subsequent electrostatic calculations have rationalized the mutational effect in terms of the change of the effective electrostatic interaction with the transferred proton in both paths. This trend is reproduced in a much more quantitative way in the present work, where the free energy of the transferred proton is modified due to the change in the local electrostatic energy, the PMF of moving the group involved and to some extent by the long range interaction with D139 (see discussion in section 6).

As clarified in our previous work it is very hard to get convergent results for the electrostatic energy of the different states by microscopic simulations. This is the reason for our use of the PDLD/S-LRA model. Although our semi-macroscopic calculations are based on many years of validations, the situation in the intirior of CcO is very complex and may involve some traps (in particular the identity of the pump site). Here, we do know the pKapp of E286 but it is harder to obtain unique results from activation barriers. In this respect, focusing on the barrier for PTR from E286 to B in the N139D mutant (where we are confident that we are dealing with this specific barrier) offers a unique opportunity to validate the PDLD/S-LRA (EVB) calculations.

It is important to mention other attempts to rationalize the effect of mutations in the D pathway. Significant attention has been placed on the idea of delocalization of the proton (e.g.) in the D pathway. While delocalization effects are interesting and potentially important, it has been demonstrated that the activation barriers for biological PT are determined by the corresponding electrostatic profiles. In fact, even in bacteriorhodopsin, where the effort in studies of a delocalized proton has been quite successful, there is no correlation between the traps of the proton and the overall barrier for PT in the rate-limiting steps (see). In our case this is probably the situation as long as the transfer through the D pathway is not rate limiting. Interestingly, even in cases of mutant studies (e.g. Y33H) where this transfer is rate limiting, the rate is probably determined by the highest point on the energy profile and not the depth of the minima.

The work of Ref. and the much more sophisticated subsequent work of Ref. suggested that E286 serves as a valve where the rotations of the type considered here play a major role in preventing leaks and in supporting proton pumping. Although this is an appealing idea, it seems to us that as long as the barriers for rotation are small relative to the other barriers it is hard to use this element as a valve, since it is likely that the system will be able to relax to the lowest energy state of E286. In this case it may be useful to simulate the kinetics by the approach used here and to explore the resulting populations.

The considerations of the possible role of low barriers and possible conflicts with microscopic reversibility should address concerns of those who are not familiar with the power of quasi-equilibrium treatments and the ability to assess the kinetics of the system once the relaxation times are relatively clear. More specifically, if we have two states (say proton on Prd and proton on E286) we can estimate the barrier for a process that occurs at about milliseconds ignoring processes that would occur at much longer time scales (e.g. protein unfolding) and safely assume that the forward and backward rates are fully determined by the free energy barriers evaluated by considering all the configurations available in the millisecond time scale. This is related to our considerations of the time dependence of the energetics of proton transfer in bacterial reaction center, where we know the rate of change in the energy of the charge-separated state between picoseconds to nanoseconds and they would be relevant for processes that occurs as long as the back reaction is within the nanosecond range. That is, after the initial relaxation from the P* to The P+H- state (see) if we look for the competition between the forward rate (electron transfer to the quinone) and the back reaction after a nanosecond we can use the energy of the relaxed (P+H- state). Overall perhaps the best current way for handling the relaxation and quasi-equilibrium problem is to use a time dependent dielectric constant as in the approach proposed by us. Using this approach with any reasonable set of dielectric relaxation times it is very hard to construct a model that would support a time dependant gate, which does not involve a formation of another state. In other words, possible proposals of water rearrangement that would block the back reaction are likely to violate quasi-equilibrium conditions.

The above comments are also relevant to the results described in Ref., where it was assumed that the changes of the PT barrier in the D pathway can affect the pumping and explain the effect of the N139 mutation. In the case of the mutations where the reprotonation of E286 is significantly faster than the other steps (see Table 1), the barriers in the D pathway are not relevant (this is even more true if one only considers the energetics of the unprotonated water molecules as was done in ref).

The present work places significant emphasis on quantifying the role of the rotation of E286. Related efforts have been made in Refs. (see the discussion above). Here we would like to point out that while this rotation is important, the rotational barrier obtained by current calculations is too low to provide a gate under any quasi-equilibrium considerations. Stable gates and valves probably require features of the type reproduced in and considered in).

In view of the rather long discussion above, we would like to emphasize that the main unique point of the current work is in using structural information and calculations of the relevant activation barriers to explore the effect of mutations on proton pumping by CcO. Doing so provides interesting insights and allows one to fully exploit the available structural information.


This work was supported by NIH grant GM40283 (to AW) and a grant from the Swedish Research Council (to PB). We gratefully acknowledge the University of Southern California’s High Performance Computing and Communications Center for computer time.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Wikström MKF. Proton pump coupled to cytochrome-c oxidase in mitochondria. Nature. 1977;266(5599):271–273. [PubMed]
2. Michel H, et al. Cytochrome c oxidase: structure and spectroscopy. Annu Rev Biophys Biomol Struct. 1998;27:329–56. [PubMed]
3. Hosler JP, Ferguson-Miller S, Mills DA. Energy transduction: Proton transfer through the respiratory complexes. Ann Rev Biochem. 2006;75:165–187. [PMC free article] [PubMed]
4. Brzezinski P, Larsson G. Redox-driven proton pumping by heme-copper oxidases. Biochim Biophys Acta - Bioenergetics. 2003;1605(1–3):1–13. [PubMed]
5. Brzezinski P, Gennis RB. Cytochrome c oxidase: exciting progress and remaining mysteries. Journal of Bioenergetics and Biomembranes. 2008;40(5):521–531. [PubMed]
6. Ostermeier C, et al. Structure at 2.7 Å resolution of the Paracoccus denitrificans two-subunit cytochrome c oxidase complexed with an antibody FV fragment. Proc Natl Acad Sci USA. 1997;94:10547–10553. [PubMed]
7. Yoshikawa S, et al. Redox-coupled crystal structural changes in bovine heart cytochrome c oxidase. Science. 1998;280:1723–1729. [PubMed]
8. Qin L, et al. Redox-Dependent Conformational Changes in Cytochrome c Oxidase Suggest a Gating Mechanism for Proton Uptake. Biochemistry. 2009;48(23):5121–5130. [PMC free article] [PubMed]
9. Qin L, et al. Identification of conserved lipid/detergent-binding sites in a high-resolution structure of the membrane protein cytochrome c oxidase. Proceedings of the National Academy of Sciences. 2006;103(44):16117–16122. [PubMed]
10. Svensson-Ek M, et al. The X-ray crystal structures of wild-type and EQ(I-286) mutant cytochrome c oxidases from Rhodobacter sphaeroides. J Mol Biol. 2002;321(2):329–339. [PubMed]
11. Luna VM, et al. Crystallographic Studies of Xe and Kr Binding within the Large Internal Cavity of Cytochrome ba3 from Thermus thermophilus: Structural Analysis and Role of Oxygen Transport Channels in the Heme–Cu Oxidases†,‡ Biochemistry. 2008;47(16):4657–4665. [PubMed]
12. Soulimane T, et al. Structure and mechanism of the aberrant ba3-cytochrome c oxidase from Thermus thermophilus. EMBO J. 2000;19(8):1766–1776. [PubMed]
13. Pawate AS, et al. A mutation in subunit I of cytochrome oxidase from Rhodobacter sphaeroides results in an increase in steady-state activity but completely eliminates proton pumping. Biochemistry. 2002;41(45):13417–13423. [PubMed]
14. Pfitzner U, et al. Tracing the D-pathway in reconstituted site-directed mutants of cytochrome c oxidase from Paracoccus denitrificans. Biochemistry. 2000;39:6756–62. [PubMed]
15. Han D, Morgan JE, Gennis RB. G204D, a mutation that blocks the proton-conducting D-channel of the aa(3)-type cytochrome c oxidase from Rhodobacter sphaeroides. Biochemistry. 2005;44(38):12767–12774. [PubMed]
16. Kim YC, Wikstrom M, Hummer G. Kinetic models of redox-coupled proton pumping. Proc Natl Acad Sci USA. 2007;104(7):2169–2174. [PubMed]
17. Siegbahn PEM, Blomberg MRA. Energy diagrams and mechanism for proton pumping in cytochrome c oxidase. Biochim Biophys Acta - Bioenergetics. 2007;1767(9):1143–1156. [PubMed]
18. Olsson MHM, et al. Exploring pathways and barriers for coupled ET/PT in cytochrome c oxidase: A general framework for examining energetics and mechanistic alternatives. Biochim Biophys Acta-Bioenergetics. 2007;1767(3):244–260. [PMC free article] [PubMed]
19. Olkhova E, et al. Dynamic water networks in cytochrome c oxidase from Paracoccus denitrificans investigated by molecular dynamics simulations. Biophys J. 2004;86(4):1873–1889. [PubMed]
20. Siegbahn PEM, Blomberg MRA, Blomberg ML. Theoretical study of the energetics of proton pumping and oxygen reduction in cytochrome oxidase. J Phys Chem B. 2003;107(39):10946–10955.
21. Song YF, Mao JJ, Gunner MR. Electrostatic environment of hemes in proteins: pKas of hydroxyl ligands. Biochemistry. 2006;45(26):7949–7958. [PMC free article] [PubMed]
22. Olsson MHM, Sharma PK, Warshel A. Simulating redox coupled proton transfer in cytochrome c oxidase: Looking for the proton bottleneck. FEBS Lett. 2005;579(10):2026–2034. [PubMed]
23. Olsson MHM, Warshel A. Monte Carlo simulations of proton pumps; On the working principles of the biological valve that controls proton pumping in cytochrome c oxidase. Proc Natl Acad Sci USA. 2006;103(17):6500–6505. [PubMed]
24. Wikstrom M, Verkhovsky MI, Hummer G. Water-gated mechanism of proton translocation by cytochrome c oxidase. Biochim Biophys Acta - Bioenergetics. 2003;1604(2):61–65. [PubMed]
25. Popovic DM, Stuchebrukhov AA. Electrostatic study of the proton pumping mechanism in bovine heart cytochrome c oxidase. J Am Chem Soc. 2004;126(6):1858–1871. [PubMed]
26. Xu JC, Voth GA. Computer simulation of explicit proton translocation in cytochrome c oxidase: The D-pathway. Proc Natl Acad Sci USA. 2005;102(19):6795–6800. [PubMed]
27. Warshel A. Conversion of light energy to electrostatic energy in the proton pump of halobacterium halobium. Photochem Photobiol. 1979;30:285–290. [PubMed]
28. Warshel A, Chu ZT. Nature of the surface crossing process in bacteriorhodopsin: Computer simulations of the quantum dynamics of the primary photochemical event. J Phys Chem B. 2001;105(40):9857–9871.
29. Braun-Sand S, et al. The energetics of the primary proton transfer in bacteriorhodopsin revisited: it is a sequential light-induced charge separation after all. Biochim Biophys Acta. 2008;1777(5):441–52. [PMC free article] [PubMed]
30. Warshel A, Parson WW. Computer Simulations of Electron Transfer Reactions in Solution and Photosynthetic Reaction Centers. Annu Rev Phys Chem. 1991;42:279–309. [PubMed]
31. Warshel A, Schlosser DW. Electrostatic control of the efficiency of light-induced electron transfer across membranes. Proc Natl Acad Sci USA. 1981;78:5564–5568. [PubMed]
32. Ädelroth P, Ek M, Brzezinski P. Factors determining electron-transfer rates in cytochrome c oxidase: investigation of the oxygen reaction in the R. sphaeroides enzyme. Biochimica et Biophysica Acta (BBA) - Bioenergetics. 1998;1367(1–3):107–117. [PubMed]
33. Kaila VR, et al. Glutamic acid 242 is a valve in the proton pump of cytochrome c oxidase. Proc Natl Acad Sci U S A. 2008;105(17):6255–9. [PubMed]
34. Pisliakov AV, et al. Electrostatic basis for the unidirectionality of the primary proton transfer in cytochrome c oxidase. Proc Natl Acad Sci U S A. 2008;105(22):7726–31. [PubMed]
35. Kaila VRI, et al. Mechanism and energetics by which glutamic acid 242 prevents leaks in cytochrome c oxidase. Biochimica Et Biophysica Acta-Bioenergetics. 2009;1787(10):1205–1214. [PubMed]
36. Namslauer A, et al. Redox-coupled proton translocation in biological systems: Proton shuttling in cytochrome c oxidase. Proc Natl Acad Sci USA. 2003;100(26):15543–15547. [PubMed]
37. Gorbikova EA, et al. Time-resolved ATR-FTIR spectroscopy of the oxygen reaction in the D124N mutant of cytochrome C oxidase from Paracoccus denitrificans. Biochemistry. 2007;46(45):13141–13148. [PubMed]
38. Zaslavsky D, et al. Observation of a novel transient ferryl complex with reduced CuB in cytochrome c oxidase. Biochemistry. 1999;38:2307–11. [PubMed]
39. Smirnova IA, et al. Aspartate-132 in Cytochrome c Oxidase from Rhodobacter sphaeroides Is Involved in a Two-Step Proton Transfer during Oxo-Ferryl Formation† Biochemistry. 1999;38(21):6826–6833. [PubMed]
40. Lee HJ, et al. Properties of Arg481 Mutants of the aa3-Type Cytochrome c Oxidase from Rhodobacter sphaeroides Suggest That neither R481 nor the Nearby D-Propionate of Heme a3 Is Likely To Be the Proton Loading Site of the Proton Pump. Biochemistry. 2009;48(30):7123–7131. [PMC free article] [PubMed]
41. Namslauer A, et al. Intramolecular Proton-Transfer Reactions in a Membrane-Bound Proton Pump: The Effect of pH on the Peroxy to Ferryl Transition in Cytochrome c Oxidase†,[perpendicular] Biochemistry. 2003;42(6):1488–1498. [PubMed]
42. Durr KL, et al. A D-pathway mutation decouples the Paracoccus denitrificans cytochrome c oxidase by altering the side-chain orientation of a distant conserved glutamate. J Mol Biol. 2008;384(4):865–77. [PubMed]
43. Brzezinski P, Johansson AL. Variable proton-pumping stoichiometry in structural variants of cytochrome c oxidase. Biochimica et Biophysica Acta (BBA) - Bioenergetics. 2010;1797(6–7):710–723. [PubMed]
44. Warshel A, et al. Modeling electrostatic effects in proteins. Biochim Biophys Acta. 2006;1764(11):1647–1676. [PubMed]
45. Warshel A. Computer Modeling of Chemical Reactions in Enzymes and Solutions. New York: John Wiley & Sons; 1991.
46. Aqvist J, Warshel A. Simulation of Enzyme Reactions Using Valence Bond Force Fields and Other Hybrid Quantum/Classical Approaches. Chem Rev. 1993;93:2523–2544.
47. Koepke J, et al. High resolution crystal structure of Paracoccus denitrificans cytochrome c oxidase: New insights into the active site and the proton transfer pathways. Biochimica et Biophysica Acta (BBA) - Bioenergetics. 2009;1787(6):635–645. [PubMed]
48. Kato M, Pisliakov AV, Warshel A. The barrier for proton transport in aquaporins as a challenge for electrostatic models: The role of protein relaxation in mutational calculations. Proteins: Struct Funct Bioinf. 2006;64(4):829–844. [PubMed]
49. Lee FS, Chu ZT, Warshel A. Microscopic and semimicroscopic calculations of electrostatic energies in proteins by the POLARIS and ENZYMIX programs. J Comp Chem. 1993;14:161–185.
50. Schutz CN, Warshel A. What are the dielectric ‘constants’ of proteins and how to validate electrostatic models. Proteins: Struct Funct Genet. 2001;44:400–417. [PubMed]
51. Sham YY, Muegge I, Warshel A. Simulating proton translocations in proteins: probing proton transfer pathways in the Rhodobacter sphaeroides reaction center. Proteins: Struct Funct Genet. 1999;36(4):484–500. [PubMed]
52. Sham YY, Muegge I, Warshel A. The effect of protein relaxation on charge-charge interactions and dielectric constants of proteins. Biophys J. 1998;74:1744–1753. [PubMed]
53. Blomberg MRA, Siegbahn PEM. A quantum chemical study of the mechanism for proton-coupled electron transfer leading to proton pumping in cytochrome c oxidase. Molecular Physics: An International Journal at the Interface Between Chemistry and Physics. 2010;108(19):2733–2743.
54. Salomonsson L, et al. The timing of proton migration in membrane constituted cytochrome c oxidase. Proc Natl Acad Sci USA. 2005;102(49):17624–17629. [PubMed]
55. Salomonsson L, Brändén G, Brzezinski P. Deuterium isotope effect of proton pumping in cytochrome c oxidase. Biochimica et Biophysica Acta (BBA) - Bioenergetics. 2008;1777(4):343–350. [PubMed]
56. Barbara Schowen K, Schowen RL. [29] Solvent isotope effects on enzyme systems. In: Daniel LP, editor. Methods in Enzymology. Academic Press; 1982. pp. 551–606. [PubMed]
57. Fothergill M, et al. Structure-Energy Analysis of the Role of Metal Ions in Phosphodiester Bond Hydrolysis by DNA Polymerase I. Journal of the American Chemical Society. 1995;117(47):11619–11627.
58. Xu JC, et al. Storage of an excess proton in the hydrogen-bonded network of the D-pathway of cytochrome c oxidase: Identification of a protonated water cluster. Journal of the American Chemical Society. 2007;129(10):2910–2913. [PMC free article] [PubMed]
59. Xu J, Voth GA. Redox-coupled proton pumping in cytochrome c oxidase: Further insights from computer simulation. Biochim Biophys Acta - Bioenergetics. 2008;1777:196–201. [PMC free article] [PubMed]
60. Rousseau R, et al. Modeling protonated water networks in bacteriorhodopsin. Physical Chemistry Chemical Physics. 2004;6(8):1848–1859.
61. Garczarek F, Gerwert K. Functional waters in intraprotein proton transfer monitored by FTIR difference spectroscopy. Nature. 2006;439(7072):109–112. [PubMed]
62. Namslauer I, et al. A pathogenic mutation in cytochrome c oxidase results in impaired proton pumping while retaining O(2)-reduction activity. Biochim Biophys Acta. 2010;1797(5):550–6. [PubMed]
63. Parson WW, Chu ZT, Warshel A. Reorganization Energy of the Initial Electron-Transfer Step in Photosynthetic Bacterial Reaction Centers. Biophys J. 1998;74:182–191. [PubMed]
64. Warshel A, Parson WW. Dynamics of biochemical and biophysical reactions: insight from computer simulations. Quart Rev Biophys. 2001;34:563–670. [PubMed]
65. Warshel A. Correlation between structure and efficiency of light-induced proton pumps. In: Packer L, editor. Methods in Enzymology. Academic Press Inc; London: 1986. pp. 578–587.
66. Henry RM, et al. Functional Hydration and Conformational Gating of Proton Uptake in Cytochrome c Oxidase. Journal of Molecular Biology. 2009;387(5):1165–1185. [PubMed]
67. Lepp H, et al. Impaired proton pumping in cytochrome c oxidase upon structural alteration of the D pathway. Biochimica et Biophysica Acta (BBA) - Bioenergetics. 2008;1777(7–8):897–903. [PubMed]