Skip to main content


Springer Nature is making Coronavirus research free. View research | View latest news | Sign up for updates

Modular co-evolution of metabolic networks



The architecture of biological networks has been reported to exhibit high level of modularity, and to some extent, topological modules of networks overlap with known functional modules. However, how the modular topology of the molecular network affects the evolution of its member proteins remains unclear.


In this work, the functional and evolutionary modularity of Homo sapiens (H. sapiens) metabolic network were investigated from a topological point of view. Network decomposition shows that the metabolic network is organized in a highly modular core-periphery way, in which the core modules are tightly linked together and perform basic metabolism functions, whereas the periphery modules only interact with few modules and accomplish relatively independent and specialized functions. Moreover, over half of the modules exhibit co-evolutionary feature and belong to specific evolutionary ages. Peripheral modules tend to evolve more cohesively and faster than core modules do.


The correlation between functional, evolutionary and topological modularity suggests that the evolutionary history and functional requirements of metabolic systems have been imprinted in the architecture of metabolic networks. Such systems level analysis could demonstrate how the evolution of genes may be placed in a genome-scale network context, giving a novel perspective on molecular evolution.


Cellular functions are carried out in a modular way, and functional modules are basic building blocks of cellular organization [1]. From the perspective of molecular biology, a functional module is regarded as a group of spatially isolated or chemically specific biological components that work together for a discrete biological function. Various functional modules such as protein complexes [24], signalling/metabolic pathways [58] and transcriptional clusters [9, 10] have been detected from functional genomic techniques or bioinformatics analyses of genomic data. Recent studies suggest, to varying degrees, functional modules correlate with evolutionary modules [11], the latter being defined as cohesive evolutionary blocks in cellular systems [12, 13]. It was found that genes within functional modules tend to evolve in a coordinated way [1215], while some fraction of evolutionary modules (or phylogenetic modules) agree well with known functional modules [1619].

On the other hand, purely topological analysis by graph-theoretic methods has revealed that molecular networks, such as protein interaction [2022], gene regulatory [23, 24] and metabolic networks [2530], consist of topological modules – densely connected sub-networks within which there is a high density of edges, and between which there is a lower density of edges [31]. Since graph-theoretic methods analyze networks from topological point of view using minimal prior knowledge about biological function or evolution, they have the potential to shed new light on biological systems based on the unbiased structural information [32]. Actually, numerous studies have demonstrated, to some extent, topological modules in molecular networks tend to be functionally modularized [2030]. Furthermore, studying molecular evolution from the viewpoint of network architecture is becoming a subject of current interest. Some recent studies suggested that the node degrees of molecular networks may constrain the evolution of proteins [3338], and protein interaction hubs situated within modules are more evolutionarily constrained than those bridging different modules [39, 40]. However, little has been known about how network modularity affects protein evolution. Thus more studies are expected to reveal the possible correlation between topological modules and evolutionary modules in molecular networks.

In this study, we ask to which extent the identified topological modules of metabolic networks co-evolve. We explore this question by analysing the metabolic network of H. sapiens (hsa) reconstructed from the KEGG database [4143]. We first break up the metabolic network into modules by the simulated annealing algorithm proposed in [28] and study the linkage pattern between modules. Then we investigate the evidence for co-evolution of modules by analysing the phylogenetic profiles, evolutionary ages and evolutionary rates of enzyme genes within modules. To mine the inherent relations between structure, function and evolution of metabolic networks, the features from H. sapiens network were then compared with those of the properly randomized counterparts, here, topological null model and biological null model, respectively.

Results and discussion

Identifying topological modules and their functions

We reconstructed the metabolic network of H. sapiens from the KEGG database [4143] and represented the network by a directed substrate graph in such a way that the nodes correspond to metabolites and arcs correspond to enzyme-catalyzed reactions between these metabolites [44]. The metabolic network of H. sapiens consists of 1378 metabolites and 666 enzymes, with the biggest connected cluster includes 948 metabolites and 614 enzymes.

The simulated annealing algorithm [28] was utilized to decompose the metabolic network. Totally 25 topologically compact modules were obtained. In Figure 1 we display the decomposition result as a network of modules, in which each node corresponds to a module and is represented by the functional cartography [28].

Figure 1

Cartographic representation of the metabolic network for H. sapiens. Each circle represents a module and is coloured according to the KEGG pathway classification of the reactions belonging to it, while the arcs reflect the connection between clusters. The area of each colour in one circle is proportional to the number of reactions that belong to the corresponding metabolism. The width of an arc is proportional to the number of reactions between the two corresponding modules. For simplicity, bi-directed arcs are presented by grey edges. The modularity metric of this decomposition is 0.8868.

Figure 1 exhibits a global view of interactions between modules, suggesting that the modules are linked in a core-periphery organized pattern [45, 46]. Some modules interact frequently and are interconnected densely to form a core, while others such as module 3, 6, 7 and 14 communicate with only one or two other modules and reside in the periphery of the network. We define the inter-module degree of one module as the number of its links with other modules, where a link by a bi-directed arc is counted as degree 2. Hence core modules have high inter-module degrees, while periphery modules have low inter-module degrees.

