Modeling the evolution of a classic genetic switch

Background The regulatory network underlying the yeast galactose-use pathway has emerged as a model system for the study of regulatory network evolution. Evidence has recently been provided for adaptive evolution in this network following a whole genome duplication event. An ancestral gene encoding a bi-functional galactokinase and co-inducer protein molecule has become subfunctionalized as paralogous genes (GAL1 and GAL3) in Saccharomyces cerevisiae, with most fitness gains being attributable to changes in cis-regulatory elements. However, the quantitative functional implications of the evolutionary changes in this regulatory network remain unexplored. Results We develop a modeling framework to examine the evolution of the GAL regulatory network. This enables us to translate molecular changes in the regulatory network to changes in quantitative network function. We computationally reconstruct an inferred ancestral version of the network and trace the evolutionary paths in the lineage leading to S. cerevisiae. We explore the evolutionary landscape of possible regulatory networks and find that the operation of intermediate networks leading to S. cerevisiae differs substantially depending on the order in which evolutionary changes accumulate; in particular, we systematically explore evolutionary paths and find that some network features cannot be optimized simultaneously. Conclusions We find that a computational modeling approach can be used to analyze the evolution of a well-studied regulatory network. Our results are consistent with several experimental studies of the evolutionary of the GAL regulatory network, including increased fitness in Saccharomyces due to duplication and adaptive regulatory divergence. The conceptual and computational tools that we have developed may be applicable in further studies of regulatory network evolution.


Background
Regulatory networks are known to underlie many biological processes, and therefore their characterization and analysis forms a central focus of systems biology [1][2][3][4]. Despite their importance, relatively little is known about how regulatory networks are formed during evolution and shaped by natural selection.
One of the best studied regulatory networks in molecular biology is the "GAL network", which is responsible for the inducible metabolism of galactose in budding yeast. In addition to being extremely well-characterized in S. cerevisiae [5][6][7] it has also been the subject of a number of quantitative modeling efforts [8][9][10][11] and evolutionary studies, which have revealed many interesting patterns of regulatory network evolution [12][13][14][15].
Perhaps most general of these evolutionary paradigms is the duplication and divergence of function of an ancestral bi-functional gene. Before the whole genome duplication (WGD), which occurred along the lineage leading to S. cerevisiae [16,17], the GAL network is thought to have employed the Gal1/3 protein to perform both an enzymatic and regulatory function. This bi-functional protein has been retained in species that diverged before the WGD [18]. Following the WGD, the two gene copies of GAL1/3 were converted to a highly inducible enzyme (GAL1) [19] and a weakly inducible regulatory protein (GAL3) [20] in post-WGD species. Recent work has demonstrated that these changes confer a growth advantage in galactose [21], and suggested that gene duplication and subsequent neofunctionalization represented an example of how a regulatory network can overcome an adaptive conflict, whereby the network cannot improve in one aspect without impairing function in another.
Here we set out to explore the consequences of evolutionary changes in the GAL network using a model of the regulatory network ( Figure 1) to relate specific sequence changes to changes in quantitative function. Starting with experimentally characterized regulatory differences between the pre-and post-WGD GAL genes (see Methods), we inferred an ancestral organization of the regulatory network (see Results). Simulations of the ancestral and S. cerevisiae networks show that gene duplication and specialization lead to elevated gene expression in the presence of galactose, and decreased gene expression in the absence of galactose, thus leading to an improved switch-like system in S. cerevisiae. We then use the model to explore the significance of the order in which a set of maximally parsimonious evolutionary events separating a post-WGD ancestor and S. cerevisiae occur, and find important consequences for the function and evolution of the switch. We introduce the idea of evolutionary paths in the space of possible regulatory networks and develop quantitative measures to compare paths. We find that there are evolutionary paths in the lineage to S. cerevisiae that optimize particular features of their constituent network intermediates, some of which have been shown to be directly related to fitness. Perhaps more importantly, we find that it does not seem possible to optimize all network features in any single path.

