Modeling the Energy Landscape of Side Reactions in the Cytochrome bc1 Complex

Much of the metabolic molecular machinery responsible for energy transduction processes in living organisms revolves around a series of electron and proton transfer processes. The highly redox active enzymes can, however, also pose a risk of unwanted side reactions leading to reactive oxygen species, which are harmful to cells and are a factor in aging and age-related diseases. Using extensive quantum and classical computational modeling, we here show evidence of a particular superoxide production mechanism through stray reactions between molecular oxygen and a semiquinone reaction intermediate bound in the mitochondrial complex III of the electron transport chain, also known as the cytochrome bc1 complex. Free energy calculations indicate a favorable electron transfer from semiquinone occurring at low rates under normal circumstances. Furthermore, simulations of the product state reveal that superoxide formed at the Qo-site exclusively leaves the bc1 complex at the positive side of the membrane and escapes into the intermembrane space of mitochondria, providing a critical clue in further studies of the harmful effects of mitochondrial superoxide production.

summarizes the molecular dynamics simulations used to analyze the reorganization energy related to electron transfer and superoxide unbinding. Fig. S1 summarizes the set of simulation carried out in the free energy perturbation calculations.  Table S1. The MD simulations used for the analysis in the present investigation. A combination of earlier and new MD simulations of various states of the bc 1 complex is used. a) The model of the membrane embedded bc 1 complex was initially constructed and equilibrated in an earlier study Barragan et al. (2015) focusing on the binding of QH 2 in the Qo-site. b) Building upon the equilibrated structure form a, the dynamics and binding of O 2 in the bc 1 complex in a state with Q •− bound in the Qo-site was modeled in another study Husen and Solov'yov (2016), c) In the present investigation, segments of the simulated trajectories from b with O 2 bound in a pocket near the Qo-site were extracted and used for analysis and for setting up a number of new simulations studying the unbinding of O •− 2 after a putative electron transfer according to model I and II, respectively.   Figure S1. Diagram of the 24 performed FEP simulations: Two models for the electron transfer process, see Fig. 2, were modeled by dividing each of them into a discharging (FEP1) and a charging (FEP2) process. Each subprocess was studied through both a forward and a backward simulation. All of this was repeated for three settings of the simulation time per coupling parameter window.

Geometric comparison of QC snapshots
In order to elucidate the effect of the geometry of the Q o -site and its ligands on possible electron to O 2 , a comparison was made between initial geometries of quantum chemically modeled region in snapshots that were observed to lead to different outcomes in terms of electron transfer (see Fig. S2). No obvious effect of geometric features of the protein or the semiquinone anion can be seen, but electron transfer to O 2 appear to require the O 2 molecule to be confined to a small region near Fe 2 S 2 and the head group of of Q •− .