Further investigation revealed that most periphery modules correspond to a well-defined single pathway and perform relatively independent function. On the contrary, the core modules are mixtures of several conventional biochemical pathways, thus are difficult to be assigned a simple function. Table 1 shows that the periphery modules carry out two categories of functions: catabolism of essential human nutrients taken from the diet (essential lipids, amino acids and vitamins), biosynthesis and metabolism for complex molecules (hormones, glycans and cofactors), whereas the core modules generally perform basic metabolisms of sugars, amino acids, lipids and nucleotides. Since the H. sapiens metabolic network does not include biosynthetic pathways for essential human nutrients, the pathways for processing these substrates may have functioned in characterizing the species during evolutionary diversification. On the other hand, owing to the function of hormones and glycans in mediating cell recognition, communication, and signal transduction, the ability to produce and utilize these complex molecules may have led the species to develop an advanced neural system for physiological regulation along the evolutionary process [47]. Thus the functions performed by periphery modules could be regarded as some high-level metabolisms specialized for H. sapiens, whereas the core modules are integrated to accomplish housekeeping processes of life. Especially, the three central pathways – Embden-Meyerhof-Parnas (EMP), tricarboxylic acid (TCA) and pentose phosphate pathway (PPP) are broken up into several core modules, mainly in module 8, 11 and 13. One possible explanation is that the metabolites in these central pathways are used as common precursors for biosynthesis of universal building blocks [48] and are thus placed in different modules. In addition, these metabolites have two or more functional classifications and that they may get assigned to different modules based on features other than the central pathway they belong to. This result agrees with earlier observations concerning the high diversity of the TCA and EMP pathway [49, 50], as well as the clustering results for Escherichia coli (E. coli) metabolic network obtained by other algorithms [19, 25, 26, 30].

Table 1 Inter-module degrees (d), number of nodes (N), and main functions of the topological modules for H. sapiens network

The similarity between the phylogenetic profiles of enzymes within modules

