Combining graph and flux-based structures to decipher phenotypic essential metabolites within metabolic networks

Metabolic network

A metabolic network is commonly represented as a directed bipartite graph (R∪M, E), where R and M are sets of nodes standing for reactions and metabolites, respectively (see Fig. 1 for illustration). For any r ∈ R, we define rcts(r) = {m ∈ M∣(m, r) ∈ E} and prds(r) = {m ∈ M∣(r, m) ∈ E}. In other words, when (m, r) ∈ E or (r, m) ∈ E for m ∈ M and r ∈ R, the metabolite m is called a reactant or product of reaction r, respectively. Quantitatively, reactions are subject to stoichiometry for balancing the relative quantities of reactants and products. This can be captured by an edge labeling giving the stoichiometric coefficient of a reaction’s reactants and products, viz. s:E → ℚ, respectively. We designate (R∪M, E, s) a stoichiometric metabolic network.

We denote by $\mathcal{R}ev\left(M\right)$ the set of reversible reactions of the system. It corresponds to the set of reactions which can be associated with another reaction of the system with reverted reactants, products and stoichiometry. Formally, it is defined as follows: $\mathcal{R}ev\left(M\right)=\left\{r\in R\mid \exists r_1\in R,\forall m\in \mathit{rcts}\left(r\right),\left(r_1,m\right)\in E\text{with}s\left(m,r\right)=s\left(r_1,m\right)\text{and}\forall m\in \mathit{prds}\left(r\right),\left(m,r_1\right)\in E\text{with}s\left(r,m\right)=s\left(m,r_1\right)\right\}$. See Fig. 1 for a reversible reaction illustration.

Each reaction r ∈ R is associated with a metabolic flux value (or activity rate) v_r, a real variable confined by a upper bound ub_r ∈ ℝ⁺.

(1)${\displaystyle 0\le v_r\le {\mathit{ub}}_r\forall r\in R.}$

Inputs and outputs of a metabolic network: set of seeds, boundary seeds and targeted reaction

We introduce the set S⊆M to model initiation seeds or nutrients, a set of metabolic compounds that are observed to be present in the system in its initial state. In addition, the graph structure includes a specific set of compounds whose production is intrinsically assumed to be activated by default. These boundary compounds or boundary seeds are defined as follows: S_b(G) = {m ∈ prds(r)|∃r ∈ R, rcts(r) = ∅}⊆M. For the sake of clarity, we assume herein all boundary compounds to be seeds: S_b(G) ⊂ S.

Its targeted metabolic compounds are hence assumed to be produced by the system. These are defined to be the reactants rcts(r_T) of the target reaction (see Fig. 1).

Stoichiometry-based production of a targeted compound

The concept of activated products (or reactions) can be modeled on the basis of different paradigms. According to the state-of-the-art, the most widely used formalism for producibility is Flux Balance Analysis (FBA) (Orth, Thiele & Palsson, 2010).

In this paradigm, each reaction r is associated with a metabolic flux value v_r, a real variable confined by Eq. (1). Flux distributions are formalized in terms of a system of equations relying on the stoichiometric coefficients of reactions. Reaction rates are governed by the law of mass conservation which assumes a steady state, i.e., input and output rates of reactions consuming and producing a metabolite are balanced. (2)${\displaystyle \sum _{\begin{array}{c}\left(r,m\right)\in E\hfill \end{array}}s\left(r,m\right)\cdot v_r-\sum _{\begin{array}{c}\left(m,r\right)\in E\hfill \end{array}}s\left(m,r\right)\cdot v_r=0\forall m\in M.}$

Worth noticing that from this definition, stoichiometrically activated reactions depend heavily on boundary compounds for which (2) is always satisfied and allows the values of fluxes to be initiated. For illustration, in Fig. 1, assuming that ub(r_S1) > 0, ub(r_S2) > 0 and ub(r_SH) > 0 implies that all reactions, including r_T, are stoichiometrically activated. On the contrary, in the absence of S₂ (i.e, ub(r_S2) = 0), all reactions are stoichiometrically inactivated except r₇ and r₈. Indeed, applying the mass conservation law to B and S₂, all flux distributions must satisfy v₀ = v₁ + v₇, v₁ = v₀ + v₈ with v₇, v₈ ≥ 0, which implies v₇ = v₈ = 0. In addition, the stoichiometry associated with the edge from r₄ ensures that the cycle r₂- r₃- r₄ is activated and that it produces extra quantities of C to feed the reaction r_T. This stoichiometry over the reaction r₄ is crucial to ensure that the cycle r₂, r₃, r₄ is not thermodynamically infeasible (Maranas & Zomorrodi, 2016), that is, a set of reactions that do not loop wth other entering or leaving metabolites. In addition, the role of compounds H and J in the reaction r₂ is to provide an external input of matter to the cycle and avoid a second type of thermodynamic infeasibility.

