The Essential Role of Thermodynamics in metabolic network modeling: physical insights and computational challenges

TL;DR

This paper emphasizes the importance of thermodynamic constraints in metabolic network modeling, discussing computational challenges and physical insights to improve the accuracy of steady-state analyses.

Contribution

It reviews the integration of thermodynamics into metabolic models, highlighting computational challenges and proposing variational principles for better constraint implementation.

Findings

Thermodynamic constraints prevent unfeasible loops in models.

Gibbs inequalities are key to ensuring physical feasibility.

Exploration of non-convex flux spaces requires advanced algorithms.

Abstract

Quantitative studies of cell metabolism are often based on large chemical reaction network models. A steady state approach is suited to analyze phenomena on the timescale of cell growth and circumvents the problem of incomplete experimental knowledge on kinetic laws and parameters, but it shall be supported by a correct implementation of thermodynamic constraints. In this article we review the latter aspect highlighting its computational challenges and physical insights. The simple introduction of Gibbs inequalities avoids the presence of unfeasible loops allowing for correct timescale analysis but leads to possibly non-convex feasible flux spaces, whose exploration needs efficient algorithms. We shorty review on the implementation of thermodynamics through variational principles in constraints based models of metabolic networks.

Tables

Table 1. Table 1: Turnover times of selected compounds in E Coli metabolism simulated for a genome scale model in aerobic glucose limited minimal environment, from net fluxes and corrected for fluctuations.

Metabolic	Turnover time	Turnover time
compound	estimated from net flux (s)	corrected for fluctuations (s)
ATP	$2.0\pm 0.1$	$0.4\pm 0.2$
ADP	$0.120\pm 0.005$	$0.02\pm 0.01$
AMP	$0.5\pm 0.1$	$0.11\pm 0.06$
NAD	$1.1\pm 0.1$	$0.3\pm 0.1$
NADH	$3.5\pm 0.2\cdot {10}^{-2}$	$1.0\pm 0.4\cdot {10}^{-2}$
NADP	$1.6\pm 0.2\cdot {10}^{-3}$	$2\pm 1\cdot {10}^{-4}$
NADPH	$9\pm 1\cdot {10}^{-2}$	$2\pm 1\cdot {10}^{-2}$
Glutammate	$90\pm 20$	$16\pm 8$
3-Phosphoglycerate	$2.0\pm 0.2$	$0.12\pm 0.6$

Equations

$$\dot{\mathbf{c}}=\mathbf{S}\cdot \mathbf{v}$$

$$\mathbf{S}\cdot \mathbf{v}=0$$

$$\mathbf{S}\cdot \mathbf{v}=0,$$

$$v_r\in [v_r^{\min },v_r^{\max }]$$

$$\sum _{\mu }{\xi }_{i\mu }{g}_{\mu }>0\forall i$$

$$\sum _i{S}_{i\mu }{k}_i=0k_i\ge 0\forall \mu$$

$$G\equiv E-PV-TS,$$

$${\Delta }_i\equiv -{\delta }_i\sum _{\alpha }{S}_{j,\mu }g\mu \ge 0\forall i.$$

$$P_0\left(\mathbf{g}\right)=\prod _{\alpha =1}^M{P}_0^{\mu }\left(\mathbf{g}\right),$$

$$i_0=\text{mbox}arg\underset{i}{\min }{\Delta }_i,$$

$$\mathbf{g}⟶\mathbf{g}-\lambda {\delta }_{i_0}{\mathbf{S}}_{i_0};$$

$$\sum _{i\in \mathcal{I}}{k}_i{\delta }_i{S}_{\mu ,i}=0\forall \mu .$$

$$\dot{c}=P-D=0,\tau =\frac{c}{P}$$

$$P_{\mu }=\sum _i\theta (S_{i\mu }{v}_i)S_{i\mu }{v}_i$$

$$\frac{v^{+}}{v^{-}}=e^{-\Delta G/RT}$$

