The Q-Cycle Mechanism of the - ACS Publications - American

Feb 28, 2017 - STAMPEDE processors. Simulations involving lipids use the latest CHARMM36 force field in which the problem with consistency of lipid de...
0 downloads 0 Views 7MB Size
Article pubs.acs.org/JPCB

The Q‑Cycle Mechanism of the bc1 Complex: A Biologist’s Perspective on Atomistic Studies Antony R. Crofts,*,†,‡ Stuart W. Rose,‡ Rodney L. Burton,†,⊥ Amit V. Desai,§ Paul J. A. Kenis,§ and Sergei A. Dikanov∥

Downloaded via UNIV OF TOLEDO on September 25, 2018 at 13:21:46 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



Department of Biochemistry, University of Illinois at Urbana−Champaign, 419 Roger Adams Lab, 600 South Mathews Avenue, Urbana, Illinois 61801, United States ‡ Center for Biophysics and Quantitative Biology, University of Illinois at Urbana−Champaign, 179 Loomis, 1110 West Green Street, Urbana, Illinois 61801, United States § Department of Chemical & Biomolecular Engineering, University of Illinois at Urbana−Champaign, 600 South Mathews Avenue, Urbana, Illinois 61801, United States ∥ Department of Veterinary Clinical Medicine, University of Illinois at Urbana−Champaign, 1008 West Hazelwood Drive, Urbana, Illinois 61801, United States S Supporting Information *

ABSTRACT: The Q-cycle mechanism of the bc1 complex is now well enough understood to allow application of advanced computational approaches to the study of atomistic processes. In addition to the main features of the mechanism, these include control and gating of the bifurcated reaction at the Qosite, through which generation of damaging reactive oxygen species is minimized. We report a new molecular dynamics model of the Rhodobacter sphaeroides bc1 complex implemented in a native membrane, and constructed so as to eliminate blemishes apparent in earlier Rhodobacter models. Unconstrained MD simulations after equilibration with ubiquinol and ubiquinone respectively at Qo- and Qi-sites show that substrate binding configurations at both sites are different in important details from earlier models. We also demonstrate a new Qo-site intermediate, formed in the sub-ms time range, in which semiquinone remains complexed with the reduced iron sulfur protein. We discuss this, and a spring-loaded mechanism for modulating interactions of the iron sulfur protein with occupants of the Qo-site, in the context of control and gating roles. Such atomistic features of the mechanism can usefully be explored through simulation, but we stress the importance of constraints from physical chemistry and biology, both in setting up a simulation and in interpreting results.



INTRODUCTION The central role played by enzymes of the bc1 complex family in the energy metabolism of the biosphere has been supported by half a century of progress summarized in the Q-cycle mechanism.1−8 Crystallographic structures show that the complex in both mitochondria and bacteria is a homodimer, with a catalytic core of three subunits in each monomer, cytochrome (cyt) b, cyt c1, and the Rieske iron−sulfur protein (ISP).9−14 Each monomer has four metal centers, hemes bH and bL in cyt b, the 2Fe2S cluster in ISP, and heme c1 in cyt c1. In some bacterial complexes, no other subunits are structurally defined, but in Rhodobacter sphaeroides, the complex expressed under photosynthetic growth has an additional subunit IV of uncertain function.15,16 Mitochondrial complexes have up to eight additional subunits. For some the function is known and for many it is uncertain,9,12,17,18 but since none contain redox centers, or impinge on the catalytic sites, they are unlikely to be involved directly in the mechanism. © 2017 American Chemical Society

The interest of the Q-cycle (Figure 1) is that it doubles the energy saved from each photon compared to that if the ubiquinol (QH2) generated in the photochemistry directly reduced the cytochrome (cyt) c2 oxidized. Instead, QH2 is oxidized in a bifurcated reaction at the ubiquinol oxidizing site (Qo-site) that directs one electron via a high-potential chain, 2Fe2S cluster and heme c1, to cyt c2, leaving an intermediate semiquinone (SQ) at the Qo-site (SQo). This is oxidized by a different, low-potential chain, hemes bL and bH, which delivers the electron to reduce ubiquinone (Q) at the quinone-reducing site (Qi-site). Since this takes two electrons, the Qo-site reaction must turn over twice, with the electron from the first stored as SQi, which is then reduced to QH2 on the second turnover. Special Issue: Klaus Schulten Memorial Issue Received: October 18, 2016 Revised: February 26, 2017 Published: February 28, 2017 3701

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

Figure 1. Consensus version of the Q-cycle. Left, the redox centers of a monomer of the Rb. sphaeroides bc1 complex (taken from the MD). The reactions of the Q-cycle (equations to the right) are shown by magenta arrows (H-transfers), green arrows (Q, QH2 exchange at catalytic sites), cyan arrows (electron transfers), and red arrows (H+ transfers). QH2 bound at the Qo-site shows the ES complex of the bifurcated reaction. The mobile head of ISP provides the binding domain for the 2Fe2S cluster (space-filled lower left), the site of redox exchange. This moves from the position shown here (where the electron is transferred from QH2 at the Qo-site to ISPox) to H-bond with heme c1 (where the ISPH is oxidized), an ∼30 Å rotational displacement.