The phylogenetic profiles of any pair of enzymes can be compared to determine whether these two enzymes exhibit significant co-occurrence, so as to result in similar phylogenetic profiles. We adopted the Jaccard coefficient (JC) to measure the similarity between the phylogenetic profiles of enzyme pairs. Figure 2(A) shows the distribution of JC for all pairwise enzymes in the H. sapiens metabolic network. It can be seen that most enzyme pairs exhibit low extent of similarity on their phylogenetic profiles. The average JC of the global network is 0.28. The significant level P (0.05) of the distribution is 0.66, i.e., the probability of enzyme pairs with JC bigger than 0.66 is only 0.05. Thus we set JC of 0.66 as a threshold of similarity between phylogenetic profiles and regard enzyme pairs whose JC are bigger than 0.66 as having similar phylogenetic profiles. Figure 2(B) illustrates the relationship between the average JC for enzyme pairs within modules and the inter-module degree of module (Spearman's rank correlation is r = -0.3814, P-value is 0.059), suggesting a statistically negative correlation tendency between these two variables. Since modules with high inter-module degrees are core modules, and those with low inter-module degrees are periphery modules, this result means that enzymes within periphery modules have higher extent of similar phylogenetic profiles than those within core modules.

Figure 2

The distribution of Jaccard coefficient (JC) for enzyme pairs in the global metabolic network and its modules. (A) The distribution of JC for all pairwise enzymes in the H. sapiens metabolic network. (B) The relationship between the average JC for enzyme pairs in modules and the inter-module degree of module.

We compare the similarity of phylogenetic profiles for enzyme pairs within each module with that within the global network. Figure 3(A) shows that enzyme pairs within most modules have higher average JC than those within the global network, and Figure 3(B) demonstrates that most modules include high percentage of enzymes with similar phylogenetic profiles. We used hypergeometric cumulative distribution [51] to quantitatively measure whether a module is more enriched with enzymes of similar phylogenetic profiles (JC ≥ 0.66) than that by chance. Given significance level α = 0.05, a P-value smaller than α demonstrates low probability that the enzymes of similar phylogenetic profiles appear in the same module by chance.

Figure 3

Comparison of the similar extent of phylogenetic profiles for enzyme pairs within each module with that within the global metabolic network of H. sapiens. (A) Average JC of enzyme pairs within modules. The red column represents the global network. The modules are ordered according to their average JC in a decreasing way. (B) Percentage of enzyme pairs within modules with JC ≥ 0.66 (threshold definition). The red column represents the global network. The modules are drawn in the same order as in (A).

Integrating these three measures, we regard a module as an evolutionary module enriched with enzymes of similar phylogenetic profiles, if it satisfies all of the following three criteria:

  1. (1)

    Average JC of the module is bigger than that of the global network.

  2. (2)

    The fraction of enzyme pairs with JC ≥ 0.66 (definition of a threshold) in the module is significantly bigger than 0.05, i.e., the fraction of that in the global network. We set the cutoff to 0.1.

  3. (3)

    The P-value is smaller than α = 0.05.

As shown in Figure 3, all modules except module 2 satisfy the criteria (1), in which 13 modules (module 7, 3, 25, 9, 16, 4, 6, 22, 12, 15, 19, 21 and 20) also satisfy the criteria (2). Of the 13 modules, only the P-value of module 20, which equals to 0.50067, is bigger than 0.05. That is to say, although the average JC of enzyme pairs in module 20 is big enough, and this module also includes a high fraction of enzymes with similar phylogenetic profiles, this case has a high probability to occur by chance. Thus module 20 could not be regarded as an evolutionary module by our criteria.

In summary, a total of 12 modules out of 25 (module 7, 3, 25, 9, 16, 4, 6, 22, 12, 15, 19 and 21) were found to be evolutionary modules, most of which are periphery modules. The inter-modules degree of eight modules are less than 5, suggesting that periphery modules behave more cohesively in evolution than core modules. That is to say, enzyme genes within periphery modules have higher tendency to be gained/lost together than those within core modules.

The evolutionary ages of modules

We classified enzyme genes in the H. sapiens network into seven evolutionary ages: Prokaryota, Protists, Fungi, Nematodes, Arthropods, Mammalian and Human. Hypergeometric cumulative distribution [51] was used to measure whether a module is more enriched with enzymes from a particular evolutionary age than would be expected by chance. Given significance level α = 0.05, a P-value smaller than α demonstrates low probability that the enzyme genes of a particular evolution age have appeared by chance.

We defined the evolutionary age of a module as the biggest value of the evolutionary age of enzymes included in this module, which satisfies all of the two criteria:

  1. (1)

    More than 1/3 enzymes of this module belong to the evolutionary age;

  2. (2)

    The corresponding P-value is smaller than α = 0.05.

According to this definition, a total of 16 modules appear to belong to one specific evolutionary age. See Table 2 for the evolutionary age, percentage of enzymes in the specific age, P-value and inter-module degree for each module. As can be seen in Table 2, all of the age-1 modules except for module 12 are core modules, while the modules with later evolutionary ages are periphery modules. Figure 4 shows the statistically significant negative correlation between evolutionary age of modules and average inter-module degrees. The Spearman's rank correlation between evolutionary age and inter-module degree is r = -0.8207 (P-value = 9.78 × 10-5).

Table 2 Evolutionary ages of topological modules for H. sapiens network
Figure 4

Relationship between evolutionary age of modules and average inter-module degrees.

The distribution feature of core and periphery modules in different evolutionary ages provides some evidence for the evolutionary history of metabolic networks. Matching the evolutionary ages of modules with their functions suggests how ancient cellular functions may have evolved to gain new phenotypes with improved adaptation: the core modules appeared earlier in evolution and communicate frequently to perform the basic functions, while the periphery modules "sprout" from the compactly inter-connected core modules later via sparse linkages to carry out some novel, specialized functions. In this way, the housekeeping functions are conserved in core modules while functional specialization is achieved by extending periphery modules.

Evolutionary rates of constituent enzyme genes of modules

We adopted the evolutionary rate to estimate the evolutionary constraints on metabolic enzymes. A small value of evolutionary rate suggests a smaller fraction of accepted amino acid substitutions (see Methods part), hence a higher evolutionary constraint on the enzyme. We extracted the evolutionary rate of each enzyme gene included in the H. sapiens metabolic network from the HomoloGene database [52], computed by the approaches in [53]. Then the evolutionary rates were averaged over all enzyme genes within the same module.

The relationship between average evolutionary rate of enzyme genes within modules and inter-module degrees is shown in Figure 5. Figure 5(A) indicates a statistically significant negative correlation between the average evolutionary rate of module and inter-module degree (Spearman's rank correlation is r = -0.4983, P-value = 0.011.). The histogram of average evolutionary rate of modules versus binned inter-module degree is displayed in Figure 5(B), indicating that enzymes in core modules evolve more slowly than those in periphery modules. That is to say, core modules have higher evolutionary constraints, thus are more evolutionarily conserved than periphery modules. This observation is in agreement with the housekeeping functions of core modules.

Figure 5

Relationship between average evolutionary rate of enzyme genes within modules and inter-module degrees. (A) The average evolutionary rate of module is negative correlation with the inter-module degree. (B) Core modules (inter-module degree > 5) are evolutionarily constrained at higher extent than periphery modules (inter-module degree ≤ 5) do.

Comparison of the H. sapiens network with its randomised counterparts

In order to investigate whether the topological, functional and evolutionary features of modularity we presented above are intrinsic for metabolic networks, we constructed two versions of randomised counterparts for the H. sapiens metabolic network and compared them with the real network.

The first version of randomized network, called topological null model, was constructed by rewiring the links between nodes while preserving some low-level topological properties of the H. sapiens metabolic network. When shuffling the links, considering that the distribution of bi-directed arcs is an inherent topological feature of metabolic networks [54], we kept not only the degree of each node [55, 56], but also the total number of directed and bi-directed arcs of the metabolic network (see the algorithm in [30]). Totally 50 random networks were generated and then decomposed by simulated annealing, resulting in mean modularity metric 0.7804 and the standard deviation 0.0056. The modularity metric of the H. sapiens metabolic network is 0.8868, 19 multiples of standard deviation above that of the randomized counterparts, suggesting that the higher modularity of H. sapiens network is unlikely to arise by chance. Moreover, as shown in Figure 6 and Table A2 of additional file [see Additional file 1], the modules of the random network are so densely connected that their linkage pattern does not show a clear-cut core-periphery dichotomy. These comparisons indicate that the highly modular core-periphery organization of H. sapiens network could be an intrinsic topological character of metabolic networks, rather than a random phenomenon.

Figure 6

Decomposition of one randomized network of the H. sapiens network, shown as a network of modules. Each circle represents a module, while the arcs reflect the connection between clusters. The width of an arc is proportional to the number of links between the two corresponding modules. For simplicity, bi-directed arcs are presented by grey edges. The modularity metric of this decomposition is 0.7845.

We generated the second version of randomized network, called biological null model, through shuffling the enzymes that catalyze the reactions while preserving the network topology. At each step of randomization, we randomly chose two reactions that are catalyzed by different enzymes and then exchanged the enzymes catalyzing them. For one randomized network, we performed the same analysis about the functional distributions, phylogenic profiles, evolutionary ages, and evolutionary rates of its topological modules as we did for the H. sapiens network. Although the randomized network has the same topology as the H. sapiens network, its topological modules are high heterogeneity of reactions from different pathways [see Figure A1 of Additional file 1], thus could not be functional specific modules. From evolutionary view, enzymes within the same modules do not exhibit similar phylogenic profiles [see Figure A2 of Additional file 1]; neither the modules have specific evolutionary ages. In addition, the correlation between the average evolutionary rates of modules and the inter-module degrees is not statistically significant (Spearman's rank correlation between evolutionary rate and inter-module degrees was calculated as r = -0.0553, P-value = 0.793). These results demonstrate that the topological modules of the biological null model show absence of evolutionary modularity. Therefore, the functional and evolutionary modularity of topological modules could be inherent for metabolic networks.


In this study we have conducted a system-level survey about how evolution of enzyme genes is related to the structure and function of metabolic networks. From topological point of view, metabolic networks exhibit highly modular core-periphery organization pattern. Furthermore, the core modules are more evolutionarily conserved and perform some housekeeping metabolism functions, while the periphery modules appear later in evolution history and accomplish relatively specific functions. Our results suggest that the core-periphery modularity organization reflects the functional and evolutionary requirements of metabolic systems. The denser inter-connections between core modules may offer effectual protections to the basic metabolic process and keep the robustness of the metabolic system. On the other hand, the looser inter-linkages of periphery modules could function in favour of the flexibility and evolutionary ability of the system, so that the mutation or evolution of these parts may generate new phenotypes with improved adaptation while not significantly affecting other modules or even causing malfunction of the whole system. Our observation may shed light on a more global understanding of the topology, function and evolution for metabolic networks.


Data preparation and network reconstruction

In this study, the metabolic network of H. sapiens (hsa) was reconstructed using data downloaded (in May 2006) from the FTP service of KEGG (Kyoto Encyclopedia of Genes and Genomes). The "hsa_enzyme.list" file in the GENOME section of the KEGG database includes a list of the known enzymes encoded by H. sapiens's genome and the corresponding genes. The "reaction" file in the LIGAND section was first scanned for all reactions catalyzed by enzymes present in H. sapiens's genome, and totally 1492 reactions were determined. Then, the resulting reactions were matched to the "reaction_mapformula.lst" file, which includes direction information and the main metabolites for each reaction.

Some small molecules, such as adenosine triphosphate (ATP), adenosine diphosphate (ADP), nicotinamide adenine dinucleotide (NAD) and CO2, are normally used as carriers for transferring electrons or certain functional groups and participate in many reactions, while typically not participating in product formation. Therefore, in order to reflect biologically relevant transformations of substrates, we excluded these kinds of small molecules whose list is as follows [57, 58]:

ATP, ADP, AMP, NAD, NADH, NADP, NADPH, NH3, CoA, O2, CO2, Glu, Pyrophosphate, H+.

To construct the metabolic network, substrates and products were extracted from each of the enzyme-catalyzed reactions and all resulting substrate-product pairs were listed to specify the connections between the substances. A metabolic network is represented by a directed graph whose nodes correspond to metabolites and whose arcs correspond to reactions between these metabolites, in which irreversible reactions are presented as directed arcs while reversible ones as bi-directed arcs. For example, the irreversible reaction,

(S)-2-Methylmalate → Acetate + Pyruvate

corresponds to two directed arcs, i.e., (S)-2-Methylmalate → Acetate and (S)-2-Methylmalate → Pyruvate. The resulting metabolic network of H. sapiens embraces 1378 metabolites and 666 enzymes, in which 948 metabolites and 614 enzymes are included in the biggest connected cluster, while the other 52 enzymes and 430 metabolites scatter in 124 small clusters. Like previous studies of metabolic networks based on topology [2530], we only analyze the biggest connected cluster.

Modularity and network decomposition

In this study, we applied the simulated annealing algorithm developed by Guimera and Amaral [28] to break up the metabolic network of H. sapiens into modules (The software Modul-w was kindly obtained from Guimera and Amaral). This algorithm identifies topological modules by maximizing the network's modularity metric through an exhaustive search, thus may generate the best decomposition of the network.

For a given decomposition of a network, the modularity metric is defined as the gap between the fraction of arcs within clusters and the expected fraction of arcs if the arcs are wired with no structural bias [59]:

M = i = 1 r [ e i i ( j e i j ) 2 ] MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGnbqtcqGH9aqpdaaeWbqaaiabcUfaBjabdwgaLnaaBaaaleaacqWGPbqAcqWGPbqAaeqaaaqaaiabdMgaPjabg2da9iabigdaXaqaaiabdkhaYbqdcqGHris5aOGaeyOeI0IaeiikaGYaaabuaeaacqWGLbqzdaWgaaWcbaGaemyAaKMaemOAaOgabeaaaeaacqWGQbGAaeqaniabggHiLdGccqGGPaqkdaahaaWcbeqaaiabikdaYaaakiabc2faDbaa@4807@

where r is the number of clusters, e ij is the fraction of arcs that leads between vertices of cluster i and j. The maximum modularity metric corresponds to the partition that comprises as many as within-module links and as few as possible inter-module links.

Since the simulated annealing algorithm is stochastic, different runs may yield different decompositions of the network because of different initial conditions. The initial condition is the seed for the random number generator, which must be a negative integer. To verify the robustness of the decomposition, we first randomly generated 30 different negative integers, and then applying each of them as the initial condition to perform 30 times of independent simulated annealing and generate 30 partitions of the network. Figure 7 shows the fraction of times that each pair of nodes in the network is clustered into the same module. It is thus confirmed that the topological modules are robustly and consistently identified.

Figure 7

Accuracy of simulated annealing algorithm to identify topological modules. We obtain 30 independent decompositions of the H. sapiens metabolic network and plot the fraction of times that each pair of nodes is clustered into the same module.

Hypergeometric distribution and P-value

If we randomly draw n samples from a finite population, the probability of getting i samples with the desired feature by chance obeys hypergeometric distribution:

f ( i ) = ( K i ) ( N K n i ) ( N n ) , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGMbGzcqGGOaakcqWGPbqAcqGGPaqkcqGH9aqpdaWcaaqaamaabmaabaqbaeqabiqaaaqaaiabdUealbqaaiabdMgaPbaaaiaawIcacaGLPaaadaqadaqaauaabeqaceaaaeaacqWGobGtcqGHsislcqWGlbWsaeaacqWGUbGBcqGHsislcqWGPbqAaaaacaGLOaGaayzkaaaabaWaaeWaaeaafaqabeGabaaabaGaemOta4eabaGaemOBa4gaaaGaayjkaiaawMcaaaaacqGGSaalaaa@43A8@

where N is the size of the population, K is the number of items with the desired feature in the population. Then the probability of getting at least k samples with the desired feature by chance can be represented by hypergeometric cumulative distribution defined as P-value:

P = 1 i = 0 k 1 f ( i ) = 1 i = 0 k 1 ( K i ) ( N K n i ) ( N n ) MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGqbaucqGH9aqpcqaIXaqmcqGHsisldaaeWbqaaiabdAgaMjabcIcaOiabdMgaPjabcMcaPiabg2da9aWcbaGaemyAaKMaeyypa0JaeGimaadabaGaem4AaSMaeyOeI0IaeGymaedaniabggHiLdGccqaIXaqmcqGHsisldaaeWbqaamaalaaabaWaaeWaaeaafaqabeGabaaabaGaem4saSeabaGaemyAaKgaaaGaayjkaiaawMcaamaabmaabaqbaeqabiqaaaqaaiabd6eaojabgkHiTiabdUealbqaaiabd6gaUjabgkHiTiabdMgaPbaaaiaawIcacaGLPaaaaeaadaqadaqaauaabeqaceaaaeaacqWGobGtaeaacqWGUbGBaaaacaGLOaGaayzkaaaaaaWcbaGaemyAaKMaeyypa0JaeGimaadabaGaem4AaSMaeyOeI0IaeGymaedaniabggHiLdaaaa@5A55@

Given significance level α, which is usually set as 0.05, a P-value smaller than α demonstrates low probability that the items with the desired feature are chosen by chance. Hence P-value can be used to measure whether the n samples drawn from the population is more enriched with items of the desired feature than would be expected by chance [51].

Phylogenetic profiles

For every enzyme, its phylogenetic profile is defined as a binary vector that encodes its absence (0) or presence (1) in the reference genomes. In this study, we applied the method in [16, 17] to construct the phylogenetic profiles of enzymes. We first chose 115 organisms (16 eukaryote, 83 bacteria and 16 archaea) as reference genomes. To reduce the effect of bias in the organism distribution, we then merged these 115 genomes into 54 taxa according to the NCBI taxonomy [16, 17, 52] [see Table A1 in Additional file 1]. At last, the ENZYME section of the KEGG LIGAND database [43] was utilized to determine whether the corresponding enzyme is coded in the reference genomes or not. The resulting phylogenetic profile of each enzyme is a 54-dimention binary vector.

We applied the Jaccard coefficient to measure the similarity between the phylogenetic profiles of any pair of enzymes. Jaccard coefficient is defined as follows [60],

J C i j = n i j n i + n j n i j , MathType@MTEF@5@5@+=feaafiart1ev1aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacH8akY=wiFfYdH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8kuc9pgc9s8qqaq=dirpe0xb9q8qiLsFr0=vr0=vr0dc8meaabaqaciaacaGaaeqabaqabeGadaaakeaacqWGkbGscqWGdbWqdaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabg2da9maalaaabaGaemOBa42aaSbaaSqaaiabdMgaPjabdQgaQbqabaaakeaacqWGUbGBdaWgaaWcbaGaemyAaKgabeaakiabgUcaRiabd6gaUnaaBaaaleaacqWGQbGAaeqaaOGaeyOeI0IaemOBa42aaSbaaSqaaiabdMgaPjabdQgaQbqabaaaaOGaeiilaWcaaa@441F@

where n ij is the number of taxa that encode both enzyme i and j; n i , n j are the number of taxa that encode enzyme i and j, respectively.

Evolutionary age

Also based on the ENZYME section of the KEGG LIGAND database, we classified the enzymes in the metabolic networks into the following seven evolutionary ages according to their inferred first appearance during evolution,

  1. 1)

    Prokaryota: E. coli group [see Table A1 of Additional file 1]

  2. 2)

    Protists: Plasmodium falciparum, Trypanosoma brucei, Entamoeba histolytica

  3. 3)

    Fungi: Saccharomyces cerevisiae; Schizosaccharomyces pombe

  4. 4)

    Nematodes: Caenorhabditis elegans

  5. 5)

    Arthropods: Drosophila melanogaster

  6. 6)

    Mammalian: Mus musculus; Rattus norvegicus; Canis familiaris

  7. 7)

    Human: H. sapiens

