Optimal efficiency of the Q-cycle mechanism around physiological temperatures from an open quantum systems approach

The Q-cycle mechanism entering the electron and proton transport chain in oxygenic photosynthesis is an example of how biological processes can be efficiently investigated with elementary microscopic models. Here we address the problem of energy transport across the cellular membrane from an open quantum system theoretical perspective. We model the cytochrome \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b}_{6}\,f$$\end{document}b6f protein complex under cyclic electron flow conditions starting from a simplified kinetic model, which is hereby revisited in terms of a Markovian quantum master equation formulation and spin-boson Hamiltonian treatment. We apply this model to theoretically demonstrate an optimal thermodynamic efficiency of the Q-cycle around ambient and physiologically relevant temperature conditions. Furthermore, we determine the quantum yield of this complex biochemical process after setting the electrochemical potentials to values well established in the literature. The present work suggests that the theory of quantum open systems can successfully push forward our theoretical understanding of complex biological systems working close to the quantum/classical boundary.

Approaching the field of biological sciences from the perspective of microscopic physical processes is extremely intriguing and yet partly unexplored in current research. Presently, photosynthesis occupies a particularly privileged position in this quest: the molecular basis of its fundamental mechanisms are actively investigated at the biological, biochemical, and biophysical levels. Additionally, the nature of the process itself is well suited to be investigated with the tools of either classical or quantum statistical physics. In fact, it involves analyzing the light-matter interaction within the cell as well as excitation and particle transport through molecular networks, together with compelling energetic considerations. It has been more than ten years now since the attention of quantum physicists was triggered by remarkable spectroscopic results reporting possible evidence for quantum coherent processes in light harvesting complexes 1 . We have henceforth witnessed the rapid development of the field of quantum biology 2,3 , a growth that continues nowadays towards many directions at the forefront of pure and applied scientific research 4,5 .
While the significance and interpretation of many initial findings in quantum biology, most notably the origin and the effect of quantum coherences, is still actively debated 6-10 , one of the most relevant theoretical achievements in the field lies in the application of the theory of open quantum systems to model biological processes featuring a single or a few excitations embedded in the complex landscape represented by a protein, a membrane, or even a cell 3 . Following a similar approach, we hereby focus on the electron transport chain, which represents the second stage of the photosynthetic machinery, immediately connected to the primary light-absorbing processes. Several early applications of open quantum systems (OQS) theory to such crucial particle and charge transport processes at the molecular level are known in the literature [11][12][13][14][15][16][17][18] , including investigations of coherence 11 , correlations 15 and non-Markovian effects 14 . Building on an OQS-based kinetic model originally proposed by Smirnov and Nori 17 , here we develop a simple but accurate Markovian theory of the Q-cycle. The latter is a remarkable biochemical mechanism that efficiently conserves energy stored in excited electrons to increase the number of protons per electron that are pumped against the concentration gradient across the cellular membrane.
In deriving the relevant rate equations describing the dynamics of this complex biochemical process, we exploit the Lindblad formalism to describe single-particle transfer reactions on a microscopic level. Consistently with previous derivations 17 , our approach connects charge transport within the Q-cycle to the reservoir dynamics of the surrounding environment. After benchmarking our alternative model by calculating the so-called quantum yield (i.e., the number of transferred protons per electron during each cycle) and comparing our results to the original ones 17 , we then extend the original study and provide new predictions of biological relevance: in particular, we determine the thermodynamic efficiency of the Q-cycle as a function of the external temperature, and find optimal performances around ambient temperatures compatible with physiological conditions for living cells. We further strengthen our analysis by introducing a new figure of merit combining information from the quantum yield and the transferred charge, which again shows optimality around room temperatures.
We emphasize that even though no genuine quantum effect such as entanglement or coherence-enhanced transport is expected to play any relevant role in the Q-cycle itself, our work confirms that the standard tools of OQS theory and quantum statistical mechanics prove to be very well suited for the physical study of microscopic biochemical transport processes, giving adequate control over the relevant parameters and leading to very natural results with potential biological and experimental relevance. Complemented by alternative mathematical approaches 19,20 , this theoretical work potentially provides a deeper and more complete understanding of many bioenergetic aspects at the molecular level.

