Genotype networks, innovation, and robustness in sulfur metabolism

Abstract

Background

A metabolism is a complex network of chemical reactions. This network synthesizes multiple small precursor molecules of biomass from chemicals that occur in the environment. The metabolic network of any one organism is encoded by a metabolic genotype, defined as the set of enzyme-coding genes whose products catalyze the network's reactions. Each metabolic genotype has a metabolic phenotype. We define this metabolic phenotype as the spectrum of different sources of a chemical element that a metabolism can use to synthesize biomass. We here focus on the element sulfur. We study properties of the space of all possible metabolic genotypes in sulfur metabolism by analyzing random metabolic genotypes that are viable on different numbers of sulfur sources.

Results

We show that metabolic genotypes with the same phenotype form large connected genotype networks - networks of metabolic networks - that extend far through metabolic genotype space. How far they reach through this space depends linearly on the number of super-essential reactions. A super-essential reaction is an essential reaction that occurs in all networks viable in a given environment. Metabolic networks can differ in how robust their phenotype is to the removal of individual reactions. We find that this robustness depends on metabolic network size, and on other variables, such as the size of minimal metabolic networks whose reactions are all essential in a specific environment. We show that different neighborhoods of any genotype network harbor very different novel phenotypes, metabolic innovations that can sustain life on novel sulfur sources. We also analyze the ability of evolving populations of metabolic networks to explore novel metabolic phenotypes. This ability is facilitated by the existence of genotype networks, because different neighborhoods of these networks contain very different novel phenotypes.

Conclusions

We show that the space of metabolic genotypes involved in sulfur metabolism is organized similarly to that of carbon metabolism. We demonstrate that the maximum genotype distance and robustness of metabolic networks can be explained by the number of superessential reactions and by the sizes of minimal metabolic networks viable in an environment. In contrast to the genotype space of macromolecules, where phenotypic robustness may facilitate phenotypic innovation, we show that here the ability to access novel phenotypes does not monotonically increase with robustness.

Background

In any biological system, genotypes contain the information needed to make phenotypes. The relationship between genotype and phenotype is also known as a genotype-phenotype map [1]. The ability to analyze different kinds of biological systems computationally has allowed a detailed characterization of genotype-phenotype maps for different systems. One common feature of genotype-phenotype maps is the existence of genotype networks, connected sets of genotypes that adopt the same phenotype. They exist in systems as different as model proteins [2], RNA secondary structures [3], regulatory circuits [4], and metabolic networks [5, 6]. Another feature is the large phenotypic diversity that is found in different neighborhoods of a genotype network [3–6]. These two properties facilitate the exploration of novel and potentially beneficial phenotypes in genotype space. By analyzing genotype-phenotype maps of different systems, one can identify general features of genotype maps, as well as features that are specific to a system.

In this work we concentrate on the genotype-phenotype maps of metabolic networks involved in the utilization of sulfur. We have two aims. Aim 1 is to examine how general earlier observations about the genotype-phenotype map of carbon metabolism are [5, 6]. We do so by examining if these observations also apply to sulfur metabolism. In particular, we investigate the existence of genotype networks whose members share the same phenotype, and the amount of phenotypic diversity in their neighborhoods. Aim 2 is to study how rapidly evolving populations of networks "discover" metabolic innovations in metabolic genotype space. Specifically, we are interested in how the rate of discovery depends on the robustness of a metabolic system. This robustness indicates a metabolic network's ability to preserve its biosynthetic capacity upon random removal of reactions. Previous work on macromolecules showed that the robustness of a molecule's phenotype to mutations can accelerate the rate at which evolving populations discover new phenotypes [7, 8]. We will ask whether the same holds for metabolic systems.

Carbon metabolism comprises so many reactions that the computational demands of studying population processes in its genotype space are too high for current computational technology. Sulfur metabolism, in contrast, comprises a smaller number of chemical reactions, which renders the computational analysis of population processes more tractable. Despite being involved in fewer reactions, sulfur is no less essential to biological organisms than other elements, such as carbon or nitrogen. Sulfur is a versatile and integral element in the biochemistry of organisms [9, 10]. Its presence in biological organisms ranges from 0.5% to 50% of dry weight [9]. It occurs in multiple oxidation states, ranging from the highly oxidized S⁴⁺ to the reduced state S^2-. This versatility in oxidation state may explain the diversity of sulfur metabolism and why it is involved in both anabolism as well as catabolism. In catabolism, depending on the environment, sulfur can be used as an electron acceptor or an electron donor, and in some cases even both as donor and acceptor. In anabolism, sulfur must first be reduced in a sequence of energetically expensive steps before being incorporated into biomass [9].

Sulfur is present in two major constituents of biomass, the amino-acid cysteine, which confers stability to proteins through disulfide bonds, and the amino acid methionine, which is the first amino acid of many proteins. Sulfur is also a part of S-adenosylmethionine (also known as AdoMet or SAM). This compound is a cysteine metabolite that is a major methyl donor to the methyl carrier metabolite tetrahydrofolate, which is indispensable for amino acid synthesis, and for the methylation of biomolecules. Furthermore, sulfur is the active element in coenzyme-A, an acyl carrier metabolite involved in the calvin cycle and in lipid synthesis. Sulfur is also present in the active core of iron-sulfur proteins, which are involved in a number of important reactions. Examples include nitrogenase, which enables the fixation of nitrogen, and hemoglobin, which enables the transport of oxygen. Another prominent molecule involving sulfur is glutathione, a peptide responsible for protection against oxidative stress in cells.

