Topological Biosignatures: Large-Scale Structure of Chemical Networks from Biology and Astrochemistry

Abstract

The chemical basis of life involves more than simply the presence of biological molecules; biochemical systems embody a complex network of reactions with characteristic topological features. At the same time, chemical complexity is also present in nonbiological contexts, inviting us to clarify the relationship between chemistry and life through comparative studies. This study examines chemical networks from biology (the metabolism of E. coli) and astronomy (gas-phase reactions in dark molecular clouds) to establish common topological features that may be generic for any complex chemical system, as well as clear differences that may be topological signatures of life. The biological and astrochemical networks exhibit different scaling behaviors, and the network motifs found in the two systems show similarities as well as significant differences. The PageRank algorithm was used to quantify the degree to which individual species act primarily as products or reactants; in the metabolic network, these two roles are nearly identical for most species, whereas the astrochemical network shows a clearer partitioning into reactants and products. Key Words: Network theory—Metabolism—Interstellar medium—Astrochemistry. Astrobiology 12, 29–39.

1. Introduction

L ife on Earth can be defined largely in terms of its chemistry. A well-defined set of building blocks—amino acids, nucleic acids, lipids, sugars, and a variety of other small molecules—forms the material basis for life as we know it. In addition, this set of building blocks is dynamically interconnected by a complex network of metabolic pathways that allow cells to synthesize complex biomolecules out of simple precursors and, conversely, degrade biomass back into chemically simple components (Ma and Zeng, 2003; Tanaka, 2005). On the other hand, complex chemical networks are not unique to biology—planetary atmospheres (Solé and Munteanu, 2004; Wilson and Atreya, 2004) and dark molecular clouds (Millar, 2004; Woodall et al., 2007; Herbst and van Dishoeck, 2009) exhibit complex chemistry that takes place in the absence of biology as we know it. This raises a natural question: Which features of complex biochemical networks are inherently biological, as opposed to those that are generic features of complex chemical networks?

This distinction is of fundamental relevance to our understanding of life on Earth, as well as the search for life elsewhere in the Universe. To imagine what might be possible for biology in other planetary contexts, it can be helpful to cast the features of terrestrial biological systems in more abstract mathematical terms that may still be applicable even if the chemical details of life were dramatically different. In reality, only some of the features of living systems are outcomes of evolution by natural selection, whereas others are constrained by the physics of complex systems. The aim of the present study is to compare a complex biochemical network (E. coli metabolism) with a nonbiological chemical network of similar size (gas-phase reactions in dark molecular clouds) in a search for contrasts and commonalities. These two chemical systems exist in environments that are unimaginably different—one operates at around 310 K in condensed aqueous medium, whereas the other is active at 10 K at lower densities than can be achieved in a laboratory vacuum chamber—and topological elements that are shared between them may be generic features of complex chemical networks over which evolution has little control.

2. Theory

Throughout this study, chemical networks will be represented as directed bipartite graphs (Fig. 1H), meaning that they consist of two different types of nodes (representing reactions and species) and that edges are directional and exist only between nodes of different types. An edge pointing from a species node to a reaction node indicates that the species is a reactant; an edge pointing from a reaction node to a species node designates that species as one product of the reaction. Many previous network-theoretic studies of chemical networks have relied on projection-based methods in which all nodes represent species and edges represent participation in a common reaction (Jeong et al., 2000; Alves et al., 2002; Ma and Zeng, 2003; Arita, 2004; Zhao et al., 2007), or nodes represent reactions and edges denote a shared reactant or product (Wagner and Fell, 2001; Light and Kraulis, 2004; Light et al., 2005; Spirin et al., 2006; Kreimer et al., 2008). It has been shown, however (Montañez et al., 2010; Zhou and Nakhleh, 2011), that such projection schemes introduce statistical artifacts and discard potentially important information about the system's topology.

FIG. 1.

Illustration of graph properties. (A) Random network with a power-law degree distribution. Nodes are shaded according to their degree, from low (white) to high (black). Most nodes have low degree; a few centrally located nodes have very high degree. (B) Random network with an exponential degree distribution; shading scale is the same as in (A). Note the lack of extremely high-degree nodes; these tend to be much scarcer in exponential distributions. (C) Random network with a transitivity of zero. (D) Random network with a transitivity of 0.6. Note that interactions in (D) are much more “local” than those in C, because the connections made by each node are strongly correlated with those made by its neighbors. (E–G) Changes in network structure as the modularity increases from 0.3 (E) to 0.5 (F) to 0.7 (G). Each node is assigned to a module and shaded accordingly; in (G) these modules are clearly visible in the graph structure, whereas in (E) they are hardly significant. (H) A bipartite network, with nodes colored either white or black according to their partition. Note that edges exist only between white and black nodes, not within either group. The black node indicated by the light gray arrow has two white neighbors (i.e., a degree of two) and three black second-nearest neighbors. Its strength, however, is four, because one of those neighbors can be reached by two distinct paths. The black node indicated by the dark gray arrow has a degree of two, four second-nearest neighbors, and a strength of four.

One important traditional measure of the large-scale structure of networks is the degree distribution, which measures the number of connections that a randomly chosen node makes with others. Many biological networks, including metabolic networks, have been shown to have a power-law degree distribution in which the probability of a node's having degree k scales as P(k)∼k ^−γ, with (typically) 2<γ<3 (Barabási and Oltvai, 2004). Power-law distributions (Fig. 1A) have certain very special properties that are not shared by other types of distributions. For example, in an exponential system (Fig. 1B) where the degree probability distribution scales as P(k)∼e ^−αk, rescaling all the k values by a positive constant would require a modification of the decay constant α, whereas the power-law scaling P(k)∼k ^−γ is unaffected by a rescaling of k. Only power laws exhibit this scale-free behavior. The origin of this power-law scaling is somewhat more controversial—many authors borrow from the physics-oriented literature of self-organized criticality and refer to such networks as “scale-free” (Albert and Barabási, 2002; Barabási, 2009), while others (largely from an engineering background) refer to “scale-rich” networks that are highly optimized for tolerance to defined types of perturbations (Carlson and Doyle, 1999; Tanaka, 2005). It has been observed (Solé and Munteanu, 2004; Jolley and Douglas, 2010) that nonbiological chemical networks, such as those in dark molecular clouds and planetary atmospheres, often follow exponential degree scaling, which leads to a relative scarcity of highly connected nodes relative to the power-law case often seen in biology. One exception to this rule is that models of Earth's atmosphere, in contrast to atmospheric models for other planets in the Solar System, show a power-law degree distribution, which suggests that the same non-equilibrium physics that shapes biological systems may be shaping the planet as a whole (Solé and Munteanu, 2004).

