Correlation-enhanced viable core in metabolic networks

Given our goal of characterizing a stable subset of metabolism by a structural criterion, we exploit the fact that a compound can maintain constant concentration only if at least one reaction produces it and another consumes it, which can support adjusting generation and consumption rates respectively, as introduced in Refs. Lemke et al. (2004); Smart et al. (2008). Translating this condition for our undirected networks, we define a compound $c$ as viable in a given network if its degree $k_{c}$, the number of reactions involving it and thus connected to it, is equal to or larger than two ($k_{c}\geq 2$); one reaction can generate and another can consume the compound under the assumption that the reactions are reversible. If a compound has no or just a single reaction as its neighbor, then it is considered as inviable. A reaction node $r$ is defined as viable if all of its $q_{r}$ connected compound nodes are viable; A reaction can occur only if none of its substrates and products are inviable. We present the examples of viable and inviable nodes in Fig. 1(a).

Then one can obtain the maximal induced-subgraph consisting of such viable compounds and viable reactions, which we define as the viable core. One can consider the viable core as the subset of metabolism which can maintain constant fluxes and concentrations of the reactions and compounds involved once they are established; Whether the viable core will reach such a steady state depends on the initial condition and the environment.

The viable core of a network can be obtained by removing iteratively all inviable nodes as follows. We first remove inviable compounds, i.e., the compound nodes of degree smaller than two. Their connected reaction nodes then become inviable and are removed. After this removal of inviable compounds and reactions in the first step, there can appear new inviable compounds, which had degree $k\geq 2$ but now have $k<2$ due to the loss of some of their connected reactions in the previous step. These new inviable compound nodes and their connected reaction nodes are removed in the second step. Such pruning is repeated until there are no more inviable nodes left, when the remaining compound and reaction nodes constitute the viable core. An example is presented in Fig. 1(a). This pruning process is almost the same as the greedy leaf removal (GLR) process considered in Ref. Liu et al. (2012), in which the nodes of degree smaller than two and all of their neighbor nodes are removed iteratively in unipartite networks without distinguishing the type of nodes.

Applying the procedure described in the previous subsection, we have obtained the viable core of the metabolic network for each species. As an example, the viable core of Klebsiella pneumoniae is shown in Fig. 1(b), which includes 18 viable reactions and 29 viable compounds among a total of 62 reactions and 91 compounds. To understand how large or small the viable core of a species $s$ is, we consider its relative size $v^{s}$, the ratio of the number of nodes in the viable core to the total number of nodes. We compare it with the average relative size $\langle v^{s}\rangle_{\rm ran}$ of the viable cores of the randomized networks, obtained by rewiring links while preserving the degrees of individual nodes; Each instance of these randomized networks, shown in Fig. 1(c) as an example, has been created by swapping $L^{s}$ times the end nodes of two randomly chosen links ($L^{s}$ is the total number of links for species $s$) and we have generated one hundred such instances for each species. One can see that the viable core of the real network of Klebsiella pneumoniae [Fig. 1(b)] is much larger than that of a randomized network [Fig. 1(c)]. This holds for all species; The empirical core size ranges between 0.1 and 0.5 across species, and $v^{s}>\langle v^{s}\rangle_{\rm ran}$ (P-value less than $6\times 10^{-6}$) for every species $s$ as shown in Fig. 1(d). This result means that the core is larger than expected based on the given degree sequences when others are random, and suggests that the network structure beyond the degree sequence plays a role in forming the viable core. Such enhanced stability has been reported in different contexts for a few selected species, e.g., it has been shown that the cascades of failures are more constrained in the real metabolic networks than in the randomized networks Smart et al. (2008); Güell et al. (2014). Also, as we will show below, highly heterogeneous degrees of compounds Jeong et al. (2000) strongly influence the size of the viable core.