We next outline the order of our analyses in the Results section. First we introduce two concepts that allow us to estimate some properties of genotype space organization. The first is that of a minimal metabolic network. This is a metabolic network from which no reactions can be removed without destroying its viability in a given environment, that is, its ability to synthesize all essential biomass molecules. The second concept is that of a superessential reaction. For our purpose, a superessential reaction is an essential chemical reaction that occurs in all minimal networks. After these preliminary analyses, we demonstrate the existence of long phenotype-preserving paths through metabolic genotypes space that allow exploring this space through many single phenotype-preserving mutations (aim 1). The maximum length of these paths and metabolic network size can be estimated and varies linearly with the number of superessential reactions. We show that the robustness of metabolic networks depends both on their size and on the average size of minimal metabolic networks viable on a given number of sulfur sources. Next, we show that the existence of neutral paths allows evolving metabolic networks to encounter an increasing number of novel phenotypes (aim 2). We finally explore the relationship between robustness and a population's ability to access novel phenotypes through changes in a network's reactions. In contrast to macromolecules, where robustness may facilitate phenotypic innovation [7, 8], we find that the ability to find novel phenotypes in our system peaks at intermediate robustness.

Results

The model

We follow an approach taken in a previous study of large-scale metabolic networks [5]. We define a metabolic genotype as the set of biochemical reactions that may take place in an organism, and that are catalyzed by gene-encoded enzymes. The set of all reactions used in this work is a subset of 1221 reactions out of 5871 reactions we curated previously [5] from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [11]. These reactions comprise all elementally-balanced reactions that involve sulfur containing metabolites (see methods for details).

A metabolic genotype can be represented in at least 2 different ways (Figure 1). The first views it as a metabolic network graph whose nodes are metabolites. Reactions are represented as directed links from substrate metabolites to product metabolites (Figure 1A). The second views it as a list of reactions (Figure 1C), or, equivalently a binary vector whose length - 1221 reactions in our case - corresponds to the number of reactions in a known reaction "universe". Each position i in this vector corresponds to a reaction. Its values ('0') or ('1') at position i indicate the inability or ability of the organism to catalyze the corresponding reaction (Figure 1C).

Figure 1

Genotype-phenotype map of metabolic networks. Different representations of a hypothetical metabolic network (A), as a node in a genotype network (B), or as a binary vector (C) listing the reactions in the network. Each genotype (circles) on the genotype network in (B) has 1221 neighbors (not all edges are drawn) that differ by a single mutation. Neighbors in (B) are connected by edges. The colors of the genotypes represent different phenotypes. The phenotypes of the metabolic networks are computed using flux balance analysis applied to 124 environments with different sulfur sources. Two hypothetical phenotypes are represented in (D) as binary vectors listing the environments a genotype is viable in (D). Random evolutionary walks can be seen as paths on a genotype network that stay on genotypes with the same phenotype (represented as the genotype color). "Mutations" correspond to additions or deletions of individual reactions from the metabolic network. The number of genotypes in the genotype space is 2¹²²¹.

We define the phenotype of a metabolic network as the subset of sulfur sources (out of 124 possible sources we consider, see Methods) that allow the network (metabolic genotype) to synthesize all biomass components, if one of the sulfur sources is provided as the sole sulfur source to the organism. We represent this phenotype as a binary vector of length 124 whose entry at position i indicates viability if sulfur source i is the sole sulfur source (Figure 1D). This is not the only way to define a metabolic phenotype, but it is appropriate for our purpose. An obvious alternative phenotype definition would count the number of biomass metabolites that a network can produce in a given environment. However, because all these metabolites are essential for survival of an organism, networks that cannot synthesize some of them are of limited biological relevance. Additional advantages of the phenotype definition we chose are that it allows a straightforward and systematic comparison of phenotypes, and enables us to study metabolic innovation in a biologically sensible way. Using this phenotype definition, a metabolic innovation is the ability to synthesize biomass metabolites from a new, previously unusable sulfur source.

To determine metabolic phenotypes from genotypes, we use flux balance analysis [12], a computational method that finds a growth-maximizing steady-state metabolic flux through all reactions in a metabolic network. This method requires information about the stoichiometry of every metabolic reaction, a maximally allowed flux of each metabolite in and out of the environment, and information about an organism's biomass composition (see Methods for details). We focus on a metabolic network's qualitative ability to produce all sulfur-containing biomass precursors. We will study networks that are able to do so from each one of a specific set of sole sulfur sources. For brevity, we call such networks viable. We will also refer to the number S of sulfur sources that a metabolic network must be viable on as the environmental demand imposed on the network.

We next introduce the concept of a genotype network for metabolic networks (Figure 1B) [5]. The nodes in this network correspond to individual genotypes (metabolic networks) with the same phenotype. Two genotypes are linked - they are neighbors - if they differ in a single reaction. A genotype network thus is a network of metabolic networks. This concept is useful when we examine the evolution of metabolic networks through the addition and elimination of metabolic reactions, which can occur, for example, by horizontal gene transfer [13, 14], or through loss-of-function mutations in enzyme-coding genes. Consider the metabolic network genotype G ₁ of some organism. This genotype is a node on the genotype network associated with this genotype's phenotype. If some variant G ₂ of this network - obtained through an addition or a deletion of a reaction - has the same phenotype as G ₁ , it will be a neighbor of G ₁ on the same genotype network. In this manner, one can envision phenotype-preserving evolutionary change of metabolic genotypes as a path through a genotype network. Such paths correspond to successive hops from genotype to genotype, by way of the edges connecting neighboring genotypes (Figure 1B). For our analysis, it will be useful to define a distance >D between two metabolic network genotypes as the fraction of reactions in which two metabolic networks differ, or

