Bistability in feedback circuits as a byproduct of evolution of evolvability

Abstract

Noisy bistable dynamics in gene regulation can underlie stochastic switching and is demonstrated to be beneficial under fluctuating environments. It is not known, however, if fluctuating selection alone can result in bistable dynamics. Using a stochastic model of simple feedback networks, we apply fluctuating selection on gene expression and run in silico evolutionary simulations. We find that independent of the specific nature of the environment-fitness relationship, the main outcome of fluctuating selection is the evolution of increased evolvability in the network; system parameters evolve toward a nonlinear regime where phenotypic diversity is increased and small changes in genotype cause large changes in expression level. In the presence of noise, the evolution of increased nonlinearity results in the emergence and maintenance of bistability. Our results provide the first direct evidence that bistability and stochastic switching in a gene regulatory network can emerge as a mechanism to cope with fluctuating environments. They strongly suggest that such emergence occurs as a byproduct of evolution of evolvability and exploitation of noise by evolution.

Figure 1

(A) Cartoon representation of the model incorporating a single gene auto-regulatory network. (B) Steady-state analysis of the system showing the production (f_p, in black) and degradation (f_d, in blue) curves for G. The dashed and solid black lines correspond to parameter values a=1.0, b=1, K_D=50 and N=0 (the starting condition for evolutionary simulations) and a=0.31, b=4.7, K_D=52 and N=5 (mean parameters resulting form one of the simulations under v=0.05), respectively. The solid and open circles indicate the stable and unstable steady states of the system.

Figure 2

Mean fitness of the population averaged over the last 5000 generations of the evolutionary simulations, μ_w, versus its variance, σ_w. Open circles and triangles represent results from simulations embedding a deterministic and stochastic phenotype, respectively. Simulations with different environmental fluctuation rates, v, are color coded as indicated by the legend. Results from additional simulations with smaller and higher values of v are omitted to achieve clarity. (A–C) Results from simulations run with mutation rate set to 0.01, 0.001 and 0.0001, respectively.

Figure 3

The steady-state behavior of the system with the mean values of parameters as averaged over the last 5000 generations of the evolutionary simulations. The different panels give the phase plane of the system (i.e., steady-state level of G versus system parameters N, a and b) for different values of the parameter K_D as shown on each panel. On each panel, the level of G is color coded, where locations with split colors indicate bistable regime with two distinct steady-states levels of G. Solid circles and triangles correspond to results from simulations embedding a deterministic and stochastic phenotype, respectively. Results shown are from simulations with v=0.05 and mutation rate set to 0.01. Note that to achieve the mapping of the results on the phase plane, the mean value of K_D obtained from each simulation is rounded to the nearest tenth. Only results with a K_D value rounded to the interval 40–60 are shown for clarity and are representative of results with other K_D values (in particular, with respect to values of N).

Figure 4

(A) Evolvability and population mean of N over environmental epochs for a sample simulation with stochastic phenotype and implementing deterministic environmental fluctuations. The environment switches every 20 generations and mutation rate is set to 0.01. Evolvability is defined as the ratio of the normalized fitness change (i.e., adaptation) over the sum of relative changes in parameters (i.e., genetic shift). The red points show the actual evolvability data, calculated for each epoch, while the black line gives its moving average over epochs. The blue dotted line gives population average of N over epochs. (B) The distribution of individual parameters, fitness and expression level of G over the population and over generations. Each row on each panel encodes a distribution for a specific quantity (as indicated at the bottom of the panel) and for the specific generation shown on the y axis. The distributions are shown as a heat map ranging from red (highest density) to blue (lowest density). Environments E_low and E_high are indicated as black and white bars on the y axis of the left-most panel.

Figure 5

Boxplots showing a summary of the distribution of evolvability from 10 simulation runs for each fixed value of N as shown on the x axis. Note that at N=2.5, we have a bistable system. In these simulations, the other parameters of the system were free to evolve as before. The environment is deterministic and switches every 20 generations. Mutation rate is set to 0.01. Evolvability is defined as the ratio of the normalized fitness change (i.e., adaptation) over the sum of relative changes in parameters (i.e., genetic shift). (A) Results from simulations with stochastic phenotype. (B) Results from simulations with deterministic phenotype.

Figure 6

Analysis of mutational effects of parameters on the phenotype (i.e., level of protein G) in genotypes with increasing value of N. This analysis shows that as N increases, so does (i) the diversity in the expression level of G and (ii) the phenotypic effects of a given mutational increment in other parameters. Rows from top to bottom show phenotypic effects of mutations in parameters b, a and K_D, respectively. Columns from left to right show genotypes with increasing value of N as indicated at the top. The x axis on each panel gives the level of protein G, while the y axis gives the deviation in a particular system parameter based on its mutational increment in the evolutionary simulations (see Materials and methods). The initial parameter is the middle point of the y axis (i.e., zero deviation) and the upper and lower halves of each panel correspond to mutational deviations from that value. The basal values used on each panel are the same; b=1.8, a=0.8 and K_D=50 (different basal values gave similar qualitative results as those shown). For each value of the parameter given on the y axis, we run 1000 independent simulations for 10 generations and encoded in color the number of simulations that converged to the corresponding level of G indicated on the x axis. Solid dots correspond to simulations run with the deterministic phenotype.