$$\dot{c}=(P^{+}+D^{-})-(D^{+}+P^{-})=0$$

$$\tau =\frac{c}{P^{+}+D^{-}}$$

$$\frac{{\tau }_0-\tau }{\tau }\simeq \frac{1}{e^{-\frac{\Delta G_2}{RT}}-1}\simeq 45\%$$

$$Q(\mathbf{v})=\frac{1}{N}\sum _{i=1}^N{v}_i^2,$$

$$\sum _i{\Omega }_i^{\mu }{k}_i=0,$$

$$w_i=v_i+\alpha k_i{v}_i,$$

$$Q(\mathbf{v})=Q(\mathbf{w})+\frac{{\left({\sum }_i{k}_i{v}_i^2\right)}^2}{{\sum }_i{k}_i^2{v}_i^2}.$$

$$\mathcal{F}({\mathbf{v}}_F,{\mathbf{v}}_R)={\mathbf{v}}_F\cdot \left[\log ({\mathbf{v}}_F+\mathbf{c}-1)\right]+{\mathbf{v}}_R\cdot \left[\log ({\mathbf{v}}_R+\mathbf{c}-1)\right]$$

$${\mathbf{S}}_I({\mathbf{v}}_F-{\mathbf{v}}_R)+{\mathbf{S}}_E{\mathbf{v}}_E^{\star }=0,$$

$$\mathbf{g}=-2RT\mathbf{\lambda }.$$

$$H=\sum _{\mu \in \mathrm{Ext}}\frac{(u^{\mu }{)}^2}{c_{\mathrm{ext}}^{\mu }}$$

$$\mathbf{Sv}=\mathbf{u}$$

Full text

Chapter 0 The Essential Role of Thermodynamics in metabolic network modeling: physical insights and computational challenges

\body

Introduction

Because of its uniquely universal nature, thermodynamics has been linked to physiology from its very inception, both to rationalize observations and to elucidate fundamental limits to physiological functions. With the advent of genetics and molecular biology, the discovery of the molecular mechanisms underlying physiology became the primary challenge. However, as the molecular actors and their interactions were mapped out at increasingly fine resolution, the focus gradually shifted on understanding their system-level organization [1]. And, perhaps unsurprisingly, it has become more and more clear that thermodynamic aspects are crucial for the emergent large-scale behaviour of these systems. Currently, renewed interest has flourished around the thermodynamics of cellular processes, only this time with the possibility of relying on a host of data at various scales for quantitative analyses [2].

In no area of physiology is thermodynamic analysis more central than in metabolic network modeling [3]. In brief, metabolic networks encode for the set of chemical reactions that, in any cell, break down nutrients and harvest free energy to synthetize the macro-molecular building blocks essential to life (amino acids, nucleotides, fatty acids, etc.) and, ultimately, biomass. Their structures can be inferred by combining gene-enzyme-reactions associations with regulatory information and transcriptional data. The availability of detailed metabolic network reconstructions for a large number of organisms and cell types is, in our view, among the most significant successes that computational methods have reaped in biology to date [4].

Building reliable and predictive dynamical models of metabolism based on this information is however challenging, mainly due to our vastly incomplete knowledge about intracellular enzyme kinetics, transport mechanisms and rate constants. On the other hand, non-equilibrium steady state approaches appear to be more feasible. Such methods are perhaps best represented by the broad class of computational schemes known as ‘constraint-based models’ [5, 6]. From a physical viewpoint, such models should essentially rely on two “Kirchhoff-type” assumptions regarding (a) mass balance for chemical species (i.e. metabolic homeostasis) and (b) energy balance for reactions (i.e. thermodynamic feasibility of material fluxes) [7]. In the most basic setup, energy balance simply requires reaction fluxes at steady state to proceed downhill in free energy, in accordance with the second law of thermodynamics. Unfortunately, implementing this constraint in genome-scale models is drastically harder than enforcing the stationarity of metabolite concentrations. Inclusion of thermodynamic constraints is however essential not only to obtain physically viable flux patterns, but also to highlight timescales and turnover rates and to allow for the estimation of metabolite concentrations. Needless to say, the range of applications of such results, from biotechnology to pharmacology, would be enormous.