$$D=1-\frac{2R_c}{N_1+N_2},$$

where R _c is the number of reactions shared by both networks, and N ₁ and N ₂ are the total numbers of reactions in networks G ₁ and G ₂ , respectively. This formula simplifies to D = 1-R_C/N when both networks have the same size N..

Minimal viable metabolic networks can be diverse and contain many superessential reactions

We begin with an analysis of minimal viable networks, which provides insights into the reactions that are essential to utilize a specific set of sulfur sources. A minimal metabolic network is a network in which all reactions are essential and none can be removed without rendering it inviable. For any one given phenotype P, there may be multiple viable minimal networks. Random minimal networks can be generated by starting from a network comprised of all 1221 reactions - which is viable on all sulfur sources - and eliminating randomly chosen reactions one-by-one, until no reactions can be further removed without rendering the network inviable on the sulfur sources defined by P. We note that a minimal network is not the same as the network with the smallest possible number of reactions with a given phenotype, which could be very difficult to find in a vast metabolic genotype space.

We generated 1000 random minimal metabolic networks viable on a given number S of sulfur sources (see methods). Specifically, we generated 100 minimal networks for 10 random sets of sulfur sources with the same number S - but not necessarily identity - of sources. We note in passing that such networks often also happen to be viable on additional sulfur sources that we did not require them to be viable on (Additional file 1). Figure 2A shows the distribution of genotype distances for pairs of minimal metabolic networks viable on S = 1, 20, or 60 sulfur sources. The figure demonstrates that, first, random minimal metabolic networks can be very different from one another. Their genotype distance may exceed D = 0.8, meaning that they may share fewer than 20 percent of reactions. Second, their average distance depends on the number of sulfur sources a network needs to be viable on. Specifically, the average genotype distance is largest D _avg = 0.6 for minimal metabolic networks viable on S = 1 sulfur source, and decreases to D _avg = 0.3 for networks viable on S = 60 sulfur sources. Third, the distribution of genotype distances is much wider for metabolic networks subject to few environmental demands (S = 1) where it ranges from D _avg = 0.2 to D _avg = 0.8, than for metabolic networks subject to many environmental demands (S = 60) where it ranges from D _avg = 0.2 to D _avg = 0.4.

Figure 2

Distances, sizes and superessential reactions found in random minimal metabolic networks subject to different environmental demands.(A) Distribution of genotype distance between pairs of minimal metabolic networks viable under the same environmental demands. (B) Average size (closed circles) and average number of superessential reactions (open circles) of 1000 minimal metabolic networks as a function of environmental demands S on a network. The number of superessential reactions was obtained by counting the number of reactions common to 100 minimal metabolic networks generated with the same set of sulfur sources. For each data point, we used 10 different sets of sulfur sources with size S. Error bars represent the standard deviations of the distributions.

Figure 2B (filled circles) shows the average size of minimal networks N _min as a function of the number of sulfur sources they are viable on. It ranges from N _min = 14 reactions for S = 1 to N _min = 87 reactions for S = 60. By definition, all reactions in a minimal network are essential, but some of these reactions are special because they occur in all minimal networks viable on a given set of sulfur sources. We call these reactions superessential reactions [6]. The open circles in Figure 2B shows the number of superessential reactions as a function of the environmental demands S on a network. The number of superessential reactions R _SE increases with S, but it is generally much lower than the total number of reactions in a minimal metabolic network. For example, at S = 1, 4 out of 14 reactions are superessential. At S = 60, 44 out of 87 reactions are superessential. We will show that the number of superessential reactions plays an important role in one of our analyses below.

Many viable sulfur metabolic network genotypes are connected via paths that lead far through metabolic genotype space

We next extended our previous work on carbon metabolism to ask about the existence of genotype networks in the space of sulfur-involving reactions, and of neutral paths that traverse such networks while preserving a metabolic phenotype. We define a neutral path as a series of mutations (reaction additions or deletions) that leave a phenotype intact (Figure 1B). We emphasize that we do not use the term neutrality in its meaning of unchanged fitness in the field of molecular evolution [15], but merely for brevity, in the sense of preserving viability on a specific set of sulfur sources. Changes to metabolic networks such as additions or deletions of chemical reactions can potentially have a positive or negative effect on fitness. The addition of a chemical reaction may have a beneficial effect if it increases the rate at which biomass is synthesized, or it may have a deleterious effect if it generates metabolites that interfere with cell physiology. Similarly, the deletion of a reaction may have both beneficial and deleterious fitness effects [16]. Studies on compensatory mutations in macromolecules, mutations that compensate for the fitness effects of previous mutations with negative effects, show that fitness neutrality is not a prerequisite for a population's genetic change on a genotype network [17, 18].

We are especially interested in two questions. How far does a neutral path typically lead through genotype space? And how does this distance depend on the number N of reactions in a network, and on the environmental demands on the network? To answer these questions, we performed 200 random walks of 10 000 mutations each for metabolic networks of various sizes, and for various environmental demands. Specifically, for networks of each size we performed 20 random walks for each of 10 different sets of S of sole sulfur sources that we required the network to be viable on. Each random walk started from a random initial viable metabolic network comprising N reactions (see methods for details). We allowed N to vary by no more than one reaction during the random walk. Moreover, each step in the random walk had to preserve viability. Finally, none of the steps was allowed to decrease the distance to the starting network, in order to maximize the distance from this network (see methods for details).

