Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Phys Chem B. Author manuscript; available in PMC 2017 December 8.
Published in final edited form as:
PMCID: PMC5721205

Mechanism of the Primary Charge Transfer Reaction in the Cytochrome bc1 Complex


The bc1 complex is a critical enzyme for the ATP production in photosynthesis and cellular respiration. Its biochemical function relies on the so-called Q-cycle, which is well established and operates via quinol substrates that bind inside the protein complex. Despite decades of research, the quinol-protein interaction, which initiates the Q-cycle, has not yet been completely described. Furthermore, the initial charge transfer reactions of the Q-cycle lack a physical description. The present investigation utilizes classical molecular dynamics simulations in tandem with quantum density functional theory calculations, to provide a complete and consistent quantitative description of the primary events that occur within the bc1 complex upon quinol binding. In particular, the electron and proton transfer reactions that trigger the Q-cycle in the bc1 complex from Rhodobacter capsulatus are studied. The coupled nature of these charge transfer reactions was revealed by obtaining the transition energy path connecting configurations of the Qo-site prior and after the transfers. The analysis of orbitals and partial charge distribution of the different states of the Qo-site has further supported the conclusion. Finally, key structural elements of the bc1 complex that trigger the charge transfer reactions were established, manifesting the importance of the environment in the process, which is furthermore evidenced by free energy calculations.

Graphical abstract

An external file that holds a picture, illustration, etc.
Object name is nihms907621u1.jpg


Cellular respiration and photosynthesis constitute the most fundamental energy conversion processes for sustaining living cells. These two processes rely on a series of energy transport events, where the electron and proton transfer reactions are the primarily tools for energy transfer across the cellular energetic apparatus. One of the key elements involved in bacterial photosynthesis and mitochondrial respiration, which utilizes such electron and proton transfer reactions to effectively create a transmembrane proton gradient, is the bc1 complex.13

The bc1 complex is a catalytic transmembrane protein that, by a series of proton and electron transfer reactions, oxidizes quinol (QH2) cofactors at the so called Qo active site and reduces quinone (Q) cofactors at the Qi active site, in an overall process referred to as the Q-cycle.4 The initial step of the Q-cycle corresponds to the binding of a QH2 molecule to the Qo-site followed by two electron transfer reactions, taking place in a bifurcated manner, towards different prosthetic groups of the bc1 complex subunits. As depicted in Fig. 1, one electron is transferred to the heme c1 of the cytochrome (cyt.) c1 subunit via Fe2S2 cluster of the iron-sulfur protein (ISP) subunit, while a second electron to heme bH via heme bL of the cyt. b subunit.3 During each QH2 oxidation, the bifurcated electron transfer occurs alongside two proton transfers from the QH2 to the positive side of the membrane, in order to create a transmembrane potential that is necessary for ATP synthesis.4

Figure 1
Rhodobacter capsulatus bc1 complex

The mechanism of the proton and electron transfers at the Qo-site of the bc1 complex still remains elusive. However, it is believed that the bifurcation of the electrons is accompanied by the proton transfers, i.e., that electrons and protons are transferred simultaneously from QH2 to the bc1 complex in a coupled fashion.5,6 In fact, such a proton-coupled electron transfer (PCET) reactions, are common in various biological systems,7,8 but it has not been proven for the bc1 complex yet. In this study, we investigate the possibility of coupled charge transfer reactions at the Qo-site the bc1 complex from Rhodobacter capsulatus,9 and demonstrate the feasibility of the PCET process.

Previous investigations of the bc1 complex based on molecular dynamics (MD) simulations and quantum chemistry (QC) calculations have revealed two feasible QH2 binding motifs at the Qo-site,10 giving rise to the possibility of two different charge transfer models. The two models, Model I and Model II, differ in the protonation state of the Fe2S2-bound residue H156 (numbering is consistent with crystal structure of the ISP subunit of the Rhodobacter capsulatus bc1 complex9), which is one of the key elements involved in the QH2 binding and, as it will be evidenced in this study, in the subsequent charge transfer reactions at the Qo-site. The Qo-site for both models is featured in Fig. 2. This figure depicts the differences between the two models; in particular, the varied protonation state of the H156 residue and the impact it has on the charge transfer pathways.11

Figure 2
Primary charge transfer reactions at the Qo-site of the bc1 complex

The primary charge transfer reactions at the Qo-site involve one electron and one proton transfer from QH2 to their bc1 complex acceptor sites12,13 and, in general, these two processes can occur either sequentially or simultaneously. The possible scenarios in the deprotonated-H156 Model I and protonated-H156 Model II are schematically illustrated bellow.

Model I

equation image

Model II

equation image

Since the residue H156 is considered deprotonated in Model I, the proton from QH2 can be transfered directly to H156; however, in Model II, the protonated H156 suggests a proton transfer to a different residue, which in previous investigations11,14 was thought to be an H2O molecule. The redox changes of the QH2 and ISP fragments are depicted in Schemes (1) and (2) as orange and magenta arrows, corresponding to electron and proton transfers respectively. The schemes illustrate sequential charge transfer reactions, leading to intermediate states with a charged semiquinone and ISP. In the case of Model I, only QH2 and ISP change their redox states during the transfers, while in Model II an additional water molecule, which acts as an intermediate proton acceptor, is involved. This water molecule is, therefore, included in the reaction scheme of Model II. The diagonal arrows in Schemes (1) and (2) illustrate the possibility of simultaneous electron and proton transfers, corresponding to the PCET regime, that undergoes through a transition state (TS). Such TS corresponds to a state in which both charges have started to be transferred from QH2 to their corresponding acceptors, but these are not yet considered in a final state.

The present investigation supports the hypothesis of a possible PCET reaction as the primary reaction occurring upon quinol binding at the Qo-site of the bc1 complex: exploration of all possible states portrayed in Schemes (1) and (2) allowed to reveal the coupled character of the charge transfer processes. Through an accurate calculation of the reaction energy profiles, the initial state (reactant), TS, and final state (product) of the bc1 complex Qo-site were established. The analysis of molecular orbitals, charge delocalization and electrostatic properties for the reactant and product states evidenced furthermore the coupled nature of the proton and electron transfers and allowed to determine the driving force that stimulates the charge transfer reactions at the Qo-site in the two different binding scenarios of the quinol substrate.


