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 [Formula: see text] 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.

Introduction

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-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 () 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 complex, which will be the object of our analysis in the following of this work.

Figure 1

Model for electron transport chain and cytochrome complex. (a) The electron transport chain of oxygenic photosynthesis. Solid purple arrows show the ideal linear electron flow (LEF) in which electrons go through the whole chain and are finally transferred to a biological redox agent called NADPH. Dashed purple arrows show the alternative path that makes cyclic electron flow (CEF) possible. (b) Molecular structure of the cytochrome complex (protein data bank, accession code PDB-1VF5) with a description of the Q-cycle mechanism under CEF conditions. The electron recycling chain (heme and heme prosthetic groups) is marked with blue arrows, and biological counterparts for the elements considered in the model are reported. (c) 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 models the electron re-injection site from the Fd pool to , while collectively represents the ISP-cytochrome 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 and , reported in blue consistently with the color code in panel (b), represent the heme and groups respectively. The black arrows show the ideal path of the electrons in the system.

The cytochrome 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-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 = 1/2 formulation, in which the occupied Fock state 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 ), while we will use capital letters for protons (e.g. annihilation operators ). Annihilation operators obey anticommutation rules and . For simplicity of notation, we will assume in the following.

In Fig. 1b,c we provide a schematic description of the working principles of the cytochrome 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-31]. Under CEF, electrons transferred on the P-side to PSI are reinjected into the 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 complex exist in the literature[32-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 towards the interior of the thylakoid (lumen) and on the other side (stroma), where specific mobile electron carriers, called plastoquinone/plastoquinol (/) can fit one per site at a time. Each 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 . In principle, a whole pool of / 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 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 and 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 (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 reduction inside the complex is effectively balanced by a 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 enters the reaction site, it interacts with two electron binding sites, namely heme- and the Iron-Sulfur (ISP) domain ( in Fig. 1c). One electron leaves the and is transferred through the ISP via cytochrome 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 and then to heme 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 site close to heme . The electron that traversed the - 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 complex. Two protons are taken up from the stromal side, resulting in a fully reduced , ready to diffuse back into the membrane and towards the P-side. The hemes and that constitute an electron-recycling chain are again described by a composite Hamiltonian containing free energy terms, , accompanied by a Coulomb repulsive interaction, , which reflects the redox anticooperativity of the two hemes, as measured in molecular complexes structurally and functionally similar to the [37]. Referring to Fig. 1c, the site collectively represents the ISP-cytochrome chain that transfers electrons to the Pc pool. On the other hand, site models the electron re-injection site from the Fd pool to . Each of the and 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 / 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. . We hereby describe and sites interacting with the electron source and drain in terms of typical Lindblad terms (with )[38] where the Lindblad dissipator is , and the average occupation numbers of the reservoirs are described by Fermi-Dirac distributions at the and energies where is the temperature and 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, , , and add the following contributions to the master equation to describe the transition from state to state via the tunneling connection from site to site . 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 is the energy eigenvalue of state as an eigenstate of the free Hamiltonian of the system, and 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 , and is the absolute value of the energy difference between the states and . If we now define we can write the full evolution of the system in the form of a single master equation where we impose if sites and 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 where is a vector with components describing the populations (diagonal elements of ), and the time-independent relaxation matrix, , receives contributions from all the components of Eq. (9). Indeed, the latter can in general be recast in the form[38] from which we see that, in our particular case, . 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 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 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, ) and represents a gaussian noise source with zero average and correlation function . The potential is added in order to confine the shuttle inside the membrane, while 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 , , , and 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 and the upper sign for . Here  nm, and 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 where the upper sign refers to and the lower to , 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 complex, which is taken into account by adding an internal potential between the membrane boundaries . Such a potential is linearly varying between two temperature-dependent extreme values  T and  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 s of 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], in such a way to guarantee a consistent comparison of the results for benchmarking our theoretical method. All such values are reported in the Methods section. First, the electrochemical gradient difference for protons, , was gradually varied while symmetrically keeping . Since the following relation generally holds[41] at a fixed temperature the scan over is equivalent to a scan over pH, namely the concentration gradient of . 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, . 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 was varied at constant and . At room temperature ( K) we have and with . Following Smirnov and Nori, we varied while keeping 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 even at the extremes of the parameter range.

Figure 2

Benchmark simulations. (a) Quantum Yield as a function of the electrochemical gradient of protons. (b) Quantum Yield as a function of the electrostatic surface potential. For comparison, in both panels we report the original data extracted graphically from the plots published in ref.[17] using the WebPlotDigitizer software[52].

Figure 3

Quantum Yield as a function of temperature. (a) Quantum Yield as a function of for fixed of protons and (scheme I) (b) Quantum Yield as a function of for fixed pH of protons and (scheme II).

Figure 4