much of the current interest lies in how the mechanism has evolved to minimize such processes. The second new set of information relating to gating mechanisms comes from development in collaboration with the Kenis lab of a rapid mix/freeze quench apparatus based on a microfluidic mixer. The approach has allowed us to characterize an intermediate product of the QH2 oxidation reaction, proposed to involve SQo in complex with the reduced ISP (ISPH•), a novel state with interesting properties.24 These include spectroscopic features, the sub-ms kinetics, thermodynamic parameters, and incorporation into the kinetic scheme. High resolution EPR approaches should allow us to establish molecular parameters for the spin coupling that can be directly compared to outcomes from quantum chemical (QC) calculation. In this paper, we review recent progress on understanding partial processes of the bc1 complex mechanism at the atomistic level, discuss what remains to be achieved, and focus attention on those aspects of importance from the perspective of a biologist. In several recent studies, the interest has been in the application of state-of-the-art computational approaches and the biology has sometimes taken second place. These studies have arrived at important conclusions, but these might be suspect if the starting computational models have defects. The main impetus in developing a new model comes from this obvious starting point; modeling of mechanism through simulation is necessarily conjectural, and success in relating the outcome to the experimental base through which computation can be validated must depend on how realistic the model is. The Role of Spring Loading in Control. The mutations in C. elegans above have focused attention on the role of ISP and movement of its extrinsic cluster binding domain in control and gating of ROS production. From previous work in bacterial systems,25−29 the primary mutation in isp-1(qm150) (Pro-225 → Ser, in the cluster binding domain, equivalent to P146 (bovine), P166 (yeast), P136 (Rb. sphaeroides)) likely introduces distortion of structure at the surface of the cluster domain that interfaces with cyt b. This would result in a steric hindrance to impede formation of the enzyme−substrate (ES) complex, and hence limit ROS production. Then, the suppressor mutants (eight different point mutations in the tether) would relieve the inhibition. The fact that all of the

This QH2 is recycled to replace one of those oxidized at the Qosite. Since protons are released on one side and taken up on the other, and electrons cross the insulating phase through the lowpotential chain (Figure 1), the Q-cycle pumps 2H+/QH2 oxidized across the membrane, with the following overall reaction: QH 2 + 2ferricyt c2 + 2H+ N ⇌ Q + 2ferrocyt c 2 + 4H+P (1)

Despite its “textbook” status, many features of the Q-cycle model remain controversial. The sites at which ubiquinone is reduced and oxidized involve four different proton-coupled electron transfer (PCET) reactions, in which the different properties of the semiquinone intermediate confer very different patterns of reaction. These are exploited both in the reaction mechanism and in control and gating functions. Our recent kinetic model for the reaction at the Qo-site in the antimycin-inhibited complex includes kinetic and thermodynamic parameters for 14 partial processes, most of them determined directly using conventional physicochemical approaches.4,6 The model takes account of existing thermodynamic constraints but is already under revision to include control processes revealed in recent work in two related areas. First, collaboration with the Kaeberlein lab on understanding the longevity in C. elegans following mutation of ISP (strain isp1(qm150)), and suppressor strains in the same subunit, has suggested a spring-loaded mechanism.19,20 To understand why this is interesting, note in Figure 1 the large distance between the 2Fe2S cluster, to which the electron from QH2 is transferred, and heme c1, to which the electron is next delivered. Structures show that the ISP extrinsic domain containing the cluster can rotate on a tethering span to close the distance.14,21,22 The spring-loading involves the counteraction of forces associated with formation of the enzyme− substrate (ES) complex, which pull on the tether span to extend a helical segment, and the forces trying to remake the helix. These involve additional partial processes, and their thermodynamic parameters need to be accommodated by redistribution of work terms among existing equilibrium constants. Longevity relates to a darker side of the enzymeits production of reactive oxygen species (ROS)which makes it an important player in the aging process.2,19,20,23 The SQo intermediate can reduce O2 to generate superoxide, the precursor of ROS, and 3702

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

the ES complex, ISPox is not paramagnetic, so this “handle” cannot be used, but binding information can be determined from similar displacements through kinetic characterization.3,4,31,32 In considering the role of spring-loading, the critical point of interest is that the binding to form the ES complex involves similar counteracting forces, but on electron transfer to the intermediate products, ISPH• and SQo, these must dissociate for further progress. The force stored in the tether on extension of the “spring” when the ES complex is formed is available to drive the dissociation of ISPH• from SQo. This would favor the forward chemistry, minimize SQo occupancy, and thus help to explain protection against ROS production. However, the new SQo.ISPH• complex introduced above does not dissociate on a ms time scale, as would be necessary for rapid forward chemistry. The intermediate state is trapped, a feature of obvious interest to be discussed further below. Examples of Some Problems in Computational Simulation of the bc1 Complex. Most molecular dynamics (MD) simulations start from the crystallographic structures. The atoms are frozen in a lattice constrained by contacts with neighboring components of the unit cell, and information on dynamic features is limited to the mobility shown by B-factors. For membrane proteins, the prison is even more unnatural. There is no membrane, but instead, ancillary lipids, detergent molecules, waters, etc., fill the interstices. The MD simulation is set up to liberate the native structure from this prison, an essential preliminary to mechanistic exploration. Crystallographic models of the bc1 complex from vertebrate mitochondria or bacteria show a dimeric structure with a substantial empty volume in the dimer interface to the N-side of the closely packed protein interface between the bL hemes. It is unlikely that this void represents a vacuum. In support of this, in higher resolution structures of the yeast mitochondrial complex, electron density in the site has been resolved, and specific phospholipids identified.33 One cardiolipin molecule occupies the central cavity, and tails from other lipids contribute more peripherally.33,34 The volume is defined by a protein scaffolding, symmetrical about the vertical axis of the dimer, consisting of transmembrane helices, and transverse (amphipathic) helices at the level of phospholipid head groups. The scaffolding restricts access at this level but allows access from the hydrophobic lipid phase. The former restriction might be expected to impede ready diffusion of phospholipid from the membrane into the volume. The scaffolding helices screen the two Qi-site volumes from the “void”, but substrate access is from a contiguous volume so that the quinone tails mingle with the lipids in the void. These features can be seen in our model, with the “void” volume now partly occupied by two phosphatidylglycerol (PG) molecules (Figure 3, Figure SI-1). In two previous MD simulations of Rhodobacter complexes accessible to the authors, it was assumed that the problem of the void would be addressed by “the physics”, as explored in MD simulation. In practice, the MD eliminated the void but in one case by partial collapse of the protein. Significantly, during 350 ns of simulation, although lipid tails explored and partly filled the volume, no phospholipid molecule diffused in. In another case, the void was filled by flooding with waters. Neither of these physical solutions is likely to be natural. In the former case, a partial unfolding from the crystallographic configuration of one of the transverse scaffolding helices occurred, which also disrupted the volume around one Qi-site of the complex, precluding realistic simulations involving that

