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 $b{c}_1$ 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 Q_o-site exclusively leaves the $b{c}_1$ 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.

1 Introduction

The cytochrome $b{c}_1$ complex is a transmembrane protein complex in the inner mitochondrial membrane of eukaryotes or the plasma membrane of photosynthetically active bacteria. Through its sophisticated reaction cycle, the Q-cycle (Mitchell, 1975; Lhee et al., 2010; Crofts et al., 2017; Barragan et al., 2021) depicted in Figure 1, it serves as a crucial energy transducer in the ATP synthesis pathway. It is, however, also suspected to be a source of harmful superoxide production through stray reactions with molecular oxygen, O₂ (Cape et al., 2007; Dröse and Brandt, 2008; Rottenberg et al., 2009; Brand, 2010; Guillaud et al., 2014).

FIGURE 1

FIGURE 1. The $b{c}_1$ complex and its reaction cycle. During the Q-cycle, two substrate quinol molecules (QH₂) are oxidized to quinone (Q) at the Q_o-site, while one Q is reduced to QH₂ at the Q_i-site of the $b{c}_1$ complex. In this process, two protons are absorbed from the negative side of the membrane, and four are released to the positive side, hence maintaining the transmembrane electrochemical gradient. An oxygen molecule may occasionally bind in a pocket near the Q_o-site (Husen and Solov’yov, 2016), which could lead to superoxide production.

As part of the Q-cycle, an ubiquinol (QH₂) molecule from the membrane is oxidized to ubiquinone (Q) in a bifurcated reaction at one of the two substrate binding sites, the Q_o-site, of the $b{c}_1$ complex. It was previously shown (Husen and Solov’yov, 2016, Husen and Solov’yov, 2017) that O₂ molecules can occasionally enter the protein complex through the membrane and become trapped close to the iron-sulfur cluster (Fe₂S₂) at the Q_o binding site of the $b{c}_1$ complex during a short-lived state of the Q-cycle, in which a bound QH₂ has been partly oxidized to a semiquinone anion (${\text{Q}}^{•-}$) at the Q_o-site (Crofts et al., 2003; Barragan et al., 2016). As this reaction intermediate is a radical, it is conceivable that it could react with a nearby oxygen molecule to produce a potentially harmful superoxide anion, ${\mathrm{O}}_2^{•-}$.

Earlier computational studies (Salo et al., 2017) have indicated two candidate mechanisms of superoxide production in the $b{c}_1$ complex initiated from a state with a radical semiquinone anion bound at the Q_o-site of the complex and an oxygen molecule trapped at the previously identified binding pocket near Fe₂S₂Husen and Solov’yov (2016). In the present investigation, these are studied separately as model I, in which an electron is transferred from semiquinone to the O₂ molecule,

${\text{Q}}^{•-}+{\text{O}}_2\to \text{Q}+{\text{O}}_2^{•-},(1)$

and model II, in which the electron is instead transferred from the iron-sulfur cluster (Fe₂S₂) in the iron-sulfur protein (ISP) subunit of the $b{c}_1$ complex:

${\text{Fe}}_2{\text{S}}_2^{-}+{\text{O}}_2\to {\text{Fe}}_2{\text{S}}_2+{\text{O}}_2^{•-}.(2)$

The two models are illustrated in Figure 2. The mechanism described in model II would appear to be prohibited due to the high midpoint redox potential of Fe₂S₂ and can essentially be ruled out from the experimental evidence (Sun and Trumpower, 2003; Bleier and Dröse, 2013). It has, however, been included in the present investigation for completeness, as the reaction is consistently observed in both earlier quantum chemical calculations (Salo et al., 2017) and new extended calculations presented here. The inclusion of this implausible reaction will serve as a test for the model and to shed light on the importance of the local environment, when modeling electron transfer processes.

FIGURE 2

FIGURE 2. Two charge transfer mechanisms possibly leading to superoxide production. In model I, an electron is transferred from the radical semiquinone anion to the O₂ molecule, oxidizing the semiquinone to a neutral quinone. In model II, an electron is instead transferred from Fe₂S₂⁻ to O₂ to yield a superoxide radical.

Studying complex biomolecular processes computationally generally require a range of computational modeling tools covering different time and length scales (Korol et al., 2020). The present work employs several computational modeling techniques at both the quantum and classical level to characterize the possible reactions leading to superoxide formation at the Q_o-site of the $b{c}_1$ complex. Essentially, this study aims to complete an endeavor to computationally model superoxide production at the Q_o-site of the $b{c}_1$ complex from O₂ binding to ${\mathrm{O}}_2^{•-}$ unbinding, specifically with anionic semiquinone at the Q_o-site (Husen and Solov’yov, 2016, Husen and Solov’yov, 2017; Salo et al., 2017) to gauge the importance of this mechanism as a contributor to metabolic superoxide production. Firstly, previously reported (Salo et al., 2017) quantum chemical (QC) calculations are repeated with significantly extended sampling over an ensemble of molecular dynamics (MD) simulation snapshots from a previous study (Husen and Solov’yov, 2016) with O₂ bound near Fe₂S₂, leading to a much more clear classification of electron transfer reactions at the Q_o-site into the two main reactions, models I and II described above, with superoxide as a product. The same calculations are also repeated with a modified cluster model to gauge the sensitivity to the local protein environment of the QC calculations. Secondly, extended MD simulation of the product state of the modeled reactions are carried out to study the binding time and unbinding of superoxide. Lastly, the free energy landscapes of the two reactions are studied using two separate MD-based approaches in order to gauge the feasibility of the reactions and estimate possible rates.

The protonation state of semiquinone at the Q_o-site is debated (Pietras et al., 2016), and the choice of an anionic semiquinone in our model is based on the assumption that it would be the most reactive with oxygen, and some experimental indications exist of this form of semiquinone (Mulkidjanian, 2005; Cape et al., 2007). The anionic form is also interesting for superoxide production, as a proposed protective complex (Pietras et al., 2016; Bujnowicz et al., 2019) with a hydrogen bond between ${\text{QH}}^{•}$ and the ISP could not be formed in this state, so if both are possible during the Q-cycle, the anionic form would likely be the most relevant for superoxide production. Similarly, the semiquinone at the Q_o-site could alternatively be a result of a semireverse mechanism, where QH₂ has rapidly donated its second electron to heme b_L in the cytochrome b subunit, in accordance with the Q-cycle, but the electron is then later transferred back to the oxidized quinone (Dröse and Brandt, 2008; Sarewicz et al., 2010). The semireverse reaction is in principle also consistent with our model, but it places fewer constraints on the state and position of the Fe₂S₂-cluster, as it can happen with more delay. Separate computational investigations using the same methodology would be warranted to study mechanisms involving the semireverse reaction or a neutral semiquinone at the Q_o-site as possible contributors to superoxide production.

The combined investigations include all levels of modeling of an important biomolecular process: Diffusion of O₂, the quantum chemistry of electron transfer and the free energy budget of the complex protein-ligand enviroment of the hypothetical reaction. Experimental methods can only address the problem much more indirectly, so the details we can get out of modeling are essential to make sure we understand the mechanism and for this, the completeness of the model is essential to allow comparison with experiment.

2 Computational Methods

The computational protocols for the quantum chemical modeling and all-atom MD simulations employed in the present work are described in the following sections.

2.1 Quantum Chemical Calculations