An enzyme was assigned to evolutionary age of "Prokaryota" if an ortholog was detected in any organism of the 15th taxonomy in Table A1 of Additional file 1 (the E. coli group); "Protists" if an ortholog was detected in P. falciparum, or T. brucei, or E. histolytica but not in E. coli group; "Fungi" if an ortholog was detected in S. cerevisiae or S. pombe but not in E. coli group, P. falciparum, T. brucei, and E. histolytica, and so forth.

Evolutionary rate

The evolutionary rate of an enzyme gene is defined as the normalized ratio of non-synonymous substitutions per nucleotide site (K a ) to synonymous substitutions per nucleotide site (K s ) that occurred in this enzyme gene [53]. We extracted the evolutionary rate of every enzyme gene in the H. sapiens metabolic network based on the complete genome of Pan troglodytes from HomoloGene database [52]. Nei and Gojobori's method [53] was used to calculate the synonymous and nonsynonymous substitution rates in HomoloGene database.


  1. 1.

    Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402: C47-C52. 10.1038/35011540.

  2. 2.

    Gavin AC, Bosche M, Krause R, Grandi P, Marzioch M, Bauer A, Schultz J, Rick JM, Michon AM, Cruciat CM: Functional organization of the yeast proteome by systematic analysis of protein complexes. Nature. 2002, 415: 141-147. 10.1038/415141a.

  3. 3.

    Ho Y, Gruhler A, Heilbut A, Bader GD, Moore L, Adams SL, Millar A, Taylor P, Bennett K, Boutilier K: Systematic identification of protein complexes in Saccharomyces cerevisiae by mass spectrometry. Nature. 2002, 415: 180-183. 10.1038/415180a.

  4. 4.

    Mewes HW, Frishman D, Guldener U, Mannhaupt G, Mayer K, Mokrejs M, Morgenstern B, Munsterkotter M, Rudd S, Weil B: MIPS: A database for genomes and protein sequences. Nucleic Acids Res. 2002, 30: 31-34. 10.1093/nar/30.1.31.

  5. 5.

    Nakao M, Bono H, Kawashima S, Kamiya T, Sato K, Goto S, Kanehisa M: Genome-scale gene expression analysis and pathway reconstruction in KEGG. Genome Informatics. 1999, 10: 94-103.

  6. 6.

    Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000, 28 (1): 27-30. 10.1093/nar/28.1.27.

  7. 7.

    Karp PD: Pathway databases: A case study in computational symbolic theories. Science. 2001, 293: 2040-2044. 10.1126/science.1064621.

  8. 8.

    Overbeek R, Larsen N, Pusch GD, D'Souza M, Jr ES, Kyrpides N, Fonstein M, Maltsev N, Selkov E: WIT: integrated system for high-throughput genome sequence analysis and metabolic reconstruction. Nucleic Acids Res. 2000, 28 (1): 123-125. 10.1093/nar/28.1.123.

  9. 9.

    Salgado H, Santos-Zavaleta A, Gama-Castro S, Millan-Zarate D, Diaz-Peredo E, Sanchez-Solano F, Perez-Rueda E, Bonavides-Martinez C, Collado-Vides J: RegulonDB (version 3.2): Transcriptional regulation and operon organization in Escherichia coli K-12. Nucleic Acids Res. 2001, 29: 72-74. 10.1093/nar/29.1.72.

  10. 10.

    van Nimwegen E, Zavolan M, Rajewsky N, Siggia ED: Probabilistic clustering of sequences: Inferring new bacterial regulons by comparative genomics. Proc Natl Acad Sci. 2002, 99: 7323-7328. 10.1073/pnas.112690399.

  11. 11.

    Fraser HB: Coevolution, modularity and human disease. Current Opinion in Genetics & Development. 2006, 16: 637-644. 10.1016/j.gde.2006.09.001.

  12. 12.

    Snel B, Huynen MA: Quantifying modularity in the evolution of biomolecular systems. Genome Res. 2004 , 14: 391-397. 10.1101/gr.1969504.

  13. 13.

    Campillos M, von Mering C, Jensen LJ, Bork P: Identification and analysis of evolutionarily cohesive functional modules in protein networks. Genome Res. 2006, 16: 374-382. 10.1101/gr.4336406.

  14. 14.

    Chen Y, Dokholyan NV: The coordinated evolution of yeast proteins is constrained by functional modularity . Trends in Genetics. 2006, 22 (8): 416-419. 10.1016/j.tig.2006.06.008.

  15. 15.

    Pereira-Leal JB, Levy ED, Teichmann SA: The origins and evolution of functional modules:lessons from protein complexes. Phil Trans R Soc B. 2006, 361: 507-517. 10.1098/rstb.2005.1807.

  16. 16.

    Yamada T, Kanehisa M, Goto S: Extraction of phylogenetic network modules from the metabolic network. BMC Bioinformatics. 2006, 7 (1): 130-10.1186/1471-2105-7-130.

  17. 17.

    Yamada T, Goto S, Kanehisa M: Extraction of phylogenetic network modules from prokayrote metabolic pathways. Genome Informatics. 2004, 15: 249-258.

  18. 18.

    von Mering C, Zdobnov EM, Tsoka S, Ciccarelli FD, Pereira-Leal JB, Ouzounis CA, Bork P: Genome evolution reveals biochemical networks and functional modules. Proc Natl Acad Sci. 2003, 100 (26): 15428-15433. 10.1073/pnas.2136809100.

  19. 19.

    Spirin V, Gelfand MS, Mironov AA, Mirny LA: A metabolic network in the evolutionary context: Multiscale structure and modularity. Proc Natl Acad Sci. 2006, 103 (23): 8774-8779. 10.1073/pnas.0510258103.

  20. 20.

    Spirin V, Mirny LA: Protein complexes and functional modules in molecular networks. Proc Natl Acad Sci. 2003, 100: 12 123-12 128. 10.1073/pnas.2032324100.

  21. 21.

    Rives AW, Galitski T: Modular organization of cellular networks. Proc Natl Acad Sci. 2003, 100: 1128-1133. 10.1073/pnas.0237338100.

  22. 22.

    Gagneur J, Krause R, Bouwmeester T, Casari G: Modular decomposition of protein-protein interaction networks. Genome Biol. 2004, 5: R57-10.1186/gb-2004-5-8-r57.

  23. 23.

    Ihmels J, Friedlander G, Bergmann S, Sarig O, Ziv Y, Barkai N: Revealing modular organization in the yeast transcriptional network. Nature Genetics. 2002, 31: 370 -3377.

  24. 24.

    Ma HW, Buer J, Zeng AP: Hierarchical structure and modules in the Escherichia coli transcriptional regulatory network revealed by a new top-down approach. BMC Bioinformatics. 2004, 5: 199-10.1186/1471-2105-5-199.

  25. 25.

    Ma HW, Zhao XM, Yuan YJ, Zeng AP: Decomposition of metabolic network into functional modules based on the global connectivity structure of reaction graph. Bioinformatics. 2004, 20 (12): 1870-1876. 10.1093/bioinformatics/bth167.

  26. 26.

    Gagneur J, Jackson DB, Casari G: Hierarchical analysis of dependency in metabolic networks. Bioinformatics. 2003, 19 (8): 1027-1034. 10.1093/bioinformatics/btg115.

  27. 27.

    Holme P, Huss M, Jeong H: Subnetwork hierarchies of biochemical pathways. Bioinformatics. 2003, 19 (4): 532-538. 10.1093/bioinformatics/btg033.

  28. 28.

    Guimera R, Amaral LAN: Functional cartography of complex metabolic networks. Nature. 2005, 433 (7028): 895-900. 10.1038/nature03288.

  29. 29.

    Ravasz E, Somera AL, Mongru DA, Oltvai ZN, Barabasi AL: Hierarchical Organization of Modularity in Metabolic Networks. Science. 2002, 297 (5586): 1551-1555. 10.1126/science.1073374.

  30. 30.

    Zhao J, Yu H, Luo J, Cao Z, Li Y: Hierarchical modularity of nested bow-ties in metabolic networks. BMC Bioinformatics. 2006, 7: 386-10.1186/1471-2105-7-386.

  31. 31.

    Girvan M, Newman MEJ: Community structure in social and biological networks. Proc Natl Acad Sci. 2002, 99 (12): 7821-7826. 10.1073/pnas.122653799.

  32. 32.

    Papin JA, Reed JL, Palsson BO: Hierarchical thinking in network biology: the unbiased modularization of biochemical networks. Trends in Biochemical Sciences. 2004, 29 (12): 641-647. 10.1016/j.tibs.2004.10.001.

  33. 33.

    Fraser HB, Hirsh AE, Steinmetz LM, Scharfe C, Feldman MW: Evolutionary Rate in the Protein Interaction Network. Science. 2002, 296: 750-752. 10.1126/science.1068696.

  34. 34.

    Wuchty S, Barabási AL, Ferdig MT: Stable evolutionary signal in a Yeast protein interaction network. BMC Evolutionary Biology. 2006, 6: 8-10.1186/1471-2148-6-8.

  35. 35.

    Vitkup D, Kharchenko P, Wagner A: Influence of metabolic network structure and function on enzyme evolution. Genome Biology. 2006, 7 (5): R39-10.1186/gb-2006-7-5-r39.

  36. 36.

    Light S, Kraulis P, Elofsson A: Preferential attachment in the evolution of metabolic networks. BMC Genomics. 2005, 6: 159-10.1186/1471-2105-6-159.

  37. 37.

    Fraser H, Wall D, Hirsh A: A simple dependence between protein evolution rate and the number of protein-protein interactions. BMC Evolutionary Biology. 2003, 3 (1): 11-10.1186/1471-2148-3-11.

  38. 38.

    Saeed R, Deane C: Protein protein interactions, evolutionary rate, abundance and age. BMC Bioinformatics. 2006, 7 (1): 128-10.1186/1471-2105-7-128.

  39. 39.

    Fraser HB: Modularity and evolutionary constraint on proteins. Nature Genetics. 2005, 37 (4): 351-352. 10.1038/ng1530.

  40. 40.

    Ekman D, Light S, Bjorklund A, Elofsson A: What properties characterize the hub proteins of the protein-protein interaction network of Saccharomyces cerevisiae?. Genome Biology. 2006, 7 (6): R45-10.1186/gb-2006-7-6-r45.

  41. 41.

    Goto S, Nishioka T, Kanehisa M: LIGAND: chemical database for enzyme reactions. Bioinformatics. 1998 , 14: 591-599. 10.1093/bioinformatics/14.7.591.

  42. 42.

    Goto S, Nishioka T, Kanehisa M: LIGAND: chemical database of enzyme reactions. Nucleic Acids Res. 2000, 28 (1): 380-382. 10.1093/nar/28.1.380.

  43. 43.

    Goto S, Okuno Y, Hattori M, Nishioka T, Kanehisa M: LIGAND: database of chemical compounds and reactions in biological pathways. Nucleic Acids Res. 2002, 30 (1): 402-404. 10.1093/nar/30.1.402.

  44. 44.

    Zhao J, Yu H, Luo J, Cao Z, Li Y: Complex networks theory for analyzing metabolic networks. Chinese Science Bulletin. 2006, 51 (13): 1529-1537. 10.1007/s11434-006-2015-2.

  45. 45.

    Borgatti SP, Everett MG: Models of core/periphery structures. Social Networks. 1999, 21: 375-395. 10.1016/S0378-8733(99)00019-2.

  46. 46.

    Holme P: Core-periphery organization of complex networks. Phys Rev E. 2005, 72: 46111-10.1103/PhysRevE.72.046111.

  47. 47.

    Tanaka T, Ikeo K, Gojobori T: Evolution of metabolic networks by gain and loss of enzymatic reaction in eukaryotes. Gene. 2006, 365: 88-94. 10.1016/j.gene.2005.09.030.

  48. 48.

    Csete M, Doyle J: Bow ties, metabolism and disease. Trends in Biotechnology. 2004, 22 (9): 446-450. 10.1016/j.tibtech.2004.07.007.

  49. 49.

    Dandekar T, Schuster S, Snel B, Huynen M, Bork P: Pathway alignment: application to the comparative analysis of glycolytic enzymes. Biochemical Journal. 1999, 343: 115-124. 10.1042/0264-6021:3430115.

  50. 50.

    Huynen MA, Dandekar T, Bork P: Variation and evolution of the citric-acid cycle: a genomic perspective . Trends in Microbiology. 1999, 7 (7): 281-291. 10.1016/S0966-842X(99)01539-5.

  51. 51.

    Tavazoie S, Hughes JD, Campbell MJ, Cho RJ, Church GM: Systematic determination of genetic network architecture. Nat Genet. 1999, 22 (3): 281-285. 10.1038/10343.

  52. 52.

    Wheeler DL, Barrett T, Benson DA, Bryant SH, Canese K, Chetvernin V, Church DM, DiCuccio M, Edgar R, Federhen S, Geer LY, Helmberg W, Kapustin Y, Kenton DL, Khovayko O, Lipman DJ, Madden TL, Maglott DR, Ostell J, Pruitt KD, Schuler GD, Schriml LM, Sequeira E, Sherry ST, Sirotkin K, Souvorov A, Starchenko G, Suzek TO, Tatusov R, Tatusova TA, Wagner L, Yaschenko E: Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2006, 34 (suppl_1): D173-180. 10.1093/nar/gkj158.

  53. 53.

    Nei M, Gojobori T: Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986, 3 (5): 418-426.

  54. 54.

    Zhao J, Tao L, Yu H, Luo JH, Cao ZW, Li Y: Bow-tie topological features of metabolic networks and the functional significance. Chinese Science Bulletin. 2007, 52: 1036-1045. 10.1007/s11434-007-0143-y.

  55. 55.

    Maslov S, Sneppen K: Specificity and stability in topology of protein networks. Science. 2002, 296 (5569): 910-913. 10.1126/science.1065103.

  56. 56.

    Maslov S, Sneppen K, Zaliznyak A: Detection of topological patterns in complex networks: correlation profile of the internet. Physica A: Statistical and Theoretical Physics. 2004, 333: 529-540. 10.1016/j.physa.2003.06.002.

  57. 57.

    Ma H, Zeng AP: Reconstruction of metabolic networks from genome data and analysis of their global structure for various organisms. Bioinformatics. 2003, 19 (2): 270-277. 10.1093/bioinformatics/19.2.270.

  58. 58.

    Huss M, Holme P: Currency and commodity metabolites: Their identification and relation to the modularity of metabolic networks. IET Systems Biology. 2007, 1: 280-285. 10.1049/iet-syb:20060077.

  59. 59.

    Newman MEJ, Girvan M: Finding and evaluating community structure in networks. Physical Review E. 2004, 69: 26113-10.1103/PhysRevE.69.026113.

  60. 60.

    Lamboy WF: Computing genetic similarity coefficients from RAPD data: the effects of PCR artifacts. PCR Methods Appl. 1994, 4 (1): 31-37.