For bipartite graphs, degree may not always be a useful measure. In chemical networks, the degree distribution of species is likely to be very diverse, as some species will participate in many reactions, and others will participate in far fewer. The degree distribution of reactions will be far less interesting—most will have a degree of 2–6, depending on their number of reactants and products. A more useful quantity in bipartite graphs is the strength (Fig. 1H), calculated as the number of paths from a given node to its second-nearest neighbors of the same type (Montañez et al., 2010). Note that this quantity can be larger than the number of second-nearest neighbors (which often corresponds to the degree in projected unipartite graphs) because a node could be connected to one of its second-nearest neighbors by paths through more than one nearest neighbor. This is one example of information that is discarded in the transition from bipartite to projected unipartite graphs.

Another unipartite measure that does not generalize well to bipartite graphs is the transitivity, defined as the probability that, if edges exist between nodes A and B and between nodes B and C, nodes A and C will also be connected. Transitivity can provide a measure of the “cliquishness” of a network (Fig. 1C, 1D) and has been used extensively in the sociological literature (Newman, 2003). Geometrically, a calculation of the transitivity involves counting up the number of triangles in the graph—in a highly transitive network, triangles could be considered a “network motif” (Alon, 2007) that is overrepresented relative to what would be expected in a similar random graph. Bipartite graphs have a transitivity of zero by definition; the presence of triangles would require a connection between two nodes of the same type. By analogy, however, one could count up the number of quadrilaterals in a graph to define a “bipartite transitivity” measure. Furthermore, the directionality of edges allows for the detection of seven different types of quadrilaterals (Fig. 2A), allowing for a more nuanced view of the manner in which a particular network is transitive.

FIG. 2.

(A) The seven different types of quadrilaterals present in a bipartite-directed network. Light gray circles represent species nodes, and dark gray circles represent reaction nodes. Note that most reactions have additional reactants and products beyond what is included in the motif; see Fig. 6 for examples of these motifs embedded in the network. (B) Frequency of directed quadrilaterals in real and randomized chemical networks. For each network, the solid bar represents the native (or modified) network, whereas the transparent bars represent networks that have been shuffled by 10n _e edge-swapping moves, where n _e is the number of edges. In all the chemical networks considered, quadrilateral Types 0 and 4 are strong network motifs; Types 5 and 6 may represent weak antimotifs. Note that quadrilateral Types 1–3 may also be network motifs for EcoCyc, but not for UDfA.

Degree distribution and transitivity are both calculated solely from local interactions within the graph; other measures give a more direct measurement of the graph's large-scale structure. One example is the modularity (Newman, 2006), defined as the degree to which the graph can be broken into disjoint modules for which intramodule edges are more prevalent than intermodule edges (Fig. 1E–G). Modularity is traditionally calculated by using community detection algorithms, which attempt to break the graph into successively smaller units until no more viable modules can be identified. Community detection does not require any modifications to work with bipartite chemical graphs; the modules returned are collections of species and reactions that form a self-contained subnetwork that is weakly connected to the rest of the system.

As a final example of a method that yields information about the global structure of a graph, bipartite chemical graphs can also be analyzed with the PageRank algorithm (Page et al., 1999; Langville and Meyer, 2005), which forms the algorithmic basis for the Google search engine. PageRank was initially developed as a method of assessing the “importance” of web pages by simulating the motion of an idealized random walker on a directed graph and ranking nodes based on the amount of time that the random walker spends there. To use a more physical analogy, nodes begin with an equal amount of PageRank and distribute it evenly among their outbound edges, until the distribution of PageRank becomes self-consistent. There are two ways for a node to have a high PageRank—one would be to have a large number of inbound connections (especially from highly ranked nodes), and the other would be for relationships with connecting nodes to be fairly exclusive, so that they do not share their donated PageRank broadly. It should be noted that, although PageRank is an attribute of individual nodes, it is strongly context-dependent—nodes with a similar local environment (i.e., similar degree and strength) can have dramatically different PageRank values, based on their position within the network as a whole (Fig. 3).

FIG. 3.

Illustration of results of the PageRank algorithm. All graphs are shaded such that nodes with the highest PageRank are black and those with the lowest are white. The direction of edges is indicated by a thicker line at the end toward which the edge is directed. (A) In a simple cycle graph, all nodes have equal PageRank—an idealized random walker is equally likely to visit each. (B) In a directed tree graph, an idealized random walker will always travel toward the periphery rather than remaining at the center, resulting in much higher PageRank values near these terminal nodes. (C) In this two-dimensional grid graph, the node with the highest PageRank (black, top row) has three incoming edges, one of which comes from the node with the second-highest PageRank (dark gray, top row). Other nodes with three or four incoming edges can be found, but their PageRank is lower because the PageRank of their neighbors is lower. (D) In this pseudofractal Dorogovtsev-Goltsev-Mendes graph (Dorogovtsev et al., 2002), almost all nodes have two incoming edges, but most have outgoing edges as well, causing their PageRank to “flow” toward the periphery. The high-PageRank nodes are the terminal states that do not provide any outbound links.

Because of the inherent directionality of the PageRank algorithm, it can be evaluated in two different directions to provide different types of information. In the forward direction, nodes receive PageRank from their inbound connections; highly ranked species would be important products (meaning that they are formed in important reactions), and highly ranked reactions would be important consumers (meaning that they consume important species). In the reverse direction, highly ranked species would be important reactants (consumed in important reactions), and highly ranked reactions would be important producers (generating important species).

