Identification of ubiquinol binding motifs at the Qo-site of the cytochrome bc1 complex.

Enzymes of the bc1 complex family power the biosphere through their central role in respiration and photosynthesis. These enzymes couple the oxidation of quinol molecules by cytochrome c to the transfer of protons across the membrane, to generate a proton-motive force that drives ATP synthesis. Key for the function of the bc1 complex is the initial redox process that involves a bifurcated electron transfer in which the two electrons from a quinol substrate are passed to different electron acceptors in the bc1 complex. The electron transfer is coupled to proton transfer. The overall mechanism of quinol oxidation by the bc1 complex is well enough characterized to allow exploration at the atomistic level, but details are still highly controversial. The controversy stems from the uncertain binding motifs of quinol at the so-called Qo active site of the bc1 complex. Here we employ a combination of classical all atom molecular dynamics and quantum chemical calculations to reveal the binding modes of quinol at the Qo-site of the bc1 complex from Rhodobacter capsulatus. The calculations suggest a novel configuration of amino acid residues responsible for quinol binding and support a mechanism for proton-coupled electron transfer from quinol to iron-sulfur cluster through a bridging hydrogen bond from histidine that stabilizes the reaction complex.

Introduction

The biosphere is sustained by chemiosmotic circuits, driven by light (in photosynthesis) or chemical potential difference (in respiration), through which electron transfer reactions coupled to proton transfer across a cellular membrane maintain the proton gradient that drives cellular metabolism through synthesis of adenosine triphosphate (ATP), transport of ions and metabolites, and other processes.[1−4] Specialized proteins and protein complexes facilitate electron transfer between distinct electron donor and acceptor sites. In the case of cellular respiration or photosynthesis, a multiprotein enzyme that drives proton translocation across the cellular membrane while performing electron transfer is the cytochrome bc1 complex.[5,6] The overall reaction of the bc1 complex, also referred to as the Q-cycle, results in the net oxidation of two quinol (QH2) molecules with the release of four protons to the positive side of the membrane (or intermembrane space in the case of mitochondria), and reduction of one quinone (Q) with uptake of two protons from the negative side.[7]

Numerous bc1 complex structures have been resolved through X-ray crystallographic methods, namely bovine,[8] chicken,[4] yeast,[9] and bacterial[10] ones. The functional core of bacterial and mitochondrial bc1 complexes consists of a dimeric protein arrangement, where each monomer contains one cytochrome b (cyt. b), one cytochrome c1 (cyt. c1), and one iron–sulfur protein (ISP) subunit, as illustrated in Figure 1a. These subunits bind several prosthetic groups that are involved in electron transfer through the bc1 complex and, therefore, are crucial for the Q-cycle. In particular, the cyt. b subunit contains two iron-containing heme b groups, the cyt. c1 subunit contains a heme c group, and the ISP subunit contains an Fe2S2 cluster, as shown in Figure 1b.

Figure 1

bc1 complex from Rhodobacter capsulatus. (a) The studied molecular system consists of a lipid bilayer membrane, water molecules, ions, and the six protein subunits forming the homodimeric bc1 complex. The bc1 complex features two monomers (A and B), each consisting of one cytochrome c1 (cyt. c1), one cytochrome b (cyt. b), and one iron–sulfur protein (ISP) subunit. (b) Each monomer binds four metal centers, heme c1 in the cyt. c1 subunit, and hemes bL and bH in the cyt. b subunit, while the ISP binds an iron–sulfur (Fe2S2) cluster. The quinol (QH2) and the quinone (Q) substrate molecules interact with the hemes and the Fe2S2 cluster to facilitate electron and proton transfers through the complex (A or B) at two distinct binding sites (Q and Q). The arrows indicate schematically pathways of electrons and protons at the initial phase of the Q-cycle. (c) The QH2 substrate molecule at the Q-site of the bc1 complex interacts closely with the H156 residue of the ISP and several other residues of the bc1 complex. The exact binding mode is addressed in the present study.

Each monomer of the bc1 complex (see Figure 1a) includes two active sites that govern the functional capabilities of the entire complex and are referred to as the Q- and Q-sites (see Figure 1b). Prior to and during the Q-cycle, the substrate molecules QH2 and Q bind to the bc1 complex at two distinct sites: the Q-site captures two QH2 and transforms them to Q, while the Q-site captures one Q and transforms it to QH2. These transformations involve electron and proton exchange between the substrate molecules and the proteins of the bc1 complex.

Even though the mechanism that governs the charge transfer reactions in the Q-cycle has been investigated for decades,[11] the limitations of the experimental techniques and the complexity of the entire process make it difficult to resolve the physicochemical mechanism underlying the Q-cycle. In fact, the individual electron and proton transfer pathways in the bc1 complex are still controversial.[7]

The controversy follows from the following difficulty. The rate-limiting reaction proceeds from an initial complex (the enzyme–substrate complex) formed only under metastable conditions and this complex is by its nature inaccessible to crystallography; indeed, no structures have shown any quinone species in the Q-site. However, a substantial literature has explored physicochemical parameters for kinetics of the reaction in wildtype and mutant strains, from which a picture of the molecular architecture of the site can be derived.[12−14] In addition, spectroscopic evidence has demonstrated a relatively stable complex between quinone and the reduced ISP. On the other hand, many structures with quinone analogues bound are available, which show several different configurations for the site, and likely coordinating residues.[15−17]

The electron and proton transfer reactions taking place in the bc1 complex can be resolved into partial processes including intermediate redox and protolytic reactions and chemical states.[18] Technical limitations and the metastable nature of the reaction complex leave many dynamical features inaccessible to direct measurement. For a complete understanding experiment needs to be complemented by computational modeling to identify physicochemical details of the mechanism at the atomistic level. A first attempt to combine crystallographic data with large scale classical molecular dynamics simulations was made early on,[19] but the quantum calculations needed to address the chemistry of catalysis were not feasible at that time. However, computational power is now much improved, and given the central role of the bc1 complex in all respiratory-based bioenergetics, a combined classical and quantum chemical description is a task that should be undertaken with urgency.

Naturally, the experimental–computational studies need to first focus on the initial state with QH2 bound to the Q-site and Q bound to the Q-site. The redox states of the bc1 complex residues at the Q- and Q-sites coordinating initially the substrates are crucial for the Q-cycle, as the coordination and redox states largely influence the rates of electron and proton transfers to and from the substrate molecules. Several coordinating residues of the substrate molecules have been proposed experimentally[20] and have been recently studied through molecular dynamics (MD) simulations.[21] In the case of the Q-site, the key residue for the binding of QH2 is a histidine (H156 in the case of bc1 complex from Rhodobacter capsulatus), covalently bonded to the Fe2S2 cluster of the ISP, as illustrated in Figure 1c. It is assumed that this histidine forms a hydrogen bond with the QH2 molecule upon its docking to the Q-site,[13] thereby keeping the substrate ready for the primary electron and proton transfers.

