PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. May 2011; 7(5): e1001137.
Published online May 19, 2011. doi:  10.1371/journal.pcbi.1001137
PMCID: PMC3098196
A Toolbox Model of Evolution of Metabolic Pathways on Networks of Arbitrary Topology
Tin Yau Pang1,2 and Sergei Maslov1*
1Department of Condensed Matter Physics and Materials Science, Brookhaven National Laboratory, Upton, New York, United States of America
2Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York, United States of America
Eugene I. Shakhnovich, Editor
Harvard University, United States of America
* E-mail: maslov/at/bnl.gov
Conceived and designed the experiments: SM. Performed the experiments: TYP. Analyzed the data: TYP SM. Wrote the paper: TYP SM. Derived mathematical solutions: SM.
Received September 20, 2010; Accepted April 14, 2011.
In prokaryotic genomes the number of transcriptional regulators is known to be proportional to the square of the total number of protein-coding genes. A toolbox model of evolution was recently proposed to explain this empirical scaling for metabolic enzymes and their regulators. According to its rules, the metabolic network of an organism evolves by horizontal transfer of pathways from other species. These pathways are part of a larger “universal” network formed by the union of all species-specific networks. It remained to be understood, however, how the topological properties of this universal network influence the scaling law of functional content of genomes in the toolbox model. Here we answer this question by first analyzing the scaling properties of the toolbox model on arbitrary tree-like universal networks. We prove that critical branching topology, in which the average number of upstream neighbors of a node is equal to one, is both necessary and sufficient for quadratic scaling. We further generalize the rules of the model to incorporate reactions with multiple substrates/products as well as branched and cyclic metabolic pathways. To achieve its metabolic tasks, the new model employs evolutionary optimized pathways with minimal number of reactions. Numerical simulations of this realistic model on the universal network of all reactions in the KEGG database produced approximately quadratic scaling between the number of regulated pathways and the size of the metabolic network. To quantify the geometrical structure of individual pathways, we investigated the relationship between their number of reactions, byproducts, intermediate, and feedback metabolites. Our results validate and explain the ubiquitous appearance of the quadratic scaling for a broad spectrum of topologies of underlying universal metabolic networks. They also demonstrate why, in spite of “small-world” topology, real-life metabolic networks are characterized by a broad distribution of pathway lengths and sizes of metabolic regulons in regulatory networks.
It has been previously reported that in prokaryotic genomes the number of transcriptional regulators is proportional to the square of the total number of genes. We recently offered a general explanation of this empirical powerlaw scaling in terms of the “toolbox” model in which metabolic and regulatory networks co-evolve together. This evolution is driven by horizontal gene transfer of co-regulated metabolic pathways from other species. These pathways are part of a larger “universal” network formed by the union of all species-specific networks. In the present work we address the question of how topological properties of this universal network influence the powerlaw scaling of regulators in the toolbox model. We also generalize its rules to include reactions with multiple substrates and products, branched and cyclic metabolic pathways, and to account for optimality of metabolic pathways. The main conclusion of our analytical and numerical modeling efforts is that the quadratic scaling is the robust feature of the toolbox model in a broad range of universal network topologies. They also demonstrate why, in spite of “small-world” topology, real-life metabolic networks are characterized by a broad distribution of pathway lengths and sizes of regulons in regulatory networks.
In prokaryotic genomes the number of transcriptional regulators is known to quadratically scale with the total number of protein-coding genes [1]. The toolbox model of co-evolution of metabolic and regulatory networks was recently proposed [2] to explain this scaling in parts of the genome responsible for metabolic functions. In this model prokaryotes acquire new metabolic capabilities by horizontal transfer of entire metabolic pathways from other organisms. One can conveniently think of these new pathways as coming from some “universal network” formed by the union of metabolic repertoires of all potential donor organisms. The essence of the toolbox argument [2] can be summarized as follows: as the non-regulatory part of the genome of an organism (its “toolbox”) grows, it typically needs to acquire fewer and fewer extra new genes (“tools”) in a pathway offering it some new metabolic capability (e.g. the ability to utilize a new nutrient or synthesize a new metabolic product). As a consequence, the number of pathways and by extension the number of their transcriptional regulators grows faster than linearly with the number of non-regulatory genes in the genome. While this qualitative explanation is rather general and therefore does not depend on specific details such as topology of the universal network, the exact value of the exponent α connecting the number of transcription factors (equal to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e001.jpg- the number of pathways or leaves of the network) to the number of metabolites in the metabolic network of an organism An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e002.jpg, as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e003.jpg, is in general model-dependent. In [2] we mathematically derived the quadratic scaling (An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e004.jpg) for the toolbox model with linear pathways on a fully connected graph in which any pair of metabolites can in principle be converted to each other in just one step via a single metabolic reaction. While this situation is obviously unrealistic from biological standpoint, before present study it remained the only mathematically treatable variant of the toolbox model. The universality of the exponent An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e005.jpg was then corroborated [2] by numerical simulations of the toolbox model with linearized pathways on the universal network formed by the union of all metabolic reactions in the KEGG database. While the agreement between the values of the exponent An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e006.jpg in these two cases hinted at underlying general principles at work, the detailed understanding of these principles remained elusive.
The question we address in this study is how the topology of the universal network determines this scaling exponent. To answer this question we first consider and solve a more realistic (yet still mathematically treatable) case in which the universal metabolic network is a directed tree of arbitrary topology. While being closer to reality than previously solved [2] case of fully connected network, the toolbox model on a tree-like universal network still retains a number of simplifications such as strictly linear pathways and one substrate → one product reactions.
To make our approach even more realistic we propose and numerically study a completely new version of the toolbox model incorporating metabolic reactions with multiple substrates and products as well as branched and cyclic metabolic pathways. Furthermore, unlike random linear pathways on a universal network [2] that can be long and therefore suboptimal from an evolutionary standpoint, the new model uses evolutionarily optimized pathways with the smallest number of reactions from the KEGG database sufficient to achieve a given metabolic task.
The toolbox model on a tree-like universal network: General mathematical description
We will first consider the case where the universal metabolic network is a directed tree. For simplicity in this section we will consider the case of catabolic pathways, while identical arguments (albeit with opposite direction of all reactions) apply to anabolic pathways. The root of the tree corresponds to the central metabolic core of the organism responsible for biomass production. Peripheral catabolic pathways (branches of the tree) convert external nutrients (leaves) to this core, while the internal nodes of the tree represent intermediate metabolites. Each of metabolites is characterized by its distance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e007.jpg from the root of the network. The universal network has An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e008.jpg metabolites at distance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e009.jpg from the root that included An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e010.jpg leaves (nutrients used in the first step of catabolic pathways) and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e011.jpg branching points corresponding to intermediate metabolites generated by more than one metabolic reaction at the next level (see Figure 1). An organism-specific network (filled circles and thick edges in Figure 1) at distance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e012.jpg from the root contains An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e013.jpg metabolites composed of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e014.jpg leaves, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e015.jpg branching points, and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e016.jpg metabolites inside linear branches (“one reaction in-one reaction out”) . For simplicity we assume that in the universal network (and thus also in any of its organism-specific subnetworks) no more than two reaction edges can combine at any node (metabolite), while the most general case of an arbitrary distribution of branching numbers can be treated in a very similar fashion.
Figure 1
Figure 1
An example of organism-specific metabolic network and the corresponding universal network.
The toolbox model specifies rules by which organism acquires new pathways in the course of its evolution. It consists of the following steps: 1) randomly pick a new nutrient metabolite (a leaf node of the universal network that currently does not belong to the metabolic network of the organism) 2) use the universal network to identify the unique linear pathway which connects the new nutrient to the root of the tree (the metabolic core) and finally 3) add the reactions and intermediate metabolites in the new pathway to the metabolic network of the organism (filled circles and thick edges in Figure 1). One needs to only add those enzymes that are not yet present in the “genome” of the organism. Graphically it means that the new branch of the universal network is extended until it first intersects the existing metabolic network of the organism.
Consider an organism capable of utilizing An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e025.jpg nutrients represented by leaves in the universal network, where An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e026.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e027.jpg. Since we assume that each nutrient utilization pathway is controlled by a dedicated transcriptional regulator sensing its presence or absence in the environment (e.g. LacR for lactose), the corresponding regulatory network would also have An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e028.jpg transcription factors (in the model we ignore transcription factors controlling non-metabolic functions). The non-regulatory part of the genome consists of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e029.jpg enzymes catalyzing metabolic reactions (strictly speaking An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e030.jpg is the number of metabolites/nodes so that the number of enzymes/edges is An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e031.jpg). Quadratic scaling plots [1] shows the number of transcriptional regulators An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e032.jpg vs. the total number of genes in the genome (both regulatory and non-regulatory) An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e033.jpg. However, since in all organism-specific networks NM [dbl greater-than sign] NL, the quadratic scaling between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e034.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e035.jpg is equivalent to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e036.jpg.
We further assume that due to random selection An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e037.jpg nutrients are expected to be uniformly distributed among all d levels. Therefore, the expected number of leaves at a given level is given by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e038.jpg where the fraction An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e039.jpg is the same at all levels. On the other hand the fraction An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e040.jpg varies from level to level. It usually tends to increase as one gets closer towards the root of the tree reaching 1 for d = 0 (the root node itself). To derive the equation for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e041.jpg, one first notices that each of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e042.jpg metabolites at level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e043.jpg is converted to another intermediate metabolite at level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e044.jpg. Due to merging of pathways at An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e045.jpg branching points the number of unique intermediate metabolites at the level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e046.jpg is actually smaller: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e047.jpg. To calculate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e048.jpg one uses the fact that each of the two nodes downstream of a branching point in the universal network is present in the organism-specific network with probability An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e049.jpg. The probability that they are both present is An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e050.jpg and thus the number of branching points at level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e051.jpg of the organism-specific metabolic network is An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e052.jpg. The intermediate metabolites together with new nutrients An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e053.jpg entering at the level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e054.jpg add up to the total number of metabolites at level An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e055.jpg:
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e056.jpg
(1)
This equation allows one to iteratively calculate An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e057.jpg for all d starting from An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e058.jpg. We will use this equation to derive the relationship between the number of leaves and the total number of nodes first for a critical branching tree and then for a supercritical one.
The toolbox model on a critical tree
The Galton-Watson branching process [3] is a simple stochastic process generating random trees, and we will consider its version where a node can have two, one, or zero neighbors (parents) at the previous level with probabilities p2, p1 and p0 correspondingly. If the average number of parents An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e059.jpg equals one, then the process is referred to as critical, and if An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e060.jpg is greater than one then the process is supercritical. More generally critical and supercritical branching trees can be generated by a variety of random processes such as e.g. directed percolation [4]. While for simplicity we used the Galton-Watson branching process in our derivation below, it can be readily extended to this more general case.
The principal geometric difference between supercritical and critical trees is that in the former case the number of nodes in a layer An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e061.jpg exponentially grows with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e062.jpg [3], while in a critical tree it grows at most algebraically (for the Galton-Watson critical process An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e063.jpg [3]). The other difference is that while the critical branching process always stops on its own at a certain finite height An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e064.jpg, a supercritical process will go on forever so that to generate a tree one has to manually terminate it at a predefined layer An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e065.jpg. The most significant feature of a critical tree is that it has much longer branches than a supercritical one of the same size. Indeed, the diameter (the maximal height) of a random critical tree with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e066.jpg nodes is An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e067.jpg while in a supercritical tree it is much shorter: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e068.jpg. Thus supercritical trees (unlike their critical counterparts) have the small world property.
A random critical network where each node has at most has two parents in the previous layer is defined by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e069.jpg. Indeed, in this case An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e070.jpg. In such network An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e071.jpg and hence the Eq. (1) can be rewritten as
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e072.jpg
(2)
A critical branching process that has not terminated by level d satisfies An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e073.jpg or An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e074.jpg. More generally if An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e075.jpg algebraically increases with An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e076.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e077.jpg asymptotically approaches 1 as
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e078.jpg
(3)
Here An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e079.jpg as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e080.jpg, thus for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e081.jpg An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e082.jpg remains approximately constant and according to Eq. (2) this constant ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e083.jpg is defined by
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e084.jpg
(4)
This quadratic relation is exact in a critical branching tree where each node can branch out into at most two nodes at the next layer, and it is still correct to a leading order in An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e085.jpg for a critical branching tree with arbitrary branching ratios (see “Quadratic relation between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e086.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e087.jpg for general critical branching processes” of Text S1). Furthermore, one can show (see “Calculation of the average An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e088.jpg in the toolbox model on a critical tree” of Text S1) that in large critical networks the overall fraction of metabolites present in organism-specific metabolic network is very close to this stationary limit of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e089.jpg: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e090.jpg.
As was explained in the previous section the ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e091.jpg between the total number An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e092.jpg of metabolic-related genes in the genome of an organism and its theoretical maximal value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e093.jpg for a genome containing the entire universal network is also given by An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e094.jpg. Furthermore, in our model the number of leaves is equal to the number of nutrient-utilizing pathways or, alternatively, their transcriptional regulators An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e095.jpg . Thus like in a much simpler model of Ref. [2] the toolbox model on any critical tree-like universal network gives rise to quadratic scaling of the number of transcription factors with the total number of genes:
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e096.jpg
(5)
The geometrical properties of the universal network such as its total number of nodes/edges An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e097.jpg and number of leaves/branches An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e098.jpg determine the prefactor of this scaling law. Simulation of the toolbox model on the critical tree (Figure 2) verified our mathematical predictions with the best fit to binned datapoints in Figure 2 giving the exponent α = 1.9±0.1.
Figure 2
Figure 2
An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e099.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e100.jpg.
The toolbox model on a supercritical tree
For a supercritical branching process An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e106.jpg and according to Eq. (1) (See SI for the derivation) the steady state value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e107.jpg of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e108.jpg satisfies
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e109.jpg
(6)
Here An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e110.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e111.jpg. Notice that for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e112.jpg one has two solutions for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e113.jpg: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e114.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e115.jpg. This indicates transition in which for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e116.jpg exactly at zero one has An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e117.jpg, while for an arbitrary small yet positive An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e118.jpg the value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e119.jpg asymptotically converges to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e120.jpg for An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e121.jpg. This transition resembles the first order phase transition, e.g., liquid-gas transition, where right at the transition point very small variation of the external parameter such as temperature (analogous to An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e122.jpg in this model) results in a large jump of the order parameter such as density (analogous to our An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e123.jpg). (See [5] for details), The number of layers over which this conversion is taking place is itself a function of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e124.jpg and for small An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e125.jpg it is large. For exponentially growing supercritical networks and for small An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e126.jpg, the network average value of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e127.jpg defined as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e128.jpg satisfies
A mathematical equation, expression, or formula.
 Object name is pcbi.1001137.e129.jpg