simulation suppressor mutations were in the same subunit as the primary lesion but were well-separated from that site in the sequence and in the structure. This localizes the underlying mechanisms for both longevity and its suppression to this subunit and the involvement of these specific regions. When the extrinsic head of ISP docks on cyt b, the driving force is associated with binding (predominantly, the H-bond between the occupant and His-152) and the protein interfaces involved. The work from binding contributes one set of the forces in the spring-loaded scenario. The other set of counteracting forces is associated with a change from helical to elongated chain in the tether span where the suppressor mutations are located; the binding force pulls on the tether to extend it. The loss of binding energy for formation of the ES complex from the primary mutation is likely overcome in the suppressor mutations by weakening the helical forces in the tether (Figure 2). Values of direct relevance to the spring-loaded scenar-

Figure 2. The ISP subunit from structures of mitochondrial bc1 complexes showing conformations that change with the occupancy of the Qo-site (adapted from refs 21, 22, 30, and 83). The structures show ISP with the fully extended tether (left) on binding at the Qo-site (here with stigmatellin) and with the fully relaxed tether (right) when no bond is formed with a Qo-site occupant. In the spring-loaded mechanism,19,20,25 the experimentally determined binding free energy (∼ −6 kJ/mol) is the difference between the work involved in binding (−23 kJ/mol) and the work needed to extend the tether (−17 kJ/ mol), referred to the relaxed state (ΔG° = 0 for both). The values shown are estimated from work on ISP mutants in Rb. capsulatus and Rb. sphaeroides.25−27,29 The curved arrow shows the first electron transfer reaction after formation of the ES complex. When SQo.ISPH• separates by spatial displacements to allow interaction with acceptors, the binding force is lost, and the energy stored in the extended spring is available to drive displacement.

io19,20,25 have been estimated for both forces from studies of Rhodobacter mutants, though these were originally constructed to test the role of the tether in movement.26−29 Complexes of reduced ISP (ISPH•) are formed with several different inhibitors, or with quinone (Q). In binding of ISPH• with different occupants, the paramagnetic property allows characterization of the interaction through changes in the EPR spectrum. Displacement of values for Em, ISP when ISPH• is pulled out of the equilibrium mix allows estimation of binding free energies,21,22,30 which depend on occupant and vary with mutation. Simple thermodynamic considerations then provide values for changes in the counteracting forces, as detailed in the Supporting Information for ref 19, summarized in Figure 2. In 3703

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

Figure 3. Cross section through the MD model of the Rb. sphaeroides bc1 complex after 31 ns of a production run. (A) The initial occupants were replaced by parametrized molecules, with ubiquinone (Qi) in the Qi-site and ubiquinol (QoH2) in the Qo-site. A slice through the protein, shown as a cartoon, reveals the prosthetic groups colored by chain, embedded in the membrane between aqueous phases, with lipids and waters represented by lines for the bonds. The redox centers are shown by VDW spheres (for hemes and 2Fe2S-cluster) or by licorice bonds, colored as below. The “void” is the ∇-shaped space defined by a scaffold of membrane spanning and transverse amphipathic helices (center, top of protein), here occupied by lipid. (Stereo pair for crossed-eye viewing.) (B) The protein stripped away to show the redox centers, to facilitate identity: all redox groups except 2Fe2S are shown by licorice bonds, the hemes are in CPK colors; QH2 at the Qo-site is yellow; Q at the Qi-site is cyan; the 2Fe2S cluster is shown by VDW spheres. (C) The same view of the protein but zoomed to highlight two lipids (PG, shown by van der Waals spheres) occupying the void, with the b-type hemes (licorice bonds, CPK colors) for reference, and the scaffolding helices, showing how exchange of phospholipids would be impeded at the headgroup level. Structure taken from the trajectory exploring formation of the ES complex, at a frame of ∼31 ns, when the bond to H152 had stabilized.