The chemical specificity of the QH2···H156 hydrogen bonding is, however, still debated and the actual protonation state of H156 greatly impacts the charge transfer reactions at the Q-site;[14,22] the protonation state active in normal forward chemistry can be tested by varying environmental pH in the range of the pK of H156 (see also the Supporting Information for a more specific discussion). To establish the binding mode of QH2 at the Q-site, it is, therefore, necessary to consider both protonation states of H156. In addition, one needs to identify other key residues that contribute to the binding of the QH2 substrate molecule within the bc1 complex. Particularly important, in this respect, is the inclusion of all charged and polar residues that can potentially contribute to substrate binding at the Q-site, as these residues can impact critically the rate of electron and proton transfers.

In the present study we investigate, through atomistic molecular dynamics (MD) simulations supported by quantum chemical density functional theory (DFT) calculations, the binding of QH2 at the Q-site of the bc1 complex from Rhodobacter capsulatus.[23] The MD simulations were performed for two states of the Q-site differing in the protonation state of histidine H156. Quantum chemical analysis allowed us to obtain a close view of QH2 binding at the Q-site. The results from the combination of classical and quantum chemical methods provide new insights into the role of the amino acid residues that hold the QH2 molecule at the Q-site of the bc1 complex. The reported findings are consistent with some features of earlier MD simulations,[21] but they identify also rearrangements of binding site residues not discussed previously and investigate two, rather than only one, H156 protonation state. The present study is thus a systematic key first step toward an atomic level investigation of the entire Q-cycle of the bc1 complex.

Methods

We have investigated by means of MD simulation and quantum chemistry calculations two possible binding modes of the QH2 substrate at the Q-site of the bc1 complex. The two respective simulations differ in the protonation state of the H156 residue of the ISP (see Figure 1c) that holds QH2 at the Q-site; H156 was assumed to be either in its protonated (Model I) or in its deprotonated (Model II) form.

MD simulations for the two computational models were performed through NAMD 2.9[24] assuming the CHARMM36 force field with CMAP corrections[25] for the proteins; for lipids and cofactors supplementary force fields[26] were employed, as discussed below. The quantum chemical calculations were carried out with the Gaussian 09 package[27] by using the DFT model chemistry. Analysis of results and snapshots of molecular structures were accomplished with VMD 1.9.1.[28]

System Preparation

The simulated systems considered in Models I and II were constructed, using VMD 1.9.1, from the X-ray crystal structure of the bc1 complex of Rhodobacter capsulatus (PDB ID: 1ZRT),[23] embedding the latter in a bilayer membrane, solvating protein and lipids within a TIP3P water box at a salt (NaCl) strength of 0.05 mol/L, and neutralizing the entire system with the salt ions. The bc1 complex forms a dimeric arrangement of six catalytic subunits,[29] each including cofactors that in the simulations were considered in the oxidized states, as summarized in Table 1.

Table 1

Oxidation States of the bc1 Complex Cofactors Assumed in the MD Simulationsa

subunit	cofactor	formal charge	oxidation state
cyt. b	heme bL	–1	oxidized
	heme bH	–1	oxidized
cyt. c1	heme c	–1	oxidized
ISP	Fe2S2	0	oxidized

Cofactors of monomers A and B of the bc1 complex, shown in Figure 1b, were simulated in identical oxidation states listed for all MD simulations performed.

The assumed oxidized form of all cofactors corresponds to the initial state of the bc1 complex prior to any charge transfer reaction. Charges and topologies of the bc1 complex proteins were assumed standard according to the CHARMM36 force field, while parameters for the cofactors were adopted to be consistent with an earlier study[21] and were adopted from an earlier investigation.[26]

Monomers A and B of the bc1 complex contain a Q-site and a Q-site, where QH2 and Q substrates become oxidized and reduced, respectively, during the Q-cycle. The 1ZRT crystal structure[23] of the bc1 complex includes bound stigmatellin molecules at the Q-sites of monomers A and B that were replaced in the performed simulations by QH2 molecules, aligning for this purpose the QH2 head groups with the respective head groups of stigmatellin, an approach also used in an earlier study.[21] The two Q molecules at the Q-sites were placed in the positions of antimycin molecules from the bc1 complex X-ray crystal structure of wild type Rhodobacter sphaeroides (PDB ID: 2QJP), antimycin being added, instead of Q, by the crystallographers for its inhibiting property.[30] Charges and topology of the QH2 and Q substrates were taken for the present study from an earlier investigation.[26]

The lipid bilayer was modeled as a random distribution of cardiolipin (CL 18:2/18:2/18:2/18:2), phosphatidylcholine (PC 18:2/18:2), and phosphatidylethanolamine (PE 18:2/18:2) lipids, with the concentration matching an earlier simulation;[21] the studied membrane patch included 102 CL, 406 PC, and 342 PE lipids. Since standard CHARMM36 parameters for CL are not available, the force field parameters from a prior study[31] were used for modeling the CL head group, while the parameters for the lipid tails were taken from the standard CHARMM36 force field. For modeling the PE and PC lipids, the standard CHARMM36 force field was employed.[32]

The missing hydrogen atoms from the X-ray structure of the bc1 complex were reconstructed by using the VMD plugin psfgen.[28] All histidine residues of the bc1 complex were considered as δ-protonated except for H156, which has been assumed ϵ-protonated in Model I and deprotonated 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 that were included in the computational models. The simulated system consisted of 502 165 atoms in Model I and 500 791 atoms in Model II, including proteins with cofactors, substrate molecules, lipids, water molecules, and ions.

The H156 residue in the Model II simulation was considered in its deprotonated form and needed to be specifically parametrized as no parameters were available for this residue protonation state in the context of the other residues in the binding site. For this purpose the complex of H156 and the Fe2S2 cluster together with the ligating residues C133, C153, and H135 was optimized using Gaussian 09, employing the B3LYP/6-31G(d) model chemistry.[33] The optimized structure was used to obtain the partial charges, which were determined by means of the electrostatic potential (ESP) fitting procedure;[34] the charges are given in the Supporting Information. The parameters of the bonding, angular, and dihedral interactions for deprotonated H156 were taken from the analogous parameters of the histidine residue in its standard protonation state.

Molecular Dynamics Simulations

MD simulations were performed with a time step of 2 fs. Electrostatic and van der Waals interactions were treated with a smooth cutoff of 12 Å. Long-range electrostatic interactions were calculated using the PME method, employing periodic boundary conditions.[35] The equilibration of the system was performed in the NPT ensemble, where the temperature was kept at 310 K by applying to all heavy atoms in the system Langevin forces with a damping coefficient of 5 ps–1. Pressure control was achieved during the equilibration simulations through Nosé–Anderson–Langevin piston pressure control[36] at 1 atm, using a piston oscillation period of 200 fs and a damping time scale of 50 fs. The production simulations were carried out in the NVT ensemble.

The protocol of the simulations performed, listed in Table 2, can be subdivided into two parts: (i) system equilibration and (ii) MD simulation. The equilibration was carried out in several steps. After energy minimization of the initial bc1 structure, lipids, water molecules, and ions were simulated for 60 ns, keeping atoms of the bc1 complex harmonically constrained and employing a combination of CHARMM22 and CHARMM27 force fields. The combination of force fields employed was the same as used in prior MD studies.[21] Following the first 60 ns, a 90 ns simulation with the CHARMM36 force field was performed, still keeping the entire bc1 complex constrained.