Download references


We thank Dr. Luís A. Nunes Amaral and Dr. Roger Guimerà for kindly providing us the software Modul-w for network decomposition; Dr. Jingchu Luo and Dr. Mikael Huss for suggestive comments on the manuscript; and the anonymous reviewers for their helpful comments. JZ thanks Dr. Petter Holme and Dr. Mikael Huss for stimulating discussions about modularity. This work was supported in part by grants from Ministry of Science and Technology China (2006AA02Z317, 2004CB720103, 2003CB715901, 2006AA02312), National High Technology Research and Development Program of China (2006AA020805), National Natural Science Foundation of China (30500107, 30670953, 30670574), International Cooperation Project of Science and Technology Commission of Shanghai Municipality (06RS07109), and Grant from Science and Technology Commission of Shanghai Municipality (04DZ19850, 06PJ14072, 04DZ14005).

Author information

Correspondence to Zhi-Wei Cao or Yi-Xue Li.

Additional information

Authors' contributions

JZ conceived of the study, designed the analysis, implemented the analysis and prepared the manuscript. GHD implemented part of the analysis and helped to revise the manuscript. LT, HY and ZHY helped JZ to implement the analysis. JHL managed the project. ZWC and YXL helped JZ to design the analysis, provided guidance, coordinated and participated in the biological and theoretical analyses, and revised the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Zhao, J., Ding, G., Tao, L. et al. Modular co-evolution of metabolic networks. BMC Bioinformatics 8, 311 (2007).

Download citation


  • Metabolic Network
  • Simulated Annealing Algorithm
  • Enzyme Gene
  • Core Module
  • Phylogenetic Profile