A variational principle for computing nonequilibrium fluxes and potentials in genome-scale biochemical networks
Highlights
► We derive a convex optimization problem that simultaneously enforces many constraints in genome-scale biochemical networks. ► Constraints enforced are steady state mass conservation, energy conservation and the second law of thermo-dynamics. ► We establish, in an exact manner, the duality relationship between reaction rates and chemical potentials. ► Efficient polynomial-time algorithms exist for solving such convex optimization problems based on interior point methods.
Introduction
The biochemical system of any organism can be represented mathematically by a network of chemicals (nodes) and reactions (edges). To analyze these networks at genome scale, systems biologists often use a linear optimization technique called flux balance analysis (FBA) (Savinell and Palsson, 1992). Flux balance requires that the sum of fluxes into and out of each node in the network be zero. This is equivalent to Kirchhoff׳s current law in an electrical network. Recent work has sought to augment flux balance analysis with Kirchhoff׳s loop law for energy conservation as well as the second law of thermodynamics (Beard et al., 2002, Price et al., 2006, Yang et al., 2005, Nagrath et al., 2007, Fleming et al., 2010, Schellenberger et al., 2011).
The incorporation of thermodynamic constraints into genome-scale models has produced models that are biologically more realistic and reveal greater insight into the control mechanisms operating in these complex biological systems (Beard et al., 2002, Kummel et al., 2006, Xu et al., 2008, Garg et al., 2010, Liebermeister et al., 2010). See Soh and Hatzimanikatis (2010) for a recent broad review of the application of thermodynamic constraints to biochemical networks. The second law of thermodynamics may be applied to each reversible elementary reaction by constraining net reaction flux to be in the direction of a negative change in chemical potential (Burton et al., 1957). Constraints on net reaction flux can be incorporated within flux balance analysis as additional linear inequality constraints (Fleming et al., 2009). Software packages to quantitatively assign reaction directionality for genome-scale models of metabolism are available (Fleming and Thiele, 2011, Schellenberger et al., 2011).
In contrast to the addition of constraints arising from the second law of thermodynamics, the addition of energy conservation constraints has been problematic because the resulting equations are nonlinear and/or nonconvex. Previous attempts required computing the global solution of a nonconvex continuous optimization problem (Beard et al., 2002, Nagrath et al., 2007, Fleming et al., 2010), solving an NP-hard problem (Yang et al., 2005), or solving a mixed integer linear program (Henry et al., 2007, Schellenberger et al., 2011). Mixed integer programs have unpredictable computational complexity.
The purpose of this work is to show that Kirchhoff׳s loop law and the second law of thermodynamics arise naturally from the optimality conditions of a convex optimization problem with flux balance constraints. Furthermore, every set of reaction fluxes that satisfies Kirchhoff׳s loop law and the second law of thermodynamics must be optimal for some instance of this problem. This suggests that there is an underlying variational principle operating in biochemical networks. Our convex optimization formulation leads to polynomial-time algorithms (Ye, 1997) for computing steady state fluxes that also satisfy energy conservation and the second law of thermodynamics.
Section snippets
Linear resistive networks
Consider a simple electrical circuit consisting of current sources, batteries, and resistors, as illustrated in Fig. 1. This is a linear resistive network with m nodes and n edges, where the node variables represent potentials and the edge variables represent flows (or currents) in the network. The circuit topology is defined by a node-edge incidence matrix , and properties of the network are encoded in a set of data vectors: is a vector of current sources, is a vector
Biochemical networks
The mathematical representation of a biochemical network is the stoichiometric matrix S. Like A above, is a sparse incidence matrix that encodes the network topology. However, biochemical networks, unlike linear resistive networks, are nonlinear networks or hypergraphs. That is, a single edge may link many nodes to many nodes, and the entries in S, which are integer stoichiometric coefficients, are not confined to the set . Each row of S corresponds to an individual chemical
Thermodynamic constraints
Flux balance analysis predicts fluxes that satisfy steady state mass conservation but not necessarily energy conservation or the second law of thermodynamics. Whilst the domain of steady state mass conserved fluxes includes those that are thermodynamically feasible, additional constraints are required in order to guarantee a flux that additionally satisfies energy conservation and the second law of thermodynamics. Recent work has tried to add constraints based on Kirchhoff׳s loop law and the
A variational principle
We now present the main theorem of this paper. The theorem introduces a new convex optimization problem with the same flux balance constraints as problem (FBA) but with a negative entropy objective function. It states that the thermodynamic constraints (4), (5) hold at its unique solution. We use e to denote a vector of ones. Theorem 1 Let be any set of optimal exchange fluxes from problem (FBA). Define , and let c be any vector in . The convex equality-constrained problem
Dual variational principle
Given the primal optimization problem (EP) it is instructive to consider the equivalent dual optimization problem, as its objective function lends insight into the properties of the optimal dual variables for the primal problem.
By Lemma 1, the constraints of Problem (EP) have feasible solutions that are strictly positive. Similarly, there must be feasible solutions that are finite, and these have finite objective values. Hence the optimum of the strictly convex objective function of problem
Thermodynamic feasibility and mass action kinetics
From the optimality conditions for problem (EP), observe that forward flux is a function of substrate and product chemical potentials. For an elementary reaction, one would expect that forward flux should depend only on substrate chemical potential(s) and reverse flux should depend only on product chemical potential(s) (Ederer and Gilles, 2007). According to mass action kinetics, the rate of an elementary forward reaction only depends on a forward kinetic parameter and substrate
Discussion
With problem (EP) and Theorem 1 we provide a new variational principle that enforces steady state mass conservation, energy conservation and the second law of thermodynamics for genome-scale biochemical networks. Moreover, we prove that all thermodynamically feasible steady state fluxes are instances of problem (EP) for a different free parameter vector (Theorem 2). The values of the free parameters influence the optimal forward and reverse fluxes through the relation .
Acknowledgments
We would like to acknowledge invaluable assistance from Ines Thiele and some helpful thermodynamic discussions with Hong Qian. We would also like to thank two anonymous reviewers for their constructive comments and for directing us to the 1951 paper by W. Millar (Millar, 1951). This project was supported by the U.S. Department of Energy (Office of Advanced Scientific Computing Research and Office of Biological and Environmental Research) as part of the Scientific Discovery Through Advanced
References (52)
- et al.
Energy balance for analysis of complex metabolic networks
Biophys. J.
(2002) - et al.
Thermodynamically feasible kinetic models of reaction networks
Biophys. J.
(2007) - et al.
Quantitative assignment of reaction directionality in constraint-based models of metabolism: application to Escherichia coli
Biophys. Chem.
(2009) - et al.
Integrated stoichiometric, thermodynamic and kinetic modeling of steady state metabolism
J. Theor. Biol.
(2010) - et al.
Thermodynamics-based metabolic flux analysis
Biophys. J.
(2007) - et al.
Genome-scale thermodynamic analysis of Escherichia coli metabolism
Biophys. J.
(2006) - et al.
Group contribution method for thermodynamic analysis of complex metabolic networks
Biophys. J.
(2008) - et al.
Tellegen׳s theorem and thermodynamic inequalities
J. Theor. Biol.
(1971) - et al.
Candidate states of Helicobacter pylori׳s genome-scale metabolic network upon application of “loop law” thermodynamic constraints
Biophys. J.
(2006) - et al.
Thermodynamics of stoichiometric biochemical networks in living systems far from equilibrium
Biophys. Chem.
(2005)
Network analysis of intermediary metabolism using linear optimization. I. Development of mathematical formalism
J. Theor. Biol.
Elimination of thermodynamically infeasible loops in steady-state metabolic models
Biophys. J.
Network thermodynamics in the post-genomic era
Curr. Opin. Microbiol.
Ab initio prediction of thermodynamically feasible reaction directions from biochemical network stoichiometry
Metab. Eng.
Thermodynamics of Biochemical Reactions
Biochemical Thermodynamics: Applications of Mathematica
A computational study of the homogeneous algorithm for large-scale convex optimization
Comput. Optim. Appl.
Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli
Nat. Chem. Biol.
Nonlinear Programming
Convex Optimization
Energy Transformations in Living Matter
Fundamentals of Enzyme Kinetics
Enzymes in context: kinetic characterization of enzymes for systems biology
Biochemist
Primal-dual solution perturbations in convex optimization
Set-Valued Anal.
Nonlinear networks IIa
Bull. Am. Math. Soc.
Cited by (29)
Structural conserved moiety splitting of a stoichiometric matrix
2020, Journal of Theoretical BiologyCitation Excerpt :To this end, the specific structure of stoichiometric matrices should be taken into account to design novel context-specific algorithmic methodologies for such big data problems. Several classes of optimisation models can be developed for responding to different problems arising in system biology such as linear (Ma et al., 2017), quadratic (Segré et al., 2002), convex (Fleming et al., 2012), duplomonotone (Artacho and Fleming, 2014) and other nonlinear problems (Ahookhosh et al., 2019, Ahookhosh et al. 2020.). For each class of optimisation models, the specific structure of the stoichiometric matrices might be used to decompose the original problem to several computationally tractable subproblems using the moiety matrix splitting described in Section 7.
Quantitative systems pharmacology and the personalized drug-microbiota-diet axis
2017, Current Opinion in Systems BiologyCitation Excerpt :To reliably and efficiently integrate PBPK and COBRA at large scale, using the available resources (Table1, Figure 2), the available hybrid modeling approaches needs to be refined. For instance, to ensure that the changes in calculated flux vectors between two consecutive time steps respond as well-behaved functions (single-valued and Lipschitz continuous) of perturbations in the input data, the flux balance analysis problem can be regularized to ensure that it is a strictly convex optimization problem [68,69] by minimizing the Euclidean distance between the current and next optimal flux vector. In addition, one could mathematically reformulate the problem by adapting techniques from discrete mechanics to pose the dynamic trajectory of an integrated pharmacokinetic and metabolic system as the optimal solution to a variational integration problem [70].
Estimation of flux distribution in metabolic networks accounting for thermodynamic constraints: The effect of equilibrium vs. blocked reactions
2016, Biochemical Engineering JournalThe Influence of Crowding Conditions on the Thermodynamic Feasibility of Metabolic Pathways
2015, Biophysical JournalGeneric flux coupling analysis
2015, Mathematical BiosciencesCitation Excerpt :Loop-law thermodynamic constraints have the advantage that they do not require data on equilibrium constants and metabolite concentrations. Therefore, they have also been widely studied in the literature [3,13,16,20,28,30,32,34,35,38,47,49] The key property of lattices is that metabolic pathways can be combined by taking the union of the supports.