Table 2

Protocol for bc1 Complex Simulations Carried out in the Present Studya

	time interval (ns)
process	Model I	Model II
1. equilibration
structure minimization	50 000 NAMD steps
lipid bilayer, water molecules and ions released; rest constrained	60 (CH22 + CH27); 90 (CH36)
protein side chain released	70
turns, bridges, and coils motifs released	—	30
all atoms released	60	150
2. MD simulation
NVT ensemble	360

Simulations for Model I (502,165 atoms) and Model II (500,791 atoms) of the bc1 complex included an equilibration simulation and a MD simulation for the analysis of QH2 binding to the Q-site. CH22 + CH27 indicates a combined CHARMM22 and CHARMM27 force field, that was employed initially, while all consecutive simulations were done with the CHARMM36 force field. “NVT ensemble” denotes the canonical ensemble that was employed in MD simulations.

Next, the side chains of the bc1 complex were released and the system was equilibrated additionally for 70 ns. Finally, for the Model I simulation, all atoms were released and equilibrated for 60 ns, while for Model II the more flexible motifs of the secondary structure were kept constrained additionally for 30 ns prior to releasing all atoms and performing a 150 ns equilibration of a constraint-free system. After the equilibration, the MD simulation was carried out for 360 ns, for both computational models, in the NVT ensemble.

The duration of the equilibration was guided by monitoring the area of the membrane patch and the root-mean-square deviation (RMSD) of the bc1 complex proteins as these quantities needed to relax to constant values prior to equilibration. The relaxation of membrane patch area and RMSD are shown in the Supporting Information, Figures S2 and S3, respectively, for Model I, and Figures S5 and S6, respectively, for Model II.

Quantum Chemistry Calculations

The motif of the Q-site with bound QH2 was studied using the quantum chemistry package Gaussian 09,[27] employing the UB3LYP DFT method[33] for both Models I and II. This method has been widely used previously in optimizations of iron–sulfur containing systems.[37−43] Two standard 6-311G(d) and 6-311+G(d) basis sets were employed to expand the electronic wave functions. Both methods are of triple-ζ accuracy, while the latter includes additional diffuse functions.[44] For both Models I and II the quantum chemistry geometry optimizations included the QH2 head group; pre-equilibrated side chains of residues Y147, I292, E295, and Y302 of cyt. b; residues C133, C153, C155, H156, and H135 of ISP; and the Fe2S2 cluster of ISP, thereby, taking into account the major environmental effects that other hybrid methods have included likewise, but through a dielectric model rather than the explicit treatment of nearby side groups, for studying iron–sulfur cluster containing systems.[22] The structure of the Q-site was studied through quantum chemical energy minimization, where the positions of the Cα atoms were fixed to positions taken form the pre-equilibrated structure, to avoid an unphysical collapse of the Q-site model. For the quantum chemical calculations the Cα atoms of the side chain residues were actually replaced by CH3 groups, employing for this purpose the MOLEFACTURE plugin of VMD.[28]

For Model I, a water molecule has been suggested to play a key role in the binding of the quinol molecule to the Q-site.[21] and was included also in the present quantum chemical calculations; this water molecule is not stably bound in Model II and, hence, not included in this model. The quantum chemical studies for both models included all polar and charged residues at the Q-site and contained approximately 150 atoms.

Results

The binding of QH2 to the Q-site of the bc1 complex was characterized first through classical MD simulations. For this purpose the hydrogen bonds that QH2 forms at the binding site in the course of an MD trajectory were analyzed. This analysis was complemented through calculation of the interaction energy between the QH2 substrate and bc1 complex residues. The binding of QH2 at the Q-site was further investigated through quantum chemical calculations accounting for the effect of wide range electron polarization.

Quinol Binding Motif at the Q-Site

Based on early crystallographic results,[45] the key residue holding QH2 at the Q-site is histidine H156 of ISP, shown in Figure 1c. On oxidation of the Fe2S2 cluster, this residue undergoes a dramatic change in pK, from ∼12.5 to 7.6, and this allows it to serve as an H-bond acceptor in the bond from QH2, thus stabilizing the reaction complex, and initiating the Q-cycle.[14,46−52] Since the protonation state of H156 is still debated,[13,21] two MD simulations were performed, as described in the Methods section, to address QH2 binding for two suggested[13,21] protonation states of H156.

Figure 3

Analysis of quinol bonding. Time evolution of the key hydrogen bond lengths stabilizing QH2 binding at the Q-sites of the bc1 complex in monomer A (left plots) and B (right plots). Blue lines show the bond lengths calculated for Model I (see Figure 2a), while red lines show the bond lengths for Model II (see Figure 2b). The insets in each panel illustrate the corresponding hydrogen bonding motifs, with lengths labeled d1, d2, and d3. Bond length d1 is defined differently in the case of Models I and II.

Figure 2 shows the two binding motifs of QH2, at the Q-site of the simulated protein–membrane–solvent system. In Model I, shown in Figure 2a, H156 is ϵ-protonated and forms a hydrogen bond with its HE2 hydrogen atom to QH2, in which the histidine is the H-donor; this hydrogen bond is characterized through the distance d1. ϵ-protonated H156 also had been assumed in an earlier MD study.[21] In Model II, shown in Figure 2b, H156 is deprotonated and hydrogen bonded to QH2 through its NE2 nitrogen atom, acting as an H-acceptor. This initial state was also assumed in earlier studies,[53] though with a more truncated quantum chemical model than considered presently.

Figure 2

Quinol binding at the Q-site of the bc1 complex. Shown are binding and coordination of QH2 at the Q-site. Dashed lines represent key hydrogen bonds that coordinate QH2 to residues H156 and Y147. The labels next to these lines indicate the corresponding bond lengths that are shown in Figure 3 and discussed in the text. QH2 binding is primarily coordinated through H156, which is in its ϵ-protonated form in Model I (a), and in its deprotonated form in Model II (b). In the case of Model I, QH2 binding is additionally stabilized through a water molecule.

The simulations show that in Model I a water molecule assumes a stable interaction with QH2 and links the substrate to the C155 and I292 backbone atoms of the ISP and cyt. b subunits, respectively, as illustrated in Figure 2a. The water molecule is only experiencing a stable binding position in the case of the protonated H156 residue; in the case of a deprotonated H156 the water is not present as no stable binding position exists. An earlier MD study[21] demonstrated also stable binding of a water molecule in the case of protonated H156 and proposed that this water molecule is key for proton transfer to the positive side of the membrane.

The QH2 molecule employs both of its hydroxyl groups in binding to the Q-site; one hydroxyl group binds to H156, while the other group forms a hydrogen bond with other residues, e.g., with the Y147 and E295 residues of cyt. b as shown in Figure 2. The binding of QH2 to cyt. b residues has been extensively studied through the effects of mutation on kinetic and thermodynamic properties, but because the second electron transfer from the QH2 molecule is not rate limiting, these approaches are less informative than when applied to the first electron transfer. As a consequence, the chemistry of the second electron transfer from the QH2 is even more controversial.[54] Both QH2 bindings are important for electron and proton transfers occurring at the Q-site of the bc1 complex and are, therefore, considered in the present study.