To avoid becoming caught in traps of cyclically referencing nodes, the idealized random walker of the PageRank algorithm has a finite probability at each step of “teleporting” to a randomly chosen node and resuming her walk from that point. This behavior can be modified somewhat through the usage of a “personalization vector” that defines each node's probability of being chosen as a teleportation destination. When the algorithm is run on a chemical network in the forward direction, this vector should be used to specify an initial set of reactant species—nebular abundances of heavy atoms and H₂ in the case of the astrochemical network, or growth media components in the case of the biochemical network. When the algorithm is run in the reverse direction, the personalization vector should specify the final state of the system—either the equilibrium abundances obtained at the end of an astrochemical simulation or the measured components of E. coli biomass.

3. Data Sets

Gas-phase chemical reactions in dark molecular clouds were derived from the UMIST Database for Astrochemistry (UDfA), which contains 2937 gas-phase reactions between 225 distinct molecular species (Woodall et al., 2007). Grain-surface processes were not included, and reactions not viable at 10 K were removed. The ultimate basis for databases such as UDfA lies in studies of dark molecular clouds in which radio and IR spectroscopy are used to identify prevalent chemical species. The success of a chemical model depends on its ability to reproduce estimated abundances of observed chemical species based on the inferred age and initial composition of the cloud. UDfA is designed to be used in (pseudo-) time-dependent simulations and therefore contains rate constant information for all reactions; this rate information is crucial to the correct outcome of astrochemical simulations but was not used in this study.

Reaction rates in the UMIST database have been based, to the extent possible, on experimental measurements. Woodall et al. (2007) reported that ∼30% of rate constants in the 2006 release of UMIST are measured. About 66% are curated from the literature, including in particular the NIST Chemistry Webbook (http://webbook.nist.gov/chemistry), the IUPAC Subcommittee on Gas Kinetic Data Evaluation for Atmospheric Chemistry (http://www.iupac-kinetic.ch.cam.ac.uk), and the extensive database of gas-phase cation-molecule reactions by Anicich (2003). For the purposes of this study, the presence of a reaction in the chemical network is binary—either a reaction occurs at a non-negligible rate and involves species known or suspected to be present in the interstellar medium (ISM), or it does not. To a first approximation, we can regard the development of an astrochemical model as a process in which species are detected in the ISM by radio or IR observations and these species are then connected by reactions whose plausibility is established by laboratory measurements or theoretical considerations. Such a model can then be subjected to uncertainty propagation and sensitivity analysis (Wakelam et al., 2010), and parameters with high uncertainty can be adjusted to improve agreement with observations.

When the PageRank algorithm was evaluated with a personalization vector, initial abundances were chosen as specified by Woodall et al. (2007), and final abundances were obtained from the steady-state values at the end of a 10⁸ y simulation by using the UDfA reaction set and physical parameters appropriate for dark molecular clouds (10 K, 10⁴ cm⁻³, 10 mag visual extinction). A table with relative abundances at the end of the simulation is included in the Supplementary Information (Supplementary Information is available online at www.liebertonline.com/ast).

The metabolism of E. coli K-12 was modeled based on the information available in the online E co li encyclopedia (EcoCyc) database (Keseler et al., 2009). Compared to UDfA, EcoCyc was, in a sense, constructed in the opposite direction. Instead of observing abundances of chemical species and using prior chemical knowledge to infer a reaction system that would reproduce observed abundances, the EcoCyc developers began with the annotated genome of E. coli K-12, used the genome sequence to determine a set of enzyme-catalyzed reactions, and then added the species known to participate in these reactions. This different approach has led to two crucial differences between UDfA and EcoCyc. One is that EcoCyc does not contain rate information for the vast majority of reactions.

The second (and for the purposes of this study far more important) difference is that EcoCyc contains a significant number of generic reactions, in which the substrate of a particular enzymatic reaction is either poorly characterized or inherently nonspecific. This means that the EcoCyc database could not be used to naïvely build a network model; instead, these nonspecific reactions had to be identified and, in most cases, removed from the network. Many of these nonspecific reactions were in lipid metabolism, where enzymes will act on lipids or precursors without regard to the number of carbon atoms in their fatty acid chains. Because lipids make a significant contribution to the biomass of a bacterial cell, these nonspecific reactions were replaced by a series of specific reaction steps that lead to the synthesis of dipalmitoyl diphosphatidylethanolamine, a major constituent of E. coli's cell membrane (Ganong and Raetz, 1982); a diagram showing the added reactions is included in the Supplementary Information accompanying this article. In addition to species involved in generic reactions, species were removed if they were either completely disconnected from the largest (undirected) connected subgraph component (i.e., not integrated with the rest of the network) or if they were connected only via generic species or currency metabolites (ATP, NADPH, etc.). A list of all species removed from the network and the rationale for their removal are provided in the Supplementary Information accompanying this article. In addition, reaction nodes that were nearest neighbors of removed species (i.e., generic or non-integrated reactions) were also removed. The final EcoCyc-derived metabolic network included 821 species participating in 1289 reactions.

When the EcoCyc network was analyzed with the PageRank algorithm in the forward direction, initial abundances were estimated by using molar ratios in minimal growth media (59% \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm PO}_4^{3 -}$$ \end{document} , 15.8% NH₃, 1.7% \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm SO}_4^{2 -}$$ \end{document} , 23.5% glucose, although similar results were obtained for glycerol or acetate carbon sources). For final abundances, cellular composition was estimated from published sources (Neidhardt and Curtiss, 1996; Weickert and Apostol, 1998; Peterson et al., 2001; Bennett et al., 2009; Milo et al., 2009), assuming that proteins, nucleic acids, and polysaccharides had been hydrolyzed into monomers (see Supplementary Information).

All network calculations were performed with the NetworkX library in Python (Hagberg et al., 2008).

4. Results