Our main goal here is to present this problem and its multiple ramifications, which span from basic biochemistry to some fundamental algorithmic challenges, under a statistical physics lens. We shall discuss what, in practice, makes it so hard to solve efficiently, and review some of the alternative approaches that have been attempted. Finally, we will point to some recent developments that may hold some promising keys to finally unlock the puzzle of metabolic network thermodynamics at genome resolution.

Background

In the most simple setting metabolism can be modeled in terms of the dynamics of the chemical compounds concentration levels [8]. Upon assuming well-mixing and neglecting noise, we still have a large possibly non-linear dynamical system whose parameters could be not known in their entirety. For a chemical reaction network in which $M$ metabolites participate in $N$ reactions with the stoichiometry encoded in a matrix $\mathbf{S}=\{S_{\mu i}\}$, the concentrations $c_{\mu}$ change in time according to mass-balance equations

[TABLE]

where a component of the vector ${\bf v}$, $v_{i}$, is the flux of the reaction $i$ that is in turn a (possibly unknown) function of the concentration levels $v_{i}(\mathbf{c})$ (and several other parameters, like enzyme copy number, etc). On the other hand, in order to analyze phenomena with timescales longer than diffusion and typical turnover times (like cell growth) it is possible to assume a steady state, i.e. a flux configuration satisfying

[TABLE]

In so called constraints-based modeling, apart from mass balance, fluxes are bounded in certain ranges $v_{r}\in[v_{r}^{{\rm min}},v_{r}^{{\rm max}}]$ that take into account thermodynamic irreversibility, kinetic limits and physiological constraints. The set of constraints

[TABLE]

defines a convex closed set in the space of reaction fluxes: the polytope of feasible steady states. The productive capabilities of the network can be investigated computationally by maximizing suited linear objective functions in the aforementioned space [9], in particular the biomass growth itself (flux balance analysis [10], based on linear programming). On the other hand more generic inference problems can be afforded quite efficiently with Monte Carlo methods given the convexity of the space [11]. It shall be noted that thermodynamics is implemented in a very simple way, i.e. by setting reaction reversibility, i.e. $f_{i}\geq 0$ for some reactions. When flux bounds are not provided, it is customary to set them to an arbitrary large number, $f_{i}\in\left[-C,C\right]$ (typically $C=10^{3},10^{4}$), that, for a meaningful model, shall not influence the results. On the other hand a more rigorous yet simple approach consists in postulating that fluxes shall follow a free energy gradient. If $f_{i}\neq 0$, then $f_{i}\Delta G_{i}<0$ (Gibbs inequality) where $\Delta G_{i}$ is the free energy change of reaction $i$. The $\Delta G$’s can be written as the difference between the chemical potentials $g_{\mu}$ of products and substrates through the stoichiometric matrix $\Delta G_{i}=\sum_{\mu}S_{i\mu}g_{\mu}$. In terms of chemical potentials we have thus a system of linear inequalities ($\xi_{i\mu}=-sign(f_{i})S_{i\mu}$)

[TABLE]

whose feasibility (existence of a solution $g_{\mu}$) is necessary for the thermodynamical feasibility of the flux configuration. Even in this basic approach, the addition of free energy variables ($g_{\mu}$ and $\Delta G_{i}$) makes the problem non-linear, in particular quadratic and possibly non-convex. On the other hand, upon conditioning on flux variables we can get useful hints from duality theorems of the alternative that characterize thermodynamic feasibility in terms of the unfeasibility of particular flux configurations [12]. Specifically, according to the Gordan theorem we have that system (4) has a feasible solution if and only if

[TABLE]

has no non-trivial solutions (unfeasible loops). Such a duality can be exploited in order to define efficient algorithms as we discuss in the following section.