Analysis of Quinol Bonding at the Q-Site

QH2 binding at the Q-site was analyzed by monitoring the lengths of hydrogen bonds formed by the hydroxyl groups of QH2 with the polar residues of the bc1 complex. Figure 3 shows the time evolution of these bond lengths calculated for the two monomers of the bc1 complex for both Model I (blue lines) and Model II (red lines). The insets to the panels introduce atoms participating in the particular bonding.

The time evolution of distance d1 between the hydroxyl group of QH2 and the H156 residue is shown in the upper panel of Figure 3. In the case of Model I the bond length d1 is short throughout the entire 360 ns long simulation, while in the case of Model II this bond length spontaneously increases during the simulation, thereby indicating that QH2 is more mobile at the Q-site in the latter case. This behavior is observed for both monomers A and B of the bc1 complex.

The time evolution of distance d2, the hydrogen bond length between hydroxyl groups of QH2 and the Y147 residue of cyt. b, is presented in the middle panels of Figure 3. Our simulations reveal that there arises a hydrogen bond between QH2 and the Y147 residue in the case of Model I, the bond being less stable in the case of Model II. Indeed, for Model I, the distance d2 fluctuates around 2.3 Å, with few increases up to ∼5 Å. The simulation of Model II also features formation of this hydrogen bond, which, however, is seen to become broken more frequently. The difference in d2 observed in the two simulations is more profound in monomer A and is very small in monomer B, suggesting, therefore, that spontaneous formation and breaking of the hydrogen bond between the QH2 molecule and the Y147 residues occurs to a similar degree for Models I and II.

In the course of our MD simulations, the Y147 and E295 side chains rotate such that Y147 comes to lie between QH2 and E295 and forms hydrogen bonds with the −OH of the QH2 molecule. Formation of such a hydrogen bond between QH2 and the Y147 residue was unexpected, since this configuration had not been observed in any of the structures. However, because the reaction complex is formed under metastable conditions, a set of forces comes into play that are not explored under crystallographic conditions. Previous studies of QH2 binding through MD simulation suggested involvement of some water molecules instead.[19,21] In the original X-ray structure of the bc1 complex the Y147 residue is located far from the Q-site and only in the course of the present MD simulations is this side chain seen to turn and to form a hydrogen bond with QH2. Interestingly, in the course of MD equilibration of both Models I and II, the Y147 rotation occurs rapidly after releasing the side chains of the bc1 complex proteins, while still keeping backbone atoms constrained; the displacements observed require little movement of the backbone. Formation of the hydrogen bond between QH2 and the Y147 residue has significant implications on the proton transfer path at the Q-site. Based on the H-bond formed with stigmatellin in all structures with this inhibitor, E295 has previously been considered to be the most likely acceptor of the second proton in the Q-site reaction.[19]

Snapshots from Models I and II MD simulations, featuring hydrogen bonding between Y147 and E295, are shown in Figure 2. In the simulations, the QH2 molecule forms a hydrogen bond with the Y147 residue, thereby preventing hydrogen bonding between QH2 and the E295 residue of cyt. b, as suggested earlier.[13] Despite the lack of this permanent hydrogen bond, the E295 residue remains located in close proximity to Y147 and our MD simulations reveal that E295 spontaneously forms hydrogen bonds to Y147. The lower panels of Figure 3 show the time dependence of the distance d3 between the side chains of Y147 and E295. The figure illustrates that for both Models I and II, in either monomer of the bc1 complex, the distance d3 is ∼3.9 Å, occasionally going down to ∼1.8 Å, i.e., to a value typical for a hydrogen bond formed.

Quinol Stabilization through a Water Molecule

An important attribute of QH2 binding to the Q-site in Model I is a water molecule as illustrated in Figure 2a. The water molecule binds to QH2 and keeps it attached to the I292 and C155 backbone oxygen atoms of cyt. b and ISP, respectively. This binding, however, is stable only in the case of Model I, where H156 is protonated.

Stabilization of QH2 binding at the Q-site through a water molecule can be characterized through hydrogen bonds that the latter forms with QH2 and surrounding residues. The time evolution of the corresponding bond lengths is shown in Figure 4. The recorded distances, d, d, and d, in Figure 4 are defined in Figure 2a for Model I. Although not shown in Figure 2b, a water in a similar configuration can bind when H156 is deprotonated (Figure 4b for Model II), but with different H-bonding characteristics.

Figure 4

Quinol binding at the Q-site coordinated by a water molecule. Shown is the evolution of the lengths of hydrogen bonds formed between QH2 and a water molecule trapped within monomer A and monomer B of the bc1 complex. Green lines represent the length d of the bond between the H1 atom of QH2 and the OH2 atom of the trapped water molecule, while red and blue lines correspond to the lenghts d and d of hydrogen bonds formed between the H1 and H2 atoms of the water molecule and the O atoms of C155 and I292, respectively. The hydrogen bonding network along with d, d and d are shown in Figure 2a. A water molecule is bound to the QH2 molecule throughout the entire simulation in the case of Model I (a) and is seen to bind only sporadically in the case of Model II (b).

Figure 4a shows that in the case of Model I, in both monomers of the bc1 complex a water molecule is spontaneously bound to QH2 as all three distances, d, d, and d, fluctuate around 2 Å. Only for monomer A d is seen to fluctuate around 2.5 Å, the other two bonds remaining at lower bond length values. Our simulations show that water molecules at the Q-site interchange on a time scale of 15 ps.

The presence of a water molecule at the Q-site, seen in the present simulations, is consistent with an earlier MD simulation.[21] In that study this water molecule was suggested to be important not only because it stabilizes QH2 binding at the Q-site, but also because it can act as the initial acceptor of a proton from QH2, triggering the Q-cycle. The analysis of the hydrogen bond distances performed earlier[21] is fully consistent with the present investigation, thereby validating the present simulation. The analysis of hydrogen bonds associated with a water molecule in Model I indicates that this water molecule acts as a key that fits the space between the QH2 molecule and the protein’s backbone atoms.

In case of Model II, the water molecule is not fitting as a key; the distances d, d, and d turn out to be significantly different from each other and are not seen to go below 2 Å simultaneously, so that the triple H-bonded configuration of Model I is not seen; the distance d between the QH2 molecule and a nearby water molecule rarely gets below 2.5 Å and can increase up to 8 Å and even higher. Water molecules occasionally form hydrogen bonds with the backbone oxygen of I292, as illustrated in Figure 4b through the distance d. This observation is, however, likely irrelevant for QH2 binding as the distances d and d are ∼4–5 Å in all cases when d goes below 2 Å. The water molecule, thus, is expected to float around the I292 residue, but not to participate in functional chemistry of the Q-site as it does in the case of Model I.

Interaction Energy between Quinol and the bc1 Complex

The interaction energy of the QH2 molecule with the bc1 complex is an important characteristic for QH2 binding at the Q-site. Figure 5 shows the interaction energy of the QH2 head group with the rest of the system for Models I and II and for monomers A and B.

Figure 5