The primary charge transfer reactions occurring at the Qo-site of the bc1 complex were characterized through an in-depth analysis of MD simulations combined with QC calculations of the two different configurations of the complex, corresponding to Model I and II introduced in Fig. 2. The calculations, performed for Model I and Model II, were divided in three main stages, summarized in Table 1. First, MD simulations of the entire system in the reactant state10 allowed to obtain an equilibrated configuration of the bc1 complex that describes its initial state prior charge transfer reactions and, therefore, is considered optimal for QC analysis. Second, QC calculations were performed for a selected fragment of the bc1 complex, composed of the residues largely involved in the charge transfer reactions and the QH2 head group, as shown in Fig. 2. Finally, by using the atomic charges obtained from the QC calculations, a refined reactant, as well as a product state of the entire system, were simulated dynamically to acquire sufficient conformational statistics for describing the PCET free energy calculations.

Table 1
Summary of performed MD simulations and QC calculations. Calculations are listed according to the redox state of the system. MD simulations were performed for all atoms (including protein, cofactors, water, lipids and ions), while QC calculations were ...

All the MD simulations were performed employing NAMD 2.1115 assuming the CHARMM36 force field with CMAP corrections16 for the proteins. The QC calculations were carried out with the Gaussian 09 package,17 employing the UB3LYP DFT method,18 widely used previously in iron-sulfur containing systems optimizations.1925 All images of the bc1 complex, including molecular orbitals and electrostatic potentials, were obtained with VMD Technical details of the methods employed in the calculations are described below.

MD simulations prior charge transfer reactions

In an earlier study,10 360 ns long MD simulations were performed for the bc1 complex represented through Model I and Model II, allowing the bc1 complex, with a bound QH2 at the Qo-site, to equilibrate first and then reach a stable conformation prior to the PCET reaction.

For modeling of the system in VMD 1.9.2,26 X-ray crystal structure of the bc1 complex of Rhodobacter capsulatus (PDB ID: 1ZRT)9 was embedded in a bilayer membrane, composed of 102 cardiolipin (CL 18:2/18:2/18:2/18:2), 406 phosphatidylcholine (PC 18:2/18:2), and 342 phosphatidylethanolamine (PE 18:2/18:2) lipids to represent a mitochondrial membrane.14 The lipid membrane with the embedded protein were solvated within a TIP3P water box at a salt (NaCl) concentration of 0.05 mol/L, and neutralized with salt ions. The Rhodobacter capsulatus crystal structure originally contained stigmatellin and antimycin molecules bound at the Qo- and Qi-site respectively, while the substrate molecules for the Q-cycle are QH2 and Q. The Q and QH2 molecules were thus aligned to the original antimycin and stigmatellin positions. The total simulation system consisted of 500,791 atoms in Model I and 502,165 atoms in Model II, including proteins with cofactors, substrate molecules, lipids, water molecules, and ions.

Addition of the hydrogen atoms that are missing in the crystal structure of the bc1 complex was performed with the VMD plugin psfgen.26 Standard charges and topologies of the bc1 complex proteins were assumed, in accordance with the CHARMM36 force field. However, parameters for the prosthetic groups, hemes and Fe2S2 cluster were adopted to be consistent with earlier investigations,14,27 in which the groups are considered prior PCET, i.e., in the oxidized form. Charges and topology of the QH2 and Q cofactors were taken from an earlier study.27 The standard CHARMM36 force field was employed28 for the PE and PC lipids, as well as for lipid tails of CL. The CL head group charges and parameters were adopted from an earlier investigation.29

All histidine residues of the bc1 complex were considered as δ-protonated except for H156, which has been assumed deprotonated in Model I and ε-protonated in Model II. Inspection of the bc1 complex crystal structure suggested disulfide bonds between the C144 and C167 residues from cyt. c1, and between C138 and C155 residues from ISP. Both disulfide bonds were included in both computational models.

The simulations were performed in the NVT ensemble, where the temperature was kept at 310 K. The long-range electrostatic interactions were calculated, by employing periodic boundary conditions using the PME method,28 with a smooth cutoff of 12 Å; the same cutoff was used for van der Waals interactions. All MD simulations were performed with a time step of 2 fs, and following a simulation protocol in which the system was equilibrated while keeping constrains on selected atoms: (i) first all protein backbone, (ii) then highly movable non-transmembrane segments of the ISP and cyt. c2 subunits, and (iii) finally releasing all the atoms.

The trajectories obtained in an earlier investigation10 were utilized to analyze the reactant state of the bc1 complex and determine the driving force that the environment exerts on the transferring charges at the Qo-site. For this purpose, the analysis of the electrostatic potential was performed on a smoothed electrostatic potential grid, calculated by using the PMEPOT plugin30 in VMD, using the entire 360 ns long MD trajectories for Model I and Model II.

QC calculations

The QC calculations included the residues and cofactors involved in the charge transfer reactions, selected from the model systems previously equilibrated during the 360 ns MD simulations. The set of residues, shown in Fig. 2, constitute the computational models of the Qo-site and were selected based on criteria such as proximity to the quinol head group, proximity to the Fe2S2 cluster, and the residue charge. These selection criteria guaranteed the inclusion of the charge donor and charge acceptor residues, as well as all charged or polar residues that would largely contribute to the driving forces during the charge transfer process.

The computational model of the Qo-site consisted of 168 atoms for Model I and 172 atoms for Model II, including the QH2 head group, Fe2S2 cluster, and pre-equilibrated side chains of cyt. b residues Y147, I292, E295, and Y302, and of ISP residues C133, H135, C138, C153, C155, and H156. In Model II, an H2O molecule was additionally included in the QC calculations as it corresponded to the most probable intermediate proton acceptor10 in this particular case. The Cα atoms of the side chain residues were replaced by CH3 groups, employing for this purpose the MOLEFACTURE plugin of VMD.26