The quantum chemical calculations followed the protocol of an earlier study (Salo et al., 2017), but with a significantly increased number of included snapshots to allow improved statistics following the aim of the present investigation. A total of 200 snapshots from an earlier MD simulation (Husen and Solov’yov, 2016) of the membrane embedded $b{c}_1$ complex from Rhodobacter capsulatus (PDB ID: 1ZRT (Berry et al., 2004)) in presence of molecular oxygen were extracted at random among the parts of the simulated trajectory with O₂ present in the binding pocket near Fe₂S₂ (see Supplementary Figure S5 in the supplementary material (SM). A region, illustrated in Figure 3, including the bound semiquinone and oxygen molecule, the iron-sulfur cluster and a number of nearby amino acids was extracted from each snapshots and used for quantum chemical calculations. Water molecules were also observed in the vicinity of the O₂ binding position in the earlier simulations (Husen and Solov’yov, 2016), so for each extracted snapshot, the four water molecules closest to the midpoint between TYR302 of the cytochrome b subunit and HIS135 of the ISP subunit were also included in the studied quantum region. The unsaturated bonds in the amino acids were capped with hydrogen. A separate set of calculations was also carried out with the quantum region extended to included the residues V293, P294 and E295 of the cythocrome b subunit. These results of the extended calculations are presented in the SM.

FIGURE 3

FIGURE 3. The quantum region used for QC calculations. The region consists of the head group of the semiquinone anion, the iron-sulfur cluster, 16 amino acid residues from the protein, the bound O₂ molecule and four water molecules. The two iron atoms are modeled as having four and five unpaired electrons, respectively, with anti-parallel spin orientation between the two atoms as an initial guess for the QC calculations. Inset: the naming of the iron and sulfur atoms of Fe₂S₂ used in the model.

Quantum chemical geometry optimization calculations were carried out in Gaussian 09 (Frisch et al., 2013) using the CAM-B3LYP density functional theory method (Yanai et al., 2004) with the 6-31G basis set (Rassolov et al., 2001) on the extracted quantum region from each sampled MD snapshot. Backbone carbons of the amino acids as well as the carbon atom connecting the ring and the tail of ${\mathrm{O}}_2^{•-}$ were held fixed during the optimization in order to impose the overall conformation sampled from MD, while still allowing the system to relax according to the quantum model. The quantum region was divided into fragments to allow setting the charge and multiplicity for the different molecules and amino acids separately in the initial wavefunction guess. Furthermore, the iron-sulfur cluster was split into separate fragments for each of the four atoms to allow modeling the anti-ferromagnetic coupling between the two iron atoms (Szilagyi and Winslow, 2006). The multiplicities for the two iron atoms were set to five and six, respectively, with anti-parallel spins between the two. The multiplicity of five for one of the irons was due to the fact that this iron atom has accepted the first electron from the bound QH₂.

The relatively simple 6-31G basis set with no polarization or diffuse functions was chosen despite the inclusion of iron atoms in the model due to the size of the quantum region with a total of 263 atoms, which would have made the calculations unrealistic with a more advanced basis set. The 6-31G basis set has, however, been shown to be able to model the anti-ferromagnetic coupling in the Fe₂S₂-cluster (Szilagyi and Winslow, 2006) and the pKa of its coordinating histidine residues of the ISP (Barragan et al., 2015) correctly, indicating that a reasonable quantum model of the Fe₂S₂-cluster and its surroundings can be achieved using the basis set.

2.2 Molecular Dynamics Simulations

The MD simulation protocol from earlier studies (Barragan et al., 2015; Husen and Solov’yov, 2016) was adopted, such that equilibrated structures and simulated trajectories from those studies were employed for the further modeling and analysis in the present investigation. Briefly, the X-ray crystal structure of the $b{c}_1$ complex from Rhodobacter capsulatus (PDB ID: 1ZRT (Berry et al., 2004)) was embedded in a membrane patch consisting of a mixture of phosphatidylcholine (PC 18:2/18:2), phosphatidylethanolamine (PE 18:2/18:2) and cardiolipin (CL 18:2/18:2/18:2) with a total of 102 CL, 406 PC and 342 PE lipids. The histidines of the protein were all considered δ-protonated, except for H135 and H156 of the ISP, which are coordinating the Fe₂S₂-cluster through the N_δ-atom. These are instead modeled as ϵ-protonated. Two disulfide bonds were introduced based on inspection of the crystal structure: between C144 and C167 of the cytochrome c₁ subunit and between C138 and C155 of the ISP. The protonation state of H156 upon quinol binding is debated (Crofts et al., 1999; Postila et al., 2013), but as we are modeling the semiquinone state, the histidine would be protonated in either case. Similarly, E295 of cytochrome b was considered as protonated as it is expected to have received the second proton from oxidation of QH₂. Finally, ASP252 of cytochrome b was considered protonated as it appears to form a hydrogen bond with PHE248 in the crystal structure. Both protonation states of ASP252 could be relevant in a proposed switching mechanism (Postila et al., 2016) involved in proton transport at the Q_i-site.

The system was solvated in a water box of $197\text{Å}\times 177\text{Å}\times 142\text{Å}$ size, and a concentration of 0.05 mol/L NaCl was added, resulting in a total of 497,562 atoms. In the adopted equilibrated structures (Barragan et al., 2015; Husen and Solov’yov, 2016), the membrane embedded $b{c}_1$ complex exhibits a stable conformation, where the transmembrane helical scaffold is filled by the Q_i-site quinones and lipids from the membrane. No water molecules are observed in the transmembrane region. While interactions between cardiolipin and the Q_i-site of the $b{c}_1$ complex are reported in the literature (Catucci et al., 2012; Postila et al., 2016), no cardiolipin was found near the Q_i-site in the equilibrated structure. This is most likely due to the fact that lipids where randomly placed in the original simulations (Barragan et al., 2015; Husen and Solov’yov, 2016), and the timescale of lipid diffusion in the membrane is too great to observe proper mixing in realistic MD simulations. However, as the present study focuses on the Q_o-site activity, we are mostly concerned with the overall structural stability in the transmembrane region.

In the present study, the $b{c}_1$ complex was modeled using MD simulations in three different redox states: The reactant state with semiquinone at the Q_o-site and O₂ bound and the product state in model I and II, respectively, with a newly formed superoxide bound. The position of the bound O₂ or ${\mathrm{O}}_2^{•-}$ molecule is depicted in Supplementary Figure S5 in the SM. All earlier and new simulations have quinone at the Q_i-site. The simulated trajectories for the reactant state were taken from an earlier study (Husen and Solov’yov, 2016), while the simulations of the product state in the two models were carried out in the present investigation. The overall simulation protocol including both previous and new simulations is shown in Supplementary Table S1.

The reactant state simulations (Husen and Solov’yov, 2016) included a concentration of molecular oxygen added initially to the water phase in the simulation box to study the dynamics O₂ and identify its potential binding sites in the $b{c}_1$ complex. The simulated trajectories were analyzed to identify events of O₂ binding at a particular site near ${\text{Q}}^{•-}$ and Fe₂S₂ at the Q_o-site (see Supplementary Figure S5), and the parts of the trajectories with a bound O₂ molecule were extracted and used for the energy sampling in the present investigation and for setting up quantum chemical calculations and product state simulations. A binding event was identified as starting when an oxygen molecule comes within a distance of 9 Å from the Fe₂S₂-cluster and ending when it is last observed within 13 Å of Fe₂S₂ before reaching a threshold distance of 20 Å. This definition allows O₂ to temporarily fluctuate more than 13 Å away, as long as it returns again rather than leaving entirely.

The product state simulations were set up by choosing a simulation snapshot with a trapped O₂ molecule in the binding pocket near Fe₂S₂ and modifying the atomic partial charges and force field parameters to model a product state according to model I and II, respectively, based on earlier parametrizations (Kaszuba et al., 2013; Husen and Solov’yov, 2016) of the involved molecular constituents in different redox states. In model I, the Q_o-site was modeled with a neutral quinone and a negatively charged iron-sulfur cluster, while in model II, a semiquinone anion and a neutral Fe₂S₂ was assumed, as specified in Supplementary Table S1. In both cases, the bound O₂ was changed to ${\mathrm{O}}_2^{•-}$. For each model, 100 replicate MD simulations with ${\mathrm{O}}_2^{•-}$ initially bound were then carried out for long enough to observe it unbind, following an approach previously used to model O₂ binding in proteins (Husen et al., 2019).

MD simulations were performed using NAMD (Phillips et al., 2005) with the CHARMM 36 (MacKerell et al., 1998, MacKerell et al., 2004; Best et al., 2012) force field, and VMD (Humphrey et al., 1996) was used extensively for system construction, data analysis and visualization. The force field parameters and partial charges for the heme groups, iron-sulfur clusters and bound quinones and semiquinones were taken from an earlier study (Kaszuba et al., 2013; Postila et al., 2013), except for the partial charges of Fe₂S₂⁻ and ${\text{Q}}^{•-}$ which were obtained separately (Husen and Solov’yov, 2016).

2.3 Free Energy Perturbation Simulations

The free energy perturbation (FEP) simulations (Dixit and Chipot, 2001; Woo and Roux, 2005) follow the same protocol as the MD simulations described above, except that constant pressure, rather than constant volume, is employed in the simulations, and electrostatic interactions were modulated using the alchemical transformation method in NAMD. FEP transformations, where electrostatic interactions were gradually turned on or off through a coupling parameter, λ, were performed to measure the free energy change due to going from reactant to product state, when modeling the reactions in Eqs. 1, 2 for model I and II, respectively. The van der Waals and bonded interactions were assumed to be essentially unchanged after electron transfer, so the same force-field parameters were used for the product and reactant states, and the van der Waals and bonded interactions were kept fully coupled during the transformations, which simplifies the method considerably. The FEP transformations were carried out in two parts – first turning off electrostatic interactions in the reactant configuration with O₂ bound at the binding pocket near Fe₂S₂, FEP1, and then turning the interactions back on in the product configuration with ${\mathrm{O}}_2^{•-}$ bound in place of O₂ FEP2:

${\text{Q}}^{•-}+{\text{O}}_2\stackrel{\text{FEP}1}{\underset{\text{backward}}{\overset{\text{forward}}{\rightleftharpoons }}}{\text{Q}}^{\left(0\right)}+{\text{O}}_2\stackrel{\text{FEP}2}{\underset{\text{backward}}{\overset{\text{forward}}{\rightleftharpoons }}}\text{Q}+{\text{O}}_2^{•-}\left(\text{model I}\right),(3)$

${\text{Fe}}_2{\text{S}}_2^{-}+{\text{O}}_2\stackrel{\text{FEP}1}{\underset{\text{backward}}{\overset{\text{forward}}{\rightleftharpoons }}}{\text{Fe}}_2{\text{S}}_2^{\left(0\right)}+{\text{O}}_2\stackrel{\text{FEP}2}{\underset{\text{backward}}{\overset{\text{forward}}{\rightleftharpoons }}}{\text{Fe}}_2{\text{S}}_2+{\text{O}}_2^{•-}\left(\text{model II}\right),(4)$

where the “(0)” superscript indicates that all electrostatic interactions are turned off, equivalent to setting all partial charges of atoms to zero.

For model I, the part of the system undergoing the FEP transformation included ${\text{Q}}^{•-}$ at the Q_o-site and the nearby bound O₂ molecule, while for model II, the FEP region included O₂, the iron-sulfur cluster and its coordinating cysteine and histidine residues as well as SER158 of the ISP (implicitly included with Fe₂S₂ in Eq. 2), as the partial charges on these residues differ between the modeled charge states of the Fe₂S₂ cluster. Each FEP transformation was carried out in both forward and backward mode, i.e. stepping the coupling parameter, λ, from 0 to 1 and then back again to 0 in 16 steps for each transformation. The simulations were carried out with three settings of the simulation time per λ-window, 10 ps, 1 and 2 ns, to test the sensitivity of the results to the simulation length. The resulting set of 24 FEP calculations is illustrated in Supplementary Figure S1: two model reactions, each modeled through two partial transformations (FEP1 and FEP2) and all transformations carried out in both forward and backward mode and for three choices of simulation length. Backward transformations were carried out starting from the atomic coordinates resulting from the corresponding forward transformations, and the FEP2 forward transformations were continued from the end of the FEP1 forward transformations as also illustrated in the figure. The O₂ molecule was restrained to move within a sphere of radius 8 Å and 10 Å for model I and II, respectively, to prevent it from leaving the binding site. For model I, the center of the confined region was defined 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 Fe2 iron atom of Fe₂S₂ (see Figure 3) and the nearest of the two unprotonated oxygens of ${\text{Q}}^{•-}$ was used.

The ParseFEP plugin (Liu et al., 2012) of VMD was used to analyze the results of the FEP simulations and produce free energy curves using the Bennet acceptance ratio (BAR) estimator (Bennett, 1976) based on the combined results from forward and backward simulations.

3 Results and Discussion

The starting point of the present investigation was an earlier computational study (Husen and Solov’yov, 2016) of the dynamics and binding of molecular oxygen inside the $b{c}_1$ complex (PDB ID: 1ZRT (Berry et al., 2004)). From these earlier simulations, parts of simulated trajectories with O₂ bound in a pocket near the Q_o-site were extracted and used for the simulations and analysis presented here. First, we discuss the results from quantum chemical modeling of a region around the Q_o-site with O₂ bound, where we generate statistics of local spin densities and identify possible chemical reactions leading to superoxide production. Next, we utilize MD-based approaches to estimate the free energy due to such electron transfer processes. Finally, we estimate the possible superoxide production rate based on the free energy results.

3.1 Quantum Chemical Modeling of the Q_o-Site

A total of 200 snapshots from MD simulations with O₂ trapped close to the Q_o-site of the $b{c}_1$ complex were randomly selected, and a 263 atom “quantum region” consisting of the local environment at the Q_o-site (see Figure 3) was extracted and analyzed in QC calculations using Gaussian 09 (Frisch et al., 2013; Yanai et al., 2004; Rassolov et al., 2001) following a protocol described previously (Salo et al., 2017). The large number of snapshots was employed to extract ensemble statistics from the earlier MD simulations Husen and Solov’yov (2016) in a local minimum characterized by having O₂ bound in the pocket near Fe₂S₂ depicted in Supplementary Figure S5 in the SM, while ${\text{Q}}^{•-}$ is bound at the Q_o-site. The need for increased sampling is emphasized by the large variability in outcome from quantum chemical modeling observed earlier (Salo et al., 2017). The QC calculations presented here essentially consist of a geometry optimization using the CAM-B3LYP DFT method using the 6-31G basis set (Rassolov et al., 2001; Yanai et al., 2004) carried out for each extracted snapshot with backbone carbon atoms kept at fixed positions and thereby carrying the conformational variation sampled from MD simulations. The relatively limited 6-31G basis set has earlier been shown to model important properties of iron-sulfur clusters correctly (Szilagyi and Winslow, 2006; Barragan et al., 2015) and was chosen to support the inclusion of a rather large region at the Q_o-site in the quantum model. The limitations of the basis set, lacking diffuse and polarization functions, should however be considered as a source of uncertainty in the results as also discussed below.

In its ground state, O₂ is a triplet diradical, and the two unpaired electrons were modeled in both spin up and spin down configurations, where the unpaired electrons of the iron atom Fe2 of the Fe₂S₂-cluster (see Figure 3) are considered as spin up. The calculations successfully converged for 127 and 114 snapshots, when modeled with O₂ spin up and O₂ spin down, respectively. After geometry optimization, a spontaneous electron transfer to the bound O₂ molecule was observed in some, but not all cases: as can be seen in Figure 4A, the spin density of O₂ is either 2 (no charge transfer) or 1 (charge transfer occurred) for the spin up case and −2 and −1, respectively, for the spin down case. The sign of the spin densities in Figure 4A indicate the spin orientation of the fragment, so it can e.g. be observed that the anti-parallel spin between the iron atoms of Fe₂S₂ is maintained after the QC geometry optimization.

FIGURE 4

FIGURE 4. Local spin densities resulting from QC calculations on MD snapshots. (A) The local spin densities obtained from QC calculations on MD snapshots with O₂ bound near the Q_o-site. The snapshots where electron transfer to O₂ has occurred, as determined by a change in O₂ spin density, are indicated with orange color, while snapshots with no electron transfer are shown in blue. (B) The difference $\text{Δ}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₂ in the spin up state, while the botton shows results for the spin down state.

To identify the source of the electron transferred to the O₂ molecule, Figure 4B shows the difference in local spin density between each of the studied snapshots, where electron transfer to O₂ was observed, and the average of all snapshots with no electron transfer, i.e. the quantity:

$\text{Δ}{s}_i=s_i-{\langle s\rangle }_{\text{no ET}},i\in \left\{\text{snapshots with ET}\right\}.(5)$

Here, $s_i$ is the local spin density calculated for a snapshot i, where electron transfer was observed, and ${\langle s\rangle }_{\text{no ET}}$ is the average of the local spin density among all snapshots, where no electron transfer was observed. In the case of O₂ spin up, the electron is consistently transferred from the Fe₂S₂-cluster (primarily atoms Fe2, S1 and S2, see Figure 3), while in the spin down case, the electron comes either fully from the semiquinone anion or again from the Fe₂S₂-cluster, but this time primarily from the two sulfur atoms. These results confirm earlier findings that an electron can be transferred to O₂ either from ${\text{Q}}^{•-}$ (model I) or from Fe₂S₂ (model II) and its coordinating amino acids (Salo et al., 2017). The extended statistics provided here clearly separates the possible reactions into two principal scenarios, corresponding to the two models studied through MD modeling in the following. A comparison of the initial geometries leading to different electronic configurations after QC calculations, shown in Supplementary Figure S2 in the SM, revealed that the position of the bound O₂ molecule is important for electron transfer to happen, but no clear deciding geometric features of ${\text{Q}}^{•-}$ or the protein are immediately apparent.

The distributions of total energies of the modeled Q_o-site of the $b{c}_1$ complex separated into cases with and without observed charge transfer events are shown in Figure 5. The mean energies and standard deviations of the separate distributions are listed in Table 1 and were used to draw the Gaussian fit functions presented in Figure 5. A small number of outlier cases, determined based on an energy threshold, fall well outside the main identified distributions. Various alternative electron transfers, mostly not involving the O₂ molecule, were observed in the outlier cases. In the snapshots leading to electron transfer from Fe₂S₂⁻ to O₂ (model II), the energy is noticeably lower than for snapshots with no electron transfer, while there is virtually no energy difference in case of electron transfer from ${\text{Q}}^{•-}$ (model I). The calculations indicated that both reactions are clearly allowed, when considering the quantum energetics of the Q_o-site alone. In the following, classical MD methods are used to model effects of the local environment in the two model reactions, including binding of ${\mathrm{O}}_2^{•-}$, $b{c}_1$ complex reorganization and the free energy of the combined reaction complex. This is in poor agreement with the experimental evidence, which essentially rules out Fe₂S₂ as an electron donor (Sun and Trumpower, 2003; Bleier and Dröse, 2013). The relatively small quantum region has, however, been isolated in our model, while the effects of the surroundings, i.e. the solvent or, in this case, the protein environment, is likely to play a significant role in electron transfer processes (Moser and Dutton, 1992; Moser et al., 1992; Miyashita et al., 2003).

FIGURE 5

FIGURE 5. Distribution of total energy of the quantum region after QC geometry optimization. The calculations were performed separately with O₂ assumed to be in the (A) spin up and (B) spin down state, respectively. The snapshots have been classified into cases with no electron transfer and electron transfer from either the semiquinone (model I) or the Fe₂S₂-cluster (model II) to O₂. The energies are shown relative to the mean energy of the snapshots with no charge transfer.

TABLE 1

TABLE 1. Energy statistics from the QC geometry optimizations. The table shows the numbers N of snapshots along with mean energies $\langle E\rangle$ and standard deviations σ, where no electron transfer or electron transfer according to model I, Eq. 1, or model II, Eq. 2, was observed following QC geometry optimization of randomly selected MD simulation snapshots. The energies are computed relative to the mean energy of snapshots with no electron transfer observed.

The size of the studied quantum region and the choice of basis set needed to be balanced due the high computational complexity of QC modeling methods. The modeled region was chosen with proximity to the O₂ molecule bound near the Q_o-site as the main criterion and matches the model studied in an earlier investigation (Salo et al., 2017). To assess the significance of the choice of cluster model on the results from QC modeling, all calculations were repeated using the same MD snapshots, but with the cluster model extended to include E295 from cytochrome b, which is likely to be an important interaction partner with QH₂/ ${\text{Q}}^{•-}$, along with its immediate neighbors, V293 and P294 (Lee et al., 2011; Crofts et al., 2017). The results are shown in Supplementary Figures S3, S4 as well as Supplementary Table S2 of the SM. In this case, only three snapshots out of 109 led to the charge transfer described by model I, while the model II reaction remains prevalent in the extended model. The high sensitivity to the choice of included amino acid residues indicates that a limited cluster model is insufficient to reliably and quantitatively model electron transfer processes in a highly complex protein environment. Furthermore, the use of the limited basis set, 6-31G, also introduce uncertainties, but these are likely overshadowed by the uncertainties due to the choice of quantum region. In the present investigation, the results from quantum chemical modeling are ultimately used in an exploratory fashion to identify possible reactions leading to superoxide production: while the initial motivation was the hypothesis that a highly reactive semiquinone anion could react with an oxygen molecule, the frequent observation of electron transfer from Fe₂S₂⁻ led us to include both models in the further free energy calculations.

Aside from the redox states of the potential reaction partners at the Q_o-site, the protonation of the bound semiquinone was also tracked in the extended quantum calculations to test for possible proton transfer back to ${\text{Q}}^{•-}$. The second proton transfer after binding of QH₂ at the Q_o-site is generally believed to go to E295 (Crofts et al., 2017), either directly or through Y147 as an intermediate (Barragan et al., 2016), so in the extended model, which includes E295, it is conceivable that the proton would be transferred back, if the anionic semiquinone at the Q_o-site is unstable. None of the studied snapshots led to proton transfer from E295 or Y147, but in around 4% of cases, proton transfer from Y302 to ${\text{Q}}^{•-}$ was observed. In the majority of cases, ${\text{Q}}^{•-}$ thus remained fully deprotonated, so the existence of a local minimum with anionic semiquinone bound at the Q_o-site remains plausible.

3.2 Unbinding of Superoxide

To study the binding and dynamics of superoxide after an electron transfer event at the Q_o-site, a set of 100 MD simulations starting with ${\mathrm{O}}_2^{•-}$ bound in place of O₂ in the pocket depicted on Supplementary Figure S5 in the SM were carried out for each of the two models following an approach previously used to model O₂ binding and unbinding in proteins (Husen et al., 2019). All simulations were continued until ${\mathrm{O}}_2^{•-}$ unbinding was observed. ${\mathrm{O}}_2^{•-}$ is indeed able to bind in the same pocket as O₂ near the Q_o-site, where it stays for a while before unbinding as illustrated in Figure 6 and Supplementary Figure S5C for model I. While the neutral O₂ enters the binding pocket from the membrane (Husen and Solov’yov, 2016), ${\mathrm{O}}_2^{•-}$ was found to consistently escape directly into the water phase on the positive side of the membrane, i.e. the intermembrane space in case of mitochondria, in all 100 simulations for each model. Two example trajectories of ${\mathrm{O}}_2^{•-}$ unbinding are shown in Figure 6A with red and blue lines, and the red surface is an isosurface of the ${\mathrm{O}}_2^{•-}$ density averaged over all 100 unbinding simulations, i.e. indicating the general localization of the anion during the simulations.

FIGURE 6

FIGURE 6. The dynamics and binding of a newly formed superoxide (model I). ${\text{O}}_{\text{2}}^{•-}$ stays bound for a while in the binding pocket near the Q_o-site, but eventually escapes into the bulk water. (A) superoxide escapes into the intermembrane space. The red regions indicate the averaged localization of O₂ during all 100 simulated unbinding events. Two example trajectories are shown with red and blue colors. (B) a histogram of the number $N_b\left(t\right)$ of simulations that still has ${\text{O}}_{\text{2}}^{•-}$ bound after a specified simulation time along with a fit of an exponential model, Eq. 6, with a characteristic binding time of $\tau =4.5\text{ns}$.

The observed occupancy of ${\mathrm{O}}_2^{•-}$ in the binding position near the Q_o-site of the $b{c}_1$ complex, defined as the fraction of simulations that still has ${\mathrm{O}}_2^{•-}$ bound, is shown in Figure 6B as a function of simulation time for model I. A fit of an exponential model,

$N_b\left(t\right)=N_0{e}^{-\frac{t}{\tau }},(6)$

reveals a characteristic binding time of $\tau =4.5$ ns for model I and $\tau =1.6$ ns for model II (see Supplementary Figure S6 in the SM). The fact that ${\mathrm{O}}_2^{•-}$ remains briefly bound points to the existence of a local minimum in the product state of the electron transfer reaction, which puts the free energy analysis in the following on a stronger footing. For both models, the superoxide anion consistently leaves into the positive side of the membrane (the intermembrane space, in case of mitochondria), when it unbinds. This is unlike the initial binding of the neutral (and non-polar) O₂ molecule, which arrives through the membrane (Husen and Solov’yov, 2016). The finding clearly indicates that mitochondrial superoxide production originating from the $b{c}_1$-complex will primarily emit superoxide into the intermembrane space.

3.3 Reorganization Energy From Classical MD Simulations

In order to construct free energy diagrams of possible electron transfer processes that could lead to superoxide production at the Q_o-site of the $b{c}_1$ complex, a number of MD simulations with O₂ bound in the binding pocket near Fe₂S₂ were carried out to sample the potential energy surfaces in the reactant and product states; in particular, the potential energy difference, $\text{Δ}U$ , for transferring the system from the reactant to the product charge state for a given set of nuclear coordinates is sampled. The sampled energy histograms are then used to calculate the free energy curves with $\text{Δ}U$ serving as the reaction coordinate. The simulations and analysis were carried out for both of the two models of superoxide formation, differing in the source of the electron transferred to O₂. For both models, simulations from an earlier study (Husen and Solov’yov, 2016) were used in the analysis as the reactant state with O₂ bound near the Q_o-site in the presence of ${\text{Q}}^{•-}$, i.e. after the initial electron transfer from QH₂ in the Q-cycle. Simulations of the product states were then set up with O₂ converted to ${\mathrm{O}}_2^{•-}$ and either ${\text{Q}}^{•-}$ oxidized to Q (model I) or Fe₂S₂⁻ oxidized to Fe₂S₂ (model II) at the Q_o-site.

Figure 7 shows histograms of the energy difference $\text{Δ}U$ sampled from MD simulations. Specifically, $\text{Δ}U$ is obtained by calculating the potential energy in a given simulation snapshot, which defines a set of atomic coordinates, ${\overrightarrow{x}}_1,\dots ,{\overrightarrow{x}}_N$ , by separately assigning atomic partial charges and model parameters describing the system in the reactant and product states:

$\text{Δ}U\left({\overrightarrow{x}}_1,\dots {\overrightarrow{x}}_N\right)=U_{\text{product}}\left({\overrightarrow{x}}_1,\dots {\overrightarrow{x}}_N\right)-U_{\text{reactant}}\left({\overrightarrow{x}}_1,\dots {\overrightarrow{x}}_N\right).(7)$

FIGURE 7

FIGURE 7. Distributions of the potential energy shift from reactant to product charge state. The potential energy difference $\text{Δ}U$ for going from the reactant to product charge state for a given set of atomic coordinates was sampled through MD simulations separately in the product and reactant states for the two studied reactions, Eqs 1, 2, where the electron received by O₂ originates from (A) semiquinone (model I) or from (B) the Fe₂S₂-cluster (model II), respectively.

The smooth curves in Figure 7 indicate normal distributions based on the statistical mean values, $\langle \text{Δ}U\rangle$, and standard deviations, σ, from the sampled $\text{Δ}U$ distributions:

$\rho \left(\text{Δ}U\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{\left(\text{Δ}U-\langle \text{Δ}U\rangle \right)}^2}{2{\sigma }^2}}.(8)$

