### Simulated evolution of metabolic pathways

To evaluate the role of mutation-selection-drift balance in biochemical pathway evolution, a population of cells with a key metabolic pathway was evolved under different selective schemes. The simplified kinetic model designed to capture features of glycolysis [35] and methylglyoxal metabolism [36] is shown in Fig. 1. The glycolysis-like aspects of the pathway include the feedback loop (as an approximation to glycolysis regulation) and the synthesis of final metabolite F as analogous to pyruvate in a linear pathway. The methylglyxoxal-like pathway elements include the toxic intermediate (B) as analogous to methylglyoxal (a highly toxic intermediate) and again, the synthesis of the final metabolite (F) is analogous to pyruvate.

This model is expressed in terms of a system of ordinary differential equations where reactions are described by reversible Michaelis-Menten kinetics. Each enzyme has parameters for enzyme concentration [Enzyme] (mmol/l), the catalytic constant (k_{cat}) (mmol/l/s), the Michaelis constant for the substrate (K_{M}) (mmol/l), the reversible catalytic constant (k_{catr}) (mmol/l/s), and the Michaelis constant for the product (K_{Mr}) (mmol/l). The kinetic model has a single inhibitory reaction that is described in the system by the inhibition constant K_{I} (mmol/l). The COPASI [37] modeling environment is used to solve this system of equations. The steady state solution is used, with a constantly replenishing concentration of A and mass action to utilize F, as described in Additional file 1: Table S1.

In order to model the evolutionary process, a forward time simulation with discrete generations is employed. In general, the simulation represents each individual in the population as an instance of the described model, subjecting this model individual to mutations which may elicit fitness effects, and then sampling individuals based on fitness to populate the next generation using weighted sampling with replacement. The pathway architecture remains unchanged during the course of the simulations. These simulations proceed by establishing an initial population of 100 homogenous individuals with parameter values given in Additional file 1: Table S1. Because a set of differential equations must be solved for each individual in each generation, the population size was limited by computational capacity. It is not expected that the results obtained in this study are driven by the size of the population. Each forward simulation was repeated 5 times. Mutations were introduced with a probability of 3*10^{−3} per parameter per individual per generation. K_{M} and k_{cat} were treated as evolving independently, although there is a mechanistically unpredictable degree of dependence (and link to protein stability) in their evolution in nature from current understanding, as described below. The mutational effect on the catalytic rate constant and enzyme concentration (both indicated by p below) are derived from a normal distribution with variable mean \( {\mu}_{n_1} \), where

$$ {\mu}_{n_1}=-0.01{e}^{c*\cdot {p}_{n_1-1}} $$

The mutational effects on the binding constants (K) are described by a standard normal distribution with a variable mean \( {\mu}_{n_2} \),

$$ {\mu}_{n_2}=\frac{1}{-0.01{e}^{c\cdot {K}_{n_2-1}}} $$

The index value c is used to scale the mutational effects, with the following values for each constant:

$$ c=\left\{\begin{array}{c}\hfill 2.5\times {10}^{-2},\wedge enzymeconcentration\hfill \\ {}\hfill 2.5\times {10}^{-2},\wedge inhibitionconstant\hfill \\ {}\hfill 1.0\times {10}^{-2},\wedge catalyticconstant\hfill \\ {}\hfill 3.{3}^{\prime}\times {10}^{-4},\wedge reversiblecatalyticconstant\hfill \\ {}\hfill 1,\wedge productconstant\hfill \\ {}\hfill 3.{3}^{\prime}\times {10}^{-2},\wedge reversableproductconstant\hfill \end{array}\right. $$

This mutational scheme allows for scaling across orders of magnitude in kinetic parameters and generates a distribution of mutational effects with a bias towards slightly degrading change that is dependent upon the activity and expression level of the protein. The mutational scheme is consistent with current thought in molecular evolution, where the range and distribution of mutational effects are influenced by the current state [22]. Most of the mutations are slightly deleterious or neutral, while advantageous mutations are rare, although slightly less so as the activity of the molecule decreases. Intuitively, as a sequence decreases in fitness contribution, the number of sequences with higher potential fitness contribution increases and as it increases in fitness contribution, the number of sequences with a higher potential fitness contribution decreases, as expected by Fisher’s geometric model.

Five different selection schemes were employed to examine the influence of various factors on pathway evolution. The first scheme involved selection on steady state flux alone, where the fitness of an individual is described below:

$$ {F}_1=\frac{1}{1+{e}^{-{\left( flux-650\right)}^{0.07}}} $$

Values in this logistic function control the asymptotic fitness and the gradient of the flux to fitness relationship. As enzymes reach limits of adaptation because of the ability to utilize products, so do pathways, where the end products are also subjected to the rules of binding and catalysis [23, 38, 39]. The asymptote of 650 and slope of 0.07 are arbitrary, but are chosen to reflect the ultimate utilizable flux. Changing them would be expected to alter the distribution of fitness effects (fraction of deleterious changes at equilibrium), but not the overall evolutionary dynamics of the system. A second (negative control) scheme was implemented to examine mutational opportunity and mutational pressure. In this experiment individuals were sampled at random from the population and only the mutational process acted.

Another control was used to examine the evolutionary stability of rate limiting steps, by implementing a scheme with selection on the first reaction rate to become rate limiting by preventing the buildup of the intermediate after the reaction, and using a neutral mutational distribution (with respect to molecular phenotype) that eliminated mutational pressure. We used the multiplicative fitness function,

where

$$ {F}_2=\frac{1}{e^{s\cdot \left[B\right]}} $$

Here, [B] is the concentration of the deleterious metabolite and s (9.4 × 10^{−4}) is a scalar chosen to control the flux and the intersection point of the two curves. As indicated, the mean of the mutational distribution is set at 0, and the distribution is parameter-independent.

A fourth experiment was implemented to examine the role of preventing the buildup of the deleterious intermediate on pathway evolution, resulting in the same multiplicative fitness function above. This experiment used the biological (parameter dependent) mutational distribution as previously outlined. Finally, the cost of protein production was also considered using another multiplicative fitness function where *s* is a normalizing constant (1.0 x 10^{−6}), cost_{AA} (30.3) and cost_{nuc} (49.2) reflect the per unit costs of synthesis [24], and enzyme lengths are given in Additional file 1: Table S2,

where

$$ {F}_3=\frac{1}{1+s\cdot \left( cos{t}_{protein} + cos{t}_{mRNA}\right)} $$

and

$$ \begin{array}{l} cos{t}_{protein}= cos{t}_{AA}\left\{\left[ EnzymeA\right] lengt{h}_A+\left[ EnzymeB\right] lengt{h}_B+\left[ EnzymeC\right] lengt{h}_C\right.\\ {}\left.+\left[ EnzymeD\right] lengt{h}_D+\left[ EnzymeE\right] lengt{h}_E\right\}\end{array} $$

$$ \begin{array}{l} cos{t}_{mRNA}= cos{t}_{nuc}\left\{\frac{3\cdot lengt{h}_A\cdot \left[ EnzymeA\right]}{1000}\right.+\frac{3\cdot lengt{h}_B\cdot \left[ EnzymeB\right]}{1000}\\ {}\left.+\frac{3\cdot lengt{h}_C\cdot \left[ EnzymeC\right]}{1000}+\frac{3\cdot lengt{h}_D\cdot \left[ EnzymeD\right]}{1000}+\frac{3\cdot lengt{h}_E\cdot \left[ EnzymeE\right]}{1000}\right\}\end{array} $$

Each of these simulations was run for 22,000 generations and the point of mutation-selection balance was reached by generation 20,000 under each of these selective schemes (the scheme with no selection did not reach equilibrium because there is no mutation-selection balance without selection). The point of mutation-selection balance was determined by the stability of the fitness of the median individual across generational time as assessed by observation of approximately equal rates of positive and negative changes (Additional file 1: Figures S20–S21). The point of balance was confirmed for the experiment with selection on flux alone by replicate experiments approaching the same point from lower fitnesses that were reached from higher fitnesses (Additional file 1: Figure S1).

### Identification of rate limiting steps

The sensitivity of each of the reactions across the last 2000 generations was determined by reducing each reaction rate of the median individual by 10 % while fixing the rest of the reaction rates. The difference in flux between the perturbed and unperturbed systems was used a measure of sensitivity, and the most sensitive step was determined by the reaction for which this value was the largest.

### Examination of evolution and coevolution

Examination of parameter evolution and co-evolution was based upon the values in the median individual at each generation for the 2000 generations after equilibrium was reached. Since the reversible and inhibitory reaction constants have minimal impact on the system, they were removed from the analysis. Five replicates of the same experiment were analyzed together and the rate of change of each parameter was calculated for every generation. In order to control for directional change within enzyme concentrations, catalytic, and binding constants, the average amount of change is calculated for each group and removed from each parameter within the group. 10,000 replicates were bootstrapped from this dataset by random re-sampling within each replicate and complete linkage clustering was performed using absolute correlations as a measure of relatedness between the rates of change (Additional file 1: Figures S22–S26) [40]. The largest clusters significant at the 0.05 level are used to identify co-evolving parameters.

### Simulations with variable population sizes

In order to evaluate the effects of population size on the evolutionary stability of rate-limiting steps, experiments with two different population sizes (10^{2} and 10^{6}) were performed, where the small population size was a control for the results of the previous set of experiments. For this purpose, several simplifications to the procedure were made for computational tractability. In each generation, a single mutation was proposed per generation, with mutational effects and fitness as above for the experiment with selection on pathway flux. The Kimura fixation probability was used to evaluate the fixation of proposed mutations, eliminating an explicit population and any probability of multiple segregating changes. We have

$$ \psi =\frac{1-{e}^{-2c{N}_esp}}{1-{e}^{-2c{N}_es}} $$

representing the fixation probability, where N_{e} is the population size, c is the ploidy (haploid, c = 1), s is the selective coefficient (f’/f_{0}-1, where f’ is the fitness after mutation and f_{0} before) and p is the initial frequency of the allele in a population. The initial frequency p was set to ½ rather than 1/N for computational efficiency, giving the property that a neutral mutation has a 50 % chance of fixation, which scales the selective coefficient. The effects of population size played out in rising from a 0.5 frequency to fixation and the introduction of new mutations was independent of population size.

The population scheme was run for 200,000 generations per experimental replicate and the rate-limiting step length was calculated as was previously described. Both population sizes were run for 30 replicates.

### Simulations with more thermodynamic realism

To evaluate the effects of biophysical constraints on the reaction landscape, simulations where mutations were constrained by Haldane’s relationship were performed for the 10^{6} population size. Although more degrees of freedom are possible with regulation, multiple substrates and products, and the involvement of cofactors, the simplest expression shows that four parameters which are non-independent as,

$$ Keq = \frac{k_{cat}*{K}_{Mr}}{k_{cat r}*{K}_M} $$

Here, K_{eq} is the equilibrium constant driven by the thermodynamics of the reaction. For this experiment, kinetic parameter initial values were set according to Haldane’s relationship (Additional file 1: Table S3)_{.} To maintain the ratio, the mutational scheme was modified from that for other experiments described above. Mutations for K_{Mr} and K_{M} are drawn from a normal distribution with a mean at −1 %. The mutational effect for k_{cat}is also drawn from a normal distribution with a mean at −1 % and has a modifier that is dependent on the original ratio of k_{cat} and k_{catr} and the ratio of mutated K_{M} and K_{Mr}. K_{catr} is calculated from Haldane’s relationship with the mutated k_{cat}, K_{M}, K_{Mr.} This experiment was replicated for 30 times. Rate-limiting step lengths were evaluated as was previously described.

### Characterizing allele segregation with explicit populations

In order to estimate the observed allele segregation within each population, the number of alleles for each parameter was calculated. Parameter numbers were calculated within each population for 2000 generations every 10 generations (total of 200 data points for each parameter) as the mean of total number of alleles per generation (for all parameters and for forward reaction-only parameters). Mean and standard deviation per 2000 generations were retrieved for each parameter as well as the minimum and maximum of the dataset. These values were compared with the expected allele segregation number as calculated for the population with a selectively neutral regime as previously described by Kimura and Crow [25],

Here n is the number of alleles for particular parameter, μ is mutation rate, and Ne the effective population size.

### Statistical tests and bootstrap confidence intervals

For the simulation experiments exploring the evolutionary stability of the rate limiting step, we performed permutation tests under the null hypothesis of no stability. Stability was measured by the number of consecutive generations that a reaction remained rate limiting, once it became the rate limiting step, and under the null hypothesis, each reaction should have the same average number of consecutive generations. For each of the replicates, the correspondence between each reaction and its average time spent as rate limiting was permuted, and the average absolute deviation of each reaction from the overall mean was calculated. In this manner, a null distribution was generated through 100,000 permutation replicates, and an empirical p-value was found by comparing the average absolute deviation in the actual data as compared to this null distribution.

Confidence intervals for each average number of consecutive generations that a reaction was rate limiting were constructed by first bootstrapping the replicates, and then bootstrapping the consecutive runs within each selected replicate. In this manner, a Monte Carlo sampling distribution for each average was generated, and 95 % confidence intervals were generated by taking the 2.5th and 97.5th percentiles from each bootstrap sampling distribution. Error bars in the corresponding figures reflect these confidence intervals.

To test the question of whether selection has an effect on expression cost, we examine enzyme concentration, k_{cat} and K_{M} against the length of the enzyme. Using the 2000 generations at equilibrium across the five replications and comparing the selection against flux experiment to the selection against flux and protein expression, we ran a linear mixed-effects model with random effects for the replicates, and an interaction term for enzyme length and experiment. The null hypothesis is that the interaction term is equal to 0, or in other words, that the different selection regimes have the same effect on expression cost. Due to computational burden within the mixed-effects model when attempting to account for the correlation structure induced by the Markovian nature of the simulations, only the unique values of each corresponding outcome variable were selected to be fit in the model, and standard linear mixed-effects models were run with a random effect for each replicate. While this loss of information is suboptimal, this should result in conservative inference. Assumptions of homoskedasticity and linearity of the relationship between enzyme length and parameter observations appear to be satisfied. The small effective sample size may be a concern, but should also lead to conservative inference. Statistical significance of the interaction term was assessed via likelihood ratio tests comparing the interaction model to the null model, which did not contain the interaction term.