the former. The problem of eliminating the void has been averted in simulations from the Róg group by introduction of a cardiolipin molecule to occupy the void,35,36 effectively simulating in a Rhodobacter capsulatus bc1 complex the yeast configuration. A second set of problems lies in representation of the native membrane. Early MD simulations of the mitochondrial bc1 complex had used a simple 1-palmitoyl-2-oleoyl-sn-glycero-3-

volume. In the latter, the waters would have introduced a high dielectric phase in a volume lined by hydrophobic residues suitable for lipid interactions. Because of the proximity of the Qi-sites, occupation by waters would significantly change the physical chemistry of the site. However, the focus in both papers was the Qo-site reaction, and since this is on the other side of the protein from the Qi-site, it was supposed that disruption of structure at the latter would have little effect on 3704

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B phosphocholine (POPC) membrane,37 still preferred in a recent effort using a Rhodobacter complex.38 A more natural membrane model was introduced by Postila et al.,39 with a composition based on that of mitochondrial membranes, also adopted in our own recent simulation.8 Unfortunately, in equilibrating the membrane model, forces came into play that converted the cis fatty acid side chains to the unnatural trans configuration, thereby substantially altering the membrane properties. Whether this was important or not to simulation of function is not clear. The Rhodobacter are versatile bacteria, like many others, and can adapt their membrane composition to cope with environmental stress. For example, under phosphatelimited growth, much of the phospholipid component of the Rb. sphaeroides membrane (though not cardiolipin or PG) was substituted by nonphosphorus glycolipids.40 Strains can be engineered by mutation to eliminate synthesis of cardiolipin from PG (CL− strain41), and these conditions have been combined to explore the dependence of growth, expression of cytochromes, and activity of respiratory and photosynthetic chains on membrane composition. The tested parameters were not attenuated, even when the CL− strain was grown under phosphate limiting conditions.40,41 The only phospholipid present under the latter conditions was PG, from which it is obvious that no other common phospholipid could be required. Although it is not obvious that modifications in lipid content would alter the protein behavior, a natural membrane is obviously preferable. Of special interest is the high ubiquinone content of the native membrane.42−44 Any natural simulation should at least include this component, even though the second order processes involved in substrate or product binding are not readily accessible on the MD time scale. It seems likely that the void is accounted for by disorder in lipids naturally incorporated on assembly of the complex which were not detected by X-ray diffraction. Indeed, as noted above, a more stable structure has been achieved in Rb. capsulatus models by populating the void with a cardiolipin,35,36 as seen in the yeast complex. These and other possibilities are readily explored by simulation. The Limitations of Conventional Approaches to Atomistic Understanding of Mechanism. In a complex mechanism like the Q-cycle, it is widely recognized that the parameters determined experimentally (rate and equilibrium constants, driving forces, reorganization energy, etc.) reflect interplay of many “invisible” partial processes at the atomistic level, even when a reaction can be well-defined by a chemical equation. For these processes, direct measurements are not available, so feasible parameters must be estimated within constraints from conventional measurements.4,6 We have pioneered exploration of some of these atomistic features more quantitatively through previous collaborations with the Schulten group, applying MD and QC approaches to several different mechanistic problems using MD models based on mitochondrial or Rhodobacter bc1 complexes.4,8,37 Building on this experience, we here introduce a new MD model using the Rb. sphaeroides bc1 complex, in which the protein environment has been modeled in a native membrane, including ubiquinone, and are currently running unconstrained MD simulations (production runs) to explore the mechanism at the atomistic level under XSEDE support, discussed below.

available MD program, with demonstrated scalability on many platforms. Currently, the molecular systems presented here can be modeled for 10−15 ns/day when run at 2 ns/step on 400 STAMPEDE processors. Simulations involving lipids use the latest CHARMM36 force field in which the problem with consistency of lipid density has been resolved.46 The simulations with proteins and prosthetic groups use the CHARMM36 force field with CMAP corrections,47 supplemented by custom built topologies from collaborators. Water molecules are represented explicitly by the TIP3P model.48 In all simulations, the temperature was maintained constant at 310 K using Langevin dynamics with a damping coefficient of 1 ps−1 and the pressure at 1 atm using the Langevin Nosé−Hoover method.49,50 Long-range electrostatic forces will be calculated without truncating using the particle mesh Ewald (PME) method.51 For all simulations, 2 fs will be used as the integration time step. Umbrella sampling (US)52 with replica exchange molecular dynamics53 will be used in some projects to calculate the potentials of mean force (PMF) associated with the studied processes. The weighted histogram analysis method (WHAM)54−56 will be used to construct the PMF from the calculated trajectory files. All of these methods are standard and fully implemented in NAMD. Collaborations will allow us to extend the range to include application of more specialized approaches. Rapid Mix/Freeze Quench Protocol. To poise bc1 complex with heme bH reduced and bL oxidized, 11 μM bc1 complex in a buffer of 50 mM 3-(N-morpholino)propanesulfonic acid (MOPS) at pH 8.5, 100 mM KCl, with 0.1% dodecylmaltoside (DDM) and 10% glycerol (buffer A) was incubated in the presence of 65 μM 6-(10-hydroxydecyl)2,3-dimethoxy-5-methyl-1,4-benzoquinone (OH-decylUQ, idebenone), 70 nM NDH-2 (NADH.ubiquinone reductase from E. coli), and saturating antimycin A. The mixture was flushed with argon for 2 h in a Thunberg cuvette before mixing with 0.5 mM NADH held in the bulb, and then allowed to incubate for 30 min. The redox state of the b-hemes and NADH was monitored spectrophotometrically. Heme bH became reduced with a lag on reduction of the Q-pool (t1/2 ∼ 6 min), but somewhat remarkably, heme bL remained oxidized (at least out to 90 min), despite the strongly reducing conditions provided by the excess NADH (Eh ∼ −350 mV). In a separate anaerobic chamber, 800 μM equine ferricyt c in buffer A was flushed with argon to keep it anaerobic. Other mixtures could be used to explore different starting states. The two mixtures were transferred to the reaction syringes under strictly anaerobic conditions before mixing to start the reaction. The reaction progress was assayed using a rapid-mix, freezequenching approach.57−59 An ultrafast microfluidic mixer was used to rapidly mix the reactants and eject the mixed solution in the form of a jet onto supercooled rotating copper wheels to quench the reaction within ∼30 μs. The Z-shaped mixing channel was designed to enable sub-millisecond mixing,60 and the operating conditions were optimized (on the basis of the physics of fluid mechanics and surface tension) to eject the mixed solution in a jet as opposed to droplets. A myoglobin/ sodium azide reaction was used to calibrate our apparatus as in Appleyard et al.,61 calibrated to an instrumental dead time of ∼30 μs, similar to the 50 μs observed in ref 60.