Results
A hybrid stochastic-deterministic modeling framework for examining the evolution of the GAL regulatory network In order to explore the functional consequences of evolutionary changes in the GAL network, we sought a quantitative approach to relate these molecular changes to changes in network operation. Because the evolutionary events of interest include variations in gene copynumber, changes in protein function, as well as changes in the number and spacing of cis-regulatory sequences, we required a modeling framework in which we could incorporate all these features and relate them to network behavior.
To model transcription/translation we modified and implemented a general physicochemical model of transcription regulation [22,23] in a stochastic context. This allowed us to vary the regulatory features of promoters driving the GAL genes and hence to predict the effects of evolutionary changes in the cis-regulatory sequence of a particular promoter. To model changes in gene copy number, we simply changed the corresponding probability of transcription, for example, to model a gene duplication event, we multiplied the probability of that gene's transcription by two. Because of the difference in time scales and species numbers between transcription/translation reactions and protein-protein/galactose interactions [24], we modeled the latter using a system of ordinary differential equations. In order to model the effect of evolutionary changes in protein function, we assigned different sets of reaction possibilities to new molecular species created during the evolution of the network.
Hybrid approaches [25][26][27][28] such as this permit the inclusion of gene expression noise [29,30], multiple time scales, and variable system mass in biological models while making simulations computationally feasible. Briefly, each step of the hybrid simulation algorithm consists of a stochastic part for the slow transcription/ translation/degradation processes, and a deterministic part for the fast molecular interactions which is solved to equilibrium (Figure 2a, see Methods section for a detailed description of the simulation algorithm).
Gal4p -Gal80p Figure 1 Schematic of the GAL regulatory network model. A transcriptional activator protein, Gal4p, binds a consensus upstream activating sequence (UAS) in the promoter region of structural and regulatory genes promoting their expression. In the absence of galactose, however, Gal4p forms a complex with the Gal80p repressor molecule, which masks its C-terminal activation domain, preventing the recruitment of transcriptional machinery. Galactose is imported from the environment by non-specific hexose transporters, as well as by highly-specific Gal2p permeases, and activates the co-inducer Gal3p molecule or the bi-functional Gal1/3p if it is present in the network. The activated form of Gal3p (and/or Gal1/3p) can then bind and sequester Gal80p, relieving the inhibition of gene expression and permitting induction of the GAL cluster. GAL1 has a core of three UASs arranged in an in-phase helical manner. This configuration allows cooperative interactions (black connecting lines in figure) between adjacently bound activating (Gal4p) and inhibiting (Gal4p-Gal80p complex, not shown) protein pairs, which enhance gene expression or inhibition in induced or uninduced conditions respectively. The promoter region of GAL1/3 also contains this core, but, unlike GAL1, these sites are arranged in an anti-phase helical configuration which does not permit cooperative interactions. GAL3 has retained only one UAS, resulting in decreased inhibition in uninduced conditions and a relatively small induction in the presence of galactose. In addition to the changes in the number and arrangement of cis-regulatory elements, differences also exist between the intrinsic strengths of ancestral and extant promoters (see Table 4). Note that cooperative interactions between the two UASs on GAL2 are permitted [5], and that the 4 th UAS on GAL1 (and GAL1/3) is not included in the model after indications that it may not be as important for gene expression as the other UASs [42].   Model parameters were chosen to reproduce the quantitative operation of the GAL network in S. cerevisiae ( [5,6,9,21], see Methods for details of parameter choice and parameter sensitivity). Typical traces from simulations are shown in Figure 2b. Protein distributions were obtained at equilibrium from 100 independent simulations of each variant network at constant extracellular galactose concentrations in the range [10 -8 , 10 -1 ] M in the absence of glucose (Figure 2c). The induction-response curve for S. cerevisiae in increasing concentrations of galactose displays the characteristic switch-like response ( Figure 2d) and fold-inductions in gene expression are consistent with previously published experimental data ( Table 1 in Methods).
Reconstructing the operation of an inferred ancestral network GAL1 and GAL3 are paralogous genes which have diverged after a WGD event in the lineage leading to S. cerevisiae from a common bi-functional ancestor, GAL1/3, which serves as both the galactokinase and coinducer molecule. It has recently been shown that the fitness of cells growing in galactose is related to the phenotype and operation of the GAL network. The cisregulatory evolution of GAL1 and GAL3 has led to a more efficient genetic switch by subfunctionalizing the processes and regulation of galactose phosphorylation and induction [21].
To study the effects of these changes on the network's quantitative function we inferred an ancestral form of the network by removing the GAL1 and GAL3 genes from the S. cerevisiae network and substituting a bi-functional GAL1/3 gene driven by the ancestral promoter for Klac-GAL1 [21]. We have assumed that the other ancestral regulatory genes have remained unchanged in function and copy-number compared to S. cerevisiae.
We find that, relative to the inferred ancestor, the number of Gal1p (galactokinase) molecules in S. cerevisiae is increased in high galactose and reduced in very low galactose, which implies a more switch-like response in (105 fold in S. cerevisiae vs. 6 fold in the ancestor). At the same time, the number of Gal3p (inducer) molecules is decreased in all conditions relative to the ancestral protein (e.g., at very low galactose, there are~5000 molecules in S. cerevisiae vs.~27000 in the ancestor, Figure 2d). Thus, a single bi-functional gene seems to perform poorly in controlling the enzymatic and regulatory aspects of the network, consistent with experimental results [21].
The effect of an evolutionary change on the quantitative function of the GAL network depends on the network's evolutionary history Having established the quantitative differences in network function between the ancestral and extant GAL network in S. cerevisiae, we next sought to explore the significance of the order of evolutionary events after the WGD on network operation.
Five events were considered: the regulatory and functional divergence of the two GAL1/3 copies, and the loss of the GAL2, GAL4, and GAL80 gene duplicates, according to the principle of maximum parsimony [31]. Since we know that there was an ancestral genome duplication event, we expect each gene to have been at two copies in some ancestor. We then make the maximally parsimonious assumption that each other gene was lost once. This is the minimum number of evolutionary changes that could have occurred, though not necessarily the scenario that was actually followed. Similarly, we know that there are three active promoter binding sites in the extant GAL1 promoter, and three in the ancestral GAL1/3 promoter, so we assume that were no additional binding site gains and losses during the specialization of GAL1/3 to GAL1. Likewise, we know that the GAL3 promoter has a single binding site, so we assume that two sites in GAL1/3 were lost without any additional intervening binding site gains or losses. Finally, since there is little data from other species to estimate parameter changes in the GAL network, we have assumed that these either remained unchanged or have incurred small changes (see Methods for parameter perturbation experiments).
We simulated all 33 of the possible network configurations that can arise as combinations of the regulatory Mean numbers of molecules +/-standard deviation is reported from equilibrium protein distributions of 100 simulations in very low galactose (uninduced) or very high galactose (induced).
changes required for the ancestral network configuration to evolve to the extant configuration in S. cerevisiae (1 pre-WGD ancestor + 1 post-WGD ancestor + 5 configurations with one event + 10 configurations with two events + 10 configurations with three events + 5 configurations with four events + 1 S. cerevisiae. See Table 2). These networks link S. cerevisiae and its ancestor via 120 unique evolutionary paths. Each path consists of 7 networks which are connected via the WGD event and a unique sequence of five evolutionary changes ( Figure 3a). To evaluate the quantitative function of a network we considered the number of protein molecules at equilibrium in very low galactose (10 -8 M, the "uninduced state") and conditions of very high galactose (0.1 M, the "induced state"). We defined three network features of interest for further examination: Feature 1 "repression strength", is defined as the number of galactokinase molecules (Gal1p, Gal1/3p) in the uninduced state (10 -8 M galactose), such that a smaller number indicates better repression; Feature 2 "induction strength", is defined as the number of galactokinase molecules in the induced state (0.1 M galactose), such that a larger number indicates stronger induction; Feature 3 "switch effectiveness" is defined as the number of co-inducers (Gal3p, Gal1/3p) in the uninduced state (10 -8 M galactose), such that a smaller number indicates a more effective switch.
Changes in gene copy-number, protein function, and the number and arrangement of cis-regulatory sequences often resulted in significant changes in the network's quantitative function, as indicated by the differences in the steady-state galactokinase and co-inducer concentrations between the 33 network configurations (Figure 3b). The magnitude and direction of these changes, however, depended both on the particular evolutionary event being effected as well as on the configuration of the network being changed -that is, on the history of evolutionary changes that a network has accumulated. For example, we found that loss of the GAL80 repressor gene duplicate has a significantly smaller impact on repression strength if it occurs prior to the specialization of both copies of GAL1/3 (data not shown). This is because Gal80pmediated repression of GAL1 is much stronger than that of GAL1/3 due to differences in the promoter regions of the two genes.
Quantitative assessment of the order of evolutionary changes reveals that there are evolutionary paths that optimize specific network features All of the evolutionary paths in the space of possible regulatory networks terminate at S. cerevisiae and, as a consequence, eventually accrue identical changes in function relative to the inferred ancestral network. Nevertheless, because the effect of an evolutionary change depends on the configuration of the network at that point in evolution, intermediate GAL networks leading to S. cerevisiae occupy markedly different regions in functional space depending on the sequence in which evolutionary changes accumulate.
For each of the 120 evolutionary paths we applied a scoring scheme (Figure 4a) which penalizes an Table 2 Configurations for the GAL network variants modeled in this study Gene copy-number Evolutionary changes, such as the WGD, gene duplicate loss, and specialization of GAL1/3 to GAL1 and GAL3 connect these networks along evolutionary trajectories to S. cerevisiae (see Figure 4).  The path that optimizes repression strength (0 penalty) consists of the sequence: specialization of GAL3, specialization of GAL1, loss of GAL4 duplicate, loss of GAL2 duplicate, loss of GAL80 duplicate (red path in Figure 3a, b). The sequence: specialization of GAL1, specialization of GAL3, loss of GAL4 duplicate, loss of GAL2 duplicate, loss of GAL80 duplicate, optimizes switch effectiveness. The sequence: specialization of GAL1, loss of GAL2, loss of GAL80, specialization of GAL3, loss of GAL4 (green path in Figure 3a, b) optimizes induction strength.
We found that evolutionary paths where the two specialization events of GAL1/3 to GAL1 and GAL3 occur in close sequence perform well in maximizing repression strength and switch effectiveness for intermediate GAL networks. Path scores in these two features are negatively correlated (R 2 = 0.85) with the number of evolutionary events separating the specialization of GAL1 and GAL3 (Figure 5a). In addition, paths where the second copy of GAL80 is lost after all other evolutionary changes show significantly higher repression strength compared to other paths (-0.34+/-0.19 vs. -1.96 +/-1.05, P = 0.001, two-sample t-test, Figure 5b). It was also found that the group of paths where specialization of GAL1/3 takes place before any gene duplicates are lost has a higher average switch effectiveness score than paths where duplicates are lost prior to specialization (-0.16+/-0.12 vs. -1.24+/-0.57, P < 0.001, two-sample t-test, data not shown). Finally, we find a correlation (R 2 = 0.72) between how early GAL1 specializes and the strength of induction (Figure 5c), as well as a weaker correlation with how late GAL3 specializes and the induction strength scores (R 2 = 0.34, data not shown).
We next explored the relationship between the scores of the different features. Evolutionary paths tend to maximize switch effectiveness scores and repression strength scores together (R 2 = 0.72, data not shown). Interestingly, however, paths cannot optimize induction strength and switch effectiveness/repression strength simultaneously (Figure 5d). As a result, evolutionary paths may optimize switch effectiveness/repression or strength of induction, but not both features at the same time. Score for 1st observed transition x to be maximized: x to be minimized: In the preceding sections, we have made the assumption that the parameters of the network remain the same as those observed or inferred in S. cerevisiae. However, in order to test whether the results presented above are sensitive to parameters, we have performed a series of perturbation experiments in which parameters were randomly modified during evolution and found that the observations presented were all recapitulated in the perturbed networks (see Methods).

