Research article
Open Access

A proof for loop-law constraints in stoichiometric metabolic networks

BMC Systems Biology20126:140

https://doi.org/10.1186/1752-0509-6-140

©  Noor et al.; licensee BioMed Central Ltd. 2012

Received: 5 July 2012

Accepted: 30 October 2012

Published: 12 November 2012

Abstract

Background

Constraint-based modeling is increasingly employed for metabolic network analysis. Its underlying assumption is that natural metabolic phenotypes can be predicted by adding physicochemical constraints to remove unrealistic metabolic flux solutions. The loopless-COBRA approach provides an additional constraint that eliminates thermodynamically infeasible internal cycles (or loops) from the space of solutions. This allows the prediction of flux solutions that are more consistent with experimental data. However, it is not clear if this approach over-constrains the models by removing non-loop solutions as well.

Results

Here we apply Gordan’s theorem from linear algebra to prove for the first time that the constraints added in loopless-COBRA do not over-constrain the problem beyond the elimination of the loops themselves.

Conclusions

The loopless-COBRA constraints can be reliably applied. Furthermore, this proof may be adapted to evaluate the theoretical soundness for other methods in constraint-based modeling.

Keywords

Thermodynamics Gordan’s theorem Metabolism

Background

Constraint-based modeling has become a successful framework for the analysis of large and complex stoichiometric biochemical networks[1]. The underlying concept of this framework is that one can use the stoichiometry of each reaction in a reconstructed metabolic network[2, 3] and known bounds on reaction fluxes to compute metabolic flux for each reaction. These predictions represent allowable steady-state metabolic flux distributions in a cell under a given growth condition. Some example constraints include mass balance and metabolite uptake rates. One set of constraints that has been more challenging to implement are those associated with thermodynamic limitations. Without thermodynamic constraints, non-physical fluxes can be computed for some metabolic reactions if they produce an internal cycle (Figure1a). Such cycles of reactions violate a “loop law” that is analogous to Kirchhoff’s second law for electrical circuits, as discussed previously by Beard et al.[4]. Many approaches have successfully constrained these loops using known flux directionality[5], energy-balance equations[4], and known[6–11] or predicted[12] thermodynamic parameters. Loops have also been indirectly removed by minimizing network flux[6, 13–15], or by coupling flux to enzyme synthesis costs[16].
Figure 1

Loop law constraints on metabolic networks. (a) Metabolic network reconstructions frequently have sets of reactions that cycle all metabolites internally. The fluxes of these reactions are therefore unconstrained. (b) Metabolic network solutions are found within a convex space, which is enclosed by known constraints on metabolite inputs, outputs, and known fluxes. Loops result in unconstrained dimensions in the solution space (blue). By implementing loopless-COBRA constraints, all loop-containing solutions are removed, leaving only solutions that do not contain loops (orange).

A new approach, called loopless-COBRA, was recently presented[17]. Unlike previous loop-removal algorithms, this method does not necessitate extra inputs or data, such as metabolite concentrations or thermodynamic parameters. Basically, this method imposes the second law of thermodynamics by using a mixed-integer linear programming (MILP) approach to constrain flux solutions so that they obey the loop law. Thus, flux solutions from this method are all within the portion of the flux space that is devoid of loops (Figure1b). Excitingly, the loop removal improved the consistency of the simulation results[17] with experimental data[18]. Specifically, it provides more realistic flux values for reactions that normally contribute to loops in a model. Otherwise, flux predictions for such reactions would usually have to be ignored in any subsequent analysis.

The loopless-COBRA method was shown to work in various scenarios, and the paper that presented the approach provides an explanation for why this method works. However, there is no mathematical proof for its formulation as an optimization problem. Specifically, it does not demonstrate that the additional MILP constraints do not over-constrain the problem and eliminate some non-loop containing solutions. Since constraint-based methods attempt to only eliminate impossible in silico phenotypes (i.e., steady-state flux distributions that the cell cannot maintain), it is important to verify that solutions representing real phenotypes are not removed by accidentally over-constraining the problem.

Here, we address this issue by presenting a mathematical proof for the completeness and soundness of the loopless-COBRA method, thereby adding fundamental support and rigorous proof for the constraints presented by Schellenberger et al.[17].

Results and discussion

Formal definition of loop-law constraints

In loopless-COBRA, the constraints added to the linear problem are:
G i < 0 for all v i > 0 G i > 0 for all v i < 0 G i ∈ R for all v i = 0 N int G = 0
where v i are the flux variables and N int is the null-space matrix of S int (the stoichiometric matrix of internal reactions). The third constraint ( G i ∈ R ) is not actually a constraint, but a way to say that G i can have any value if v i  = 0. We can rewrite all these constraints succinctly as:
∃ G ∈ R n s.t. null ( S int ) · G = 0 ∧ ∀ i sign ( G i ) = − sign ( v i ) ∨ v i = 0
(1)