Figure 3A shows the maximum genotype distance D _max obtained in such random walks for networks up to 300 reactions, where we required viability on S= 1, 5, 10, 20, 40, or 60 different sole sulfur sources. This distance is in general large. For example, D _max is greater than 0.7 for all metabolic networks with more than 200 reactions. For each value of S, the data point at the smallest value of N (horizontal axis) corresponds to the minimal metabolic networks we discussed earlier. Perhaps surprisingly, these minimal networks can not only be very diverse, as we saw earlier, but neutral paths starting from any one such network can also reach far through genotype space. For example, the maximal length of neutral paths is D _max = 0.65 for minimal metabolic networks viable on S = 1 sulfur source, and still a sizeable D _max = 0.38 for metabolic networks viable on S = 60 sulfur sources. To provide a point of reference, the E. coli metabolic network has 142 reactions involving sulfur. Random viable metabolic networks of this size would have maximum genotype distances between D _max = 0.96 (for S = 1) and D _max = 0.60 (for S = 60).

Figure 3

Maximum genotype distances and robustness depend strongly on the number of superessential reactions and the sizes of minimal metabolic networks. (A) Average maximum genotype distance for metabolic networks of different sizes and subject to different environmental demands after random walks of 10 000 accepted reaction changes. (B) Robustness of random metabolic networks with different sizes and subject to different environmental demands. Each data point is an average over 200 random walks (20 random walks for 10 different sets of environmental demands with the same number S of sulfur sources).

Maximal genotype distance and robustness of metabolic networks are well approximated by simple properties of minimal networks

We asked whether the maximal genotype distance of networks of a given size, as well as their robustness to reaction removal can be predicted from properties of the underlying minimal networks. The answer is yes.

Figure 3A shows that the maximal possible genotype distance D _max between metabolic networks of the same size increases with metabolic network size N. The solid lines show the relationship between the maximal genotype distance D _max and metabolic network size N as predicted by the equation

$$\rho =1+m-\frac{N_{\min }}{N}(1+m),$$

(1)

Here, R _SE is the number of reactions that are super-essential for a given environmental demand S. We had estimated this number in our previous analysis of minimal networks (Figure 2B). The simple relationship of equation (1) fits our numerical data (Figure 3A) remarkably well and corresponds exactly to our distance function when the number of superessential reactions R _SE replaces the number of common reactions R _C between networks of the same size. This implies that the only common reactions between maximally distant networks are the superessential reactions. Therefore, as the size of a network increases, more phenotype-preserving changes become possible. For networks of the smallest size, D _max , systematically overestimates the maximal genotype distance, but it does so by no more than 10% percent. We note that our estimates of maximum genotype distances are only lower bounds, such that this discrepancy may result from our limited ability to find maximal genotype distances accurately. In sum, a simple, linear function of the number of superessential reactions at any one environmental demand S approximates the maximal genotype distance between networks well.

Next we examined how network robustness depends on the size of metabolic networks and on environmental demands. We define robustness as the fraction of non-essential reactions in a metabolic network. Figure 3B shows the robustness of metabolic networks as a function of network size N and varying environmental demands S on a network. Large metabolic networks with 200 reactions or more have a robustness ρ > 0.6 for all values of S. For smaller metabolic networks (N < 200), robustness ranges from ρ = 0.9 under low environmental demands (S = 1) to ρ = 0.2 under high environmental demands (S = 60). The relationship between ρ and network size N can be explained by noting that ρ = 1-R_ess/N where R _ess is the number of essential reactions. We find that R _ess decreases linearly with increasing metabolic network size (Additional file 2) and is described by the function R_ess = N_min(1+m)-Nm, In this equation, N _min is the average size of minimal metabolic networks (estimated above for given S) and m is the rate at which the number of essential reactions decreases with increasing metabolic network size (estimated from data in Additional file 2). The question why the number of essential reactions R _ess decreases with increasing size has a simple answer. As one increases the size of a metabolic network by adding reactions and entire pathways to minimal metabolic networks, some reactions may become non-essential because the added reactions create alternative pathways for biomass metabolite synthesis. Describing ρ in terms of N, N _min and m, we arrive at the following relation

$$\rho =1+m-\frac{N_{\min }}{N}(1+m),$$

(2)

which is plotted as the solid lines in Figure 3B and fits the data very well.

This relationship means that network robustness is a linear function of the ratio N _min /N, whose inverse indicates how much larger a given network is than a minimal network for a given S, and of the rate at which reaction essentiality declines (robustness increases) with increasing N.

The diversity of phenotypes found in the neighborhood of two metabolic networks changes rapidly with their genotype distance

Thus far, we have concentrated on the characteristics of individual sets of genotypes viable on a given number S of sulfur sources, and on the genotype networks they form. Long paths through a genotype network can contribute to evolutionary innovation in metabolic phenotypes, if many novel phenotypes can be encountered near such a path. We next asked whether this is the case, and how this number of novel phenotypes depends on environmental demands on a network. We consider a phenotype to be novel if it confers viability in a set of new sulfur sources, in addition to those required by the environmental demands imposed on the metabolic network.