Relaxational Algorithms

Reconstruction of complete metabolic networks is becoming, thanks to new, modern and accurate experiments, a possible option for many simple organisms. What is important in the present context is that this reconstruction needs to be compatible with thermodynamic principles, and that implementing thermodynamic requirements can help in an accurate reconstruction of the network. This is sometimes a difficult computational problem, and we describe here an useful algorithm to implement thermodynamic consistency and use this to help in a (correct) network reconstruction [13].

We will in this way gather information about Gibbs free energy and about their landscape, that will have to be compatible with the selected vector of reaction directions. We use stoichiometric information via a constructive algorithm inspired by perceptron [14] approach learning. In the method we use the a preliminary reconstruction of the network structure to iteratively build up correlations between the chemical potentials of the chemical species, until we reach a thermodynamically consistent profile. The algorithm is nicely scalable, and it can allow the crucial result of guaranteeing the feasibility of flux configurations, or of identifying and removing unfeasible cycles. The algorithm can also be useful to get an estimate of reaction affinities, and it can be used to derive bounds for concentrations.

We consider the Gibbs energy at temperature $T$ and volume $V$

[TABLE]

where $E$ is the internal energy of the systems, $P$ its pressure and $S$ the entropy. Let us call $\delta_{i}$ the direction of chemical reaction $i$: $\delta_{i}=\pm 1$, i.e. the reaction can proceed in the ”forward” direction or in the reverse direction. If $\Delta G_{i}$ is the Gibbs energy difference induced by reaction $i$ one needs that $\delta_{i}\,\Delta G_{i}\leq 0\;\forall i$ (the Gibbs energy cannot increase in the direction where the chemical reaction operates). Let us consider now the stoichiometric coefficients, that are $S_{\mu,j}<0$ for substrates and $S_{\mu,j}>0$ for products. Let us also define the vector $\boldsymbol{g}\equiv\{g_{\mu}\}$ of the Gibbs energies per mole of species $\mu$. In terms of $\boldsymbol{g}$ we have that $\Delta G_{i}=S^{T}_{j,\mu}\;g_{\mu}$. Now we can discuss the following problem. Given a set of reaction directions $\{\delta_{i}\}$ (that have been inferred in a first step of the procedure) determine, if it exists, $\boldsymbol{g}$ such that

[TABLE]

For fixed $\delta_{i}$ the solution space is convex. Relaxation methods are a typical and potentially effective choice for solving a problem of this kind. If a solution does not exist the reconstructed network is not consistent and it has, at best, to be cured.

To solve this problem we have introduced [13] an algorithm based on the so called MinOver approach [14] which was originally developed for neural network learning. The algorithm starts from a configuration of the $g_{\mu}$, that is extracted under the probability distribution $P_{0}(\boldsymbol{g})$. $P_{0}$ is selected a priori after phenomenological considerations. All the experimental input to the algorithm is indeed in the choice of $P_{0}$. We assume the simple ansatz

[TABLE]

where $P_{0}^{\mu}$ is uniform around the values estimated from experiments with a range also suggested from experimental data. We now generate a random vector under $P_{0}\left(\boldsymbol{g}\right)$ and we compute the vector $\boldsymbol{\Delta}=\{Delta_{i}\}$ from 7. Let us call $i_{0}$ the index of the most broken constraint, i.e. let us set

[TABLE]

Now if $\Delta_{i_{0}}\geq 0$ $\boldsymbol{g}$ is a thermodynamically consistent chemical potential, and solves our problem. We can accept it and exit, or look for more solutions (including the one already found) by restarting the algorithm from a different seed. If instead $\Delta_{i_{0}}<0$ $\boldsymbol{g}$ we do not have a solution and we update $\boldsymbol{g}$ by setting

[TABLE]

with $\lambda$ an appropriate constant. One iterates till convergence, that is guaranteed (maybe after a very long, unpractical time) if a solution exists.