METHODS Classical Molecular Dynamics Simulations. MD simulations are performed using NAMD,45 a highly parallel, publicly

RESULTS AND DISCUSSION The Bifurcated Reaction at the Qo-Site. The finer features of control and gating cannot be appreciated without





3705

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

for the ES complex. However, this involves QH2 and the oxidized ISP (ISPox), for which the polarity would demand that the ring >COH acted as a donor to the H-bond to Nε of the dissociated H152 of ISPox 22,67,68 (pKox1 ∼ 7.6). This configuration provides an obvious pathway for both proton and electron transfer, through the H-bond, to histidine and the cluster, respectively. A second clue then came from the redox chemistry of ISP, which showed that, over the physiological range, the Em value varied with pH, indicating that reduction involved both an electron and a proton.7,68−71 We suggested that the changes in pKox1 in the physiological range seen on redox titration of mutant strains, and reflected in the pH dependence of electron transfer,72,73 could be well explained as resulting from the changes in pKs of the two histidine ligands (H131 and H152 in Rb. sphaeroides) to Fe2 of the cluster.1,69,71 This hypothesis has now been validated by direct measurement by NMR spectroscopy on the Thermus thermophilus ISP.74 A third clue came from work in the Nocera group75 on the rate of electron transfer between H-bonded partners in a model system, in which the groups involved (a diamine on one side Hbonded to a carboxylate on the other) could be exchanged without much effect on the photochemically generated driving force. Surprisingly, the observed rate in aprotic solutions was strongly dependent on the donor−acceptor configuration of the H-bond through which the electron and proton had to pass. They recognized that the proton-coupled electron transfer depended on the difference between the pK values of the groups forming the H-bond, which defined a Brønsted term specifying the probability that the H+ would be favorably situated for electron transfer. By combining these hints, we could propose not only a simple model for the ES complex but also a simple explanation for the pH dependence of electron transfer in the physiological range, and the simple explanation for why the reaction was so slow, summarized in eq 2 and discussed extensively elsewhere.1,32,69,73 Previous Computational Studies. In the recent collaborative work8 with the Schulten and Solov’yov groups, we have used molecular dynamics, supplemented by snapshots of critical states treated through QC calculation to examine the ES complex of the bifurcated reaction in the bc1 complex from Rb. capsulatus. Modeling by MD simulation suggested a novel configuration, consistent with previously unexplained outcomes of mutagenesis studies of the key residue, Tyr-147,76 and in marked contrast to the direct participation of E295 previously expected.77 In subsequent studies,78−80 polar residues contributing to binding of QH2 and the 2Fe2S cluster were approximated through density function theory (DFT) to provide a QM model in which QC calculations were used to explore this state. In extending these calculations to an investigation of the first electron transfer, they used an approach based on the intrinsic reaction coordinate to explore potential pathways to product states, starting with different protonation states of H156 (equivalent to H152 in Rb. sphaeroides).78 Although configurations with His-156 either protonated or not can be represented by competent QC models, only the complex in which ISPox has His-156 initially dissociated at Nε provided a productive forward chemistry.78 Several other recent papers8,38,39,77 have included QC calculations of the first step of the bifurcated reaction at the Qo-site. This is rate limiting under substrate saturation, and therefore of special mechanistic interest. Critical to the reaction profile is the configuration of the ES complex; conclusions about the PCET reactions involved depended on the model