All the QC calculations employed two standard 6-31G(d) and 6-311G(d) basis sets to expand the electronic wave functions. The simpler 6-31G(d) basis set was used to efficiently find optimized reactant and product states of the Qo-site and then a more refined calculations were performed with the 6-311G(d) basis set as it describes more accurately systems with heavy atoms, such as Fe and S, due to its triple-zeta accuracy and additional diffuse functions.31 The comparison of results obtained with the two basis sets is provided in the Supporting Information (SI). Initial QC geometry optimizations of the Qo-site model in the reactant state were performed, and in order to avoid an unphysical collapse of the atoms, the Cα atom positions, taken form the pre-equilibrated structure, were kept fixed during the optimization calculations.

In order to consider the possible sequential and concerted pathways of the electron and proton transfers during the primary charge transfer reactions at the Qo-site, see Scheme (1) and (2), the system was set up in all possible redox states that it could populate during single or coupled charge transfers from QH2 to their respective acceptors (see Table 1). QC geometry optimizations were performed for all the single-transfer states. The transition states (TSs) were obtained through the Synchronous Transit-Guided Quasi-Newton (STQN) Method32,33 in Gaussian 09.17 Additional confirmation of the TS for the PCET was made throughout vibrational frequency calculations followed by intrinsic reaction coordinate (IRC) integration,34,35 also performed in Gaussian 09.

Once the TS, for the two models were established, QC geometry optimizations were carried out from these states towards the reactant and product configurations of the Qo-site, generating a reaction path (energy profile) for the PCET reactions. 31 configurations of the Qo-site were selected from these paths for each model, including the TS, product and reactant states for further analysis, such as the calculations of molecular orbitals and atomic charges, derived following the ESP Merz-Singh-Kollman scheme.36,37

The ESP fitted charges were used for re-defining the atomic charges of Qo-site residues, allowing to perform consequent MD simulations for the bc1 complex in the reactant and product states. For this purpose, Qo-site residues were redefined, allowing to update the topology files needed for MD simulations with refined charges, corresponding to the Qo-site in the reactant and product states. During the charge fitting process, charge symmetries, were taken into account and only the atomic charges of the residue side-chains were reassigned, while the charges of polypeptide backbone atoms were kept at the standard CHARMM36 force field values.38 The modified atomic charges, used in the MD simulations, are provided in the SI.

To stress the influence of Y147 and E295 on the PCET reaction, additional QC calculations of the system were performed for two alternate configurations where (i) Y147 was replaced by a H2O molecule and (ii) where both residues Y147 and E295 were removed.

MD simulations after PCET

MD simulations were performed for the bc1 complex in the reactant and product states, utilizing the refined topology files obtained by reassigning the atomic charges of the Qo-site residues with the ESP charges taken from the QC calculations (see Table 1). By employing the modified charges 80 ns long MD simulations, in the NVT ensemble, were performed using NAMD 2.11, for the bc1 complex Model I. All the simulation parameters, such as temperature, interaction cutoff distances, ions concentration, and others, were the same as in the initial 360 ns long MD simulations. The 80 ns long MD trajectories were used to sample the total energy of the system in the different redox states and to carry out a statistical analysis of energy differences of the system in the reactant and product states. The energy differences were evaluated using Mathematica 10.339 and allowed to establish the free energy profile for the charge transfer reactions in Model I.

Results and Discussion

Earlier studies10,11,14 suggest that two QH2 binding motifs at the Qo-site of the bc1 complex from Rhodobacter capsulatus are plausible, differing in the protonation state of the residue H156 of the ISP, as shown in Fig. 2. For both binding motifs, referred to as Model I for H156 deprotonated and Model II for H156 protonated, the primary electron and proton transfer reactions at the Qo-site were studied through an in-depth MD and QC analysis, allowing to establish the energetics of the reactions and the nature of the proton and electron transfers. For this purpose TSs, corresponding to different charge transfer pathways, were obtained for Model I and Model II. From the calculated TSs, the charge transfer reaction pathways were revealed by scanning the potential energy surfaces along the charge transfer reaction coordinate. For a particular case of PCET reaction, analysis of the charge delocalization, electrostatic potential at the Qo-site and free energy calculations allowed to identify the key residues involved in the charge transfer reaction and to establish the corresponding driving forces.

TS of the charge transfer reaction

The TS for a molecular system separates the final state (product) from the initial state (reactant) on the potential energy landscape of a chemical reaction. Thus, once the TS of a reaction is identified, it is possible to follow a reaction coordinate in two opposite directions, starting from the TS and, thereby, to reconstruct the energy profile of the reaction. Finding the reactant, transition and product states, is key in the description of the reaction mechanism as it allows to establish the activation energy EA and the reaction energy ΔE, and hence to characterize the reaction itself.

By performing QC geometry optimizations over all the states depicted in Schemes (1) and (2), investigations of possible sequential and simultaneous electron and proton transfers at the Qo-site of the bc1 complex were carried out. Despite all the efforts, intermediate states, in which either the electron or the proton had been transferred, could not be established since the system always relaxed to either the initial state, in which neither charge has been transferred, or to the final state, in which both charges have been transferred simultaneously; i.e., only reactant and product states could be established. These results hint strongly that only a simultaneous reaction, in which the electron and the proton are transferred in a concerted manner, is feasible at the Qo-site, as indicated by diagonal arrows in Schemes (1) and (2), and that the reaction coordinate corresponds to the diagonal route in the Schemes.

In contrast, the TSs for the concerted reactions were established for Model I and Model II, thus allowing to reconstruct the reactant and product states and revealing the energy landscapes of PCET reactions. Figure 3 shows the energy profiles calculated for the two models of the Qo-site and indicates the TSs, which corresponds to the highest energy values of the profile in each model, as well as the initial and final electronic configuration of the Qo-site.