Quinol interaction with bc1 complex. Shown is the time evolution of the interaction energy for the QH2 head group and the rest of the simulated system, including water molecules, lipids, and the bc1 complex proteins. Blue and red lines show the energies calculated for Models I and II, respectively. The energies calculated for each step of the simulation are shown in shaded colors, while intense color shows a gliding average with energies averaged over a gliding window of 100 ps. Vertical arrows indicate the time instances for monomers A and B, when the QH2 molecule unbinds from H156 as seen in the Model II simulation.

The interaction energy of QH2 with the bc1 complex in Models I and II fluctuates by about ±5 kcal/mol (the value expected for a hydrogen bond) around −43 kcal/mol until the QH2 molecule unbinds from the H156 residue. The time instance at which this unbinding occurs is seen in Figure 5 for both monomers of the bc1 complex. QH2 unbinding from H156 is characterized through the distance d1, which spontaneously increases from 2.3 Å up to 4.0 Å in the case of Model II, as illustrated in the upper panels of Figure 3.

The calculated interaction energy of H156 with QH2 differs between Models I and II; Figure 5 shows that the energy in Model I is lower than the energy in Model II (accordingly, the binding is stronger), the difference amounting to 7.8 kcal/mol for monomer A and 4.3 kcal/mol for monomer B. The lower binding affinity of QH2 in the case of Model II supports our earlier conclusion that QH2 is more mobile at the Q-site in this case.

We note finally, that in the case of monomer A of Model I there is a short disruption in QH2 binding to H156 occurring at ∼240–300 ns. This disruption is reflected by an increase of the interaction energy of the QH2 molecule with the bc1 complex and is correlated also with a slight increase of the distance d1 between QH2 and H156, as seen in the upper panels of Figure 3.

Quantum Chemistry Study of Q-Site Quinol Binding

The description of QH2 binding at the Q-site of the bc1 complex by means of classical MD simulations cannot account for polarization of electrons across the binding complex that is expected to contribute to binding strength;[55,56] the participation of electronic degrees of freedom can only be accounted for through quantum chemical calculations. The quantum effects also contribute to the geometry of the binding motif of the QH2 molecule at the Q-site; this contribution was studied for both Models I and II through quantum chemical structure optimization of the Q-site binding complex. This complex included in our calculation a QH2 molecule, the Fe2S2 cluster of ISP, and several surrounding polar residues that are expected to play a role in primary electron and proton transfer of the Q-cycle as highlighted in Figure 6. We note that a similar methodology was successfully applied before to different protein systems.[57−59]

Figure 6

Quantum chemistry model of QH2 binding at the Q-site of the bc1 complex. Included in the quantum chemical description are the components shown here in licorice representation, namely the QH2 head group and all residues within 10 Å from the head group. The coloring of the bc1 complex secondary structure illustrates the charge state of the protein amino acids: negative (red), polar uncharged (green), and hydrophobic uncharged (yellow). Side chain groups of polar and negatively charged amino acids surrounding the QH2 head group (Y147, H135, H156, E295, Y302) were included into the computational model to describe environmental effects on the QH2 binding, while distant charged side chains that point away from the QH2 head group and E295 were not included in the quantum chemical description. The Fe2S2 cluster and all its coordinating amino acids (licorice) have also been included in the computational model.

The present choice of the quantum chemistry models is dictated by the electrostatic potential at the Q-site. Figure 6 shows a characteristic snapshot from the MD simulations of Model I in which the charged and polar amino acid residues can be distinguished from the neutral ones. Figure 6 shows that within a range of about 10 Å from the QH2 head group only the amino acids H156, H135, E295, Y147, and Y302 should have an impact on QH2 binding. Therefore, the side chains of these amino acids were used to construct the quantum mechanical model of the Q-site. The initial structures of the Q-site, employed for the quantum chemical studies, were taken from MD simulations for Models I and II, as shown in Figures 7a and 8a, respectively. The details of the optimization protocol are provided in the Methods section.

Figure 7

Quantum chemical optimization of the binding complex Q-site for Model I. Shown are residues of the bc1 complex directly involved in the binding of QH2 at the Q-site and covered in our quantum chemistry analysis of Model I. The quantum chemical Q-site for Model I includes residues Y147, I292, E295, and Y302 from the cyt. b subunit, as well as the Fe2S2 cluster together with residues C133, C153, C155, H156, and H135 of the ISP subunit. The initial configuration (a) used in the quantum chemical calculations is a conformational average from a bc1 complex before equilibration. The optimized structure, shown in part b, features rearrangements of residues Y147 and E295 accompanied by spontaneous proton transfer from Y147 to E295. The Cα atoms of the amino acid residues (cyan spheres) were fixed during the quantum chemical optimization process.

Figure 8

Quantum chemical optimization of the Q-site binding complex for Model II. Shown are residues of the bc1 complex directly involved in the binding of the QH2 molecule at the Q-site and covered in our quantum chemical analysis of Model II. The Q-site for Model II is constructed similarly to the Q-site for Model I introduced in Figure 7. The key difference here is that residue H156 is deprotonated and that there is no water molecule linking QH2 with residues I292 and C155. As in case of Model I, the initial configuration of the Q-site used for quantum chemical optimization (a) is a conformational average from a bc1 complex before equilibration. The optimized structure, shown in part b, features rearrangements of residues Y147 and E295. The Cα atoms of the amino acids residues (cyan spheres) were fixed during the quantum chemical optimization process.

For the initial geometries of the Q-site for Models I and II that started the quantum chemical optimizations we selected average postequilibration configurations arising in our MD simulations, as done before.[60,61] Figures 7 and 8 show the initial and optimized geometries of the binding complexes at the Q-site for Models I and II.

In the course of the quantum chemical optimization of Model I the hydrogen bond network involving the QH2 molecues, the H156 residue, and the water molecule remains largely intact, as a comparison of parts a and b of Figure 7 shows. Indeed, the quantum chemically optimized distance d1 is 1.78 Å, while the average ⟨d1⟩ from the MD simulation is 2.18 Å. The distances are indicated in Figure 7 and labeled in Figure 2. The distances relevant for QH2 binding at the Q-site are summarized in Table 3. One can see that the distances d, d, and d for Model I, characterizing positioning of the water molecule discussed above, are also consistent between MD simulations and QM calculations.

Table 3

Hydrogen Bond Lengths at the Q-Sitea

	bond distance averages (errors) from MD trajectories		bond distances from QMO
distances (Å)	Model I	Model II	Model I	Model II
d1	2.18 (0.13)	2.31 (0.40), 0–200 ns	1.78	1.96
		3.96 (0.51), 200–360 ns
d2	2.52 (0.53)	3.28 (1.26)	1.59	1.91
d3	4.10 (0.67)	3.82 (1.09)	1.02	1.54
dw1	2.03 (0.33)		1.75
dw1	2.46 (0.36)		1.80
dw1	2.65 (0.97)		1.73

Listed are average hydrogen bond lengths calculated from MD trajectories and from QM optimizations (QMO) for Models I and II. The values in parentheses give the standard deviations of the corresponding distances calculated for the 360 ns trajectories (see Table 2). Due to lack of stability a water molecule binding between QH2 and residues I292 and C155 was not included in the quantum chemical optimization of Model II, and, therefore, the distances d, d, and d are not given in that case. Since the distance d1 in Model II experiences a step-wise change (see Figure 3), the average distance value in that case was calculated for the first 200 ns of the MD trajectory. The average value of d1 for the 200–360 ns interval is 3.96 Å.

