An approach to learn regulation to maximize growth and entropy production rates in metabolism

Elucidating cell regulation remains a challenging task due to the complexity of metabolism and the difficulty of experimental measurements. Here we present a method for prediction of cell regulation to maximize cell growth rate while maintaining the solvent capacity of the cell. Prediction is formulated as an optimization problem using a thermodynamic framework that can leverage experimental data. We develop a formulation and variable initialization procedure that allows for computing solutions of the optimization with an interior point method. The approach is applied to photoheterotrophic growth of Rhodospirilium rubrum using ethanol as a carbon source, which has applications to biosynthesis of ethylene production. Growth is captured as the rate of synthesis of amino acids into proteins, and synthesis of nucleotide triphoshaptes into RNA and DNA. The method predicts regulation that produces a high rate of protein and RNA synthesis while DNA synthesis is reduced close to zero in agreement with production of DNA being turned off for much of the cell cycle.

1 Introduction

Biological systems can be understood as dissipative systems analogous to tornadoes and hurricanes. Dissipative systems act by taking the most probable path to reduce the energy difference in the environment. These most probable paths result in cyclical patterns of material movement that act to transport the material from regions of high energy to regions of lower energy. In tornadoes, air movement becomes correlated and cyclical in which hot air is driven up and cool air down. Whether one considers entropy from a thermodynamic or information standpoint, the concept is the same—the most probable path allowed by physical constraints is taken. It is in this respect that we expect the dynamics of biological systems act to maximize their entropy production rates to the extent possible. In this regard, biological processes are highly constrained by physical and biological limitations—constraints that are critical for their function. Moreover, for biological systems, the adaptation of their dynamics to maximize their entropy production rates occurs on the time scale of natural selection. This is to say, through natural selection their genome structure and metabolism become updated such that they act to move to the most probable state possible by taking the least action path. On the shorter timescale of a cell’s lifetime, those adaptations also include complex mechansims for adapting to changing environmental conditions. These adaptation mechanisms serve to control the dynamics in such a way that their ability to carry out auto-catalysis (Hinshelwood, 1952) is preserved. Adaptation is accomplished through regulation of metabolic reaction pathways such that organisms remain viable in a physico-chemical sense despite rapid or dramatic shifts in the environment.

Because regulation is at the heart of understanding biological processes, discovery and understanding of the basis for regulation is critical to the field of biology. Regulation may be in the form of just-in-time regulation, graded regulation, switch-like regulation, pre-programmed regulation such as the circadian clock system, or even the choice to not regulate when it is too costly to alter the respective enzyme expression but to instead constitutively turn on an activity (Sivak and Thomson, 2014). Regulation is an important aspect of fitness in that organisms and even individual cells must regulate themselves so that they act as efficiently and as quickly as possible, otherwise they are outcompeted in their environment. Consequently, the environment over the time of evolution shapes regulation. Those cell phenotypes that operate in the most efficient manner in nutrient poor and dramatically shifting environments will tend towards fixation in the population. In rich environments in which selection pressure is not strong, dysregulation of cells can result in uncontrolled growth, manifest as cancer in metazoans.

Regulation in biological systems has been historically uncovered using experiments carefully designed to test hypotheses, typically employing isotope labeling methods to track reaction fluxes. Such meticulous but labor intensive studies have been quite successful at uncovering metabolic regulation (Larsson-Raznikiewicz, 1967; Newsholme and Start, 1973; Waygood and Sanwal, 1974; Waygood et al., 1975; 1976; Aithal et al., 1985; Jitrapakdee and Wallace, 1999; Holness and Sugden, 2003; Lehninger et al., 2005; Hallows et al., 2012; Saha et al., 2014; Wang et al., 2014), yet we are far from a complete understanding of regulation in even model organisms. High-throughput methods to uncover regulation are highly desirable, and recent work in this area has made progress (Hackett et al., 2016; Reznik et al., 2017), but we are still a long way from a routine method for determining regulation, either by experiment or prediction.

In a recent study by Hackett et al. (2016), Michealis-Menten steady-state kinetic models with regulation were used to predict reaction fluxes and concentrations. The predictions were compared to reaction fluxes inferred from 13C isotope experiments and concentrations derived from mass spectrometry and NMR measurements (Hackett et al., 2016). The correlation between simulation-predicted fluxes and experimentally-inferred fluxes were evaluated with and without regulation in the simulation. If the match was better with regulation, then regulation was assumed. The work was a tour de force in that 25 chemostat studies were used to carefully measure both absolute and relative metabolomics data while at the same time cover as much of the proteome as possible.

An approach used by Reznik et al. (2017), has less reliance on multimodal experimental designs, and instead used sophisticated informatics to develop a model of small molecule regulatory networks from curated databases of enzymes. They integrated the regulatory network with a metabolic model of Escherichia coli, and distilled information on how substrates and inhibitors contribute to metabolic flux regulation (Reznik et al., 2017).

More recently, Britton et al. (2020) have developed optimization and reinforcement learning methods that predict which enzymes need to be regulated to maintain metabolite concentrations at a level such that the diffusion process within the cell remains viable. These methods build on work by Cannon et al. (2018) that takes advantage of a maximum path-entropy/caliber (Jaynes, 1985; Dewar, 2009; Dixit et al., 2018) formulation of the law of mass action in order to predict likely metabolite concentrations. The advantage of these approaches is that minimal measurements and few experimentally-derived parameters are required to predict regulation.

Maximum entropy methods are attractive because they use the most likely values of the parameters needed to model a process (Jaynes, 1985; Dewar, 2009; Dixit et al., 2018). That is, rather than a large search through parameter space to determine the most likely kinetic parameters, a maximum entropy formulation can provide the solution directly, then if needed, reaction parameters can be back calculated. Strictly speaking, these methods maximize the path entropy of a system, meaning that they maximize the entropy subject to the stoichiometric constraints imposed by the reactions and the constraints due to the system boundary conditions. Without the stoichiometric constraints, each molecular species would move to its natural abundance determined by its standard chemical potential. Without non-equilibrium boundary constraints, the reactions would all go to the equilibrium state.

In the study by Britton et al. (2020), regulation was inferred only with regard to which reactions needed to be controlled to keep metabolites at concentrations that are not so high that they would impede diffusion (Britton et al., 2020). As mentioned above, in an environmental niche with finite nutrient resources for growth, for cells to persist they must replicate in a timely and efficient manner in order to compete with other species. Therefore in the context of evolution we expect to see regulation in cells that maximizes growth, and only produces those enzyme catalysts that contribute to growth in a timely and efficient manner.

In this study, we develop an optimization approach that allows for extension of the work by Britton et al. (2020) to be applied for more general constraint based modeling, and demonstrate our method for prediction of regulation that maximizes growth pathways while still constraining metabolite levels to physiological values. While flux balance type analyses optimize a similar objective and may include information about regulation (Covert et al., 2001) or even thermodynamic feasibility constraints (Henry et al., 2007), the method we present is fundamentally different in that we use a thermodynamic perspective to directly infer the most likely distribution of both fluxes and metabolite concentrations subject to given constraints.