If the problem has no solutions the assignment of the directions is not consistent. This happens if and only if there is at least one unfeasible loop. The main problem is at this point that is not easy to find a loop that can be, in principle, also very long. We can phrase better the problem by saying that there is an unfeasible loop if there is a set $\mathcal{I}$ of reactions such that a set of positive constants $k_{i}>0$ exists such that

[TABLE]

If there is a loop the algorithm does not converge, since the least satisfied constraint rotates on the loop and does not get fixed. If this is happening we can start from our algorithm to localize and kill not too long loops in the following way. After discarding the first part of the iteration steps (where we typically are in a transient region) we start storing the values of $i_{0}(t)$ (where $t$ labels the iterations of the procedure), i.e. the value of the most broken constraint at iteration $t$ of the MinOver procedure. Now we look among the reactions appearing in the set of the $\{i_{0}(t)\}$ for $t>\tilde{T}$, where we can vary the minimum iteration $\tilde{T}$, and we search for loops of length $\mathcal{L}$, starting from $\mathcal{L}=3$, and increasing it if needed.When we find a loop we change one of the directions and try the MinOver procedure again. The loops in the dual space can be searched exhaustively or with aid of Monte Carlo methods [15]

The possibility of improving the reconstruction of a network, making it compatible with thermodynamical basic principles, is important, and our algorithm helps in this direction.

Flux scales

It is interesting to notice that scales analysis in metabolic networks need implementation of thermodynamics constraints beyond reversibility assignment. This leads to a geometrically more complex picture of the flux space, but its lack leads possibly to wrong conclusions. Previous work on sampling the flux space seemed to show scale free distributions [16], in contradiction with the existence of physical limiting factors, e.g. resources availability [17] or maximum ribosome elongation rate [18]. As we stressed in the background, models of metabolic networks come with arbitrary bounds on the fluxes, upon which the solution space could depend, in turn hampering scaling analysis, in particular in presence of the aforementioned unfeasible loops. Suppose in fact that $f_{i,0}$ is a feasible state and $k_{i}$ is a solution of (5), the line $f_{i}(L)=f_{i,0}+k_{i}s_{i}L$ verifies the steady state mass balance constraints by construction and it will be inside the polytope till $L$ possibly reaches the arbitrary bounds.

A simple illustration is depicted in Fig. 1: $A$ is injected with rate $u$, it can be transformed either in $B$ and subsequently in $C$ with rate $x$ or directly in $C$ with rate $u-x$, that is consumed with rate $u$. (if $u<0$, $C$ is injected and $A$ consumed). The variable $x$ is unbounded unless closed loops are forbidden, leading to the non linear constraints $ux\geq 0$, $|x|\leq|u|$. In general thermodynamics forbids the orthants in the flux space that include closed loops: the remaining feasible space is not convex anymore (see Fig. 1), but its scales now reflect true physical constraints. These issues have been studied in a genome scale metabolic network model, specifically the typical steady states of the E Coli metabolic network iJR904 [19] in a glucose limited minimal medium in aerobic conditions [20]. Flux configurations have been uniformly sampled and corrected from unfeasible loops with the methods described in the previous section. Results for the distribution of flux intensities $|f_{i}|$ are shown in Fig. 2, before and after correcting for unfeasible loops.

The long tail corresponding to the uniform sampling depends on the arbitrary constant fixed for flux bounds. After thermodynamic correction the flux intensity distribution has a cut-off that scales simply with the glucose input [20].

Enhanced turnover

What is the overall turnover time of metabolism? This question receives useful hints from thermodynamic information even if the knowledge of reactions kinetic is lacking. In particular fluctuation analysis upon application of the fluctuation theorem returns a faster picture of metabolism with respect to standard turnover estimates [21]. In general and formally relaxation times are calculated from a linear stability analysis of the steady states, but this requires detailed knowledge of reactions kinetics. On the other hand if at least fluxes and concentrations are experimentally known, it is possible to calculate the metabolites turnover times $\tau$, i.e. the ratio between the concentration $c$ and the net flux of production $P$ (or equivalently consumption $D$, given the steady state), schematically