Strength distributions: The strength was calculated for each node in the two chemical networks, and rank-strength graphs are shown in Fig. 4. The x axis of each plot is the node strength, normalized so that the maximum strength of each curve is 1. The y axis is the fraction of nodes with a (normalized) strength greater than the value on the x axis, equivalent to 1 minus the cumulative distribution function (CDF). A steep slope indicates a heterogeneous distribution in which higher-ranking nodes have much higher strength than lower-ranking nodes; a shallow slope indicates a more homogeneous distribution.

FIG. 4.

Strength distributions in the metabolic (EcoCyc) and astrochemical (UDfA) reaction networks. The strength axis has been normalized to a [0,1] scale for each network to highlight the similar shapes of the strength distributions for the networks. The vertical axis shows the fraction of nodes with a (normalized) strength greater than the given value on the horizontal axis; this can be thought of as 1 minus the CDF. The EcoCyc species nodes show a region that is approximately a straight line on the log-log plot, indicative of power-law scaling over a limited range, whereas the other three sets of nodes show a shape more characteristic of exponential scaling. Non-normalized plots along with power-law and exponential fits for each curve are included in the Supplementary Information accompanying this article.

Quadrilateral motifs: Figure 2B shows the results of counting the number of directed quadrilaterals of each type in the network, normalized to the total number of edges. To determine whether these directed quadrilaterals represent significant network motifs, both graphs were shuffled in order to obtain randomized graphs. The shuffling procedure consists of randomly selecting two directed edges: (n_A,n_B) points from node A to node B and (n_C,n_D) points from C to D. As long as nodes A and C are of the same type (as are B and D), these two edges can be removed from the network and replaced by (n_A,n_D) and (n_C,n_B). This swap procedure is repeated 10n _e times, where n _e is the number of edges in the network. This procedure preserves the in and out degree of each node and leaves the strength distribution virtually unchanged while randomizing any larger-scale structure. Any measurement of structure beyond the single-edge scale is significant only if it is not preserved when the graph is shuffled.

Modularity: Community-detection calculations show a qualitative difference between the metabolic and astrochemical reaction networks. The EcoCyc network showed a clustering coefficient of 0.492 [comparable to those reported elsewhere (Newman, 2006)], which decreases to 0 upon shuffling of the network. The UDfA network, on the other hand, shows a clustering value of 0, which is unchanged upon shuffling. The EcoCyc network splits into 12 distinct subnetworks; the largest contains 334 nodes and the smallest 8.

PageRank: The top-ranking species for forward and reverse PageRank evaluations on both networks are shown in Table 1. The identities of the top-ranking species are not surprising; in each case they correspond to highly abundant species that participate in a large number of reactions and form “hubs” in the network. Closer inspection, however, shows an important difference—for the EcoCyc network, the rankings for the forward and reverse algorithms are nearly identical, whereas for UDfA they show larger differences. This intuition can be made more systematic by plotting the forward versus reverse PageRank of each species (Fig. 5). EcoCyc shows a few conjugate pairs of currency metabolites that are located far from the diagonal—ATP/ADP, NAD⁺/NADH, NADP⁺/NADPH—but most other species are clustered close to the diagonal (and the origin). UDfA shows several species located far from the diagonal—He⁺, He, and e⁻ are located above the diagonal and are therefore important reactants; CO, H, H₂, O, and OH are located below the diagonal and are stable leaving groups.

FIG. 5.

PageRank (PR) calculations for chemical networks. (A) Histogram of differences between forward and reverse PR for each node in the EcoCyc and UDfA networks. The distribution is significantly broader for the UDfA network, reflecting a partition of the network into nodes that function primarily as products and nodes that function primarily as reactants. This partitioning is less clear in the EcoCyc network, where the differences between forward and reverse PR are more tightly clustered around zero. (B, C) Scatter plot of forward vs. reverse PR for both chemical networks. Equality of the two PageRank values is indicated with a solid line. Species that are more important as reactants than as products lie above this line; important products lie below it.

Table 1.

Top-Ranking Species for Forward and Reverse PageRank for EcoCyc and UDfA Networks

EcoCyc forward	EcoCyc reverse	UDfA forward	UDfA reverse
ATP	ATP	H	He⁺
ADP	L-glutamate	HeH⁺	He
D-glucono-1,5-lactone	2-oxoglutarate	e⁻	e⁻
Diphosphate	NAD⁺	OH	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm H}_3^+$$ \end{document}
Deoxyribose-1-P	ADP	CO	H
NADPH	NADPH	O	CH₄
Ubuquinol	NADH	H₂	\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm H}_2^+$$ \end{document}

A node's forward PageRank tends to correlate strongly with its inbound degree and its reverse PageRank with its outbound degree. In fact, plots similar to Fig. 5 can be made comparing the in and out degree of nodes in the UDfA and EcoCyc networks rather than the forward and reverse PageRank values. The results look qualitatively similar (see Supplementary Information), except that construction of a histogram similar to Fig. 5A is complicated by the fact that the in and out degree, unlike PageRank, must take on integer values. One striking difference is the location of He in the UDfA plot; with an in degree much higher than its out degree, He would be located well below the diagonal on an in versus out degree plot but is located well above the diagonal in Fig. 5.

The results described above were initially obtained for the PageRank algorithm in the absence of a non-uniform personalization vector. If a chemically realistic personalization vector is provided as described above, the first major change that is observed is that a strong correlation emerges between initial (final) abundance and forward (reverse) PageRank. In this situation, more information about the network structure can be derived from the “excess PageRank,” in which the contribution to the PageRank that is attributable simply to initial (or final) abundance is removed by subtracting out a linear fit. Once this correction is applied, the distribution of PageRank asymmetry described above is unchanged. There is a positive correlation between PageRank values obtained for individual nodes with and without the personalization vectors; this correlation is fairly weak (R ²∼0.59) for EcoCyc and stronger (R ²∼0.79) for UDfA.

5. Discussion