The Q-site geometry optimization for Model I as shown in Figure 7 involves significant rearrangements of residues Y147 and E295 of cyt. b. The side chain of residue Y147 turns toward the hydroxyl group of QH2 to form a hydrogen bond with it. This turn is accompanied by rearrangement of residue E295, the side chain of which also interacts strongly with Y147. The latter interaction turns out to be so strong that in the course of the quantum chemical optimization procedure the Y147 residue loses its proton and donates it to E295, as seen in Figure 7b.

The rearrangements involving Y147 and E295 are observed both in the quantum chemical optimization and in the classical MD simulations. Table 3 summarizes the mean values of the studied distances, calculated for the entire MD simulation trajectories, as well as values obtained from the corresponding quantum chemistry optimizations. In the case of Model I the distance d2 becomes as small as 2.34 Å, while it fluctuates around a mean value of 2.52 Å. The distance d3 varies around its mean value of 4.1 Å in Figure 3, decreasing occasionally to 1.79 Å. The spontaneous proton transfer from Y147 to E295, however, cannot arise in an MD simulation, as this process involves breaking of a chemical bond and, therefore, involves electronic degrees of freedom. This explains the discrepancy in the values of distances d2 and d3 between the MD simulations and the QM optimization, as listed in Table 3; the quantum chemically optimized structure of the Q-site reveals that a proton has been transferred from Y147 to E295, while MD simulations show only a rearrangement that makes such a proton transfer possible.

The optimized geometry of the Q-site of Model II does not change significantly in the course of the quantum chemical optimization; however, the side chain of residue E295 undergoes a rearrangement. This rearrangement is illustrated in Figure 8 and results in the formation of a hydrogen bond between residues E295 and Y147. The turn of the side chain of E295 for Model II is similar to an analogous rearrangement in Model I; however, no spontaneous proton transfer occurs in the former case. Nevertheless, the H-bonded network is certainly appropriate for proton transfer by a Grotthus-type mechanism, as the bond length of the hydroxyl group of the Y147 side chain has increased by 0.06 Å during the quantum chemical optimization, while d2 and d3 decreased by 0.09 and 2.08 Å, respectively. Rearrangement of Y147 and E295 suggests a probable QH2 → Y147 → E295 path for proton transfer at the Q-site. Such a path had not been discussed earlier, but the strong interaction between the involved residues and the spontaneous proton transfer between Y147 and E295 suggest it.

The comparison of the hydrogen bond lengths in Table3 reveals that the distance d1 values, calculated for the Q-site through quantum chemical energy optimization and through MD simulations, agree well with each other while the QH2 molecule is bound to H156. The distances d2 and d3 computed from MD simulations show irregular behavior (see middle and lower panels of Figure 3), and, therefore, their average values are somewhat larger than in the case of a single quantum chemical optimization.

Partial charges and spin density distributions were also analyzed for the quantum chemistry models in order to provide a more accurate description of the QH2 binding. For that purpose, both quantum chemistry models were split into 11 fragments as defined in Figure 9. Each atom of the Fe2S2 cluster was defined as a fragment itself as the iron ions of the cluster should have an antiferromagnetic coupling,[37,44] and it was defined so in the calculations. To focus on this coupling, the total difference spin density of the α-electron spins (spin up) and β-electron spins (spin down) was computed. The spin density difference is illustrated by the distributions obtained from the quantum chemistry calculations and shown as symmetric transparent surfaces in red and blue in Figure 9.

Figure 9

Fragments of the Q-site and spin densities. The Q-sites for Models I (a) and II (b) have been subdivided into 11 fragments, whose total charges were analyzed separately and summarized in Table 4. Atoms belonging to a certain fragment are highlighted with the same color. Transparent surfaces around the Fe-atoms show the difference of the total spin density calculated between all α-electrons (spin up) and all β-electrons (spin down) in the system. The surfaces are shown for the contour values 0.01 (blue) and −0.01 (red) and indicate antiferromagnetic coupling of the two iron ions.

Table 4 summarizes the partial charges of all fragments (calculated through the Mulliken and ESP-fitted schemes), and the total spin density; the data were obtained by using the B3LYP/6-311G(d) and B3LYP/6-311+G(d) (values in parentheses) methods. It can be noticed that there are some discrepancies between charges calculated with the two methods, mostly for the iron ions, that can be attributed to the additional diffuse functions which are present in the B3LYP/6-311+G(d) computational method. The antiferromagnetic coupling between Fe and Fe ions, however, is still evidenced in both computational approaches; for both Models I and II as seen in Table 4, the spin densities of the iron ions turn out to be around ±3.8 (B3LYP/6-311G(d)) and ±3.4 (B3LYP/6-311+G(d)).

Table 4

Charges and Spin Densities of Q-Site Fragmentsa

	Mulliken charges		ESP charges		spin densities
fragment	Model I	Model II	Model I	Model II	Model I	Model II
1: FeA	1.34 (−0.17)	1.37 (−0.25)	0.29 (−0.06)	0.37 (−0.01)	3.82 (3.42)	3.85 (3.40)
2: FeB	1.30 (1.42)	1.30 (1.23)	1.01 (0.74)	1.10 (0.77)	–3.78 (−3.41)	–3.78 (−3.36)
3: SA	–1.04 (−0.75)	–1.05 (−0.65)	–0.57 (−0.35)	–0.60 (−0.35)	0.19 (0.27)	0.13 (0.20)
4: SB	–0.91 (−0.58)	–0.95 (−0.57)	–0.57 (−0.38)	–0.65 (−0.43)	0.22 (0.33)	0.13 (0.24)
5: C133	–0.48 (−0.48)	–0.54 (−0.53)	–0.39 (−0.32)	–0.47 (−0.40)	–0.31 (−0.38)	–0.28(−0.33)
6: C153	–0.58 (−0.81)	–0.60 (−0.86)	–0.52 (−0.42)	–0.58 (−0.49)	–0.27 (−0.37)	–0.26(−0.38)
7: H135	0.13 (0.62)	0.11 (0.56)	0.35 (0.37)	0.35 (0.36)	0.07 (0.08)	0.06(0.09)
8: H156 + C155	0.18 (0.75)	–0.54 (0.11)	0.34 (0.39)	–0.47 (−0.39)	0.07 (0.06)	0.14(0.15)
9: QH2	–0.09 (−0.11)	–0.11 (−0.04)	–0.08 (−0.21)	–0.05 (−0.06)	0.00 (0.00)	0.00 (0.00)
10: I292 + E295 + Y302 + H2O + H+	–0.12 (0.01)	–0.41 (−0.37)	–0.04 (0.05)	–0.32 (−0.31)	0.00 (0.00)	0.00 (0.00)
11: Y147–	–0.74 (−0.90)	–0.58 (−0.62)	–0.81 (−0.82)	–0.66 (−0.67)	0.00 (0.00)	0.00 (0.00)

