Optimizing enzymatic catalysts for rapid turnover of substrates with low enzyme sequestration

TL;DR

This paper investigates how to optimize enzyme-substrate binding energies to minimize enzyme sequestration while maintaining a required reaction flux, revealing a trade-off that favors moderate binding affinities.

Contribution

It provides a theoretical analysis of enzyme binding free energies, identifying optimal values that balance sequestration and reaction rate, which was previously not well understood.

Findings

Optimal binding free energies are moderate, not too high or low.

Both substrate and product binding energies are constrained by thermodynamic considerations.

Strategies to drive reactions via binding energy differences have limited efficacy.

Abstract

We analyse the mechanism of enzyme-substrate catalysis from the perspective of minimizing the load on the enzymes through sequestration, whilst maintaining at least a minimum reaction flux. In particular, we ask: which binding free energies of the enzyme-substrate and enzyme-product reaction intermediates minimize the fraction of enzymes sequestered in complexes, while sustaining at a certain minimal flux? Under reasonable biophysical assumptions, we find that the optimal design will saturate the bound on the minimal flux, and reflects a basic trade-off in catalytic operation. If both binding free energies are too high, there is low sequestration, but the effective progress of the reaction is hampered. If both binding free energies are too low, there is high sequestration, and the reaction flux may also be suppressed in extreme cases. The optimal binding free energies are therefore neither too high nor too low, but in fact moderate. Moreover, the optimal difference in substrate and product binding free energies, which contributes to the thermodynamic driving force of the reaction, is in general strongly constrained by the intrinsic free-energy difference between products and reactants. Both the strategies of using a negative binding free-energy difference to drive the catalyst-bound reaction forward, and of using a positive binding free-energy difference to enhance detachment of the product, are limited in their efficacy.

Figures

Equations

$$\mathrm{E}+\mathrm{S}\underset{k_{-0}}{\overset{k_{+0}}{\rightleftharpoons }}\mathrm{ES}\stackrel{k_{+\mathrm{cat}}}{\to }\mathrm{E}+\mathrm{P}.$$

$$\Delta G_j=-\ln \frac{{\nu }_j}{{\nu }_{-j}}.$$

$$\Delta G=-\ln \prod _j\frac{{\nu }_j}{{\nu }_{-j}}.$$

$$\mathrm{E}+\mathrm{S}\underset{k_{-0}}{\overset{k_{+0}}{\rightleftharpoons }}\mathrm{ES}\underset{k_{-\mathrm{cat}}}{\overset{k_{+\mathrm{cat}}}{\rightleftharpoons }}\mathrm{EP}\underset{k_{-1}}{\overset{k_{+1}}{\rightleftharpoons }}\mathrm{E}+\mathrm{P},$$

$$\mathcal{R}=1-{\pi }_E={\pi }_{ES}+{\pi }_{EP},$$

$$\Psi =k_{+0}[S]{\pi }_E-k_{-0}{\pi }_{ES}.$$

$$\frac{k_{+0}[S]k_{+\mathrm{cat}}{k}_{+1}}{k_{-0}{k}_{-\mathrm{cat}}{k}_{-1}[P]}=\exp (-\Delta G)=const.$$

$$\exp (-\Delta G_{ES}/k_{\mathrm{B}}T)=\frac{k_{+0}}{k_{-0}},$$

$$\exp (-\Delta G_{EP}/k_{\mathrm{B}}T)=\frac{k_{-1}}{k_{+1}},$$

$$\exp \left((-\Delta G_{EP}+\Delta G_{ES}+\Delta \mu )/k_{\mathrm{B}}T\right)=\frac{k_{+\mathrm{cat}}}{k_{-\mathrm{cat}}}.$$

$$k_{+\mathrm{cat}}=k_{\text{cat}}\min (1,e^{-\Delta G_{EP}+\Delta G_{ES}+\Delta \mu +\Delta G_c}),$$

$$k_{-\mathrm{cat}}=k_{\text{cat}}{e}^{\Delta G_c}\min (1,e^{-\Delta G_{ES}+\Delta G_{EP}-\Delta \mu -\Delta G_c}),$$

$$\mathcal{R}=1-{\pi }_E=1-\frac{{\tau }_E}{{\tau }_{E\to E}}.$$

$$\mathcal{R}=1-\frac{1}{1+k_0[S]{\tau }_{ES\to E}+k_1[P]{\tau }_{EP\to E}}.$$

$${\tau }_{ES\to E}={\tau }_{ES}+P_{ES\to EP}{\tau }_{EP\to E},$$

$${\tau }_{EP\to E}={\tau }_{EP}+P_{EP\to ES}{\tau }_{ES\to E},$$

$${\tau }_{ES\to E}=\frac{{\tau }_{ES}+P_{ES\to EP}{\tau }_{EP}}{1-P_{EP\to ES}{P}_{ES\to EP}},$$

$${\tau }_{EP\to E}=\frac{{\tau }_{EP}+P_{EP\to ES}{\tau }_{ES}}{1-P_{EP\to ES}{P}_{ES\to EP}},$$

$$\mathcal{R}=1-$$

$$\displaystyle\frac{1}{1+k_{0}[S]\Big{(}\frac{\tau_{ES}+P_{ES\rightarrow EP}\tau_{EP}}{1-P_{EP\rightarrow ES}P_{ES\rightarrow EP}}\Big{)}+k_{1}[P]\Big{(}\frac{\tau_{EP}+P_{EP\rightarrow ES}\tau_{ES}}{1-P_{EP\rightarrow ES}P_{ES\rightarrow EP}}\Big{)}}.$$

$${\Psi }_b=\frac{(1-\mathcal{R})k_0{k}_1{k}_{\text{cat}}[E_{\text{tot}}]([S]e^{\Delta \mu +\Delta G_c}-[P])}{k_0{k}_1{e}^{\Delta \mu +\Delta G_c+\Delta G_{ES}}+k_0{k}_{\text{cat}}+k_1{k}_{\mathrm{cat}}{e}^{\Delta \mu +\Delta G_c}}$$

$${\Psi }_f=\frac{(1-\mathcal{R})k_0{k}_1{k}_{\text{cat}}[E_{\text{tot}}]([S]e^{\Delta \mu +\Delta G_c}-[P])}{k_0{k}_1{e}^{\Delta G_{EP}}+k_0{k}_{\text{cat}}+k_1{k}_{\mathrm{cat}}{e}^{\Delta \mu +\Delta G_c}}$$

$$P_{ES\to EP}=\frac{k_{\mathrm{cat}}}{k_0{e}^{\Delta G_{ES}}+k_{\mathrm{cat}}},$$

$$P_{EP\to ES}=\frac{k_{\text{cat}}{e}^{-\Delta G_{ES}-\Delta \mu }}{k_1+k_{\text{cat}}{e}^{-\Delta G_{ES}-\Delta \mu }}.$$

$$P_{EP\to ES}=\frac{k_{\mathrm{cat}}{e}^{\Delta G_c}}{k_1{e}^{\Delta G_{EP}}+k_{\mathrm{cat}}{e}^{\Delta G_c}},$$

$$P_{ES\to EP}=\frac{k_{\text{cat}}{e}^{-\Delta G_{EP}+\Delta \mu +\Delta G_c}}{k_0+k_{\text{cat}}{e}^{-\Delta G_{EP}+\Delta \mu +\Delta G_c}}$$