Discussion
Using a quantitative model of the regulatory network we were able integrate the wealth of experimental knowledge about the S. cerevisiae GAL network with promoter swap experiments and to infer and simulate the operation of an ancestral form of the GAL network. Two of our observations are consistent with the report of increased fitness after promoter specialization [21]. First, our observation of a higher expression of GAL1 relative to the ancestor could lead to an elevated metabolic capacity (galactose flux through the galactokinases) for S. cerevisiae in high galactose, and therefore perhaps to increased growth rate. Second, we observed a decrease in the number of protein molecules in the absence of galactose, which could lead to increased growth in the absence of galactose through a decrease in the energetic cost of protein synthesis (smaller number of protein molecules) (see [32] for an account of cost-benefit analysis of gene expression in a different system). Despite the apparent consistency between our results and the report of adaptive specialization along the S. cerevisiae lineage, we found it difficult to establish a rigorous relationship between the quantitative function of the network and the network's fitness landscape. For example, in order to understand the potential benefit of increased levels of GAL1 expression, it is necessary to consider the impact of the galactokinase on the entire metabolic network [33]. In particular, consideration of the toxicity of the metabolic intermediates is important and greatly complicates this analysis [34]. Perhaps even more difficult is that we do not know the relative importance to the population of the induced or uninduced states of the GAL network. Clearly, the amount of time that S. cerevisiae may spend exposed to galactose depends on the environment and history of the population. It is unclear if this information will ever be available to us.

