Estimating the size of the solution space of metabolic networks

Abstract

Background: Cellular metabolism is one of the most investigated system of biological interactions. While the topological nature of individual reactions and pathways in the network is quite well understood there is still a lack of comprehension regarding the global functional behavior of the system. In the last few years flux-balance analysis (FBA) has been the most successful and widely used technique for studying metabolism at system level. This method strongly relies on the hypothesis that the organism maximizes an objective function. However only under very specific biological conditions (e.g. maximization of biomass for E. coli in reach nutrient medium) the cell seems to obey such optimization law. A more refined analysis not assuming extremization remains an elusive task for large metabolic systems due to algorithmic limitations.

Results: In this work we propose a novel algorithmic strategy that provides an efficient characterization of the whole set of stable fluxes compatible with the metabolic constraints. Using a technique derived from the fields of statistical physics and information theory we designed a message-passing algorithm to estimate the size of the affine space containing all possible steady-state flux distributions of metabolic networks. The algorithm, based on the well known Bethe approximation, can be used to approximately compute the volume of a non full-dimensional convex polytope in high dimensions. We first compare the accuracy of the predictions with an exact algorithm on small random metabolic networks. We also verify that the predictions of the algorithm match closely those of Monte Carlo based methods in the case of the Red Blood Cell metabolic network. Then we test the effect of gene knock-outs on the size of the solution space in the case of E. coli central metabolism. Finally we analyze the statistical properties of the average fluxes of the reactions in the E. coli metabolic network.

Conclusion: We propose a novel efficient distributed algorithmic strategy to estimate the size and shape of the affine space of a non full-dimensional convex polytope in high dimensions. The method is shown to obtain, quantitatively and qualitatively compatible results with the ones of standard algorithms (where this comparison is possible) being still efficient on the analysis of large biological systems, where exact deterministic methods experience an explosion in algorithmic time. The algorithm we propose can be considered as an alternative to Monte Carlo sampling methods.

Figure 1

Embedding. Sketch of the the parameterization ϕ = ℝ^{N - M} → V ⊂ ℝ^N, in a simple case where N = 3 and M = 1. The mass balance equations defined in Equation 2 define the hyperplane V and the set of inequalities defined in Equation 3 restrict V to polytope Π indicated as the grey triangle in the left part of the figure. The application ϕ map the full-dimensional counter-image of Π (the grey triangle on the right indicate with U) onto Π.

Figure 2

Marginal probability. Geometrical interpretation of the marginal probability distributions: the marginal probability at point ν_i = $\overline{\nu }$ is proportional to length of the blue segment which is the intersection between the polytope Π (the grey triangle) and the plane ν_i = $\overline{\nu }$.

Figure 3

Tiling. Discretization of the volume: in the left part of the figure we display the the regular orthogonal grid Λ_ε of side ε partitioning ℝ^N. The counter-image of Λ_ε via ϕ is given by the the grid Γ_ε . Let number of Γ_ε squares intersecting ϕ^-1(Π) is proportional to the number of Λ_ε cubes intersecting Π. The smaller the ε the better the intersection of the grid Λ_ε with the polytope Π will approximate Π.

Figure 4

Discrepance histogram. Histograms of δS = (S_BP - S_LRS)/S_LRS, where S_LRS and S_BP are the estimates of the logarithm of the volume computed using respectively the program LRS [8], and our message-passing algorithm. The measurement if δS are taken over a set of 1000 random realizations of a N = 12 and M = 4 stoichiometric matrix where each of the metabolites participates to K = 3 different reactions. The four histograms represent the distribution of the experimental frequencies of the measured δS. Measurements were taken at different values of of the discretization parameter q^max = 16, 64, 256, 1024. For larger values of q the distributions of the measured δ_S peak around zero.

Figure 5

Running time. Logarithm of the running time τ expressed in seconds vs. N for LRS algorithm, and BP algorithm. Averages are taken over 5000 realizations of random stoichiometric matrices at each value of N except in the N = 12 case where we analyzed 500 realizations. LRS outperforms BP up to sizes N ~12 where the running time of LRS explodes exponentially while BP maintains a modest almost-linear behavior.

Figure 6

Distribution of fluxes in Blood Cells. Marginal probability distributions of the flux values for each of the 46 reaction in the red blood cell network computed using our message-passing algorithm (filled area) and the MC method (lines). Panels are arranged following the same sequence of figure 5 in [10]. The maximum values of the fluxes are all identical, ${\nu }_i^{\max }$ = 1.

Figure 7

Distribution of fluxes in Blood Cells. Marginal probability distributions of the flux values for each of the 46 reaction in the red blood cell network computed using our message-passing algorithm (filled area) and the MC method (lines). Panels are arranged following the same sequence of figure 5 in [10]. The maximum values of the fluxes were taken from [10].

Figure 8

Flux knock-out curves. Flux knock-out impact on the volume of the space of solutions in E. coli central metabolism. On the x-axis we display the percentage of reduction of any given flux, and on the y-axis the relative volume difference with respect to the unperturbed system. We use upper-case keys for internal fluxes, and lower-case to exchange fluxes. We indicate keys only for the 20 fluxes having larger impact on the volume (dots) and we display the rest fluxes with thin scattered lines.

Figure 9

Top impact flux knock-outs. Impact of the 20 fluxes which have a larger impact on the volume of the space of solution. Here S_KO = S₀ - S_KO (${\nu }_i^{\max }$ = 0)measures the difference of the volume of the unperturbed system and that modified by setting ${\nu }_i^{\max }$ = 0. In the x-axis we indicate the relative flux name using upper-case keys for internal fluxes, and lower-case to exchange fluxes.

Figure 10

Correlation of impact and value. Change in entropy S₀ - S_KO after reaction knockouts and the average as a function of the average flux value ⟨ν⟩ of the unperturbed system. Red dots are relative to 75% knock-out and blue one to complete knock-out. The six red dots with average flux ~0.3 and entropy change larger than 0.15 are G1PP, GLCP, HEX1, CS, PDH, glycogen which are key check-points of central metabolism.

Figure 11

Distribution of average fluxes. Distribution of average fluxes (⟨ν⟩) of the reactions for the E. coli metabolism (points) obtained from the marginal probability distributions computed using our message-passing algorithm. We also display a power-law fit function P_fit(ν) = A/(ν + ν 0)γ (solid line) with ν₀ = 0.0002 and γ = 1.48.

Figure 12

Integrated distribution of average fluxes.Integrated distribution P_<of average fluxes ⟨ν⟩ of reactions in E. coli metabolism obtained from the marginal probability distributions computed using our message-passing algorithm. Note the jumps for ⟨ν⟩ = 0.4, ⟨ν⟩ = 0.5 and ⟨ν⟩ = 0.6.

Figure 13

Linear cavity summation. An example showing the summation strategy for a metabolite with 8 fluxes. Once all elements of the table are computed (6 convolutions are needed), one needs only 3 more convolutions to compute each one of the 8 cavity messages. In the figure the table elements that are needed to compute ⊗_{j ≠ 5} h_j = (h₁ ⊗ h₂ ⊗ h₃ ⊗ h₄)⊗ (h₇ ⊗ h₈) ⊗ h₆ are colored in blue.