Robustness of the attractor landscape.

<p>Plots of the probability that a percentage <i>q</i> of the network attractors is conserved after a gene the knockout of one gene for critical networks constructed <i>de novo</i> (A) and the critical networks that result from the evolutionary process (B). The differently colored distributions in B correspond to populations that started in different dynamical regimes (as in <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002669#pcbi-1002669-g001" target="_blank">Fig. 1</a>). Note the high probability for the <i>de novo</i> networks to lose all their attractors by a gene knockout () which does not happen for the evolved networks ( in all cases). Conversely the probability to conserve all attractors is considerably larger for the latter than for the former ( and, respectively). These data were computed from populations of 1000 networks and 500 attractors per network.</p

Criticality Is an Emergent Property of Genetic Networks that Exhibit Evolvability

<div><p>Accumulating experimental evidence suggests that the gene regulatory networks of living organisms operate in the critical phase, namely, at the transition between ordered and chaotic dynamics. Such critical dynamics of the network permits the coexistence of robustness and flexibility which are necessary to ensure homeostatic stability (of a given phenotype) while allowing for switching between multiple phenotypes (network states) as occurs in development and in response to environmental change. However, the mechanisms through which genetic networks evolve such critical behavior have remained elusive. Here we present an evolutionary model in which criticality naturally emerges from the need to balance between the two essential components of evolvability: phenotype conservation and phenotype innovation under mutations. We simulated the Darwinian evolution of random Boolean networks that mutate gene regulatory interactions and grow by gene duplication. The mutating networks were subjected to selection for networks that both (i) preserve all the already acquired phenotypes (dynamical attractor states) and (ii) generate new ones. Our results show that this interplay between extending the phenotypic landscape (innovation) while conserving the existing phenotypes (conservation) suffices to cause the evolution of all the networks in a population towards criticality. Furthermore, the networks produced by this evolutionary process exhibit structures with hubs (global regulators) similar to the observed topology of real gene regulatory networks. Thus, dynamical criticality and certain elementary topological properties of gene regulatory networks can emerge as a byproduct of the evolvability of the phenotypic landscape.</p> </div

Gene expression variability <i>α</i>.

<p>This parameter, which measures whether the genes are frozen in one state, either 0 or 1, or if they more or less switch back and forth between these two states, can be computed in two distinct ways. The first way (horizontal variability ) is to measure along each attractor state, as shown in A, and then average over all the attractor states and over all the attractors in the attractor landscape. The second way (variability ) shown in B, is to measure the variability for each gene throughout time along the attractor cycle, then average over all the genes in the network and over all the attractors. This is illustrated by the particular example for the same attractor where , whereas .</p

Survival times of the different strains in the population.

<p>(<b>A</b>) Plot of the network labels (strains) that are present in the population at a given generation. Each horizontal line indicates the survival time of a particular strain. The vertical lines indicate the fixation events in which all the networks in the population are relabeled after only one strain was left in the entire population. (B) Distribution of survival times computed during generations (black curve). This distribution was computed using logarithmic bins. Only data for are presented because we checked the existence of strains every 20 generations. The red dashed line is the best fit which corresponds to the power-law . The inset shows the corresponding cumulative distribution , which better reveals the goodness of the power-law fit. The fact that has a more or less power-law behavior implies that almost all the strains disappear from the population very quickly, whereas only very few networks are able to survive the Darwinian selection mechanism given by the ACC, the AIC and the <i>α</i>-fitness criterion.</p

Mutations in the regulatory region of a gene.

<p>This second type of mutation affects the way in which a given gene regulates its targets. (A) If is one of the regulators of a target gene . Then mutation in the coding region of may afford gene the capacity to bind to a new site in the regulatory region of the same target gene (B), or abrogate the capacity to bind to an existing binding site of that target gene (C). Conversely a change in the coding region of gene may provide new binding capacity for a site in the regulatory region of a new target, as it is shown in (D) and (E).</p

Efflux Pump Regulatory Network.

<p>Arrows indicate positive regulation. Blunt arrows indicate repression. A) Literature base reconstruction of the AcrAB-TolC efflux pump regulatory network of <i>Escherichia coli</i> as reported on [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0118464#pone.0118464.ref009" target="_blank">9</a>]. B) Simplified version of the AcrAB-TolC efflux pump regulatory network (EPRN). The activator (Act) and repressor (Rep) are two Transcriptional Factors that belong to the same transcriptional unit (EPRN operon, indicated by the dashed line). When the repressor occupies its DNA binding site, the expression of the operon is restrained. Nonetheless, when the antibiotic (or inducer, <i>Ind</i>) enters the cell, it inactivates the repressor by binding to it, allowing the operon to be actively transcribed, promoting the production of pumps and decreasing the synthesis of porins (this last process is known to occur through an intermediary). Both food and inducer are expelled by the efflux pump system. In the population model, a reduction in food concentration implies an increase in the division time.</p

Network structure generated by the evolutionary process.

<p>The top-left network shows the structure of the giant component of the transcription factor interaction network of <i>E. coli</i> according to the RegulonDB <a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1002669#pcbi.1002669-GamaCastro1" target="_blank">[3]</a>. This network has <i>N</i> = 101 nodes and average connectivity <i>K</i> = 2.46. The structure on the top-right corresponds to the typical network that results from our evolutionary algorithm, which in this particular case has <i>N</i> = 100 nodes and average connectivity <i>K</i> = 1.85. Note the existence of global regulators, i.e. nodes with a great number of output connections. The bottom panel presents in a log-log plot the out-degree distribution of these two networks to illustrate their remarkable similarity.</p

Evolution towards criticality.

<p>(A) Evolution of the average network sensitivity for four different populations, each initially composed of networks in one of the three dynamical regimes: ordered (, black), critical (, red), and chaotic (, green; and , blue). Under the Darwinian selection given by the ACC and AIC, all the populations quickly become critical (), regardless of their initial dynamical regime. The inset shows that convergence towards criticality occurs during the first 10000 generation steps. The control curves (in light gray) were obtained by evolving populations without selection, and show that the mutagenic method alone drives the networks into the chaotic regime (). (B) Distribution of sensitivities at two different generations for the population that started with chaotic networks. In early generations is quite broad (dashed line), reflecting a great diversity of networks. However, through evolution all the surviving networks approach criticality and the distribution narrows down (solid line). The distribution shown here at generation has .</p

Evolution of the network topology.

<p>(A) The common ancestor network has 10 nodes and one of them (node number 9) is a global regulator that regulates 8 other nodes. (B) Diagram of strain survival times showing the first fixation event at generation <i>g</i> = 6411 (indicated by the red arrow). The common ancestor network is the one that gives rise to the population of the first fixation event. (C) Structure of a randomly chosen network in the final population (generation <i>g</i> = 250000). The initial hub (node number 9) is the one marked with the red circle. Note that at the end this is not a hub anymore, but just another ordinary node of the network. (D) Distribution of the link persistence for the 10 connections of the common ancestor. The black and red histograms represent the populations at the first fixation event and at the end of the simulation, respectively. Even after the first fixation event, the links , and almost disappear from the population. Furthermore, in the final population none of the links of the initial hub occur at significant frequency. By contrast, link is present in all the networks of the final population because node 2 became a hub throughout the evolutionary processes. Link is indicated with the blue bold arrow in (A).</p

Reversibility of the resistance phenotype.

<p>(A) This tracking plot shows that the expression of the activator increases while the antibiotic shocks are applied as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0118464#pone.0118464.g002" target="_blank">Fig. 2</a>. Then, when the antibiotic is removed (indicated by the tilted black arrow), the expression of the activator decreases abruptly and eventually reaches its basal level. (B) Size of the population as a function of time for the same simulation as in A. After each antibiotic shock (small black arrows) the population size decreases exponentially and the recovery time becomes longer with each shock. After the antibiotic is removed (tilted black arrow) the population comes back again to its wild-type (WT) growth rate. To carry out the simulations in the antibiotic-free phase, every time the population reached the maximum size N = 5000, we took a random sample of 10% of the cells and made them grow without antibiotic, until the population reached again this maximum size, and so on. (C) Average transcription rate μ_β = ⟨β₀⟩ in the population as a function of time. Note that the average increases while the shocks are applied and then gradually comes back to small values when the antibiotic is removed. Error bars indicate the standard deviation. It can be observed that the standard deviation increases with the antibiotic stress. The panels below show the full distribution G (μ_β, σ_β) at three different times: before any antibiotic is introduced (circle); after several antibiotic shocks (star); after a long period of time without antibiotic (line). Time is measured generations, being one generation the time it takes for a cell with β₀ = 1 to reach θ_F starting from F = 0.</p