(7)
Note that this equation connecting An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e130.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e131.jpg (see SI for detailed derivation) is markedly different from Eq. (6) for steady state value An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e132.jpg in middle layers.
In conclusion, while the toolbox model on a critical universal network is characterized by a quadratic scaling between An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e133.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e134.jpg (see Eq. (4)), the same model on a supercritical, exponentially expanding universal network gives rise to a linear scaling of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e135.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e136.jpg albeit with logarithmic corrections (see Eq. (7)). This difference in exponent equally applies to the scaling of the number of regulators An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e137.jpg vs. the total number of genes An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e138.jpg in the toolbox model on critical and supercritical universal network.
Simulation of the toolbox model on the KEGG network with linearized pathways
To test our mathematical results for a more realistic version of the universal tree we linearized pathways and reactions in the network formed by the union of all reactions in the KEGG database [6]. To this end we generated a random spanning tree on the KEGG network by the following algorithm: the metabolite pyruvate was selected as the root of the tree. We then randomly picked a metabolite located upstream of it and generated a linear pathway (tree branch) as a self-avoiding random walk on the KEGG network extended until it either merges with another pathway or reaches the root of the tree. This step was repeated until all upstream metabolites were covered. The resulting spanning tree was then used as the universal network on which we simulated the toolbox model by gradually increasing the number of pathways An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e139.jpg and recording the total number of metabolites An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e140.jpg in organism-specific metabolic networks. Our numerical simulations generated approximately quadratic scaling An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e141.jpg (see Ref. [2]).
To better understand the origins of this scaling we investigated the topology of the underlying universal tree. The criticality of a tree is defined by the asymptotic value of the ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e142.jpg for large An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e143.jpg: for supercritical trees it reaches a limit An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e144.jpg, while for critical ones it converges to 1 as described in Eq. (3). Figure 3 showing An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e145.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e146.jpg in the linearized KEGG network convincingly demonstrates its criticality. Thus the quadratic scaling between the number of transcriptional regulators and the number of metabolites in the toolbox model simulated on the linearized KEGG network is explained by the mathematical formalism described in previous sections.
Figure 3
Figure 3
An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e147.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e148.jpg for KEGG-based universal network with linearized pathways.
In addition to using a random spanning tree to linearize the KEGG network we also tried a version using minimal paths. In this version the universal network is generated by randomly picking a metabolite and connecting it to the root of the tree (pyruvate) by the shortest path. At a first glance such “minimal path” selection appears to be reasonable from evolutionary standpoint. Indeed, evolution would favor simpler and shorter pathways in order to minimize the expenditure of resources to achieve a given metabolic goal [7] . However, the minimal paths version of linearization of the KEGG resulted in a supercritical universal network with logarithmically short branches An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e152.jpg. As predicted for supercritical trees (Eq. (7)) the toolbox model in this case had an approximately linear scaling of the number of transcriptional regulators (leaves of branches on the network) with the total number of metabolites: the measured best fit exponent was only An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e153.jpg.
How do we reconcile the evolutionary pressure apparently selecting for minimal pathways with dramatically wrong scaling properties of this model? We believe that most of the ultra-short “small world” pathways generated by minimal paths on the KEGG network are unrealistic from biochemical standpoint. Indeed, highly connected co-factors often position metabolites with very different chemical formulas in close proximity to each other. For example, the KEGG reaction R00134: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e154.jpg would appear as a miraculous “one-step” conversion of carbon dioxide into formate, while the reaction R03546: An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e155.jpg would artificially link carbon dioxide and cyanate. The combination of these two reactions gives rise to equally impossible two-step path: formate → CO2 → cyanate. As a consequence of such artificial shortcuts branches of the universal network linearized by minimal paths are much shorter than they are in reality. .This problem is at least partially alleviated by 1) removing unusually high-degree nodes corresponding to common co-factors such as H2O, ATP, NAD in the metabolic network so that some unrealistic paths are eliminated, and also 2) using random spanning tree instead of the shortest paths. In Ref. [2] we followed both of these recipes to successfully reproduce the quadratic scaling in real-life data. Still no linearization procedure could completely avoid biochemically meaningless shortcuts. In the next section we introduce and study a new considerably more realistic version of the toolbox model operating on branched and interconnected universal networks. Pathways in this version of the toolbox model satisfy the evolutionary requirements for minimal size. Proper treatment of metabolic reactions with multiple substrates prevents biochemically meaningless shortcuts and as a consequence restores the quadratic scaling.
The toolbox model on KEGG network with branched pathways and multi-substrate reactions
Real metabolic reactions routinely include multiple inputs (substrates) and multiple outputs (products) (see Table 1 and Table 2 for statistics in the KEGG database). Furthermore, metabolic networks often have two or more alternative pathways generating the same set of end-products from the same set of nutrients. Both these factors result in metabolic networks that are branched and interconnected. Here we propose and simulate a more realistic version of the toolbox model. The most prominent feature of the new model of pathways is the “AND” function acting on inputs of multi-substrate reactions. It reflects the constraint that a reaction cannot be carried out unless all its substrates are present.
Table 1
Table 1
The distribution of irreversible reactions classified by their numbers of substrates and products.
Table 2
Table 2
The distribution of reversible reactions classified by their numbers substrates/products.
The new version of the toolbox model simulates addition of anabolic pathways aimed at production of new metabolites from those the model organism can currently synthesize (its current metabolic core). The new pathways are optimal in the sense that they contain the smallest number of reactions necessary to synthesize the desired end-product. As for previous versions of the toolbox model, one can modify the rules of this model to apply to catabolic pathways but for simplicity we will limit the following discussion to anabolic pathways. The rules of the new model are:
  • At the beginning of the simulation, the model organism starts with a “seed” metabolic network consisting of 40 metabolites classified by the KEGG as parts of central carbohydrate metabolism, plus a number of “currency” metabolites such as water, ATP and NAD (see the section “Seed metabolites of the scope expansion” of Text S1 for additional details). It is assumed that our organism is able to generate all of these metabolites by some unspecified catabolic pathways.
  • At each step a new metabolite that cannot yet be synthesized by the organism is randomly selected from the “scope” [8] of our seed metabolites. This scope consists of all metabolites that in principle could be synthesized from the seed metabolites using all reactions listed in the KEGG database (see Ref. [8] for details).
  • To search for the minimal pathway that converts core metabolites to this target we first perform the “scope expansion” [8] of the core until it first reaches the target. In the course of this expansion reactions and metabolites are added step by step (or layer by layer). Each layer consists of all KEGG reactions that have all their substrates among the metabolites in the current metabolic core of the organism (light blue area in Figure 4) and those generated by reactions in all the previous layers. (See Figure 4 for an illustration).
    Figure 4
    Figure 4
    Diagram of a new pathway added to the metabolic network of the organism.
  • Next we trace back added reactions starting from the target and progressively moving to lower levels. One starts by finding the reaction responsible for fabrication of the target metabolite and adding it to the new pathway (if several such reactions exist in the last layer we randomly choose one of them). In case of multi-layer expansion process some substrates of this reaction are not among the core metabolites (otherwise this reaction would be in the first layer). One then goes down one layer and adds the reactions fabricating these missing substrates. This is repeated all the way down to the first level of the original expansion. The resulting pathway includes the minimal (or nearly minimal) set of reactions needed to generate the target metabolite from the current metabolic core of the organism. Starting from the next step of the model the target and all intermediate metabolites become part of the metabolic core. Genes for enzymes catalyzing these new reactions are assumed to be horizontally transferred to the genome of the organism. The newly added metabolic pathway is assumed to have a dedicated transcriptional regulator so that the number of transcription factors in our model is always equal to the number of pathways or their target metabolites.
  • Steps 1–5 are repeated until metabolic network of the organism reaches its maximal size. At this stage it includes the entire scope [8] of the starting set of metabolites in step 1.