Graph-based (topological) production of a targeted reaction

An alternative to stoichiometry-based modeling for activated reactions is graph-based modeling, by figuring out how metabolites can be transformed by chains of reactions using the graph topology. In this area a fundamental issue is to elucidate how cycles are activated according to the graph topology (De Figueiredo et al., 2009; Schaub & Thiele, 2009; Acuna et al., 2012). In the following, we focus on the most stringent semantics for graph-based production, originally introduced in Handorf, Ebenhoh & Heinrich (2005). These semantics were first shown to be of interest for identifying important metabolites in the context of network evolution (Raymond & Segrè, 2006; Goldford et al., 2017), and for inferring minimal nutrient requirements (Handorf et al., 2008). They have also also been used in the reconstruction of metabolic networks for non-model organisms (Romero & Karp, 2001; Handorf, Ebenhoh & Heinrich, 2005; Prigent et al., 2017; Frioux et al., 2017).

Because the semantics associated with reaction activation are based on graph topology, their definition is recursive. Given a metabolic network G, a reaction r ∈ R is topologically activated from a set of seeds S if all reactants in rcts(r), that is, the predecessors of r in the graph, are reachable from S. To deploy such a recursive definition, we say that, a metabolite m ∈ M is topologically activated from S if m ∈ S or if m ∈ prds(r) for some reaction r ∈ R, where all m′ ∈ rcts(r) are topologically activated from S. In other words, at least one predecessor of m in the graph is itself activated from S. The scope of S corresponds to the set of metabolites activated from the initial seeds of the network. Therefore, it is nothing other than the closure of the set of seeds on the hypergraph, as shown in several theoretical studies (Dittrich & Di Fenizio, 2007; Cottret et al., 2008).

More formally, the concept of topological activation is defined as follows.

According to the toy example in Fig. 1, all metabolites are topologically activated from the set of seeds {S₁, S₂, H}, implying that r_T is topologically activated. Notice, however, that this property is not valid any more if we consider the seeds S₁ and S₂ independently. The set of reactions that are topologically activated from S₁ is {r₆, r₉, r_e} and the set of metabolites that are topologically activated from S₁ is {S₁, A, G, F}. Indeed, the only reaction topologically activated from the seed S₁ is r₆. Therefore, G and A are topologically activated from S₁. By recursivity, we get that r₉, hence F and r_e are activated. On the contrary, C cannot be activated from S₁ only, implying that r_T is not activated. Similarly, the sets of reactions and metabolites that are topologically activated from {S₂, H} are {r₀, r₁, r₂, r₃, r₄, r₇, r₈} and {S₂, B, C, D, E, H, J}. Therefore, r_T is not topologically activated from {S₂}. Assuming that all S₁, S₂ and H are seeds, all reactions and compounds are topologically activated so that r_T is topologically activated from {S₁, S₂, H}.

Structural properties of the set of phenotypic essential metabolites

All highly connected compounds are either PEMs, biomass components or seeds

Generic graph-based approaches usually point out the importance of hubs, that is, nodes with a high level of connectivity (Liberal & Pinney, 2013). In our context, the degree of connectivity of a metabolic compound m in a metabolic network (R∪M, E) was defined as connectivity(m) = card(r ∈ R∣m ∈ rcts(r)∪prds(r)). It depicts the number of reactions which either consume or produce the metabolite m. However, we notice that this definition does not take into account the functioning role of the considered compound with respect to the metabolic network, that is, its impact on the production of targeted metabolites. As a first analysis, we ordered the compounds of each metabolic network according to both their degree of connectivity and their role in the functioning of the network classified as PEM, biomass component, seed (that is, roughly, general system inputs from the extracellular compartment) or generic compounds. Our analysis confirmed that, as expected, regardless the considered network, all highly connected compounds are either a PEM, a biomass component or a seed (see Supplemental Information 1). Indeed, in the iJR904, iJO1366 and iAF1260 networks, the 14 compounds with the highest degrees of connectivity (from 70 to 1040) have such a role. The only exception is the periplasmic H2O of the iJO1366 network, a component which appears to play a redundant role in the network with cytoplasmic H2O. In the A. ferrooxidans and T. lutea networks, the 31 compounds with the highest degrees of connectivity (from 10 to 310) have such a role. We noticed however that the PEM concept encompasses a much larger set of compounds than the single family of highly connected compounds. This emphasizes that, in addition to the expected highly connected nodes, the PEM concept sheds light on metabolic compounds with a possibly impactful role in the production pathways of target compounds regardless their connectivity.

