Viable flux distribution in metabolic networks

The metabolic networks are very well characterized for a large set of organisms, a unique case in within the large-scale biological networks. For this reason they provide a a very interesting framework for the construction of analytically tractable statistical mechanics models. In this paper we introduce a solvable model for the distribution of fluxes in the metabolic network. We show that the effect of the topology on the distribution of fluxes is to allow for large fluctuations of their values, a fact that should have implications on the robustness of the system.


Introduction
Dynamical models on networks have attracted a large interest because of the non-trivial effects of network structure [1,2,3,4] on the dynamics defined on them [5]. Important examples of the dynamics on networks with relevant applications are the Ising model [6,7,8], the spreading of a disease [9] and the synchronization models [10,11]. In this paper introduce a solvable model for the distribution of fluxes in the metabolic network. While motivations come from the study of the metabolic network, the problem is quite general and can be applied to supply networks and to many other linear problems [12] of constraint satisfaction on continuous variables on a network.
Metabolic networks describe the stoichiometric relations between substrates in biochemical reactions inside the cell. They have been mapped [13] for a large number of organisms in the three different domains of life (archea, bacteria and eukaryotes). They provide the biomass needed for cell duplication, and the rate of biomass production (growth rate) can be identified with a fitness of the cell. The structure of the metabolic network can be represented as a factor graph with nodes that are chemical reactions and function nodes that are chemical metabolites. The projection of the network on the metabolites has a power-law degree distribution and a hierarchical structure [14,15,16]. To each factor node, which describes a chemical reactions, it is associated an enzyme which itself is produced by a regulated gene network. Important aspect of the functioning of these very complex systems include dynamical considerations. Flux-balance-analysis [17,18,19] make a major semplification in the problem. In fact it considers only the steady state of the dynamics and includes all the dynamical terms inside the definition of the flux of a reaction. For this reason it was able to predict with sufficient accuracy the fluxes of the reactions in the graph for a given environmentand it consitute a real break-through in the field. Special interest has been addressed to the perturbation of the distribution of the fluxes after knockout of a gene or in different environments [20,21]. The problem of identifying the flux distribution in Escherichia coli was studied experimentally [22] and by means of Flux-Balance-Analysis [23] . A fat tail in their distribution with different power-law exponents α < 2 was found.
Metabolic networks provide a very interesting framework for the construction of analytically tractable models using tools of statistical mechanics of disordered systems. In this paper we will discuss the impact of the network structure (degree distributions) on the steady state distribution distribution of the fluxes. We shall consider random networks with the same degree distribution as the real ones i.e. networks in the the hidden-variable ensemble [24,25,26] with same expected degree distribution as the metabolic factor graphs. Formally the problem is resolved with replica calculations on diluted networks [7] extended to the case of continuous variables. Due the simplicity of the Hamiltonian the problem is solved with an expansion of the order parameter in terms of Gaussians. The problem shares some similarity with other problems in statistical mechanics of disordered systems [28,29]. In a recent paper [30] a similar problem was considered in the framework of a different model where the steady state of the fluxes is not a priori considered and the positive fluxes don't have any upper limit.

The model
The metabolic network has a bow tie structure [16] : the metabolites are divided into: (i) input metabolites which are provided by the environment, (ii) the output metabolites which provide the biomass and (iii) the intermediate metabolites. The stochiometric matrix is given by ((ξ µ,i )) where µ = 1, . . . , M indicates the metabolite and i = 1, . . . , N the reaction and the sign of ξ µ,i indicates if the metabolite µ is an input or output metabolite of the reaction i. As in the Flux-Balance-Analysis method we assume that each intermediate metabolite has a concentration c µ which is consumed/produced by a reaction i at a rate f i . At steady state, we have where f i is the flux of the metabolic reaction i. For the metabolites present in the environment and the metabolites giving rise to the biomass production we can fix the incoming flux given by The fluxes in Eqs.(1-2) can vary inside a fixed volume Ω. We assume for simplicity that this volume is an hypercube Ω = [0, 2L] N . Changing the variables f i in the variables s i = f i − L and the equations that the fluxes s i must satisfy are given by where g µ = a µ − L i ξ µ,i . The volume of solutions V , given the constraints (3), is proportional to the quantitỹ where we have used the heterogeneous spherical constraints and integrated over L ′ in the interval [0, L] in order to allow analytical treatment of the problem.

