Elsevier

Journal of Theoretical Biology

Volume 292, 7 January 2012, Pages 71-77
Journal of Theoretical Biology

A variational principle for computing nonequilibrium fluxes and potentials in genome-scale biochemical networks

https://doi.org/10.1016/j.jtbi.2011.09.029Get rights and content

Abstract

We derive a convex optimization problem on a steady-state nonequilibrium network of biochemical reactions, with the property that energy conservation and the second law of thermodynamics both hold at the problem solution. This suggests a new variational principle for biochemical networks that can be implemented in a computationally tractable manner. We derive the Lagrange dual of the optimization problem and use strong duality to demonstrate that a biochemical analogue of Tellegen׳s theorem holds at optimality. Each optimal flux is dependent on a free parameter that we relate to an elementary kinetic parameter when mass action kinetics is assumed.

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 yRm represent potentials and the edge variables xRn represent flows (or currents) in the network. The circuit topology is defined by a node-edge incidence matrix ARm×n, and properties of the network are encoded in a set of data vectors: fRm is a vector of current sources, bRn is a vector

Biochemical networks

The mathematical representation of a biochemical network is the stoichiometric matrix S. Like A above, SRm×n 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 {1,0,1}. 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 ve be any set of optimal exchange fluxes from problem (FBA). Define b=Seve, and let c be any vector in Rn. The convex equality-constrained problemminimizevf,vr>0ϕvf

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 cRn (Theorem 2). The values of the free parameters influence the optimal forward and reverse fluxes through the relation vf.vr=exp(2c).

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)

  • J.M. Savinell et al.

    Network analysis of intermediary metabolism using linear optimization. I. Development of mathematical formalism

    J. Theor. Biol.

    (1992)
  • J. Schellenberger et al.

    Elimination of thermodynamically infeasible loops in steady-state metabolic models

    Biophys. J.

    (2011)
  • K. Soh et al.

    Network thermodynamics in the post-genomic era

    Curr. Opin. Microbiol.

    (2010)
  • F. Yang et al.

    Ab initio prediction of thermodynamically feasible reaction directions from biochemical network stoichiometry

    Metab. Eng.

    (2005)
  • Akle, S., Dalal, O., Fleming, R.M.T., Saunders, M.A., Taheri, N.A., Ye, Y., 2011. Existence of Positive Equilibria for...
  • R.A. Alberty

    Thermodynamics of Biochemical Reactions

    (2003)
  • R.A. Alberty

    Biochemical Thermodynamics: Applications of Mathematica

    (2006)
  • E.D. Andersen et al.

    A computational study of the homogeneous algorithm for large-scale convex optimization

    Comput. Optim. Appl.

    (1998)
  • B.D. Bennett et al.

    Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli

    Nat. Chem. Biol.

    (2009)
  • D.P. Bertsekas

    Nonlinear Programming

    (1999)
  • S. Boyd et al.

    Convex Optimization

    (2004)
  • K. Burton et al.

    Energy Transformations in Living Matter

    (1957)
  • A. Cornish-Bowden

    Fundamentals of Enzyme Kinetics

    (2004)
  • A. Cornish-Bowden et al.

    Enzymes in context: kinetic characterization of enzymes for systems biology

    Biochemist

    (2005)
  • A.L. Dontchev et al.

    Primal-dual solution perturbations in convex optimization

    Set-Valued Anal.

    (2001)
  • R.J. Duffin

    Nonlinear networks IIa

    Bull. Am. Math. Soc.

    (1947)
  • Cited by (29)

    • Structural conserved moiety splitting of a stoichiometric matrix

      2020, Journal of Theoretical Biology
      Citation 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 Biology
      Citation 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].

    • Generic flux coupling analysis

      2015, Mathematical Biosciences
      Citation 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.

    View all citing articles on Scopus
    View full text