Model and Theoretical Description
In Fig. 1a we show a concise diagram of the full photosynthetic electron transport chain. Light is absorbed and stored in the form of chemical energy inside large proteins called photosystem I (PSI) and photosystem II (PSII), respectively 21 . These photosystems are typically assisted by light harvesting complexes (LHC) during the photon absorption and exciton migration stage. The two photosystems are actually successive stages of the same global energy transfer chain in which the energy contained in high energy chemical bonds fuels the transfer of protons ( + H ) against their concentration gradient across the membrane of the photosynthetic cellular organelles, the "thylakoids". This leads to a non-equilibrium distribution of protons, where the proton concentrations on the P-side ("lumen") is orders of magnitude larger than on the N-side ("stroma"). Such out-of-equilibrium unbalance results mainly from two processes: first, the electrons radiatively excited in PSII are extracted from water, leading to the release of oxygen and protons in the lumen; second, the chemical energy of excited electrons is used to transfer protons against their concentration gradient, a process orchestrated by a third membrane-spanning molecular complex, the so-called cytochrome b f 6 complex, which will be the object of our analysis in the following of this work.
The cytochrome b f 6 complex works essentially as an electron-fueled proton pump implementing the Q-cycle mechanism, as originally proposed by P. Mitchell 22,23 and successively discussed and adapted [24][25][26][27] . An OQS-based kinetic model of the Q-cycle mechanism was originally introduced in ref. 17 . The full details of such implementation are reported in the Methods section, for completeness. We now present our original formulation, which is built on the same elementary ingredients but makes explicit use of the Lindblad form of the master equation. By introducing the evolution of the total density matrix of the system, we can describe consistently the microscopic origin of all the contributions to the dynamics, including the Markovian interaction with external reservoirs. The model describes the proton and electron motions through different binding sites, which are individually treated as two-level systems. The proposed use of quantum theory allows for a straightforward derivation of the time evolution equations, directly arising from the definition of all the possible elementary states of the system and their mutual interactions. At the same time, thanks to the intrinsically open system character of our formulation, the interaction between the microscopic details of single-particle transfer reactions can be very naturally and operatively connected to the mesoscopic and reservoir dynamics of the surrounding degrees of freedom. We make use of the pseudospin s = 1/2 formulation, in which the occupied Fock state 1 indicates the presence of a particle on the specific site. Throughout the paper, quantum mechanical operators for electron binding sites will be written with small letters (e.g. annihilation operators A i ), while we will use capital letters for protons (e.g. annihilation operators A i ). Annihilation operators obey anticommutation rules ij . For simplicity of notation, we will assume  = = k 1 B in the following. In Fig. 1b,c we provide a schematic description of the working principles of the cytochrome b f 6 complex under the so called cyclic electron flow (CEF) conditions, a situation that turns out to be simpler to model and to simulate, while keeping all the interesting features of the Q-cycle mechanism [28][29][30][31] . Under CEF, electrons transferred on the P-side to PSI are reinjected into the b f 6 complex from the N-side, thus effectively moving in a closed loop contrary to the conventional Linear Electron Flow (LEF) picture, in which electrons are eventually loaded on a NADPH molecule. Several detailed structural descriptions of the cytochrome b f 6 complex exist in the literature [32][33][34][35][36] , including all the large heme and other prosthetic groups embedded inside the protein scaffold acting as electron binding sites. In particular, it is currently accepted that two reaction sites are present at the two opposite sides of this membrane-spanning protein complex, respectively called Q o towards the interior of the thylakoid (lumen) and Q i on the other side (stroma), where specific mobile electron carriers, called plastoquinone/plastoquinol (PQ/PQH 2 ) can fit one per site at a time. Each PQ molecule (purple in Fig. 1c) can be either discharged or charged with up to two electrons (mainly coming from PSII) and two protons, thus transforming into a PQH 2 . In principle, a whole pool of PQ/PQH 2 molecules is present inside the membrane, to which these hydrophobic species are confined while being free to diffuse, but only few of them can be present at the same time inside the body of the b f 6 complex. For simplicity, a single plastoquinone/plastoquinol molecule is included in our model, described as a carrier of up to two protons and two electrons, which can spatially commute between the N-and P-side of the cellular membrane and alternatively interacting with the Q i and Q o reaction sites. We will often refer to this component as the PQ-shuttle. The Hamiltonian describing the microscopic degrees of freedom of such plastoquinone contains multiple contributions, namely e p (details on the individual terms and a summary of the relevant parameters can be found in the Methods section): the first two terms describe two pairs of independent binding sites for electrons and protons, respectively, while the other terms describe the mutual electrostatic interactions between the charged particles simultaneously present on the shuttle. These are intended as an effective way to take into account the differences in free energy (i.e. in standard redox potential) between different molecular species such as quinone, semiquinone, and quinol.
Under CEF conditions, each PQ reduction inside the b f 6 complex is effectively balanced by a PQH 2 oxidation, contrary to the 1:2 ratio between reductions and oxidations which is assumed in the conventional linear picture of the Q-cycle. When a PQH 2 enters the Q o reaction site, it interacts with two electron binding sites, namely hemeb L and the Iron-Sulfur (ISP) domain (B in Fig. 1c). One electron leaves the PQH 2 and is transferred through the ISP via cytochrome f to a Plastocyanin (Pc) molecule, a water-soluble single electron carrier acting as a connector with PSI. Hence, the original redox energy of such electron is completely consumed. One proton is released to the lumen side (also called P-side for proton rich), leaving the original plastoquinol in a semiquinone state. The second electron is transferred to heme b L and then to heme b H across the membrane, while the second proton is released to the P-side: the plastoquinol is now in the fully oxidized plastoquinone state. On the other side of the membrane, such plastoquinone molecules can bind to the Q i site close to heme b H . The electron that traversed the Pictorial representation of the model described in this article. Single electron binding sites are represented as small black circles, and proton ones as small white circles. The PQ purple element represents the plastoquinone/ plastoquinol molecule, which is able to commute between N-and P-side (stromal and lumenal, respectively) reaction sites, as indicated by purple arrows. Ferredoxin (Fd) and Plastocyanin (Pc) pools on the stromal and lumenal sides, namely collections of water-soluble single electron carriers, are modeled as fermionic reservoirs. Protons in the bulk aqueous phase on the P-and N-side are described in a similar way. Site A models the electron re-injection site from the Fd pool to PQ, while B collectively represents the ISP-cytochrome f chain that transfers electrons to the Pc pool. We assume that PSI transfers electrons from the Plastocyanin pool to the Ferredoxin pool with the help of external energy coming from light absorption. Sites L and H, reported in blue consistently with the color code in panel (b), represent the heme b L and b H groups respectively. The black arrows show the ideal path of the electrons in the system. www.nature.com/scientificreports www.nature.com/scientificreports/ L-H chain reduces this plastoquinone to a semiquinone state. A second fresh electron is provided directly by PSI through a water soluble single-electron carrier called Ferredoxin (Fd) which can bind to the N-side of the b f 6 complex. Two protons are taken up from the stromal side, resulting in a fully reduced PQH 2 , ready to diffuse back into the membrane and towards the P-side. The hemes b L and b H that constitute an electron-recycling chain are again described by a composite Hamiltonian containing free energy terms, H LH free , , accompanied by a Coulomb repulsive interaction, H LH int , , which reflects the redox anticooperativity of the two hemes, as measured in molecular complexes structurally and functionally similar to the b f 6 37 . Referring to Fig. 1c, the B site collectively represents the ISP-cytochrome f chain that transfers electrons to the Pc pool. On the other hand, site A models the electron re-injection site from the Fd pool to PQ. Each of the A and B sites can bind a single electron, while the Ferredoxin and Plastocyanin pools on the stromal and lumenal sides, respectively, are modeled as collections of fermionic oscillators. We assume that PSI transfers electrons from the Plastocyanin pool to the Ferredoxin pool with the help of external energy (i.e., from light absorption), and that this mechanism is at a steady state, such that pools can be described with time-independent parameters. The protons in the bulk aqueous phase on the P-and N-side can be treated in a similar way. With a circuital analogy, the Fd and Pc pools act as source and drain leads for the proton pump, while the PQ/PQH 2 molecules are used as mobile parts. It is easy to see that, as a result of the set of reactions described above, the overall process translocates two protons to the P-side per electron passed to the Pc pool. Therefore, under ideal conditions, the quantum yield (QY) of the reaction, defined as the ratio between the number of protons released on the P-side and the number of electrons whose redox energy is consumed approaches the theoretical value This is particularly relevant in view of calculating the thermodynamic efficiency of the proton translocation process, defined as the ratio between energetic outputs and inputs, which is proportional to QY 17 : where µ ∆ denotes the change in electrochemical potential for electrons and protons. We now solve for the dynamics of this model by explicitly resorting to the theory of OQS applied to the full density matrix of the complete system, ρ. Notice that the total dimension of the Hilbert space corresponds to a collection of 8 independent two-level systems, i.e. = = d 2 256 8 . We hereby describe A and B sites interacting with the electron source and drain in terms of typical Lindblad terms (with { , }, and the average occupation numbers of the reservoirs are described by Fermi-Dirac distributions at the A and B energies  j where T is the temperature and µ i are the source and drain electrochemical potentials, respectively. The electron transfer reactions can be phenomenologically modeled through incoherent Lindblad-type terms including forward and backward Marcus rates 39,40 . A complication arises here, since the single Marcus rates from two individual states explicitly depend on their energy difference. Therefore, we shall distinguish between the cases in which other charged molecules are present or not, for example in all the transitions involving shuttle electrons. This is precisely the effect of the Coulomb interaction terms, which contribute to make some specific transitions more or less favorable. We thus introduce projectors on the states of the system, = P i i i , = … i 1, , 256, and add the following contributions to the master equation i j x y xy j x y i  to describe the transition from state i to state j via the tunneling connection from site x to site y. For simplicity, we assume that the reorganization energies only depend on the sites and not on the states of the system: we can therefore write where ω i is the energy eigenvalue of state i as an eigenstate of the free Hamiltonian of the system, and ∆ xy is a tunneling matrix element. Finally, a similar problem also occurs in defining the interaction between the proton sites on the PQ-shuttle and the N-and P-side reservoirs. In this case, it is the Fermi-Dirac distribution describing the average occupation number that depends on the energy difference of each specific transition. The solution comes again with the help of state projectors: ( where α = N,P , = k 1, 2 and Ω ij is the absolute value of the energy difference between the states i and j. If we now define we can write the full evolution of the system in the form of a single master equation Fd Pc x y Elec if sites x and y are not directly connected in the model. As already anticipated, from Eq. (9) it clearly appears that the off-diagonal entries of the density matrix evolve independently from the populations. The effect of quantum coherence is thus dynamically ineffective here, and it can safely be neglected when considering the electron transfer dynamics. Since we are dealing with a fully incoherent picture, we can then recast the relevant part of Eq. (9) in the form of a Pauli master equation 38 : indeed, the dynamics of the diagonal elements of the density matrix, namely the electron and proton populations in the possible basis states, evolve according to a relaxation equation of the form Relax where P is a vector with 256 components describing the populations (diagonal elements of ρ), and the time-independent relaxation matrix, Λ x ( )

Relax
, receives contributions from all the components of Eq. (9). Indeed, the latter can in general be recast in the form 38 ab ij abij ij from which we see that, in our particular case, Λ =R ( )

Relax ai
aaii . The explicit position dependence of the relaxation matrix comes from the fact that, as we discuss below, coupling strengths and tunneling rates to and from the PQ shuttle depend on its current location inside the membrane. Notice that, at difference with the original approach 17 , the present formalism is based on the complete density matrix of the system and retains the linear structure typical of quantum mechanics, even at the stage where it is reduced to a set of rate equations. Therefore, it fully captures all possible correlation effects (in principle, both at the quantum and semiclassical level) without further assumptions. In particular, it conserves the total number of particles in electron and proton transfer reactions. Moreover, since the full master equation gives in principle access to the dynamics of the quantum mechanical coherences, the formalism could be easily extended to include a richer phenomenology when applied to different contexts.
A very peculiar element of the description of the b f 6 complex under CEF conditions is the stochastic modeling of the PQ shuttle motion inside the membrane between the N-and the P-side, as already introduced 17 . Here we keep such an approach, thus including a stochastic differential equation (SDE) to describe the position of the PQ shuttle, whose diffusive commuting inside the rather dense lipidic membrane is modeled as an overdamped Langevin-Brownian motion where ζ is the drag coefficient (i.e. the inverse of the diffusion coefficient times the temperature, ζ = T /  ) and ξ represents a gaussian noise source with zero average and correlation function ξ ξ The potential U w is added in order to confine the shuttle inside the membrane, while U ch prevents an electrically charged shuttle from crossing the lipid core of the membrane. The tunneling rates between the electron sites on the shuttle and the sites A, B, L, and H will then depend on the position of the shuttle. This can be described with an exponential decay of the coupling rate 21 with the lower sign for = i L,B and the upper sign for = i H,A . Here = x 2 0 nm, and l e is a characteristic length. Since the shuttle can only interact with the aqueous phase of the N-or P-side when it is spatially close to the membrane border, we assume position-dependent rates www.nature.com/scientificreports www.nature.com/scientificreports/ where the upper sign refers to = i N and the lower to = i P, respectively, thus distinguishing the two sides of a membrane with a total thickness of 4 nm. Finally, the distribution of charged aminoacid residues is at the origin of a structural asymmetry in the electrical surface potential of the b f 6 complex, which is taken into account by adding an internal potential between the membrane boundaries ±x 0 . Such a potential is linearly varying between two temperature-dependent extreme values = . V 4 6 N T and = . V 5 4 P T, and its effect is to rescale the free excitation energies of the electron and proton binding sites (see the Methods section), as already described 17 .

Results
In this section, we present some of the most significant results we obtained by performing extensive numerical simulations of the hybrid quantum and stochastic model discussed above (see also details in the Methods section). Notice that all the simulations presented in this work can be performed on a standard Personal Computer. The dynamical behavior of the system was investigated by varying some physically meaningful external parameter, such as the electrochemical gradient of protons across the membrane, the static surface potential, or the temperature. For each parameter under individual scan, the relevant range of variation was empirically identified and divided into a number of steps. For each one of these steps, i.e. for each individual value of the parameter under scan, six independent stochastic diffusion trajectories were simulated, all with the same set of parameters and initial conditions. Each trajectory covered µ 100 s of b f 6 complex activity. The mean value and standard deviation of the quantum yield and the total number of transferred electrons and protons (i.e. the translocation rate) were then computed for each sampling point from the outputs of the six stochastic realizations.
A few benchmark simulations were performed in order to validate our formalism with respect to previous results 17 . In all the simulations, any parameter not involved in the scan was kept fixed at a default value chosen consistently with the relevant literature 17 at a fixed temperature the scan over µ ∆ is equivalent to a scan over ∆pH, namely the concentration gradient of + H . The results for the simulated QY are shown in Fig. 2a. As it can easily be seen, the system is stable and keeps good performances over a wide range of values, with QY well above 1 and close to the ideal value, = QY 2. The quantum yield decreases when the concentration gradient against which the protons are pumped becomes very large. These results show good qualitative agreement with previously reported data 17 , which have been extracted from Figs 3 and 4 of the original reference and plotted on the same scale as ours for the reader's convenience. Notice that our Markovian method predicts slightly higher values of the QY on average, while keeping full consistency for the response of the system under study over the whole range of protonic electrochemical potential. As a second benchmark simulation, the electrostatic surface potential difference ∆ = + V V V P N was varied at constant µ ∆ and T . At room temperature ( = T 298 . Following Smirnov and Nori, we varied ∆V while keeping − V V P N constant. Our results are presented in Fig. 2b. As detailed in the previous section, technical differences in the mathematical modeling and the solution of the quantum stochastic equations, such as the use of a Lindblad approach and a particle-conserving description based on a single global density matrix, inevitably lead to some quantitative discrepancies. Nevertheless, here we report again good qualitative agreement with the original data 17 , also explicitly shown in Fig. 2b. It is worth noticing that our results remain safely below the ideal limiting value = QY 2 even at the extremes of the parameter range.
We now go beyond previous studies by analyzing the effects of temperature variations. Temperature is directly included in the definition of many quantities affecting the dynamical behavior of the model, and influences both electron and proton transfer reactions as well as the stochastic motion of the PQ shuttle and other structural parameters. In our simulations, the temperature was scanned according to three different schemes, summarized in the following.
In scheme I, the temperature was varied while keeping fixed µ ∆ for protons and the surface potential. In such a case, temperature variations only affect the reservoir populations given by Fermi-Dirac distributions and the Marcus electron-transfer rates. Moreover, the diffusion coefficient for the stochastic motion of the PQ shuttle is also affected. According to Eq. (15), the pH gradient also depends on T at constant µ ∆ . Within this formulation, we analyzed the response of the system mainly for what intrinsically concerns the electron and proton transfer dynamics, all the other conditions and properties being unchanged. The results for this scan are presented in Fig. 3a.
In scheme II, we held the pH gradient fixed while changing the temperature. This time, the µ ∆ for protons changes at different sampling points, again according to Eq. (15). This situation is of great biological interest, since realizes the case in which the pH is under external control. At each point, given the temperature and ∆ = . pH 2 5, µ ∆ is obtained from Eq. (15) and then distributed as µ µ = ∆ P /2 and µ µ = −∆ N /2. The results of this scan, as shown in Fig. 3b, are not radically different from the previous ones in Fig. 3a: we can thus infer that the impact of a temperature change on the particle dynamics inside the system dominates over other factors, such as the µ ∆ change. As the third and more complete stage (scheme III) we also took into account the temperature dependence of the surface electrostatic potential. As already shown 17 Fig. 5a for the data of simulation scheme III. The maximum of Q is found in a temperature range that broadly lies around physiological conditions. Finally, we can apply the formal general definition for thermodynamic efficiency, as given in Eq. (2), to the CEF case if we recognize that the input energy is provided by the electrochemical potential gradient of electrons from source (Fd pool) to drain (Pc pool), and that work is performed to move protons from the N-to the P-side. We thus have It is particularly interesting to compute such quantity when both the quantum yield and the electron-to-proton energy conversion ratio vary simultaneously. This is precisely the situation, for example, of the T scan at constant pH and T -dependent surface potential, whose efficiency is plotted in Fig. 5b. Also in such a case where a purely energetic figure of merit is considered, we find a peak of the efficiency around the physiologically relevant temperatures of living organisms.

Discussion
The Q-cycle mechanism is a relevant energy conversion process entering the electron transport chain in oxygenic photosynthesis, but its complex biochemical nature, difficult to capture with theoretical approaches, has to date limited the development of analytic and numerical models. Building on a microscopic description based on an open quantum system formalism, here we have started to bridge the gap between a first principles dynamical evolution of elementary charges and physiologically relevant conclusions on a macroscopic level. First, we were able to reproduce one of the most interesting effects associated with the Q-cycle, namely the bifurcation of electron flows: upon net transfer of two protons from one side (stroma) to the other (lumen) of the cellular membrane, only one of the two electrons simultaneously present on the charge carrier gets consumed, while the other goes through a recycling path within the b f 6 complex ready to trigger another cycle. Moreover, as an original result of this work we were able to capture a few biologically relevant conclusions regarding the efficiency of the mechanism itself and the physiological optimality of its design. In particular, we found that few biochemically relevant figures of merit show maximal values around ambient temperatures for typical values of the electrochemical potentials.
In general terms, it is remarkable that the Q-cycle mechanism is highly conserved, since not only it is found in the thylakoid membranes of photosynthetic organisms, but also in the inner mitochondrial matrix, where it serves as an energy-conserving proton-pumping mechanism in the bc 1 complex as part of the electron transport chain of oxidative phosphorylation. This points not only at the common evolutionary origin of these complexes and the Q-cycle mechanism, but also indicates that the energy conservation that they implement is highly beneficial to a wide variety of organisms. Our results further suggest that such systems might also be examined under the light of selective environmental adaptation. www.nature.com/scientificreports www.nature.com/scientificreports/ In perspective, the present model is intended as a first step towards a full description of the b f 6 complex, and may serve as a starting point to investigate either the detailed pathway of CEF or its behavior in response to specific external conditions. Indeed, while it has been shown that the redox chemistry alone is sufficient to justify the electron bifurcation 19 , some phenomenological observations, mainly obtained while studying the response of the mitochondrial bc 1 complex to the introduction of specific inhibitors, pointed towards the necessity of explicitly formulating some gating mechanisms [42][43][44] . In principle, these hypotheses could be theoretically tested by extending the present model to include an external control over the allowed transitions, and better detailing the individual chemical mechanisms of such processes.
On the technical side, we highlight that our work clearly demonstrates that a Markovian open quantum systems formulation is sufficient to explain the Q-cycle as a naturally emerging process, given the general structural and physical parameters defining the protein complex environment. Moreover, such formalism is certainly suited for the exploration of more genuinely quantum effects at the mesoscale. Indeed, it is worth reminding that quantum mechanical features are naturally embedded into the description at a fundamental level: this is reflected in the linear structure of the resulting set of rate equations, which in principle preserve all possible quantum and semiclassical correlation effects and could readily be extended to situations in which the dynamics of quantum coherences plays a non-trivial role. Given the close similarity of the working principles of the b f 6 complex with other molecular structures, e.g., in the mitochondrial respiratory chain and, more in general, with a large family of particle and energy transfer processes [11][12][13][14][15][16][17][18] , it is certainly worth pointing out that the Lindblad approach could easily be applied to many other biologically relevant examples. On one hand, the conceptual simplicity and modularity of the resulting model could guide the design of experimental investigations in which, e.g., the stoichiometry of the electron and proton transport can be assessed and compared with predictions 45 , keeping all the relevant parameters and external conditions under control. On the other hand, the markovianity and weak-coupling assumptions, which are naturally embedded into the Lindblad form of the master equation, could be violated in some cases, e.g. when the interplay between particle or energy transfer events and the evolution of the protein scaffold is stronger 14 , or when the relaxation time of the environment is comparable to the relevant quantum dynamics 46 : under such conditions, our proposed formalism may serve as a guide towards non-trivial extensions, e.g., to non-Markovian quantum master equations 47 or hierarchical approaches [46][47][48] . Finally, we cannot neglect that the interplay between biology and its formal analysis is bidirectional: in this respect, the elementary CEF formulation adopted here might inspire new possible routes to artificially engineer coupled electron-proton translocations.

Methods
Here we provide further details on the mathematical structure of the model described in the main text, together with the techniques that were used to derive and solve the resulting set of equations.
Structure of the microscopic model. Here we explicitly report all the Hamiltonian terms entering the microscopic model of the Q-cycle discussed in the main text and further details concerning the stochastic motion of the PQ shuttle. As previously stated, these structural elements are chosen consistently with the original kinetic model in ref. 17    while other terms describe the electrostatic interaction between the charged particles on the shuttle. Each site interacts independently with the others and the interaction is diagonal on the Fock basis, thus affecting only the energies of the states without inducing any transition:

 
Equivalent expressions can be given for the the P-and N-side proton reservoirs. The asymmetric electrostatic surface potential, varying linearly with position inside the membrane, can be modeled as  while for the PQ-shuttle this means that the free energies become position-dependent We assume that the surface potential does not affect the Brownian motion: indeed, it could only act on electrically non-neutral states of the shuttle and, in such case, the motion across the membrane is strongly suppressed by the action of the hydrophobic energy barrier U ch . In other words, the stochastic random force originating from molecular collision is dominant in governing the diffusion of the PQ shuttle. The explicit analytic form of the confining potentials for the diffusion are the following: Open quantum systems treatment. The master equation for the quantum evolution was derived by using techniques from the theory of open quantum systems. The source (Fd) and drain (Pc) terms are realizations of the well known case of a single two-level system in contact with a thermal reservoir whose creation and annihilation operators obey anti-commutation rules. The proton reservoir terms are similar in spirit but distinguish the possible configurations of electrons and protons on the PQ shuttle, essentially recognizing that the difference in redox potential is well resolved despite the line-broadening induced by the interaction with the reservoirs.  17 . Here e.c.p. is a short-hand notation for electrochemical potential.
www.nature.com/scientificreports www.nature.com/scientificreports/ which are essentially an adapted version of the well known spin-boson Hamiltonian 49 featuring the two binding sites coupled by a tunneling interaction and connected to the surrounding molecular vibrational environment. Marcus transition rates are obtained by applying perturbation theory and Fermi's golden rule while assuming a Debye form for the bath spectral function ω ω ∝ J( ) / ω τ + (1 ) 2 2 (see ref. 50 for details).
Numerical solution. The dynamics of the system was numerically simulated with an original Python code.
The coupled quantum and stochastic system equations (see Eqs (10) and (12)  display two separate dominant time-scales, namely the fast (on the order of picoseconds) quantum dynamics and the slower (on the order of fractions of microseconds) mesoscopic diffusion of the PQ shuttle. At every time step µ = − dt 10 3 s, we computed the local energy values and we evolved the populations from the previous state. Notice that, at difference with the formally similar set of rate equations reported in ref. 17 , here the fully linear structure of Eq. (10) makes it possible to employ, as e.g. in ref. 11 , matrix exponentiation to compute the exact time evolution for arbitrary times x t 0 ( ) 0

Relax
We used the above expression at each step for a time dt, and we updated all the relevant observables, including the net charge on the shuttle, making use of the standard quantum mechanical formalism for the expectation values. At the end of such step, the position was updated using an Euler scheme for the integration of the SDE. The default parameters that we used in all the simulations, unless otherwise specified in the text, are summarized in Table 1: all of them are chosen in agreement with the values already known from the relevant literature and with the original work (ref. 17 ) to allow for a direct comparison. In particular, electron and proton binding energies and electrochemical potential ranges are consistent with typical redox potentials for the relevant electron transport chains and PQ/PQH 2 conversions 34 and with experimental measurements of membrane potential and pH gradient in molecular structures closely related to the b f 6 complex, e.g. the bc 1 complex in the mitochondrial resipratory chain 37 . The values and temperature dependence of the surface potential V x ( ) correspond to those reported for M. laminosus 34 , while other structural parameters such as the overall typical dimensions are retrieved from biological and crystallographic characterizations of the b f 6 complex 32,34,36 . The effective drag coefficient ζ takes into account the possible presence of other PQ/PQH 2 inside the membrane, impairing the motion of a single shuttle via random collisions 51 . Finally, the remaining parameters that are specific to our model, such as tunneling matrix elements (∆ xy ), transfer rates (γ i and Γ α ), Coulomb and reorganization energies are directly obtained from ref. 17 , where they were optimized for the Q-cycle performances at fixed temperature to reproduce the typical reaction speed and yield.
As an explicit example, we show in Fig. 6 part of a typical stochastic trajectory of the system dynamical evolution, as obtained with the default parameters. The corresponding time evolution of the electron and proton populations is also reported, for completeness. We notice that the overall dynamics is consistent with similar results reported in the literature and describing biological transport processes at the molecular level 12,13,16 .

Data availability
All the data and simulations that support the findings of this study are available from the corresponding author upon reasonable request.