[TABLE]

where e.g. the flux $P$ can be calculated from the network and the flux configuration for a given metabolite $\mu$ ($\theta$ is the Heaviside step function):

[TABLE]

Such turnover time is the typical time it takes to fully replenish a given metabolic pool. On the other hand net fluxes result from the difference between forward and backward contributions $\nu=\nu^{+}-\nu^{-}$, and the latters can be estimated by the fluctuation theorem if the free energy $\Delta G$ change is known (cit):

[TABLE]

Upon taking into account the backward contribution the turnover time can be shorter, i.e. schematically

[TABLE]

For instance consider Glucose-6-phosphate in the human red blood cell. This is produced the Hexokinase enzyme ($\Delta G_{1}\simeq-29$ KJ/mol), and consumed by the phosphoglucoisomerase enzyme ($\Delta G_{2}\simeq-2.9$ KJ/mol). At $RT=2.5$KJ/mol the turnover time $\tau_{0}$ estimated only from net fluxes overestimates the one $\tau$ that takes into account backward contribution by a factor

[TABLE]

Such analysis has been performed for the genome scale E.Coli metabolic network iJR904 [19] in a glucose limited minimal medium in aerobic conditions [21], returning a faster picture of intermediate metabolism, that we summarize in table 1 reporting the turnover times estimate from net and total fluxes of the metabolites ruling the energetics of the network.

Variational principles

An alternative route to implementing thermodynamic feasibility at network scale consists in devising global variational principles ensuring that optimal mass-balanced flux patterns are void of cycles . The simplest such principle is perhaps given by the minimization of the total flux

[TABLE]

where it is understood that flux vectors $\mathbf{v}=\{v_{i}\}$ satisfy the mass balance conditions $\mathbf{Sv=0}$ with pre-defined ranges of variability for fluxes (as well as the additional constraints that may be required on a case by case basis). An argument proving that $Q$ is minimum for flux configurations that are thermodynamically viable (assuming such configurations exist for the network under study) is as follows. Consider a mass-balanced flux configuration $\mathbf{v}$ and assume it contains an infeasible cycle. Such a cycle must be described by a non-zero solution $\mathbf{k}$ to the system

[TABLE]

where $\Omega_{i}^{\mu}=-v_{i}S_{i}^{\mu}$. Now consider the flux configuration $\mathbf{w}$ defined by

[TABLE]

with $\alpha$ a constant. Clearly, if $\mathbf{v}$ is mass-balanced, so is $\mathbf{w}$. However $Q(\mathbf{v})>Q(\mathbf{w})$ provided $\alpha$ is such that $\frac{\partial Q(\mathbf{w})}{\partial\alpha}=0$. In particular, one finds

[TABLE]

In other terms, given a thermodynamically infeasible mass-balanced flux configuration it is always possible to construct another mass-balanced flux configuration whose total flux is lower. In turn, the resulting flux pattern has to be thermodynamically feasible when $Q$ is minimized.

This idea, originally put forward in [22], has been applied in various computational schemes for genome-scale metabolism, such as pFBA (parsimonious Flux Balance Analysis) [23], CycleFreeFlux [24], the global method to remove infeasible cycles from NESS flux configurations introduced in [15]. Clearly, it is useful in practice whenever one is interested in finding a single feasible flux pattern as long as the minimum of $Q$ lies within the solution space defined by mass balance constraints. A more integrative principle has been proposed in [25]. It is most easily expressed by distinguishing forward ($F$) and reverse ($R$) directions for each flux, so that the net flux $v_{i}$ can be written as $v_{i}=v_{i,F}-v_{i,R}$ (with $v_{i,F}\geq 0$ and $v_{i,R}\geq 0$), as well as exchange fluxes corresponding to sources or sinks of the reaction network ($v_{i,E}$). In brief, it states that the triplet $(\mathbf{v}_{F}^{\star},\mathbf{v}_{R}^{\star},\mathbf{v}_{E}^{\star})$ satisfies a Flux Balance Analysis-like problem (i.e. is mass balanced and maximizes a given linear objective function) thermodynamically, provided it minimizes the functional

