On the role of conserved moieties in shaping the robustness and production capabilities of reaction networks
We study a simplified, solvable model of a fully-connected metabolic network with constrained quenched disorder to mimic the conservation laws imposed by stoichiometry on chemical reactions. Within a spin-glass type of approach, we show that in presence of a conserved metabolic pool the flux state corresponding to maximal growth is stationary independently of the pool size. In addition, and at odds with the case of unconstrained networks, the volume of optimal flux configurations remains finite, indicating that the frustration imposed by stoichiometric constraints, while reducing growth capabilities, confers robustness and flexibility to the system. These results have a clear biological interpretation and provide a basic, fully analytical explanation to features recently observed in real metabolic networks.
Biological complexity Spin-glass and other random models Systems biology
Metabolic networks represent the biochemical machinery by which cells dispose of the nutrients found in their surrounding environment to produce the macromolecules needed for their survival, including nucleic acids, membranes, cell walls, proteins and the energy carrier ATP. Understanding how functionality emerges from such a complex system of biochemical reactions is a major issue with important implications in bioengineering and pharmacology. Unfortunately, detailed kinetic models of genome-wide metabolic networks are computationally prohibitive and fine tuning of parameters would require an amount of biological information that is not available. It is therefore important to develop methods to predict the organization of metabolic fluxes from the sole knowledge of the stoichiometry, which represents the full available information about the network topology and about the proportions in which different reagents are interconnected as a set of input-output relations. This problem is central in systems biology. A key question in this context is whether the physiological states observed in real cells are optimal from a metabolic viewpoint, i.e. if reaction fluxes self-organize so as to maximize the production capabilities of the network (or of a smaller set of metabolites) [1, 2, 3, 4].
It is quite intuitive that, in a generic situation where reagents flow in a heterogeneous input-output network, the ways in which fluxes can be chosen will be heavily constrained by the need to meet prescribed production objectives, while the dependence of the global properties on the specific realization of input and output coefficients will become weaker and weaker for larger networks. Ultimately, it is reasonable to think that the feasibility of flux configurations meeting all constraints, the robustness of the solutions against localized flux variations (or reaction knock-outs), and the productive efficiency of solutions will depend crucially on structural parameters like the ratio between the number of reactions (variables) and that of reagents (constraints).
But in order to address the question of optimality in a realistic biochemical network it is necessary to take into account the peculiar nature of stoichiometric coefficients, which enforce mass balance conditions at each reaction node in the network. Such relations in turn define a family of conservation laws reflecting the existence of dynamical invariants of purely topological origin that are able to affect the kinetics of the system as a whole . To be explicit, let (respectively, ) denote the output (respectively, input) stoichiometric coefficient of metabolite in reaction . A group of metabolites satisfying
is such that the total number of molecules in the pool, or the aggregate concentration, does not change in time. Such pools are ubiquitous in real metabolic networks: an important example is the total quantity of ATP and ADP (the ‘adenylate moiety’) which remain constant during metabolic activity, as ATP is continuously discharged to ADP and vice-versa .
In this Letter we study how the presence of conserved metabolic pools, reflecting the underlying mass-balance conditions of stoichiometric origin, affect the global growth capabilities and the volume of the solution space in a toy metabolic network where reagents are fully interconnected and stoichiometric coefficients (which play the role of a quenched disorder) satisfy laws such as (1). To be as general as possible and to highlight the ingredient of conserved pools, we adopt the framework of Von Neumann’s expanding model , under which one aims at computing the maximum value of for which the system of inequalities
admits solution(s) in the form of non-zero flux vectors . From a physical point of view, the conditions (2) ensure that the total output of each metabolite is at least times the total input, so that if the maximum achievable (which we shall denote by ) exceeds (respectively, is smaller than) one, the network is in an expanding (respectively, contracting) state. signals instead a stationary network. (See  for a more thorough description of Von Neumann’s setup and a dynamical derivation of (2).)
Following , we assume that every reaction produces and consumes (in different proportions) every metabolite. In such a fully-connected framework, when and the input and output matrix are independent, identically distributed quenched random variables, it has been shown that depends only on in such a way that the system undergoes a transition from a contracting to an expanding phase at . Moreover, (2) admits a unique feasible flux configuation when . The advantage of using an unrealistic fully-connected setup lies mainly in its analytic tractability. We shall see that indeed the fully connected model provides an excellent proxy for a real metabolic network, at least as long as production capabilities and solution space are concerned. Graphical versions of the model yield essentially the same scenario .
To keep things mathematically simple, we account for conserved metabolic pools by constraining each disorder sample (i.e., each realization of the stoichiometric coefficients) to embed a pool formed by a finite fraction of metabolites. In other terms, we request that
where is a quenched random variable that equals with probability and zero otherwise. In this way, we include a single conserved pool in the system. It will become clear that both the number of pools and their size do not affect the growth properties in the framework we consider.
For a start, note that multiplying each term in (2) by and summing over one easily obtains
If the ’s are such that (3) holds, then either or all ’s connected to a metabolite in the conserved pool must vanish. The latter condition however leads to the null solution in a fully connected model and must be discarded. Hence necessarily . One then sees that an expanding regime cannot occur in a system with a conserved metabolic pool, independently of its size. This suggests a radically different and more realistic scenario than the unconstrained case. In addition, it is easy to understand that the inequalities (2) must become equalities at for all metabolites belonging to a conserved pool. In other words: if a metabolite belongs to a conserved pool there can be no net production or consumption of it in the optimal state.
To get a more thorough insight, it is necessary to compute exactly as a function of . As in , the calculation can be carried out in two steps. First, we compute the typical volume of flux configurations compatible with (2) and (3). This requires the calculation of an average over the constrained quenched disorder , for which we shall resort to the replica trick . Next, following Gardner , we impose that the average distance between solutions vanishes (i.e. that the typical volume reduces to a single point). This condition intuitively marks , if an increase of reduces the set of flux vectors satisfying (2) until no more solutions are found. We will see that this Ansatz describes the system correctly only up to a critical value of . Above this point, the solution is no longer unique. The breakdown of the Ansatz is directly linked to the existence of a conserved metabolic pool.
The volume of flux vectors satisfying (2) is given by
where the extra factor (here and in (7)) is added for convenience. By analogy with known systems with quenched disorder, we expect the typical volume of solutions at fixed to be given by , where (with over-bar denoting the average over the quenched disorder)
is a self-averaging quantity. Note that the disorder-averaging includes the constraints (3), so
where is now an average over the free, un-constrained stoichiometric coefficients . In turn, the quenched average can be computed via the replica trick:
As in , it is convenient to write the stoichiometric coefficients as and , with and zero-average Gaussian random variables. Inserting these in (2) one immediately sees that, to leading order in , the optimal growth rate depends only on the average input and output coefficients: . In turn, the corrections to the leading order can be captured by re-writing as
(3), however, requires that the average input and output coefficients are the same (this can be easily seen by direct substitution). Hence we shall set and shift the focus to : (resp. ) now signals an expanding (resp. contracting) phase. The calculation in the limit can be carried out along the lines of , except that the new constraints (3) and the fact that the average flux is now a free variable (in  it was conveniently fixed to because of the invariance of (2) under re-scalings ) lead to the introduction of extra order parameters. The key one is however still the overlap
between different solutions and (at fixed ). In the replica-symmetric approximation (which is putatively exact in the present case due to the convexity of the space of solutions),
and one is lead to consider the following saddle-point problem:
Physically, the order parameter measures the distance between solutions (it is easy to see that indeed ). It is reasonable to expect that as increases, the typical volume of solutions shrinks as it gets more and more difficult to satisfy (2). Assuming that a single flux state satisfies all constraints when is then equivalent to studying the solution of (12) in the limit . This can be done by introducing a proper re-scaling of the order parameters in terms of , similarly to , so as to allow the integrals in (13) and (14) to be evaluated by steepest descent. In the present case, it is convenient to set
With these definitions, one obtains the following set of saddle point conditions:
and denote, respectively, the fraction of “intermediate” metabolites for which the total output equals the total input and the fraction of reactions with zero flux: in particular,
Equations (17) can be solved numerically to obtain (and all order parameters) as a function of . Note that the dependence of on the disorder variables is embodied in the combination , meaning that the optimal growth rate is larger the bigger is the spread in the difference between input and output coefficients (in other words, the system maximizes growth by taking advantage of imbalances between inputs and outputs in the stoichiometric coefficients). Results for as a function of are shown in Fig. 1.
The critical point where (now a function of ) becomes larger as increases and ultimately, when all metabolites form a conserved pool, becomes equal to two. This signals that conserved metabolic pools reduce the optimal growth capabilities, e.g. twice as many reactions are required to ensure a steady state with when than when . Similarly, the quantities and defined in (20) and (21) are displayed in Fig. 2.
Notice that as increases a larger and larger fraction of reagents is fully re-cycled into production, while the fraction of zero fluxes becomes consistently smaller as increases (as more active reactions are required to keep a steady growth rate).
These results can be confirmed numerically by calculating via the MinOver algorithm introduced in  based on  (see Fig. 3). Its main idea is to rotate the vector iteratively in the direction of the worst-satisfied constraint, until all constraints are satisfied at fixed . Then one can increase by a quantity and iterate the procedure, until no more solutions are found.
Simulations suggest that as the solution is unique for (where ), while for the assumption of a unique solution that underlies the limit ceases to be correct and multiple solutions must exist. This can be seen directly in numerical experiments by measuring the overlap
between different solutions and obtained by MinOver. If , then . In general, . (Note that the definition (22) must be corrected to deal with null fluxes, since the overlap of the null solution with itself must obviously equal one.) In particular, we are interested in the average of (averaged over many solutions), denoted by . This is reported in Fig. 4.
It is natural to expect that as increases, the volume of solutions decreases and correspondingly increases. We see that indeed the volume initially shrinks but it stops contracting before the optimal growth rate is reached. Ultimately, at , confirming that many flux states satisfy (2). By contrast, when no conserved pools are present keeps growing and tends to as , as a unique solution exists.
From a biological viewpoint, one sees that stoichiometric constraints on one hand reduce the production capabilities limiting the system to a stationary state. On the other hand, they increase robustness, since many (microscopic) flux states are compatible with the optimal (macroscopic) growth performance. Clearly, increased flexibility is crucial for biological stability, since environmental changes or localized flux variations are more likely to be sustained by the system by re-arranging fluxes while keeping maximal production rates. This trade-off between stability and optimal growth has been recently observed in the metabolic network of the bacterium E. coli : a finite volume of flux vectors is indeed compatible with a maximal growth assumption for E. coli and the predicted flux range is in good agreement with the measured experimental fluxes in different environments.
It is easy to extend these results to the case in which the reaction network presents dispersion at reaction nodes. The simplest possibility is to consider, instead of (3), the condition
with ’s real constants. A positive (respectively, negative) says that for reaction more (respectively, less) reagents are produced than they are consumed. While the typical properties of large, random instances with random ’s can be characterized analytically along the lines described above, it is more useful to concentrate on the case where , which provides for an immediate interpretation of the limiting cases and . It is easy to understand that (23) implies that the average input and output stoichiometric coefficients must be linked by
This in turn yields
Hence, the influence of is essentially to allow, in a finite system, for an expanding phase if . If , instead, the system is necessarily confined to a contracting regime. In particular, the introduction of does not affect . Analytical and numerical studies fully confirm this prediction (see Fig. 5).
To summarize, the existence of conserved metabolic pools (groups of reagents whose aggregate concentration is invariant in time) stems directly from stoichiometric balance and characterizes all biochemical reaction networks. The model we addressed aims at understanding how these invariant structures affect the global growth properties within a simple exactly solvable setting. We have seen that stoichiometry frustrates the system by reducing its optimal production capability. At the same time, a finite volume of flux configurations is compatible with the optimal conditions, implying an increase of stability and flexibility. Both ingredients are crucial at the biological level since, reasonably, metabolic networks should be robust to local flux variations. The basic mechanism by which local conservation laws build up to regulate the trade-off between growth and stability in the present model is likely to lie at the core of recently observed properties of real metabolic networks. An important aspect that cannot be analyzed theoretically at the level of a fully-connected model is how dynamical invariants affect the way in which reaction networks re-arrange fluxes in response to knock-outs or environmental changes. It is reasonable to expect that the graphical version of the problem, introduced in , will be useful to shed some light on this issue.
Acknowledgements.We are indebted with S. Franz and T. Jörg for useful comments and questions.
- Present address: Department of Molecular, Cellular and Developmental Biology, Yale University, KBT 1048, P.O. Box 208103 New Haven, CT 06520 (USA)
- Present address: Department of Mathematics, King’s College London, Strand, London WC2R 2LS (UK)
- \NameA. Varma B. O. Palsson \REVIEWNature Biotechnology121994994.
- \NameD. Segrè, D. Vitkup G. Church \REVIEWProc. Nat. Acad. Sci. USA99200215112.
- \NameE. Almaas, B. Kovacs, T. Vicsek, Z. N. Oltvai A.-L. Barabasi \REVIEWNature4272004839.
- \NameR. Schuetz, L. Kuepfer U. Sauer \REVIEWMolecular Systems Biology32007119.
- \NameI. Famili B. O. Palsson \REVIEWBiophys. J.85200316.
- \NameJ. Von Neumann \REVIEWRev. Econ. Stud.1319451.
- \NameA. De Martino M. Marsili \REVIEWJ. Stat. Mech.2005L09003.
- \NameA. De Martino, C. Martelli, R. Monasson I. Perez Castillo \REVIEWJ. Stat. Mech.2007P05012.
- \NameM. Mezard, G. Parisi M. A. Virasoro \BookSpin glass theory and beyond \PublWorld Scientific, Singapore \Year1987.
- \NameE. Gardner \REVIEWJ. Phys. A: Math. Gen.211988257.
- \NameW. Krauth M. Mezard \REVIEWJ. Phys. A: Math. Gen.201987L745.
- \NameC. Martelli, A. De Martino, E. Marinari, M. Marsili I. Perez Castillo \REVIEWsubmitted2008.