The table summarizes Mulliken charges, ESP-fitted charges, as well as the spin densities of the Q-site fragments calculated with the B3LYP/6-311G(d) and B3LYP/6-311+G(d) methods for Models I and II. The fragments (first column) are defined in Figure 9. Charges and spin densities for each fragment are shown which correspond to the calculations done with the B3LYP/6-311G(d) and B3LYP/6-311+G(d) methods; the values of the B3LYP/6-311+G(d) calculation are given in parentheses. The electronic spin density is defined as the total electron density of electrons of one spin minus the total electron density of the electrons of the opposite spin. The fragments are defined similarly in both models (see Figure 9); however, a water molecule is present in fragment 10 of Model I only, while H156 in fragment 8 is protonated in Model I and deprotonated in Model II. The proton from Y147 (fragment 11) has intentionally been included in fragment 10 to elucidate the charge migration between Y147 and E295.

Although the diffuse functions significantly impact the partial charges of the two iron ions of the Fe2S2 cluster, both DFT methods (B3LYP/6-311G(d) and B3LYP/6-311+G(d)) can still be used to describe QH2 binding, as follows from the analysis of total charges of the QH2 and the ISP fragments only. Table 5 summarizes the total charges of the relevant subsystems for the QH2 binding: the ISP (Fe2S2 cluster and its covalently bonded amino acids), the QH2 head group, and the environment (involving in the present quantum mechanical calculations all other surrounding amino acids). The comparison of results for both B3LYP/6-311G(d) and B3LYP/6-311+G(d) shows that redistribution of charges in the ISP due to the diffuse functions does not affect the total charge of the ISP subsystems and, therefore, there is no charge delocalization between ISP and QH2. Thus, one concludes that deviation of fragment charges seen in Table 4 for the two employed computational methods shows small sensitivity of the Q-site models to the choice of the computational method.

Table 5

Total Charges of the Charge Transfer Subsystems of the Q-Sitea

	Model I		Model II
fragment	6-311G(d)	6-311+G(d)	6-311G(d)	6-311+G(d)
ISP	–0.06 (−0.06)	0.00 (−0.03)	–0.90 (−0.95)	–0.96 (−0.94)
QH2	–0.09 (−0.08)	–0.11 (−0.21)	–0.11 (−0.05)	–0.04 (−0.06)
environment	–0.86 (−0.85)	–0.89 (−0.77)	–0.95 (−0.98)	–0.99 (−0.98)

Total Mulliken and ESP-fitted charges of the quinol QH2, ISP part of the Q-site, and of the remaining system included in the calculations (E292 + E295 + Y302 + Y147). For each fragment the first numbers correspond to the Mulliken charges, while the numbers in parentheses are the ESP-fitted charges. The charges were calculated with the B3LYP/6-311G(d) and B3LYP/6-311+G(d) methods for Models I and II, as indicated.

It is important to note that the calculated partial charges of fragments containing Y147 (fragment 11) and E295 (fragment 10), shown in Table 4, clearly indicate that Y147 loses a proton, as the partial charge of fragment 11 (Y147 without its hydroxyl hydrogen) is highly negative. Therefore, one concludes that such transfer occurring during the quantum chemical optimization (see Figure 7) is not a proton coupled electron transfer between these two amino acids. The proton is likely shifted away from Y147 due to the flatness of the potential energy landscape between Y147 and E295. This idea is supported by the quantum chemical optimization results of Model II, in which no proton shift from Y147 to E295 has been observed, despite the similarity of Models I and II. Since the proton can move rather freely between Y147 and E295, it is natural to expect that its exact localization should largely be irrelevant for QH2 binding to the Q-site.

To study the binding strength of the QH2 substrate at the Q-site, its interaction energy with the rest of the system was computed. For the subsystems, as defined above, the total energy and QH2 binding energy are summarized in Table 6. It is important to note that the QH2 binding in Model I is stronger than in Model II, as also indicated by MD simulations (see Figure 5), and, therefore, in qualitative agreement with the bonding analysis shown in Figure 3. Table 6 also shows that utilization of diffuse functions in the calculations (B3LYP/6-311+G(d) method) consistently lowers the binding energy in both models by approximately 0.3 eV.

Table 6

Total and Interaction Energies of Q-Site Fragmentsa

	energies
fragment	Model I	Model II
total energy (au)	–8224.507921 (−8224.633022)	–8147.459089 (−8147.584573)
QH2 energy (au)	–690.504046 (−690.521193)	–690.506375 (−690.523010)
empty Q0 energy (au)	–7533.923715 (−7534.043413)	–7456.915233 (−7457.032775)
ΔE (eV)	–2.181 (−1.862)	–1.020 (−0.783)

The energies of the Q-site with the quinol, of the individual quinol QH2, and of the empty Q-site were calculated quantum chemically employing the B3LYP/6-311G(d) and B3LYP/6-311+G(d) methods for Models I and II. The energy difference ΔE is the binding energy of the quinol at the Q-site. The B3LYP/6-311+G(d) values are indicated in parentheses.

Discussion and Outlook

The bc1 complex converts, in the photosynthetic apparatus, the energy available from light harvesting to a proton gradient by using work stored on reduction of Q (quinone) to QH2 (quinol). Critical for reaching the high efficiency observed in the energy transformation is a bifurcated electron transfer at the Q-site that sends the first electron of a QH2 substrate to the Fe2S2 center (down in Figure 1a,b) and the second to the Q-site (up in Figure 1a,b), releasing both QH2 protons to the periplasmic space (in the down-direction of the photosynthetic membrane in Figure 1). The Q thus formed at the Q-site is released, a new QH2 bound, and the bifurcated reaction undergoes a second cycle. As a result, two electrons are transferred across the membrane up and two electrons are passed via the Fe2S2 center down to cyt. c2, the latter shuttling the electrons one-by-one back to the reaction center. Simultaneously, four protons are released on one side (the down-side) of the membrane. The electrons transferred to the Q-site reduce a Q bound there to QH2, with uptake of two protons from the up-side of the membrane, in effect recycling one of the two QH2 processed at the Q-site.

Such bifurcated electron transfer along with monodirectional proton transfer requires a reaction complex in which QH2 binds in a geometrically highly structured site, in which the configuration constrains the reaction coordinate so as to enable productive forward chemistry but minimize the nonproductive bypasses. The structural and kinetic detail needed to establish the mechanism underlying initiation of the Q-cycle requires an approach in which experiment and computational modeling complement each other. Because of the central role of electron and proton transfer processes in the Q-cycle computational modeling has to combine classical molecular dynamics and quantum chemical calculations. The suggestion of the Q-cycle was made already 39 years ago,[62] but only today do we have the computational means available to carry out the quantum chemical calculations needed to demonstrate its detailed structure–function characteristics. In this regard the present study opens a critical new chapter in the field of bioenergetics. Clearly, the starting point of the research needs to be establishment of the QH2 bound state formed in the Q binding site of the bc1 complex.