Analysis of evolutionary paths
Our scoring of evolutionary paths in the functional space indicates that most potential phenotypic divergence takes place in the uninduced state (repression strength and switch effectiveness, Figure 5b). This is consistent with experimental findings that report a greater effect of cis-regulatory variation on gene repression than on induction [35]. We also found that no evolutionary path takes optimal steps (steps which incur no penalty) in every functional dimension at every network intermediate. In particular, we found that paths that maximized repression strength and switch effectiveness were sub-optimal with respect to induction strength.
Interestingly, we note that GAL networks in other extant species have specialized their GAL1/3 homologues, but have retained various duplicates. For example, S. bayanus and S. castellii have specialized forms of GAL1 and GAL3, but S. bayanus has kept two GAL80 copies and reacquired a tandem GAL2 duplicate [15], while S. castellii has kept two GAL4 and GAL80 copies [36]. According to our evolutionary path assessment, the specialization of GAL1 and GAL3 before the loss of regulatory gene duplicates, as well as the close proximity of these specialization events, improves switch effectiveness and repression strength. Retention of the GAL80 duplicate in both species is also consistent with better repression strength scores (Figure 4a, b).
We speculate a possible explanation for these observations. The GAL network allows for a larger phenotypic divergence in the uninduced state (consistent with our observations of larger variance in path scores pertaining to switch effectiveness and strength of repression). The network's capacity to improve and, therefore, the growth advantage, is greater if networks follow the path of uninduced state optimization.