Extended quantum region
The QC calculations presented in Figs. 4 and 5 and Table 1 were repeated with the quantum cluster model (Fig. 3) extended to include the residues V293, P294 and E295 of cytochrome b. The results of the extended QC model are shown in Figs S3 and S4 and Table S2. With the extended model, the electron transfer Figure S2. Geometric comparison of snapshots depending on outcome in QC calculations. The figure shows overlaid initial geometries of different parts of the cluster model used in the primary quantum calculations colored by the outcome in terms of observed electron transfer to O 2 . a and c show the variation in the semiquinone (Q •− ) and O 2 geometries in caclulations with O 2 modeled as spin up or spin down, respectively, while b and d show the same for the amino acid residues included in the quantum model. Figure S3. Local spin densities resulting from QC calculations on MD snapshots using the extended cluster model. a: The local spin densities obtained from QC calculations on MD snapshots with O 2 bound near the Q o -site. The snapshots where electron transfer to O 2 has occurred, as determined by a change in O 2 spin density, are indicated with orange color, while snapshots with no electron transfer are shown in blue. b: The difference ∆s between the local spin density in each snapshot where electron transfer is observed and the average of all snapshots with no electron transfer. The top half of the figure shows the results of QC calculations with O 2 in the spin up state, while the botton shows results for the spin down state. from semiquinone to O 2 (model I) was only observed in three snapshots. This shows that the method is highly sensitive to the choice of model, although the inclusion of more residues does not change the overall conclusion that pure QC calculations indicate two possibilities for electron transfer. In the main paper, we include the results of the smaller quantum region, as this model has better statistics of the electron transfer events.
As E295 is believed to receive the second proton from QH 2 , either directly or indirectly through Y147 (Crofts et al., 1999;Postila et al., 2013;Barragan et al., 2016), including the residue in quantum chemical modeling can also serve as a test of the stability of the anionic semiquinone in the Q o -site. Out of 118 (109) completed calculations, Q •− was observed to receive a proton in 4 (6) cases, when O 2 was modeled in the spin-up (spin-down) state, i.e. it remained in anionic form in at least 95% of cases. This result indicates that the semiquinone anion exists at least as a metastable state at the Q o -site.  (2), was observed following QC geometry optimization of randomly selected MD simulation snapshots using the quantum cluster model extended to include V293, P294 and E295 of cyt b. The energies are computed relative to the mean energy of snapshots with no electron transfer observed.
The O 2 /O •− 2 binding pocket The QC calculations in the present work are started from snapshots of MD simulations from an earlier study (Husen and Solov'yov, 2016) that identified a number of O 2 binding sites within the bc 1 complex by studying the diffusion of O 2 molecules initially added to the bulk water. Fig. S5a shows an O 2 molecule trapped in a binding pocket near the semiquinone ring and Fe 2 S 2 after migrating from the membrane through the narrow channel occupied by the tail of Q •− (Husen and Solov'yov, 2016)).  (Husen and Solov'yov, 2016), which also traps superoxide for a while after electron transfer. a: A snapshot of O 2 bound in the reactant state along with a transparent rendering of the protein and Q •− in a small subset of other trajectory frames to illustrate the conformational variation (Husen and Solov'yov, 2016). The red surface indicates the localization of O 2 during a binding event. b: Close-up of the same snapshot with distances (inÅ) between O 2 and nearby amino acids indicated. c: A snapshot from the product state simulations (model I) with O •− 2 bound in place of O 2 and Q •− converted to Q. intermembrane space in mitochondria, when it unbinds. The dynamics of O •− 2 on the way out appear to be more complex than for model I, but it generally unbinds a lot faster. Fig. S6b shows a characteristic binding time of 1.6 ns, when a threshold distance of 20Å away from Fe 2 S 2 is used to characterize unbinding. The results for the two sub-transformations (FEP1 and FEP2) in each model have been combined, such that the rightmost value in the diagram indicates the total free energy change, ∆G 0 , for the studied reaction.

Free energy perturbation simulations
In FEP calculations, the O 2 molecule was kept within a small region with a radius of 8Å and 10Å for model I and II, respectively, by means of an artificial confining potential to prevent it from leaving the position near Fe 2 S 2 entirely. For model I, the center of the confined region was chosen as the midpoint between the carbonyl carbon atoms of TYR302 of cytochrome b and HIS135 of the ISP, and for model II, the midpoint between the closest oxygen atom of Q •− and iron atom of Fe 2 S 2 was used. The confining potential is constantly zero within the boundary and takes the shape of a half-harmonic potential on the outside. In other words, the O 2 molecule only feels the confining potential, if it pushes against the boundary.  If the O 2 molecule remains naturally within the bounded region most of the time, the effect of the artificial potential on the measured free energy of the system should be minimal. Fig. S8 shows the distribution of the bounded distance keeping O 2 at the binding site. For most simulations, the O 2 molecule remains naturally within the bounded region, and the constraint only serves to prevent rare events of O 2 leaving the binding site. In this case, we assume the free energy penalty due to the constraint as negligible. In the FEP2 transformation for model I, however, the O 2 molecule spends a significant Figure S8. Distribution of the constrained distance preventing O 2 from leaving the bound position near the Q o -site. In model I (a), the constrained distance is between O 2 and the midpoint between TYR302 of cytochrome b and HIS135 of the ISP, and the upper distance limit is 8Å, indicated by a dashed line. For model II (b), the constrained distance is defined as the distance between O 2 and the midpoint between the Fe2 iron atom (see Fig. 3) of Fe 2 S 2 and the closest oxygen atom of Q •− , and the upper distance limit is 10Å. The colors indicate different durations of the λ-windows used in the simulations. amount of time near edge of the boundary, which could lead to an error in the free energy calculation. Fig. S9 shows the mean and standard deviation of the constraint distance at each λ-window of the FEP2 transformation for model I. During the backward simulations, the O 2 molecule appears to jump back and forth between a local minimum near its initial position and an alternative position near the boundary of the confined region. However, the O 2 molecule primarily pushes against the boundary at the unphysical intermediate stages of the FEP transition and less so at the final fully coupled product state, indicating that the constraint should not significantly disturb the estimated end-to-end free energy difference between the reactant and the product states.