Except for a small skewness for the product state in model II, the sampled $\text{Δ}U$ appears to be normal distributed, indicating that the sampling is carried out within a single local minimum of the free energy landscape.

The potential energy difference is a function of atomic coordinates alone and is therefore suitable as a reaction coordinate describing the mechanical reorganization accompanying the modeled electron transfer process (Miyashita and Go, 2000; Miyashita et al., 2003). By sampling $\text{Δ}U$ fluctuations in local minima around both product and reactant states, a set of free energy curves can be constructed by taking the logarithm of the densities of states (Warshel, 1982; Miyashita et al., 2003) shown in Figure 7:

$G=-k_{B}T\text{ln}\left[\rho \left(\text{Δ}U\right)\right]+\text{constant}.(9)$

Assuming a normal distribution of $\text{Δ}U$, Eq. 8, the free energy curve acquires a parabolic trend:

$G=\frac{k_{B}T}{2{\sigma }^2}{\left(\text{Δ}U-\langle \text{Δ}U\rangle \right)}^2+G_0.(10)$

The resulting free energy curves for the studied electron transfer processes are depicted in Figure 8. The figure shows the free energy calculated from the sampled histograms in Figure 7 by using Eq. 9 with thick lines, while extrapolated parabolas obtained through Eq. 10 are shown with thin lines. The free energy of the product state was shifted to make the two curves intersect at $\text{Δ}U=0$ (Miyashita et al., 2003). $\text{Δ}{G}^0$ shown in the figure is the difference between the minima of the two free energy curves, and the reorganization energy, ${\text{Λ}}_p$ , is the energy required to reorganize the atoms of the local environment from their equilibrium position in the product state, i.e. after electron transfer, to correspond to the equilibrium positions in the reactant state. The reorganization energy ${\text{Λ}}_p$ and the similarly defined ${\text{Λ}}_r$ for moving atoms in the reactant state to the equilibrium positions in the product state were observed to differ (see Table 2), indicating that the medium fluctuations do not globally follow Gaussian statistics (LeBard and Matyushov, 2010; Martin and Matyushov, 2012; Martin et al., 2013; Matyushov, 2013). The product state value, ${\text{Λ}}_p$ , was used, as the two parabolas intersect at the product state equilibrium when ${\text{Λ}}_p=-\text{Δ}{G}^0$ , leading to maximal overlap integral between the reactant and the product harmonic oscillator states (Moser et al., 1992; Moser and Dutton, 1992).

