Pathway Thermodynamics Highlights Kinetic Obstacles in Central Metabolism

In metabolism research, thermodynamics is usually used to determine the directionality of a reaction or the feasibility of a pathway. However, the relationship between thermodynamic potentials and fluxes is not limited to questions of directionality: thermodynamics also affects the kinetics of reactions through the flux-force relationship, which states that the logarithm of the ratio between the forward and reverse fluxes is directly proportional to the change in Gibbs energy due to a reaction (ΔrG′). Accordingly, if an enzyme catalyzes a reaction with a ΔrG′ of -5.7 kJ/mol then the forward flux will be roughly ten times the reverse flux. As ΔrG′ approaches equilibrium (ΔrG′ = 0 kJ/mol), exponentially more enzyme counterproductively catalyzes the reverse reaction, reducing the net rate at which the reaction proceeds. Thus, the enzyme level required to achieve a given flux increases dramatically near equilibrium. Here, we develop a framework for quantifying the degree to which pathways suffer these thermodynamic limitations on flux. For each pathway, we calculate a single thermodynamically-derived metric (the Max-min Driving Force, MDF), which enables objective ranking of pathways by the degree to which their flux is constrained by low thermodynamic driving force. Our framework accounts for the effect of pH, ionic strength and metabolite concentration ranges and allows us to quantify how alterations to the pathway structure affect the pathway's thermodynamics. Applying this methodology to pathways of central metabolism sheds light on some of their features, including metabolic bypasses (e.g., fermentation pathways bypassing substrate-level phosphorylation), substrate channeling (e.g., of oxaloacetate from malate dehydrogenase to citrate synthase), and use of alternative cofactors (e.g., quinone as an electron acceptor instead of NAD). The methods presented here place another arrow in metabolic engineers' quiver, providing a simple means of evaluating the thermodynamic and kinetic quality of different pathway chemistries that produce the same molecules.


.1 Original problem
We formulate the optimization of minimal driving forces as a maximization problem (1) where S ∈ R n×m ; x, x min , x max ∈ R n ; G • , 0 (m) , 1 (m) ∈ R m ; and B ∈ R. x represents the natural logarithm of the compound concentrations divided by their standard (1M). G • represent the standard reaction Gibbs energies (in units of RT ). 0 (m) , 1 (m) are column vectors with only 0 or 1 values, respectively. In order to convert the problem to standard form, we define: Then the system in equation 1 is equivalent to the following linear problem: x is free (2)

Dual problem
According to duality theory, equation 2 is also equivalent to: which translates to the following system: where w ∈ R m , z, u ∈ R n and y = (w, z, u).
From the first constraint we get that u = Sw + z. Then, substituting u in the objective yields: The right-hand side of this objective is a weighted sum of the different dimensions of z, with positive coefficients (since x max,j − x min,j > 0). Also, the only other constraints on z are the two inequality constraints, z ≥ 0 (n) and Sw + z = u ≥ 0 (n) . We can write it in short as: z ≥ max(0 (n) , −Sw).
Since there are no other constraints and the part of the objective which depends on z is to minimize (x max − x min ) z, we get that this inequality is actually and equality: And we can also solve for u: So, we reach this simplified formulation of the dual:

Interpretation -linear pathway example
Assume that we have a linear pathway, i.e. S ∈ R m×(m−1) is a bi-diagonal matrix such as: and that x min,i and x max,i are the same for all compounds. Note that −Sw = ∆w, which is the difference between w and the same vector shifted by 1 (so that ∆w i = w i − w i−1 and we pad the edges with zeros). The dual problem then becomes very simple: In the following sections, we will prove that the solution for this problem will always be of the form w = (0, . . . , 0, φ −1 , . . . , φ −1 , 0, . . . 0), where φ is the number of non-zero elements. In other words, w will contain a single block of non-zero values where all the weights are the same. To show this, we will first prove that any series of consecutive non-zero weights must all have the same value (lemmas 1-2). Then we will conclude that only one such block can exist. Lemma 1. Let w be a non-degenerate feasible optimal solution to the linear problem in (6). Suppose that there is an index i such that the only non-zero values are w i and w i+1 . Then it must be that w i = w i+1 .
x Proof. We will assume that w i = w i+1 and show that this contradicts the optimality of w.
Define D ≡ min(w i , w i+1 ) > 0. Now, for any |δ| < D let the vector w be changed tow ≡ w+δ(e i −e i+1 ), where e i is the ith unit vector. This new vector satisfies the constraints of the linear system (all elements are non-negative and sum up to 1).
Importantly, if δ < |w i − w i+1 |, then sign(∆w) = sign(∆w). This means that the objective of the linear problem f (w + δ(e i − e i+1 )) will be a linear function of δ within this range. This fact is a contradiction to the optimality (and non-degeneracy) of f (w), i.e. when δ = 0, since the optimum of a linear function can only be at the boundary of its range.
Lemma 2. Let w be a non-degenerate feasible optimal solution to the linear problem in 6. If there exist p < q where w p−1 = w q+1 = 0 and ∀i = p, . . . , q : w i > 0 then all these non-zero values are the same, i.e. w p = w p+1 = . . . = w q .
Proof. We will again assume that these values are not all the same and then reach a contradiction. Assume that there is p ≤ i ≤ q such that w i = w i+1 . Then using the same idea as in the proof of Lemma 1, we can define a vector d that will be −x for all indexes between p and i and +y for all indexes between i + 1 and q, and that will sum up to 0. Then, similarly, we can definew ≡ w + δd and find a small enough δ that will not change the sign of ∆w. This would mean that f (w + δd) will be linear in δ and this would again lead to a contradiction.

Conclusion
Lemma 2 means that if w is a non-degenerate feasible optimal solution then it is a series of "blocks" with zeros between them. Then, using a similar proof to Lemma 2, we can also show that only one such block can exist.

Solving the simplified dual for linear pathways
By assuming w is of the form described in the above conclusion, the only two parameters defining it are the first and last index of the plateau (e.g. p and q). The height is just 1 over the length, i.e. 1/(q − p + 1). The x max,q penalty is added only once (for the first index of the plateau) and the −x min,p+1 is added once for the last index. Therefore the dual for linear pathways is equivalent to solving: 2 Max-min driving force and minimal enzyme cost 2.1 Max-min driving force calculation is equivalent to cost minimization of substratesaturated enzymes Consider a linear pathway with N uni-uni reactions, where each product is the substrate of the consecutive reaction. All N enzymes have the exact same k + cat , are substrate saturated, and are not allosterically regulated. The standard Gibbs energies of the reactions (in units of RT ) are stored in the vector G • ∈ R N . The log-scale concentrations of the metabolites are stored in the vector x ∈ R N +1 satisfying the constraints x min ≤ x ≤ x max . The enzyme abundances are given by the vector E ∈ R N . Now, the vector of reaction Gibbs free energies (g ∈ R N +1 , also in units of RT ) reads g i = G • i +x i+1 −x i , and according to the second law of thermodynamics, we require that g < 0. Using the separable rate law for simple MM reactions and assuming full substrate saturation (as in section 2.3.2 in [1], i.e. [S] K S and . We assume a predefined, stationary flux v 0 through the pathway. To realize this flux, the enzyme abundances must satisfy and the total enzyme cost will be One can assume that evolution can adjust the enzyme abundances, and thereby the metabolite levels, in a way that minimizes the total enzyme cost per unit flux. To solve for this, we simply need to minimize E with respect to c while fixing v 0 , k + cat and G • . This is a convex optimization problem. After dropping the constants v 0 and k + cat ), it can be formulated as: The minimal driving force achieved by the minimal-enzyme cost algorithm can never be better than the value obtained by the MDF algorithm, by definition. This is trivial. However, shortly we will prove that it will also not be worse.