Putative role of sustainability, producibility and optimal-efficiency-PEMs

The distribution of sustainability, producibility and optimal-efficiency-PEMs is described in Table 1. According to this study, the most frequent PEM is the optimal-efficiency-PEM: depending on the metabolic network studied, from 77% to 99% of PEMs are optimal-efficiency-PEMs. Accordingly, for all networks, a small proportion of PEMs are not optimal-efficiency-PEMs. The number of sustainability and producibility-PEMs is fairly comparable: the ratio of producibility-PEMs ranges from 44.7% (iJO1366 network) to 99% (A. ferrooxidans network), whereas the ratio of sustainability-PEMs ranges from 46% (iJO1366 network) to 86% (T. lutea network). Similarly as before, for all networks, some PEMs are neither sustainability nor producibility-PEMs. This analysis confirms that the three classes of PEMs are complementary and must be considered together to capture the full complexity of the functioning of the network.

Overlaps between the different classes of compounds for the six metabolic networks are shown in Fig. 3. At least 26% (and up to 93%) of the network PEMs are sustainability, producibility and optimal-efficiency-PEMs simultaneously. Comparing compounds which are sustainability, producibility and optimal-efficiency-PEMs in the iJR904, iJO1366 and iAF1260 networks highlights the fact that 20 metabolic compounds satisfy these three properties in all the three networks, including the highly connected nodes udp, h, phosphate and phosphoenolpyruvate. These four compounds can therefore be viewed as the major biomass production regulators of E. coli networks in the graph, stoichiometry and optimal flux-based formalisms. Surprisingly, the other 16 PEMs in the common skeleton of the iJR904, iJO1366 and iAF1260 networks have a low degree of connectivity in the network (from 2 to 7). These compounds are Chorismate, DTDP, ahdt, Dihydroneopterin, Dephospho-CoA, 1,2-Dihydronaphthalene-1,2-diol, 6-hydroxymethyl-dihydropterin pyrophosphate, 4-amino-4-deoxychorismate, Dihydroneopterin monophosphate, 6-hydroxymethyl dihydropterin, N-((R)-4-Phosphopantothenoyl)-L-cysteine, 6-hydroxymethyl dihydropterin, 5-O-(1-Carboxyvinyl)-3-phosphoshikimate, D-4′-Phosphopantothenate, Pantetheine 4′-phosphate. One interpretation is that these metabolites may refer to structural compounds of the network necessary to produce the biomass and may constitute a skeleton of the network structure regarding the biomass production.

Network redundancies

Sustainability and producibility-PEMs that are not optimal-efficiency-PEMs depict network redundancies.

This situation occurs when a compound is required to produce a biomass component from the stoichiometry and graph-based viewpoints. All the reactions consuming this component can be removed without impacting the optimal biomass growth rate. One interpretation of this is that this type of compounds is the substrate of two pathways, both equally capable of producing the targeted biomass component. An example is shown in Fig. 4. In the A. ferrooxidans network, the compound g1p is a precursor to glycogen and lps_AFE, two biomass components. Their production is ensured by a pathway starting from β-f6p and β-g6p. They can be transformed into either α-g6p or in β-g1p. No other pathway can produce these targets. Therefore, β-g6p is a both a sustainability and a producibility-PEM, but it is not an optimal-efficiency PEM.

Among the six metabolic networks studied, this situation is not frequently observed in the iJO1366, A. ferrooxidans and T. lutea networks. In the Synecchocystis, iAF 1260 and iJR904 networks this situation occurs from eight to 22 times. In particular, the relatively high number of such PEMs in the Synecchocystis network suggests that its optimal functioning is not highly constrained. This could be explained by the presence of a relatively high number of export fluxes (this has already been observed when studying PEMs that are producibility-PEMs only) which generate redundancies with intracellular pathways.

Optimal-efficiency-PEMs that are neither sustainability nor producibility-PEMs shed light on optimal flux-distributions among alternative pathways.