We define the viability of each reaction $v_{r}$ as the fraction of the species having the reaction in their viable cores. It can quantify the likelihood that a reaction belongs to the viable core and is thus stable. In Fig. 1(e), it is shown that a significant number of reactions ($5549/12077\simeq 46\%$) do not belong to the viable core in any species, while their viability in the randomized networks is mostly non-zero and distributed broadly—5530 reactions have $v_{r}=0$ but $\langle v_{r}\rangle_{\rm ran}>0$. On the other hand, the other 6528 reactions have non-zero viability in the empirical networks and tend to have lower viability in the degree-preserving randomized networks, i.e., $v_{r}>\langle v_{r}\rangle_{\rm ran}$. These findings indicate the stronger heterogeneity of reaction’s viability in empirical networks than in randomized networks and imply that the viable cores of real metabolic networks are primarily formed by selected high-viability reactions.

We have shown that the real-world metabolic networks have larger stable components than expected for the randomized networks. Will it hold also even with perturbations? In terms of the viable core, we can seek the answer by measuring the viable core in the damaged networks.

To be specific, we remove a fraction $Q$ of randomly chosen reaction nodes to mimic their malfunction due to e.g., the inactivation or mutation of the enzymes catalyzing them, and then compute the viable core. For Escherichia coli [Fig. 2(a)], the real damaged networks have larger cores than the damaged degree-preserving randomized networks in the whole range of $Q$. For all species $s$, we measure $v^{s}(Q)$ and $\langle v^{s}\rangle_{\rm ran}(Q)$ as functions of $Q$ in the original and degree-preserving randomized networks, respectively. We find that the characteristic core size $\widetilde{v^{s}}\equiv\int_{0}^{1}v^{s}(Q)dQ$, an average of $v$ in the range $0\leq Q\leq 1$, of the real damaged networks is much larger than that of the damaged randomized networks $\widetilde{\langle v^{s}\rangle}_{\rm ran}$ as shown by comparing their distributions over species [Fig. 2(b)].

Since the degree-preserving randomized networks are obtained by shuffling connection among nodes while preserving nodes’ degrees, they have no or very weak correlation between the degrees of connected nodes Park and Newman (2003); Catanzaro et al. (2005). For such uncorrelated networks, the size of the viable core can be obtained analytically by solving self-consistent equations formulated in terms of the degree distributions as done in Refs. Dorogovtsev et al. (2006); Liu et al. (2012); Lee et al. (2023) or the adjacency matrix, called the message-passing method, as done recently in Ref. Bianconi and Dorogovtsev . Here we take the former approach to understand the influence of the degree distribution on the core size in uncorrelated networks.

Following Ref. Liu et al. (2012), let us first consider a compound node reached by following a link from a viable reaction node and denote the probability of the compound node to be inviable by $\alpha$. The compound node will be inviable if all the remaining neighboring reaction nodes, except the viable one at the end of the followed link, are inviable. Next, let us consider a reaction node reached by following a link from a viable compound node, and denote the probability of the reaction node to be inviable by $\beta$. The reaction node will be viable if (i) it is not removed by the random removal of $Q$ fraction of reaction nodes and (ii) all of its neighboring compound nodes are viable. These reasonings lead us to $\alpha=\sum_{k}\frac{kP_{\rm cmp}(k)}{\bar{k}}\beta^{k-1}$ and $1-\beta=(1-Q)\sum_{q}\frac{qP_{\rm rxn}(q)}{\bar{q}}(1-\alpha)^{q-1}$, where $P_{\rm cmp}(k)\,[P_{\rm rxn}(q)]$ and $\bar{k}\,(\bar{q})$ are the degree distribution and the mean degree of compound (reaction) nodes, respectively. Introducing the generating functions $G_{\rm cmp}(z)\equiv\sum_{k=0}^{\infty}P_{\rm cmp}(k)z^{k}$ and $G_{\rm rxn}(z)\equiv\sum_{q=0}^{\infty}P_{\rm rxn}(q)z^{q}$, one can represent the previous equations as