Temperature scan at constant pH and varying surface potential. (a) Quantum Yield as a function of for fixed pH of protons and with varying with temperature (scheme III). (b,c) Average total number of translocated electrons () and protons () during a full trajectory under scheme III.

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 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 shuttle is also affected. According to Eq. (15), the pH gradient also depends on 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 , is obtained from Eq. (15) and then distributed as /2 and /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], this dependence takes the form and . Again, we kept pH fixed and plotted the results in Fig. 4. The connection between and pH was not modified, assuming that ion balance across the thylakoid membrane always provides a net zero electrostatic component to the protonic electrochemical potential. The -scan simulations show that the quantum yield decreases while increasing the temperature, while the number of translocated protons increases. If achieving a quantum yield as high as possible is certainly beneficial, this comes at the cost of a reduced number of transferred protons per unit time: in the biological optimization perspective, a trade-off should then be found at the point where the activity is both sufficiently quick and efficient. A possible figure of merit is the product which we expect to show a maximum in the intermediate range of temperatures. This is indeed the case, as shown in Fig. 5a for the data of simulation scheme III. The maximum of is found in a temperature range that broadly lies around physiological conditions.

Figure 5

Figures of merit under scheme III. (a) Quantum Yield times the number of translocated protons. (b) Thermodynamic efficiency.

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 scan at constant pH and -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 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 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.

In perspective, the present model is intended as a first step towards a full description of the 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 complex to the introduction of specific inhibitors, pointed towards the necessity of explicitly formulating some gating mechanisms[42-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 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-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-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 shuttle. As previously stated, these structural elements are chosen consistently with the original kinetic model in ref.[17] and constitute the basis for our original developments based on an open quantum systems formalism. The explicit expression of the full system Hamiltonian can be written as For the shuttle, we have where the electron and proton Hamiltonians are (here the indeces refer to the two different binding sites for electrons, with operators , and protons, operators ) 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: where and . The Hamiltonian for the and sites also contains free energy terms and a Coulomb repulsive interaction Finally, and are treated as single-electron binding sites: The Ferredoxin and Plastocyanin pool on the stromal and lumenal side respectively are modeled as collections of fermionic oscillators 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 The effect of is to rescale the free excitation energies of the electron and proton binding sites. For sites , , and this is just a static shift 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 . In other words, the stochastic random force originating from molecular collision is dominant in governing the diffusion of the 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 shuttle, essentially recognizing that the difference in redox potential is well resolved despite the line-broadening induced by the interaction with the reservoirs. Finally, the Marcus transition rates can be derived starting from Hamiltonian terms of the form

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 / (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) in the main text) 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 shuttle. At every time step 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 We used the above expression at each step for a time , 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 / conversions[34] and with experimental measurements of membrane potential and pH gradient in molecular structures closely related to the complex, e.g. the complex in the mitochondrial resipratory chain[37]. The values and temperature dependence of the surface potential 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 complex[32,34,36]. The effective drag coefficient takes into account the possible presence of other / 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 (), transfer rates ( 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.

Table 1

Summary of the relevant parameters and of the corresponding default values assumed for the numerical simulations performed in this work, chosen consistently with ref.[17]. Here e.c.p. is a short-hand notation for electrochemical potential.

Parameter	Description	Value	Parameter	Description	Value
ϵ_Q0	electron binding energy on PQ	280meV	E_Q0	proton binding energy on PQ	822meV
μ_Fd	e.c.p. of Fd pool	410meV	μ_Pc	e.c.p. of Pc pool	−440meV
μ_N	e.c.p. of N-side protons	−75meV	μ_P	e.c.p. of P-side protons	75meV
ϵL′	Site L electron binding energy	360meV	ϵH′	Site H electron binding energy	220meV
ϵA′	Site A electron binding energy	465meV	ϵB′	Site B electron binding energy	−495meV
Δ_AQ	A-PQtunneling rate	0.10meV	Δ_BQ	B-PQtunneling rate	0.10meV
Δ_LQ	L-PQtunneling rate	0.06meV	Δ_HQ	H-PQtunneling rate	0.06meV
Δ_LH	L-H tunneling rate	0.10meV	λ_AQ	A-PQreorganization energy	100meV
λ_BQ	B-PQreorganization energy	100meV	λ_LQ	L-PQreorganization energy	100meV
λ_HQ	H-PQreorganization energy	100meV	λ_LH	L-H reorganization energy	250meV
V_P	surface potential on P-side	140meV	V_N	surface potential on N-side	120meV
U_ee	electron repulsion energy on PQ	305meV	U_pp	proton repulsion energy on PQ	76.30meV
U_ep	electron-proton attraction energy	610meV	U_LH	L-H electron repulsion energy	240meV
γ_Fd	Fd pool electron transfer rate	10^−4meV	γ_Pc	Pc pool electron transfer rate	10^−4meV
Γ_P	P-side proton transfer rate	0.002meV	Γ_N	N-side proton transfer rate	0.002meV
U_ch0	hydrophobic confining energy	770meV	U_w0	membrane confining energy	500meV
x_ch	hydrophobic confining position	1.70nm	x_w	membrane confining position	2.70nm
l_ch	hydrophobic confining length	0.05nm	l_w	membrane confining length	0.10nm
l_p	proton coupling decay length	0.25nm	l_e	electron coupling decay length	0.25nm
T	temperature	25meV	ζ	drag coefficient	8.55meVμsnm^2

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].

Figure 6

Typical stochastic trajectory of the model system. (a) Position of the shuttle inside the membrane. (b–d) Electron () and proton () average population on the shuttle and on the , , and sites. Notice the spikes corresponding to the charging and discharging of the shuttle when it approaches the N- and P-side, as well as the remarkable coupling between electron and proton motion (the latter are translocated against their electrochemical gradient). (c) Total number of electrons transferred to the Pc pool () and total number of protons () brought from the N- to the P-side.