Strength distributions: Rank-strength graphs for the reaction and species nodes of each chemical network are shown in Fig. 4. Each curve was fitted to both a power-law and an exponential decay curve; the EcoCyc species were more consistent with a power-law distribution with γ=2.06 (p _power-law=0.00015 vs. p _exponential=0.0030), the UDfA species fit somewhere between the two distribution types (p _power-law=0.0013 vs. p _exponential=0.0018), and the reaction nodes for both systems clearly fit an exponential distribution more closely (p=0.00094 vs. 0.00010 for EcoCyc; p=0.0002 vs. 4.4×10⁻⁵ for UDfA). Note that even the EcoCyc species do not follow a straight line on the log-log plot for the entire distribution; the power law holds up well at an intermediate range but shows deviations at the high and low extremes. Note also that, although the power-law scaling e ^−γ can take on any real value, strength (and degree) take on discrete integer values; this makes attempts to verify power-law scaling with statistical tests based on continuous distributions (such as a Kolmogorov-Smirnov goodness-of-fit test) inherently problematic. Most studies reporting power-law scaling in networks rely on less-sensitive measures such as the linear regression fits described above.

The scaling observed for species strength is consistent with the published results described in the Theory section above for degree distributions, which suggests that strength distributions can be used to measure features of a bipartite graph similar to the unipartite graph features measured by degree distributions. This result also underscores the importance of using the full bipartite graph when dealing with chemical systems; a projection scheme would have obscured the fact that reactions and species in EcoCyc show qualitatively different strength distributions and possibly imposed a power-law scaling upon a projected reaction subgraph (Montañez et al., 2010).

Quadrilateral motifs: Figure 2B shows clearly that quadrilateral Types 0 and 4 are clear motifs in both the EcoCyc and the UDfA networks. Type 0 corresponds to a catalytic cycle (Fig. 6); Type 4 corresponds to two parallel reactions that share a common substrate and a common product. Type 0 is a much stronger motif in EcoCyc than in UDfA; this may result from the fact that biological systems have means for driving unfavorable reactions at room temperature, whereas the cold, near-vacuum conditions of the ISM require that all non-photochemical reactions be energetically favorable and barrierless, making the formation of closed loops much more difficult. In addition, quadrilateral Types 5 and 6 are antimotifs for both EcoCyc and UDfA, meaning that they are less prevalent in the real chemical networks than in their randomized counterparts. The Type 5 antimotif can be seen as a two-step linear pathway that is folded back on itself through re-use of a by-product. For example, the metabolic example in Fig. 6 is a two-step reduction of methylglyoxal to L-1,2-propanediol via acetol. Both reductions use NADH as a reductant, and the “folding back” of the pathway occurs because of the re-use of NAD⁺ as a by-product. Likewise, the astrochemical example involves the removal of two hydrogen atoms from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm SiH}_2^+$$ \end{document} , first by reaction with a free electron and second by reaction with atomic oxygen. One reason that this antimotif may be rare in real chemical networks is that it requires a starting compound (e.g., methylglyoxal or \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\rm SiH}_2^+$$ \end{document} ) that can undergo two similar reactions in succession, so that the two reactions share a product (e.g., NAD⁺ or H). The Type 6 antimotif can be thought of as two reactions branching from a single substrate (ATP or NS in Fig. 6), where a product from one of these reactions (5-phospho-alpha-D-ribose-1-bisphoshate or H, respectively) is re-used as a reactant in the other reaction. In both of the cases presented in Fig. 6, this branching substrate needs to be able to participate in two qualitatively different types of chemistry—ATP is involved in both phosphate transfer to yield ADP and bond formation to its adenine base to yield phosphoribosyl-ATP. Although NS is much simpler, it is also involved in two different types of chemistry—electron transfer to H⁺ and replacement of S by H. Such re-use of species in branching (Type 6) or linear (Type 5) pathways appears to be unusual in real chemical networks, which leads these quadrilateral types to appear as network antimotifs.

FIG. 6.

Specific examples of quadrilaterals from the EcoCyc and UDfA networks. The shading scheme is the same as in Fig. 1; in each case the two central species nodes and the reaction nodes attached to them form a directed bipartite quadrilateral, whereas the four species nodes at the edges are included for convenience in understanding the full reactions. The Type 0 motif is a catalytic cycle in which the two central species are constantly recycled, and the Type 4 motif is the result of a single substrate-product pair (ATP/ADP or SO/O) involved in driving multiple reactions. The Type 5 and 6 antimotifs, in contrast, require the species in the quadrilateral to participate in more than one type of reaction and are rare in real chemical networks.

Modularity: As described above, community-detection calculations revealed a significantly modular structure for the biological EcoCyc network and a pronounced lack of modularity for the astrochemical UDfA network. One way to gain insight about the nature of the subnetworks returned during community detection is to calculate the betweenness centrality for each node—the fraction of shortest paths between pairs of nodes that passes through a given node. Most of the subnetworks, particularly the larger ones, contain a single highly connected metabolite (P_i, diphosphate, ammonia, etc.) or a conjugate pair of metabolites (ADP/ATP, NAD⁺/NADH, NADP⁺/NADPH, SAM/SAH) with extremely high betweenness centrality values; often more than 0.4. The modularity of the bipartite EcoCyc network is not, however, solely a result of the highly connected tail of the strength distribution—the fact that shuffling abolishes the network's modularity indicates that this feature results from larger-scale structure. It is probably significant that many of the highly connected subnetworks involve conjugate pairs of currency metabolites; ATP and ADP (for example) will still have high degree (and therefore high strength) in a randomly shuffled network but are far less likely to be connected to the same reaction (as reactant and product) than in the native network. This diminishes their ability to function as the nucleus of a highly connected clique, effectively eliminating the network's modularity. If the primary currency metabolites ATP, ADP, NAD⁺, NADH, NADP⁺, NADPH, and P_i were removed from the network, the community detection algorithm returned a modularity of 0, consistent with a dominant role for these metabolites in creating a modular network.

