The free lunch of a scale-free metabolism

TL;DR

This paper demonstrates that apparent scale-free patterns in metabolic fluxes are artifacts caused by thermodynamically unfeasible cycles, and correcting for these reveals a more meaningful, constrained metabolic space with complex structure.

Contribution

It identifies the source of scale-free tails in metabolic flux data as thermodynamic artifacts and proposes a correction method revealing the true constrained metabolic space.

Findings

Scale-free tails are artifacts of unfeasible cycles.

Correcting thermodynamics reveals a constrained, multimodal metabolic space.

Metabolic space scales with physical limits after correction.

Abstract

In this work it is shown that scale free tails in metabolic flux distributions inferred from realistic large scale models can be simply an artefact due to reactions involved in thermodynamically unfeasible cycles, that are unbounded by physical constraints and would be able to perform work without expenditure of free energy. After correcting for thermodynamics, the metabolic space scales meaningfully with the physical limiting factors, acquiring in turn a richer multimodal structure potentially leading to symmetry breaking while optimizing for objective functions.

Equations

$${\dot{c}}_{\mu }=\sum _i{S}_{i\mu }{f}_i$$

$$\sum _i{S}_{i\mu }{f}_i=0$$

$$f_i\in \left[f_{i,min},f_{i,max}\right]$$

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

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

$$u_i\ge 0$$

$$\sum _i{u}_i\le U$$

$$u_i=y_i+x_{i+1}-x_i$$

$$y_i\ge 0$$

$$V(\{u_i\})={\int }_{x_{i+1}-x_i\le u_i}\prod d{x}_i(\text{Not all the same sign})$$

$$-r_1-r_2\le u_1$$

$$r_2-r_3\le u_2$$

$$r_3+r_1\le u_3$$

$$I={\int }_0^{u_3}(u_3-r_3)(u_2+r_3)d{r}_3=\frac{u_3^2}{2}(u_2+u_3/3)$$

$$2V(u_1,u_2,u_3)=U\sum _i{u}_i^2-\frac{1}{3}\sum _i{u}_i^3$$

Full text

Key words: Thermodynamics, Metabolism, Physical Scales.

Introduction

How physical constraints shape the space of feasible living beings? Could we have a bacterium 1 meter long and/or men as tall as buildings? Typically natural systems are characterized by certain scales, e.g. it is well known from statistical physics that we can break a magnet and get two downto a certain size, known as correlation length[1]. Scale-free systems however are present in nature, the proven example being equilibrium thermodynamic systems at criticality, ruled by a beautiful and deep physical theory[2]. The wide and rich response to perturbations that these systems endow as well as intriguing observations[3] have raised the quest for scale free processes and structures in biology[4] with interesting ongoing discussions, e.g. on biological neural networks structure and dynamics[5, 6, 7]. A key biological process with deep physical roots ideally at the core of emergence of scales in living systems is metabolism, i.e. the fundamental set of reactions building up cell components and transducing free energy. The energetic constraints of metabolism have been considered to play a role in many important issues, from eucaryotes evolution to aging, all the way down to the origin of life itself[8]. Nowadays it is possible to keep track for a given organism the entire set of enzymes devoted to metabolic functions upto to reconstruct the entire chemical reaction network at genome scale level. The analysis of such complex systems is difficult, simplifying assumption are maded into a framework known as constraint based modelling and/or flux balance analysis[9] with wide applications e.g. in metabolic engineering[10]. Apart from practical applications are there general laws ruling the large scale organization of metabolic networks? Based on general inferences on these models, it has been argued within the framework of complex network theory[11] that metabolic fluxes are globally organized in a scale free manner around a backbone of reactions carrying high-intensity fluxes[12]. While it is generally acknowledged the existence and importance of central core pathways in metabolism, it is generally believed that such core has however scales that are simply given by the physical limiting factors, being them extrinsic (resources avalability[13]) or intrinsic (maximum ribosome elongation rate[14]). Further, it is generally known in the field of flux balance analysis that high intensity fluxes retrieved in models can be a spurious effect due to an erroneous handling of thermodynamics constraints. The way these constraints impact and shape the space of feasible states[15] is thus very important and attracted a lot of attention in recent years[16] upto the deeper clarification of its theoretical basis from a mathematical[17] and physical point of view[18]. In this work it is showed that once thermodynamic constraints are implemented the long tail in flux intensity distribution is truncated and the whole space scales meaningfully with the physical limiting factors, e.g. the availability of resources and/or the dilution rate. On the other hand it is pointed out that the space is not anymore convex thus acquiring a rich multimodal structure possibly leading to symmetry breaking upon optimizing for objective functions. After illustrating the background, we analyze in particular a genome scale reconstruction of the metabolism of the bacterium E. Coli and illustrate the mechanisms in simple toy models.