Assumptions of the approach
Our method of exploring the evolution of the GAL network is based on several assumptions. First, in constructing the evolutionary paths leading to S. cerevisiae, we assume that the evolutionary distance between two GAL network intermediates that are connected by a regulatory change is always constant -that is, we treat all evolutionary changes as equally likely to occur at each step of a path. Second, we make the assumption that regulatory elements not present in S. cerevisiae were also not present in the ancestor. While network features not present in S. cerevisiae, such as a unique UAS configuration on a gene promoter, could have existed in the ancestor or intermediate networks, we have constrained the network's regulatory space to include only those evolutionary events required under the assumption of maximum parsimony -which states that, when confronted with multiple evolutionary scenarios that explain some data, the scenario that requires the fewest evolutionary events is the most probable to have been followed [31]. Third, while we attempt to model parameter perturbations during network evolution, we cannot exclude the possibility that mutations can occur that have very large effects on parameter values. Potentially, experimental characterization of the GAL network in other species [37] could give an indication of such events, which could then be modeled as distinct evolutionary changes in the history of the GAL network.

Conclusions
A hybrid stochastic-deterministic modeling framework has been used to explore the effect of regulatory divergence in the GAL network following the whole genome duplication. We have shown that some evolutionary paths optimize distinct features of network operation and that intermediate GAL networks in the lineage to S. cerevisiae could not have optimized all network features at the same time.