Figure 3
Energetics of the primary charge transfer reactions at the Qo-site of the bc1 complex

As indicated in Fig. 3, Model I reactants correspond to a bound QH2 and ISP residues and cofactors (labeled as ISP…QH2), while Model II reactants correspond to a bound QH2 to ISP(H+), representing the ISP subunit with a protonated H156 residue and, additionally, a quinol-binding H2O molecule (labeled as ISP(H+)…QH2+H2O). Products in both models correspond to a radical semiquinone QH, an ISP(H) and, only in the case of Model II, an additional hydronium ion H3O+ is present.

Energetics of the charge transfer reactions

The intrinsic differences in the molecular structure of the charge acceptors between the two studied computational models, determine the pathway of the proton transfer from QH2 to the bc1 complex. In Model I the proton is transferred to H156 of the ISP, leading to the formation of a radical ISP(H); however, as per a protonated H156 in Model II, the proton is transferred to an acceptor H2O molecule, leading to the creation of a hydronium ion H3O+, in addition to the ISP(H) radical that is being formed by the electron transfer to the initial ISP(H)+. Such differences in the reaction mechanism manifest as differences in the reaction energy profiles, demonstrated by the activation energy, EA, and reaction energy, ΔE, in Fig. 3.

In both considered models, the energy of the reactant state is lower than the energy of the product state, indicating that the associated charge transfer reactions are uphill processes. This means that the reactions can occur backwards, bringing the charges back to the initial QH2 donor, even after the initial transfers have occurred. The probability of this back transfer to happen is, however, higher in the case of Model II than in the case of Model I, were the difference EA−ΔE is considerably higher (see Table 2).

Table 2
Activation energy and reaction energies. Energies, indicated in Fig. 3, were computed for Model I and Model II through QC calculations carried out with the B3LYP/6-311G(d) and B3LYP/6-31G(d) methods. The B3LYP/6-31G(d) values are shown in the brackets. ...

Table 2 summarizes the activation and reaction energies obtained from the QC calculations, as well as the difference EA−ΔE for both studied models of the Qo-site. The activation energies are 14.20 kcal/mol and 8.94 kcal/mol for Model I and Model II, respectively. This comparison suggests that the PCET is more likely to occur in a configuration of the system described by Model II, with a protonated H156 residue. Such difference has a direct impact on the kinetics of the charge transfer reactions,40 leading to a higher PCET rate constant in the case of Model II than in the case of Model I. However, the rate constants of the reverse charge transfer reaction have to be taken into account in order to establish the stability of the PCET process. The rate constant for the reverse reaction in the case of Model II also depends of the diffusivity of H3O+, which is the proton carrier in this case. The diffusion of the hydronium ion away from the Qo-site could assist the prevention of the proton transfer in the reverse reaction, even though the low energy barrier of the reverse reaction in the case of Model II makes it favorable to occur.

As listed in Table 2, the barriers for the back reactions, EA−ΔE, differ for both models, being equal to 11.69 kcal/mol for Model I and 0.52 kcal/mol for Model II, respectively, strongly affecting the reaction kinetics, namely, the stability of the product state after the PCET. In Model I, EA−ΔE is considerably larger than this in the case of Model II: a small energy barrier in the case of Model II reveals that the system is equally likely to populate the TS and the product state, making the latter a rather unstable state of the system, and allowing a proton transfer back towards its initial donor QH2.

The observed difference between the two models is primarily attributed to the fact that the proton acceptor in the case of Model II is a solvent H2O molecule, as opposed to Model I where the acceptor is the H156 residue. Even though Model I seems to represent a more stable product state, it is possible that the hydronium is necessarily formed as an initial proton transporter to be translocated across the membrane. Previous studies10,14 suggest this possibility and here it is evidenced to still be a plausible scenario.

Table 2 summarizes the energies calculated by using the B3LYP/6-31G(d) method with a double-zeta precision basis set. All calculations were carried out using the triple-zeta B3LYP/6-311G(d) and double-zeta B3LYP/6-31G(d) basis sets for the purpose of testing the accuracy of different methods. A comparison of the calculations for the two basis sets is given in the SI. A relatively small difference in energies, shown in Table 2, computed with different methods, suggests that both methods could actually be used to describe the studied PCET reactions.

Reactant and product states of the PCET

Structurally, the TS of the PCET occurring at the Qo-site, corresponds to a quinol-donor hydrogen H2 shift towards its acceptor atom, NE2 of H156 in the case of Model I or OH2 of the H2O molecule in the case of Model II, accompanied by a redistribution of the transferring electron between the donor, QH2, and acceptor, Fe2S2. However, in the reactant and product states, the transferring charges are well localized in their specific donor and acceptor residues. In the reactant state, the electron is localized at the donor, QH2 site, and the proton (H2 atom) is bonded to it, while in the product state, the proton is bonded to the acceptor, H156 or H2O (in Model I and Model II, respectively), and the electron is localized at the Fe2S2 cluster site. Snapshots of the product, transition and reactant states in Model I are depicted in Fig. 4, where the hydrogen atom, highlighted as a magenta sphere, depicts the transferring proton, and the calculated highest occupied molecular orbital (HOMO), in orange and blue, illustrates the transferring electron.

Figure 4
Electron and proton transfers at the Qo-site of the bc1 complex

The hydrogen atom is localized equidistant from its donor, QH2, and acceptor, H156, in the TS (middle panel), and appears bonded to the donor (upper panel) and the acceptor (lower panel) in the reactant and product states, respectively. The orange arrow shown atop the hydrogen atom, corresponds to the imaginary normal vibration mode, obtained through QC calculations for the TS. A single imaginary frequency in the normal vibration spectrum was revealed by the B3LYP/6-311G(d) calculation, indicating a first order TS which hints on the proper reaction coordinate, and indicates a bond breakage of the donor-proton bond, as expected in a TS.41 The vibration mode involves the motion of the H-atom along the line connecting the donor and acceptor sites. The delocalization of the HOMO evidences a partially transferred electron in the TS, while a well localized HOMO at the donor and the acceptor sites can be seen for the reactant and product state, respectively.