Proof
Outline of the proof If we optimize for the total enzyme cost, the pathway will contain some stretches of reactions with the following properties: • The first metabolite hits the upper concentration bound.
• The last metabolite hits the lower concentration bound.
• All reactions in the stretch have the same reaction Gibbs free energy (GFE).
• The reaction before the stretch and the reaction following the stretch have lower reaction GFEs.
• One of these stretches will have the highest GFEs, and this GFE cannot be decreased any further.
Therefore, any minimal enzyme cost solution is a solution of the MDF problem.
Lemma 3. Let x * be the solution to equation (9). The following holds ∀i = 2 . . . N : Proof. Since x * is a minimum of E, then ∀i ∈ {2 . . . N } the partial derivative ∂E ∂xi must either be equal to 0, or x i must be bound such that moving it in the direction of the negative slope is impossible. Explicitly, Since x i affects only g i−1 and g i , the partial derivative of E with respect x i is Next, we note that the function h(x) = e x /(1 − e x ) 2 is monotonically increasing (in the range x < 0). Therefore, for any x, y < 0 we conclude that sign(h(x) − h(y)) = sign(x − y), and in our case The conclusion of the lemma follows directly.
Rest of the proof Using Lemma 3, we can show that the x * consists of stretches of reactions with equal g i s where the concentration of the first substrate is x max,q and the concentration of the last product is x min,p+1 . If there is a reaction right before or after this stretch, its Gibbs energy is more negative. The sum of all the Gibbs energies of the reactions in this stretch is and is a constant. Therefore, there is no way to decrease the Gibbs energy of all the reactions within the stretch simultaneously. Furthermore, since all g i values are equal, we can calculate it using: and it means that this solution also minimizes the value of the maximal g i within this stretch of reactions. Therefore, the driving force in the thermodynamic bottleneck reactions, min(−g i ), will be equal to the MDF.
3 Reactions close to equilibrium exert little flux control in simple metabolic pathways In this section, we use Metabolic Control Analysis to show a relationship between driving forces and flux control coefficients for a simple but prototypical case. Consider a linear chain of uni-uni reactions in steady state, all supporting the same net flux J. As analyzed in the main article, a reaction with a low driving force is expected to require a high level of its enzyme (assuming all kinetic parameters are similar). Furthermore, the forward and backward fluxes in such a reaction are significantly higher than its net flux J. Thus, if the enzyme level of a low driving force reaction is increased, the forward and backward conversion fluxes will increase almost by the same proportion. Hence, the concentrations of the reactants will further approach their equilibrium ratio and therefore the driving force of the reaction will go even lower. The final outcome is that the potential increase in the net flux is dampened considerably and therefore the effect on J will be minimal. We therefore claim that the reactions with the smallest driving force exert relatively little control over the pathway flux and are thus not expected to be used as regulation points.

Control coefficients and reaction elasticities in a linear pathway
To show this more rigorously, we use the flux connectivity theorem from Metabolic Control Theory [2]. For the scaled flux control coefficients C Now we consider a linear chain of uni-uni reactions. The first reaction converts metabolite 1 into metabolite 2, the second reaction converts metabolite 2 into metabolite 3, and so on. Since each intermediate metabolite (index i) affects only the previous and the following reaction (with indices i and i + 1), it has two non-zero elasticities, and the connectivity theorem for the pathway flux J and this metabolite reads so Thus, the ratio of two subsequent flux control coefficients is given by the negative ratio of the reaction elasticities. The same formula also holds for unscaled control and elasticity coefficients.

The elasticities in separable rate laws
The reversible Michaelis-Menten rate law can be rewritten as a product of three terms [1]: where V + = E k + cat is the maximal rate capacity of the enzyme, the binding term κ (ranging between 0 and 1) describes the effects of incomplete substrate saturation, and the thermodynamic term γ (also between 0 and 1) describes the effect of microscopic backward fluxes. The maximal rate capacity V + does not depend on metabolite concentrations. Together with the kinetic term, it yields the microscopic forward rate v + i = V + i · κ i of the reaction. The flux-force relationship [3] states that for every reaction, the ratio of forward and backward fluxes is determined by the Gibbs free energy: v + i /v − i = e −gi , and so the ratio between net flux and forward flux reads This is exactly our thermodynamic term γ i (see [1]), which is therefore strictly a function of the reaction's Gibbs energy change g i (in units of RT ). Its derivative with respect to g i reads and the derivative with respect to the logarithmic reactant levels (x j ≡ ln c j ) can be written as From the definition of g i we get that where the n ji are the stoichiometric coefficients for substrates (negative) and products (positive) in reaction i. We can now have a general expression for the scaled elasticities: Taking everything together, the scaled elasticities can therefore be written as