some deeper insights into the Q-cycle mechanism through which the complex operates (Figure 1). Our current picture for the Qo-site reaction involves mobile dancers (ISP extrinsic domain, SQo, several side chain rotations) in a molecular ballet of reactions linked to the chemistry, driven by overall work terms well characterized by conventional approaches. However, the ballet is choreographed by atomistic forces. Recent evidence suggests that, at the Qo-site, control and gating of the first electron transfer involves the spring-loaded mechanism mentioned above,19,20 while that of the second electron transfer involves Coulombic couplings that come into play as H+ and electron separate on oxidation of SQo (initially formed as the neutral SQ (QH•) by heme bL). Such Coulombic effects are complex, especially so when electrostatic forces from the neighboring protein and solvent perturb electronic orbitals, and vice versa. Additional complexity is introduced by equilibration between dimers. Kinetic data suggest that electron transfer across the dimer interface (between the bL hemes) does not contribute a significant fraction of the normal forward flux,3,62,63 but certainly, Coulombic effects and slow electronic equilibration can occur.3,64 To explore such systems requires application of several different areas of physicochemical expertise, and the field has recently attracted such diverse interest. The Rate Limiting Process. The first electron transfer at the Qo-site is rate-limiting overall,1,4,32 with the rate constant determined by the unique properties of the proton-coupled electron transfer mechanism, which can be summarized in the following Marcus−Brønsted form: log10 k = 13 −