$$\begin{array}{rl}\Psi & =k_{+0}{\pi }_E\times \text{P}(EP\to E\text{occurs before}ES\to E)\\ & -k_{-1}{\pi }_E\text{P}(ES\to E\text{occurs before}EP\to E).\end{array}$$

$$\underset{\Delta G_{ES}\to \infty }{\lim }\Psi =k_0[S]P_{ES\to EP}-k_1[P]P_{EP\to ES}.$$

$$\underset{\Delta G_{ES}\to \infty }{\lim }{P}_{ES\to EP}=\underset{\Delta G_{ES}\to \infty }{\lim }\mathrm{const}\exp (-\Delta G_{ES})=0,$$

$$\underset{\Delta G_{ES}\to \infty }{\lim }{P}_{EP\to ES}=\underset{\Delta G_{ES}\to \infty }{\lim }\mathrm{const}\exp (-\Delta G_{ES})=0,$$

$$k_{+0}=k_0[S]\min (1,e^{-(\Delta G_{ES}-\Delta G_c^{ES})}),$$

$$k_{-0}=k_0{e}^{\Delta G_c^{ES}}\min (1,e^{\Delta G_{ES}-\Delta G_c^{ES}}),$$

$$k_{+1}=k_1{e}^{\Delta G_c^{EP}}\min (1,e^{\Delta G_{EP}-\Delta G_c^{EP}}),$$

$$k_{-1}=k_1[P]\min (1,e^{-(\Delta G_{EP}-\Delta G_c^{EP})}),$$

$$k_{+\mathrm{cat}}=k_{\text{cat}}\min (1,e^{-\Delta G_{EP}+\Delta G_{ES}+\Delta \mu +\Delta G_c}),$$

$$k_{-\mathrm{cat}}=k_{\text{cat}}{e}^{\Delta G_c}\min (1,e^{-\Delta G_{ES}+\Delta G_{EP}-\Delta \mu -\Delta G_c}).$$

$$\underset{\Delta G_{ES}\to \infty }{\lim }\Psi =k_0[S]e^{-\Delta G_{ES}+\Delta G_c^{ES}}\frac{P_{ES\to EP}{P}_{EP\to E}}{1-P_{EP\to ES}{P}_{ES\to EP}}$$

Full text

1 Introduction

Enzymatic catalysts are ubiquitous in biology, forming crucial parts of the networks that implement metabolism [1], signalling [2, 3], and the central dogma of molecular biology [4]. Analysing the mechanism by which they function is fundamental to understanding the exquisite behaviour of natural networks, to engineering existing systems [5, 6], and to developing synthetic analogs de novo [7, 8, 9].

A catalytic enzyme enhances the overall rate of a molecular process by participating in reaction intermediates. The enzyme, however, is recovered unscathed at the end of the process. This mechanism of action allows a single enzyme to turn over a large number of metabolites, but it has equally fundamental consequences for signalling and information processing systems.

Catalytic action allows an enzyme to modify its substrates in a persistent way, so that the products do not immediately convert back once they detach from the enzyme – unlike in simpler mechanisms of direct allosteric action [10]. The importance of this persistence is exemplified by kinase signalling networks, which are central to signal transduction and information processing in eukaryotes [11, 12]. In these systems, kinase enzymes catalyse phosphorylation of specific amino acid residues within protein substrates. The phosphorylated products demonstrate a change in activity, becoming, for example, activated transcription factors or kinases that in turn activate further downstream species. Crucially, the activated products do not need to remain bound to their upstream enzymes to be functional - their activation persists beyond the timescale of the catalytic interaction.

In this way, catalytic mechanisms allow the formation of (meta)stable memories in the states of downstream molecules, which retain information on their prior interactions [13]. These memories can suppress sensing noise through time-integration of signals [14]; underlie the calculation of time derivatives of an input in the context of adaptive chemotaxis [15]; and are central to the measurement and feedback cycles of recently-proposed molecular integral feedback controllers [16, 17, 18]. Furthermore, since a single catalyst can activate multiple downstream substrates, catalyst motifs can be used to split and amplify signals [19].

Although catalysts are not consumed by reactions, they must transiently participate in the intermediate complexes. In so doing they face a central paradox: catalysts must bind strongly enough to participate in the reaction, but weakly enough to be recovered at the end. Moreover, whilst participating in a reaction, enzymatic catalysts are generally sequestered by their binding partners and cannot act on additional substrates. In some cases this sequestration may allow novel, advantageous behaviour. For example, sequestration is necessary for the mechanism of “zero-order ultrasensitivity” that allows for sharp responses of the output substrate to small changes in concentration of the input catalyst [20, 21, 22]. Frequently, however, sequestration is a potential disadvantage [23]: it can limit the maximal rate of substrate turnover, and can cause “retroactive” loading effects in signalling networks that lead to breakdowns in assumptions of modularity [24, 25, 26, 27, 28, 29]. The latter is a particular concern in the rational design of synthetic systems [30].

The simple Michaelis-Menten model of enzymatic kinetics illustrates the problem of sequestration for substrate turnover [31, 32]. Consider the following reaction scheme:

[TABLE]

Here E is the catalytic enzyme, S is the substrate, ES is the enzyme-substrate complex, and P is the product. The bimolecular rate constant $k_{+0}$ describes the speed with which substrates and enzymes bind; the unimolecular rate constants $k_{-0}$ and $k_{+\rm cat}$ describe the rates of unbinding and product generation and release, respectively. Assuming the concentration of the ES complex, $[ES]$, reaches a quasi-steady state rapidly relative to any depletion of the substrate [33, 34], one obtains $k_{+0}[E][S]=(k_{-0}+k_{\rm+cat})[ES]$. The conservation law for the total enzyme concentration gives $[E]+[ES]=[E_{\rm tot}]$, implying that the concentration of enzymes in complexes with the substrate is $[ES]=\frac{k_{+0}[E_{\rm tot}][S]}{k_{+0}[S]+k_{-0}+k_{\rm+cat}}$. Therefore, the overall flux of substrates through the reaction is $r=k_{\rm+cat}[ES]=k_{\rm+cat}\frac{k_{+0}[E_{\rm tot}][S]}{k_{+0}[S]+k_{-0}+k_{\rm+cat}}=r_{\rm max}\frac{[S]}{[S]+K_{M}}$, where $r_{\rm max}=k_{\rm+cat}[E_{\rm tot}]$ is the limiting possible rate of the reaction and $K_{M}=k_{-0}/(k_{+0}+k_{\rm+cat})$. The quantity $\frac{[S]}{[S]+K_{M}}\leq 1$ reflects the effect of sequestration in quasi-steady state. It is also known as the efficiency of the enzyme catalyzed reaction [35, 36]. At low substrate concentration $[S]\rightarrow 0$, the reaction flux is proportional to $[S]$; plenty of enzymes are available to process additional substrates at the same rate per substrate. For $[S]\gtrsim K_{M}$, however, the reaction flux plateaus because a substantial fraction of the enzymes become sequestered, and fewer are available to process additional substrates.

In this article, we analyse how properties of an enzyme might be tuned, either by evolution or bioengineers, to achieve the goal of minimising sequestration while maintaining a certain reaction flux. As shown by the pioneering work of Terrell Hill [37, 38], and emphasised in the more recent works of Beard, Qian and coworkers [39, 40, 41], catalysts operate out of equilibrium and understanding the fundamental constraints on their behaviour requires a thermodynamic perspective. Numerous authors have considered the implications of thermodynamics for catalytic function in various contexts, from a generic relationship between free energy and one-way reaction fluxes [42] to inherent trade-offs between thermodynamic cost and performance in sensing and signalling [14, 19] and insulating motifs that are designed to suppress retroactive effects [29, 27]. In our case, we will assume a fixed overall thermodynamic drive to a catalytic reaction, and instead explore how the thermodynamic stability of reaction intermediates determine whether a catalyst can achieve high substrate turnover rates at low levels of sequestration.