We first introduce the notion of a (1-mutant) neighborhood around a metabolic network genotype, which comprises all networks that differ from the genotype by a single reaction (Figure 1B). Because our genotype space has 1221 metabolic reactions, each metabolic network has 1221 neighbors. Of all these neighbors, some will be inviable in any given environment (these are the mutants that have lost an essential reaction), some will maintain the same phenotype, and some will have a novel phenotype while being viable in this environment. That is, they will have gained viability on a new sulfur source. We focus on the latter class of neighbors in this section.

We asked how different are the novel phenotypes in the neighborhood of two metabolic networks G and G _k on the same genotype network, where G _k is a metabolic network derived from G through k random reaction changes. That is, we determined the fraction of novel phenotypes that occurred in the neighborhood of only one but not the other network. Below we refer to it as the fraction of novel phenotypes unique to one neighborhood. If this fraction is very small even for large k, then networks in different regions of a genotype space will have mutational access to similar novel phenotypes. Figure 4 shows that the opposite is the case. We obtained the data shown during phenotype-preserving random walks starting from an initial network, by recording the fraction of novel phenotypes that occur in the neighborhood of the changing metabolic network, but not of the initial network. Every data point is an average over 20 random walks each for 10 different initial metabolic networks (thus, 200 random walks in total) at every value of S. Figure 4 shows that the fraction of unique novel phenotypes reaches high values for modest distance between two metabolic networks - small compared to the maximum genotype distance - and does not depend much on the number of sulfur sources S on which viability is required. It also does not depend strongly on metabolic network size (results not shown). In sum, the neighborhood of moderately different metabolic networks contains very different novel phenotypes.

Figure 4

Fraction of phenotypes unique to the neighborhood of an evolving metabolic network G _k (vertical axis) when compared to the neighborhood of a starting metabolic network G after k "mutations", i.e., reaction changes (horizontal axis). Each data point represents an average over 200 evolving metabolic networks of size N = 200 (20 random walks for 10 different sets of environmental demands with the same number S of required sulfur sources).

The ability of metabolic networks to encounter novel phenotypes does not depend monotonically on their phenotypic robustness

The question of how robustness relates to phenotypic variability has raised considerable interest in recent years [19, 20]. Macromolecules - RNA and protein - whose phenotypes are more robust to mutations can access more novel phenotypes than less robust phenotypes [7, 8]. This holds for both large and small evolving populations of such molecules, at least in RNA [8]. We next asked whether these observations are specific to macromolecules, or whether they hold more generally, that is, also for the genotype-phenotype map of metabolic networks.

Above we considered the robustness of a metabolic genotype as its fraction of non-essential reactions. Analogously, we now consider the robustness of a metabolic phenotype as the average fraction of non-essential reactions of all networks with this phenotype [8]. We showed that robustness decreases as networks are required to be viable on more and more sulfur sources (Figure 3B). That is, for networks of any given size, the number S of sulfur sources on which they are viable can serve as a proxy for phenotypic robustness. The greater a phenotype's S is, the smaller is its robustness.

When analyzing how evolving populations explore a genotype network, we need to distinguish between two different kinds of populations. The first are populations where the product of population size and mutation rate is much smaller than one. For brevity, we refer to such populations here as small populations. The second are populations where this product is much greater than one. We refer to these as large populations.

Small populations are genotypically monomorphic most of the time [15], and effectively explore a genotype network much like a single changing network would, i.e., through a random walk on the genotype network. During such a random walk, the changing network encounters different phenotypes in its neighborhood. We determined the cumulative number of different novel phenotypes found in the neighborhood of a random walker. That is, if a phenotype was encountered twice, either in the same neighborhood, or in a neighborhood encountered during an earlier step, we counted it only once. We did so for networks of varying size N and number of sulfur sources S. Specifically, for each N and S, we carried out 200 random walks of 10 000 mutations each (20 walks for 10 different sets of sulfur sources at each S). Figure 5A shows the resulting data. The cumulative number of novel phenotypes is a unimodal function of S, indicating that metabolic networks under few and many environmental demands encounter fewer novel phenotypes than under an intermediate number of environmental demands (S≈20). The cumulative number of novel phenotypes depends strongly on metabolic network size for S < 20, where larger metabolic networks encounter more novel phenotypes throughout the random walk. It is not sensitive to N for larger values of S.

Figure 5

Metabolic networks have diverse phenotypes in their neighborhood. (A) Cumulative number of novel phenotypes encountered in the neighborhoods of evolving metabolic networks of different sizes and subject to different environmental demands. (B) The average pairwise genotype distance found in populations of evolving metabolic networks. Each population consists of 100 individual metabolic networks. (C) Number of novel phenotypes found in the (1-mutant) neighborhood of random metabolic networks of different size N and subject to different environmental demands S. Each data point is an average over 200 metabolic networks (20 random walks for 10 different sets of environmental demands, with the same number S of required sulfur sources).

We next examine large evolving populations. Such populations are polymorphic most of the time. To model their evolutionary dynamics, one needs to track every individual in the population, unlike for monomorphic populations. We determined the cumulative number of novel phenotypes that are mutationally accessible to a population of metabolic networks evolving on (and restricted to) a specific genotype network. This number can be determined by examining, for each generation, the neighborhood of each individual in the population, and by counting the total number of different novel phenotypes encountered. We simulated populations of 100 individuals evolving for 2000 generations (see Methods for details). Additional file 3 shows the average number of cumulative unique novel phenotypes accessible to a population through generation 2000. Each data point represents an average and standard deviation over 200 simulations (20 simulations for 10 random sets of sulfur sources at a given S). Qualitatively, the figure resembles our observations for a single random walk (Figure 5A), except that the absolute number of cumulative unique phenotypes encountered is higher in evolving populations.