Hybrid stochastic-deterministic simulation algorithm
Our model of the GAL regulatory network contains 18 molecular species with extracellular galactose acting as the input to the system ( Table 3). Simulations of the GAL network were performed using a hybrid stochasticdeterministic algorithm where each iteration of the algorithm consists of a stochastic part and a deterministic part (the latter of which may be skipped as an approximation method).

Deterministic part
The first part of the algorithm begins with an algebraic step which models galactose transport to the cell interior: where we used a = 0.001 and b = 0.01 as transport coefficients for the non-specific hexose transporters and galactose permeases (Gal2p) respectively. In the above equation, and in equations that follow, molecular names in square brackets refer to the concentration of that species, whereas names without brackets refer to that species' number of molecules.
Protein-galactose and protein-protein interactions involved in signal transduction are modeled as a system of ODEs: where f P i is the promoter model assigned to gene GAL i (see Promoter models section). Our model accounts for changes in gene copy-number because the probability of transcript production is proportional to the number of gene copies for that molecular species in the network. This is a result of the stochastic simulation algorithm, where multiplying the promoter function by gene number is equivalent to introducing additional identical reaction channels. The propensity function for protein degradation, d, is where g is the degradation rate. The Gillespie algorithm samples the joint probability distribution of reaction events and times at each iteration of the stochastic part, selects a reaction, and generates a time interval containing no reactions, τ, which is used to advance the simulation time.

Promoter models
The propensity function of protein synthesis in the stochastic part of the simulation algorithm is driven by a physicochemical promoter model. An upstream activation sequence (UAS) on a promoter may be empty, occupied by a Gal4p molecule, or occupied by a Gal4-80p molecule. Sites occupied by Gal4p interact with and recruit RNA polymerase and the transcriptional machinery with an affinity that is intrinsic to the particular promoter (Table 4). Empty sites and sites occupied by Gal4-80p may also allow some basal transcription to occur. The promoter function, f PI , for transcription is the sum of the contribution of all sites to transcription initiation in each occupancy configuration, weighted by the configuration's probability: Where c index the occupancy configurations for the promoter driving gene GAL i . K c is the sum of contributions to transcription initiation by all sites in configuration c. n4 c and n480 c are the numbers of Gal4p-and Gal480p-bound sites in configuration c respectively. b (1.68 kcal/mol) is the product of the gas constant and the system's temperature (300 K). ΔG c is the sum of the Gibbs' free energy changes from all DNA-protein binding events in configuration c, as well as any cooperative interactions between adjacently bound activator and/or repressor proteins which stabilize the DNA-protein complex. Taking the configuration where all sites are empty as the reference state of each promoter, we used -13.86 kcal/mol as the energy change for the binding of either a Gal4p or Gal4-80p to a UAS [39]. We chose -2.00 kcal/mol and -3.00 kcal/mol for energy changes due to a single cooperative interaction between appropriately positioned Gal4p and Gal480p respectively because they led to proper quantitative promoter function in S. cerevisiae, and we could not identify these parameters in the literature. In the absence of any UAS on a promoter, as in the case for the GAL4 promoter, the propensity function evaluates to a constant, Kc, giving constitutive expression.