Raw free energy perturbation results
In this section, the results from the FEP calculations are presented in more detail. Figs. S10-S21 shows the free energy results from the forward and backward transformations as well as the distribution of the sampled potential energy difference ∆U and the convergence of the estimated free energy change for each λ-window. The 12 figures cover the two models, the two partial FEP tranformations (FEP1 and FEP2) for each model and three choices of the simulation length per λ-window (see Fig. S1 in the main text). Tables S3 and S4 summarize the results. The sum ∆G forward + ∆G backward indicates the discrepancy between results from the forward and the backward transformations: if perfect micro-reversibility is achieved, the sum should be zero Dixit and Chipot (2001). ∆G BAR is calculated from a combination of the forward and backward data using the Bennet acceptance-ratio (BAR) estimator Bennett (1976). The ParseFEP plugin   Figure S10. Results of the FEP1 calculations for model I with 10 ps/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S11. Results of the FEP2 calculations for model I with 10 ps/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S12. Results of the FEP1 calculations for model I with 1 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S13. Results of the FEP2 calculations for model I with 1 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S14. Results of the FEP1 calculations for model I with 2 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S15. Results of the FEP2 calculations for model I with 2 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S3). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations.  -183.63 -193.68 -191.53 Combined ∆G 0 = ∆G BAR,FEP1 + ∆G BAR,FEP2 171.49 160.58 162.72 Table S4. Results of free energy perturbation calculations for model II. The free energy change during the forward transformation, ∆G forward , the forward/backward discrepancy, ∆G forward + ∆G backward , as well as estimates based on combined data from the forward and backward transformations using the BAR estimator, ∆G BAR , are shown in the table separately for the FEP1 and FEP2 transformations and for three different settings of the simulation length. ∆G 0 is the combined estimate of the reactant-product free energy change. All energies are in kcal/mol. Figure S16. Results of the FEP1 calculations for model II with 10 ps/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S17. Results of the FEP2 calculations for model II with 10 ps/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S18. Results of the FEP1 calculations for model II with 1 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S19. Results of the FEP2 calculations for model II with 1 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S20. Results of the FEP1 calculations for model II with 2 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations. Figure S21. Results of the FEP2 calculations for model II with 2 ns/λ-window. Top: the accumulated free energy change ∆G is shown as a function of λ in the forward (black curve) and backward (red curve) transformations, respectively. The two curves have been connected at λ = 1, such that the value of the red curve at λ = 0 indicates the discrepancy between the results from the forward and backward transformations (see Table S4). Bottom left: the sampled probability distributions, P frwd (∆U ) and P rev (∆U ), of the potential energy change ∆U for going from reactant to product state at a given set of atomic coordinates is shown for each λ window for the forward (black) and backward (red) transformations, respectively. Bottom right: the convergence of the estimated free energy change within each simulated λ-window in the forward, ∆G frwd , and backward, ∆G rev , transformations.