FIGURE 8

FIGURE 8. Free energy diagrams for the two model reactions. The free energy was calculated as a function of the potential energy shift $\text{Δ}U$ sampled in the MD simulations (see Figure 7) in the reactant and product states. (A) results for model I, (B) results for model II. The thick curves indicate the free energies sampled directly from the simulations, while the thinner curves are the extrapolated parabolas, obtained using Eq. 10. The free energy difference $\text{Δ}{G}^0$ between the reactant and product state minima and the reorganization energy ${\text{Λ}}_p$ calculated in the product state are indicated in the figure, both values in kcal/mol.

TABLE 2

TABLE 2. Free energy parameters and estimated rates obtained from simulations. $\text{Δ}{G}_{\text{reorg}}^0$ is the free energy difference between the reactant and product state equilibria estimated from the two-state reorganization energy calculations, and ${\text{Λ}}_p$ and ${\text{Λ}}_r$ are the corresponding reorganization energies estimated in the product and reactant state, respectively (see Figure 8). $\text{Δ}{G}_{\text{FEP}}^0$ is the free energy change estimated through FEP simulations. $k_{\text{et}}$ is the electron transfer rate estimate obtained using $\text{Δ}{G}_{\text{FEP}}^0$, ${\text{Λ}}_p$ and Eq. 13.