Results

Background.

The assumption behind constraints based modelling is essentially to consider metabolism as a well mixed chemical reaction network in the steady state. Concentration variations of metabolites are linearly related by mass balance to reactions fluxes through the matrix of stoichiometric coefficents

[TABLE]

Once a steady state is assumed and reaction fluxes are bounded by physical limits (reversibility, medium, kinetics) the reactions fluxes span a convex polytope

[TABLE]

When reaction bounds are not known/provided, it is costumary to set them to an arbitrary large number, e.g. $f_{i}\in\left[-10^{4},10^{4}\right]$, that, for a meaningful model, shall not influence the results. In the most simple way thermodynamics is implemented by setting reaction reversibility, i.e. $f_{i}\geq 0$ for some reactions. On the other hand a more rigorous yet simple approach consists in postulating that fluxes shall follow at least a free energy gradient[19], i.e. if $f_{i}\neq 0$, then $f_{i}\Delta G_{i}<0$ where $\Delta G_{i}$ is the free energy drop of reaction $i$. Upon decomposing the reaction free energy in terms of chemical potentials we have a system of linear inequalities ($\xi_{i\mu}=-sign(f_{i})S_{i\mu}$)

[TABLE]

whose feasibility is necessary for the thermodynamic soundness of the flux configuration. The feasibility of such a system is ruled by the Gordan theorem by which the dual system

[TABLE]

shall have no non-trivial solution(s). It is thus seen that any solution of the latter system, that has the form of a closed cycle, it provides an unfeasible “gauge invariance”[20] for the involved fluxes such that, if $f_{i,0}$ is a feasible state and $k_{i}$ is such a solution, the line $f_{i}(L)=f_{i,0}+k_{i}s_{i}L$ is still inside the polytope till $L$ reaches the (arbitrary) bounds. From a physical point of view it corresponds to a “perpetuum mobile” from which we could get work for free. A simple illustration is provided in the toy network depicted in Fig 1, where we see that the variable $x$ reaches unfeasibly high values without taking into account properly thermodynamics, whose implementation leads on the other hand to the non linear constraints $ux\geq 0$, $|x|\leq|u|$. In essence thermodynamics forbid certain flux directions cutting in turn orthants within the feasible space that get decomposed in a connected star of convex polytopes, that is overall not convex anymore (see Fig 1), but whose scales now shall reflect the truly physical constraints of the given model. If we suppose that the flux $x$ is performing work we would have perpetual motion outside the thermodynamically feasible space (a “realistic” example for ATP production has been shown in [19]).

E. Coli genome scale model analysis.

We consider here the typical steady states of the genome scale model of the E Coli metabolic network iJR904[21] in a glucose limited minimal medium in aerobic conditions with maximum glucose uptake $u=6$mmol/gDWh (the default case provided with the model). After removing leaves this space has dimension $D=233$ for $N=666$ reactions. The uniform sampling of the underlying convex polytope has been performed by an hit-and-run Montecarlo Markov chain[22], whose convergence is guaranteed in feasible times[23] upon handling ill-conditioning[24] (that can be severe exactly because of thermodynamic unfeasibilities). Thermodynamically unfeasible cycles have been enumerated for this network[19] and they can be detected and removed in many ways[25, 26, 27]. They have been corrected in a minimal way by adjusting the gauge parameter $L$ in turn to silence the less intense reaction flux that amounts to a projection onto the thermodynamically feasible space[20]. Results for the overall distribution of flux instensities $|f_{i}|$ are shown in Fig 2, before and after correcting for thermodynamic unfeasibilities.