Algorithm execution
The model was implemented and simulated in MATLAB ® (The MathWorks). Simulations were run for 8,000 seconds, which was found to be much longer than the typical time for equilibrium (data not shown). To decrease computation time -particularly in systems with high protein concentrations -the deterministic step was executed only every n iterations, where n scaled linearly with the total number of molecules participating in the deterministic step. The rationale behind this approximation is that stochastic synthesis or degradation of a few molecules will not significantly change the equilibrium of interacting proteins when large concentrations of proteins are already reacting. Trial simulations run with and without this approximation procedure followed similar induction dynamics, attained  [5]). We also do not consider the 4 th UAS on GAL1 (and GAL1/3) after indications that it may not be as important for gene expression as the other UAS [42].
equilibrium at the same times, and had similar equilibrium behavior (data not shown). Conversions between number of molecules and concentrations are performed at each step as required using 1.66 × 10 -15 and 23 × 10 -15 L as nuclear and cellular volumes respectively [40,41].
The mean and variance for the number of molecules for each species was obtained from the equilibrium distributions of 100 independent simulations at 8 extracellular galactose concentrations (10 -8 , 10 -7 , 10 -6 , 10 -5 , 10 -4 , 10 -3 , 10 -2 , and 10 -1 M). This was done for each of the 33 GAL networks to construct their response curves to increasing galactose concentrations.

Scoring scheme for evolutionary paths
Each evolutionary event is scored according to its effect on some feature -for example, the number of molecules with galactokinase activity at a certain extracellular galactose concentration. We begin scoring evolutionary events after the genome duplication and consider the 120 unique combinations of the five events required for the post-genome duplication ancestor to evolve to S. cerevisiae. For each feature to be scored, we defined either the largest fold-increase or fold-decrease to be the optimal change (see Results section). Events are scored relative to the possible alternative events at each intermediate network along the path. An observed event is penalized if it does not result in the best possible change in the feature from amongst the possible events. To keep the scoring scheme consistent, if the feature is to be maximized, we used: where x a,i , is the fold-change in the network feature resulting from the i-th evolutionary step after the genome duplication, for the a-th alternative network. We indicate the observed evolutionary choice in the path at step i, as a = o. See Figure 4a for an illustration of the scoring scheme.
The scoring scheme provides an objective means of ranking evolutionary paths in regulatory space that optimize some network feature. The path with the highest possible score (which is 0) consists of the order of events which results in the best possible sequence of changes in the scored feature. Conversely, the worstscoring path follows none of the best alternatives (except for the last evolutionary event). 4. Binding energies (promoter models): Additive factor in [-1, 1] from U(-1,1).
where U is the uniform distribution. These random parameter perturbations can be thought of as minor mutations in the network components which do not form part of the parsimonious set of evolutionary changes along the lineage to S. cerevisiae, but which may nevertheless have occurred and subsequently been lost.
Under parameter perturbations, evolutionary path score distributions had larger variance, often skewed toward worse scores, with occasional perturbations completely abolishing the switch-like behavior of the system (Additional File 1, Figures S1-S3). Most perturbations, however, did not have a large impact on perturbed path scores.
Specifically, we find that the number of evolutionary events separating the specialization of ScGAL1 and ScGAL3 again correlates negatively with repression strength scores (R 2 = 0.85, Additional File 1, Figure  S4a). The number of events preceding the specialization of ScGAL1 correlates negatively with induction strength scores (R 2 = 0.64, Additional File 1, Figure S4b). A significant difference in repression strength scores was again found between evolutionary paths were GAL80 duplicate loss is the last event and all other paths (-0.40 +/-0.23 versus -2.08 +/-1.03, P < 0.01, Additional File 1, Figure S4c). Perturbed evolutionary paths tend to maximize switch effectiveness scores and repression strength scores together (R 2 = 0.73, data not shown). Finally, we report that, even with a perturbed parameter set, there exists no single path which optimizes all three aspects of quantitative network behavior (Additional File 1, Figure S4d).

Additional material
Additional file 1: Path score distribution comparisons and assessment of evolutionary paths under parameter perturbations. Box plot comparisons of the three evolutionary path score distributions between unperturbed and perturbed paths (Figures S1-3). Correlations between evolutionary path scores and classes of paths under parameter perturbations ( Figure S4).