Molecular orbitals redistribution upon charge transfer

The coupled nature of the primary charge transfer reaction at the Qo-site of the bc1 complex, can be described through a quantum mechanical analysis of the electronic structure and atomic spatial distribution studied along the reaction path.

For this purpose, the electronic distribution at the Qo-site was calculated for the 31 selected configurations, following the reaction profile in Fig. 3. The HOMO, computed for each configuration, allowed to visualize the distribution of the valence electrons at the donor and acceptor sites during the charge transfer reaction. Most of the remaining occupied molecular orbitals exhibit minor perturbations. The HOMO calculated for the 31 configurations, is displayed in a movie (see bc1.mpg in SI), which features the reaction pathway dynamically. Figure 4 shows the HOMO, calculated for the initial, transition, and final states of the Qo-site in Model I. Initially, the HOMO is localized around the QH2 substrate head group, while it shifts towards the ISP upon the charge transfer reaction; TS features electron delocalization between QH2 and ISP, as it is expected during an electron transfer.42

HOMO delocalization in the TS of the two computational models, suggests that in the course of the charge transfer reaction the system features charge delocalization throughout the different residues of the Qo-site. An analysis of the atomic charge assignments provides a more quantitative description of such delocalization and hence a more accurate reaction mechanism characterization.

Charge distribution at the Qo-site

Through QC calculations, the atomic charges were obtained for each of the 31 selected configurations. Figure 5 shows the total ESP-fitted charges computed for the residues of the Qo-site in the case of Model I and Model II. With exception of the transferred proton charge, the atomic charges are summed into fragments to illustrate the effect of charge exchange at the Qo-site during the charge transfer reactions. The residues at the Qo-site are colored by fragments (upper panel), for which the total charge is calculated as a function of the reaction coordinate (lower panels). In both models charge delocalization is evidenced given a specific molecular structure: in Model I a negative charge is transferred from the QH2 substrate (blue) to the ISP (magenta), while keeping a nearly constant charge distribution in the Y147 (red) and the E295 (green) residues. The transference of charges during the reaction is described through step-wise changes in the total charge of the donor and acceptor molecules upon the change of the reaction coordinate.

Figure 5
Charge delocalization at the Qo-site of the bc1 complex

In Model I, one notes a symmetrical change of charge of the total charge of the QH and the ISP fragments, see blue and magenta and lines in Fig. 5, which provides the evidence of the electron transfer reaction between the QH2 and the ISP fragments. In Model II, however, the existence of a coordinating H2O molecule in the Qo-site affects the charge transfer between QH2 and ISP, showing a slight asymmetry in the total charge of these fragments. Furthermore, a slight increase of the orange line indicates that the total charge of the H2O fragment (proton acceptor) increases simultaneously as the charge of the ISP fragment (electron acceptor) decreases.

The difference in the initial total charge of the ISP fragment for both computational models is due to the difference in the protonation state of the residue H156, which makes the initial total charge of the ISP in Model I ~ −1.0 e, and ~ 0 e in Model II. In both models, the Y147 and E295 residues maintain a highly conserved charge for all the reaction coordinate values, indicating that these residues do not act as charge donor or acceptors in the course of the reaction.

The described charge exchange between donor and acceptor residues at the Qo-site, furthermore evidences a PCET process in the Qo-site of the bc1 complex which is expected to be driven by electrostatic and thermodynamic effects of the environment that surrounds the Qo-site.

Effect of environment on PCET at the Qo-site

The total energy of the system measured for different redox states of the Qo-site allows to establish the free energy of the PCET reaction. For this purpose, the total energy of the system in the reactant state must be established, and compared to the energy of the system in the product state. In the present investigations, MD simulations were performed, considering the system as it resembles the reactant and the product state, and a reaction coordinate, ΔE, was defined as the energy difference between the reactant and the product state. This definition of the reaction coordinate is in general used for free energy calculations of charge transfer reactions43 as it allows to describe the free energy in terms of energy differences between reactant and product states.

Since the PCET reaction is essentially a quantum mechanical process taking place at the Qo-site, one needs to differentiate the Qo-site from its environment in the different redox states of the bc1 complex. The Qo-site model is defined here according to its involvement in the quantum mechanical process and it includes the atoms that have been considered in the QC calculations, as shown in Fig. 2. The environment surrounding the Qo-site, on the other hand, is composed of the remaining protein, lipids, water, and ion atoms that have been considered in the MD simulations, see Fig. 6a. The total energy of the system can, therefore, be described as the sum of the energy of the Qo-site EQo, the energy of the environment EEnv, and the interaction energy between the Qo-site and the environment, Eint, as:

Figure 6Figure 6Figure 6
Free energy differences of the PCET

The total energy of the system is then calculated from the classical MD simulation, using the force field approximation, for every frame of the MD trajectories. However, an additional correction to the energy has to be taken into account as the PCET could not be described by classical mechanics, but instead, quantum mechanically. This means that the energy of the Qo-site in Eq. (3), obtained from MD simulations, should be replaced by the QC energy of the Qo-site, and thus the total energy of the system is:


Here, the superscript (MD) indicates that the energy was calculated from MD trajectories, while the superscript (QC) indicates that this was obtained from QC calculations.

In order to obtain the reaction coordinate for the free energy calculation, i.e., the energy differences ΔE between the reactant and the product states, the calculation of the total energy in Eq. (4), is performed for the reactant and the product states of the system for every frame in the refined 80 ns long MD trajectories (see Table 1).

For every frame of the MD trajectory in which the system is in the reactant state, one considers two redox states of the Qo-site, while keeping the environment in the reactant state. The energy difference, thus, corresponds to the total energy difference between the two configurations of the system where the Qo-site has been changed upon the PCET reaction, while the environment did not have enough time to respond to this changes. The energy difference thus reads as:


where the first letter in the superscript indicates the redox state of the environment (R[equivalent]reactant), and the second superscript corresponds to the redox state of the Qo-site (R[equivalent]reactant; P[equivalent]product). The energy EtotalRR is readily obtained from MD simulations corrected through QC calculations employing Eq. (4). To calculate the energy EtotalRP it is required to set the atomic charges of the atoms in the Qo-site to the values that correspond to the product state, while preserving the positions and charges of the atoms of the environment as in the reactant state simulations. In other words, EtotalRP could be obtained once the system resembles an environment of the reactant state and the Qo-site of the product state. Analogous calculations are carried out for the product state, in which case the energy difference is:


Once the energy differences ΔER and ΔEP are obtained for the 80 ns long MD trajectories of the reactant and product states, it is possible to compare the probability distribution of the energy differences for both states of the bc1 complex, as shown in Fig. 6b. The distributions pER) and pEP) are expected to follow the Gaussian profile


where the superindex i corresponds to R (reactant) or P (product) states, σi is the width and μi the average energy difference of the distribution pEi). The free energy could be readily calculated once pER) and pEP) are known:43


Here the superscript i stands for P (product) and R (reactant) states. Fig. 6c shows the resulting free energy curves. The definition of the free energy in Eq. (7) implies that the energy profiles for the reactant and product states cross at ΔE = 0; this is achieved by shifting the energy of one state in reference to the other. The resulting free energy profiles presented as a function of the energy difference indicate that the rate of the backwards PCET process differs from the rate of the forward PCET process, which could be concluded from the asymmetry of the two energy curves in Fig. 6c. Moreover, the relative position of the curves minima indicates that the PCET appears to be energetically a downhill process. This result should be differentiated from the energy profile described by the pure QC calculations in Fig. 3, as the later only describes a single conformation of the system, while the free energy calculations take into account the environment of the Qo-site and provide a statistical averaging over a considerable number of possible configurations, that the system could populate.

The free energies obtained in this study are intended to further characterize the primary PCET at the Qo-site of the bc1 complex, and can be used to obtain the rate constants of the underlaying charge transfer processes. However, further calculations of the PCET rate constants require additional information such as establishing the adiabaticity regime in which the reactions occur, the coupling of the quantum states that describe the system in the reactant and product states, and possibly, extending MD simulations as well as performing multiple calculations for the different configurations of the Qo-site.

Key role of E295 and Y147

In a previous study,10 based on MD simulations and QC computations of the Qo-site of the bc1 complex, it was demonstrated that the residues E295 and Y147 feature rearrangements to form a hydrogen bonding network with QH2. In the present MD simulations Y147 residue occasionally turns away from the QH2 head group, letting a water molecule to occupy its place instead. In order to study the specific role of E295 and Y147 in the PCET process, the Qo-site Model I was modified such that (i) both residues (E295 and Y147) were removed, and (ii) the Y147 residue was replaced by a water molecule.

The TS of the PCET reaction could not be established once the Qo-site is missing the E295 residue, indicating that in this case the charge transfer reaction is energetically unfavorable, or even impossible. On the contrary, once E295 is present, but Y147 is replaced by H2O molecule the TS could be found from QC calculations. The energy profile of the corresponding charge transfer reaction is shown in Fig. 7, where the optimized structures for the reactant and product states are indicated in the upper panels.

Figure 7
Minor impact of Y147 on the primary charge transfer reaction

The findings provide further support for previous investigations in which it was suggested that E295 acts as a proton acceptor, while Y147 does not play a fundamental role in the QH2 binding or proton transfer from QH2. However, all calculations where Y147 is present indicate a mediation of this residue in the QH2 binding as well as in the subsequent proton transfer. The findings thus strongly suggest that Y147 acts as an intermediate bridge for the proton transfer between QH2 and E295.

It is remarkable that the energy profile in Fig. 7 is largely similar to the energy profile calculated for the complete Qo-site Model I shown in Fig. 3, which indicates that the charge transfer reactions at the Qo-site are possible even in the absence of Y147 at the Qo-site. The E295 residue, however, is essential for the reaction to occur, as its mutation could lead to the prevention of the primary PCET reaction.

Electrostatic potential distribution at the Qo-site

Since the primary charge transfer reaction at the Qo-site of the bc1 complex seems to be affected by the E295 residue, one should expect that its electrostatic properties should contribute greatly to the generation of a proper electrostatic environment for rendering the primary charge transfer reactions from QH2 to the ISP.

Earlier 360 ns long MD simulations10 were employed here to compute the time averaged electrostatic potential of the Qo-site for the two studied models. Figure 8 shows two well defined positive and negative electrostatic potential regions around the QH2 head group: the negative potential (red surface) is strategically centered around the side chain of E295 and the positive potential (blue surface) embraces the iron-sulfur cluster. This particular electrostatic potential distribution can produce a driving force for the valence electron of the QH2 head group, indicating that it may be the key element that enables its transfer to the ISP.

Figure 8
Electron transfer driving force

As the repulsive force generated by the negative electrostatic potential, localized near the E295 residue, drives the valence electron towards the acceptor Fe2S2, it would seem that the PCET is initiated by the electron transfer that consequently drives a proton transfer. However, a more detailed analysis is called for at this point, in which the adiabaticity of the electron and proton transfers are accurately evaluated considering the present results.


The present study evidences the coupled nature of the primary proton and electron transfer (PCET) reactions that initiate the Q-cycle at the Qo-site of the Rhodobacter capsulatus bc1 complex. This PCET process was established through the detailed analysis of reaction energetics, computed quantum mechanically for possible charge transfer reactions from the QH2 to the Fe2S2 cluster of the ISP subunit in the bc1 complex, featuring two different models of the Qo-site, which differ in the protonation state of the key H156 residue. Particularly, for the deprotonated H156 model, which seems to support a more stable reaction, the charge delocalization at the Qo-site indicates that electron and proton from the QH2 molecule are transferred in tandem, driven by a specific electrostatic potential distribution at the Qo-site.