PageRank: The high reverse PageRank of He in the UDfA network (Fig. 5, Table 1) deserves some additional comment—one would assume that such a chemically inert species would be unlikely to function as an important reactant. He is so stable, in fact, that He⁺ is a very powerful oxidant that will take an electron away from a variety of different species in order to fill up the unoccupied valence and form He. Although He⁺ is a highly reactive species, UDfA only contains one reaction that forms it—the cosmic ray ionization of He. He⁺ receives a high reverse PageRank because of the large number of reactions in which it participates, whereas He receives a high reverse PageRank because of its exclusive connection to a single reaction that produces a highly ranked species, which gives it a crucial role in connecting the overall network.

Taking the difference between the forward and reverse PageRanks for each species gives a measure of the degree to which a given species is more important as a reactant than as a product (Fig. 5). This shows a particularly striking difference between the EcoCyc and UDfA networks; this difference measure is tightly distributed around zero for the EcoCyc network (σ=0.00122), whereas UDfA gives a significantly broader distribution (σ=0.00484; 2-sample Kolmogorov-Smirnov p=2.5×10⁻⁹). This is suggestive of the approach to equilibrium favored in astrochemical systems, which would lead to a clearer partition of reactant and product species as opposed to the far-from-equilibrium nature of biochemical systems in which most metabolites can function in both catabolic and anabolic pathways. Removing some of the distinctive features of metabolic networks increases σ only slightly—to 0.00189 if reversible reactions are made irreversible and 0.00188 if all currency metabolites are removed. This suggests that the differences in PageRank asymmetry between the two networks reflect more subtle features of the network structure.

Although these results are most readily interpreted in terms of kinetics and thermodynamics, it is important to bear in mind that the graph-based calculations do not explicitly include this type of information. Instead, kinetics enters indirectly through the graph topology, because the only reactions that are included in the network are those that were considered kinetically feasible, either based on general physical-chemical principles (UDfA) or because an enzyme exists that is able to catalyze them (EcoCyc), possibly by coupling unfavorable reactions to favorable processes such as ATP hydrolysis.

6. Conclusions

Despite the tremendous differences in their physical environments in terms of density, temperature, solvation, and energy sources, the reaction networks present in a bacterial cell and in dark molecular clouds show a number of striking similarities. Both systems show reaction node strength distributions that follow exponential scaling, which suggests topological similarities in the large-scale ordering of the networks. The same two strong directed quadrilateral motifs can be found in both networks, which suggests that the same fundamental principles of chemical reactivity shape the local connectivity structure of both networks. Examination of the PageRank asymmetries depicted in Fig. 5 shows that both systems contain species that are more important as reactants than as products: He, He⁺, and e⁻ in UDfA, and ATP and NADPH in EcoCyc. Both networks also contain stable leaving groups that are more important as products than as reactants: CO, H, and H₂ in UDfA, and ADP, P_i, CO₂, and NADP⁺ in EcoCyc.

The most interesting features of these two chemical networks, however, are probably their differences. The ultimate explanation for the power-law scaling of EcoCyc's species degree distribution is beyond the scope of this study, and extensive discussions of power-law behavior have been given elsewhere (Newman, 2005). Whatever the reason for the power-law scaling, it remains a telling difference between biological and nonbiological chemical networks. Although the two quadrilateral network motifs revealed by Fig. 2B are common between EcoCyc and UDfA, it should be noted that Type 0, which corresponds to catalytic cycles, plays a larger role in EcoCyc, whereas UDfA is more completely dominated by the parallel pathways represented in Type 4 motifs. Another interesting difference is that, while Types 1, 2, and 3 in UDfA are nearly indifferent to shuffling, which indicates that they are present in roughly the same levels as would be expected in a similar random network, every type of directed quadrilateral is either a motif or an antimotif in EcoCyc, with Types 5 and 6 representing stronger antimotifs than they do in UDfA. This suggests that, although the fundamental physics of complex chemical networks dictate a commanding role for Types 0 and 4, evolutionary pressures led to a less-random character for the distribution of the other quadrilateral types in EcoCyc. Compared to UDfA, EcoCyc also shows more balance between the forward and reverse PageRank values for most nodes (Fig. 5). This could be a result of evolutionary pressure to maintain a highly dynamic system in which metabolic fluxes can be reconfigured by enzymatic regulation and any metabolite can participate in a variety of different pathways, depending on the changing needs of the cell. A system in which most metabolites can clearly be classified as preferentially reactants or products will tend to be more static and to progress toward equilibrium, as is generally the case for nonbiological reaction systems.

The presence of significant network motifs and antimotifs in both the astrochemical and metabolic networks emphasizes the fact that chemical reaction networks are decidedly nonrandom. The randomized chemical networks used to assess the significance of the quadrilateral network motifs exhibit the same degree and strength distributions as the real networks but are dramatically different at larger scales of organization. This has clear implications for network-based models for the origin of life (Kauffman, 1986), which often rely on the presence of phase transitions in random graph models. It is well known that Erdős–Rényi random graphs are not a good approximation of metabolic reaction networks (Jeong et al., 2000) and that models with a power-law distribution are likely to provide a better approximation. The results in this study make it clear that constructing a plausible chemical network will likely require more than simply reproducing the correct degree (or strength) distribution. Chemically realistic random networks will, for example, need to exhibit the types of quadrilateral motif distributions shown in Fig. 2. Networks based on “artificial chemistry,” where in silico “molecules” are related to each other by chemistry-like rules, may be a valuable approach for the generation of hypothetical chemical networks with plausible structure (Benkö et al., 2003; Hintze and Adami, 2008).

It should be emphasized that neither of the chemical networks examined here provides a complete enough description of the physical environment to fully explain why these systems function the way they do in nature. In the case of astrochemical networks, the timescale for chemical equilibration is uncomfortably close to the timescale for gravitational collapse. This is what is meant when simulations based on chemical networks are referred to as “pseudo-time dependent”—true time dependence would require coupling between the physical and chemical evolution of the cloud (Lee et al., 1996). The biological case is even more complicated. Biological reactions occur in the crowded confines of the cell, where local concentrations can be high, often permitting enzymes to operate in both directions under different intracellular conditions and leading to the large number of reversible reactions observed in biochemical networks. An even greater source of additional complexity comes from regulatory mechanisms, by which the cell is able to control both the quantity and activity of catalytic enzymes in order to meet changing physiological needs. The binary nature of the networks examined in this study, where reactions are either present or absent, does not capture this level of regulatory nuance.