We apply our method, that we term pathway-controlled optimization (PCO), to the metabolism and enzymatic activities of Rhodospirilium rubrum, a purple non-sulfur photosynthetic bacterium that is being used as a synthetic biology organism for the purpose of ethylene production (North et al., 2020). We compare the regulation and reaction fluxes predicted by this new method to the approach presented in (Britton et al., 2020) that adjust regulation only to maintain metabolite concentrations at physiological levels.

2 Methods

We begin this section with a brief introduction to the maximum path entropy approach developed in (Cannon et al., 2018) and (Britton et al., 2020). Then we present how we extend this approach to account for regulation in cells to optimize growth.

2.1 Maximum path entropy solution without regulation

For a set of molecular species $\mathcal{I}=\{A,B,C,D\}$ having respective concentrations n_i for each species $i\in \mathcal{I}$, participating in a set of reactions $\mathcal{J}=\{1,-1\}$, with unsigned stoichiometric coefficients ν_i,j for each molecular species i and each reaction $j\in \mathcal{J}$, the reversible reaction is described by the chemical equation,

${\nu }_{\mathrm{A},1}{n}_A+{\nu }_{\mathrm{B},1}{n}_B\underset{k_{-1}}{\overset{k_1}{\rightleftharpoons }}{\nu }_{\mathrm{C},1}{n}_C+{\nu }_{\mathrm{D},1}{n}_D.(1)$

Here, k₁ and k₋₁ are the reaction rate parameters.

Let $|\mathcal{I}|$ denote the size of the set $\mathcal{I}$ and let $n\in {\mathbb{R}}^{|\mathcal{I}|}$ be the vector of molecular counts with elements n_i. If S is the stoichiometric matrix of elements S_i,j = γ_i,j where γ_i,j are the signed stoichiometric coefficients such that ν_i,j = |γ_i,j|, we have according to the law of mass action,

$\frac{dn}{dt}=SJ\left(n\right)(2)$

where for a given n, J(n) is the net flux of the set of forward and reverse reactions. The chemical species occurring on the left hand side of the equation are known as reactants and belong to the subset ${\mathcal{I}}_{R_j}\subset \mathcal{I}$, while those on the right hand side are known as products and belong to the subset ${\mathcal{I}}_{P_j}\subset \mathcal{I}$.

For any reversible reaction with forward and reverse reactions + j and −j, the net flux is given by,

$J_j\left(n\right)=k_j\prod _{i\in {\mathcal{I}}_{R_j}}{n}_i^{{\nu }_{i,j}}-k_{-j}\prod _{i\in {\mathcal{I}}_{R_{-j}}}{n}_i^{{\nu }_{i,-j}}.(3)$

Eq. 3 is a purely kinetic description of the reaction flux. For elementary reactions, thermodynamic terms are introduced into the law of mass action by a simple factorization

$J_j\left(n\right)=k_{-j}\prod _{i\in {\mathcal{I}}_{R_{-j}}}{n}_i^{{\nu }_{i,-j}}\left(\frac{k_j{\prod }_{i\in {\mathcal{I}}_{R_j}}{n}_i^{{\nu }_{i,j}}}{k_{-j}{\prod }_{i\in {\mathcal{I}}_{R_{-j}}}{n}_i^{{\nu }_{i,-j}}}\right)-k_j\prod _{i\in {\mathcal{I}}_{R_j}}{n}_i^{{\nu }_{i,j}}\left(\frac{k_{-j}{\prod }_{i\in {\mathcal{I}}_{R_{-j}}}{n}_i^{{\nu }_{i,-j}}}{k_j{\prod }_{i\in {\mathcal{I}}_{R_j}}{n}_i^{{\nu }_{i,j}}}\right)(4)$

$=k_{-j}\prod _{i\in {\mathcal{I}}_{R_{-j}}}{n}_i^{{\nu }_{i,-j}}\left(K_j\prod _i{n}_i^{{\gamma }_{i,j}}\right)-k_j\prod _{i\in {\mathcal{I}}_{R_j}}{n}_i^{{\nu }_{i,j}}\left(K_{-j}\prod _i{n}_i^{{\gamma }_{i,-j}}\right),(5)$

where $K_j=\frac{k_j}{k_{-j}}$ is the equilibrium constant for reaction j. In Eq. 5, the terms in parentheses are related to the thermodynamic forces on the reaction while the terms outside the parentheses are the time-dependent kinetic terms. Setting the latter terms to a constant gives the Marcelin Equation (Marcelin, 1910), in which each forward and reverse reaction occurs on the same timescale by postulate. While this postulate is incorrect for determining the true dynamics of the system, the Marcelin formulation of Eq. 5 does give a way to the maximum path entropy solution (Cannon et al., 2018; Britton et al., 2020) by assuming that all reactions occur on the same timescale and that their fluxes are proportional to the thermodynamic forces on them. Once these fluxes and concentrations are known from the maximum path entropy solution, they can be used to back-calculate the rate constants that will give the corresponding true dynamics (Cannon et al., 2018; Britton et al., 2020).

Using the maximum path entropy solution, and representing the thermodynamic reaction forces for each reaction $j\in \mathcal{J}$ as,

$f_j\left(n\right)=K_j\prod _{i\in \mathcal{I}}{n}_i^{-{\gamma }_{i,j}},(6)$

the flux is then given by

$J_j\left(n\right)=f_j\left(n\right)-\frac{1}{f_j\left(n\right)}(7)$

and we can express the time dependence of each metabolite $i\in \mathcal{I}$ as follows

$\frac{d{n}_i}{dt}=\sum _{j\in \mathcal{J}}{\gamma }_{i,j}\left(f_j\left(n\right)-\frac{1}{f_j\left(n\right)}\right).(8)$

2.2 Steady state

At the maximum path entropy configuration the system is in a steady state, where a set of metabolites, denoted ${\mathcal{I}}_f$, are assumed to be held at fixed concentrations as boundary conditions for the system. For the rest of the metabolites, ${\mathcal{I}}_v=\mathcal{I}\backslash {\mathcal{I}}_f$, their steady state concentration is free and their rate of change is zero. For all $i\in {\mathcal{I}}_f$, let ${\bar{n}}_i$ be the fixed concentration of the metabolite. A steady state is defined as a solution to the following system of equations,

$\begin{array}{cc}\hfill & \frac{d{n}_i}{dt}=0\forall i\in {\mathcal{I}}_v,\hfill \\ \hfill & n_i={\bar{n}}_i\forall i\in {\mathcal{I}}_f.\hfill \end{array}(9)$

2.3 Controlling metabolite concentrations

Metabolite concentrations predicted using the maximum path entropy approach without regulation will produce values that are much too large to be physiologically reasonable—concentrations may approach the limit of their solubility, causing the cytoplasm to become glass-like (Parry et al., 2014; Heimlicher et al., 2019).

That regulation is needed to control concentrations in vivo was proposed early in the field of enzymology (Atkinson, 1969;Atkinson, 1977). In a previously reported method to control concentrations (Britton et al., 2020), Britton, et al. applied regulation to reactions using a scalar valued activity coefficient α_j ∈ [0.0, 1.0] for each reaction j, that linearly scales the reaction rate such that the time dependence of metabolite concentration n_i is given by,

$\frac{d{n}_i}{dt}=\sum _{j\in \mathcal{J}}{\gamma }_{i,j}{J}_j\left(n,{\alpha }_j\right),(10)$

where