° + λET)2 (ΔG ET β (R − 3.6) − γ λET 2.303

− (pK QH2 − pK app)

(2)

Here β/2.303 is the slope of dependence of log10 k on electron transfer distance65 with β = 1.4 Å−1, γ is F/(4 × 2.303RT) = 4.23 in this treatment,66 and λET is the Marcus reorganization energy for the electron transfer component. The terms to the right can be seen as log10 probabilities whittling away the limiting rate constant (k ∼ 1013 s−1 at bond vibrational rates). The observed rate constant is k ∼ 103 s−1, but models based on structures suggested that the distance for electron transfer is ∼7 Å, for which the expected value would be very much faster (∼108 s−1 using typical values for reorganization energy, λ ∼ 0.7 V). The rate constant is lowered from the value expected for simple electron transfer (the Marcus value, the first three terms on the right in eq 2) by the improbability of finding the proton in the favorable position along the H-bond through which both electrons and protons are transferred (the rightmost Brønsted term). This explanation for the rate constant discrepancy came from clues provided by structures. Although none of them have shown a quinone species in the catalytic site, structures with some types of inhibitor (stigmatellin, 5-n-undecyl-6-hydroxy4,7-dioxobenzothiazole (UHDBT)) all showed specific Hbonding from a ring >CO (or >CO−) to ISPH•, involving one of the histidine ligands to the 2Fe2S cluster (His-152 in Rb. sphaeroides). With pKred ∼ 12.5, this side chain is protonated at physiological pH, and would serve as a H-bond donor. Since the inhibitors were thought of as quinone analogues, it was natural to suggest a similar pattern for binding of ubiquinone, involving ISPH• and one of the redox active >CO groups of the quinone ring in the gx = 1.800 complex.21 Extrapolating from this to the quinol form then provided an obvious model 3706

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

other strains when additional mutants were included, suggesting that Glu-295 had little if any effect on binding of QH2.6 In addition, at least one structure (PDB ID 1sqv) has been published in which the equivalent tyrosine (Y131 in beef) H-bonds to the inhibitor UHDBT at the Qo-site.10 This arrangement has been discussed in detail by Berry et al.9 They were also able to identify such a configuration in electron densities for some of their own UHDBT-occupied structures. These results are of obvious importance in supporting the new configuration. QC Simulations of the First Electron Transfer. Both the ES complex and its transition through the first electron transfer have recently been simulated in MD/QC calculations in a continuation of our collaboration with the Schulten and Solov’yov laboratories.8,78 The outcomes of these calculations show a gratifying agreement with the mechanism suggested above but with a new wealth of atomistic detail. Free-energy profiles of the reaction coordinate from ref 78 are similar to that expected from eq 2, including energy levels, and endergonic nature. In the physicochemical profile (Figure 4), the

chosen. In all models, the ES complex involved QH2 in the site, H-bonded to the same histidine of ISPox docked at its interface with cyt b. In Martin et al.,38 the electron transfer was treated in a somewhat truncated QC model, with parallel studies of the dielectric response of the protein and environment. In determination of the latter, the reaction was considered as starting with QH−.ISPH+•, unencumbered by the Brønsted term, and treated as if the proton transfer was separate, so that the electron was transferred with full Coulombic consequence. However, the reaction is intrinsically neutral: QH2.ISPox ⇌ QH• + ISPH•. In the PCET mechanism above,1,7,32,69,71,81 with the ES complex formed with QH2 H-bonded with Nε of His152 of ISPox dissociated (pKox1 ∼ 7.6 in the isolated protein70), formation of the complex (and hence the rate of reaction) depends on pH, as shown by the strong correlation between change in pKapp for electron transfer rate and change in pKox1 in ISP mutant strains.7,8 Since the ⩾C−OH of QH2 (with pK ∼ 11.5) donates the H atom to the H-bond, electrons and protons would follow the same path on oxidation by ISPox. Most of dielectric response needed would then be accommodated by electronic redistribution in the cluster and its ligands (as seen in the pK change of the histidines on reduction7,69−71,82). In the model considered by Postila et al.,39 the electron transfer also occurred though the H-bond to histidine (H156 in Rb. capsulatus, H152 in Rb. sphaeroides) but with the histidine Nε initially protonated, so that the histidine was not able to serve as a proton acceptor. Release of the proton therefore had to involve the other quinol ⩾C−OH, and a pathway, separate from that of the electron, via a bound water molecule, to which they assigned a key role. The MD model used by Barragan et al.8 was initially parametrized following that of Postila et al.39 but then was adapted so as to allow a testing of the Crofts et al. model,1,7,32,69,71,81 with the histidine Nε initially dissociated. This simulation revealed unexpected features of the ES complex configuration, already noted above. Previous guesses based on the binding of stigmatellin had assumed that Glu-295 would provide a H-bond to the ⩾C−OH of the quinol distal from ISP. Since E295 mutations had a profound effect on function, a critical role had been assumed as acceptor for the second H+, released on oxidation of the intermediate SQo. This was reinforced by the striking rotational displacement of E295, seen in structures with myxothiazol (or E-β-methoxyacrylate (MOA) type inhibitors), that both opened a volume for diffusion of SQo and brought the carboxylate group into contact with a water chain to the aqueous exterior, thus providing a pathway for proton release.21,83 Barragan et al.8 found instead that Tyr-147 spontaneously rotated from the position in the stigmatellin structure to serve the H-bonding function to QH2. In their new model, Y147 also H-bonded to E295, and this allowed it to serve a relay function with the glutamate, so that the latter could still act as a proton acceptor. Reinterpretation of results from earlier literature on mutations at this residue,76 and our recent detailed analysis6 of additional mutations at Glu295, seemed to provide support for this novel configuration. In Y147F,76 in which the H-bonding potential was lost, the apparent Km determined from the dependence of the rate of QH2 oxidation on the degree of reduction of the Q pool indicated a substantially reduced affinity for the Qo-site, consistent with a binding function in the native strain. On the other hand, the small increase in Km on mutation of E295 to D or G observed in earlier studies,83 confirmed in later work, was found to be spread to either side of the wild-type value in

Figure 4. Marcus−Brønsted reaction profile of the Qo-site reaction (updated from refs 1 and 7). Values for ΔG are calculated from experimental data. The vertical dotted arrows show limits on estimated values. The reaction follows the path expected from the Marcus− Brønsted treatment. The profile is similar to that generated by QC simulation in Barragan et al.,78 set up with the ES complex (ISP··· QH2) as the starting state, in which QH2 is the H+ donor to H152 Nε in the H-bond, so that the reaction would be expected to follow the same path as that from the Marcus−Brønsted approach.

contribution of the Brønsted barrier is shown as a distinct state, but in the equivalent linear form of eq 2, the last two terms would represent the product of probabilities, and, with the rapid rate constants for H+ transfer in H-bonds,84 would therefore behave as a composite term. The approach through a QM superposition of states simulates not only the profile but also this composite feature, in the landscape from QC calculations.78 In an interesting investigation of the diffusion of O2 into the Qo-site, and its reduction to generate super oxide, Husen and Solov’yov80 have also extended the model from ref 8 to approximate the product state in which the first electron transfer has generated Q•−.ISPH•. In another extension, Barragan et al.78 have explored the QC trajectory between ES 3707

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

(occupied in the initial crystallographic model (PDB ID 2qjy) respectively by stigmatellin and ubiquinone), when occupied by UQ, respectively, in reduced or oxidized form at Qo- and Qisites (Figure 3) (this redox poise would be appropriate to the pool initially 90% oxidized). The stable configuration of the Qisite confirms that in ref 35, in which features of the Qi-site mechanism were simulated in a Rb. capsulatus model of the bc1 complex with cardiolipin in the void. However, our initial unconstrained production run, exploring a continuous trajectory of ∼125 ns, revealed significant differences from those previously reported in configuration of reaction complexes at both sites (Figures 3 and Figure SI-1). The basis of these differences obviously needs to be resolved. In the Barragan et al.78 complex leading to productive forward chemistry, three H-bonds stabilized the quinol: from H156 Nε to QH2 −OH (O2), from Y147 −OH to QH2 −OH (O5), and from Y147 −OH to E295 −COO−. In our simulations, these three residues were all found to participate in H-bond pairings of the Barragan et al. model, suggesting that drastic revision of previous work might not be necessary. Figure 5 shows (in part A) the state of the Qo-site in frame 153 of the trajectory (at ∼30 ns), in which important H-bonds are highlighted, and (in part B) the distances for the H-bonds above, read from the trajectory as it evolves. In the time course shown, ISPox and QH2 are initially separate, but over the first 90 frames (18 ns), the H-bond from QH2 −OH to ISPox Nε of H152 forms, dissociates, and then reforms and stabilizes the ES complex. During the entire trajectory, Y147, E295, and N279 explore configurations in which Y147 rarely visits QH2. Mostly, E295 is busy swapping its association between the other two residues. However, no configuration in which all three H-bonds were simultaneously engaged (which formed the basis of their QC calculations) has yet been reached in our simulation. This volume of the protein also includes several exchangeable waters which are involved in H-bonding with the polar residues, and connecting to the heme bL propionates and Arg-94, likely providing H+ conducting pathways, including one to the Pphase water. Over the remaining time captured in this trajectory, Y147 occasionally (about only 10% of the time; see time course in part B) H-bonds with QH2, as shown in part A. In the frame captured in part A, the two primary H-bonds to QH2 stabilizing the ES complex (from H152 and Y147) are both present, but E295 is distant. In a continuation of this run, with the sampling interval decreased to 1 fs, the pattern persists at longer times (currently to ∼125 ns). If this pattern is confirmed in additional runs, it would require at least a modest change in interpretation of the previous result, that the release of a H+ from QH• involves collisional exchange with Y147, E295, and N279 via stochastic H-bonding, rather than the direct relay previously suggested.78 The Barragan et al. model78 included Y147 and E295 in the QC calculation, and they found that, without E295 included, the forward chemistry was thwarted. This may capture an important feature of the mechanism. The same role could be implemented in the stochastic mechanism above but, unfortunately, not in a selfcontained QC calculation. On the positive side, a smaller model, omitting Y147 and E295, would be more tractable in QC calculations. More difficult to implement would be QC evolution along the reaction profile using snapshots, interrupted by MD simulation allowing equilibration with the stochastic possibilities from structural dynamics. More extensive revision could be justified; it is possible that the earlier modeling21,22,30,83 of the ES complex as involving

complex and this product state. We will discuss these models later in the context of new insights from the SQo.ISPH• state introduced above. Two additional recent papers have considered the role of cardiolipin in stabilizing the structure,36 and in facilitating delivery of protons to the Qi -site mechanism,35 also discussed later. A New MD Model Using the Rb. sphaeroides bc1 Complex. Biochemical approaches can validate (or otherwise) quantum calculations, and insights from quantum chemistry can guide mutational and spectroscopic efforts. For the new model discussed above, the bc1 complex from Rb. sphaeroides has been embedded in a native membrane with the physiological complement of lipids, including ubiquinone, to fit the computational model better to the experimental reality. In designing this model, we have taken note of the two problems introduced above: (1) the void in crystallographic structures (Figure 3) and (2) the non-native membrane used in previous work. Since interpretations of the outcome in these previous simulations might have been distorted by these architectural blemishes, an important preliminary task has been to test alternative scenarios for dealing with the void. We have stabilized the complex by insertion of two phospholipid molecules (PG in Figure 3c), simulating what might be expected on population of the void from the membrane on assembly of the complex. The Róg group has shown in their simulation of the Rb. capsulatus complex35,36 a similar improvement in stability when cardiolipin occupies the void. We chose two PG instead of the cardiolipin because function of the bc1 complex in Rb. sphaeroides does not require cardiolipin (see below). These are important improvements, because they minimize structural displacements of protein or water in response to the void. No population of the hydrophobic volume by waters occurred in the model containing the two lipids. As already noted, although the POPC membranes of earlier models have been replaced by more natural ones, for example, the composition of the mitochondrial membrane,8,39 in setting up these models, the membrane was relaxed by procedures that allowed the natural cis fatty acid side chains to isomerize to the unnatural trans form. In our present model, we have used the published composition (including ubiquinone) of the Rb. sphaeroides chromatophore membrane determined under conditions of growth similar to our own,85 to implement a more natural simulation environment, and taken care to avoid artificial isomerization. Preliminary simulations after relaxation of the model with all constraints removed (short runs ∼5 ns) showed that all three phases (protein, lipid, and water) behaved naturally. In longer trajectories, free diffusion of lipids was observed in the bulk membrane; diffusion of phospholipids was limited to 2D within the constraints of the polar heads, but ubiquinone and ubiquinol also explored the lipid phase through 3D excursion, except when bound at the catalytic sites. The two PGs placed in the initially void volume were retained for at least 125 ns. The behavior in the membrane will allow direct estimation of diffusion coefficients over the time course of a trajectory (to be discussed elsewhere). Water chains previously observed crystallographically were populated by waters, but these exchanged rapidly with the bulk. No waters were found in the “void” volume, except at the headgroup/water interface. These results provide some confidence that modifications to correct previous defects have been effectively implemented. Both the overall structure and layout of the Qo- and Qi-sites retained configurations close to those seen in the starting model 3708

DOI: 10.1021/acs.jpcb.6b10524 J. Phys. Chem. B 2017, 121, 3701−3717

Article

The Journal of Physical Chemistry B

two models indistinguishable experimentally from the earlier model. Possible scenarios are discussed further below. One consideration in deciding between these scenarios is the need to constrain the SQ in the Qo-site. Preference might depend on a simple physical principle; the low probability (high energy cost) of solvating a charged species favors mechanisms involving Q•− as the liberated form, rather than QH•, since the former would have a much lower probability of escape into the lipid phase via the hydrophobic entrance channel. In parallel, we have also run a preliminary simulation with antimycin occupying the Qi-site, a commonly used experimental situation, since it allows ready measurement of turnover of the Qo-site through reduction of heme bH. The measured rate of QH2 oxidation (in the first turnover) when heme bH is initially oxidized is the same in the absence or presence of antimycin.5,86 In line with this, the configuration of the ES complex examined followed essentially the same pattern as that in the absence of inhibitor (not shown). Control and Gating in the Oxidation of SQo. The need for control/gating of the reaction was first noted explicitly by Oscyzka et al.,87 who pointed out that, if rate constants for electron transfer were determined by distance, similar values would apply to oxidation of SQo in forward chemistry, and bypass reactions involving reduction of SQo by heme bL, if both reactions occurred from the location at which it was generated from QH2. Consequently, the bypass reactions, and those associated with ROS production, could occur at rates comparable to the normal forward reaction. To account for the observed slow rates, the reaction would have to be gated. Quite how this happens is under debate. Although turnover is constrained by the limiting process, it is widely recognized that the second electron transfer must occur with a rate constant intrinsically much more rapid than the first. In normal forward chemistry, intermediate product states of the Qo-site reaction must have a low occupancy (