The free energy diagrams in Figure 8 and the parameters $\text{Λ}$ and $\text{Δ}{G}^0$ could be used to estimate the rate of electron transfers leading to superoxide production using the Marcus theory (Marcus and Sutin, 1985; Page et al., 1999). However, as the sampled $\text{Δ}U$ distributions are rather far from directly overlapping, quite extensive extrapolation is involved, which leads to a high degree of uncertainty in the estimated values of $\text{Δ}{G}^0$ and ${\text{Λ}}_p$. To alleviate this uncertainty, a different approach of free energy estimation was employed.

3.4 Free Energy Perturbation Calculations

Instead of connecting the product and reactant states in a single step, an alchemical free energy perturbation (Dixit and Chipot, 2001; Miyashita et al., 2003) (FEP) approach was employed as a second attempt of identifying the free energy change due to an electron transfer to O₂. In this approach, the electrostatic interactions involving O₂/${\mathrm{O}}_2^{•-}$ and ${\text{Q}}^{•-}$ /Q for model I and Fe₂S₂⁻/Fe₂S₂ (including the coordinating amino acids and SER158 of the ISP) for model II were gradually modulated by a parameter λ during MD simulations, where $\lambda ={\lambda }_1,{\lambda }_2,\dots {\lambda }_N$ steps between 0 and 1 in a series of simulation windows. First, the electrostatic interactions were turned off in the reactant state. Then, they were turned back on in the product state charge configuration. Only electrostatic interactions were modulated, as no atoms or bonds are introduced or removed in either of the studied reactions, and the van der Waals interactions were assumed to be the same for the two redox states. The total free energy difference between reactant and product state was then calculated by adding up the contributions from each λ-step:

$\text{ΔΔ}{G}_i=-k_{B}T\text{ln}\langle e^{-\frac{\text{Δ}{U}_i}{k_{B}T}}\rangle ,\text{Δ}{G}^0={\displaystyle \sum }_i\text{ΔΔ}{G}_i.(11)$