$J_j\left(n,{\alpha }_j\right)={\alpha }_j\left(f_j\left(n\right)-\frac{1}{f_j\left(n\right)}\right).(11)$

When α_j = 1.0, the reaction is fully active and when α_j = 0.0, the activity, and hence reaction rate J_j (n, α_j), is zero. Activities are adjusted in a deterministic manner by characterizing the sensitivity of a metabolite concentration to an enzyme activity using Metabolic Control Analysis (MCA) (Sauro, 2018). In MCA, the concentration control coefficient C_i,j measures the sensitivity of a concentration n_j to an activity coefficient α_j,

$C_{i,j}=\frac{\partial \log n_i}{\partial \log {\alpha }_j}.(12)$

Using this approach, metabolites whose predicted concentrations are furthest from their experimentally observed values, as measured by the log ratio of the predicted to observed concentration, are reduced first. The concentration of a metabolite n_i is reduced by adjusting the activity coefficient α_j that has the largest influence on n_i and all other metabolites that exceed their experimentally observed concentrations, as determined by the sensitivity analysis. This process is iteratively carried out until all concentrations are at or below the experimentally observed concentrations. Details of the approach are provided in the study by Britton et al. (2020). In the cases discussed herein, the maximum value of the experimentally observed metabolite concentrations are taken to be n_i = 1.0 mM for $i\in {\mathcal{I}}_v$ and for $i\in {\mathcal{I}}_f$ the values are taken from mass spectrometry observed concentrations for bacteria (Bennett et al., 2009; Park et al., 2016).

2.4 Controlling metabolite concentrations and maximizing growth

While controlling metabolite concentrations may be a primary role of metabolic regulation, natural selection also requires that organisms be regulated to grow fast and efficiently—using the available energy from the environment to ensure survival and compete with others. Here we develop a novel approach, we call pathway-controlled optimization (PCO), whose goal is to obtain this biological objective, and derive a formulation that allows for tractable numerical solutions.

Let $\mathcal{G}\subset \mathcal{J}$ be the set of reactions corresponding to production of biomass. We formulate the steady state with maximum biomass production as the solution of the optimization problem.

$\max \sum _{j\in \mathcal{G}}{J}_j\left(n,{\alpha }_j\right)(\mathrm{13a})$

$\begin{array}{ll}\hfill & \text{subject\hspace{0.17em}to:}\hfill \\ \hfill & \frac{d{n}_i}{dt}=0\forall i\in {\mathcal{I}}_v,\hfill \end{array}(\mathrm{13b})$

$n_i={\bar{n}}_i\forall i\in {\mathcal{I}}_f,(\mathrm{13c})$

$0\le n_i\le n_{\max }\forall i\in {\mathcal{I}}_v,(\mathrm{13d})$

$0\le {\alpha }_j\le 1\forall j\in \mathcal{J}.(13e)$

The objective seeks to maximize the flux through the growth reactions $\mathcal{G}$ while the constraints (13b) - (13c) ensure that the steady state (9) is satisfied. The activity coefficients and metabolite concentrations are further restricted to physiologically meaningful values with the constraints (13e) and (13d).

The formulation of the PCO problem is simple to express but difficult to solve. The steady state constraints (13b) are non-linear and non-convex presenting significant challenges to optimization. Values for the flux, activity coefficients, and metabolite concentrations can also vary over many orders of magnitude, which introduces additional difficulty in employing numerical methods to compute solutions. In this work we present a more computationally tractable reformulation of the constraints and present numerical solutions from an interior point solver.

2.4.1 Representing the steady state condition

To simplify solving the optimization, we reformulate the steady state constraint (13b) to be more numerically tractable. Let $S_v\in {\mathbb{R}}^{|{\mathcal{I}}_v|,|\mathcal{J}|}$ be the sub-matrix of S with the rows corresponding to only the variable metabolites. Then at a steady state the reaction fluxes must satisfy

$S_{v}J\left(n,\alpha \right)=0,(14)$

that is $J(n,\alpha )\in \mathcal{N}(S_v)$, where $\mathcal{N}(S_v)$ is the nullspace of S_v defined

$\mathcal{N}\left(S_v\right)=\left\{x\in {\mathbb{R}}^{|\mathcal{J}|}:S_{v}x=0\right\}.(15)$

We separate (14) into the identification of a vector of fluxes $y\in {\mathbb{R}}^{|\mathcal{J}|}$ which satisfy steady state and the construction of a vector of metabolites $n\in {\mathbb{R}}^{|\mathcal{I}|}$ which achieve those fluxes with the two conditions,

$\begin{array}{cc}\hfill y& \in \mathcal{N}\left(S_v\right),\hfill \\ \hfill y_j& =J_j\left(n,{\alpha }_j\right)\forall j\in \mathcal{J}.\hfill \end{array}(16)$

The expression y_j = J_j (n, α_j) can be made more tractable by considering the computation of the reaction flux in terms of the log of the metabolites. For $\eta \in {\mathbb{R}}^{|\mathcal{I}|}$ let

${\eta }_i=\log \left(n_i\right)\forall i\in \mathcal{I},(17)$

then we have from (7) that the flux as a function of η in the maximum path entropy formulation is given by

$J_j\left(\eta ,{\alpha }_j\right)={\alpha }_j{K}_j{e}^{-\langle {\left(S\right)}_j,\eta \rangle }-{\alpha }_j\frac{1}{K_j}{e}^{\langle {\left(S\right)}_j,\eta \rangle },(18)$

where (S)_j is the jth column of S, and ⟨⋅, ⋅⟩ is the standard inner product. It follows that in terms of the log of the metabolite counts, conditions (16) are equivalent to.

$y\in \mathcal{N}\left(S_v\right),(\mathrm{19a})$

$S^T\eta =\widehat{y}\left(y,\alpha \right),(\mathrm{19b})$

${\widehat{y}}_j\left(y,\alpha \right)=\log \left(\frac{K_j}{2{\alpha }_j}\left(-y_j+\sqrt{y_j^2+4{\alpha }_j^2}\right)\right)\forall j\in \mathcal{J}.(\mathrm{19c})$

Computation of $\widehat{y}(y,\alpha )$ is prone to error. For example, in the instance where y_j is positive, rounding may result in evaluating the log of zero. To avoid these issues we use the following equivalent expression to (19c)

${\widehat{y}}_j\left(y,\alpha \right)=\log \left(K_j\right)+\mathrm{sgn}\left(y_j\right)\left(\log \left(2{\alpha }_j\right)-\log \left(|y_j|+\sqrt{y_j^2+4{\alpha }_j^2}\right)\right),(20)$

with sgn the signum function given by