The electrostatic properties of residues forming the Q-site of the bc1 complex provide an optimal environment for the binding of QH2 and are responsible for the success of the further redox processes. Thus, the interaction network between QH2 and some key residues of the bc1 complex determine the QH2 binding and the subsequent initiation of the Q-cycle. The results from the MD simulations and the quantum chemistry optimizations in the present investigation have revealed that two binding motifs of the substrate molecule are feasible, namely Models I and II, which differ in the protonation state of ISP residue H156. Experimental results, including pH dependence of the electronic turnover rate,[46,47,63,64] site specific mutagenesis,[14,47,65,66] studies of the thermodynamic cycle,[14,48,67−69] structural,[49,70−72] and spectroscopic studies,[50−52,73,74] have shown that the protonation state of H156 plays a key role in determining the rate of the initial proton-coupled electron transfer reaction at the Q-site. The maximal rates are observed for deprotonated H156, where the oxidized ISP forms a hydrogen bond with the hydroxyl group of the quinol through the Nϵ atom (Model II in the present paper). However, the electron turnover is also possible at lower rates,[46,47,63,64] at pH values well below pKA, demonstrating that the deprotonated histidine configuration is not essential. In the present study, both protonation states of H156 were considered as they both are supported in the literature.[20,48,74]

The performed analysis of QH2 bonding at the Q-site revealed that the substrate molecule forms a hydrogen bonding network, which bridges residues suitable for charge transfer reactions. In this respect, the central role in accommodating QH2 to the Q-site is played by the H156 residue, as it is involved in the first charge transfer reactions occurring between QH2 and cyt. c1. This residue links the QH2 substrate molecule and the Fe2S2 cluster of the ISP, and, therefore, operates as an acceptor group which carries the flux of electrons in oxidation of QH2 by ISPox to generate the intermediate semiquinone. The critical question to be addressed is whether this path also carries the proton flux (Model II), or whether the two charges separate so as to reach the aqueous phase through separate channels (Model I).

Quinol binding in the two computational models (Models I and II) allowed us to identify different scenarios of the primary electron and proton transfer reactions. Thus, in the case of Model I the simulations indicate that the primary proton transfer happens directly to the positive side of the membrane, through a water channel; it follows from the simulations that a water molecule is constantly bound to one of the QH2 hydroxyl groups through a hydrogen bond, as seen in Figure 2a; the coupled electron and proton transfer reaction in this case iswhich corresponds to the configuration assumed in a previous study.[21]

However, in the case of deprotonated H156 (Model II), the QH2 molecule is expected to donate the proton directly to the H156 residue of the ISP, together with a shift of an electron to the Fe2S2 of the ISP. The reaction describing the charge transfer reaction in this case is simplywhere ISP includes both the Fe2S2 cluster and the H156 residue. As previous investigations suggest,[37,44] the two iron ions of the Fe2S2 cluster exhibit anti ferromagnetic coupling. However, such coupling is not explicitly shown in eqs 1 and 2 (spin states of the Fe2S2 cluster ions not specified), as the respective spin state is largely conserved during the proton-coupled electron transfer.

In the present study we have also identified a key role of residues Y147 and E295 in stabilizing QH2 binding and in providing a possible pathway for the deprotonation of the semiquinone intermediate QH• that sets up the reaction complex from which the second electron transfer occurs. A critical consideration in understanding the bifurcated reaction is that of how bypass reactions are minimized. Undesirable bypasses could, in principle, oxidize the semiquinone intermediate without passing its electron to the low potential chain.[75−77] Most interesting from a medical point of view (the bc1 complex arises also in the respiratory pathway common to most human cells) are reactions leading to reactive oxygen species generated by the semiquinone via reduction of O2 to the superoxide anion, as such species play a central role in many of the diseases associated with aging, arthritis, or heart disease.[78−80]

Identification of a possible role for Y147 is particularly important because such a role opens possibilities for control of the bifurcation that have not previously been explored.[81] An extensive earlier mutagenesis study[82] of Y147 had concluded that, although this residue serves an important function, the tyrosine hydroxyl group was relatively unimportant, because mutation Y → F gave a strain with about half the activity of wild type. In order to reconcile this finding with the role proposed here, we note that Y147 is not serving as the proton acceptor; that function is served by E295 as suggested earlier.[67,83] Instead, Y147 serves as an intermediary between semiquinone and E295. Although E295 plays a critical role in the second electron transfer through exit of the second proton, the suggested function as a direct ligand stabilizing quinol binding[6,67] is no longer supported.[82,84] This function of Y147 as an intermediary might alternatively be fulfilled by a water bridge, as seen in previous simulations.[21,53] The configuration with a water molecule replacing Y147 might allow the reduced flux seen in the Y → F strain; currently, QM calculations are being extended to explore the role of water in quinol binding stabilization.

As recognized in other systems,[84] in proton coupled electron transfer processes, the full gamut can be found, depending on how closely the charges are coupled, and the several different examples of such reactions in the bc1 complex cover a wide range. Although the first electron transfer of the Q-cycle has a special importance as the rate limiting reaction, this feature also means that it is the most accessible to direct kinetic assay. Indeed, the Marcus–Brønsted approach previously developed from kinetic studies[67,82,84] provides a satisfactory model that is compatible with the atomistic picture developed here. The reverse also holds true, namely that the computational approaches pioneered here will be even more essential in studies of less accessible processes.

In contrast to the first step of the Q-cycle, where the electron and proton transfer are likely tightly coupled, in the second step, the charges likely separate so that the electron and proton follow different pathways. The Coulombic consequences provide a rich area for computational investigation. Similarly, many characteristics of the two-electron gate at the Q-site, also proton-coupled electron transfers, remain to be satisfactorily explained. Add to the monomeric processes the complexities of dimeric interactions, and the challenges are exciting. MD can also play a direct role in understanding electron transfer processes.

The two computational models studied here provide insight into the binding of QH2 at the Q-site of the bc1 complex. Whereas the key residue for such QH2 binding, namely H156 of the ISP, is often considered to be in its deprotonated state, our data reveal that a protonated H156 can provide an alternate binding configuration. Because alternative QH2 binding motifs at the Q-site have now been identified, one should investigate the quantum chemistry of the bifurcated electron transfers from both reaction complexes to ascertain what determines the productive path, and should extend these studies to all reaction pathways for electrons and protons in the Q-cycle. Earlier experimental studies[14,49,51,52,65,67,73,74] already provided some support in favor of the proposed QH2 binding regimes. With the aim of testing the proposed pathway for forward redox reaction, electron paramagnetic resonance (EPR) measurements can be carried out, since significant accumulation of QH• at the Q-site is reached when the second quinol oxidation step is impeded.[83,85,86] Through mutational analysis it might be possible to examine structural associations of QH• and ISP(H)• and to determine rate constants through rapid-mix freeze-quench EPR. Furthermore, CW-EPR and pulsed EPR approaches can be utilized to measure interactions between QH• and ISP(H)•[87] if both contribute to the product state. The SQo states recently examined[53,85,880] suggest a downstream product, after dissociation. Identification of the intermediate state would allow detailed study of the suggested electron and proton transfer paths. The detailed atomistic analysis performed here can then definitely guide further investigations.

A more profound understanding of the entire bc1 complex function requires use of highly accurate quantum chemistry methods[88−90] and more extensive MD simulations for identifying possible conformational changes occurring during the Q-cycle. Using the full power of available computational tools in the context of a rich experimental background the inner mechanism of the bc1 complex, ubiquitous in bioenergetics membranes of the photosynthetic and the respiratory type, should hopefully be revealed.