Taken together, these observations show that the number of novel phenotypes accessible to a population does not increase monotonically with phenotypic robustness. It decreases with increasing robustness (decreasing S) for low values of S, and it increases with robustness at higher values of S. Additional file 4 demonstrates this relationship in a 3 D plot of the number of novel phenotypes versus robustness and network size (A) or environmental demands (B).

We next examined two candidate explanations of this pattern. The first is that environmental demands and network size affect how rapidly a population can diversify on its genotype network, and thus also how many novel phenotypes it can access. To find out whether this diversification rate matters, we examined the average pairwise genotype distance of our evolving populations. The smaller this difference is, the more slowly a population diversifies. Figure 5B shows a plot of pairwise genotype distances, averaged over an entire population, at the end of 2000 generations. One can see that populations of smaller networks are less diverse. However, environmental demand (S) influences genotypic diversity only weakly, and not in the same unimodal way as seen in Additional file 3. Thus, population dynamic processes alone cannot explain the pattern observed in Figure 5A and Additional file 3.

The second candidate explanation is that the patterns of Figure 5A and Additional file 3 may simply reflect how the number of novel phenotypes in the neighborhood of random metabolic networks varies with N and S. Figure 5C shows the number of novel phenotypes in the neighborhood of random viable metabolic networks of varying size, and with varying environmental demands on the network. This figure is based on random samples of 200 metabolic networks (see Methods) for every value of N and S (20 metabolic networks for 10 different sets of sulfur sources at each S). The vertical axis of this figure shows the mean and standard deviation of the number of unique novel phenotypes in the neighborhood of the examined networks. It shows similar unimodal characteristics as the data in Figure 5A and Additional file 3. The figure demonstrates that the number of novel phenotypes depends strongly on metabolic network size for environments with S < 20. In this regime, larger metabolic networks have more novel phenotypes in their neighborhood than smaller networks. For S > 20, the dependency on metabolic network size disappears and the number of accessible novel phenotypes declines again. In sum, the different accessibility of novel phenotypes in evolving populations, at least qualitatively, emerges from how novel phenotypes are distributed in genotype neighborhoods, and how this distribution depends on S and N.

Discussion

The genotype-phenotype map we characterized here shows both similarities and differences to previously characterized such maps [2–6]. One similarity is the existence of connected genotype networks that extend far through genotype space, and that link genotypes having the same phenotype. Connected sets of metabolic networks viable on the same set of sulfur sources exhibit large maximum genotype distances D _max . For example, networks with as few as 200 reactions can show D _max > 0.7, meaning that they share fewer than 30 percent of their reactions. A second similarity regards phenotypic innovations, genotypes whose phenotypes allow viability on novel sulfur sources. The neighborhoods of two genotypes G ₁ and G ₂ tend to contain very different phenotypic innovations, even if G ₁ and G ₂ are only moderately different. Both features, taken together, facilitate the exploration of novel phenotypes. They would allow a population of organisms (networks) to explore different regions of genotype space, preserving their phenotype while exploring many novel phenotypes.

A major difference to previously studied genotype-phenotype maps regards the relationship between a phenotype's robustness to mutation and a population's ability to explore novel phenotypes. In macromolecules, this relationship appears to be positive: Greater robustness facilitates innovation [7, 8]. Although robust molecules can access, on average, fewer novel phenotypes in their mutational neighborhoods, populations of robust molecules can spread faster through genotype space. In balance, the second process dominates and allows evolving populations to access more novel phenotypes through mutations.

In sulfur metabolism, we do not see this relationship. Robust phenotypes in this context are characterized by viability on few sulfur sources. They are less easily disrupted through eliminations of individual reactions. We found that the number of phenotypic innovations such phenotypes can access in their neighborhood - through changes of single reactions - is highest at intermediate robustness, that is, for phenotypes viable on approximately 20 out of 60 carbon sources we examined. (It can also depend on metabolic network size, being lowest for small networks.) This phenomenon cannot solely be explained by the evolutionary dynamics of evolving populations, partly because populations whose members have intermediate robustness do not spread fastest through genotype space. Instead, the phenomenon is a simple consequence of how many novel phenotypes occur in the neighborhoods of individual genotypes. This number peaks for genotypes whose phenotypes have intermediate robustness. It shows the same qualitative dependence on robustness as the number of novel phenotypes accessible to populations. Thus, in this case, population dynamics do not dominate the process of novel phenotype exploration. We note that the total number of possible novel phenotypes decreases exponentially with the number S of sulfur sources on which a network is already viable. If we took this exponential decrease into account, for example by determining the cumulative fraction instead of the number of novel phenotypes accessible to evolving population, this fraction would decrease with increasing S.

These observations raise the question whether they are unique to sulfur metabolism or whether they occur in other metabolic systems. As we stated earlier, part of our motivation to study sulfur metabolism was to avoid the much larger number of reactions of carbon metabolism, which render population approaches like ours computationally intractable. Nonetheless, very limited analyses for carbon metabolism are possible. Additional file 5 panels A and B show the results of such an analysis, based on a small number of populations of networks at moderate size. The analysis has large uncertainties, but it shows a pattern that is at least reminiscent of sulfur metabolism: Innovation peaks at intermediate robustness (the number of alternative carbon sources a phenotype is viable on).