$\mathrm{sgn}\left(x\right)=\left\{\begin{array}{cc}1\hfill & \hspace{0.17em}x>0\hfill \\ 0\hfill & \hspace{0.17em}x=0\hfill \\ -1\hfill & \hspace{0.17em}x<0\hfill \end{array}\right..(21)$

A further complication in using the steady state conditions (19) as constraints for optimization is that the value of $\widehat{y}(y,\alpha )$ can become extremely sensitive to small changes in the flux. In particular the partial derivative

$\frac{\partial {\widehat{y}}_j}{\partial y_j}=\frac{1}{\sqrt{y_j^2+4{\alpha }_j^2}}(22)$

grows unbounded as α_j and y_j go to zero. In this application we consistently found solutions with α_j and y_j on the order of 10^–9 for some reactions while others were unregulated with values of y_j on the order of 10³. These discrepancies in scale are challenging for numerical methods.

However, it is not necessary to explicitly compute the activity coefficients in order to characterize a steady state. Instead it is sufficient only to identify metabolite concentrations that produce flux values that can be reduced to a steady state. We can capture this criteria with the conditions

$\begin{array}{cc}\hfill & |y_j|\le |J_j\left(n,\mathbf{1}\right)|\forall j\in \mathcal{J},\hfill \\ \hfill & y_j{J}_j\left(n,\mathbf{1}\right)>0,\hfill \end{array}(23)$

where 1 is the vector of all ones corresponding to no regulation for all reactions. If (23) is satisfied then any reaction fluxes J_j (n, 1) with magnitude greater than the steady state value y_j can always be reduced to satisfy equality so long as they have the same sign, which is captured by the second condition in (23). This idea can also be captured in our modified formulation (19) for the steady state conditions.

For metabolites η and steady state fluxes y there exists an α such that (19) holds if the following hold

$\begin{array}{cc}\hfill & |g_j-\log \left(K_j\right)|\ge -\left(\log \left(2\right)-\log \left(|y_j|+\sqrt{y_j^2+4}\right)\right)\forall j\in \mathcal{J},\hfill \\ \hfill & \left(\log \left(K_j\right)-g_j\right)y_j\ge 0\forall j\in \mathcal{J},\hfill \end{array}(24)$

where g = S^Tη. This follows from (19b)-(19c), (20), and the inequality

$\log \left(2{\alpha }_j\right)-\log \left(|y_j|+\sqrt{y_j^2+4{\alpha }_j^2}\right)\le \log \left(2\right)-\log \left(|y_j|+\sqrt{y_j^2+4}\right),(25)$

for all $y_j\in \mathbb{R}$ and α_j ∈ [0, 1].

Therefore we can formulate steady state conditions equivalent to (19) without the need to explicitly specify the activity coefficients as follows.

$y\in \mathcal{N}\left(S_v\right),(\mathrm{26a})$

$g=S^T\eta ,(\mathrm{26b})$

$h_j=\mathrm{sgn}\left(y_j\right)\left(\log \left(2\right)-\log \left(|y_j|+\sqrt{y_j^2+4}\right)\right)\forall j\in \mathcal{J},(\mathrm{26c})$

$|g_j-\log \left(K_j\right)|\ge |h_j|\forall j\in \mathcal{J},(\mathrm{26d})$

$\left(\log \left(K_j\right)-g_j\right)y_j\ge 0\forall j\in \mathcal{J}.(26e)$

2.4.2 Numerical formulation of the steady state constraints

Implementing conditions (26a) - (26e) as constraints for numerical optimization requires further reformulation. The constraint (26d) is non-convex, and represents an ‘or’ condition depending on the sign of the flux for each reaction. This switching condition makes the feasible set non-convex and therefore challenging for optimization methods to search over. To find solutions we use a big M relaxation of (26d) with the constraints

$h_j=\mathrm{sgn}\left(y_j\right)\left(\log \left(2\right)-\log \left(|y_j|+\sqrt{y_j^2+4}\right)\right)\forall j\in \mathcal{J},(\mathrm{27a})$

$g_j-\log \left(K_j\right)\ge h_j-u_{j}M\forall j\in \mathcal{J},(\mathrm{27b})$

$g_j-\log \left(K_j\right)\le h_j+\left(1-u_j\right)M\forall j\in \mathcal{J},(\mathrm{27c})$

$2u_j-1=\mathrm{sgn}\left(y_j\right)\forall j\in \mathcal{J},(\mathrm{27d})$

where M is taken to be a large constant. For each reaction the variable u_j is a “switching” term that relaxes either the constraint (27b) or (27c) depending on the sign of y_j such that only the correct constraint is imposed for each reaction.

In several constraints we utilize the signum function of the flux which is a discrete function, whereas the optimization method we use here requires continuous functions of all variables for the constraints. Therefore we approximate the signum function with

$s\tilde{g}n\left(x\right)=\frac{\lambda x}{\lambda |x|+ϵ}(28)$

for scaling parameters λ and ϵ.

To implement the constraint $y\in \mathcal{N}(S_v)$ we construct a basis for the nullspace. For m the dimension of $\mathcal{N}(S_v)$, let $B\in {\mathbb{R}}^{|\mathcal{J}|,m}$ be a matrix with columns given by a basis of $\mathcal{N}(S_v)$. Then we have that

$\mathcal{N}\left(S_v\right)=\left\{B\beta :\beta \in {\mathbb{R}}^m\right\},(29)$

and we capture the constraint $y\in \mathcal{N}(S_v)$ with the expression y = Bβ for a $\beta \in {\mathbb{R}}^m$.

2.4.3 Maximum growth optimization numerical formulation

Using the steady state conditions as outlined in sections 2.4.2 and 2.4.1 we formulate the maximum growth pathway-controlled optimization problem as

$\underset{y,\eta }{\max }\sum _{j\in \mathcal{G}}{y}_j(\mathrm{30a})$

$\begin{array}{ll}\text{subject\hspace{0.17em}to:}\hfill & \hfill \\ y=B\beta ,\hfill \end{array}(\mathrm{30b})$

$g=S^T\eta ,(\mathrm{30c})$

$h_j=s\tilde{g}n\left(y_j\right)\left(\log \left(2\right)-\log \left(|y_j|+\sqrt{y_j^2+4}\right)\right)\forall j\in \mathcal{J},(\mathrm{30d})$

$g_j-\log \left(K_j\right)\ge h_j-uM\forall j\in \mathcal{J},(30e)$

$g_j-\log \left(K_j\right)\le h_j+\left(1-u\right)M\forall j\in \mathcal{J},(\mathrm{30f})$

$\left(\log \left(K_j\right)-g_j\right)y_j\ge 0\forall j\in \mathcal{J},(\mathrm{30g})$

$2u_j-1=s\tilde{g}n\left(y_j\right)\forall j\in \mathcal{J},(\mathrm{30h})$

${\eta }_i\le {\eta }_{{\max }_i}\forall i\in {\mathcal{I}}_v,(\mathrm{30i})$

${\eta }_i={\bar{\eta }}_i\forall i\in {\mathcal{I}}_f,(\mathrm{30j})$

where the constraints (30b)-(30h) represent the steady state constraint (13b). Note that the activity coefficients α for reactions are not explicitly included, however for a solution with optimal log metabolite counts η^* and fluxes y^*, the associated activity coefficients can be recovered as

${\alpha }_j^{\ast }=\frac{y_j^{\ast }}{J_j\left({\eta }^{\ast }\right)}\forall j\in \mathcal{J}.(31)$

It is important here to note that when (30) admits solutions they are not necessarily unique. For instance, those reactions that do not directly impact the objective may have a range of choices for the fluxes and constituent metabolite values that satisfy constraints (30b)–(30h) without changing the objective value. In particular constraints (30e) and (30f) are not neccesairly all that restrictive.

2.4.4 Initialization for optimization

The greatest difficulty in solving (30) is in identifying the correct sign for the fluxes. If the correct flux directions are known then the problem simplifies considerably. Using IPOPT, a general non-linear interior point solver (Wächter and Biegler, 2006), we found that the solution and convergence of the solver was extremely sensitive to the variable initialization, as is common for non-convex problems. In particular, starting with fluxes close to the optimal directions was crucial for convergence to a good solution. To initialize the variables we solve for metabolite values that produces fluxes close to the gradient direction of the objective. Using random or other initialization methods we found that the solver failed to converge, or converged to points with smaller objective values.

To compute initial variable values we first find the projection of the gradient of the objective onto the space of steady state solutions, with

$y_g={\mathbb{P}}_{\mathcal{N}\left(S_v\right)}\left({\nabla }_y\left[\sum _{j\in \mathcal{G}}{y}_j\right]\right)(32)$

where ${\mathbb{P}}_{\mathcal{N}(S_v)}$ is the projection operator onto $\mathcal{N}(S_v)$. We then solve the constrained linear least squares problem

$\underset{\eta }{\min }{\left|\left|\mathrm{sgn}\left(-y_g\right)◦\left(S^T\eta -\log \left(K\right)\right)-s\right|\right|}^2(\mathrm{33a})$

$\begin{array}{cc}\text{subject\hspace{0.17em}to:}\hfill & \hfill \\ s_j\ge 0,\forall j\in \mathcal{J},\hfill \end{array}(\mathrm{33b})$

${\eta }_i\le {\eta }_{{\max }_i}\forall i\in {\mathcal{I}}_v,(\mathrm{33c})$

${\eta }_i={\bar{\eta }}_i\forall i\in {\mathcal{I}}_f,(\mathrm{33d})$

where s is a vector of positive slack variables and ◦ is the Hadamard product. Consider that if the objective is zero at a solution η∗ then sgn (−y_g) = sgn (S^Tη∗ − log(K)) which implies that there exists a γ > 0 such that (η∗, γy_g) is a feasible point of (30). Therefore this is a feasible point with fluxes in the directions of y_g. In practice solutions of (33) may only be close to zero, so we choose initial fluxes y∗ as

$y^{\ast }=-\zeta {\mathbb{P}}_{\mathcal{N}\left(S_v\right)}\left(S^T{\eta }^{\ast }-\log \left(K\right)\right)(34)$

for a ζ > 0. Larger zeta values move the fluxes further from the origin and was found to improve convergence, likely as it reduces the number of fluxes that will switch sign early in the optimization.

2.5 Numerical solution

We solved (33) with the SciPy version 1.7.1 lsq_linear routine using the default configuration (Virtanen et al., 2020). A solution of (30) was then computed with IPOPT version 3.14.4 using Pyomo (Hart et al., 2017) again using the default optimization parameters except with tol = 1.0 × 10⁻⁷ and maximum iterations set to 10,000.

The parameters for $\widetilde{\mathrm{sgn}}$ in constraint (30d) were λ = 1.0 × 10⁵⁰ and ϵ = 1.0 × 10⁻⁵⁰, while in constraint (30h) they were λ = 1 and ϵ = 1 × 10⁻⁵⁰. The big M value was set to 100, and we used ζ = 10 for the initialization computation.

For all variable metabolites the upper bound was ${\eta }_{{\max }_i}=13.308$. Values for the fixed metabolite concentrations, equilibrium constants, and stochiometeric matrix are given in the Supplementary Material. Equilibrium constants were computed using eQuilibrator version 0.3.1 (Beber et al., 2021). Values for the fixed metabolites were taken from mass spectrometry observed concentrations for bacteria (Bennett et al., 2009; Park et al., 2016).

2.6 Metabolic model

The metabolic model of R. rubrum consisted of 184 reactions and 204 metabolites including central metabolism, the Calvin cycle, the oxidative and reductive TCA cycles, biosynthetic pathways for all amino acids, pyrimidines and purines, the S-adenylsyl methionine (SAM) cycles I and II, the methyl alkane reductase pathway (North et al., 2020), sulfate assimilation, and generic RNA, DNA and protein synthesis pathways. The model represents an engineered species that is predicted to not need the ethylmalonyl CoA pathway due to the inclusion of a oxaloacetate decarboxylase/malic enzyme that converts oxaloacetate to pyruvate and CO₂ and the inactivation of phospho-glycerate mutase. Consequently, all 3-phophoglycerate produced by the Calvin Cycle is stays within the Calvin cycle.

Boundary conditions (fixed species concentrations) consisted of concentrations of nutrients (ethanol, diphosphate, orthophosphate, sulfate and ammonia), cell waste products (CO₂, H₂, ethylene), ATP derivatives to set up the energy gradient due to photoheterotrophic growth (ATP, ADP, AMP), internal redox pairs (NAD+, NADH, NADP+, NADPH, reduced and oxidized ferredoxin, reduced and oxidized thioredoxin), folates (N10-formyltetrahydrofolate, 5-methyltetrahydrofolate, 5,10-methylenetetrahydrofolate, 7,8-dihydrofolate and tetrahydrofolate), and other internal metabolites for which scavenging pathways were not modeled (coenzyme A, (4S)-4,5-dihydroxypentan-2,3-dione, adenine, adenosine, adenosine-3’,5’-bisphosphate).

The standard free energy of reaction for each reaction was calculated with eQuilibrator, version 0.3.1. The free energies for the methyl alkane reductase pathway, the SAM cycles I & II, the ferredoxin oxidoreductases for 2-oxoglutarate and pyruvate, DNA, RNA and protein polymerization reactions were not available from eQuilibrator and were estimated manually from half-reactions reactions or overall reactions for the respective pathways. If estimates for some individual reactions were available for these pathways, they were adjusted accordingly, otherwise the standard free energy change for the pathways was divided equally among the constituent reactions. Details are included in the computational notebook available as supplemental material.

The model was constructed from a custom R. rubrum BioCyc database at PNNL. The draft model is available as a Jupyter notebook along with the necessary auxillary files for running the notebook in King and Cannon (2023).

3 Results

In order to place the results in the proper context, a brief overview of the metabolism of R. rubrum, the subject of the biomass optimization, is necessary. R. rubrum is capable of growing both photoautotrophically, meaning that it only uses CO₂ as a carbon source, or photoheterotrophically, meaning that it can carry out photosynthesis while simultaneously assimilating organic carbon. In this study, we are modeling photoheterotropic growth, in which energy is acquired through photosynthesis while carbon is mostly acquired by uptake and assimilation of ethanol.

During photosynthesis, R. rubrum uses the Calvin cycle to regenerate chemical precursors for the assimilation of CO₂. The Calvin cycle includes reactions shared with the non-oxidative branch of the pentose phosphate cycle, some of which can act as an alternate route in the Calvin Cycle. The two routes in the Calvin cycle differ in that one route produces erythrose-4-phosphate and the other route consumes erythrose-4-phosphate. Erythrose-4-phosphate is a key metabolite for the synthesis of vitamin B6, pyridoxal-5’-phosphate, an essential cofactor in amino acid, and hence protein, synthesis. Therefore, it is important that during growth, erythrose-4-phosphate not be entirely used as a Calvin cycle precursor.

We define the biomass pathways to be optimized as the pathways regarding the incorporation of amino acids into proteins and the incorporation of nucleotide triphosphates (NTPs) into RNA and DNA. The synthesis of amino acids and NTPs from precursors, however, are not included in the biomass pathways. During growth, proteins, RNA and DNA can be simultaneously synthesized, although RNA and DNA synthesis compete for both available energy and molecular precursors, NTPs and deoxynucleotide triphosphates (dNTPs), respectively. R. rubrum is known to have a cell cycle in which DNA can only be synthesized during the S phase, while RNA is synthesized continuously throughout the cell cycle except during mitosis. Whether DNA is synthesized depends on the redox state of the cell, specifically, the ratio of NADP to NADPH.

3.1 Controlling concentrations

Because some metabolites have highly favorable free energies of formation or standard chemical potentials, their concentrations can rise to dangerously high levels if not controlled (Atkinson, 1977). The primary role of control in a biological system must be to maintain viability. Metabolism will not operate if the concentrations of metabolites becomes so great that the cytoplasm becomes too viscous for diffusion-reaction processes to occur. We previously showed that concentrations of metabolites can be controlled by reducing key reaction fluxes through the activities of the respective enzymes identified by concentration control coefficients from Metabolic Control Analysis. In this approach, described briefly in the Methods section and in detail in reference (Britton et al., 2020), sensitivity analyses are used to characterize the influence of an enzyme’s activity on the concentration of any of the models metabolites. If a metabolite concentration is too high, the enzyme with the most influence over that metabolite is chosen to be regulated. This process occurs until all metabolite concentrations are at physiological levels. The approach can be implemented either as an optimal control problem or using reinforcement learning. Analysis of central metabolism showed that both analyses produced similar results and identified known regulators of central metabolism.

In Table 1, we compare the regulation of enzymes that were identified as necessary to be controlled to maintain physiological levels of metabolites in the MCA predictions with those from the pathway-controlled optimization (PCO) predictions. Of course, in the PCO method, reactions are controlled both to maintain physiological levels and to maximize flux to the growth reactions, so any amount of regulation predicted by the PCO method may result from a combination of the two objectives of maximizing growth while constraining metabolite concentrations to physiological levels. Nevertheless, as shown in Table 1, the PCO predicted regulation of reactions is consistent with the MCA predicted regulation of reactions.

TABLE 1

TABLE 1. Reactions regulated in the MCA approach in which the goal was to control concentrations. Ten out of twelve of these reactions are also regulated in the pathway-controlled optimization approach.

The PCO method regulated ten out of twelve of the same reactions as the MCA method, suggesting that the PCO method reduced the activity of these reactions in order to maintain physiological levels of the metabolites. However, the observation that the flux in some of these reactions is now effectively zero may indicate that this subset of reactions were additionally shutdown to redirect flux to the growth reactions.

In addition, the PCO method may need to control more reactions than the MCA method, since reaction fluxes are also redistributed due to optimizing the growth reaction rates (discussed below). Decreased or increased reaction flux may result in higher reactant or product levels, respectively, such that the increased levels need to be controlled to keep them at physiological levels. The PCO approach regulated reactions that are not regulated in MCA and yet have non-zero fluxes. These are shown in Table 2. Presumably, these additional reactions need to be regulated to maintain metabolite concentrations at physiological levels in the new regime.

TABLE 2

TABLE 2. Reactions regulated in the PCO method but not in the MCA method and that have significant reaction flux in the PCO method. Presumably these reactions are regulated to control metabolite concentrations.

3.2 Increasing growth

In order to increase flux through the growth reactions in accordance with the PCO objective function, additional reactions are regulated to redirect flux to that end. The cumulative value of all the reaction fluxes for both the MCA-regulated reactions and the PCO-regulated reactions are shown in Figure 1. Some reactions in the PCO model are indeed reduced and some are increased, with the overall effect being a large increase in absolute rates of a subset of reactions. In general, more flux is directed to biosynthetic reactions, with the net result being a redirection of flux away from central metabolism in comparing the MCA-derived flux distribution to the PCO-derived flux distribution, as indicated in the flux maps for the two models in Figure 2.

FIGURE 1

FIGURE 1. Comparison of the cumulative distribution of the reaction fluxes for the MCA model and the PCO model. Flux through the growth reactions are maximized by turning reaction fluxes down but not off.

FIGURE 2

FIGURE 2. Comparison of the fluxes through central metabolism for the MCA model (top) and the PCO model (bottom). Flux through the growth reactions are maximized by turning reaction fluxes down but not off. The decreased flux through central metabolism in the PCO model reflects the same decrease shown below the red line for the MCA model in Figure 1.

The reactions with the highest fluxes from the PCO approach are shown in Table 3. The reaction with the highest absolute flux in the PCO regulated set of reactions is the RNA synthesis reaction, which increased to 8.11 × 10³, a more than 7-fold increase as compared to the MCA regulated model for controlling only metabolite levels. Likewise, the reaction for protein synthesis is included in this group. Both of these reactions are part of the growth reaction set used in the objective function. The other ten reactions with the highest flux are all involved in uptake and processing of the carbon source, ethanol. Much of the carbon from ethanol flows into fumarate and then into the pathways for synthesis of the pyrimidine nucleosides, uridine triphosphate (UTP) and cytidine triphosphate (CTP), precursors of RNA synthesis. In contrast the purine nucleosides are readily available from the already high levels of ATP.

TABLE 3

TABLE 3. Top 12 reactions from the constrained optimization that have the highest fluxes. The reaction that has the highest flux in the constrained optimization is the growth reaction for RNA synthesis. Likewise, the reaction for protein synthesis is included in this group. Other reactions with the highest flux are all involved in uptake and processing of the carbon source, ethanol.

Reactions that were effectively shutdown are shown in Table 4. As can be seen at the bottom of the table, DNA synthesis was minimized by significantly down regulating the synthesis of dNTPs, making more NTPs available for RNA synthesis. In addition, the alternate pathway in the Calvin cycle that consumes erythrose-4-phosphate is shutdown (top three reactions in Table 4), making more erythrose-4-phosphate available for synthesis of vitamin B6, a necessary component for production of the branched chain amino acids tryptophan and tyrosine.

TABLE 4

TABLE 4. Reactions regulated in the PCO method but not in the MCA method and that have insignificant reaction flux in the PCO method. Presumably these reactions are shut down to redirect mass flow to the growth reactions.

4 Discussion

4.1 Biological discussion

The PCO method we present here agrees with the majority of the regulation imposed by the MCA method for maintaining solvent capacity of the cell, which was shown previously to agree with known points of regulation from experiments as reported in the literature Britton et al. (2020). However, the consideration of growth as an additional objective for regulation in our PCO approach produced significant additional regulation in line with adjusting flux through growth pathways. The inclusion of the consideration of growth is likely to be crucial for a complete understanding of regulation in cells.

While this study uses a limited model (184 reactions and 204 metabolites), that DNA synthesis was turned off, as it is during most of the cell cycle, is interesting, and presumably was done so in order to maximize RNA synthesis. Further evaluation with a more complete model is needed. It is not clear if under low NADP/NADPH conditions that DNA synthesis is likewise maximized by reducing RNA synthesis.

Of course, it will be critical to test the results of the model against experimental observations in order to draw firm biological conclusions. In this regard, previous studies of the metabolic dynamics of R. rubrum using Metabolic Flux Analysis (McCully et al., 2020), an experimental method to derive fluxes from observations of isotope distributions, will provide a good test of the regulation model.

If successful, the incorporation of thermodynamics into the optimization of growth with the PCO approach developed here, along with the ability to infer kinetic parameters from maximum path entropy distributions, will provide a way to further explore the relationship between natural selection and thermodynamics in detail. Ideally, these capabilities will allow the study of the expression patterns and activities of enzymes at various stages in the cell cycle.

4.2 Optimization discussion

One of the biggest challenges in solving the PCO problem is the range of the variable values. For the solutions presented here both metabolites concentrations and activity coefficients for reactions vary over many orders of magnitude. To make the problem more numerically tractable our formulation in essence seeks to do computations with respect to the log of the variables to mitigate the large discrepancies in scale. This is relatively straightforward with respect to the metabolite concentrations in part because they must be positive. However, to handle the flux we must separate out its magnitude and sign which complicates the formulation. It is important to note that our approach does not strictly adhere to the requirements of the IPOPT solver, namely, that the objective and constraints be twice continuously differentiable. In particular, both constraint (27a) and (30h) violate this condition.

For the constraint (27a), the original expression (19c) it is derived from is twice continuously differentiable but resulted in numerical errors. By introducing the form (20) and using the approximation (28) to the signum function a cusp is introduced into the first derivative at zero and error is introduced in an interval about the origin determined by the parameters λ and ϵ. We choose λ and ϵ to shrink this interval such that the true derivative is closely approximated at all points evaluated by the solver. We found this choice to produce more accurate solutions and better convergence behavior than using twice continuously differentiable approximations of the signum function or other forms for (19c). It may be possible to further improve performance by supplying the solver directly with alternative derivative computations rather than relying on automatic differentiation, particularly if solutions approach the origin.

The relaxed formulation (30) also uses the approximation (28) to the signum function in the constraint (30h) specifying the value of the ‘switching’ variable u. Here choices of λ and ϵ are made such that the right hand side of the constraint is effectively constant and partials with respect to the flux will evaluate to zero, in order to hold u constant at 1 or −1 so long as the sign of the flux remains the same. This formulation captures the problem well when the sign on the fluxes are in the optimal direction but it is unclear to what extent it constrains the domain reachable by the solver from an arbitrary initial condition. We also evaluated solutions including an explicit computation of the activity coefficients as in (20) but found that convergence of the solver required setting a lower bound on regulation of roughly 1 × 10⁻⁹, and solutions had a significantly lower objective value than were found using the implicit regulation formulation (30).

Due to the potential restrictions on the solutions that can be obtained with our approach using the relaxed formulation (30) we will consider methods to more thoroughly explore the feasible set in future work including stochastic search methods. Of course a more complete exploration of the space of solutions can come at significant computational cost and if other methods perform similarly to the approach we present here, then it may have a clear practical advantage in that the times-to-solution will likely be much faster. Additionally, as discussed in Section 2.4.3 solutions are not necessarily unique. Methods to explore and characterize the set of equivalent solutions, that is solutions that can achieve the same objective value, will allow for a greater understanding of the flexibility that may be possible in controlling reaction fluxes and metabolite concentrations while achieving high growth rates.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

WC conceived and designed the maximum entropy approach to model metabolism and developed the R. rubrum metabolic model. EK designed and implemented the optimization method for learning regulation to maximize growth. JH contributed to the optimization development. JN contributed to the development of the metabolic model of Rhodospirilium rubrum. The manuscript was written by EK and WC.

Funding

WC and JN were supported by the U.S. Department of Energy, Office of Biological and Environmental Research through contracts 78266 and award DE-SC0022091, respectively. EK and JH were supported by the Data Model Convergence Initiative, and EK additionally by the Predictive Phenomics Initiative, under the Laboratory Directed Research and Development Program at the Pacific Northwest National Laboratory. PNNL is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy under contract DE-AC05-76RLO 1830.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsysb.2023.981866/full#supplementary-material

Supplementary Data Sheet S1 | Equilibrium constants.

Supplementary Data Sheet S2 | Metabolite concentrations.

Supplementary Data Sheet S3 | Stoichiometric Matrix.

References

Aithal, H. N., Walsh-Reitz, M. M., and Toback, F. G. (1985). Regulation of glyceraldehyde-3-phosphate dehydrogenase by a cytosolic protein. Am. J. Physiol. 249, C111–C116. doi:10.1152/ajpcell.1985.249.1.C111

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkinson, D. E. (1977). Cellular energy metabolism and its regulation. New York: Academic Press.

Google Scholar

Atkinson, D. E. (1969). “Limitation of metabolite concentrations and the conservation of solvent capacity in the living cell,” in Current topics in cellular regulation. Editors B. L. Horecker, and E. R. Stadtman (Academic Press), 1, 29–43. doi:10.1016/B978-0-12-152801-0.50007-9

CrossRef Full Text | Google Scholar

Beber, M. E., Gollub, M. G., Mozaffari, D., Shebek, K. M., Flamholz, A., Milo, R., et al. (2021). eQuilibrator 3.0: a database solution for thermodynamic constant estimation. Nucleic Acids Res. 50, D603–D609. doi:10.1093/nar/gkab1106

CrossRef Full Text | Google Scholar

Bennett, B. D., Kimball, E. H., Gao, M., Osterhout, R., Van Dien, S. J., and Rabinowitz, J. D. (2009). Absolute metabolite concentrations and implied enzyme active site occupancy in escherichia coli. Nat. Chem. Biol. 5, 593–599. doi:10.1038/nchembio.186

PubMed Abstract | CrossRef Full Text | Google Scholar

Britton, S., Alber, M., and Cannon, W. R. (2020). Enzyme activities predicted by metabolite concentrations and solvent capacity in the cell. J. Roy. Soc. Interfaces 17, 20200656. doi:10.1098/rsif.2020.0656

CrossRef Full Text | Google Scholar

Cannon, W., Zucker, J., Baxter, D., Kumar, N., Baker, S., Hurley, J., et al. (2018). Prediction of metabolite concentrations, rate constants and post-translational regulation using maximum entropy-based simulations with application to central metabolism of neurospora crassa. Processes 6, 63. doi:10.3390/pr6060063

PubMed Abstract | CrossRef Full Text | Google Scholar

Covert, M. W., Schilling, C. H., and Palsson, B. (2001). Regulation of gene expression in flux balance models of metabolism. J. Theor. Biol. 213, 73–88. doi:10.1006/jtbi.2001.2405

PubMed Abstract | CrossRef Full Text | Google Scholar

Dewar, R. C. (2009). Maximum entropy production as an inference algorithm that translates physical assumptions into macroscopic predictions: Don’t shoot the messenger. Entropy 11, 931–944. doi:10.3390/e11040931

CrossRef Full Text | Google Scholar

Dixit, P. D., Wagoner, J., Weistuch, C., Presse, S., Ghosh, K., and Dill, K. A. (2018). Perspective: Maximum caliber is a general variational principle for dynamical systems. J. Chem. Phys. 148, 010901. doi:10.1063/1.5012990

PubMed Abstract | CrossRef Full Text | Google Scholar

Hackett, S. R., Zanotelli, V. R. T., Xu, W., Goya, J., Park, J. O., Perlman, D. H., et al. (2016). Systems-level analysis of mechanisms regulating yeast metabolic flux. Science 354, aaf2786. doi:10.1126/science.aaf2786

PubMed Abstract | CrossRef Full Text | Google Scholar

Hallows, W. C., Yu, W., and Denu, J. M. (2012). Regulation of glycolytic enzyme phosphoglycerate mutase-1 by sirt1 protein-mediated deacetylation. J. Biol. Chem. 287, 3850–3858. doi:10.1074/jbc.M111.317404

PubMed Abstract | CrossRef Full Text | Google Scholar

Hart, W. E., Laird, C. D., Watson, J.-P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., et al. (2017). Pyomo–optimization modeling in python second edn, 67. Springer Science & Business Media.

Google Scholar

Heimlicher, M. B., Bachler, M., Liu, M., Ibeneche-Nnewihe, C., Florin, E. L., Hoenger, A., et al. (2019). Reversible solidification of fission yeast cytoplasm after prolonged nutrient starvation. J. Cell Sci. 132, jcs231688. doi:10.1242/jcs.231688

PubMed Abstract | CrossRef Full Text | Google Scholar

Henry, C. S., Broadbelt, L. J., and Hatzimanikatis, V. (2007). Thermodynamics-based metabolic flux analysis. Biophysical J. 92, 1792–1805. doi:10.1529/biophysj.106.093138

CrossRef Full Text | Google Scholar

Hinshelwood, C. N. (1952). 136. On the chemical kinetics of autosynthetic systems. J. Chem. Soc. 745, 745. doi:10.1039/JR9520000745

CrossRef Full Text | Google Scholar

Holness, M. J., and Sugden, M. C. (2003). Regulation of pyruvate dehydrogenase complex activity by reversible phosphorylation. Biochem. Soc. Trans. 31, 1143–1151. doi:10.1042/bst0311143

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaynes, E. T. (1985). Macroscopic prediction. Berlin Heidelberg: Springer-Verlag, 254–269. doi:10.1007/978-3-642-70795-7

CrossRef Full Text | Google Scholar

Jitrapakdee, S., and Wallace, J. C. (1999). Structure, function and regulation of pyruvate carboxylase. Biochem. J. 340 (1), 1–16. doi:10.1042/bj3400001

PubMed Abstract | CrossRef Full Text | Google Scholar

King, E., and Cannon, W. (2023). Pathway controlled optimization. Available at: https://github.com/pnnl/pathway_controlled_optimization.

Google Scholar

Larsson-Raznikiewicz, M. (1967). Kinetic studies on the reaction catalyzed by phosphoglycerate kinase. ii. the kinetic relationships between 3-phosphoglycerate, mgatp2-and activating metal ion. Biochim. Biophys. Acta 132, 33–40. doi:10.1016/0005-2744(67)90189-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehninger, A. L., Nelson, D. L., and Cox, M. M. (2005). Lehninger principles of biochemistry 4th edn. New York: W. H. Freeman.

Google Scholar

Marcelin, R. (1910). The mechanics of irreversible phenomenon. Comptes Rendus Hebd. Des. Seances De. L Acad. Des. Sci. 151, 1052–1055.

Google Scholar

McCully, A. L., Onyeziri, M. C., LaSarre, B., Gliessman, J. R., and McKinlay, J. B. (2020). Reductive tricarboxylic acid cycle enzymes and reductive amino acid synthesis pathways contribute to electron balance in a rhodospirillum rubrum calvin-cycle mutant. Microbiology 166, 199–211. doi:10.1099/mic.0.000877

PubMed Abstract | CrossRef Full Text | Google Scholar

Newsholme, E. A., and Start, C. (1973). Regulation in metabolism. London: John Wiley & Sons.

Google Scholar

North, J. A., Narrowe, A. B., Xiong, W., Byerly, K. M., Zhao, G., Young, S. J., et al. (2020). A nitrogenase-like enzyme system catalyzes methionine, ethylene, and methane biogenesis. Science 369, 1094–1098. doi:10.1126/science.abb6310

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J. O., Rubin, S. A., Xu, Y. F., Amador-Noguez, D., Fan, J., Shlomi, T., et al. (2016). Metabolite concentrations, fluxes and free energies imply efficient enzyme usage. Nat. Chem. Biol. 12, 482–489. doi:10.1038/nchembio.2077

PubMed Abstract | CrossRef Full Text | Google Scholar

Parry, B. R., Surovtsev, I. V., Cabeen, M. T., O’Hern, C. S., Dufresne, E. R., and Jacobs-Wagner, C. (2014). The bacterial cytoplasm has glass-like properties and is fluidized by metabolic activity. Cell 156, 183–194. doi:10.1016/j.cell.2013.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Reznik, E., Christodoulou, D., Goldford, J. E., Briars, E., Sauer, U., Segrè, D., et al. (2017). Genome-scale architecture of small molecule regulatory networks and the fundamental trade-off between regulation and enzymatic activity. Cell Rep. 20, 2666–2677. doi:10.1016/j.celrep.2017.08.066

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, A., Connelly, S., Jiang, J., Zhuang, S., Amador, D. T., Phan, T., et al. (2014). Akt phosphorylation and regulation of transketolase is a nodal point for amino acid control of purine synthesis. Mol. Cell 55, 264–276. doi:10.1016/j.molcel.2014.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Sauro, H. M. (2018). Systems biology: An introduction to metabolic control Analysis. Ambrosius Publishing.

Google Scholar

Sivak, D. A., and Thomson, M. (2014). Environmental statistics and optimal regulation. PLoS Comput. Biol. 10, e1003826. doi:10.1371/journal.pcbi.1003826

PubMed Abstract | CrossRef Full Text | Google Scholar

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272. doi:10.1038/s41592-019-0686-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wächter, A., and Biegler, L. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program. 106, 25–57. doi:10.1007/s10107-004-0559-y

CrossRef Full Text | Google Scholar

Wang, Y. P., Zhou, L. S., Zhao, Y. Z., Wang, S. W., Chen, L. L., Liu, L. X., et al. (2014). Regulation of g6pd acetylation by sirt2 and kat9 modulates nadph homeostasis and cell survival during oxidative stress. EMBO J. 33, 1304–1320. doi:10.1002/embj.201387224

PubMed Abstract | CrossRef Full Text | Google Scholar

Waygood, E. B., Mort, J. S., and Sanwal, B. D. (1976). The control of pyruvate kinase of escherichia coli. binding of substrate and allosteric effectors to the enzyme activated by fructose 1,6-bisphosphate. Biochemistry 15, 277–282. doi:10.1021/bi00647a006

PubMed Abstract | CrossRef Full Text | Google Scholar

Waygood, E. B., Rayman, M. K., and Sanwal, B. D. (1975). The control of pyruvate kinases of escherichia coli. ii. effectors and regulatory properties of the enzyme activated by ribose 5-phosphate. Can. J. Biochem. 53, 444–454. doi:10.1139/o75-061

PubMed Abstract | CrossRef Full Text | Google Scholar

Waygood, E. B., and Sanwal, B. D. (1974). The control of pyruvate kinases of Escherichia coli. J. Biol. Chem. 249, 265–274. doi:10.1016/s0021-9258(19)43120-7

PubMed Abstract | CrossRef Full Text | Google Scholar