Recalling that a reaction node is viable if it is not removed and all of its neighboring compound nodes are viable and a compound node is viable if it has at least two neighbor reaction nodes that are viable, one can find the probability $v_{\rm rxn}\,(v_{\rm cmp})$ of a reaction (compound) node to be viable as

Finally, the relative size $v$ of the viable core is given by

with $n_{\rm cmp}\,(n_{\rm rxn})$ being the fraction of compound (reaction) nodes, satisfying $n_{\rm cmp}+n_{\rm rxn}=1$, in the considered metabolic network.

The solutions to Eqs. (1), (2), and (3) are dependent only on the degree distributions $P_{\rm cmp}(k)$ and $P_{\rm rxn}(q)$. The degree distributions averaged over species are shown in Figs. 2(c) and 2(d). The broad degree distributions for compounds are universal, which leads the viable cores to shrink gradually as $Q$ increases as shown in Fig. 2(a) Liu et al. (2012); Dorogovtsev et al. (2006). Such a gradual shrink may not be the case when the degree distribution is narrow and the mean degrees are sufficiently large. In the completely randomized networks that are obtained by assigning links to randomly selected node pairs and thus preserve the total number of links but do not preserve the degree sequences, the viable core size drops abruptly to a small value at a critical value $Q_{c}$ [Fig. 2(a)], denoted by ‘link-preserving’. With the degree distributions given in the Poissonian form, one can find that such discontinuous transitions can occur if the mean degrees are larger than $e=2.718182\ldots$, i.e., $\bar{k}>e$ and $\bar{q}>e$ Liu et al. (2012); Dorogovtsev et al. (2006). See Appendix A for derivation. We find 5463/5469 $\simeq$ 99.8% species satisfy this condition.

Our analysis shows that the degrees $q$ and $k$ of the end nodes, reaction- and compound-type respectively, of a link are positively correlated in the empirical networks. Their Pearson correlation coefficient or assortativity is defined for each species $s$ as  Newman (2002),

where $A_{cr}^{s}$ is the adjacency matrix of the metabolic network, and $k_{\rm nn}^{s}\equiv\sum_{c,r}(A_{cr}^{s}/L^{s})k_{c}^{s}$ and $q_{\rm nn}^{s}\equiv\sum_{c,r}(A_{cr}^{s}/L^{s})q_{r}^{s}$ are the mean degrees of the end nodes of a link, and $\sigma_{\rm cmp}^{s}\equiv[\sum_{c,r}(A_{cr}^{s}/L^{s})(k_{c}^{s}-k_{\rm nn}^{s})^{2}]^{1/2}$ and $\sigma_{\rm rxn}^{s}\equiv[\sum_{c,r}(A_{cr}^{s}/L^{s})(q_{r}^{s}-q_{\rm nn}^{s})^{2}]^{1/2}$ are the standard deviations. The assortativity of the real metabolic networks is given by $\rho_{\rm emp}\simeq 0.13\pm 0.02$ and significant with P values less than 0.01 for almost all species. Also we compute the probability density function (pdf) of the standardized degrees $\tilde{k}_{c}^{s}={k_{c}^{s}-k_{\rm nn}^{s}\over\sigma_{\rm cmp}^{s}}$ and $\tilde{q}_{r}^{s}={q_{r}^{s}-q_{\rm nn}^{s}\over\sigma_{\rm rxn}^{s}}$, and the conditional mean $\langle\tilde{q}^{s}\rangle(\tilde{k})=\sum_{c,r}(A_{cr}^{s}/L^{s})\tilde{q}_{r}^{s}\delta(\tilde{k}_{c}^{s}-\tilde{k})/\sum_{c,r}(A_{cr}^{s}/L^{s})\delta(\tilde{k}_{c}^{s}-\tilde{k})$. Note that the assortativity is one of the moments of the pdf, $\rho=\int d\tilde{k}d\tilde{q}P(\tilde{k},\tilde{q})\tilde{k}\tilde{q}$. In Fig. 3(a) are shown the pdf and the conditional mean averaged over species. One can see that $\langle\tilde{q}\rangle(\tilde{k})$ grows with $\tilde{k}$ in the most probable range $-1\lesssim\tilde{k}\lesssim 0$ supporting the positive correlation between $\tilde{q}$ and $\tilde{k}$.