The manuscript is structured as follows: in section 2, we justify the use of a specific thermodynamically self-consistent model of enzyme-substrate catalysis [43, 44, 45]. We then formulate our question as an optimization problem. We ask: how do we choose binding free energies of the intermediate enzyme-substrate and enzyme-product complexes so that the system minimises the number of enzyme-substrate complexes whilst maintaining a required minimum flux of reactants into products? In Section 3, we analyse this optimization problem under the assumption that association reactions are diffusion limited, finding that there is an inherent trade-off in such motifs. Choosing very high binding free energies for the intermediate complexes reduces retroactivity, but also reduces the flux through the circuit. On the other hand, choosing very low binding energies implies that the system spends a large proportion of time in the intermediate states, increasing the retroactivity of the system. We show that the optimal binding free energies are not only moderate as a consequence of this trade-off, but they are strongly related to each other. In particular, the difference between the optimal binding free energies is a constant that is related to the intrinsic free energy difference between the products and reactants. In addition, we also show that the optimal circuit saturates the bound on the flux requirement. In Section 3.2, we relax the assumption that binding rates are diffusion limited. We find that many of our observations from the diffusion-limited regime carry over qualitatively to this new regime. The optimal binding free energies are still moderate and the difference between them is confined to a value close to the intrinsic free energy difference between the products and reactants.

2 Model and methods

We now introduce basic modelling assumptions, in the process explaining why the classic Michaelis-Menten model is insufficient for our purposes. Henceforth, we will use natural units in which $k_{\rm B}T=1$, all rates are defined dimensionlessly relative to $1$ s*-1*, and all concentrations given dimensionlessly relative to 1M. We model dilute biochemical systems at the level of molecular macrostates [46]; reactions are described by mass-action kinetics with well-defined rate constants. We consider an ensemble of enzymes interacting with substrate S and product P molecules; these substrates and products are assumed to have approximately constant concentrations $[S]$ and $[P]$ on the timescale of interest [47]. In this limit, the trajectory of a single enzyme through its discrete binding states can be analysed independently as a continuous time Markov chain, with pseudo-first-order transition rates that depend on $[S]$ and $[P]$ [48]. The resultant probabilities are proportional to the expected concentrations of enzymatic states in a bulk system. The systems we consider form irreducible Markov chains, and therefore tend toward a well-defined steady-state probability distribution $\pi_{i}$ describing the occupancy of enzyme binding states $i$.

We will first illustrate our approach with the commonly-used model of Eq. 1, before arguing that it is insufficiently rigorous to allow a meaningful optimisation. We will then present the extended model that will form the basis of this work. The model of Eq. 1 has a Markov chain representation shown in Fig. 1 (a). The binding states of the enzyme are unbound (E) and substrate-bound (ES). Both the release of the product and unbinding of the substrate contribute to the same transition (ES to E) at the level of the enzyme’s binding states. The probability of the enzyme being unbound (equal to the fraction of unbound enzymes in an ensemble) is $\pi_{E}$, and the net rate of product output per enzyme is $k_{+\rm cat}\pi_{ES}$.

We consider the challenge of optimizing the enzyme properties to achieve a desired steady state rate of conversion of S into P per enzyme, $\Psi$, with a minimal steady-state fraction of sequestered enzymes, $\mathcal{R}$, at fixed concentrations $[S]$ and $[P]$. For the model in Fig. 1 (a), these quantities are given by $\Psi=k_{+\rm cat}\pi_{ES}$ and $\mathcal{R}=1-\pi_{E}$, respectively. Furthermore, since $\pi_{E}=\frac{k_{\rm+cat}+k_{-0}}{k_{+0}[S]+k_{\rm+cat}+k_{-0}}$ and $\pi_{ES}=\frac{k_{+0}[S]}{k_{+0}[S]+k_{\rm+cat}+k_{-0}}$ , it is immediately clear that the sequestration fraction $\mathcal{R}$ can be made arbitrarily small, without compromising the flux $\Psi$, by allowing the catalytic rate $k_{\rm+cat}\rightarrow\infty$. To obtain meaningful insight, it is therefore necessary to consider physically-motivated constraints on the kinetic parameters.

The most important constraint is that since the enzyme E is not consumed in the reaction, its properties cannot influence the overall free-energy change of reaction, $\Delta G=\ln\frac{[P]}{[S]}-\Delta\mu$. Here, $\Delta\mu$ is the intrinsic free-energy difference between S and P; by our sign convention, a positive $\Delta\mu$ implies P is more thermodynamically stable than S, providing a forward drive to the reaction. Note that $\Delta\mu$ can also incorporate the contribution to $\Delta G$ of the consumption of ancillary fuel molecules, such as ATP, which are treated implicitly in the model of Eq. 1. If the environment of substrates and products is fixed, an optimization over enzyme properties corresponds to optimizing at fixed $\Delta G$.

For a single reaction step $j$, the principle of detailed balance dictates that the free energy change is directly related to the forwards and backwards transition rates $\nu_{\pm j}$ [46, 49]:

[TABLE]

For a multi-step reaction, one can simply add together Eq. 2 for each step $j$, obtaining

[TABLE]

Since the catalytic step of the Michaelis-Menten model has no reverse complement, $\Delta G$ is undefined, making it impossible to perform an optimization at fixed $\Delta G$. It is therefore necessary to introduce a backwards transition, which would allow E and P to bind, and be converted into ES [50, 51, 52, 53].

Having included this reaction, it is hard to justify combining both the chemical conversion of substrate into product, and its release from the enzyme, in a single step, as in Eq. 1. If both P and S can be converted into each other by E, shouldn’t the binding and unbinding of P also be treated explicitly using a binding state EP? Indeed, ignoring EP corresponds to assuming that the enzyme-product complex is arbitrarily short-lived, yet does not present a barrier to the conversion of S into P. This assumption seems to ignore the very challenge of the optimization problem itself. We therefore use the following model [43, 44, 45],

[TABLE]

as the minimal description of catalysis in which we can meaningfully ask how enzyme properties can be adjusted to minimize sequestration $\mathcal{R}$ at fixed flux per enzyme ${\Psi}$. This molecular model can be represented as a continuous time Markov process over the enzymatic states E, ES and EP as shown in Fig. 1 (b); within this description, we obtain

[TABLE]

and

[TABLE]

The requirement of fixed $\Delta G$ equates to

[TABLE]

Even with the restriction to fixed $\Delta G$, and fixed concentrations $[S]$ and $[P]$, the optimization problem is still poorly constrained. As it stands, all rate constants in Eq. 7 could be increased by an arbitrary factor, allowing an arbitrarily high $\Psi$ whilst leaving the stationary distribution $\pi$ (and hence $\mathcal{R}$) unchanged. It would therefore be possible to obtain any flux $\Psi$ whilst ensuring $\mathcal{R}\rightarrow 0$. In practice, molecular processes cannot be arbitrarily fast due to basic constraints from physical chemistry. In any case, it is far easier to tune, whether by rational engineering or evolution, the binding free energies of the metastable substrate-bound and product-bound states than it is to further optimise the precise chemistry of the ensemble of transition states between the two. To reflect this fact, we restrict ourselves to optimizing the two standard binding free energies $\Delta G_{ES}$ and $\Delta G_{EP}$, which are related to the rate constants by