Formal definition of loops

In order to prove that this constraint eliminates loops (and only loops), we must first find a mathematical formulation for a loop, using the same notation as above. We thus define a loop as a nonzero vector x ∈ R n which satisfies the mass-balance equation for the internal reactions, i.e. S int  · x = 0. This means that although there is a nonzero net flux in some of the reactions, overall, the internal network is at steady-state (an obvious violation of the second law of thermodynamics). It is important to point out, that this equation for defining loops must not be confused with the steady-state assumption commonly used in flux balance analysis models, namely S·v = 0, where the full stoichiometric matrix (S) is used.

According to this definition, a flux distribution ( v ∈ R n ) will contain a loop if and only if there exists a vector x ∈ R n ∖ { 0 } which is consistent with the flux directions in v (i.e. x i is either zero or has the same sign as v i ) and is itself a loop (i.e. S int  · x = 0). Formally, v has a loop if and only if:
∃ x ∈ R n ∖ { 0 } s.t. ∀ i sign ( x i ) ∈ { sign ( v i ) , 0 } ∧ S int · x = 0
(2)

We have now finished laying the groundwork for our mathematical proof that loopless-COBRA is sound and complete. In order to do that, we are left only to show that Equation 1 is satisfiable if and only if Equation 2 is unsatisfiable (in other words, there are no loops).

Gordan’s theorem

We start our proof by quoting Gordan’s theorem: For all A ∈ R m × n exactly one of the following two statements is true:
( a ) ∃ x ∈ R + n ∖ { 0 } s.t. Ax = 0 ( b ) ∃ y ∈ R m s.t. A ⊤ y > 0

We will show that statement (a) in Gordan’s theorem is equivalent to having a loop (Equation 2) and statement (b) is equivalent to the MILP constraints used by loopless-COBRA (Equation 1). After doing so, we would easily reach the conclusion of the proof.

As a guidance for the following sections, one can see that statement (a) already resembles Equation 2 if we define A = S int . The only difference is that x is constrained to have only non-negative values (note the ‘+’ in R + n ), instead of being consistent with the sign of v. Corollaries 1 and 2 will show how we can overcome this discrepancy by defining A in a slightly different way.

At first glance, statement (b) might look unrelated to Equation 1. However, in the last part of our proof, we show that choosing G ∈ R n , which satisfies null(S int )G = 0, is the same as choosing y ∈ R m and then taking G = S int ⊤ · y . Only for sake of understanding the algebra, one can think of y as the vector of formation Gibbs energies, and of G as the vector of reaction Gibbs energies. The rest of the proof, like for statement (a), deals with adjusting A to fit with the non-positive values in v. The toy example in Figure2 shows how Gordan’s theorem corresponds to having or not having a loop, for a network with 3 compounds and 3 internal reactions.
Figure 2

Illustrative example for Corollary 2. This example shows a small network with 3 internal reactions (x2−4). The flux directions were chosen according to the direction of the arrows. The matrix A is the internal stoichiometric matrix. (a) A flux distribution is shown where all 3 internal reactions are active and form a loop. Therefore there is a solution (x2 = x3 = x4 = 1) for the mass balance equation Ax = 0. In this case, no solution exists for A⊤y > 0. Therefore this flux distribution will be eliminated by loopless-COBRA. (b) A loopless flux distribution, in this case x4is not active. There is no solution for Ax = 0(except for the trivial solution x = 0). Gordan’s theorem claims that there must be a solution for A⊤y > 0, e.g. the one shown in the figure. Thus, loopless-COBRA will not eliminate any such flux distributions.

Corollary 1

For all A ∈ R m × n and d ∈ {−1,0,1} n exactly one of the following two statements is true:
( 1 a ) ∃ x ∈ R n ∖ { 0 } s.t. ∀ i sign ( x i ) ∈ { d i , 0 } ∧ Ax = 0 ( 1 b ) ∃ y ∈ R m s.t. ∀ i sign ( A ⊤ y ) i = d i ∨ d i = 0

Proof

First, define a new matrix  that is the same as A, without the columns corresponding to d i  = 0 and where columns corresponding to d i  = −1 are multiplied by −1. Statement (1a) is true for A if and only if (a) is true for  . The forward direction is easily shown by removing the zeros from x, where d i  = 0, and negating values corresponding to d i  = −1 (as previously done for A). Reversing this process (i.e. taking a positive solution for  x ̂ = 0 , adding back the zeros and negating the same values) shows the other direction is true as well.