To scrutinize the influence of the degree-degree correlation on the size of the viable core, we generate an ensemble of artificial degree-preserving correlated networks that exhibit a given assortativity. To be specific, for given species $s$, we perform the Monte-Carlo sampling of different configurations—different adjacency matrices $A_{cr}$—while preserving the given degree sequence by swapping the end nodes of two links with a Hamiltonian $\mathcal{H}=-J\sum_{c,r}A_{cr}k_{c}^{s}q_{r}^{s}$. The assortativity is measured in the equilibrated networks for each $J$. Note that the degree-preserving randomized networks are generated with $J=0$. The relative size $v$ of the viable cores of these degree-preserving correlated networks, averaged over species, increases as assortativity $\rho$ increases over a wide range of $\rho$, i.e., $-0.5\lesssim\rho\lesssim 0.2$, which includes most of the empirical values of assortativity [Fig. 3(b)]. This demonstrates the contribution of the positive degree-degree correlations to the enhancement of the viable cores in the empirical networks. Yet, we also note that they are not the only cause of the core enhancement; The size of the empirical viable cores, $v=0.4\pm 0.04$, is larger than $v=0.22\pm 0.02$ of those artificial correlated networks at the empirical values of assortativity $\rho_{\rm emp}$.

To understand why a positive (negative) correlation makes a larger (smaller) viable core for a given degree sequence, we investigate the branching ratio $g$ of inviable compounds between adjacent steps of pruning; If there are $C_{1}$ inviable compounds in the current step and after removing them and their connected reactions nodes, there appear $C_{2}$ newly inviable compounds in the next step, the branching ratio will be given by $g\equiv{C_{2}\over C_{1}}$. One can expect a larger (smaller) core with a smaller (larger) branching ratio. To see the dependence of the branching ratio on assortativity, we decompose the branching ratio into the compound-to-reaction one and the reaction-to-compound one as

where $R_{1}$ is the number of the current-step inviable reaction nodes, those connected to the $C_{1}$ current-step inviable nodes, and $C_{2}$ is the number of the next-step inviable nodes that will be made inviable only in the next step by pruning in the current step.

If a network has a negative degree-degree correlation, the compound-to-reaction branching ratio $g_{{\rm cmp}\to{\rm rxn}}$ can be small; Inviable compounds may be concentrated, being connected to a few hub reaction nodes. Each of those current-step inviable reaction nodes has many neighbor compounds and thus can leave many next-step inviable compounds, resulting in large $g_{{\rm rxn}\to{\rm cmp}}$. In contrast, $g_{{\rm cmp}\to{\rm rxn}}$ cannot be so small in a network with a positive degree-degree correlation since inviable compounds are distributed, each being connected to each different reaction node of a small degree. Those current-step inviable reaction nodes have small degrees and thus each will leave a small number of next-step inviable compounds, resulting in small $g_{{\rm rxn}\to{\rm cmp}}$.

These reasonings can be summarized as follows: As assortativity increases, $g_{{\rm cmp}\to{\rm rxn}}$ increases but $g_{{\rm rxn}\to{\rm cmp}}$ decreases. This is confirmed in our artificial correlated networks. In the inset of Fig. 3(c), we show the two branching ratios obtained by computing $C_{1},R_{1}$, and $C_{2}$ in the early-stage of pruning in each correlated network and averaging over species, which behave as functions of $\rho$ as expected within a range of interest. A crucial observation is that the increase of $g_{{\rm cmp}\to{\rm rxn}}$ with $\rho$ is much weaker than the decrease of $g_{{\rm rxn}\to{\rm cmp}}$ with $\rho$. Such quantitative difference in the dependence on $\rho$ results in the decrease of $g$, the product of the two branching ratios, with increasing $\rho$; $g\simeq 0.4$ at $\rho=-0.4$ while $g\simeq 0.35$ at $\rho=0.1$.