The long tail corresponding to the uniform sampling depends on the arbitrary constant fixed for flux bounds (in this case $10^{3}$) and it can be extended thus arbitrarely. In regard to single reaction flux distributions, we show in Fig 2 (bottom) the retrieved histogram of adenylate kinase, an important enzyme involved in energy homeostasis performing the reaction $2ADP\leftrightarrow ATP+AMP$. Since this reaction is involved in an unifeasible cycle[19], its distribution results to be flat and wide upon uniform sampling, while after thermodynamic correction shows a bimodal fashion opening possible speculations about the presence of two effective states for energetic homeostasis. In order to keep track in a more rigorous way the scales of the metabolic space in dependence of physical limiting factors we performed a principle component analysis of the sampled configuration upon varying uptake conditions. The square root of the retrieved eigenvalues are proportional to the diameters of the maximum inscribed ellipsoid into the convex hull of the space and thus quantify more rigorously its scales. We show in Fig 3 (left) the $100$th largest eigenevalues before and after thermodynamic correction. At odds with the interesting case of assets cross-correlations in financial markets[28], here the largest eigenvalues for the uniform sampling correspond to unfeasible cycles that desapper after thermodynamic correction. In fig 3 (right) it is shown that the two highest eigenvalues scale linearly with the maximum glucose uptake after taking into account thermodynamics, while they do not in the infeasible case, that is indeed scale free.

Towards symmetry breaking in toy models.

We consider the toy model in Fig 4, where three units uptake some nutrients with intensity $u_{i}$: they either use it ($y_{i}$) or exchange it ($x_{i}$) and resources are globally limited by a certain constant $U$. We have the constraints (where we identify $4=1$)

[TABLE]

Because of thermodynamic constraints, the $x_{i}$ cannot have the same sign. Suppose that $\sum_{i}y_{i}=\sum u_{i}$ is maximized, this will be equal to $U$. This leads to some degeneracy and each configuration $u_{i}$ has a phase space volume (to which probability is proportional in case of a flat prior)

[TABLE]

For three nodes we can evaluate the integral dividing it in the $6$ integrals in which the sign $s_{i}$ of the $x_{i}$ is fixed (we exclude the $2$ configurations in which signs are equal). For instance consider the integral $I$ upon fixing $s_{1}=+1$ and $s_{2}=s_{3}=-1$. We replace the $x_{i}$ with positive definite variables $r_{i}\geq 0$: $x_{i}=s_{i}r_{i}$. We have the constraints (from $x_{i+1}\leq x_{i}+u_{i}$)

[TABLE]

the first constraint is trivially satisfied, we can integrate $0\leq r_{1}\leq u_{3}-r_{3}$ and $0\leq r_{2}\leq u_{2}+r_{3}$ obtaining

[TABLE]

Considering the other cases that can be obtained by symmetry we have finally after some calculations

[TABLE]

This function has a minimum for the symmetric solution $u_{i}=U/3$ and three (equal) maxima where one unit is uptaking all resources, e.g. $u_{1}=U$ and $u_{2}=u_{3}=0$.

Conclusions

Although scientists have established to be impossible under the laws of physics, perpetual motion continues to capture the imagination of inventors and to enter in a subtle way in scientific research as well. In this work it has been shown that in order to meaningfully assess scales in metabolic networks it would be necessary to take into account thermodynamics, given the free energy-transducing nature of these sytems. It has beeen shown in one realistic example that the long tails inferred in fluxes distributions in metabolic network models are a spurious effect due to reactions involved in infeasible cycles. From the positive side, such an energy balance analysis can be rewarding from several aspects, since, apart from assessing metabolites concentrations[29], free energies[30], potentially regulated sites[31], it reveals that the metabolic space has a more complex multimodal structure. The additional layer of complexity provided by thermodynamics could enrich the possibility of biophysical modelling especially in presence of heterogeneity, cell-to-cell interactions and/or metabolic shuttling, and from a theroretical viewpoint could potentially lead to new variational principles[32]. A simple toy model has been provided to illustrate a possible mechanism of symmetry breaking while maximizing for an objective function. On the other hand the inference from such constraints is challenging from a computational point of view given the non-convexity of the space[33]. Suitable extensions of fast approximate analytical methods for the convex case[34, 35] could be very promising in this respect.

Acknowledgments

The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme ($FP7/2007-2013$) under REA grant agreement $n[291734]$.