Numerical simulation of this model shows that the number of transcriptional regulators scales with the number of metabolites with power An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e156.jpg (Figure 5). This is consistent with quadratic scaling we observed and mathematically derived for a simpler model with linearized pathways composed of single-substrate reactions.
Figure 5
Figure 5
An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e157.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e158.jpg of toolbox model with branched pathways and multi-substrate reactions.
The mathematical formalism derived in the previous sections is limited to tree-like universal networks and thus does not directly apply to the new model. Nevertheless, one generally expects the quadratic scaling to be limited only to critical, “large world” networks in which organisms with small genomes initially tend to acquire sufficiently long pathways. As noted before, from purely topological standpoint the KEGG network has a “small world” property making long pathways unlikely. It is important to check if the realistic treatment of multi-substrate reactions did in fact restore the “large world” property and criticality to the KEGG universal network by increasing the minimal number of steps required for connecting target metabolites to the metabolic core. To quantify the criticality of the expansion process as before we use the ratio An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e163.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e164.jpg denotes the number of metabolites reached at step An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e165.jpg of the scope expansion starting from the initial seed subset of metabolites. As in the case of critical branching trees this ratio asymptomatically converges to 1 thus confirming the criticality of the scope expansion process. The mere existence of ~40 steps in this process (the x-axis in Figure 6) can serve as evidence in favor of “large world” character of the KEGG universal network characterized by the existence of long pathways.
Figure 6
Figure 6
An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e166.jpg vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e167.jpg for the universal network consisting of all KEGG reactions.
Geometrical properties of branched pathways in the model
Unlike linearized pathways in the original version of the toolbox model [2], branched pathways in the more realistic model from previous section are interesting objects in their own right. We identified several geometrical properties of these pathways (see Figure 4 for illustration) quantifying their position relative to the core network to which they were added: 1) nborder rxn–the number of added reactions that are connected (as a substrate or a product) with at least one metabolite in the core, 2) nbase–the number of metabolites in the core that serve as substrates to reactions in the added pathway, 3) nfeedback–the number of core metabolites that are products of reactions in the new pathway, 4) nbyproduct–the number of final metabolic products of the added pathway that are neither core metabolites nor the target, 5) length-the number of steps (layers of the scope expansion process) it takes to transform core metabolites into the target product. 4 illustrates the definition of these parameters while Figure 7 and Figure 8 plot these parameters as a function of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e170.jpg (the number of metabolites in the added pathway) or An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e171.jpg (the number of reactions in the added pathway).
Figure 7
Figure 7
nbyproduct vs.An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e172.jpg.
Figure 8
Figure 8
Various linear relationships on the individual pathways.
Approximately linear relationship between nborder rxn vs. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e178.jpg (Figure 8a) suggests that added pathways tend to be located at or near the surface of the core metabolic network of the organism. Most of reactions in these pathways use metabolites from this core network either as substrates (nbase) or as products (nfeedback). Further analysis indicates that “currency metabolites” (common co-factors that serve as substrates or products of many reactions) constitute a significant fraction (~25%) of all core metabolites involved in border reactions (see the section “Analysis of the currency metabolites in the toolbox model” of Text S1 for details). On the other hand, the fact that the number of steps in a pathway (its length) constitutes a good fraction of its overall number of reactions An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e179.jpg implies that, in spite of these numerous “shortcut” connections to the core, added pathways remain very thin and essentially linear. That is to say, these pathways tend to work as a single “conveyor belt” sequentially converting intermediate products into each other instead of having two or more parallel “processing lines” and assembling final products of these lines only at final stages of the pathway. This finding provides an intuitive reason why models with branched and linearized pathways have similar scaling properties. One can argue that this is because pathways in both models are essentially linear. Yet, in spite of their linearity and optimality (each has the smallest number of reactions to generate the target from the core) added pathways in the new version of the model are very different from shortest paths on the universal network. As illustrated in Figure 9 the average pathway length is several times longer than the geometrically shortest path between the target and the core.
Figure 9
Figure 9
Comparison of lengths of the pathways and shortest distances of the targets from the core.
As can be seen from Figure 7, most of added pathways (around 97%) do not generate any byproducts. They only produce the intended target and nfeedback metabolites in the core network of the organism to which they were added. The relative scarcity of byproducts suggests that pathways in our model satisfy the evolutionary constrains imposed on real-life organisms. Indeed, as previously proposed in Ref [9] it makes sense to assume that evolution favors pathways that achieve a given metabolic goal using the smallest number of enzymes and at the same time striving to generate the maximal possible yield. Unnecessary byproducts not only reduce the yield of the desired metabolic target, they also might become toxic in high concentrations and thus would require extra transporter proteins to pump them out.
The small world property of complex biomolecular networks has been extensively discussed in the literature during the last decade (see [10][12] for earliest reports in metabolic and protein interaction networks correspondingly). It was often assumed that the small world effect positively contributes to the robustness of the network by providing multiple redundant pathways for target production in metabolic networks or for propagation of signals along regulatory and protein interaction networks. In addition to its positive aspects the small world property in biomolecular networks also has a potentially negative side by facilitating system-wide propagation of undesirable cross-talk [13]. In the course of evolution different strategies appeared allowing organism to limit and attenuate these unwelcome side effects of global connectivity.
The extent of small world topology in metabolic networks has been recently disputed in [14]. There it was argued that many shortcuts in simple graph representations of metabolic networks are meaningless from biochemical standpoint. By taking into account additional structural information about metabolites Arita [14] dramatically increased the diameter of the metabolic network in E. coli. In our simulations of the toolbox model we also encountered limitations of the simple graph representation giving rise to small world topology of metabolic networks. Small world by definition implies very short pathways (or equivalently supercritical network branching with exponentially growing lists of neighbors at distance An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e180.jpg), which in its turn prevents the appearance of quadratic scaling in the linear toolbox model.
How to reconcile this apparent contradiction? The answer known from pioneering studies of R. Heinrich and collaborators (see e.g. [8], [15], [16] ) is to altogether abandon the simple graph representation in favor of realistic treatment of multi-substrate reactions. A metabolic reaction with two or more substrates will not proceed at any rate until all these metabolites are present in the cell. This implicit “AND” function operating on inputs of multi-substrate metabolic reactions makes reaching a given metabolic target much harder task and ultimately leads to dramatically longer pathways (Figure 9 quantifies this effect). These longer pathways in turn reinstate the quadratic scaling in the version of the toolbox model that was introduced in the previous section. This leads to the novel conclusion of our study that, when multi-substrate reactions are properly taken into account, the small world (supercritical) topology of the metabolic universe disappears in favor of the “large world” topology characteristic of critical branching networks. The increase in the effective diameter of the network due to this effect is dramatic. One goes from 3–4 steps diameter typical of a small world network of [12], [11] to ~8 steps of [14] and finally to 30–40 layers in the scope expansion process shown in Figure 6 (see also Figure 6 of [8]).
These arguments lead us to adapt the “scope expansion” algorithm by Heinrich et al [8] to pathway acquisition in the toolbox model. Not only did it restore the “large world” properties such as quadratic scaling to the model, it also made the added pathways plausible from evolutionary standpoint. Unlike linear random walk pathways on KEGG network used in [2], pathways in the new version of the toolbox model have the smallest number of KEGG reactions to achieve their metabolic task (production of the target metabolite from the set of metabolites already present in organism's network). As can be seen in Figure 7 a large fraction of these pathways also does not generate any byproducts. Accumulation of such byproducts inside a cell is potentially dangerous and would require specialized proteins to excrete them to the environment. The lack of byproducts also means that the useful yield of an added pathway is at or near its theoretical maximum. This is consistent with the fact that real biological pathways are optimized in the course of evolution to increase their yield while simultaneously reducing the number of reaction steps [7], [17], [18].
Optimality of metabolic pathways in central carbon metabolism was recently discussed in Ref. [17]. There it was shown that some (but not all) of these pathways coincide with the shortest walks in the space of possible metabolic transformations. This study also estimated a typical metabolic substrate can in principle be converted into any of the 20 different products in just one step. This quickly adds up to a very large number of biochemically feasible paths connecting metabolites to each other. However, this exponential growth does not necessarily result in a small world universal metabolic network. Indeed, evolutionary optimization leaves just a tiny fraction of these biochemically feasible reactions to be realized in any organism. The universal metabolic network formed by the union of all organism-specific metabolic networks is thus dramatically sparser than the set of all reactions allowed by the basic rules of biochemistry. Thus, as demonstrated in Ref. [8] and the present study, the number of metabolites one could generate in N steps starting from a small core network and using KEGG-listed metabolic reactions instead of expanding as An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e181.jpg grows with N much more slowly (algebraically). The overall picture consistent with both our observations and those of Ref. [17] is that exponentially large, supercritical tree of all possible biochemical transformations is first pruned to an evolutionary optimized critical universal network out of which individual organisms select a subset of reactions necessary to accomplish their metabolic goals: that is to utilize nutrients in their environment and generate metabolic targets necessary for their operation.
Simplified “toy” models based on artificial chemistry reactions have a long history of being used to reveal fundamental organizational principles of metabolic networks:
  • The recent model of Riehl et al [18] uses the simplest possible metabolites distinguished from each other only by the number of atoms of one element (e.g. carbon). All reactions in this case are of ligation/cleavage type (e.g. An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e182.jpg) constrained only by mass conservation. In spite of utmost simplicity of this artificial chemistry, the optimal pathways in this model display a surprisingly rich set of properties and bear some similarity to real-life metabolic pathways.
  • The study of Pfeiffer et el [19] emphasizes the role of different chemical groups forming metabolites. They consider another artificial chemistry where metabolites are defined by binary strings indicating presence or absence of each of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e183.jpg different chemical groups, and reactions transferring one such chemical group from one substrate which has it to another substrate which initially does not. Plausible evolutionary rules of their model give rise to complex scale-free metabolic networks emerging from the simple initial condition of An external file that holds a picture, illustration, etc.
Object name is pcbi.1001137.e184.jpg completely non-specific transferase enzymes.
  • Finally the artificial chemistry studied by Hintze et al [20] has molecules composed of three different types of atoms with different valences. Metabolites are linear molecules in which every atom is connected to others by as many bonds as specified by its valence. This model with rather complicated rules of evolution is then used to shed light on topics such as robustness and modularity of metabolic networks.
In our study we used the real-life (even if incomplete and sometimes noisy) metabolic universe of all reactions in the KEGG database. The only simplifying approximations remaining in the new most realistic version of the toolbox model is random selection of metabolic targets in the course of evolution and easy availability of any subset of KEGG reactions for horizontal transfer. Both these approximations can be relaxed in later versions of the model. Another promising direction is to extend the toolbox model to artificial chemistry universal networks of Refs. [18], [19], [20]. While taking away from the realism of the model such extensions might help to broaden our intuition about what topological properties of the universal network determine the scaling properties of its species-specific subnetworks.
The universal network used in our study consists of the union of all reactions listed in the KEGG database. The directionality of reactions and connected pairs of metabolites were inferred from the map version of the reaction formula: ftp.genome.jp/pub/kegg/ligand/reaction/reaction?mapformula.lst. The universal network with linearized pathways used to construct Figure 2 and Figure 3 consists of 1813 metabolites upstream of pyruvate and 2745 reaction edges out of which 1782 are irreversible and 963 are reversible. The metabolic network with branched and cyclic pathways used to construct Figure 599 consists of 1861metabolites located downstream from the central metabolism and reachable from it by the scope expansion algorithm of Ref. [8]. It has 2819 reactions out of which 1402 are irreversible and the remaining 1417 are reversible. Table 1 and Table 2 shows the statistics for the number of substrates and products of these reactions. The list of core metabolites is obtained from KEGG Pathways Modules in the category “central carbohydrate metabolism” and extended with “currency” metabolites including water, ATP and NAD. Simulations were done in Matlab and Octave.
Supporting Information
Text S1
Supplementary information.
(0.37 MB DOC)
Acknowledgments
Additional support of this work was provided by the DOE Systems Biology Knowledgebase project “Tools and Models for Integrating Multiple Cellular Networks”.
Footnotes
The authors have declared that no competing interests exist.
Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, US Department of Energy. Part of this study was supported by the program “Computational Biology and Bioinformatic Methods to Enable a Systems Biology Knowledgebase” ( DE-FOA-0000143) of the Office of Biological and Environmental Research, US Department of Energy. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
1. van Nimwegen E. Scaling laws in the functional content of genomes. Trends Genet. 2003;19:479–484. [PubMed]
2. Maslov S, Krishna S, Pang TY, Sneppen K. Toolbox model of evolution of prokaryotic metabolic networks and their regulation. Proc Natl Acad Sci U S A. 2009;106:9743–9748. [PubMed]
3. Athreya KB, Ney PE. New York: Dover Publications; 2004. Branching Processes.
4. Broadbent SR, Hammersley JM. Percolation Processes. Math. Proc. Camb. Phil. Soc. 1957;53:629–641.
5. Huang K. Statistical mechanics. New York: Wiley; 1987. 17.5 THE VAN DER WAALS EQUATION OF STATE.
6. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30. [PMC free article] [PubMed]
7. Meléndez-Hevia E, Torres NV. Economy of design in metabolic pathways: further remarks on the game of the pentose phosphate cycle. J. Theor. Biol. 1988;132:97–111. [PubMed]
8. Handorf T, Ebenhöh O, Heinrich R. Expanding metabolic networks: scopes of compounds, robustness, and evolution. J. Mol. Evol. 2005;61:498–512. [PubMed]
9. Meléndez-Hevia E, Torres NV. Economy of design in metabolic pathways: further remarks on the game of the pentose phosphate cycle. J. Theor. Biol. 1988;132:97–111. [PubMed]
10. Wagner A. The Yeast Protein Interaction Network Evolves Rapidly and Contains Few Redundant Duplicate Genes. Mol Biol Evol. 2001;18:1283–1292. [PubMed]
11. Wagner A, Fell DA. The small world inside large metabolic networks. Proc. Biol. Sci. 2001;268:1803–1810. [PMC free article] [PubMed]
12. Jeong H, Tombor B, Albert R, Oltvai ZN, Barabasi A. The large-scale organization of metabolic networks. Nature. 2000;407:651–654. [PubMed]
13. Maslov S, Ispolatov I. Propagation of large concentration changes in reversible protein-binding networks. Proc Natl Acad Sci U S A. 2007;104:13655–13660. [PubMed]
14. Arita M. The metabolic world of Escherichia coli is not small. Proc Natl Acad Sci U S A. 2004;101:1543–1547. [PubMed]
15. Handorf T, Christian N, Ebenhöh O, Kahn D. An environmental perspective on metabolism. J. Theor. Biol. 2008;252:530–537. [PubMed]
16. Ebenhöh O, Handorf T, Heinrich R. A cross species comparison of metabolic network functions. Genome Inform. 2005;16:203–213. [PubMed]
17. Noor E, Eden E, Milo R, Alon U. Central Carbon Metabolism as a Minimal Biochemical Walk between Precursors for Biomass and Energy. Mol. Cell. 2010;39:809–820. [PubMed]
18. Riehl WJ, Krapivsky PL, Redner S, Segrè D. Signatures of Arithmetic Simplicity in Metabolic Network Architecture. PLoS Comput Biol. 2010;6:e1000725. [PMC free article] [PubMed]
19. Pfeiffer T, Soyer OS, Bonhoeffer S. The evolution of connectivity in metabolic networks. PLoS Biol. 2005;3:e228. [PMC free article] [PubMed]
20. Hintze A, Adami C. Evolution of Complex Modular Biological Networks. PLoS Comput Biol. 2008;4:e23. [PubMed]
Articles from PLoS Computational Biology are provided here courtesy of
Public Library of Science