Case 1: Substrate-saturated enzymes
In this section, we assume that all enzymes in the linear pathway are fully substrate saturated (as in section 2.3.2 in [1], i.e. [S] K S and [P ] K P ). In this case, the binding term κ i for all of the enzymes is insensitive to any of the compound concentrations, and hence: Therefore, from Eq.(21) we get that: For a chain of uni-uni reactions, all stoichiometric coefficients are ±1, and the flux control coefficients between subsequent reactions i − 1 and i (and the metabolite c i in between) show the ratio Thus and, considering the sum theorem, we obtain As the driving force of a reaction approaches 0, so does its control coefficient. If a single reaction is completely irreversible (infinite driving force), its control coefficient vanishes. However, if several reactions are completely irreversible, the control coefficients will be undefined because slight changes of enzyme levels -which play a central role in defining the control coefficients -would create a situation in which no steady state is possible.

Case 2: Enzymes operating in linear range
As a second limiting case, we consider enzymes which have substrate concentrations much lower than the Michaelis-Menten constant, i.e. x j K M,j . In the formulation used in Eq.(16), we can approximate ∂ ln κ i ∂x j = n ji j is a substrate of reaction i 0 j is a product of reaction i .
Therefore, Eq. (21) yields and for a chain of uni-uni reactions, all stoichiometric coefficients are ±1, so we obtain The flux control coefficients between subsequent reactions i − 1 and i (and the metabolite c i in between) show the ratio Thus The formula shows that driving forces affect the control coefficients in two ways: on the one hand, according to the numerator, strongly driven reactions (large value of −g i ) should have a larger control, and reactions with small driving forces (−g i just above 0) should have little control. On the other hand, the magnitude of the denominator increases along the chain, leading to a tendency for lower control as a reaction is more downstream. In particular, a reaction with large driving force decreases the relative control of all following reactions, not just because the control coefficients must sum to 1, but also because it increases the denominator term for all following reactions (see example in Figure S3).

Consistency with existing formula for flux control coefficients
The formula (33), for flux control coefficients of enzymes in their linear range, agrees with similar, existing formulae. On the one hand, driving forces g i and mass-action ratios ρ i are related by ρ i = e −gi . Accordingly, Eq. (33) agrees with Kacser and Burns' formula from their seminal work on Metabolic Control Analysis [4]. On the other, it is related to a formula from [5] which relates the flux control coefficients to enzyme levels and kinetic constants. Consider again a chain of n reactions (external metabolites concentrations are x 1 and x n+1 , internal metabolites are x 2 . . . x n ). As before, the vector G • represent the standard reaction Gibbs energies (in units of RT ). From the exact flux formula (see [6]) follows, by direct differentiation, a formula for the scaled flux control coefficients (compare the equivalent formula in [5]). The absolute values are defined, as usually, by the summation theorem i C J vi = 1. This formula depends only on kinetic constants and the enzyme levels, and neither metabolite levels nor driving forces appear explicitly. However, it is equivalent to our formula. We first note that the term ci , we rewrite the concentration of the i + 1 th metabolite in terms of driving forces:

Summary
In the above two limiting cases, there is a direct relationship between the driving force of a reaction, −∆ r G(i), and its control coefficient on the pathway flux, C J vi . The kinetic parameters of the enzyme, as well as its abundance, do not affect the control coefficient directly, but can still indirectly affect it through changing the driving forces at steady-state. We believe that this trend holds in most kinetic models, i.e., even for non-linear pathways and for cases where some enzyme are neither saturated nor operating in the linear range.
When thinking about the cost of maintaining a certain pathway flux, in terms of enzyme abundances, we would like to focus our analysis on the effect of changing the level of one enzyme. In Metabolic Control Analysis, the sensitivity between an enzyme level and a flux is termed a response coefficient of the enzyme level. Since we constrain ourselves to kinetic models where the rate of each reaction is proportional to the amount of the enzyme catalyzing it, the scaled response coefficients are identical to the scaled control coefficients, C J vi . Therefore, the response coefficient is higher for reactions that are far from equilibrium and lower for the reactions of minimal driving force. Increasing the enzyme level in such a reaction will have a weaker effect on the pathway flux than increasing the level of an enzyme with a high driving force.