Taken together these analyses show that the organization of different phenotypes in genotype space can differ greatly among different classes of biological systems, such as proteins and metabolic networks. And these differences can affect the ability of a system to explore novel phenotypes in this genotype space.

A third class of analyses regards features that have not been studied previously, partly because they are unique to metabolic systems and our representation of them. One of them regards the analysis of networks with different sizes (numbers of reactions). Our genotype representation can accommodate and allows us to compare networks of different sizes, whereas commonly used representations of other systems - molecules or regulatory circuits - cannot. For example, proteins of different length form genotype spaces of different dimensions, making their comparison challenging [21].

When analyzing metabolic networks of different sizes, we found that populations of small networks can explore fewer novel phenotypes (Figure 5A and Additional file 3). This observation is easily explained if one considers that populations of such networks are less robust. Their genotype can thus be altered less easily. In consequence, they are genotypically less diverse (Figure 5B), which restricts their access to novel phenotypes (Figure 5B).

Another analysis focusing on network sizes is our characterization of minimal metabolic networks, networks in which all reactions are essential. While the process of genome and metabolic network reduction leading to small networks has been studied for specific biological networks [22], our approach does not start from such a network and can thus provides a more systematic exploration of genotype space. In our analysis of random minimal metabolic network viable on the same sulfur sources, we found that such networks can have large genotype distance. We can explain part of this observation through reactions that are very similar but differ in one of several highly related metabolites. For example, in many types of reactions involving the phosphorylation of a metabolite, the phosphor group donor can be any of ATP, ADP, AMP or even other phosphorylated nucleotide bases. This allows single reactions to be substituted by similar reactions that only use another group donor metabolite. Also, many alternate pathways require only the swapping of two reactions allowing metabolic networks with very little robustness to substitute some of their reactions. However, these may not be the only explanations of different network architectures, because minimal metabolic networks viable on the same sulfur sources can have dramatic pathway differences (results not shown). Whether such differences can be bridged through series of single reaction changes is a question for future exploration.

Properties of minimal networks are also useful in explaining the maximal genotype distance in a genotype network. For example, for metabolic networks of a given size N and viability on S sulfur sources, the maximum genotype distance within a genotype network is well approximated by one minus the fraction of superessential reactions in minimal metabolic networks. These are reactions found in all minimal networks viable on a given number of sulfur sources. We currently have no mechanistic explanation for this relationship and it, also, remains a subject for future work.

Studies like ours have several limitations. One of them is that we focus on biomass synthesis phenotypes, and not on other aspects of metabolism, such as secondary metabolite production. The reason is that biomass synthesis has the most immediate impact on an organism's survival. Other limitations include that the addition and deletion of reactions may have effects on fitness even if they do not affect biomass synthesis, that our knowledge of the reaction universe is limited, and that we face uncertainty about the biologically most important sulfur sources, about thermodynamic properties of individual metabolic reactions, and about the role of cellular compartmentalization in guiding sulfur metabolism. Our study, even though it uncovers generic features of genotype-phenotype maps with demonstrated relevance for evolutionary adaptation and innovation in other biological systems [7, 23], is thus best viewed as a modest beginning in characterizing a complex metabolic genotype space.

Conclusions

We demonstrate that metabolic networks in sulfur metabolism with the same phenotype form large genotype networks that reach far through metabolic genotype space. How far they reach through this space is a linear function of the number of super-essential reactions specific to the environment. We show that the robustness of metabolic networks depends on the size of a metabolic network, on the average size of minimal networks viable in a given environment, and on how rapidly the proportion of essential reactions decrease with increasing network size. The neighborhoods of two metabolic networks on the same genotype network typically contain different novel phenotypes. In evolving populations of metabolic networks, robustness facilitates the discovery of novel phenotypes only up to some modest value of robustness, beyond which populations discover fewer novel phenotypes. The difference in the role of robustness in the evolution of metabolic networks compared to its role in the evolution of macromolecules shows that phenotypic innovations may not occur according to the same principles in all biological systems.

Methods

Global set of sulfur-involving reactions

To obtain the global set of reactions involving sulfur-containing metabolites that can be present in the metabolic networks we studied, we used data from the LIGAND database of the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.ad.jp/kegg/ligand.html) [11]. The LIGAND database is a database of chemical compounds and reactions in biological pathways compiled from pathway maps of metabolism of carbohydrates, energy, lipids, nucleotides, amino acids and others. Also included in the database are the list of catalyzed reactions categorized by the Nomenclature Committee of the International Union of Biochemistry and Molecular Biology (NC-IUBMB) http://www.chem.qmul.ac.uk/iubmb/enzyme/ which includes all enzymes with known classification [24].

Specifically, we used the REACTION and COMPOUND sections of the LIGAND database to construct our global reaction set. From this dataset we pruned (i) all reactions involving general polymer metabolites of unspecified numbers of monomer units (C₂H₆(CH₂)_n), or, similarly, general polymerization reactions that were of the form A_n+ B → A_n+1, because their abstract form makes them unsuitable for stoichiometric analysis, (ii) reactions involving glycans, because of their complex structure, (iii) reactions that were not stoichiometrically or elementally balanced, and (v) reactions involving complex metabolites without chemical information about their structure.