Such dependencies of the branching ratios and the core size on assortativity help us understand how a structural correlation can influence the viable core. We should remark that the behaviors of $v(\rho)$ and $g(\rho)$ may be changed in extremely correlated networks; $v$ decreases as $\rho$ increases in the range $\rho\lesssim-0.5$ and $\rho\gtrsim 0.2$ [Fig. 3(b)] and $g$ increases with $\rho$ in the range $\rho\gtrsim 0.2$ [Fig. 3(c)]. The investigation of these deviations may be interesting and need further investigation. We also note that the empirical branching ratio $g=0.2\pm 0.03$ is smaller than the ratio $g\simeq 0.34$ from the artificial correlated networks at the empirical value of assortativity, implying again the influence of higher-order structure on the viable core and the branching ratio of pruning.

Given that a positive degree-degree correlation can enlarge the viable core and thereby enhance stability, one can suspect that an evolutionary pressure might have been imposed towards generating the positive correlations in the real metabolic networks during evolution. While it is a simple and compelling scenario, here we present evidence that a neutral evolution of metabolism, without an evolutionary pressure, over a growing tree of species on a long-time scale can generate the positive correlations.

The frequency $f_{r}$ ($f_{c}$) of a reaction $r$ (compound $c$) to be present in a species’ metabolism among species has been found to be so heterogeneous  Bernhardsson et al. (2011); Kim et al. (2015, 2019) that its distribution takes a power-law form with exponent one Kim et al. (2015, 2019); Lee and Lee (2024). It has been shown in Ref. Lee and Lee (2024) that such heterogeneous frequencies of reactions and compounds may originate from their different times of the first appearance in the metabolism of a living species and their inheritance to the descendant species during the co-evolution of metabolism and species tree. In this framework, the early-recruited compounds and reactions have high frequency and moreover, are expected to have larger degrees in the metabolic networks of the species containing them since they have had more chances to be connected to other nodes during evolution than the late-recruited ones. Therefore we can think of the empirical positive degree-degree correlations as originating from the co-recruitment of connected reactions and compounds in the evolving metabolism in a growing species tree.

The empirical data support this expectation. The Pearson correlation coefficient between the frequency $f_{c}$ and the degree $k_{c}^{s}$ of a compound at the end of a link in a species is $0.34\pm 0.03$ (P $\lesssim 4\times 10^{-6}$ for all species) with the fluctuation arising among species, demonstrating their significant correlations. In each species, the end compound node of a link is likely to have high frequency and such high-frequency compounds show clearly a positive correlation between the standardized degree $\tilde{k}_{c}$ and the frequency $f_{c}$ as shown by their joint pdf and the conditional mean $\langle\tilde{k}_{c}\rangle(f_{c})$ [Fig. 4(a)]. The correlation between $f_{r}$ and $q_{r}^{s}$ of a reaction is $0.06\pm 0.06$, which is weaker than the correlation for compounds’ degree and frequency but still significant for many species (P $\lesssim 0.05$ for $4130/5469\approx 76\%$ species). The pdf $P(f_{r},\tilde{q}_{r})$ and the conditional mean $\langle q_{r}\rangle(f_{r})$ show that the correlation is particularly significant for reactions of high frequency in Fig. 4(b).

Also, as expected, the frequencies of a connected pair of reaction and compound nodes are strongly correlated; Their Pearson correlation coefficient is $0.42\pm 0.04$ (P value $\lesssim 0.0015$ for all species) and the pdf $P(f_{c},f_{r})$, and the conditional mean $\langle f_{r}\rangle(f_{c})$ clearly show the correlation in Fig. 4(c). Collecting all these results, we find that for a connected pair of a compound node and a reaction node, their degrees are positively correlated with their frequencies respectively and their frequencies are very similar, as expected from the metabolism evolution model as sketched above. All these factors result in the positive correlations of the degrees of the pair and, in turn, contribute to the enhancement of the viable core.