Replica method
We assume that the support of our stochiometric matrix is a random network with given degree distribution, i.e. a realization of the random hidden-variable model [24,25,26]. In particular we fix the expected degree distribution of the nodes of the factor graphs to be q i for the reaction node i = 1, . . . N and q µ for the metabolite nodes µ = 1, . . . , M and we assume that the matrix elements ξ µ,i are distributed following where δ() indicates the Kronecker delta. Note that in (6) we have assumed that the elements of the stochiometric matrix have values 0, ±1 with a random sign and that the variables q i q µ are nothing else than the hidden-variables associated with metabolite µ of the reaction i of the hidden-variable network ensemble [24,25,26]. In order to evaluate the steady state distribution of the fluxes in a typical network realization we replicate the realizations of the s a i and we compute log(V ) over the different network realizations. To calculate this average we use the replica trick S = log(Z) = lim n→0 . The averaged unormalized volume of solutions <Ṽ n > can be expressed as Using the techniques coming from the field of diluted systems, we introduce the order parameters [31,7] getting for the volume The saddle point equations for evaluating Σ are given by We assume that the solution of the saddle point equation is replica symmetric, i.e. the distribution of the variables z a = λ a , s a conditioned to a vector field x are identically equal distributed, where φ(z|x) are distribution functions of z and P (x) is a probability distribution of the vector field x. For the function φ(z|x) the exponential form is usually assumed in Ising models. In our continuous variable case for our quadratic problem, we assume instead that φ(z|x) has a Gaussian form. This assumption could be in general considered as an approximate solution of the equations (9). Explicitly we assume that the functions c(λ) and c(s) can be expressed as the following, from which we derive forĉ(s) andĉ(λ) The saddle point equations (9), taking into account the expression for the order parameters (11)(12) closes and can be written as recursive equation for P (h λ , m λ ) and P (h s , m s ), i.e.
Finally S can be calculated at the saddle point as

Population Dynamics
We solved the equations (13) by a population-dynamical algorithm. We represent the effective field distributions (h s , m s ) (h λ , m λ ) by a large population of P ≫ 1 fields. Running the algorithm the population is first initialized randomly and then equations (13) are used to iteratively replace the fields inside the population until convergence is reached. Instead of fixing ω we introduce a further variable Λ fixing the value of the average flux in the network. The action of the algorithm is summarized in the following pseudo code • choose a metabolite i 0 with probability q i P (q i ); • choose a random reaction µ 0 with probability q µ P (q µ ) • draw d form a Poisson distribution (e −qµ q k µ /k!) while (not converged) return end We run the population dynamics algorithm and we measure the distribution of the average fluxes m s /h s , the distribution of the fields h s for different values of Λ. We consider as the underline network a network with the real degree distribution of the metabolic factor graph of Saccharomyces cerevisiae and on a network with the same number of metabolites and reactions as the real Saccahromyces cerevisiae network but with a fixed connectivity for each metabolite and reaction node. We consider a population of P = 3000 pair of fields (h s , m s ). A random fraction of 5% of the nodes is chosen as an input/output metabolite. The values of g µ are chosen randomly depending if the metabolite µ is an intermediate metabolite or an input/output metabolite.
For the input/output metabolites we assume that g µ is a random number uniformly distributed in the interval [−100Λ, 100Λ] mimicking high rate of production/consumption. For intermediate metabolites we choose g µ with a uniform distribution in the interval [−Λ, Λ].
The distribution of m s /h s as a function of Λ are plotted in figure 1(a) for the real metabolic network of Saccharomyces cerevisiae and in figure 1(b) for the random network with two delta function degree distribution (P (q i ) = q i , P (q µ ) = q i N/M ). We observe that the average fluxes distribution m s /h s in Saccharomyces cerevisiae for low Λ has a fatter tail for the real degree distribution than for the two delta peak degree distribution.
On the other side the distribution of h s is very different in the real and in the random case (see figure 1(c)-(d)). In particular for the real metabolic network degree distribution is broader allowing with higher probability for smaller value of h s than in the case of a two delta peak degree distribution. Therefore we have shown that the real topology of the networks has as a major consequence in allowing larger fluctuations of the fluxes in the network.  Fig. 1. Distribution of the fields ms/hs and of the fields hs for a the real degree distribution of the metabolic network of Saccharomyces cerevisiae (graphs (a) and (c))and for a graph with the same number of metabolites and reactions and the same number of nodes that the real metabolic network of Saccharomyces cerevisiae but with two delta peaks for the degree distribution (graphs (b) and (d)).

Conclusions
In this paper we have proposed a statistical mechanics approach for the study of flux-balance-analysis in a particular ensemble of metabolic networks. We have studied the impact of the topology of the networks on the distribution of the fluxes. We observe that the role of the real topological structure is to allow for larger variation of the fluxes, a fact that should have implications for the robustness of the system. In particular we found that the topology of real metabolites has an impact on the fat tail of the m s /h s distribution and on the small h field of the network, Further work is under consideration for the implementation of a message-passing algorithm able to predict the fluxes taking into account the full complexity of the real metabolic network.