[TABLE]

A schematic representation of the the free energy profile of this process, and its dependence on the parameters in Eq. 8, is given in Fig. 2.

A given set of binding free energies, along with $\Delta\mu$, therefore specifies ratios of forwards/backwards rate constants. Absolute values remain ambiguous. To make progress, we must make further assumptions about how rate constants respond to changes in the binding free energies. A similar issue was addressed in Ref. [54], in which the authors attempted to optimize the flux of trajectories through a series of states by adjusting their energies (in that work, unlike this one, there was no attempt to minimise occupancy of intermediate states). In that case, one of the two transitions in a backwards/forwards pair was assumed to be exponentially sensitive (“labile”) to the energy difference, whereas the other was assumed to be constant.

In our system, we will initially assume (for simplicity) that the binding reactions are diffusion-controlled [55, 56]; i.e. the on-rate is fixed by the diffusion time scales that are independent of the details of the enzyme’s interaction wih substrate and product. In particular, we set $k_{+0}=k_{0}={\rm const}$ and $k_{-1}=k_{1}={\rm const}$. In the language of [54], the substrate binding transition is “backwards labile” and the product release transition is “forwards labile”. There is no immediately obvious reason to make the intermediate step of chemical catalysis either forward or backward labile. We show in Section 3, however, that both choices give pathological results for the question we ask. Invoking the fact that a true chemical modification cannot happen arbitrarily fast, we then consider an alternative in which both the forward and backward catalytic rate constants have a finite upper bound:

[TABLE]

and

[TABLE]

for some $\Delta G_{C}\in\mathbb{R}$. The scheme is based on “Metropolis dynamics”, in which reactions that are downhill in free energy have a fixed rate and uphill reactions are slowed down [57]. This dynamics is equivalent to assuming that, in an Arrhenius picture, the transition state is a fixed free energy above the maximum free energy of the two metastable states on either side. The inclusion of a finite $\Delta G_{C}$ generalises this approach to allow an offset between the maximum forward and backward catalytic rate constants. The overall effect is to split the $\Delta G_{ES}-\Delta G_{EP}$ plane into two regions: region $I$ in which the interconversion is backwards labile and the backwards step is slow; and region $II$ in which the interconversion is forwards labile while the forwards reaction is slow. A graphical illustration of these different responses to the free energies of transition is given in Fig. 3. We are now able to fully state our main optimization problem:

Question 2.1**.**

How should binding free energies $\Delta G_{ES}$ and $\Delta G_{EP}$ in the model of Eq. 4 be chosen to minimize sequestration $\mathcal{R}$, whilst maintaining a product output rate of $\Psi_{0}$, given fixed $[S]$, $[P]$ and $\Delta\mu$, diffusion-limited binding reactions, and catalytic steps with fixed and finite upper bounds on their rates (Eqs. 9 and 10)?

3 Results and Discussion

3.1 Diffusion-controlled binding rates

To answer Problem 2.1, we seek the optimal binding free energies $\Delta G^{\rm opt}_{ES}$ and $\Delta G^{\rm opt}_{EP}$ that achieve a minimal sequestration $\mathcal{R}^{\rm opt}$, whilst maintaining an output flux $\Psi^{\rm opt}\geq\Psi_{0}$. Here $\Psi_{0}\in\mathbb{R}_{>0}$ is the target seady-state net rate of substrate turnover. We will first present a mathematical analysis, followed by a physical explanation and interpretation.

3.1.1 Detailed analysis of diffusion-controlled system

We describe a target flux $\Psi_{0}$ as achievable if choices of $\Delta G_{ES}$ and $\Delta G_{EP}$ exist that satisfy $\Psi\geq\Psi_{0}$. We will show that solutions to Problem 2.1 for achievable target fluxes $\Psi_{0}$ lie on the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$, with $\Delta G^{\rm opt}_{ES}$ and $\Delta G^{\rm opt}_{EP}$ taken as high as possible to saturate the flux constraint, $\Psi^{\rm opt}=\Psi_{0}$. To make these arguments, we require machinery from the theory of Markov chains. To guide this derivation, we illustrate the continuous time Markov chain corresponding to Eq. 4, with its explicit dependence on the parameters $\Delta G_{ES}$, $\Delta G_{EP}$, $\Delta\mu$ and $\Delta G_{c}$ as laid out in Eq. 8, in Fig. 4. We first solve for $\mathcal{R}$ and $\Psi$ in terms of basic properties of this Markov chain.

Let us denote the expected lifetime of state $i$ by $\tau_{i}$; and the expected time for the next arrival at $j$ given that the current state is $i$ (the mean first passage time) by $\tau_{i\rightarrow j}$. The series of states visited by the enzyme form an embedded discrete-state, discrete-time Markov chain [58]; let $P_{i\rightarrow j}$ represent the transition probabilities of this Markov chain. $P_{i\rightarrow j}$ is then the probability that $j$ is the next state visited given that the system is in state $i$. By the memoryless property of continuous time Markov chains, the lifetime of a state is exponentially distributed with parameter equal to the total rate of outward transition from that state, implying that $\tau_{i}=\frac{1}{\sum_{j}K_{ji}}$, where $K_{ji}$ is the rate of outward transition from state $i$ to $j$.

Since we have an irreducible Markov chain, using [58, Theorem 3.8.1] and Equation 5, we obtain

[TABLE]

Using $\tau_{E\rightarrow E}=\tau_{E}+P_{E\rightarrow ES}\tau_{ES\rightarrow E}+P_{E\rightarrow EP}\tau_{EP\rightarrow E}$, and noting that $P_{E\rightarrow ES}=k_{0}[S]\tau_{E}$ and $P_{E\rightarrow EP}=k_{1}[P]\tau_{E}$, we find

[TABLE]

Observing that

[TABLE]

we can solve for the average first passage times as

[TABLE]

As a result, we can re-write the sequestered fraction in Eq. 12 solely in terms of properties of single transition steps:

[TABLE]

For the flux $\Psi$, there are two cases:

Region $I:\Delta G_{EP}<\Delta G_{ES}+\Delta\mu+\Delta G_{c}$, ${\rm ES}\rightleftharpoons{\rm EP}$ backward labile,

[TABLE] 2. 2.

Region $II:\Delta G_{EP}>\Delta G_{ES}+\Delta\mu+\Delta G_{c}$, ${\rm ES}\rightleftharpoons{\rm EP}$ forward labile,

[TABLE]

The division of $\Delta G_{ES}-\Delta G_{EP}$ space into these two regions is shown schematically in Fig. 5. For the edge case of $\Delta G_{EP}=\Delta G_{ES}+\Delta\mu+\Delta G_{c}$ that divides the two regions, both $\Psi_{b}=\Psi_{f}$ are valid. We now state a few lemmas that describe the behaviour of retroactivity and flux with respect to changing binding energies in regions $I$ and $II$.

Lemma 3.1**.**

In Region $I$, increasing $\Delta G_{EP}$ decreases $\mathcal{R}$ and increases $\Psi$.

Proof.

The transition probabilities $P_{ES\rightarrow EP}$ and $P_{EP\rightarrow ES}$ are both independent of $\Delta G_{EP}$ in region $I$, since $\Delta G_{EP}$ is irrelevant to the transitions out of state ES, and contributes the same factor to all transitions out of EP:

[TABLE]

The expected lifetime of EP, $\tau_{EP}=\frac{1}{{\rm e}^{\Delta G_{EP}}(k_{1}+k_{\text{cat}}e^{-\Delta G_{ES}-\Delta\mu})}$, decreases monotonically with increasing ${\Delta G_{EP}}$, while $\tau_{E}=\frac{1}{k_{0}[S]+k_{1}[P]}$ and $\tau_{ES}=\frac{1}{k_{\text{cat}}+k_{0}e^{\Delta G_{ES}}}$ remain unchanged. Eq. 3.1.1 therefore shows that $\mathcal{R}$ must decrease as $\Delta G_{EP}$ increases within region $I$, and Eq. 16, shows flux $\Psi$ increases with $\Delta G_{EP}$ in region $I$. ∎

Lemma 3.2**.**

In Region $II$, increasing $\Delta G_{ES}$ decreases $\mathcal{R}$ and increases $\Psi$.

Proof.

The transition probabilities $P_{ES\rightarrow EP}$ and $P_{EP\rightarrow ES}$ are both independent of $\Delta G_{ES}$ in region $II$, since $\Delta G_{ES}$ is irrelevant to the transitions out of state EP, and contributes the same factor to all transitions out of ES:

[TABLE]

The expected lifetime of ES, $\tau_{ES}=\frac{1}{{\rm e}^{\Delta G_{ES}}(k_{0}+k_{\text{cat}}e^{-\Delta G_{EP}+\Delta\mu+\Delta G_{c}})}$, decreases monotonically with increasing ${\Delta G_{ES}}$, while $\tau_{E}=\frac{1}{k_{0}[S]+k_{1}[P]}$ and $\tau_{EP}=\frac{1}{k_{\text{cat}}e^{\Delta G_{c}}+k_{1}e^{\Delta G_{EP}}}$ remain unchanged. Eq. 3.1.1 therefore shows that $\mathcal{R}$ must decrease as $\Delta G_{ES}$ increases within region $I$, and Eq. 17, shows flux $\Psi$ increases with $\Delta G_{ES}$ in region $II$. ∎

Theorem 3.3**.**

$\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$.

Proof.

For contradiction, assume not. Then, the optimal solution either lies inside region $I$ or region $II$.

Consider region $I$: assume an optimal pair of binding free energies $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$, $\Delta G^{\rm opt}_{EP}<\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$, exists that gives a minimal $\mathcal{R}=\mathcal{R}^{\rm opt}$ while satisfying $\Psi\geq\Psi_{0}$. By Lemma 3.1, it is always possible to add $\delta G>0$ to $\Delta G^{\rm opt}_{EP}$ and increase $\Psi$ while reducing $\mathcal{R}$. The resulting pair $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}+\delta G\right)$ necessarily satisfies $\Psi>\Psi_{0}$, and has a sequestered fraction $\mathcal{R}<\mathcal{R}^{\rm opt}$. The optimality of $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ is therefore contradicted.

An exactly analogous argument can be made for optimal solutions in region $II$. Using Lemma 3.2, it is always possible to increase $\Delta G^{\rm opt}_{ES}$ by $\delta G>0$ and simultaneously reduce the sequestration fraction $\mathcal{R}$ and increase the flux $\Psi$. Therefore no pair of binding free energies $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ within region $II$ can be optimal. As a result, optimal solutions must lie on the line dividing the regions, $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$.

∎

We note in passing that had we modelled the catalytic reaction as uniformly forwards labile (or uniformly backwards labile) throughout the $\Delta G_{ES}-\Delta G_{EP}$ plane, then pathological results would have been obtained. This choice would correspond to setting the whole of the plane to be region $II$ (or region $I$). As a result, it would always be possible to improve the design by increasing $\Delta G_{ES}$ (or $\Delta G_{EP}$), leading to divergent solutions that require unphysical, infinitely-fast transitions between ES and EP.

We now prove in Theorem 3.5 that if the flux constraint is achievable by the system, it is saturated: $\Psi_{\rm opt}=\Psi_{0}$. To do so, we first prove the following lemma

Lemma 3.4**.**

$\mathcal{R}$* decreases if both binding free energies $\left(\Delta G_{ES},\Delta G_{EP}\right)$ are increased by the same amount $\delta G>0$.*

Proof.

First, note that if both binding free energies are increased by the same $\delta G$, the new binding energies are in the same region of the phase plane shown in Fig. 5 as the old binding energies. Transition rates between ES and EP do not change, since the free energy change between the two binding states is unchanged. Transitions from $EP$ and $ES$ to $E$ are accelerated by $\exp(\delta G)$. As a result of these changes in rates, $\tau_{ES}$, $\tau_{EP}$, $P_{ES\rightarrow EP}$ and $P_{EP\rightarrow ES}$ all necessarily decrease for $\delta G>0$. It immediately follows from Eq. 3.1.1 that $\mathcal{R}$ decreases if both binding energies $\left(\Delta G_{ES},\Delta G_{EP}\right)$ are increased by $\delta G>0$. ∎

Lemmas 3.1, 3.2, 3.4 and Theorem 3.3 are summarised in Fig. 5.

Theorem 3.5**.**

$\Psi^{\rm opt}=\Psi_{0}$.

Proof.

For contradiction, assume an optimal pair $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ with $\Psi^{\rm opt}>\Psi_{0}$ exists. By Lemma 3.4, a pair with $\left(\Delta G^{\rm opt}_{ES}+\delta G,\Delta G^{\rm opt}_{EP}+\delta G\right)$ with lower sequestration $\mathcal{R}<R^{\rm opt}$ can be found for arbitrary $\delta G>0$. Since $\Psi$ is a continuous function, it is always possible to choose a sufficiently small $\delta G$ such that the new system has $\Psi\geq\Psi_{0}$, contradicting our initial assumption. Therefore an optimal pair with finite $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ and $\Psi^{\rm opt}>\Psi_{0}$ cannot exist. However, it is still possible that divergent values of the binding free energies lead to an ever decreasing value of $\Psi$ that nonetheless does not tend towards $\Psi_{0}$. We now argue that $\Psi$ tends exponentially towards zero for sufficiently large $\left(\Delta G_{ES},\Delta G_{EP}\right)$. As a consequence, the procedure of iteratively adding $\delta G>0$ to any candidate pair $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ with $\Psi^{\rm opt}>\Psi_{0}>0$ is guaranteed to eventually reach an improved solution with reduced $\mathcal{R}$ and $\Psi=\Psi_{0}$.

In general, the net flux through the system is

[TABLE]

This result follows from multiplying the rate of binding transitions in steady state by the probability that those binding transitions actually lead to substrate/product turnover.

In the limit $\Delta G_{ES},\Delta G_{EP}\rightarrow\infty$ at an arbitrary fixed offset, $\Delta G_{EP}=\Delta G_{ES}+\Delta G_{\rm off}$, we have $\pi_{E}\rightarrow 1$ since both ES and EP states are unstable. Moreover, the probability of transitioning between the EP and ES states (rather than to E) is suppressed by a factor $\exp(-\Delta G_{ES})$. To first order in the small quantity $\exp(-\Delta G_{ES})$, trajectories that visit either EP or ES more than once before returning to E can therefore be neglected when calculating $\text{P}(EP\to E\,\text{occurs before}\,ES\to E)$ and $\text{P}(ES\to E\,\text{occurs before}\,EP\to E)$ .