Here, $\text{ΔΔ}{G}_i$ is the change in free energy when stepping the λ parameter from ${\lambda }_i$ to ${\lambda }_{i+1}$ , and $\text{Δ}{G}^0$ is the total change in free energy between the reactant and the product states. The sampled variable, $\text{Δ}{U}_i$ , is the difference in potential energy due to changing λ from ${\lambda }_i$ to ${\lambda }_{i+1}$ for a given set of atomic coordinates, sampled in a simulation with a force field corresponding to $\lambda ={\lambda }_i$. The FEP transformation is carried out both in forward and backward mode, and the results are combined to produce a single estimate of the free energy change as a function of λ by using the BAR estimator (Bennett, 1976) in the ParseFEP plugin (Liu et al., 2012) of VMD (Humphrey et al., 1996). The results of the individual forward and backward transformations and the BAR estimates are summarized in Supplementary Tables S3, S4 in the SM. The combined results of the free energy calculations are shown in Supplementary Figure S7. The free energy due to the discharging (FEP1) and charging (FEP2) transitions have been connected in the figure, such that the total end-to-end free energy difference in the diagram reflects the total free energy change, $\text{Δ}{G}^0$ , of the electron transfer reaction. The results are shown for three settings of the simulation length per λ-window. The results are further broken down in Supplementary Figures S10–S21 of the SM, where distributions of the sampled $\text{Δ}{U}_i$ for each λ step are shown along with the separate free energy curves for the forward and backward simulations. When 1 ns or 2 ns simulation time per λ step is used, the sampled $\text{Δ}{U}_i$ appear as normal distributed, and the forward and backward results are generally in good correspondence with discrepancies up to 5 kcal/mol (see Supplementary Tables S3, S4). The effects of the artificial restraints keeping the O₂ molecule from leaving its binding pocket are assumed to be minor, as the O₂ remains naturally in the bound position during most of the simulations as further discussed in the SM.

The electron transfer from ${\text{Q}}^{•-}$ (model I) leads to a significant decrease in free energy of around 20 kcal/mol, indicating a favorable reaction, while the reaction in model II has a quite high free energy cost of 160 kcal/mol. The reaction free energies are summarized in Table 2. The free energy change for the reaction in model I is roughly similar to the result from the reorganization energy sampling (Figure 8), while the discrepancy is somewhat bigger for model II. In the latter case, however, both models agree that the free energy cost is so large that the reaction is virtually impossible. The gradual change of potential in the free energy perturbation approach makes sure that the $\text{Δ}U$ histograms for a forward and backward step are always overlapping, leading to a much higher confidence in the resulting free energy difference. However, barriers in the free energy landscape and the reorganization energy cannot be obtained using the free energy perturbation method, since an alchemical and unphysical route is used to sample the energy profile.

3.5 Rate of Superoxide Production

The obtained free energy shifts and reorganization energies can be used to estimate the possible rates of the modeled electron transfer processes. Marcus theory (Marcus and Sutin, 1985; Miyashita et al., 2003) predicts the following form of the electron transfer rate constant:

$k_{\text{et}} = \frac{2\pi }{\hslash }|T_{AB}{|}^2\frac{1}{\sqrt{4\pi \text{Λ}{k}_{B}T}}\text{exp}(-\frac{{\left(\text{Λ}+\text{Δ}{G}^0\right)}^2}{4\lambda k_{B}T}),(12)$

where $\left|T_{AB}\right|$ is the electronic coupling constant between the reactant and product states, $\text{Λ}$ is the reorganization energy, and $\text{Δ}{G}^0$ is the total change in Gibbs free energy between the minima of the two states. $\text{Λ}={\text{Λ}}_p$ and $\text{Δ}{G}^0$ can be taken from the reorganization and FEP calculations, respectively. Meanwhile, Dutton and Moser (Moser and Dutton, 1992; Moser et al., 1992) found the following empirical relationship for electron transfers in biomolecules:

${\text{log}}_{10}{k}_{\text{et}} = 15-0.6R-3.1\frac{{\left(\text{Δ}G+\text{Λ}\right)}^2}{\text{Λ}},(13)$

where R is the edge-to-edge distance between donor and acceptor in Å, energies are assumed in eV, and the electron transfer rate is in inverse seconds. In the reactant state simulations, the distance between donor and acceptor in model I (the head group of ${\text{Q}}^{•-}$ and the O₂ molecule) fluctuates between 2 and 10 Å. Taking $\text{Λ}={\text{Λ}}_p$ and $\text{Δ}{G}^0=\text{Δ}{G}_{\text{FEP}}^0$ and applying Eq. 13 individually to 892 snapshots at 10 ps intervals from an O₂ binding event in a simulated trajectory in the reactant state leads to average electron transfer rates of $6.1\times {10}^5{s}^{-1}$ for model I as listed in Table 2. The results for model II indicates that the electron transfer from Fe₂S₂⁻ to O₂ is virtually impossible.

The estimated electron transfer rates correspond to the state, where semiquinone and O₂ are simultaneously present at the Q_o-site of the $b{c}_1$ complex. The total superoxide production rate, while the $b{c}_1$ complex is in the semiquinone state, taking into account binding and unbinding of oxygen molecules can be estimated as:

$k_{{\text{O}}_2^{•-}|\text{SQ}} = k_{\text{bind}}\frac{k_{\text{et}}}{k_{\text{et}}+k_{\text{unbind}}},(14)$

where $k_{\text{bind}}$ and $k_{\text{unbind}}$ are the rates of binding and unbinding of O₂ at the binding pocket near the Q_o-site. A previous investigation (Husen and Solov’yov, 2016) estimated the binding and unbinding rates under physiological oxygen concentrations as $k_{\text{bind}}=2\times {10}^5{\text{s}}^{-1}$ and $k_{\text{unbind}}={10}^8{\text{s}}^{-1}$ , which leads to a superoxide production rate in the semiquinone state for model I of:

$k_{{\text{O}}_2^{•-}|\text{SQ}} \simeq 1.2\times {10}^3{\text{s}}^{-1}.(15)$

This is a very rough order of magnitude estimate, as e.g. using instead the reorganization energy estimate from the reactant state, $\text{Λ}={\text{Λ}}_r$ , yields a two orders of magnitude lower electron transfer rate, but the results clearly indicate that the reaction described in model I is plausible, while the electron transfer from Fe₂S₂ in model II is prohibited due the local protein environment. An experimental study on rat mitochondria (Quinlan et al., 2011) estimated values of $k_{{\text{O}}_2^{•-}|\text{SQ}}$ ($k_9$ in the referenced study (Quinlan et al., 2011)) ranging from 1 to 40 ${\text{s}}^{-1}$ by fits of different kinetic models to experimental superoxide production data, which taking into account the uncertainties, is qualitatively consistent with our computational findings. It should be noted that the rate estimations shown here assume that the $b{c}_1$ complex is in a state of the Q-cycle that has semiquinone bound at the Q_o-site. This will most likely only be the case a small fraction of the time (one estimate (Crofts et al., 2006) puts this fraction at $x_{\text{SQ}}\simeq 4\times {10}^{-8}$), so in reality the superoxide production rate per $b{c}_1$ complex will be much lower. This is not so surprising, as under normal conditions the mitochondrial ROS production should be low.

4 Conclusion

Quantum chemical modeling confirmed earlier findings (Salo et al., 2017) of two main electron transfer processes leading to superoxide formation at the Q_o-site of the $b{c}_1$ complex, which were studied as models I and II, where an electron is donated by either the semiquinone anion, ${\text{Q}}^{•-}$, at the Q_o-site or the iron-sulfur cluster, respectively. The more extended statistics in the present study clearly separates the studied snapshots into cases corresponding to the two electron transfer models as well as cases, where no electron transfer takes place. The shift in energy between the cases with the electron at the donor and acceptor, respectively, indicates a notable decrease in energy for model II, while the electron transfer from ${\text{Q}}^{•-}$ in model I leads to a slight increase. The data from QC calculations alone thus suggests that both reactions are possible, and the model II reaction is rather favorable, which is in contrast to the expected high redox potential of the Fe₂S₂-cluster. However, including the effects of reorganization of the local protein environment changes this picture dramatically: the free energy contributions due to the local environment were estimated by two separate methods based on classical all-atom MD simulations, leading to values of the free energy difference between reactant and product states of the studied reactions as well as the corresponding reorganization energies. The free energy contribution obtained through modeling of the entire complex showed that the effects of the local environment is critical. Essentially, the electron transfer from Fe₂S₂ (model II) is virtually impossible, as expected, while the electron transfer from ${\text{Q}}^{•-}$ could happen at low rates. The calculated electron transfer rates are at best rough order of magnitude estimates, but are in reasonable correspondance with available experimental values (Quinlan et al., 2011) and provide evidence pointing to semiquinone at the Q_o-site of the $b{c}_1$ complex as a source of mitochondrial ROS production. The computational model could further be used to model other potential superoxide production mechanisms and to study the effects of e.g. mutations on ROS production rates to further connect the computational results with experimental evidence. Additionally, it is found that superoxide produced at the Qo-site will exclusively escape to the intermembrane space of mitochondria.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

PH preformed the simulations, analyzed the data and prepared the figures. IS supervised the project. Both authors contributed in writing the manuscript.