Likewise, statement (1b) for A is true if and only if (b) is true for  , since columns with d i  = −1 are negated in  and thus sign (  ⊤ y ) i = 1 . Columns with d i  = 0 have no other constraints in (1b) and the same goes for (b) since they are removed from  .

Therefore, Corollary 1 is directly derived from Gordan’s theorem. □

Since constraint-based models usually use a vector of real values ( v ∈ R n ) to represent the flux distribution, we subsequently change the formulation of the Corollary slightly to match.

Corollary 2

For all A ∈ R m × n and v ∈ R n , exactly one of the following two statements is true:
( 2 a ) ∃ x ∈ R n ∖ { 0 } s.t. ∀ i sign ( x i ) ∈ { sign ( v i ) , 0 } ∧ Ax = 0 ( 2 b ) ∃ y ∈ R m s.t. ∀ i sign ( A ⊤ y ) i = − sign ( v i ) ∨ v i = 0

Proof

Defining d i  ≡ sign(v i ), we get this directly from Corollary 1. Note that −sign(v i ) can be used in (2b), since the existence of y is equivalent to the existence of −y. □

This adjustment now allows us to apply Corollary 2 to constraint-based problems and show that it eliminates loops (see example in Figure2). In order to avoid any confusion, we point out that a solution for Ax = 0 is considered a loop only if A is the stoichiometric matrix of internal reactions. This should not be mistaken as the steady-state mass-balance equation which looks exactly the same, except that A contains both internal and external reactions.

Corollary 3

Adding the following constraint:
∃ y ∈ R m s.t. ∀ i sign ( A ⊤ y ) i = − sign ( v i ) ∨ v i = 0

is equivalent to eliminating all loops in a flux distribution v.

Proof

Using Corollary 2, all that must be shown is that statement (2a) is equivalent to having a loop. This is apparent since x is a vector in the null-space of A (i.e., a loop) and is consistent with the flux direction of v in each of its nonzero reactions. □

Note that the trivial case v = 0 can still be a solution and it should be explicitly avoided if necessary.

Applying Corollary 3 in loopless-COBRA

The added constraints in loopless-COBRA[17] are slightly different than in Corollary 3, namely:
∃ G ∈ R n s.t. null ( A ) · G = 0 ∧ ∀ i sign ( G i ) = − sign ( v i ) ∨ v i = 0
We claim here that both formulations are equivalent. The fundamental theorem of linear algebra states that the nullspace, null(A), is the orthogonal complement of the row space, image(A⊤). Therefore, we can say that null(A)·G = 0 if and only if G ∈ image(A⊤), so we can rewrite the constraint above as:
∃ G ∈ image ( A ⊤ ) s.t. ∀ i sign ( G i ) = − sign ( v i ) ∨ v i = 0

which is obviously equivalent to the constraint in Corollary 3.

Conclusions

Our results prove that the constraints proposed by Schellenberger, et al.[17] eliminate all flux solutions with loops and nothing more. This alleviates the concern as to if the loopless-COBRA constraints might eliminate true flux states. Furthermore, values in G are analogous to the change in Gibbs energy (Δ r G) of the reactions[17], and y values are analogous to the chemical potentials (or formation energies, Δ f G) of the compounds themselves. Since in most cases, there are fewer compounds than reactions (m < n), we believe that it is convenient and intuitive to use the new formulation.

In conclusion, this proof provides theoretical credibility for the loopless-COBRA constraint. However, as with any algorithmic MILP implementation, care must still be taken with respect to numerical limitations and the convergence of the optimization algorithm.

Lastly, we believe this proof may be extended to similar methods addressing loop elimination. We also hope that similar proofs will appear for other methods, since more rigorous mathematical treatments are needed in many published algorithms in computational biology to prove or disprove their correctness.

Declarations

Acknowledgements

We thank Or Sheffet, Bernhard Ø. Palsson and Uri Barenholz for mathematical assistance and helpful discussions. EN is grateful to the Azrieli Foundation for the award of an Azrieli Fellowship.

Authors’ Affiliations

(1)
Department of Plant Sciences, Weizmann Institute of Science
(2)
Department of Bioengineering, University of California San Diego
(3)
Wyss Institute for Biologically Inspired Engineering and Department of Genetics, Harvard Medical School