It should also be noted that there are many more differences between astrochemical and metabolic reaction networks than the mere presence or absence of biology, and it would be unwise to attribute every observed difference to biological organization. For this reason, extension of this type of analysis to other types of reaction systems could yield further insight into the constraints that shape the structure of chemical networks, both biological and nonbiological. For example, a network analysis of an abiotic, non-equilibrium reaction system such as the Belousov-Zhabotinsky reaction (Field et al., 1972; Mikhailov and Showalter, 2006) could help to determine whether some of the observed features of the biological system are a result of non-equilibrium dynamics or whether they are primarily a response to evolutionary pressures. Planetary atmospheres could provide examples of abiotic systems under somewhat more clement conditions than the ISM, and geochemical systems could provide a case of a system that is not biological in the sense of being acted upon directly by evolution but may be driven by the action of biological systems. Further insight into the topology of bipartite metabolic reaction networks could also be obtained by applying these tools to more sophisticated stoichiometric metabolic models (Feist et al., 2007) or to multiple-organism networks that could allow for the modeling of metabolism at the ecosystem level (Caspi et al., 2010). Such metabolic networks are unlikely to be qualitatively different from EcoCyc; however, they could provide a sense of the robustness of metabolic network features to the inclusion or exclusion of reactions that differ between differently curated databases. The tools of network theory are well suited to comparative studies of diverse types of chemical systems and can provide us with valuable insights into where life fits in the broader landscape of chemical complexity.

Footnotes

Acknowledgments

This research was supported by a grant from the NASA Astrobiology Institute (NNA08CN85A). The authors would also like to thank Aurélien Mazurie and Sarah Hörst for helpful discussions and references to pertinent literature and Aric Hagberg for assistance with the PageRank algorithm.

Author Disclosure Statement

C.J.: No competing financial interests exist. T.D.: No competing financial interests exist.

Abbreviations

CDF, cumulative distribution function; EcoCyc, E. coli encyclopedia; ISM, interstellar medium; UDfA, UMIST Database for Astrochemistry.

References

Albert

, Barabási

A.L.

2002. Statistical mechanics of complex networks. Rev Mod Phys, 74:47–97.

Alon

2007. Network motifs: theory and experimental approaches. Nat Rev Genet, 8:450–461.

Alves

, Chaleil

R.A.

, Sternberg

M.J.

2002. Evolution of enzymes in metabolism: a network perspective. J Mol Biol, 320:751–770.

Anicich

V.G.

2003. An index of the literature for bimolecular gas phase cation-molecule reaction kinetics. JPL-Publication-03-19, Jet Propulsion Laboratory: Pasadena, CA.

Arita

2004. The metabolic world of Escherichia coli is not small. Proc Natl Acad Sci USA, 101:1543–1547.

Barabási

A.L.

2009. Scale-free networks: a decade and beyond. Science, 325:412–413.

Barabási

A.L.

, Oltvai

Z.N.

2004. Network biology: understanding the cell's functional organization. Nat Rev Genet, 5:101–113.

Benkö

, Flamm

, Stadler

P.F.

2003. A graph-based toy model of chemistry. J Chem Inf Comput Sci, 43:1085–1093.

Bennett

B.D.

, Kimball

E.H.

, Gao

, Osterhout

, Van Dien

S.J.

, Rabinowitz

J.D.

2009. Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nat Chem Biol, 5:593–599.

10.

Carlson

J.M.

, Doyle

1999. Highly optimized tolerance: a mechanism for power laws in designed systems. Phys Rev E Stat Nonlin Soft Matter Phys, 60:1412–1427.

11.

Caspi

, Altman

, Dale

J.M.

, Dreher

, Fulcher

C.A.

, Gilham

, Kaipa

, Karthikeyan

A.S.

, Kothari

, Krummenacker

, Latendresse

, Mueller

L.A.

, Paley

, Popescu

, Pujar

, Shearer

A.G.

, Zhang

, Karp

P.D.

2010. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. Nucleic Acids Res, 38:D473–D479.

12.

Dorogovtsev

S.N.

, Goltsev

A.V.

, Mendes

J.F.F.

2002. Pseudofractal scale-free web. Phys Rev E Stat Nonlin Soft Matter Phys, 65. 10.1103/PhysRevE.65.066122.

13.

Feist

A.M.

, Henry

C.S.

, Reed

J.L.

, Krummenacker

, Joyce

A.R.

, Karp

P.D.

, Broadbelt

L.J.

, Hatzimanikatis

, Palsson

B.O.

2007. A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information. Mol Syst Biol, 3. 10.1038/msb4100155.

14.

Field

R.J.

, Noyes

R.M.

, Körös

1972. Oscillations in chemical systems. II. Thorough analysis of temporal oscillation in bromate-cerium-malonic acid system. J Am Chem Soc, 94:8649–8664.

15.

Ganong

B.R.

, Raetz

C.R.

1982. Massive accumulation of phosphatidic acid in conditionally lethal CDP-diglyceride synthetase mutants and cytidine auxotrophs of Escherichia coli. J Biol Chem, 257:389–394.

16.

Hagberg

A.A.

, Schult

D.A.

, Swart

P.J.

2008. Exploring network structure, dynamics, and function using NetworkX. Proceedings of the 7^th Python in Science Conference (SciPy2008) Varoquaux

, Vaught

, Millmans

SciPy Community: Pasadena, CA, 11–15.

17.

Herbst

, van Dishoeck

E.F.

2009. Complex organic interstellar molecules. Annu Rev Astron Astrophys, 47:427–480.

18.

Hintze

, Adami

2008. Evolution of complex modular biological networks. PLoS Comput Biol, 4:e23. 10.1371/journal.pcbi.0040023.