The involvement of key residues, such as H156 of the ISP subunit and Y147 and E295 of cyt. b subunit, in the primary PCET reaction is further elaborated, as previous MD analysis had indicated their importance earlier.10 The Y147 and E295 residues rearrange to form a hydrogen bonding with the QH2, and assist the second proton transfer from QH2, which is expected to flow, towards Y147 and subsequently E295.

Free energy calculations obtained for the PCET process at the Qo-site showed that the molecular environment of the Qo-site plays an important role in the PCET energetics and the driving force of the reaction. The performed free energy calculations illustrate the charge transfer processes at the Qo-site, but also suggest that in order to obtain an accurate estimate of the PCET reaction rate constants, it is necessary to reveal the adiabaticity of the proton and electron transfers. This can for example be achieved by exploring the Qo-site of the bc1 complex through hybrid computational methods such as QM/MM44 or a polarizable embedding approach.45

Although the present investigation reveals the coupled nature of the primary charge transfer reactions at the Qo-site of the bc1 complex, the studied physical mechanism provides only a first step towards describing fundamentally important ubiquitous mechanisms of energy transport in cellular respiration and higher photosynthetic organisms, and strives for follow up investigations.

Supplementary Material

Supplementary material

Supplementary movie

topology file


The research reported here has been supported by the National Institutes of Health through grant NIH 9P41GM104601 and the National Science Foundation through grant NSF PHY0822613 to KS. The authors acknowledge supercomputer time on Stampede provided by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin through Extreme Science and Engineering Discovery Environment (XSEDE) Grant XSEDE MCA93S028 and the DeIC National HPC Centre at SDU. AMB is grateful to Beckman Institute for financial support. IAS is grateful for financial support from the Lundbeck Foundation.


Supporting Information Available

In Supporting Information (SI) we present a movie of the PCET reaction, the supplementary topology files used in the present investigations, and additional figures.