Noting additionally that $\text{rate}(E\to ES)=k_{0}[S]$ and $\text{rate}(E\to EP)=k_{1}[P]$, Equation 20 becomes

[TABLE]

Moreover, since

[TABLE]

$\lim_{\Delta G_{ES}\rightarrow\infty}{\Psi}=0$ and the flux constraint is always saturated.

∎

3.1.2 Qualitative physical discussion of diffusion-controlled system

The general behaviour outlined in Fig. 5 is exemplified by a specific system in Fig. 6, in which we show contour plots of sequestration $\mathcal{R}$ and flux $\Psi$ as a function of $\Delta G_{ES}$ and $\Delta G_{EP}$. Numerical optimization (using the code in https://doi.org/10.5281/zenodo.2656526) is used to identify points that minimize $\mathcal{R}$ at specific target fluxes $\Psi_{0}$; these points are indicated on the contour plots. As expected from the results in Section 3.1.1, the optimization saturates the bound on the target flux, and all optimal points lie on the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$. In addition, it can be seen that as the target flux $\Psi_{0}$ is increased, the optimal points move southwest towards the maximal flux, paying for the increase in flux with an increase in sequestration.

The flux plot, shown in Fig. 6 b, illustrates the central trade-off inherent to enzymatic operation. Focussing on just the substrate binding free energy $\Delta G_{ES}$, we see that the flux $\Psi$ is non-monotonic, with a peak at moderate values of $\Delta G_{ES}$. If $\Delta G_{ES}$ is too high, $\Psi$ vanishes because ES complexes immediately dissociate before the catalytic reaction can occur. On the other hand, if $\Delta G_{ES}$ is too low, ES complexes are too stable and the reactions proceed to completion extremely slowly. By requiring the system to minimise $\mathcal{R}$, we force the system to the highest possible binding free energies that can sustain the flux, since high free energies tend to reduce binding and therefore sequestration.

This central trade-off is not apparent in the Michaelis-Menten model of Eq. 1. In that model, there is no penalty to the binding free energy of the ES complex being arbitrarily low, because it is assumed that the complex can always be converted to ${\rm E+P}$ quickly. However, in our more complete model, if ES is to be rapidly converted into EP, the EP state must also be low in free energy. If this is the case, however, it will tend to frustrate the subsequent release of the product. As a result, optimal values of $\Delta G_{ES}$ are moderate.

Beyond this moderation of $\Delta G^{\rm opt}_{ES}$, we find a linear relationship between $\Delta G^{\rm opt}_{ES}$ and $\Delta G^{\rm opt}_{EP}$: $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$. In our simple model, this relationship means that reactions can occur arbitrarily fast out of the ES and EP states, avoiding sequestration, without compromising the tendency of the reactions to proceed in the desired direction (${\rm ES\rightarrow EP}$ rather than E, ${\rm EP\rightarrow E}$ rather than ES). At a deeper level, it reflects the fact that there is no point in making the product bind arbitrarily more weakly to the catalyst than the substrate does, or the catalytic step ${\rm ES\rightarrow EP}$ will never occur. Similarly, however, there is no point in driving the ${\rm ES\rightarrow EP}$ reaction forwards using a far more favourable enzyme-product complex, since this product would never be released. In fact, in the default symmetric case in which the kinetic offset parameter $\Delta G_{c}=0$ (Fig. 6 c), the difference in standard binding free energies is given by precisely the intrinsic free energy difference of the free product and substrate, $\Delta\mu$. If $\Delta\mu>0$, the released free energy can compensate for a limited increase in $\Delta G^{\rm opt}_{EP}$ relative to $\Delta G^{\rm opt}_{ES}$, allowing a somewhat enhanced rate of product release.

Hitherto, we haven’t discussed the free-energy profile for the optimal catalyst as it converts a single substrate into a product in detail. Naïvely, one might assume that to optimise the rate at which the system moves through its states, the optimal free energies would form a nice ladder of roughly evenly-spaced states. However, in their paper, in which the sole aim was to maximise flux, Brown and Sivak noted that uneven free energy drops could “compensate for differences in bare rate constants” [54]. Whilst such an effect is doubtless also present in our system, we also see that the additional need to minimise sequestration leads to the intermediate states being systematically pushed to higher free energies as illustrated in Fig. 6 (c).

3.2 Non-diffusion controlled binding rates

With respect to Equation 4, we have assumed hitherto that the binding rate constants $k_{+0},k_{-1}$ are diffusion-controlled and therefore fixed, independent of the binding free energies $\Delta G_{ES}$ and $\Delta G_{EP}$. While this approximation may be reasonable in many cases, we have effectively assumed that dissociation reactions ${\rm ES}\rightarrow{\rm E+S}$ and ${\rm EP}\rightarrow{\rm E+P}$ and can occur arbitrarily fast if $\Delta G_{ES}$ and $\Delta G_{EP}$ are sufficiently large, which is likely to be unphysical. Fortunately, our solutions predict finite values of $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ for any target flux $\Psi_{0}$, so our optimal enzymes are not inherently pathological in this fashion. Nonetheless, it is reasonable to consider the possibility that sufficiently large, positive values of $\Delta G_{ES}$ and $\Delta G_{EP}$ cause the assumption of diffusion-controlled binding rates to break down, leading to a “reaction-limited” regime in which association rates are suppressed as bonding becomes more and more unfavourable.

We incorporate this possibility by using the same modified Metropolis dynamics applied to the ${\rm ES\rightleftharpoons EP}$ transition in Section 3.1. Specifically, we assume that the ${\rm E+S\rightleftharpoons ES}$ transition is backwards labile (binding is diffusion limited) up until a binding free energy of $\Delta G^{ES}_{c}$, and forwards labile (reaction limited) above this point. Similarly, we assume that the ${\rm EP\rightleftharpoons E+P}$ transition is forwards labile (binding is diffusion limited) up until a binding free energy of $\Delta G^{EP}_{c}$, and backwards labile (reaction limited) above this point. Transition rates in the Markov model can then be written

[TABLE]

The resultant Markov process is illustrated graphically in Fig 7. In this case, there are eight possible combinations of forwards labile and backwards labile options for the three reactions. However, for any particular set of parameters, only a maximum of seven appear on the $\Delta G_{ES}-\Delta G_{EP}$ plane. This split into seven regions, rather than two as in Section 3.1, is illustrated schematically in Figure 8a. We will analyse the resultant mathematical optimization problem in Section 3.2.1, before turning to its biophysical interpretation in Section 3.2.2.

3.2.1 Detailed analysis of non-diffusion-controlled system

Our question in this setting amounts to the following:

Question 3.6**.**

How should binding free energies $\Delta G_{ES}$ and $\Delta G_{EP}$ in the model of Eq. 4 be chosen to minimize sequestration $\mathcal{R}$, whilst maintaining a product output rate of $\Psi_{0}$, given fixed $[S]$, $[P]$ and $\Delta\mu$, non-diffusion-controlled binding reactions, and catalytic steps with fixed and finite upper bounds on their rates given by (Eqs. 23)?

The analysis of regions $I$ and $II$ is identical to that in Section 3.1.1. By similar approaches, we deduce the directions on the $\Delta G_{ES}$-$\Delta G_{EP}$ plane that are guaranteed to increase flux and/or decrease sequestration in the different regions shown in Figure 8.

• •

Region $III$: Arguing as we did for Region $II$ during the proof of Lemma 3.2, one can show that increasing $\Delta G_{ES}$ decreases retroactivity but increases flux in this region.

• •

Region $VI$: Employing the same argument that we used for Region $I$ during the proof of Lemma 3.1, one can show that increasing $\Delta G_{EP}$ decreases retroactivity but increases flux in this region.

We have not found a direction in which $\mathcal{R}$ is guaranteed to decrease, and $\Psi$ guaranteed to increase, for regions $IV$, $V$ $VII$ and $VIII$. As a consequence, we get a weaker result for the augmented system: the optimal binding free energies either satisfy $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ or lie in regions $IV$, $V$, $VII$ or $VIII$.

Lemma 3.4, however, still holds for the augmented system; if both binding free energies are increased by the same $\delta G>0$, $\mathcal{R}$ necessarily decreases. The result can be separately verified for each region; increasing both binding free energies reduces a non-zero subset of lifetimes $\tau_{ES}$ and $\tau_{EP}$, and transition probabilities $P_{ES\rightarrow EP}$ and $P_{EP\rightarrow ES}$. If adding $\delta G$ to the binding free energies causes the system to move between two regions, the net effects can simply be added together. The result then follows from Eq. 3.1.1.

Similarly, we can also show that the flux constraint $\Psi^{\rm opt}\geq\Psi_{0}$ is saturated for the augmented system. Since, by Lemma 3.4, it is always possible to reduce $\mathcal{R}$ by increasing both binding free energies by the same $\delta G>0$, it only remains to be shown that $\Psi$ tends exponentially towards zero for sufficiently large $\left(\Delta G_{ES},\Delta G_{EP}\right)$. As a consequence, the procedure of iteratively adding $\delta G>0$ to any candidate pair $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ with $\Psi>\Psi_{0}>0$ is guaranteed to eventually reach an improved solution with reduced $\mathcal{R}$ and $\Psi=\Psi_{0}$.

To perform this analysis, it is necessary to consider regions $IV$ and $V$. Within these regions, consider taking $\Delta G_{ES},\Delta G_{EP}\rightarrow\infty$ at an arbitrary fixed offset, $\Delta G_{EP}=\Delta G_{ES}+\Delta G_{\rm off}$. Then

[TABLE]

This result follows from the fact that $\pi_{E}\rightarrow 1$ as $\Delta G_{ES},\Delta G_{EP}\rightarrow\infty$. Therefore the flux can be calculated as the rate for $E\rightarrow ES$ multiplied by the probability that the transition $EP\rightarrow E$ subsequently occurs before $ES\rightarrow E$, minus the equivalent term for the conversion of $P$ into $S$.

Within regions $IV$ and $V$, the only effect of increasing $\Delta G_{ES},\Delta G_{EP}$ at a fixed offset is to reduce the rates of the binding transitions by the same factor. All probabilities $P_{i\rightarrow j}$ are unchanged. Therefore

[TABLE]

and the flux necessarily tends exponentially to zero if $\Delta G_{ES},\Delta G_{EP}\rightarrow\infty$ with an arbitrary, fixed offset. We therefore conclude that it is always possible to improve any candidate pair $\left(\Delta G^{\rm opt}_{ES},\Delta G^{\rm opt}_{EP}\right)$ with $\Psi>\Psi_{0}>0$, and that continued improvements will eventually reach a solution with reduced $\mathcal{R}$ and $\Psi=\Psi_{0}$.

3.2.2 Qualitative physical discussion of non-diffusion-controlled system

The system with a crossover to reaction-controlled binding kinetics reproduces most of the biophysics observed in the simpler system of Section 3.1. Typical examples of $\mathcal{R}$ and $\Psi$ contour plots are given in Fig. 9, along with points indicating optimal designs for a range of target fluxes $\Psi_{0}$. The central trade-off that limits $\Delta G_{ES}^{\rm opt}$ to moderate values is still present, and the need to suppress retroactivity pushes binding free energies as high as possible whilst maintaining $\Psi_{0}$, leading to a free-energy profile that is not uniformly downhill.

The major difference is that optimal solutions are no longer constrained to lie exactly on the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ if the free energies are sufficiently positive. To understand why, consider what happens when we decrease $\Delta G_{ES}$ and move away from the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ into region $IV$. Generally, decreasing binding free energies increases sequestration. However, unlike in region $II$, such a move can potentially increase the flux $\Psi$, because the rate of substrate binding now increases exponentially. Against this fact, the rate of transition for ${\rm ES\rightarrow EP}$ will also be reduced, tending to reduce $\Psi$. Moving away from the line is therefore now a potentially fruitful way to trade sequestration for flux if the increase in flux due to the binding rate outweighs the negative contribution due to the decrease in the catalytic rate.

Qualitatively, the relative strength of the two contributions depends on whether this reduction in the rate of ${\rm ES\rightarrow EP}$ has a noticeable effect on the probability to proceed to the EP state once bound in the ES state:

[TABLE]

If the first term in the denominator of Eq. 26 is larger the second, then transitions ${\rm ES\rightarrow EP}$ are sufficiently rapid relative to ${\rm ES\rightarrow E}$ that decreasing $\Delta G_{ES}$ by a certain amount affects $P_{ES\rightarrow EP}$, and hence the flux, by much less than the exponential increase in the binding rate. Thus leaving the $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ line to enter region $IV$ can be beneficial.

This strategy, however, is clearly limited. A significant decrease in $\Delta G_{ES}$ will quickly suppress the catalytic step, and $P_{ES\rightarrow EP}$ will quickly tend towards $P_{ES\rightarrow EP}\propto\exp^{\Delta G_{ES}}$. At this point, there will be negligible gains in the flux by continuing to reduce $\Delta G_{ES}$, while increases in $\mathcal{R}$ will continue. We thus expect any deviation from the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ into region $IV$ to be limited.

A similar explanation can be made for optimal solutions in region $V$; in this case, a decrease in $\Delta G_{EP}$ will tend to increase sequestration, but unlike in region $I$ this decrease will tend to slow down unwanted backwards steps ${\rm EP\rightarrow ES}$ without also slowing the release of product from the $EP$ state by the same amount. This tactic is potentially effective in increasing the flux $\Psi$ only while the probability of backwards steps,

[TABLE]

remains high. Again, however, the exponential reduction in the rate of backwards steps with $\Delta G_{EP}$ will quickly reduce the incentive to decrease $\Delta G_{EP}$ further, whereas the sequestration of the enzyme by the product will continue to increase. We thus also expect any deviation from the line $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ into region $V$ to be limited.

Despite the breakdown of the tight constraint on $\Delta G^{\rm opt}_{EP}$ and $\Delta G^{\rm opt}_{ES}$, we therefore still expect the free energies to be closely linked. Moreover, this link is due to the physical principles outlined in Section 3.1: $\Delta G^{\rm opt}_{ES}$ shouldn’t be reduced too far relative to $\Delta G^{\rm opt}_{EP}$, or it will compromise the ${\rm ES\rightarrow EP}$ transition; and $\Delta G^{\rm opt}_{EP}$ shouldn’t be reduced too far relative to $\Delta G^{\rm opt}_{ES}$ or any favourable increase in the tendency for ${\rm ES\rightarrow EP}$ to be unidirectional will be outweighed by increased product binding. Randomly-generated example systems bear out the intuition of this semi-quantitative analysis (Fig. 9). Optimal points stay relatively close to the $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$ line, and importantly do not appear to continue to move further and further away from it as the target flux $\Psi_{0}$ is taken to zero. Indeed, optimal pairs $\left(\Delta G^{\rm opt}_{EP},\Delta G^{\rm opt}_{ES}\right)$ tend toward a line parallel to $\Delta G^{\rm opt}_{EP}=\Delta G^{\rm opt}_{ES}+\Delta\mu+\Delta G_{c}$, as would be predicted if the deviation from the line was constrained by a system-specific limit on the fractions in Eq. 26 and 27.

4 Conclusions

We have addressed the question of how catalytic enzymes might be designed or evolved to achieve a target rate of substrate turnover whilst minimizing enzymatic sequestration. This goal maximises the efficiency of metabolic enzymes, and is also helpful in the context of signalling and information-processing tasks. Our results, highlighting key trade-offs in the structure of intermediate enzyme-substrate complexes, should inform the way we rationalize the operation of natural systems and design synthetic analogs. These results are particularly relevant to the design of synthetic information-processing networks – eg. [59] – since sequestration-driven retroactivity [24, 25, 26, 27, 28, 29] disrupts our ability to rationally design complex circuits from simple modular components.

Specifically, we have demonstrated that investigating this question meaningfully requires the use of a more sophisticated model than the standard Michaelis-Menten description of enzymatic kinetics. Using a three-state model, with physically reasonable assumptions on the dependence of reaction rates on binding free energies, we have shown that this challenge centres around the key trade-off in enzymatic kinetics. Namely, binding to substrates should neither be too strong, since the substrate-bound states will act as stable sinks, nor should it be too weak, since the progress of the reaction will be hampered. The optimal binding free energies are therefore moderate.

We also find that the optimal binding free energies of the substrate and the product are closely related to each other. In our simplest description, their difference $\Delta G^{\rm opt}_{ES}-\Delta G^{\rm opt}_{EP}$ is a constant, related to the intrinsic free-energy difference between substrates and products. This design principle runs counter to the intuition that either: (a) $\Delta G^{\rm opt}_{ES}-\Delta G^{\rm opt}_{EP}$ should be extremely negative to favour release of the product; or (b) $\Delta G^{\rm opt}_{ES}-\Delta G^{\rm opt}_{EP}$ should be very positive to favour the enzyme-bound catalytic turnover. In fact, to maximise flux at fixed sequestration, an intermediate value is preferred, with the intrinsic chemical free energy difference between P and S potentially compensating for a slightly negative $\Delta G^{\rm opt}_{ES}-\Delta G^{\rm opt}_{EP}$.

Interestingly, our model predicts that the difference $\Delta G^{\rm opt}_{ES}-\Delta G^{\rm opt}_{EP}$ (unlike the individual optimal binding free energies) is largely insensitive to $[P]$ and $[S]$, and is an intrinsic feature of the enzyme-substrate-product chemistry. It would therefore be a natural candidate for an experimental or bioinformatics investigation exploring the question of which, if any, natural enzymes are optimized to minimize sequestration at fixed flux. It should be noted that even for biological systems that would benefit from reduced sequestration, other constraints are present in a real environment. For example, we have not considered how minimising sequestration might affect the specificity with which an enzyme interacts with one substrate rather than alternatives.

One way to phrase the above result is that the free energy released by product formation should be exploited to weaken the binding between the product and the catalyst to optimize turnover while minimising sequestration. Although our work has focused on a simple catalytic conversion of a single substrate into a product, this idea is more general. Another important mechanism of cellular information processing is the copying of template polymer sequences, as in RNA transcription and protein translation [4]. The templates involved must act catalytically to fulfill their biological function [60, 61], but current attempts to emulate these systems in synthetic contexts have struggled to generate spontaneous detachment of the products [62]; the templates effectively become sequestered by the copy polymers formed. Our results here suggest that an effective strategy would be to channel the free energy released by polymerisation into weakening interactions between template and copy, and we have recently presented a DNA reaction motif with exactly the required functionality [63].

In Ref. [54], Brown and Sivak considered the challenge of maximizing the flux of trajectories through a 3-state system; they did not simultaneously consider minimising the occupancy of the intermediate states, as we do here. They found that the ideal division of free-energy drops was dependent on the detailed functional dependence of the transition rates on free-energy differences. To meaningfully tackle our optimisation problem, it was necessary not only to make some clear assumptions about these functional forms. We had to go beyond the simple characterisation of reactions as “forwards labile” and “backwards labile” [54], by introducing limits to the speed of any one reaction – an approach that is likely to be relevant to a broad class of optimization problems in molecular systems.

By introducing these constraints, we found tighter restrictions on the shape of free energy landscapes than Brown and Sivak - in particular, the tendency of the free energy difference between intermediate states to be closely related to the overall free energy change of reaction. Moreover, the challenge of optimizing flux alone would give only the single point of maximal flux highlighted in green in Figs. 6, 9. Demanding a certain flux with minimal sequestration instead gives a family of solutions, with lower required fluxes leading to higher binding free energies for both enzyme-substrate and enzyme-product complexes. These families of optimal solutions tend to follow lines with gradient 1 in $\Delta G_{\rm ES}-\Delta G_{\rm EP}$-space, illustrating clearly the close relationship between enzyme-substrate and enzyme-product binding free energies in the context.

Our model of a catalyst with two discrete, metastable intermediate states is highly simplified, and more sophisticated descriptions with a higher number or even a continuum of intermediate states are possible. More complex descriptions of the dependence of transition rates on system parameters could also be considered, such as a smoothed version of the current sharp transition between forwards and backwards labile as the free energy of reaction is changed. A candidate would be transition rates proportional to $1-\tanh(\Delta G/2)$, where $\Delta G$ is the free energy change assoicated with the transition. Our proofs, which rely on the specific functional forms of the rates, may not hold exactly in such a smoothed model. However, we expect that smoothed models would produce results that are consistent with the overall biophysics identified here, with the smoothing region making otherwise exact results approximate. It wold also be interesting to consider time-dependent concentrations of substrate and product (rather than a fixed representative concentration for each). As a general point, however, we have argued at our model is the simplest in which the question of minimizing sequestration at fixed flux has a meaningful answer, and our approach therefore provides insight into the problem in general. Any additional detail is probably best suited to the study of specific systems in which the myriad choices necessary can be justified by reference to the particular biochemistry.

It would perhaps be of greater interest to treat the turnover of any ancillary fuel molecules, such as ATP, explicitly. In particular, we have ignored the possibility of futile cycles in which fuel molecules are consumed without converting substrate into product. It is possible that such cycles might allow reduced retroactivity at a given flux, at the expense of addition consumption of chemical free energy, analogous to the increases in molecular recognition specificity that can be achieved through “kinetic proofreading” schemes [64, 65, 66]. It would be instructive to compare such a strategy to the effect of simply increasing the intrinsic chemical free energy difference between substrates and products. Previous work [29, 27] has shown that fuel consuming approaches can be effective at reducing retroactivity when optimising at the level of the molecular network; we raise the question of whether this strategy is also effective at the level of individual enzyme-substrate-product interactions.

5 Codes

The codes for generating the figures are available at https://doi.org/10.5281/zenodo.2656526.

6 Acknowledgement

A.D. acknowledges support from Roth scholarship during his time at the Department of Mathematics in Imperial College London and Van Vleck Visiting Assistant Professorship from the Department of Mathematics at Wisconsin Madison. T.E.O. is supported by a Royal Society University Research Fellowship.