These PEMs are required to produce optimal biomass. Removing all reactions for which this PEM is a substrate does not prevent the biomass from being stoichiometrically and topologically activated. For instance, in the Synecchocystis network, the removal of the reactions for which Cytosolic O-Phospho-L-serine is a substrate causes the maximal biomass growth rate to decrease from 47.5 to 41.1 mmol/gDW/h.

Our interpretation is that this mainly occurs when biomass components can be produced by at least two distinct pathways and when one of them is favored by the system on the basis of optimal flux-based analysis. An example is shown in Fig. 5. In the Synecchocystis network, the cytosolic putrescine is a precursor of the biomass components. In this network, two pathways produce this compound: one of them drives an import flux of extracellular putrescine and the other, an internal one, drives a pathway involving L-arginine and agmathine. CO2 is a co-product of the latter pathway. Its relationship with the production of coproduct causes the network to favor the L-arginine pathway for producing putrescine in order to optimize the biomass growth rate. Therefore, agmathine and putrescine are optimal-efficiency-PEMs. However, they are neither sustainability and producibility-PEMs because they are not required for the graph-based and the stoichiometry-based productions of putrescine.

In the six networks studied in this paper, optimal-efficiency-PEMs that are neither sustainability nor producibility-PEMs are the most frequent PEMs after PEMs that are sustainability, producibility and optimal-efficiency-PEMs simultaneously. In the iJR904, iAF1260 and iJO1366 networks, between 32% and 50% PEMs have this property. The most recent Synecchocystis and T. lutea networks includes around 10% of such PEMs. The exception is the A. ferrooxidans network with very few differences between the graph, stoichiometry and optimal flux-based formalisms.

Dependency of pathway activation on the initial state of the cell

Sustainability-PEMs that are neither producibility nor optimal-efficiency-PEMs can provide information about the dependency of pathways on initial cell state.

A PEM which satisfies only the topology-based criteria is a necessary compound for the topological activation of a biomass component. However, this compound is unnecessary for producing the corresponding biomass component from a stoichiometric viewpoint. For instance, according to graph-based criteria, in the iJO1366 network, the cytosolic thiamine is necessary for producing the thiamine diphosphate, a biomass component. Indeed, if we prune the metabolic network by removing all the reactions whose substrate is cytosolic thiamine, the thiamine diphosphate is no longer topologically activated. On the contrary, in the pruned graph, the growth rate maximal value (the biomass reaction) is unchanged compared to the growth rate of the initial network (11.747 mmol/gDW/h). This apparent paradox is explained mainly by the dynamical assumptions underlying the stoichiometric and topological activation semantics. We assert that these different semantics impact heavily on the interpretation of a cycle functioning with respect to the dependency of metabolites on their own production (Kruse & Ebenhöh, 2008).

An example is shown in Fig. 6. In the iJR904 network, the compound Peptidoglycan subunit, a biomass component, is produced from the metabolite uaagmda which is involved in a cycle through the Undecaprenyl diphosphate. At steady state, stoichiometry-based analysis suggests that the production of Peptidoglycan subunit is assured by the self-activation of the cycle, since all the reactions in the cycle are essential to the optimal production of the biomass. However, this type of analysis implicitly assumes that at least one of the component of the cycle is present on the activation of the system, and that none of the components is degraded during system functioning, guaranteeing the system producibility at steady state. On the contrary, in the graph-based framework, the analysis requires cycles to be initiated from external pathways in order to be considered activated. This corresponds in particular to situations where the cells are growing and require external input to ensure their sustainability, as pointed out in Kruse & Ebenhöh (2008). In the example shown in Fig. 6, the cycle is initiated by a linear pathway of ten reactions starting from the pyruvate metabolite. Our analysis shows that nine compounds are required this linear pathway, from a graph-based point of view, to produce the targeted compound Peptidoglycan subunit. On the contrary, they are not needed to produce Peptidoglycan subunit at steady state (stoichiometric-based accessibility), because this compound is produced by a self-induced cycle, thus the reaction flux is equal to zero for optimal flux-based analysis. This analysis suggests that the production of Peptidoglycan subunit is dependent on pyruvate, in addition to its link with a self-activated loop.

To summarize, one interpretation of sustainability-PEMs that are neither producibility nor optimal-efficiency-PEMs is that it suggests the existence of a self-activated cycle in the flux-based analysis (producibility at steady state) and provides candidates for the initiation of cycles when the cell is not at steady state (sustainability during growth phases) (Kruse & Ebenhöh, 2008). However, with the six metabolic networks studied, this phenomenon is observed very infrequently: once in the iAF1260, iJO1366 and T. lutea networks, nine times in the Synecchocystis network and eleven times in the iJR904 network. This suggests a need for curation for networks with a medium number of sustainability-PEMs that are neither producibility nor optimal-efficiency-PEMs.