References

  1. Lewis NE, Nagarajan H, Palsson BO: Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat Rev Microbiol. 2012, 10 (4): 291-305.Google Scholar
  2. Thorleifsson SG, Thiele I: rBioNet: A COBRA toolbox extension for reconstructing high-quality biochemical networks. Bioinf (Oxford, England). 2011, 27 (14): 2009-10. 10.1093/bioinformatics/btr308.View ArticleGoogle Scholar
  3. Thiele I, Palsson BO: A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat Protoc. 2010, 5: 93-121.View ArticleGoogle Scholar
  4. Beard Da, Liang Sd, Qian H: Energy balance for analysis of complex metabolic networks. Biophys j. 2002, 83: 79-86. 10.1016/S0006-3495(02)75150-3.View ArticleGoogle Scholar
  5. Varma A, Palsson BO: Metabolic capabilities of escherichia coli: I. synthesis of biosynthetic precursors and cofactors. J Theor Biol. 1993, 165 (4): 477-502. 10.1006/jtbi.1993.1202.View ArticleGoogle Scholar
  6. Holzhütter HG: The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks. Eur J Biochem / FEBS. 2004, 271 (14): 2905-22. 10.1111/j.1432-1033.2004.04213.x.View ArticleGoogle Scholar
  7. Henry CS, Broadbelt LJ, Hatzimanikatis V: Thermodynamics-based metabolic flux analysis. Biophys J. 2007, 92 (5): 1792-1805. 10.1529/biophysj.106.093138.View ArticleGoogle Scholar
  8. Fleming RMT, Thiele I, Provan G, Nasheuer HP: Integrated stoichiometric, thermodynamic and kinetic modelling of steady state metabolism. J Theor Biol. 2010, 264 (3): 683-692. 10.1016/j.jtbi.2010.02.044.View ArticleGoogle Scholar
  9. Fleming RMT, Thiele I: von Bertalanffy 1.0: a COBRA toolbox extension to thermodynamically constrain metabolic models. Bioinf (Oxford, England). 2011, 27 (1): 142-143. 10.1093/bioinformatics/btq607.View ArticleGoogle Scholar
  10. Cogne G, Rügen M, Bockmayr A, Titica M, Dussap CG, Cornet JF, Legrand J: A model-based method for investigating bioenergetic processes in autotrophically growing eukaryotic microalgae: application to the green algae Chlamydomonas reinhardtii. Biotechnol Prog. 2011, 27 (3): 631-40. 10.1002/btpr.596.View ArticleGoogle Scholar
  11. Noor E, Bar-Even A, Flamholz A, Lubling Y, Davidi D, Milo R: An integrated open framework for thermodynamics of reactions that combines accuracy and coverage. Bioinf (Oxford, England). 2012, 27 (15): 2037-2044.View ArticleGoogle Scholar
  12. De Martino D, Figliuzzi M, De Martino A, Marinari E: A scalable algorithm to explore the Gibbs energy landscape of genome-scale metabolic networks. PLoS Computat Biol. 2012, 8 (6): e1002562-10.1371/journal.pcbi.1002562.View ArticleGoogle Scholar
  13. Lewis NE, Hixson KK, Conrad TM, Lerman Ja, Charusanti P, Polpitiya AD, Adkins JN, Schramm G, Purvine SO, Lopez-Ferrer D, Weitz KK, Eils R, König R, Smith RD, Palsson BO: Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol. 2010, 6 (390): 1-13.Google Scholar
  14. Klitgord N, Segre D: The importance of compartmentalization in metabolic flux models: yeast as an ecosystem of organelles. Genome Inform Ser. 2010, 2009 (22): 41-55.Google Scholar
  15. Schuetz R, Zamboni N, Zampieri M, Heinemann M, Sauer U: Multidimensional optimality of microbial metabolism. Science. 2012, 336 (6081): 601-604. 10.1126/science.1216882.View ArticleGoogle Scholar
  16. Thiele I, Fleming RMT, Bordbar A, Schellenberger J, Palsson BO: Functional characterization of alternate optimal solutions of Escherichia coli’s transcriptional and translational machinery. Biophys J. 2010, 98 (10): 2072-81. 10.1016/j.bpj.2010.01.060.View ArticleGoogle Scholar
  17. Schellenberger J, Lewis NE, Palsson BO: Elimination of thermodynamically infeasible loops in steady-state metabolic models. Biophys J. 2011, 100 (3): 544-53. 10.1016/j.bpj.2010.12.3707.View ArticleGoogle Scholar
  18. Lewis NE, Cho BK, Knight EM, Palsson BO: Gene expression profiling and the use of genome-scale in silico models of Escherichia coli for analysis: providing context for content. J bacteriol. 2009, 191 (11): 3437-44. 10.1128/JB.00034-09.View ArticleGoogle Scholar

Copyright

© Noor et al.; licensee BioMed Central Ltd. 2012

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Advertisement