[TABLE]

subject to

[TABLE]

where $\mathbf{S}_{I}$ and $\mathbf{S}_{E}$ stand for the intracellular and exchange parts of the stoichiometric matrix, $\mathbf{c}$ is a generic vector in $\mathbb{R}^{N}$ and $\mathbf{1}$ is the vector with all entries equal to 1. In other terms, if $\mathbf{v}_{E}^{\star}$ is the vector of optimal exchange fluxes for the solution of an FBA problem, one can obtain a thermodynamically viable solution to the same FBA problem by minimizing $\mathcal{F}$. Importantly, the vector $\mathbf{g}$ of chemical potentials are related to the vector $\boldsymbol{\lambda}$ of Lagrange multipliers enforcing (21) by

[TABLE]

This result follows from standard convex analysis (see [25] for details) and benefits from the standard conceptual and computational advantages of convex optimization problems with linear constraints (uniqueness of solution, efficient computational implementation). In addition, it is fully generic, in the sense that each thermodynamically viable mass-balanced flux vectors must minimize $\mathbf{F}$ for a certain choice of $\mathbf{c}$. It therefore provides a rather transparent description of thermodynamic feasibility (as far as optimality is concerned). On the other hand, the existence of a free parameter constitutes a limitation, at least in part. The most serious drawback, in our view, however derives from the fact that the above formulation does not allow to account for explicit constraints on net fluxes (as discussed in [25])

Thermodynamic arguments have also inspired different types of variational principles that effectively extend the reach of flux-based models by allowing to account for concentrations. For instance, in Ref. [26] feasible NESS are assumed to minimize the function

[TABLE]

over exchange fluxes $\mathbf{u}$ and intracellular fluxes $\mathbf{v}$, subject to

[TABLE]

and with prescribed bounds of the form $u^{\mu}\in[u^{\mu}_{\min},u^{\mu}_{\max}]$ and $v_{i}\in[v_{i,\min},v_{i,\max}]$ for each extracellular compound $\mu$ and each intracellular reaction $i$. The explicit dependence of $H$ on external concentrations (considered as fixed parameters) makes it possible to use the above principle to infer intracellular reaction rates given the levels of a set of extracellular metabolites. On the other hand, the minimization of (23) does not ensure that the resulting flux pattern is thermodynamically viable, unless in specific cases. The reader is referred to [27] for details.

Conclusions

Integrating thermodynamics with genome-scale biochemical reconstructions is perhaps the central theoretical open challenge of metabolic network modeling. Basically, two classes of approaches are currently being attempted that can mutually benefit from each other. On one hand, novel empirical data and biochemical methods are employed to estimate more accurate standard free energies for reactions and chemical potentials for metabolites [28, 29]. Having better estimates of such quantities is crucial, especially if they cover a larger part of the reactome and of the metabolome than the ones currently available. On the other hand, optimization principles based on different physico-chemical arguments can be used to obtain approximate (but genome-scale) estimates, whose accuracy depends on the prior biochemical information as well as on the underlying assumptions. The key issue to be faced here is computational, relating both to the scalability of the algorithms required for the study of genome-scale networks and to the fact that thermodynamics may require the study of non-convex optimization problems. At the same time, ongoing work is uncovering how thermodynamics constrains flux patterns at steady state starting from the kinetics of individual enzymes, with results which suggest that the problem of computing thermodynamic potentials described here might have to be modified as more is known about individual reaction mechanisms, at least to some degree [30]. In this respect, it is the convergence of novel statistical physics [31, 32, 33, 34], biochemical and algorithmic ideas that will likely provide the tools to effectively tackle this challenge.