Funding

The authors would like to thank the Independent Research Fund Denmark, the Volkswagen Foundation (Lichtenberg Professorship to IAS), the DFG, German Research Foundation, (GRK1885—Molecular Basis of Sensory Biology and SFB 1372—Magnetoreception and Navigation in Vertebrates) and the Ministry for science and culture of Lower Saxony (Simulations meet experiments on the nanoscale: Opening up the quantum world to artificial intelligence (SMART)). Computational resources for the simulations were provided by the DeiC National HPC Center, SDU, and the CARL Cluster at the Carl-von-Ossietzky University Oldenburg, which is supported by the DFG and the ministry for science and culture of Lower Saxony. The work was supported by the North-German Supercomputing Alliance (HLRN).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.643796/full#supplementary-material

References

Barragan, A. M., Crofts, A. R., Schulten, K., and Solov’yov, I. A. (2015). Identification of Ubiquinol Binding Motifs at the Qo-Site of the Cytochromebc1Complex. J. Phys. Chem. B. 119, 433–447. doi:10.1021/jp510022w

PubMed Abstract | CrossRef Full Text | Google Scholar

Barragan, A. M., Schulten, K., and Solov’yov, I. A. (2016). Mechanism of the Primary Charge Transfer Reaction in the Cytochrome bc1 Complex. J. Phys. Chem. B. 120, 11369–11380. doi:10.1021/acs.jpcb.6b07394

PubMed Abstract | CrossRef Full Text | Google Scholar

Barragan, A. M., Soudackov, A. V., Luthey-Schulten, Z., Hammes-Schiffer, S., Schulten, K., and Solov’yov, I. A. (2021). Theoretical Description of the Primary Proton-Coupled Electron Transfer Reaction in the Cytochrome Bc1 Complex. J. Am. Chem. Soc. 143, 715–723. doi:10.1021/jacs.0c07799

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennett, C. H. (1976). Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comp. Phy. 22, 245–268. doi:10.1016/0021-9991(76)90078-4

CrossRef Full Text | Google Scholar

Berry, E. A., Huang, L.-S., Saechao, L. K., Pon, N. G., Valkova-Valchanova, M., and Daldal, F. (2004). X-Ray Structure of Rhodobacter Capsulatus Cytochrome bc1: Comparison with its Mitochondrial and Chloroplast Counterparts. Photosynth. Res. 81, 251–275. doi:10.1023/b:pres.0000036888.18223.0e

PubMed Abstract | CrossRef Full Text | Google Scholar

Best, R. B., Zhu, X., Shim, J., Lopes, P. E. M., Mittal, J., Feig, M., et al. (2012). Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 8, 3257–3273. doi:10.1021/ct300400x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bleier, L., and Dröse, S. (2013). Superoxide Generation by Complex Iii: from Mechanistic Rationales to Functional Consequences. Biochimi. Biophys. Acta. 1827, 1320–1331. doi:10.1016/j.bbabio.2012.12.002

CrossRef Full Text | Google Scholar

Brand, M. D. (2010). The Sites and Topology of Mitochondrial Superoxide Production. Exp. Gerontol. 45, 466–472. doi:10.1016/j.exger.2010.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bujnowicz, Ł., Borek, A., Kuleta, P., Sarewicz, M., and Osyczka, A. (2019). Suppression of superoxide production by a spin‐spin coupling between semiquinone and the Rieske cluster in cytochrome bc 1. FEBS Lett. 593, 3–12. doi:10.1002/1873-3468.13296

PubMed Abstract | CrossRef Full Text | Google Scholar

Cape, J. L., Bowman, M. K., and Kramer, D. M. (2007). A semiquinone intermediate generated at the Qo site of the cytochrome bc1 complex: Importance for the Q-cycle and superoxide production. Proceed. Nation. Acad. Scie. 104, 7887–7892. doi:10.1073/pnas.0702621104

PubMed Abstract | CrossRef Full Text | Google Scholar

Catucci, L., De Leo, V., Milano, F., Giotta, L., Vitale, R., Agostiano, A., et al. (2012). Oxidoreductase activity of chromatophores and purified cytochrome bc 1 complex from Rhodobacter sphaeroides: a possible role of cardiolipin. J. Bioenerg. Biomembr. 44, 487–493. doi:10.1007/s10863-012-9447-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Crofts, A. R., Barquera, B., Gennis, R. B., Kuras, R., Guergova-Kuras, M., and Berry, E. A. (1999). Mechanism of Ubiquinol Oxidation by thebc1Complex: Different Domains of the Quinol Binding Pocket and Their Role in the Mechanism and Binding of Inhibitors. Biochem. 38, 15807–15826. doi:10.1021/bi990962m

PubMed Abstract | CrossRef Full Text | Google Scholar

Crofts, A. R., Lhee, S., Crofts, S. B., Cheng, J., and Rose, S. (2006). Proton pumping in the bc1 complex: A new gating mechanism that prevents short circuits. Biochimica et Biophysica Acta (BBA) - Bioenergetics 1757, 1019–1034. doi:10.1016/j.bbabio.2006.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Crofts, A. R., Rose, S. W., Burton, R. L., Desai, A. V., Kenis, P. J. A., and Dikanov, S. A. (2017). The Q-Cycle Mechanism of the bc1 Complex: A Biologist's Perspective on Atomistic Studies. J. Phys. Chem. B. 121, 3701–3717. doi:10.1021/acs.jpcb.6b10524

PubMed Abstract | CrossRef Full Text | Google Scholar

Crofts, A. R., Shinkarev, V. P., Kolling, D. R. J., and Hong, S. (2003). The Modified Q-cycle Explains the Apparent Mismatch between the Kinetics of Reduction of Cytochromes c 1 and b H in the bc 1 Complex. J. Bio. Chem. 278, 36191–36201. doi:10.1074/jbc.m305461200

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixit, S. B., and Chipot, C. (2001). Can Absolute Free Energies of Association Be Estimated from Molecular Mechanical Simulations? The Biotin−Streptavidin System Revisited. J. Phys. Chem. A. 105, 9795–9799. doi:10.1021/jp011878v

CrossRef Full Text | Google Scholar

Dröse, S., and Brandt, U. (2008). The Mechanism of Mitochondrial Superoxide Production by the Cytochrome bc1 Complex. J. Biol. Chem. 283, 21649–21654. doi:10.1074/jbc.M803236200

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2013). Gaussian 09 Citation. Wallingford, CT: Gaussian, Inc.

Guillaud, F., Dröse, S., Kowald, A., Brandt, U., and Klipp, E. (2014). Superoxide production by cytochrome bc1 complex: A mathematical model. Biochim. Biophy. Acta. 1837, 1643–1652. doi:10.1016/j.bbabio.2014.05.358

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual molecular dynamics. J. Mol. Grap. 14, 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Husen, P., Nielsen, C., Martino, C. F., and Solov’yov, I. A. (2019). Molecular Oxygen Binding in the Mitochondrial Electron Transfer Flavoprotein. J. Chem. Inf. Model. 59, 4868–4879. doi:10.1021/acs.jcim.9b00702

PubMed Abstract | CrossRef Full Text | Google Scholar

Husen, P., and Solov’yov, I. A. (2017). Mutations at the QoSite of the Cytochromebc1Complex Strongly Affect Oxygen Binding. J. Phys. Chem. B. 121, 3308–3317. doi:10.1021/acs.jpcb.6b08226

PubMed Abstract | CrossRef Full Text | Google Scholar

Husen, P., and Solov’yov, I. A. (2016). Spontaneous Binding of Molecular Oxygen at the Qo-Site of the bc1 Complex Could Stimulate Superoxide Formation. J. Am. Chem. Soc. 138, 12150–12158. doi:10.1021/jacs.6b04849

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaszuba, K., Postila, P. A., Cramariuc, O., Sarewicz, M., Osyczka, A., Vattulainen, I., et al. (2013). Parameterization of the Prosthetic Redox Centers of the Bacterial Cytochrome Bc1 Complex for Atomistic Molecular Dynamics Simulations. Theor. Chem. Acc. 132, 1–13. doi:10.1007/s00214-013-1370-8

CrossRef Full Text | Google Scholar

Korol, V., Husen, P., Sjulstok, E., Nielsen, C., Friis, I., Frederiksen, A., et al. (2020). Introducing Viking: A Novel Online Platform for Multiscale Modeling. ACS Omega 5, 1254–1260. doi:10.1021/acsomega.9b03802