Sustainability and optimal-efficiency-PEMs that are not producibility-PEMs provide insights on internal cycles for production of non-optimal biomass production

This is a tricky situation, for when a compound has this property, it is the substrate of a reaction that is always activated when the cell produces an optimal biomass growth rate. The sustainability property means that, if this reaction is removed from the network, the targeted biomass component is no longer activated according to a graph-based criteria. Nevertheless, the biomass component still has the capability to be activated at the stoichiometric level, although with a lower growth rate. An alternative means of production of the targeted component exists, although this alternative pathway is not connected to the set of seeds from a graph-based viewpoint. For instance, in the A. ferrooxidans network, removing the reactions which consume Alpha-D-Ribose_1-phosphate causes the optimal biomass growth rate to decrease from 3.6 to 3.4 mmol/gDW/h. Therefore it is still producible although there is no longer a pathway from the set of seeds to the biomass reaction. One interpretation of this paradox is similar to the case shown in Fig. 5, concerning the self-activated loop.

For instance, Fig. 7 shows a sub-network of the iJO1366 network. In the cytosol compartment, Thiamin, a PEM, is the only precursor of Thiamin monophosphate, a biomass component. Its production is assured by a pathway initiated by the set of seeds that imports Thiamin from periplasm to cytosol. Also, when the reaction from Thiamin to Thiamin monophosphate is removed from the network, the network has the ability to activate a cycle to produce 4-Methyl-5-(2-phosphoethyl)-thiazole (4mpetz). Notice, however, that this cycle is not connected to the set of seeds since its components are not topologically activated. Therefore, its activation is heavily dependent on the presence of at least one of its components at the initiation of the cell dynamics.

Regarding the six genome-scale networks studied, this situation mainly occurs in the iJO1366 network (six times) and the T. lutea network (seven times). Therefore, both networks have the capacity to adapt themselves by activating internal mass-balanced cycles.

Control of system response through co-products equilibria

Producibility-PEMs that are neither sustainability nor optimal-efficiency PEMs may point out to the fine tuning of mass-balance equilibria.

Such PEMs have the ability to set the biomass flux rate to zero when all their associated substrate reactions are removed from the network even though none of the removed fluxes is required to produce the biomass in the initial network. In addition, the pruning operation does not remove the graph-based activation of the biomass component. An example in the iAF1260 network is the cytosolic carbon dioxide. Removing this metabolite from the network, with all the reactions that use it as a substrate, does not affect the graph-based activation of any of the target metabolites. However the growth rate is reduced from 9,021 to 0 mmol/gDW/h. The biomass component cannot be produced because of the non-balanced mass equation despite the fact that it is theoretically still able to be activated topologically.

As shown in Fig. 8, these compounds are mainly the co-products of essential reactions that need to be metabolized to ensure the biomass production. In the periplasmic compartment of the iAF1260 network (Fig. 8), the compound murein5px4p (a biomass component) is a product of a reaction whose substrate is murein5p5p. Another product of this reaction is D-Alanine. When the export reactions of this co-product are removed from the network, the compound D-Alanine accumulates and can no longer satisfy the law of conservation of mass. Biomass production, in this pruned network, is impossible: flux-based analysis gives murein5px4p a production rate of zero. Nevertheless, removing reactions that export D-Alanine has no impact on the graph-based activation of the biomass components. In addition, this compound is not an optimal-efficiency PEM because there are two reactions that can export D-Alanine.

In other words, our interpretation is that producibility PEMs that are neither sustainability nor optimal-efficiency PEMs may refer to co-products of essential reactions. They may have a very large high impact on the organism growth rate since their consumption rate is a direct controller of biomass production at steady-state.

In practice, this situation is identified infrequently in all networks except the Synecchocystis network: once in the A. ferrooxidans and the T. lutea networks and two, three and four times in the iJR904, iAF1260 and iJO1366 networks, respectively. This observation may suggest that the co-products have been extensively studied in all of these networks. On the other hand, the Synecchocystis network contains 17 producibility PEMs that are neither sustainability nor optimal-efficiency PEMs. This suggests that biomass production is under the control of many co-product export fluxes which merit sensitivity analysis and reflect a need for curation.