19.

Jeong

, Tombor

, Albert

, Oltvai

Z.N.

, Barabasi

A.L.

2000. The large-scale organization of metabolic networks. Nature, 407:651–654.

20.

Jolley

C.C.

, Douglas

2010. A network-theoretical approach to understanding interstellar chemistry. Astrophys J, 722:1921–1931.

21.

Kauffman

S.A.

1986. Autocatalytic sets of proteins. J Theor Biol, 119:1–24.

22.

Keseler

I.M.

, Bonavides-Martínez

, Collado-Vides

, Gama-Castro

, Gunsalus

R.P.

, Johnson

D.A.

, Krummenacker

, Nolan

L.M.

, Paley

, Paulsen

I.T.

, Peralta-Gil

, Santos-Zavaleta

, Shearer

A.G.

, Karp

D.P.

2009. EcoCyc: a comprehensive view of Escherichia coli biology. Nucleic Acids Res, 37:D464–D470.

23.

Kreimer

, Borenstein

, Gophna

, Ruppin

2008. The evolution of modularity in bacterial metabolic networks. Proc Natl Acad Sci USA, 105:6976–6981.

24.

Langville

A.N.

, Meyer

C.D.

2005. A survey of eigenvector methods for Web information retrieval. SIAM Rev Soc Ind Appl Math, 47:135–161.

25.

Lee

H.H.

, Bettens

R.P.A.

, Herbst

1996. Fractional abundances of molecules in dense interstellar clouds: a compendium of recent model results. Astron Astrophys Suppl Ser, 119:111–114.

26.

Light

, Kraulis

2004. Network analysis of metabolic enzyme evolution in Escherichia coli. BMC Bioinformatics, 5. 10.1186/1471-2105-5-15.

27.

Light

, Kraulis

, Elofsson

2005. Preferential attachment in the evolution of metabolic networks. BMC Genomics, 6. 10.1186/1471-2164-6-159.

28.

H.W.

, Zeng

A.P.

2003. The connectivity structure, giant strong component and centrality of metabolic networks. Bioinformatics, 19:1423–1430.

29.

Mikhailov

A.S.

, Showalter

2006. Control of waves, patterns and turbulence in chemical systems. Phys Rep, 425:79–194.

30.

Millar

T.J.

2004. Organic molecules in the interstellar medium. Astrobiology: Future Perspectives. Ehrenfreund

, Irvine

, Owen

, Becker

et al. Kluwer Academic Publishers: Dordrecht, the Netherlands, 17–31.

31.

Milo

, Jorgensen

, Moran

, Weber

, Springer

2009. BioNumbers—the database of key numbers in molecular and cell biology. Nucleic Acids Res, 38:D750–D753.

32.

Montañez

, Medina

M.A.

, Solé

R.V.

, Rodríguez-Caso

2010. When metabolism meets topology: reconciling metabolite and reaction networks. Bioessays, 32:246–256.

33.

Neidhardt

F.C.

, Curtiss

1996. Escherichia coli and Salmonella: Cellular and Molecular Biology. ASM Press: Washington DC.

34.

Newman

M.E.J.

2003. The structure and function of complex networks. SIAM Rev Soc Ind Appl Math, 45:167–256.

35.

Newman

M.E.J.

2005. Power laws, Pareto distributions and Zipf's law. Contemporary Physics, 46:323–351.

36.

Newman

M.E.J.

2006. Modularity and community structure in networks. Proc Natl Acad Sci USA, 103:8577–8582.

37.

Page

, Brin

, Motwani

, Winograd

1999. The PageRank Citation Ranking: Bringing Order to the Web. Technical Report. Stanford InfoLab: Stanford, CA.

38.

Peterson

J.D.

, Umayam

L.A.

, Dickinson

, Hickey

E.K.

, White

2001. The Comprehensive Microbial Resource. Nucleic Acids Res, 29:123–125.

39.

Solé

R.V.

, Munteanu

2004. The large-scale organization of chemical reaction networks in astrophysics. Europhys Lett, 68:170–176.

40.

Spirin

, Gelfand

M.S.

, Mironov

A.A.

, Mirny

L.A.

2006. A metabolic network in the evolutionary context: multiscale structure and modularity. Proc Natl Acad Sci USA, 103:8774–8779.

41.

Tanaka

2005. Scale-rich metabolic networks. Phys Rev Lett, 94. 10.1103/PhysRevLett.94.168101.

42.

Wagner

, Fell

D.A.

2001. The small world inside large metabolic networks. Proc Biol Sci, 268:1803–1810.

43.

Wakelam

, Smith

, Herbst

, Troe

, Geppert

, Linnartz

, Öberg

, Roueff

, Agúndez

, Pernot

, Cuppen

H.M.

, Loison

J.C.

, Talbi

2010. Reaction networks for interstellar chemical modelling: improvements and challenges. Space Sci Rev, 156:13–72.

44.

Weickert

M.J.

, Apostol

1998. High-fidelity translation of recombinant human hemoglobin in Escherichia coli. Appl Environ Microbiol, 64:1589–1593.

45.

Wilson

E.H.

, Atreya

S.K.

2004. Current state of modeling the photochemistry of Titan's mutually dependent atmosphere and ionosphere. J Geophys Res, 109. 10.1029/2003JE002181.

46.

Woodall

, Agúndez

, Markwick-Kemper

A.J.

, Millar

T.J.

2007. The UMIST database for astrochemistry 2006. Astron Astrophys, 466:1197–1204.

47.

Zhao

, Ding

G.H.

, Tao

, Yu

Z.H.

, Luo

J.H.

, Cao

Z.W.

, Li

Y.X.

2007. Modular co-evolution of metabolic networks. BMC Bioinformatics, 8. 10.1186/1471-2105-8-311.

48.

Zhou

, Nakhleh

2011. Properties of metabolic graphs: biological organization or representation artifacts? BMC Bioinformatics, 12. 10.1186/1471-2105-12-132.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.03 MB