Multiscale simulations reveal key features of the proton-pumping mechanism in cytochrome c oxidase Ruibin Lianga,b,c,d, Jessica M. J. Swansona,b,c,d, Yuxing Penge, Mårten Wikströmf, and Gregory A. Votha,b,c,d,1 a Department of Chemistry, University of Chicago, Chicago, IL 60637; bInstitute for Biophysical Dynamics, University of Chicago, Chicago, IL 60637; cJames Franck Institute, University of Chicago, Chicago, IL 60637; dComputation Institute, University of Chicago, Chicago, IL 60637; eResearch Computing Center, University of Chicago, Chicago, IL 60637; and fHelsinki Bioenergetics Group, Programme for Structural Biology and Biophysics, Institute of Biotechnology, University of Helsinki, FI-00014 Helsinki, Finland
Edited by Richard Eisenberg, University of Rochester, Rochester, NY, and approved May 25, 2016 (received for review February 9, 2016)
Cytochrome c oxidase (CcO) reduces oxygen to water and uses the released free energy to pump protons across the membrane. We have used multiscale reactive molecular dynamics simulations to explicitly characterize (with free-energy profiles and calculated rates) the internal proton transport events that enable proton pumping during first steps of oxidation of the fully reduced enzyme. Our results show that proton transport from amino acid residue E286 to both the pump loading site (PLS) and to the binuclear center (BNC) are thermodynamically driven by electron transfer from heme a to the BNC, but that the former (i.e., pumping) is kinetically favored whereas the latter (i.e., transfer of the chemical proton) is rate-limiting. The calculated rates agree with experimental measurements. The backflow of the pumped proton from the PLS to E286 and from E286 to the inside of the membrane is prevented by large free-energy barriers for the backflow reactions. Proton transport from E286 to the PLS through the hydrophobic cavity and from D132 to E286 through the D-channel are found to be strongly coupled to dynamical hydration changes in the corresponding pathways and, importantly, vice versa.
|
cytochrome c oxidase proton-pumping mechanism dynamics bioenergetics proton transport
|
|
| reactive molecular
C
ytochrome c oxidase (CcO) (Fig. 1) catalyzes the reduction of O2 to H2O and couples the free energy of this exergonic reaction to the pumping of protons across the membrane, creating a transmembrane proton electrochemical gradient that drives, for example, ATP synthesis (1–3). During each reaction cycle eight protons are taken up from the negatively charged inside (N-side) of the membrane and either react with oxygen (referred to as “chemical” protons below) or are pumped to the positively charged outside (P-side) of the membrane (referred to as “pumped” protons below). In aa3-type CcO, the D-channel is responsible for uptake of all four pumped protons and at least one out of four chemical protons (4, 5). Protons on the N-side are taken into D-channel via the amino acid residue D132 at the channel entrance, and then transferred to residue E286 in the middle of the membrane (4, 5). By Grotthuss shuttling through the water molecules in the hydrophobic cavity (HC) above E286, each proton is either transferred to react with oxygen in the binuclear center (BNC), consisting of heme a3 and the CuB complex, or transferred to the pump loading site (PLS) and then further released to the P-side of the membrane (Fig. 1). Despite decades of study, the CcO proton-pumping mechanism is still incompletely understood at the atomistic level. For example, controversy remains regarding how the internal proton transport (PT) events are coupled with electron transfer (ET) (3, 6–16), what features enable pumping while preventing proton backflow from the P-side to the N-side (3, 6, 10, 11, 15, 17, 18), and how hydration influences PT through the HC and D-channel (19–26). Recent computational work proposed a stepwise pumping mechanism in which an excess proton is first transported from E286 to the D-propionate on heme a3 (PRDa3) through a poorly hydrated HC, followed by an increase in the HC hydration (24). However, as we have recently reported (27), the migration of a hydrated excess proton can strongly and dynamically influence its solvation 7420–7425 | PNAS | July 5, 2016 | vol. 113 | no. 27
environment, thus drawing into question this prior view of the role of solvation during PT through the HC. To address these questions, we have carried out extensive multiscale reactive molecular dynamics (MS-RMD) (28–33) free-energy simulations to study the proton-pumping mechanism in CcO. Quantum-mechanical forces from targeted quantum mechanics/molecular mechanics (QM/MM) calculations are bridged, in a multiscale fashion via a variational mathematical framework, into the reactive molecular dynamics (MD) algorithm to describe the dynamics of the system nuclei, thus including chemical bond breaking and making. In this way, we explicitly simulate PT between proton binding sites, including Grotthuss shuttling of the excess proton(s) through residues and intervening water molecules. We focus in this work on the transitions where the PM state of the binuclear site (called PM′ herein, where heme a and CuA are reduced) reacts further by receiving an electron from heme a, because it is representative of how PT and ET are coupled in the native catalytic cycle, and there is extensive experimental information to which we can compare our computational results. In the flow-flash experiments, oxidation of the four-electron reduced CcO begins with O2 binding to the BNC (see refs. 3 and 34 for review). Three electrons and one proton from the BNC are then consumed during O–O bond scission (35, 36). The subsequent ET from heme a to BNC generates the PR state, followed by proton uptake at the BNC from N-side bulk to form the F state. Internal PT from E286 to the PLS occurs with the same time dependence as ET (6), and a proton is pumped to the P-side on the same timescale as the PR → F transition. We simulate the internal PT events (excluding proton uptake and release) after O–O bond scission, including the state just before ET (herein denoted PM′) Significance The long-studied proton-pumping mechanism in cytochrome c oxidase (CcO) continues to be a source of debate. This work provides a comprehensive computational characterization of the internal proton transport dynamics, while explicitly including Grotthuss shuttling, that lead to both pumping and catalysis. Focusing on the first steps of oxidation of the fully reduced enzyme, our results show that the transfer of the pumped and chemical protons are thermodynamically driven by electron transfer, and explain how proton back-leakage is avoided by large back-leak barriers and kinetic gating. This work also explicitly characterizes the coupling of proton transport with hydration changes in the hydrophobic cavity and D-channel, thus advancing our understanding of proton transport in biomolecules in general. Author contributions: R.L., J.M.J.S., Y.P., M.W., and G.A.V. designed research; R.L., J.M.J.S., and Y.P. performed research; R.L., J.M.J.S., Y.P., M.W., and G.A.V. analyzed data; and R.L., J.M.J.S., M.W., and G.A.V. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. 1
To whom correspondence should be addressed. Email:
[email protected].
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1601982113/-/DCSupplemental.
www.pnas.org/cgi/doi/10.1073/pnas.1601982113
Fig. 1. Illustration of the simulation setup for the full CcO from R. sphaeroides in a membrane and surrounded by water. The D- and K-channels, as well as metal centers, key residues, and internal water molecules are depicted.
(see SI Appendix, Table S1, for detailed definition on the states). Our results show that both PT events (E286 to the PLS and E286 to the BNC) are driven by ET from heme a to the BNC. The transfer of the pumped proton is kinetically favored while that of the chemical proton is rate limiting. These results also explain how CcO prevents the decoupling of pumping from the chemical reaction with kinetic gating. Given the computational accuracy and efficiency of the MS-RMD methodology (33), a critical component of our results is the explicit characterization of the coupling between PT and hydration changes in both the HC and D-channel. Finally, we present results that argue against the possibility of E286 being biprotonated during the pumping process. Results and Discussion Transport of the Pumped Proton and Hydration of the HC. We have
simulated PT to the PLS and the BNC focusing on coupled hydration changes. Starting with PT from E286 to the PLS, we calculated two-dimensional potentials of mean force (2D PMFs) in the PM′ state (before ET) and PR state (after ET) during the PM′ → PR → F transition (SI Appendix, Fig. S3 A and B). The collective variables used to define these 2D PMFs are (i) the progress of the excess proton center of excess charge (CEC) through the HC (horizontal axis) and (ii) the degree of hydration of the HC (vertical axis). The 2D PMFs and minimum free-energy pathways traced out from those 2D PMFs (black lines in SI Appendix, Fig. S3 A and B) verify that, as the proton moves from E286 to PRDa3, the HC becomes more hydrated (increasing from 6 to 10 waters in our chosen hydration “box,” which also includes approximately 2 waters outside of the traditionally defined interheme region). When E286 is protonated, the HC favors a low hydration state (approximately four waters in the interheme region) (SI Appendix, Fig. S6A). As the excess proton moves to the water above E286 (SI Appendix, Fig. S6B), the PRDa3 side chain rotates down to interact with the positive charge (CEC). This weakens the interactions between PRDa3 and nearby W172 and R481, allowing more water molecules to enter the HC. The hydration level in HC reaches its maximum just before (PM′), or as (PR), PRDa3 is protonated (SI Appendix, Fig. S6C). The PRDa3 then rotates from downward to upward orientation through the transition state (SI Appendix, Fig. Liang et al.
Transport of the Chemical Proton. Focusing next on the transfer of the chemical proton, we calculated similar 2D PMFs for PT from E286 to the BNC in PM′ and PR states, both with and without the PLS (PRAa3) protonated (SI Appendix, Fig. S4 A–D). The 1D curves traced along the minimum free-energy pathways (Fig. 2B) show that PT from E286 to the CuB-bound hydroxide in the BNC is thermodynamically unfavorable by more than 10 kcal/mol in the PM′ state, but favorable by ∼5 kcal/mol in the PR state. This suggests that ET from heme a to BNC also provides a thermodynamic driving force for the PT from E286 to the BNC (i.e., the formation of the F state) that leads to the chemical reaction. This conclusion is independent of the protonation state of the PLS (Fig. 2B) and the hydration level in the HC (SI Appendix, Fig. S4 A–D). The molecular mechanism of PT from E286 to the BNC in the PM′ and PR states is depicted in SI Appendix, Figs. S8 and S9. Note that, in the PM′ state, the waters between protonated E286 and the BNC (not shown in configurations A and B) are only transiently sampled, which is in line with the water gate mechanism discussed in refs. 20 and 22. First, the side chain of protonated E286 rotates toward the BNC (configurations A, B, and C). The rotation contributes to part of the reaction barrier due to intramolecular strain and disruption of the optimal solvation configuration in the HC. Then E286 deprotonates (x ’ 5–6 Å) and transfers the proton to the water molecule between E286 and BNC to form the transition state (configurations D). This is the most narrow region that the proton passes through, with a single water molecule hydrogen bonded to E286, the backbone carbonyl of G283, and CuBOH. Thus, the free-energy barrier is due to electrostatic repulsion from the BNC (much larger for PM′ than for PR), a frustrated water network, and a confined charge distribution as the proton squeezes through this narrow pass. Finally, the combination of the excess proton and CuB-OH completes the chemical reaction to form the F state. The reduced electrostatic repulsion from BNC after ET makes the chemical PT reaction faster and more favorable in PR state than in PM′ state. Rates of the Pumped and Chemical Proton Transport Events. Analysis of the rates (Table 1; reported as time constants, i.e., the inverse of rate constants) for PT to the PLS and BNC in the PM′ and PR states provides deeper mechanistic insight into the proton-pumping mechanism during the PM′ → PR → F transition. Before ET (in the PM′ state), PT from E286 to the PLS is kinetically possible within the experimentally determined A → PR transition [∼25 to 50 μs (17, 41)], but the rate for reverse PT (PLS to E286) outcompetes that of the forward PT and also the reprotonation of E286 through D-channel (Tables 1 and 2). Thus, PT to the PLS is minimal before the ET. In the same state, it is kinetically prohibitive to transfer the PNAS | July 5, 2016 | vol. 113 | no. 27 | 7421
CHEMISTRY
S6D). Subsequently, the proton moves to the A-propionate on heme a3 (PRAa3) (SI Appendix, Fig. S6E). The 1D free-energy curves in Fig. 2A are plotted from the minimum free-energy pathways traced out from the 2D PMFs. These free-energy curves describe the dominant activated reactive energetics for the water-mediated transport of the pumped proton in the PM′ and PR states, revealing several important findings. First, the free-energy minimum on PRAa3 in the pumping PMF of the PR state suggests that PRAa3 is the major PLS, in agreement with previous work (37, 38). It is interesting to note that, in ba3type CcO, PRAa3 is also suggested to be the PLS (39). However, PRDa3 is still actively involved in the pumping process because it shuttles the proton from E286 to the PLS (breaking its salt bridge to R481 and hydrogen bond with W172), as suggested by Wikström et al. (40) and Yamashita and Voth (13). Second, PT from E286 to the PLS is thermodynamically unfavorable by ∼8 kcal/mol in the PM′ state, but favorable by 3 kcal/mol in the PR state (Fig. 2A). Thus, ET to the BNC provides a thermodynamic driving force for PT to the PLS. Combined with the experimental results showing that ET is not complete without PT to the PLS in the E286Q mutant (35), one might conclude that ET to the BNC and PT to the PLS are coupled, each driving the other to its thermodynamically favored state during PM′ → PR transition.
A
B
C
D
Fig. 2. (A) One-dimensional free-energy profiles for PT from E286 to the PLS for the PM′ (red) and PR (blue) states traced out along the minimum free-energy pathways of the full 2D PMFs in SI Appendix, Fig. S3. (B) One-dimensional free-energy profiles for PT from E286 to the BNC in the PM′ state with and without a proton preloaded in the PLS (purple and red, respectively), and the PR state (black and blue, respectively) along the minimum free-energy pathways of the 2D PMFs in SI Appendix, Fig. S4. PT to the BNC with the PLS preloaded (black curve) is the most physically relevant for the PR → F transition. (C) Full 2D PMF for PT from D132 to E286 through the D-channel in the PR state, with a proton preloaded at the PLS. The 2D PMF is a function of the CEC coordinate through the D-channel and the water hydration in the asparagine gate region. The strongly coupled behavior of PT and hydration in the asparagine gate region along the 1D minimum free-energy path (black line) is clearly evident. (D) One-dimensional free-energy trace for PT through the D-channel for PR (blue) and F (red) states, along the minimum free-energy pathway in C and SI Appendix, Fig. S5. The statistical errors of the 1D free-energy profiles in A, B, and C are in the range of ∼0.1–1 kcal/mol. The statistical errors of the 2D PMF in D are in the range ∼0.1–3 kcal/mol. In all of the plots, the positions of E286, PRDa3, PRAa3, BNC, D132, N139, S200, and S201 are labeled with text boxes.
chemical proton from E286 to the BNC, which would short-circuit pumping if it were to happen before PT from E286 to the PLS (22). In contrast, after ET (in the PR state) proton backflow from the PLS to deprotonated E286 is slower than the reprotonation of E286 through the D-channel (Tables 1 and 2). This allows timely reprotonation of E286 and prevents the proton loaded at PLS from leaking back to E286 and subsequently being consumed at the BNC. Thus, full loading of the PLS is achieved only with ET. The reverse is also true considering the experiments in ref. 35. Therefore, we conclude that PT to PLS and ET to BNC gradually drive each other to completion. Moreover, the ET also facilitates the PT from E286 to the BNC (Table 1) by reducing the reaction barrier in PR state. (Fig. 2B). We note that the PT to PLS is still faster than PT to the BNC even after the ET (Table 1, PR state), which again prevents the abovementioned short-circuiting issue. Our calculated rate for transfer
of the chemical proton in PR state is in quantitative agreement with the experimental PR → F transition rate (6, 17) (Table 1, compare the calculated time constant for E286 → BNC in PR state with protonated PLS to the experimental time constant for PR → F transition), and is slower than the PT to PLS (Table 1) and PT from D132 to E286 in D-channel (Table 2). Therefore, we conclude that the PT from E286 to the CuB-bound hydroxide in the BNC is the rate-limiting step for the PR → F transition. Comparison with Proposed Mechanisms. Our thermodynamic and kinetic results build upon the previously proposed mechanism based on the orientation and connectivity of water chains in the HC (6, 15, 20, 22). In this mechanism, PT from E286 to the PLS and ET from heme a to the BNC are tightly coupled to each other and occur in concerted fashion during the PM′ → PR transition. The short-circuiting (i.e., premature PT to BNC before PT to PLS)
Table 1. Calculated and experimental time constants (inverse of rate constants) of proton transport processes E286 → BNC, μs State PM′ PR Exp (A → PR) Exp (PR → F)
E286 → PLS, μs (2.6 ± 0.2) × 10 2.1 ± 0.3
−1
Deprotonated PLS
Protonated PLS
PLS → E286, μs
(4 ± 3) × 10 (3 ± 2) × 102
(7.7 ± 0.2) × 10 (1.7 ± 0.9) × 102
(1.5 ± 0.4) × 10−5 (1.2 ± 0.6) × 103
6
8
25–50 200
Time constants for transport of the pumped proton (E286 → PLS), the chemical proton (E286 → BNC), and the back-leaked proton (PLS → E286) for the PM′ and PR states. The PT from E286 to the BNC should be compared with the experimental rate of the PR → F transition (see the main text). The rate constant calculations are explained in the SI Appendix.
7422 | www.pnas.org/cgi/doi/10.1073/pnas.1601982113
Liang et al.
State PR F Exp (PR → F)
D132 → E286, μs (5 ± 2) × 10−1 (2 ± 1) × 10−1 200
The calculated rates for PT through the D-channel (from D132 to E286) in the PR and F states are compared to the experimental rate of the PR → F transition (17). The PLS is protonated in the simulation. The rate constant calculations are described in the SI Appendix.
is avoided by the absence of a water chain leading from protonated E286 to BNC in the PM′ state. However, this mechanistic proposal was based on classical nonreactive MD simulations in the absence of a shuttling excess proton. Our results support the conclusion that PT to the BNC is unfavorable in PM′ state, but given our treatment of the explicit proton transport they further explain why PT to the PLS is preferred after ET occurs, a question raised but not fully answered in ref. 15. Moreover, our results describe the diffusive PT processes that influence the dynamic water environment in the HC, which in turn is strongly coupled to the PT. Our results are also in partial agreement with the mechanism proposed by Faxén et al. (17) in which PT from E286 to the PLS occurs during PR → F rather than PM′ → PR transition. This mechanism was further discussed later by Lepp et al. (42), who concluded that the electrometric signal during PM′ → PR transition is caused by the upward movement of K362 side chain rather than PT from E286 to PLS. However, it has been shown that this electrometric signal is large enough to include both the lysine swing and the PT from E286 to PLS (15). Moreover, the electrometric signal was shown to be lost in the E286Q mutant (35), resulting in a product that was actually no longer pure PR, but a PR/PM′ mixture with the latter dominating (see section 5.2 and figure 10 in ref. 3). Based on these experimental findings, we suggest that the coupled PT/ET process during PM′ → PR transition is at present the best explanation (15). Proton Transport Through the D-channel. The D-channel is responsible for transporting two protons from the N-side of the membrane to E286 during the PM′ → PR → F transition, one for pumping and the other for the chemical reaction. Three highly conserved asparagine residues, N139, N121, and N207, reside roughly one-third of the way into the D-channel and form a constricted region (called the asparagine gate) (43). Previous nonreactive classical MD simulations have revealed a gating motion of N139 that controls the hydration state of D-channel (25, 26). However, that work hypothesized that PT through the D-channel follows a sequence in which the proton waits on the proton donor (D132) until a “water wire” is formed, and then rapidly dissociates and transports through the preformed water wire to the proton acceptor (E286). As discussed above and in ref. 27, this historically popular depiction of PT in aqueous systems can be misleading and is inconsistent with the dynamically coupled and cooperative nature of the hydration environment and the excess proton migration. Here, as we did earlier for the PT in the HC, we investigate how the motion of the excess proton is coupled to the change of the hydration level across the asparagine gate by explicitly calculating full 2D PMFs for PT through the D-channel in the PR and F states (Fig. 2C and SI Appendix, Fig. S5 A and B). The minimum freeenergy pathways are identified on the 2D PMFs (black lines), and the corresponding 1D traces can be compared in Fig. 2D. The D-channel PT process starts at D132 (x ’ 35 Å) and ends at E286 (x ’ 0–5 Å). Initially, when D132 is protonated, the space between N139 and N121 is narrow and dehydrated, forming an effective gate for PT (SI Appendix, Fig. S7A). As the excess proton transitions to the water above D132, the N139 side chain rotates and opens a pathway for solvation and PT past the asparagine residues (SI Appendix, Fig. S7B). The transition state is reached Liang et al.
when the excess proton is in the middle of the asparagine gate (SI Appendix, Fig. S7C). Once the excess proton shuttles through, the asparagine gate gradually closes and becomes dehydrated again (SI Appendix, Fig. S7D). The curvy and nonhorizontal nature of the minimum free-energy pathways on the 2D PMFs (SI Appendix, Fig. S5 A and B) indicates that PT and hydration changes are strongly coupled processes. After traversing the gate region, the excess proton proceeds to the serine zone (S200 and S201), where it forms hydrogen bonds with the hydroxyl groups of the pore-lining serine residues (x ’ 14 Å). This downhill free-energy region importantly creates a large proton backflow barrier. Subsequently, the excess proton protonates E286 in the “down” conformation (x ’ 5 Å; Fig. 2D), reaching the global free-energy minimum, and E286 rotates up for either pumping or the chemical reaction (x ’ 1 Å; Fig. 2D). The 1D free-energy traces (Fig. 2D) reveal free-energy barriers for the proton to pass through the asparagine gate in the D-channel. For both states, the calculated rates for forward PT (D132 to E286) are much faster than the overall PR → F transition rate (Table 2), confirming that forward PT through the D-channel is not ratelimiting in the PR → F transition. The D-channel PT rates are also faster than the PT backflow rates from the PLS, as discussed above, allowing for fast reprotonation of E286 before proton backflow from the PLS to deprotonated E286. It is plausible, based on these results, that the decoupling mutants, such as N139T (42) and N139D (44), disable kinetic gating by slowing down PT through the D-channel [enabling backflow from the PLS to E286 to outcompete the reprotonation of E286 as previously suggested (3, 26, 39)], by directly slowing down the pumping rate (bypassing pumping altogether), or by lowering the backflow barriers through the D-channel. In ref. 16, an alternative biprotonated mechanism was proposed, in which the excess proton approaches a protonated (neutral) E286 through the D-channel and facilitates proton pumping through a positively charged, biprotonated E286 transition state. Our results show that the self-consistent charge density functional tight-binding (SCC-DFTB) method used in ref. 16 overestimates the pKa for aspartic acid deprotonation in bulk water and that the biprotonated E286 transition state is highly energetically unfavorable (SI Appendix). These results suggest that the conclusions reached in ref. 16 could, at least in part, be attributed to the SCC-DFTB method overestimating the proton affinity of E286 (see the SI Appendix for further discussion). Conclusions The MS-RMD simulations and quantitative analysis of the explicit PT steps in CcO presented in this work, combined with previous experimental findings (see refs. 3, 15, and 18, and references therein), lead us to the following conclusions regarding the most likely sequence of PT and ET events during the PM′ → PR → F transition (Fig. 3A). First, an electron is transferred from heme a to the BNC during the 25- to 50-μs A → PR transition, likely following the chemistry that occurs at the BNC (35, 36). Either coupled with this ET event (during the PM′ → PR transition) or immediately after it (during the PR → F transition), a proton is transferred from a singly protonated E286 to the PLS (PRAa3). Experimental findings suggest that the former (coupled PT/ET) is more likely (3, 15, 35). This PT induces and is accompanied by an increase in the HC hydration from approximately four to eight waters. Second, E286 is rapidly reprotonated through the D-channel, and the HC relaxes back to the low hydration state. Third, the uncompensated negative charge in the BNC caused by ET triggers transfer of the chemical proton from E286 to the CuB-bound hydroxide in the BNC, forming a water molecule. Fourth, E286 is reprotonated again through D-channel, accompanied by partial ET from CuA to heme a and proton ejection from PLS to the P-side. Based on the above reaction sequence, we combine the PMFs and available experimental data to estimate the free-energy diagram for the PM′ → F transition (Fig. 3B). Based on the experimental finding that ET to the BNC in the absence of PT to the PLS (E286Q mutant) results in only 30% PR formation (ref. 3, figure 10), ET from heme a to BNC is estimated to be endergonic by ∼0.5 kcal/mol. Therefore, starting from the PM′ state (state I, PNAS | July 5, 2016 | vol. 113 | no. 27 | 7423
CHEMISTRY
Table 2. Calculated and experimental time constants (inverse of rate constants) for PT in the D-channel
A I IV hem
ea
1
fast 1 1
E286
+
fast
PLS
PRD
PRAa3
w slo
a3
BNC ow a3 l s me he 2
IV V
4
V VIII
+
+
+ +
3
slow
fast
slow
D132
B
+
0 kcal/mol), the PT from E286 to the PLS coupled with ET from heme a to heme a3 leads to PR state (state II) with energy level of −2.9 + 0.5 = −2.4 kcal/mol. Following this, assuming the pKa of D132 is 4.9 [raised by 1 pH unit relative to that in bulk (45), which is ∼3.9 (46)], proton uptake to D132 from N-side bulk leads to state III (0.5 kcal/mol). (This step was not explicitly simulated in our work.) Then PT from D132 to E286 forms state IV (−5.6 kcal/mol). The proton on E286 is subsequently transferred to the BNC, forming the F state (state V, −10.8 kcal/mol). The subsequent proton uptake of D132 from N-side bulk (not explicitly simulated in our work) leads to state VI (−7.9 kcal/mol), followed by the second reprotonation of E286 that leads to state VII (−14.4 kcal/mol). If proton release from the PLS to the P-side bulk is estimated to have exergonicity of 2.8 kcal/mol [not explicitly simulated in our work, but based on the PLS pKa of 5 after F-state formation and P-side pH of 7 (37)], then completing the A → F transition leads to state VIII with an exergonicity of −17.2 kcal/mol. This reflects the overall driving force for the PM′ → F transition, and is larger than the ∼12 kcal/mol exergonicity estimated from experiments (3, 37). However, a recent computational study by Blomberg (47) indicated that the exergonicity of PM → F is ∼20 kcal/mol, in close agreement with our result. Our results show that ET from heme a to the BNC drives both PT events (E286 to the PLS and E286 to the BNC). Among all of these processes, transfer of the chemical proton from E286 to the BNC is rate-limiting, which is in agreement with previous isotope effect studies (48, 49), and our calculated rate is in quantitative agreement with experiments (6, 17). In this mechanism, the pathways that would decouple pumping from transfer of the chemical proton or damage the directionality of proton flow are avoided by kinetic gating in two ways. First, the fast PT from E286 to the PLS precedes proton transfer to the BNC during the entire PM′ → F transition. Second, the possible proton backflow from the PLS to deprotonated E286 is avoided by fast reprotonation of 7424 | www.pnas.org/cgi/doi/10.1073/pnas.1601982113
+ Fig. 3. (A) Reaction sequence during the PM′ → F transition. The redox groups and protonatable sites are labeled. The red minus sign and various colored plus signs indicate electron and proton occupancy, respectively, in the starting state, while arrows indicate charge migration. The first, second, and third proton to travel through the enzyme are colored black, blue, and orange, respectively. The potential leakage pathways that would decouple pumping from reaction are indicated by dashed gray arrows. (B) Free-energy diagram for the reaction sequence during the PM′ → F transition. The resting states are described in the main text. The transition states between them are labeled with “TS,” and the activation barriers are obtained from the free-energy profiles (1D minimum free-energy traces from the 2D PMFs) in Fig. 2 A, B, and D. The black numbers indicate the energy levels for each state in kcal/mol. The red numbers indicate the time constants for the forward transitions between the two neighboring resting states.
E286 through the D-pathway combined with a large barrier for proton backflow from the PLS to E286 after ET from heme a to BNC (PR state). Although PT in the D-channel is clearly coupled with solvation changes through the asparagine gate, it is never found to be rate-limiting in the PM′ → F transition. Although nuclear quantum effects were not included in our PT simulations, they are very unlikely to change our conclusions because the kinetic isotopic effect for PR → F transition is small (∼2.0) (48, 49), corresponding to ∼0.4 kcal/mol decrease in the chemical PT reaction barrier at room temperature. The MS-RMD simulations presented herein thus provide a comprehensive picture of the proton-pumping mechanism in CcO, the important final enzyme in the respiratory electron transfer chain. Our future work will explore proton pumping in CcO mutants and other redox-coupled proteins. Materials and Methods The full structure of CcO from Rhodobacter sphaeroides [Protein Data Bank ID code 1M56 (50)] was embedded in a dimyristoylphosphatidylcholine (DMPC) lipid bilayer and solvated by water molecules on each side of the membrane. MS-RMD simulations using metadynamics (MTD) (51) were performed to identify PT pathways (52) in both the D-channel and HC. The FitRMD method (32, 33) was used to parameterize the MS-RMD models from QM/MM data for protonatable sites in CcO. The MS-RMD umbrella sampling calculations in the HC and D-channel were carried out by restraining (i) the excess proton CEC position (see the SI Appendix for definition) along the PT pathway defined from the MTD procedure, and (ii) the water density in a predefined box (see ref. 27 for definition). Further details on simulation methods are provided in the SI Appendix. ACKNOWLEDGMENTS. We acknowledge Prof. Peter Brzezinski and Prof. Robert Gennis for helpful discussions and Dr. Vivek Sharma, Prof. Ville Kaila, and Prof. Qiang Cui for providing equilibrated structures from their published work, which were very helpful in determining the starting hydration state of the HC for our 2D PMFs. We also acknowledge Dr. J. G. Nelson
Liang et al.
Environment (XSEDE), which is supported by National Science Foundation Grant OCI-1053575, as well as from the University of Chicago Research Computing Center (RCC), the Texas Advanced Computing Center at the University of Texas at Austin, and the US Department of Defense (DOD) High-Performance Computing Modernization Program.
1. Wikstrom MKF (1977) Proton pump coupled to cytochrome c oxidase in mitochondria. Nature 266(5599):271–273. 2. Ferguson-Miller S, Babcock GT (1996) Heme/copper terminal oxidases. Chem Rev 96(7):2889–2908. 3. Kaila VR, Verkhovsky MI, Wikström M (2010) Proton-coupled electron transfer in cytochrome oxidase. Chem Rev 110(12):7062–7081. 4. Konstantinov AA, Siletsky S, Mitchell D, Kaulen A, Gennis RB (1997) The roles of the two proton input channels in cytochrome c oxidase from Rhodobacter sphaeroides probed by the effects of site-directed mutations on time-resolved electrogenic intraprotein proton transfer. Proc Natl Acad Sci USA 94(17):9085–9090. 5. Wikström M, Jasaitis A, Backgren C, Puustinen A, Verkhovsky MI (2000) The role of the D- and K-pathways of proton transfer in the function of the haem-copper oxidases. Biochim Biophys Acta 1459(2-3):514–520. 6. Belevich I, Verkhovsky MI, Wikström M (2006) Proton-coupled electron transfer drives the proton pump of cytochrome c oxidase. Nature 440(7085):829–832. 7. Kim YC, Wikström M, Hummer G (2007) Kinetic models of redox-coupled proton pumping. Proc Natl Acad Sci USA 104(7):2169–2174. 8. Belevich I, Bloch DA, Belevich N, Wikström M, Verkhovsky MI (2007) Exploring the proton pump mechanism of cytochrome c oxidase in real time. Proc Natl Acad Sci USA 104(8):2685–2690. 9. Kim YC, Wikström M, Hummer G (2009) Kinetic gating of the proton pump in cytochrome c oxidase. Proc Natl Acad Sci USA 106(33):13707–13712. 10. Siegbahn PEM, Blomberg MRA (2010) Quantum chemical studies of proton-coupled electron transfer in metalloenzymes. Chem Rev 110(12):7040–7061. 11. Hammes-Schiffer S, Stuchebrukhov AA (2010) Theory of coupled electron and proton transfer reactions. Chem Rev 110(12):6939–6960. 12. Kim YC, Hummer G (2012) Proton-pumping mechanism of cytochrome c oxidase: A kinetic master-equation approach. Biochim Biophys Acta 1817(4):526–536. 13. Yamashita T, Voth GA (2012) Insights into the mechanism of proton transport in cytochrome c oxidase. J Am Chem Soc 134(2):1147–1152. 14. Lu J, Gunner MR (2014) Characterizing the proton loading site in cytochrome c oxidase. Proc Natl Acad Sci USA 111(34):12414–12419. 15. Wikström M, Sharma V, Kaila VRI, Hosler JP, Hummer G (2015) New perspectives on proton pumping in cellular respiration. Chem Rev 115(5):2196–2221. 16. Goyal P, Yang S, Cui Q (2015) Microscopic basis for kinetic gating in cytochrome c oxidase: Insights from QM/MM analysis. Chem Sci (Camb) 6(1):826–841. 17. Faxén K, Gilderson G, Adelroth P, Brzezinski P (2005) A mechanistic principle for proton pumping by cytochrome c oxidase. Nature 437(7056):286–289. 18. Lepp H, Svahn E, Faxén K, Brzezinski P (2008) Charge transfer in the K proton pathway linked to electron transfer to the catalytic site in cytochrome c oxidase. Biochemistry 47(17):4929–4935. 19. Riistama S, et al. (1997) Bound water in the proton translocation mechanism of the haem-copper oxidases. FEBS Lett 414(2):275–280. 20. Wikström M, Verkhovsky MI, Hummer G (2003) Water-gated mechanism of proton translocation by cytochrome c oxidase. Biochim Biophys Acta 1604(2):61–65. 21. Tuukkanen A, Kaila VRI, Laakkonen L, Hummer G, Wikström M (2007) Dynamics of the glutamic acid 242 side chain in cytochrome c oxidase. Biochim Biophys Acta 1767(9):1102–1106. 22. Sharma V, Enkavi G, Vattulainen I, Róg T, Wikström M (2015) Proton-coupled electron transfer and the role of water molecules in proton pumping by cytochrome c oxidase. Proc Natl Acad Sci USA 112(7):2040–2045. 23. Ghosh N, Prat-Resina X, Gunner MR, Cui Q (2009) Microscopic pKa analysis of Glu286 in cytochrome c oxidase (Rhodobacter sphaeroides): Toward a calibrated molecular model. Biochemistry 48(11):2468–2485. 24. Goyal P, Lu J, Yang S, Gunner MR, Cui Q (2013) Changing hydration level in an internal cavity modulates the proton affinity of a key glutamate in cytochrome c oxidase. Proc Natl Acad Sci USA 110(47):18886–18891. 25. Henry RM, Yu CH, Rodinger T, Pomès R (2009) Functional hydration and conformational gating of proton uptake in cytochrome c oxidase. J Mol Biol 387(5):1165–1185. 26. Henry RM, Caplan D, Fadda E, Pomès R (2011) Molecular basis of proton uptake in single and double mutants of cytochrome c oxidase. J Phys Condens Matter 23(23): 234102. 27. Peng Y, Swanson JMJ, Kang SG, Zhou R, Voth GA (2015) Hydrated excess protons can create their own water wires. J Phys Chem B 119(29):9212–9218. 28. Knight C, Lindberg GE, Voth GA (2012) Multiscale reactive molecular dynamics. J Chem Phys 137(22):A525.
29. Swanson JMJ, et al. (2007) Proton solvation and transport in aqueous and biomolecular systems: Insights from computer simulations. J Phys Chem B 111(17): 4300–4314. 30. Liang R, Li H, Swanson JMJ, Voth GA (2014) Multiscale simulation reveals a multifaceted mechanism of proton permeation through the influenza A M2 proton channel. Proc Natl Acad Sci USA 111(26):9396–9401. 31. Yamashita T, Peng Y, Knight C, Voth GA (2012) Computationally efficient multiconfigurational reactive molecular dynamics. J Chem Theory Comput 8(12): 4863–4875. 32. Nelson JG, Peng Y, Silverstein DW, Swanson JMJ (2014) Multiscale reactive molecular dynamics for absolute pKa predictions and amino acid deprotonation. J Chem Theory Comput 10(7):2729–2737. 33. Lee S, Liang R, Voth GA, Swanson JMJ (2016) Computationally efficient multiscale reactive molecular dynamics to describe amino acid deprotonation in proteins. J Chem Theory Comput 12(2):879–891. 34. Brzezinski P, Gennis RB (2008) Cytochrome c oxidase: Exciting progress and remaining mysteries. J Bioenerg Biomembr 40(5):521–531. 35. Gorbikova EA, Belevich I, Wikström M, Verkhovsky MI (2008) The proton donor for O–O bond scission by cytochrome c oxidase. Proc Natl Acad Sci USA 105(31): 10733–10737. 36. Karpefors M, Adelroth P, Namslauer A, Zhen Y, Brzezinski P (2000) Formation of the “peroxy” intermediate in cytochrome c oxidase is associated with internal proton/ hydrogen transfer. Biochemistry 39(47):14664–14669. 37. Wikström M, Verkhovsky MI (2007) Mechanism and energetics of proton translocation by the respiratory heme-copper oxidases. Biochim Biophys Acta 1767(10): 1200–1214. 38. Lee HJ, Ojemyr L, Vakkasoglu A, Brzezinski P, Gennis RB (2009) 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 48(30):7123–7131. 39. Chang HY, et al. (2012) Exploring the proton pump and exit pathway for pumped protons in cytochrome ba3 from Thermus thermophilus. Proc Natl Acad Sci USA 109(14):5259–5264. 40. Wikström M, et al. (2005) Gating of proton and water transfer in the respiratory enzyme cytochrome c oxidase. Proc Natl Acad Sci USA 102(30):10478–10481. 41. Belevich I, et al. (2010) Initiation of the proton pump of cytochrome c oxidase. Proc Natl Acad Sci USA 107(43):18469–18474. 42. Lepp H, Salomonsson L, Zhu JP, Gennis RB, Brzezinski P (2008) Impaired proton pumping in cytochrome c oxidase upon structural alteration of the D pathway. Biochim Biophys Acta 1777(7-8):897–903. 43. Han D, et al. (2006) Replacing Asn207 by aspartate at the neck of the D channel in the aa3-type cytochrome c oxidase from Rhodobacter sphaeroides results in decoupling the proton pump. Biochemistry 45(47):14064–14074. 44. Namslauer A, Pawate AS, Gennis RB, Brzezinski P (2003) Redox-coupled proton translocation in biological systems: Proton shuttling in cytochrome c oxidase. Proc Natl Acad Sci USA 100(26):15543–15547. 45. Sharma V, Ala-Vannesluoma P, Vattulainen I, Wikström M, Róg T (2015) Role of subunit III and its lipids in the molecular mechanism of cytochrome c oxidase. Biochim Biophys Acta Bioenergetics 1847(8):690–697. 46. Lide DR (2004) CRC Handbook of Chemistry and Physics (CRC Press, Boca Raton, FL). 47. Blomberg MRA (2016) Mechanism of oxygen reduction in cytochrome c oxidase and the role of the active site tyrosine. Biochemistry 55(3):489–500. 48. Salomonsson L, Brändén G, Brzezinski P (2008) Deuterium isotope effect of proton pumping in cytochrome c oxidase. Biochim Biophys Acta 1777(4):343–350. 49. Johansson AL, et al. (2011) Proton-transport mechanisms in cytochrome c oxidase revealed by studies of kinetic isotope effects. Biochim Biophys Acta 1807(9): 1083–1094. 50. Svensson-Ek M, et al. (2002) The X-ray crystal structures of wild-type and EQ(I-286) mutant cytochrome c oxidases from Rhodobacter sphaeroides. J Mol Biol 321(2): 329–339. 51. Laio A, Parrinello M (2002) Escaping free-energy minima. Proc Natl Acad Sci USA 99(20):12562–12566. 52. Zhang Y, Voth GA (2011) A combined metadynamics and umbrella sampling method for the calculation of ion permeation free energy profiles. J Chem Theory Comput 7(7):2277–2283.
Liang et al.
PNAS | July 5, 2016 | vol. 113 | no. 27 | 7425
CHEMISTRY
for the initial simulation setup and helpful discussions in the early stages of this project. This research was supported by National Institutes of Health Grant R01-GM053148 (to G.A.V., J.M.J.S., R.L., and Y.P.) and the Sigrid Jusélius Foundation (M.W.). The researchers used computing facilities provided by the Extreme Science and Engineering Discovery