In addition, we merged all the reactions existing in the E. coli metabolic network model (iJR904) [25] that involve sulfur containing compounds. After these steps of pruning and merging, our global reaction set consisted of 1221 reactions.

Flux balance analysis

Flux balance analysis is a computational method used to find a set of fluxes through all metabolic reactions that maximize biomass production in a given metabolic network, assuming it is in a steady state [12]. This assumption means that the concentrations of internal metabolites does not change over time. To compute the maximum biomass growth using this method, one needs to know the stoichiometric coefficients of each reaction, the chemical environment of the cell (the set of upper bounds on the fluxes of external metabolites into the cell), and the biomass composition, which represents metabolite consumption during cell growth. This consumption is reflected in a "biomass growth reaction", for which we chose the reaction defined in the E. coli iJR904 metabolic model [25]. This biomass growth reaction includes all 20 proteinaceous amino acids, nucleotides, deoxynucleotides, putrescine, spermidine, 5-methyltetrahydrofolate, coenzyme-A, acetyl-CoA, succinyl-CoA, cardiolipin, FAD, NAD, NADH, NADP, NADPH, glycogen, lipopolysaccharide, phosphatidylethanolamine, peptidoglycan, phosphatidylglycerol, phosphatidylserine and UDPglucose. For the purpose of this study we concentrated only on the ability of a metabolic network to synthesize the sulfur containing biomass precursors, which are the two amino-acids cysteine and methionine, coenzyme-A, acetyl-CoA and succinyl-CoA. We thus allowed the metabolic networks to uptake any metabolites not containing sulfur. We consider a metabolic network to be viable in a given environment if it can sustain a biomass growth rate greater than 1.0 × 10^-3. In essence, the approach we take is equivalent to asking whether all the necessary sulfur containing biomass precursors are synthesizable given a metabolic network in a specified environment. Flux balance analysis relies on linear programming [26] to compute the maximum biomass production rate. We used the packages CPLEX (11.0, ILOG; http://www.ilog.com/) and CLP (1.4, Coin-OR; https://projects.coin-or.org/Clp) to solve the associated linear programming problems.

Environments and phenotypes

We here considered 124 different environments that differed in the chemical compound that could serve as the sole source of sulfur. These 124 sources were all the sulfur containing metabolites in the 1221 reactions of our global reaction set. We provided any metabolite not containing sulfur in the environment, in effect making it a rich environment limited by sulfur containing metabolites only. Also, we allowed cells to secrete all metabolites. We define a metabolic phenotype as the set of environments (each with a different sole sulfur source) in which a metabolic network is viable. The environmental demands imposed on a metabolic network correspond to the set of sulfur sources that the metabolic network must at least be viable in.

Essential and super-essential reactions

We define a reaction as essential if its removal from a metabolic network renders the metabolic network inviable on at least one of the sulfur sources that it had previously been viable on. We called a reaction super-essential if it occurred in all minimal metabolic networks generated under a given set of environmental demands.

Generating random and minimal metabolic networks

We generated random viable metabolic networks as follows. First, we generated a random environmental demand, that is, we required viability in some given number X of sulfur sources. To this end, we first created a binary vector of length 124 (each of whose entries corresponds to one sulfur source), initialized all its entries to the value zero, and then randomly changed X of these entries to one. These entries represent the set of sulfur sources on which we required our metabolic networks to be viable.

We then generated random viable metabolic network of N reactions as follows. We started from a metabolic network that contained all 1221 reactions (this networks is viable on all 124 sulfur sources) and sequentially removed randomly chosen reactions, while ensuring viability on the set of X sulfur sources chosen previously, until we had reached a network with the target number N of reactions.

We define a minimal metabolic network as a network where not a single reaction can be removed without destroying viability. To generate a (random) minimal metabolic network we used the same procedure until no reactions could be removed without destroying viability.

Metabolic network random walk maintaining viability in the environmental demands

We generated random walks for metabolic networks of given reaction numbers N and viability on a given number of sulfur sources by first generating a random metabolic network of this size, as just described. We then generated a series of steps ("mutations") in metabolic genotype space, each one either an addition or a deletion of a reaction. After each step, we recomputed the phenotype of the metabolic network. If the metabolic network was still viable on the same set of sulfur sources, we accepted the mutation and proceeded to the next mutation; if not, we rejected the mutation and repeated the process from the metabolic network prior to the mutation. We continued the resulting random walk for 10 000 accepted mutations. We kept the size of the metabolic network in the narrow interval (N, N+1) by ensuring that accepted mutations alternated between reaction additions and deletions.

In a variation on this procedure, we also carried out forced random walks through genotype space. Their aim was to obtain metabolic networks that are as different (in terms of genotype distance) as possible from the initial metabolic network. In a forced random walk, we required that any reaction addition did not involve a reaction that had been part of the initial network at the start of the walk.

Population dynamics

Populations where the product of population size and mutation rate is much greater than one are polymorphic most of the time, and show evolutionary dynamics different from those of small populations [27]. To understand their evolution, one needs to simulate them explicitly. To this end, we implemented a Fisher-Wright model of evolution [28] in populations of 100 metabolic networks. We initialized each population with 100 copies of a single viable metabolic network, and then exposed it to repeated "generations" of mutation (one reaction addition or deletion per network and generation) and selection. Specifically, for the selection procedure, we chose 100 viable individuals at random with replacement to form the next generation. A mutation that rendered a network inviable could not be chosen. Our simulations proceeded for 2000 generations.