1. Wikstrom M, Krab K, Saraste M. Proton–Translocating Cytochrome Complexes. Annu Rev Biochem. 1981;50:623–655. [PubMed]
2. Nicholls DG, Ferguson S. Bioenergetics. Academic Press; 2013.
3. Crofts AR, Berry EA. Structure and Function of the Cytochrome bc1 Complex of Mitochondria and Photosynthetic Bacteria. Curr Opin Struct Biol. 1998;8:501–509. [PubMed]
4. Cramer WA, Hasan SS, Yamashita E. The Q-cycle of Cytochrome bc Complexes: A Structure Perspective. Biochim Biophys Acta – Bioener. 2011;1807:788–802. [PMC free article] [PubMed]
5. Crofts AR. Proton-coupled Electron Transfer at the Qo-site of the bc1 Complex Controls the Rate of Ubihydroquinone Oxidation. Biochim Biophys Acta – Bioener. 2004;1655:77–92. [PubMed]
6. Hsueh K-L, Westler WM, Markley JL. NMR Investigations of the Rieske Protein from Thermus Thermophilus Support a Coupled Proton and Electron Transfer Mechanism. J Am Chem Soc. 2010;132:7908–7918. [PMC free article] [PubMed]
7. Chang CJ, Chang MC, Damrauer NH, Nocera DG. Proton-Coupled Clectron Transfer: a Unifying Mechanism for Biological Charge Transport, Amino Acid Radical Initiation and Propagation, and Bond Making/Breaking Reactions of Water and Oxygen. Biochim Biophys Acta – Bioener. 2004;1655:13–28. Special issue dedicated to Jerry Babcock. [PubMed]
8. Sjödin M, Styring S, Åkermark B, Sun L, Hammarström L. Proton-Coupled Electron Transfer From Tyrosine in a Tyrosine-Ruthenium- Tris-Bipyridine Complex: Comparison with Tyrosinez Oxidation in Photosystem II. J Am Chem Soc. 2000;122:3932–3936.
9. Berry EA, Huang L-S, Saechao LK, Pon NG, Valkova-Valchanova M, Dal-dal F. X-ray Structure of Rhodobacter Capsulatus Cytochrome bc1: Comparison with its Mitochondrial and Chloroplast Counterparts. Photosyn Res. 2004;81:251–275. [PubMed]
10. Barragan AM, Crofts AR, Schulten K, Solov’yov IA. Identification of Ubiquinol Binding Motifs at the Qo-Site of the Cytochrome bc1 Complex. J Phys Chem B. 2015;119:433–447. [PMC free article] [PubMed]
11. Crofts AR, Hong S, Ugulava N, Barquera B, Gennis R, Guergova-Kuras M, Berry EA. Pathways for Proton Release During Ubihydroquinone Oxidation by the bc1 Complex. Proc Natl Acad Sci USA. 1999;96:10021–10026. [PubMed]
12. Crofts AR, Hong S, Wilson C, Burton R, Victoria D, Harrison C, Schulten K. The Mechanism of Ubihydroquinone Oxidation at the Qo-site of the Cytochrome bc1 Complex. Biochim Biophys Acta – Bioener. 2013;1827:1362–1377. [PMC free article] [PubMed]
13. Crofts AR, Barquera B, Gennis RB, Kuras R, Guergova-Kuras M, Berry EA. Mechanism of Ubiquinol Oxidation by the bc1 Complex: Different Domains of the Quinol Binding Pocket and Their Role in the Mechanism and Binding of Inhibitors. Biochemistry. 1999;38:15807–15826. [PubMed]
14. Postila PA, Kaszuba K, Sarewicz M, Osyczka A, Vattulainen I, Róg T. Key Role of Water in Proton Transfer at the Qo-site of the Cytochrome bc1 Complex Predicted by Atomistic Molecular Dynamics Simulations. Biochim Biophys Acta – Bioener. 2013;1827:761–768. [PubMed]
15. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel RD, Kalé L, Schulten K. Scalable Molecular Dynamics with NAMD. J Comp Chem. 2005;26:1781–1802. [PMC free article] [PubMed]
16. Mackerell AD, Feig M, Brooks CL. 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 Comp Chem. 2004;25:1400–1415. [PubMed]
17. Frisch MJ, et al. Gaussian 09, revision A 1. Gaussian Inc.; Wallingford, CT: 2009.
18. Becke AD. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J Chem Phys. 1993;98:5648–5652.
19. Shoji M, Koizumi K, Kitagawa Y, Yamanaka S, Okumura M, Yamaguchi K. Theory of Chemical Bonds in Metalloenzymes IV: Hybrid-DFT Study of Rieske-type [2Fe-2S] Clusters. Int J Quant Chem. 2007;107:609–627.
20. Bassan A, Blomberg MR, Borowski T, Siegbahn PE. Oxygen Activation by Rieske Non-heme Iron Oxygenases, a Theoretical Insight. J Phys Chem B. 2004;108:13031–13041.
21. Salomon O, Reiher M, Hess BA. Assertion and Validation of the Performance of the B3LYP Functional for the First Transition Metal Row and the G2 Test Set. J Chem Phys. 2002;117:4729–4737.
22. Niu S, Nichols JA, Ichiye T. Optimization of Spin-unrestricted Density Functional Theory for Redox Properties of Rubredoxin Redox Site Analogues. J Chem The Comp. 2009;5:1361–1368. [PMC free article] [PubMed]
23. Sigfridsson E, Olsson MH, Ryde U. A Comparison of the Inner-sphere Reorganization Energies of Cytochromes, Iron-sulfur Clusters, and Blue Copper Proteins. J Phys Chem B. 2001;105:5546–5552.
24. Wilkens SJ, Xia B, Weinhold F, Markley JL, Westler WM. NMR Investigations of Clostridium Pasteurianum Rubredoxin. Origin of hyperfine 1H, 2H, 13C, and 15N NMR Chemical Shifts in Iron-sulfur Proteins as Determined by Comparison of Experimental Data with Hybrid Density Functional Calculations. J Am Chem Soc. 1998;120:4806–4814.
25. Sigfridsson E, Olsson MH, Ryde U. Inner-sphere Reorganization Energy of Iron-sulfur Clusters Studied with Theoretical Methods. Inor Chem. 2001;40:2509–2519. [PubMed]
26. Humphrey W, Dalke A, Schulten K. VMD – Visual Molecular Dynamics. J Mol Graphics. 1996;14:33–38. [PubMed]
27. Kaszuba K, Postila PA, Cramariuc O, Sarewicz M, Osyczka A, Vattulainen I, Róg T. Parameterization of the Prosthetic Redox Centers of the Bacterial Cytochrome bc 1 Complex for Atomistic Molecular Dynamics Simulations. Theor Chem Acc. 2013;132:1–13.
28. Feller SE, Yin D, Pastor RW, MacKerell A. Molecular Dynamics Simulation of Unsaturated Lipid Bilayers at Low Hydration: Parameterization and Comparison with Diffraction Studies. Biophys J. 1997;73:2269–2279. [PubMed]
29. Aguayo D, González-Nilo FD, Chipot C. Insight into the Properties of Cardiolipin Containing Bilayers from Molecular Dynamics Simulations, Using a Hybrid All-Atom/United-Atom Force Field. J Chem Theor Comp. 2012;8:1765–1773. [PubMed]
30. Aksimentiev A, Schulten K. Imaging alpha-Hemolysin with Molecular Dynamics: Ionic Conductance, Osmotic Permeability, and the Electrostatic Potential Map. Biophys J. 2005;88:3745–3761. [PubMed]
31. Szilagyi RK, Winslow MA. On the Accuracy of Density Functional Theory for Iron-sulfur Clusters. J Comp Chem. 2006;27:1385–1397. [PubMed]
32. Peng C, Bernhard Schlegel H. Combining Synchronous Transit and Quasi-Newton Methods to Find Transition States. Isr J Chem. 1993;33:449–454.
33. Peng C, Ayala PY, Schlegel HB, Frisch MJ. Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States. J Comp Chem. 1996;17:49–56.
34. Fukui K. The Path of Chemical Reactions - the IRC Approach. Acc Chem Res. 1981;14:363–368.
35. Hratchian HP, Schlegel HB. In: Theory and Applications of Computational Chemistry. Dykstra CE, Frenking G, Kim KS, Scuseria GE, editors. Elsevier; Amsterdam: 2005. pp. 195–249. Chapter 10.
36. Singh UC, Kollman PA. An Approach to Computing Electrostatic Charges for Molecules. J Comp Chem. 1984;5:129–145.
37. Besler BH, Merz KM, Kollman PA. Atomic charges derived from semiempirical methods. J Comp Chem. 1990;11:431–439.
38. MacKerell AD, et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J Phys Chem B. 1998;102:3586–3616. [PubMed]
39. Wolfram Research Inc. Mathematica, Version 10.3. Champaign, Illinois: 2015.
40. Lodish H, Berk A, Zipursky SL, Matsudaira P, Baltimore D, Darnell J. Molecular Cell Biology. 4th. WH Freeman; New York: 2000.
41. Truhlar DG, Isaacson AD, Garrett BC. In: Theory of Chemical Reaction Dynamics. Baer M, editor. CRC Press; Boca Raton, FL: 1985. pp. 65–137. Chapter 4.
42. Hammond BL, Lester WA, Jr, Reynolds PJ. Monte Carlo Methods in ab initio Quantum Chemistry. Vol. 1 World Scientific; 1994.
43. Blumberger J. Free Energies for Biological Electron Transfer from QM/MM Calculation: Method, Application and Critical Assessment. Phys Chem Chem Phys. 2008;10:5651–5667. [PubMed]
44. Warshel A, Levitt M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium ion in the Reaction of Lysozyme. J Mol Biol. 1976;103:227–249. [PubMed]
45. Olsen JM, Aidas K, Kongsted J. Excited States in Solution through Polarizable Embedding. J Chem Theor Comp. 2010;6:3721–3734.