PubMed Abstract | CrossRef Full Text | Google Scholar

LeBard, D. N., and Matyushov, D. V. (2010). Protein-water electrostatics and principles of bioenergetics. Phys. Chem. Chem. Phys. 12, 15335–15348. doi:10.1039/c0cp01004a

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D.-W., El Khoury, Y., Francia, F., Zambelli, B., Ciurli, S., Venturoli, G., et al. (2011). Zinc Inhibition of Bacterial Cytochromebc1Reveals the Role of CytochromebE295 in Proton Release at the QoSite. Biochemistry 50, 4263–4272. doi:10.1021/bi200230e

PubMed Abstract | CrossRef Full Text | Google Scholar

Lhee, S., Kolling, D. R. J., Nair, S. K., Dikanov, S. A., and Crofts, A. R. (2010). Modifications of Protein Environment of the [2Fe-2S] Cluster of the bc1 Complex. J. Biol. Chem. 285, 9233–9248. doi:10.1074/jbc.m109.043505

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, P., Dehez, F., Cai, W., and Chipot, C. (2012). A Toolkit for the Analysis of Free-Energy Perturbation Calculations. J. Chem. Theory Comput. 8, 2606–2616. doi:10.1021/ct300242f

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKerell, A. D., Feig, M., and Brooks, C. L. (2004). Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 25, 1400–1415. doi:10.1002/jcc.20065

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J., et al. (1998). All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins†. J. Phys. Chem. B. 102, 3586–3616. doi:10.1021/jp973084f

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcus, R. A., and Sutin, N. (1985). Electron transfers in chemistry and biology. Biochimica et Biophysica Acta (BBA) - Review. Bioenerg. 811, 265–322. doi:10.1016/0304-4173(85)90014-x

CrossRef Full Text | Google Scholar

Martin, D. R., LeBard, D. N., and Matyushov, D. V. (2013). Coulomb Soup of Bioenergetics: Electron Transfer in a Bacterial bc1 Complex. J. Phys. Chem. Lett. 4, 3602–3606. doi:10.1021/jz401910e

CrossRef Full Text | Google Scholar

Martin, D. R., and Matyushov, D. V. (2012). Non-gaussian Statistics and Nanosecond Dynamics of Electrostatic Fluctuations Affecting Optical Transitions in Proteins. J. Phys. Chem. B. 116, 10294–10300. doi:10.1021/jp305757t

PubMed Abstract | CrossRef Full Text | Google Scholar

Matyushov, D. V. (2013). Protein Electron Transfer: Dynamics and Statistics. J. Chem. Phys. 139, 025102. doi:10.1063/1.4812788

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitchell, P. (1975). The Protonmotive Q Cycle: a General Formulation. FEBS Lett. 59, 137–139. doi:10.1016/0014-5793(75)80359-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyashita, O., and Go, N. (2000). Reorganization Energy of Protein Electron Transfer Reaction: Study with Structural and Frequency Signature. J. Phys. Chem. B. 104, 7516–7521. doi:10.1021/jp000865z

CrossRef Full Text | Google Scholar

Miyashita, O., Okamura, M. Y., and Onuchic, J. N. (2003). Theoretical Understanding of the Interprotein Electron Transfer between Cytochrome c2 and the Photosynthetic Reaction Center. J. Phys. Chem. B. 107, 1230–1241. doi:10.1021/jp026753k

CrossRef Full Text | Google Scholar

Moser, C., and Dutton, P. (1992). Engineering protein structure for electron transfer function in photosynthetic reaction centers. Biochimica et Biophysica Acta (BBA) - Bioenerg. 1101, 171–176. doi:10.1016/s0005-2728(05)80012-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Moser, C. C., Keske, J. M., Warncke, K., Farid, R. S., and Dutton, P. L. (1992). Nature of Biological Electron Transfer. Nature 355, 796–802. doi:10.1038/355796a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mulkidjanian, A. Y. (2005). Ubiquinol oxidation in the cytochrome bc1 complex: Reaction mechanism and prevention of short-circuiting. Biochim. Biophy. Acta. 1709, 5–34. doi:10.1016/j.bbabio.2005.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Page, C. C., Moser, C. C., Chen, X., and Dutton, P. L. (1999). Natural Engineering Principles of Electron Tunnelling in Biological Oxidation-Reduction. Nature 402, 47, 52. doi:10.1038/46972

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 26, 1781–1802. doi:10.1002/jcc.20289

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietras, R., Sarewicz, M., and Osyczka, A. (2016). Distinct Properties of Semiquinone Species Detected at the Ubiquinol Oxidation Q o Site of Cytochrome bc 1 and Their Mechanistic Implications. J. R. Soc. Interface. 13, 20160133. doi:10.1098/rsif.2016.0133

PubMed Abstract | CrossRef Full Text | Google Scholar

Postila, P. A., Kaszuba, K., Kuleta, P., Vattulainen, I., Sarewicz, M., Osyczka, A., et al. (2016). Atomistic Determinants of Co-enzyme Q Reduction at the Q_i-Site of the Cytochrome Bc1 Complex. Sci. Rep. 6, 33607. doi:10.1038/srep33607

PubMed Abstract | CrossRef Full Text | Google Scholar

Postila, P. A., Kaszuba, K., Sarewicz, M., Osyczka, A., Vattulainen, I., and Róg, T. (2013). Key Role of Water in Proton Transfer at the Qo-Site of the Cytochrome bc1 Complex Predicted by Atomistic Molecular Dynamics Simulations. Biochimi. Biophy. Acta. 1827, 761–768. doi:10.1016/j.bbabio.2013.02.005

CrossRef Full Text | Google Scholar

Quinlan, C. L., Gerencser, A. A., Treberg, J. R., and Brand, M. D. (2011). The Mechanism of Superoxide Production by the Antimycin-Inhibited Mitochondrial Q-Cycle. J. Biol. Chem. 286, 31361–31372. doi:10.1074/jbc.m111.267898

PubMed Abstract | CrossRef Full Text | Google Scholar

Rassolov, V. A., Ratner, M. A., Pople, J. A., Redfern, P. C., and Curtiss, L. A. (2001). 6-31g* Basis Set for Third-Row Atoms. J. Comput. Chem. 22, 976–984. doi:10.1002/jcc.1058

CrossRef Full Text | Google Scholar

Rottenberg, H., Covian, R., and Trumpower, B. L. (2009). Membrane Potential Greatly Enhances Superoxide Generation by the Cytochrome bc1 Complex Reconstituted into Phospholipid Vesicles. J. Biol. Chem. 284, 19203–19210. doi:10.1074/jbc.m109.017376

PubMed Abstract | CrossRef Full Text | Google Scholar

Salo, A. B., Husen, P., and Solov’yov, I. A. (2017). Charge Transfer at the Qo-Site of the Cytochrome bc1 Complex Leads to Superoxide Production. J. Phys. Chem. B. 121, 1771–1782. doi:10.1021/acs.jpcb.6b10403

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarewicz, M., Borek, A., Cieluch, E., Świerczek, M., and Osyczka, A. (2010). Discrimination between two possible reaction sequences that create potential risk of generation of deleterious radicals by cytochrome bc1. Biochim. Biophy. Acta (BBA) - Bioenergetics 1797, 1820–1827. doi:10.1016/j.bbabio.2010.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J., and Trumpower, B. L. (2003). Superoxide Anion Generation by the Cytochrome Bc1 Complex. Archiv. Biochem. Biophy. 419, 198–206. doi:10.1016/j.abb.2003.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Szilagyi, R. K., and Winslow, M. A. (2006). On the accuracy of density functional theory for iron-sulfur clusters. J. Comput. Chem. 27, 1385–1397. doi:10.1002/jcc.20449

PubMed Abstract | CrossRef Full Text | Google Scholar

Warshel, A. (1982). Dynamics of Reactions in Polar Solvents. Semiclassical Trajectory Studies of Electron-Transfer and Proton-Transfer Reactions. J. Phys. Chem. 86, 2218–2224. doi:10.1021/j100209a016

CrossRef Full Text | Google Scholar

Woo, H.-J., and Roux, B. (2005). Calculation of absolute protein-ligand binding free energy from computer simulations. Proceedings of the National Academy of Sciences 102, 6825–6830. doi:10.1073/pnas.0409005102

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanai, T., Tew, D. P., and Handy, N. C. (2004). A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chemical Physics Letters 393, 51–57. doi:10.1016/j.cplett.2004.06.011

CrossRef Full Text | Google Scholar