---
abstract: 'We investigate the evolution of Boolean networks subject to a selective pressure which favors robustness against noise, as a model of evolved genetic regulatory systems. By mapping the evolutionary process into a statistical ensemble and minimizing its associated free energy, we find the structural properties which emerge as the selective pressure is increased and identify a phase transition from a random topology to a “segregated core” structure, where a smaller and more densely connected subset of the nodes is responsible for most of the regulation in the network. This segregated structure is very similar qualitatively to what is found in gene regulatory networks, where only a much smaller subset of genes — those responsible for transcription factors — is responsible for global regulation. We obtain the full phase diagram of the evolutionary process as a function of selective pressure and the average number of inputs per node. We compare the theoretical predictions with Monte Carlo simulations of evolved networks and with empirical data for *Saccharomyces cerevisiae* and *Escherichia coli*.'
author:
- 'Tiago P. Peixoto'
bibliography:
- 'bib.bib'
title: 'Emergence of robustness against noise: A structural phase transition in evolved models of gene regulatory networks'
---

Introduction
============

Many large-scale dynamical systems are composed of elementary units which are *noisy*, i.e., can behave non-deterministically, but nevertheless must behave globally with some degree of predictability. A paradigmatic example is gene regulation in the cell, which is a system of many interacting agents — genes, mRNA, and proteins — which are subject to stochastic fluctuations. What makes gene regulation particularly interesting is that it is assumed to be under evolutionary pressure to preserve its dynamic memory against stochastic fluctuations [@kitano_biological_2004; @raser_noise_2005; @maheshri_living_2007]. Many important cellular processes require such reliability, such as circadian oscillations [@raser_noise_2005]. Furthermore, in multicellular organisms, errors in signal transduction can potentially lead to catastrophic consequences, such as embryo defects or cancer [@kitano_biological_2004; @gutkind_signaling_2000]. Since the source of noise cannot be fully removed [@lestas_fundamental_2010], a gene regulation system must adopt characteristics which compensate for the unavoidably noisy nature of its elements. Since they are a product of natural selection, these characteristics must emerge from random mutations and subsequent selection based on fitness. A central question concerns the general large-scale features which are likely to emerge in this scenario that result in reliable function under noise. In this work, we study the emergence of robustness against noise in networks of Boolean elements which are subject to selective pressure, functioning as a model for evolved gene regulatory systems. We show that the system undergoes a structural phase transition at a critical value of selective pressure, from a totally random topology to a “segregated-core” structure, where a smaller and more densely connected subset of the network is responsible for the regulation of most nodes in the network. This characteristic is present to a significant degree in gene regulatory systems of organisms such as yeast and *Escherichia coli*, in which all the regulation is done by a much smaller (and denser) subset of the network, comprised of transcription factor genes.

Boolean networks (BNs) have been used extensively to model gene regulation [@kauffman_metabolic_1969; @kauffman_homeostasis_1969; @bornholdt_systems_2005; @drossel_random_2008]. The Boolean value on a given node represents the level of concentration of proteins encoded by a gene, which in the simplest approximation can be either “on” or “off.” The regulation of genes by other genes is represented by Boolean functions associated with each node, which depend on the state of other nodes called the *inputs* of the function. The dynamics on these networks serve as a model for the mutual regulation of genes which control the metabolism of cells in an organism. Gene regulation is composed of specific steps involving the production of proteins and other metabolites, which need to be carried out in specific sequences and under certain conditions. During each of these steps the dynamics is subject to stochastic fluctuations [@raser_noise_2005; @maheshri_living_2007; @raj_nature_2008], since the number of proteins involved can be very low [@raser_noise_2005; @eldar_functional_2010], and the whole process lacks an inherent synchronization mechanism. In order for the regulation process to work reliably, the network must possess some degree of robustness against these perturbations [@kollmann_design_2005]. Indeed, the investigation of real regulatory networks modelled as BNs, such as the one responsible for the yeast cell cycle [@li_yeast_2004], revealed a remarkable degree of robustness, where most trajectories in state space lead to the same attractor, regardless of the initial conditions. Similar results were also obtained for the segment polarity regulatory network in *Drosophila melanogaster* [@albert_topology_2003; @chaves_robustness_2005], which showed that wild-type attractors are significantly robust to different initial conditions and perturbations, and seem to depend only on general topological characteristics of the network, instead of specific functional details. However, the general features which make BNs robust against different types of perturbations are still being identified.

Perhaps the simplest form of perturbation one can consider is a “flip” of a single node in the network, and the propagation of flips which result from it. This corresponds to the situation where the stochastic noise is very weak, and can be modeled as a single flip event. After the perturbation, the system has an arbitrary amount of time to recover (if it recovers), and different perturbations do not build up. Many authors have considered the robustness against perturbations of this type, including Kauffmann [@kauffman_metabolic_1969; @kauffman_homeostasis_1969] who was the first to propose random Boolean networks (RBNs) — networks with fully random topology and functions — as a model of gene regulation. According to this type of perturbation, the dynamics of RBNs [@drossel_random_2008] can belong to one of two phases, depending on the number of inputs per node $K$: A frozen phase ($K<2$) where the perturbation propagates sub-linearly in time and eventually dies out; and a “chaotic” phase ($K>2$) where the perturbation grows exponentially and eventually reaches the entire system. A critical line exists at $K=2$, where the perturbation grows algebraically, and features from both phases are simultaneously observed.

Although RBNs in the frozen phase and on the critical line show features which can be interpreted as robustness in some sense, they fall short of being convincing models for gene regulation. Actual gene regulation networks are not random and show a high degree of topological [@maslov_computational_2005] and functional [@harris_model_2002] organization which are not present in simple RBNs. Conceivably, these features arise out of stringent requirements to perform specific tasks and of types of robustness which are more demanding than the containment of single flip perturbations. As an attempt at producing more complete models, many authors investigated the evolution of BN systems, where the fitness criterion is some form of robustness against perturbation which is not inherent to RBNs. The majority of authors assumed single flips as the only type of noise, but considered different types of responses as fitness criteria [@bornholdt_neutral_1998; @bornholdt_robustness_2000; @stern_emergence_1999; @bassler_evolution_2004; @aldana_robustness_2007; @szejka_evolution_2007; @mihaljev_evolution_2009; @pomerance_effect_2009], most of which are related to the capacity of the network to display the same dynamical pattern after a single-flip. In particular, in [@szejka_evolution_2007] it was found that if the fitness criterion is the ability to return to the same attractor after the perturbation, the evolved networks always achieve maximum fitness. Furthermore these networks with maximum fitness span a huge portion of configuration space, and show a high degree of variability. This means not only that this type of robustness can evolve, but also that it is not a very demanding task for the evolutionary process.

In this work, we consider the arguably more realistic situation where the perturbations are caused by transcriptional noise, which can be arbitrarily frequent [^1]. In this scenario, the effects of noise can overlap and build up in time. The appropriate fitness criterion remains whether or not the network is capable of performing some predefined dynamical pattern, but this is a task which becomes much more complicated. In fact, it can be shown that for networks which are sparse, i.e., the average number of inputs per node is some finite number, perfect robustness can never be achieved, and some amount of deviations, or “errors,” in the dynamics are always going to exist [@peixoto_redundancy_2010; @peixoto_behavior_2012]. Instead, one measures robustness not only by the amount of existing errors, but also by the ability of the system to not be overtaken by them and consequently lose all memory of its dynamical past — i.e., to become ergodic. This type of robustness is much stronger than, and not necessarily related to the ability of the system to contain single-flip perturbations. This was shown in [@peixoto_noise_2009] for RBNs subject to transcriptional noise, for which neither phase (chaotic or frozen) is robust, and both display ergodic behavior, for any non-zero value of noise.

Furthermore, unlike [@bornholdt_neutral_1998; @bornholdt_robustness_2000; @stern_emergence_1999; @bassler_evolution_2004; @szejka_evolution_2007; @mihaljev_evolution_2009], in this work we also consider the *cost* which is associated with different levels of robustness. It is generally the case that improved robustness can be obtained by introducing redundancy or some other mechanism that counteracts the effect of noise, which increases the overhead in the system. This added overhead can impact negatively on the fitness of the organism, which needs to spend more energy or more time to perform the same task. Therefore the trade-off between overhead and robustness is also driven by the evolutionary process. In this work, this overhead is controlled by fixing the average in-degree during the evolutionary process, which becomes an external parameter. By selecting the appropriate value, one automatically determines a selective pressure that yields the corresponding trade-off.

Our main result is that under transcriptional noise, the selective pressure can have a very noticeable effect on large-scale properties of the system: If it is large enough, it triggers a structural phase transition, where networks change from a random topology to a segregated-core structure, with most nodes being regulated by a smaller and denser subset of the network. This observed segregated-core topology is strikingly similar (even if qualitatively so) to what is observed in most real gene regulation networks; namely, genes are separated into two classes: target genes, and those which regulate transcription factors. Only transcription factor genes are responsible for regulation of other genes, and they are usually orders of magnitude smaller in number than target genes [@nimwegen_scaling_2006].

This work is divided as follows. We begin in Sec. \[sec:model\] by presenting the model, and in Sec. \[sec:evolution\] we define the evolutionary process and map it into an equivalent Gibbs ensemble. We then parametrize the topological characteristics of the system as a stochastic blockmodel in Sec. \[sec:block\], and obtain an expression for its entropy. In Sec. \[sec:free-energy\] we describe the technique used to minimize the free energy. We follow in Sec. \[sec:phase\] with the characterization of the existing phase transition and obtain the phase diagram. We perform comparisons with Monte Carlo simulations in Sec. \[sec:mc\] and with the gene regulatory networks of yeast and *E. coli* in Sec. \[sec:real\]. Finally, we conclude with a discussion.

The model {#sec:model}
=========

A Boolean network [@kauffman_metabolic_1969; @drossel_random_2008] is a directed graph of $N$ nodes representing Boolean variables $\mathbf{\sigma} \in \{1,0\}^N$, which are subject to a deterministic update rule, $$\label{eq:bn_dyn}
  \sigma_i(t+1) = f_i\left(\bm{\sigma}(t)\right)$$ where $f_i$ is the update function assigned to node $i$, which depends exclusively on the states of its inputs. At a given time step all nodes are updated in parallel.

We include noise in the model by introducing the probability $P$ that at each time step a given input has its value “flipped”: $\sigma_j \to 1
- \sigma_j$, before the output is computed [@peixoto_noise_2009]. This probability is independent for all inputs in the network, and many values may be flipped simultaneously. The functions on all nodes are taken to be the majority function, defined as $$\label{eq:maj}
  f_i(\{\sigma_j\}) = 
  \begin{cases}
    1  \text{ if } \sum_j\sigma_j > k_i / 2, \\
    0  \text{ otherwise, }
  \end{cases}$$ where $k_i$ is the number of inputs of node $i$. It is assumed throughout the paper that the values of $k_i$ are always odd [^2]. This is so chosen because odd-valued majority functions are optimal, since no other function performs better against noise [@evans_maximum_2003]. By using Eq. \[eq:maj\], we essentially remove the choice of functions from the evolutionary process and concentrate solely on topological aspects.

Starting from an initial configuration, the dynamics of the system evolves and eventually reaches a dynamically stable regime, where the average fraction $b_t$ of nodes with value $1$ no longer changes, except for stochastic fluctuations which vanish for a large system size [@huepe_dynamical_2002; @peixoto_behavior_2012]. In the absence of noise ($P=0$) there are only two possible attractors (if the network is sufficiently well connected) where all nodes have the same value, which can be either $0$ or $1$. We will consider these homogeneous attractors as representing the “correct” dynamics, and denote the deviations from them as “errors.” More specifically, and without loss of generality, we will name the value of $1$ as an “error,” and the value of $b_t$ as the average error on the system.

The steady-state fraction of errors $b^* \equiv \lim_{t\to\infty} b_t$ (for $b_0=0$) will increase with $P$. For any network with a finite average in-degree there will be a critical value of noise $P^*$ for which the dynamics undergoes a phase transition, and the value of $b^*$ reaches $1/2$, and remains at this value for $P>P^*$ [@peixoto_redundancy_2010]. The value $b^*=1/2$ is special, since it means that the dynamics lost the memory of its initial state, since any other initial value of $b_0$ (including $b_0>1/2$) would lead eventually to this same value of $b^*$. Therefore, the value of $P^*$ marks the transition from a nonergodic to an ergodic dynamics. Robustness against noise is synonymous with nonergodicity, since only in this regime are dynamical correlations not destroyed over time.

BNs with majority functions serve as a paradigmatic model for networks robust against noise, since they are composed of optimal elements, and they show a minimal dynamical behavior in the absence of noise, namely two homogeneous attractors with $\{\sigma_i\} = 0$ or $1$. If robustness cannot be attained for such a system, it is much less likely to be possible for a different system with a another choice of Boolean functions, or displaying a more elaborate dynamical pattern [@peixoto_behavior_2012].

In this work we will consider the value of the steady-state average error $b^*$ as the main fitness criterion governing the survival probability of an organism, since it directly measures the deviation from the situation without noise. Although the phenotype itself, i.e. the dynamics without noise, does not change during the evolutionary process considered here, its *stability*, as measured by $b^*$, does. This translates into an actual fitness criterion, since it is not enough for phenotypes to exist; they must also be stable against perturbations. If they are not, they are not viable in practice, and thus should not be observed.

Evolutionary dynamics {#sec:evolution}
=====================

We suppose that a given BN represents the genotype of a full organism, which can self-replicate and belongs to a population that is subject to an evolutionary pressure. The number of individuals in the population is assumed to be sufficiently large and constant. Individuals replicate a given number of times with a constant rate. Parents die the moment they replicate. The offspring are always initially identical to their parents, but are individually subject to point mutations represented by the matrix $\mu_{ij}$, which defines the probability of mutating from genotype (i.e. network) $i$ to $j$. The offspring survive with probability $a_i$, given by the Boltzmann selection criterion $$a_i \propto e^{\beta f_i},$$ where $f_i$ is the fitness of genotype $i$. The parameter $\beta$ controls the *selective pressure*: For large values of $\beta$ only the very best networks survive, whereas for smaller values most networks do. As mentioned previously, the fitness of a network will be given by the fraction of ones (“errors”) after a sufficiently long time, $b^*_{(i)} \equiv \lim_{t\to\infty} b_t^{(i)}$, for $b_0^{(i)}=0$, as $$f_i = -Nb^*_{(i)}.$$ Thus, the largest fitness a network can have is $f_i=0$, which should be possible only if there is no noise ($P=0$).

We suppose that the global offspring mortality rate is such that the size of the population always remains constant. If we consider that the dynamics occurs in discrete time steps, we can write the probability $\pi_i(t)$ of finding an individual in the population with genotype $i$ at time $t$ as a Markov chain, $$\label{eq:markov}
  \pi_i(t) =  a_i \sum_j \pi_j(t-1) \mu_{ji}.$$ The mutation probabilities $\mu_{ji}$ have a decisive effect on what topologies emerge. Mutations in actual biological systems may result in topological bias, such as gene duplications, which are not reversible and result in networks with broad degree distributions [@ispolatov_duplication-divergence_2005; @enemark_gene_2007]. However, the central aim here it to obtain the most likely topology that arises due to the selective pressure *alone*. For this reason we are more interested in mutations which will lead to all possible networks with equal probability in the absence of selective pressure (i.e. ergodicity). A simple and conventional choice is reversible mutations, $\mu_{ij} = \mu_{ji}$, for which the steady state $\pi_i \equiv \lim_{t\to\infty} \pi_i(t)$ obeys the detailed balance condition: $\pi_i \mu_{ij} a_j = \pi_j \mu_{ji}
a_i$. This is a sufficient condition for the desired ergodicity property, but it is not strictly necessary, since other types of mutations may also fulfill it. However, from this condition we easily obtain that the steady-state probability of finding an individual with genotype $i$ is given by its survival probability, $$\pi_i = a_i =  e^{\beta f_i} / \mathcal{Z},$$ where $\mathcal{Z} = \sum_i e^{\beta f_i} = \sum_i e^{-\beta
Nb^*_{(i)}}$. This corresponds exactly to a Gibbs ensemble of all possible genotypes, with a partition function given by $\mathcal{Z}$, where $Nb^*_{(i)}$ plays the role of the “microstate energy” and $\beta$ is the “inverse temperature” (these are of course only mathematical analogies, since these quantities do not actually represent a physical energy and temperature, respectively). The average intensive “energy” in the ensemble is thus $$\label{eq:bcan}
  b^* = \sum_i b^*_{(i)} e^{-\beta Nb^*_{(i)}}/\mathcal{Z},$$ and the canonical entropy is $$\label{eq:scan}
  \mathcal{S} = -\sum_i\pi_i\ln\pi_i = \ln\mathcal{Z} + \beta Nb^*.$$ The objective is to obtain not only $b^*$ for a given $\beta$, but also the network topologies which characterize the ensemble. Instead of considering all microstates individually (i.e., all possible networks) and computing Eqs. \[eq:bcan\] and \[eq:scan\] directly, we may parametrize the whole ensemble via some macroscopic variables $\{x_j\}$ which sufficiently describe its topological properties. These variables must be chosen so that it is possible to write both $b^*(\{x_j\})$ and $\mathcal{S}(\{x_j\})$ as functions of these variables alone. The entropy can, for instance, be obtained via the microcanonical formulation $$\mathcal{S}(\{x_j\})=\ln\Omega(\{x_j\}),$$ where $\Omega(\{x_j\})$ is the number of different networks given a macroscopic parametrization $\{x_j\}$. The values of $\{x_j\}$ which correspond to thermodynamic equilibrium \[i.e., the steady state of Eq. \[eq:markov\]\] can be obtained by minimizing the “free energy”, $$\label{eq:fe}
  \mathcal{F}(\{x_j\}) = Nb^*(\{x_j\}) - \mathcal{S}(\{x_j\})/\beta,$$ with respect to $\{x_j\}$. This stems from the principles of maximum entropy and minimum energy, for closed systems with fixed energy and entropy, respectively, which need to hold in thermodynamic equilibrium [@callen_thermodynamics_1985]. It should again be emphasized that the theory so far is only a mathematical tool, which, although exact, says nothing about the actual physical thermodynamical properties of the evolved systems, i.e., they have no relation to an actual measurable energy or temperature. Instead, the minimization of Eq. \[eq:fe\] is entirely analogous to obtaining the steady state of Eq. \[eq:markov\] by any other means. However, this approach, together with an appropriate topological parametrization, allows us to obtain the outcome of the evolutionary process on the population, without having to actually implement any dynamics. As will be described in detail in the next section, we will parametrize the ensemble as a general *stochastic blockmodel*, which allows for a wide range of topological configurations, while at the same time allowing for a tractable computation of $b^*$ and $S$, which then can be used to minimize $\mathcal{F}$.

It should also be mentioned at this point that we are interested in the properties of typical networks in the ensemble when the selective pressure $\beta$ is varied, under the restriction that the average number of inputs per node (the average in-degree) ${{\left<k\right>}}$ is always the same. As mentioned in the introduction, this restriction originates from the assumption that a larger number of inputs increases the putative cost for the organism of realizing a regulatory mechanism which depends on more elements. Thus, the value of ${{\left<k\right>}}$ should on its own impact the fitness of the organism, and should also be subject to natural selection. For simplicity, we do not describe the fitness landscape which depends on ${{\left<k\right>}}$ and its evolution, in order to emphasize the effects of robustness against noise alone. Instead, we consider ${{\left<k\right>}}$ as an external parameter, which essentially means that the fitness pressure on ${{\left<k\right>}}$ supersedes that of the other parameters, such that it cannot change during evolution. In this way, we are implicitly considering the cost associated with the robustness achieved by increasing ${{\left<k\right>}}$: the smaller is the value of ${{\left<k\right>}}$ chosen, the larger is the implied fitness penalty of having more connections.

Stochastic blockmodel {#sec:block}
=====================

Simultaneous consideration of all possible networks with a given ${{\left<k\right>}}$ is a tremendous task, due to the gigantic number of diverse configurations which are possible. For arbitrary networks the computation of $b^*$ according to Eq. \[eq:bn\_dyn\] may be very cumbersome, since it may depend on many degrees of freedom. Therefore, we narrow down the allowed subset of possible network topologies to those which can be accommodated in a *stochastic blockmodel* [@holland_stochastic_1983; @faust_blockmodels:_1992; @karrer_stochastic_2011] [^3]. As will become clear in the following, we do so without sacrificing the generality of the approach, since we can progressively add to this model as many degrees of freedom as we desire, and in this way obtain arbitrarily elaborate structures in a controlled fashion.

\[tb\]

![\[fig:blocks\](Color online) Example of a network corresponding to a blockmodel with five blocks. The vertices of each block are labeled with the same color.](block-model.pdf){width="0.5\columnwidth"}

A stochastic blockmodel assumes that the nodes in the network can be partitioned into discrete *blocks*, such that every node belonging to the same block has (on average) the same characteristics. Hence, for very large systems, we need only to describe the degrees of freedom associated with the individual blocks (see Fig. \[fig:blocks\]). By considering a system composed of many blocks, we can describe a wide array of possible topological configurations.

More precisely, a (degree-corrected [@karrer_stochastic_2011]) stochastic blockmodel is a system of $n$ blocks, where $w_i$ is the fraction of nodes in the network which belong to block $i$ (we have therefore that $\sum_iw_i=1$), and $p^i_k$ is the in-degree distribution of block $i$. Thus, the average in-degree is given by ${{\left<k\right>}} =
\sum_{k,i}kp^i_kw_i$. The matrix $w_{j\to i}$ describes the fraction of the inputs of block $i$ which belong to block $j$ (we have therefore that $\sum_jw_{j\to i}=1$). Since the out-degrees are not explicitly required to describe the dynamics (see Eq. \[eq:dynblock\] below), they will be assumed to be randomly distributed, subject only to the restrictions imposed by $w_i$ and $w_{j\to i}$.

We note that, although we have diminished the class of networks which will be accessible by the evolutionary algorithm, we still allow a very large array of possible configurations, which can in principle incorporate arbitrary in-degree distributions, degree correlations [@newman_structure_2003], assortative or disassortative mixing [@newman_mixing_2003], and community structure [@girvan_community_2002], to name only a few properties. As will become clear in the following section, this blockmodel is sufficient to characterize the most important topological property that is relevant for robustness against noise, which is the formation of densely connected central subgraphs.

The value of $b^*$ for a blockmodel
-----------------------------------

Supposing that the number of vertices $Nw_i$ belonging to each block $i$ is arbitrarily large, we can compute the value of $b^*$ using an heterogeneous version of the annealed approximation [@derrida_random_1986], by supposing that at each time step the inputs of each function are randomly chosen, such that the specified block structure given by $w_{i\to j}$ is always preserved [@peixoto_behavior_2012]. If the number of vertices in each block is large enough, we can expect this approximation to become an exact description for quenched networks as well. We can then write the average value of $b_i$ for each block over time as $$\label{eq:dynblock}
  b_i(t+1) = \sum_kp_k^im_k\left((1-2P)\sum_j w_{j\to i}b_j(t) + P\right),$$ which is a system of $n$ coupled maps, where $m_k(b)$ is the probability that the output of a majority function will be $1$, if the inputs are $1$ with probability $b$, and is given by $$\label{eq:maj_p}
  m_k(b) = \sum_{i = {{\lceil k/2\rceil}}}^k{k \choose i} b^i(1-b)^{k-i}.$$ A fixed point of Eq. \[eq:dynblock\] represents the solution of a polynomial system of arbitrary order, and therefore cannot be written in closed form. However, it can be obtained numerically by starting the system at $b_i = 0$, and iterating Eq. \[eq:dynblock\] until a fixed-point $\{b^*_i\}$ is reached. The value of $b^*$ can then be obtained as $b^* = \sum_iw_ib^*_i$.

Blockmodel entropy {#sec:entropy}
------------------

We obtain the entropy of the stochastic blockmodel ensemble [@bianconi_entropy_2009; @peixoto_entropy_2011] by enumerating all possible networks which are compatible with a given choice of $w_i$, $p_k^i$ and $w_{i\to j}$. To make the counting simpler, we ignore the difficulty of forbidding parallel edges, and consider the ensemble of *configurations*, since the occurrence of parallel edges should vanish for large network sizes (see [@peixoto_entropy_2011] for more details). Later we compare the results obtained with Monte-Carlo simulations with parallel edges forbidden, and we find very good agreement.

We begin by enumerating all possible in-degree sequences of each block which correspond to the prescribed in-degree distributions, $$\Omega_d =\prod_i\frac{(Nw_i)!}{\prod_k{(Nw_ip_k^i)!}}.$$ For a given block $i$ with a fixed in-degree sequence, we can count the number of different input choices as $$\Omega^i_e = \frac{E_i!}{\prod_jE_{j\to i}!}\prod_j(Nw_j)^{E_{j\to i}},$$ where $E_i=Nw_ik_i$ is the total number of inputs belonging to block $i$ and $E_{j\to i}=w_{j\to i}E_i$ is the total number of inputs from block $i$ which belong to block $j$. Since the set of inputs of each function is unordered, we still need to divide the whole number of input combinations by $\prod_k(k!)^{N_k}$, where $N_k=N\sum_{i,k}w_ip_k^i$ is the total number of vertices with in-degree $k$. Putting it all together we have $$\Omega = \Omega_d \frac{\prod_i \Omega^i_e}{\prod_k(k!)^{N_k}}.$$ Taking the logarithm of this expression, and the limit $N \gg 1$, and using Stirling’s approximation, we obtain the full entropy (up to a trivial constant term, which is not relevant to the minimization of the free energy), $$\begin{gathered}
  \mathcal{S}/N = {{\left<k\right>}}\ln N + \sum_iw_i\mathcal{S}_i^k \\
  - \sum_iw_ik_i\sum_jw_{j\to i}\ln \left(\frac{w_{j\to i}}{w_j}\right),\end{gathered}$$ where $\mathcal{S}_i^k$ is an entropy term associated with the degree distribution of block $i$, and is given by $$\mathcal{S}_i^k = -\sum_kp_k^i(\ln p_k^i + \ln k!).$$

Choice of single-block in-degree distribution
---------------------------------------------

We want to constrain the number of degrees of freedom in the model, such that only the average in-degree $k_i$ of each block is specified, not the entire distribution. In this way, graphs with many different global in-degree distributions are still possible by composing different blocks with different $k_i$’s, but we have a finite number of degrees of freedom per block. In order to obtain the in-degree distribution of the individual blocks, we maximize the entropy $\mathcal{S}$, with the restriction that the average in-degrees are fixed. For that, we construct the Lagrangian, $$\Lambda = \mathcal{S} - \sum_i\lambda_i'\left(\sum_kkp_k^i - k_i\right) - \sum_i\mu_i\left(\sum_kp^i_k-1\right),$$ where $\{\lambda_i\}$ and $\{\mu_i\}$ are Lagrange multipliers which keep the averages and the normalizations constant. We note that the sum over $k$ is made only over *odd* values of $k$, due to the imposed restrictions on the majority function. Obtaining the critical point $(\{\frac{\partial\Lambda}{\partial p^i_k}\},
\{\frac{\partial\Lambda}{\partial \lambda'_i}\},
\{\frac{\partial\Lambda}{\partial \mu_i}\}) = 0$, and solving for $\{p^i_k\}$ one obtains, $$\label{eq:pk}
  p_k^i = \frac{1}{\sinh \lambda_i} \frac{\lambda_i^k}{k!},$$ where $k_i = \lambda_i / \tanh \lambda_i$. Equation \[eq:pk\] is a Poisson distribution, which is defined and normalized only for odd values of $k$.

This choice of $p^i_k$ is not necessarily the optimal one. In fact, it is possible to show that single-valued distributions with zero variance tend to provide the best error resilience [@peixoto_behavior_2012]. Nevertheless, the improvement over a Poisson distribution is *very* small, and the definition of Eq. \[eq:pk\] allows for the average $k_i$ to be continuously varied, which is very practical for the optimization of the free energy.

Block splitting, decrease of entropy and the necessary number of blocks {#sec:split}
-----------------------------------------------------------------------

For the blockmodel defined in this section, there are $2n + n^2$ variables which define the topology, where $n$ is the number of blocks. In order for arbitrary topologies to be faithfully represented by the model, one would need to make $n\to\infty$, which would render this approach impractical. However, we will show that for the purpose at hand, only a minimal number of *two* blocks is sufficient to fully characterize the evolutionary process, without relying on any approximations. This is due the following two facts: 1. Any possible value of $b^*$ can be obtained with only two blocks; 2. Any other topology with the same $b^*$ will invariably have a lower entropy, and thus a larger free energy. Thus the minimum of the free energy will always lie on a two-block structure.

The first fact can be shown by construction: Consider a system of two blocks, where one of them (the “core”) is smaller and much denser, and the remaining block has an average in-degree close to the global average. The inputs of the core block belong mostly to the core itself, as well as the inputs of the remaining block. By changing the density of the core block, as well as the degree of input segregation, it can be shown [@peixoto_behavior_2012] that any possible value of $b^*$ can be achieved [^4].

The second fact can be shown by considering a system of many blocks, and selecting any two blocks, $l$ and $m$. If all other blocks are kept intact, it can be shown that the entropy will always be larger if these two blocks are merged into an effective single block. This can be shown by partially maximizing the entropy $S$ via the Lagrangian, $$\Lambda = \mathcal{S} - \mu\left(\sum_ik_iw_i - {{\left<k\right>}}\right) - \sum_i\gamma_i\left(\sum_jw_{j\to i}-1\right),$$ where $\mu$ and $\{\gamma_i\}$ are Lagrange multipliers which keep the average in-degree and the normalization of $w_{j\to i}$ fixed, respectively. Obtaining the critical point $(\frac{\partial\Lambda}{\partial w_{l\backslash m}},
\frac{\partial\Lambda}{\partial k_{l\backslash m}},
\{\frac{\partial\Lambda}{\partial w_{{l\backslash m}\to j}}\},
\{\frac{\partial\Lambda}{\partial w_{j\to {l\backslash m}}}\}) = 0$, and solving for $w_{l\backslash m}, k_{l\backslash m}, \{w_{j\to l\backslash
  m}\}, \{w_{l\backslash m\to j}\}$ one obtains, $$\begin{gathered}
  k_l = k_m \\
  \frac{w_{l\to j}}{w_l} =  \frac{w_{m\to j}}{w_m} \\
  w_{j\to l} = w_{j\to m}.\end{gathered}$$ This corresponds to the situation where the nodes from blocks $l$ and $m$ are topologically indistinguishable, i.e. the outgoing and incoming edges are randomly distributed among the nodes of both blocks. Since any arbitrary many-block structure can be converted into a single-block by successive block merges, it follows directly that any arbitrary many-block structure can be constructed by starting from a single block, and successively splitting blocks. Thus, since the merging of blocks always increases entropy, and the splitting decreases it, the entropy of the final structure must be smaller than that of both the initial single block and the succeeding two-block network.

In order for a many-block structure to have a lower free energy than the two-block structure with the same value of $b^*$, it needs to have a larger entropy. But according to the above argument, networks with a larger number of blocks tend to have *smaller* entropy. Networks with larger entropy and number of blocks would have to be more randomized than the two-block structure, which would invariably result in a larger value of $b^*$. We can therefore conclude that the global minimum of the free energy always occurs for a two-block structure, and thus we need to deal with only eight variables [^5].

Minimization of the free energy {#sec:free-energy}
===============================

Although we have an analytical expression for the entropy $\mathcal{S}$, the value of $b^*$ cannot be obtained in closed form, and thus the same is true for $\mathcal{F}$. Therefore we must resort to a numerical computation of $b^*$, via the iteration of Eq. \[eq:dynblock\], and minimize $\mathcal{F}$ with a gradient descent algorithm, using finite differences. Many of these methods work only for unconstrained optimization problems, and we need to impose several constraints: The average in-degree must be fixed, and the $w_i$ and $w_{j\to i}$ distributions must be normalized. However we can make the problem unconstrained by using the following transformations, $$\begin{gathered}
  w_i = \frac{e^{\widetilde{w}_i}}{\sum_j e^{\widetilde{w}_j}} \\
  w_{j\to i} = \frac{e^{\widetilde{w}_{j\to i}}}{\sum_l e^{\widetilde{w}_{l\to i}}},\end{gathered}$$ where $\widetilde{w}_i$ and $\widetilde{w}_{j\to i}$ are unconstrained real variables. Likewise we can transform $\lambda_i$ as $$\begin{gathered}
  \widetilde{k}_i = \frac{e^{\widetilde{\lambda}_i}}{\tanh e^{\widetilde{\lambda}_i}} \\
  k_i = c \widetilde{k}_i + \gamma \\
  \lambda_i = k_i \tanh \lambda_i  \label{eq:li},\end{gathered}$$ where $\widetilde{\lambda}_i$ is also an unconstrained real variable. To obtain $\lambda_i$ from Eq. \[eq:li\] it is simply iterated until it converges to the correct value, within some desired precision. The values of $c$ and $\gamma$ are chosen to force $k_i \ge 1$ and the average in-degree to a prescribed value ${{\left<k\right>}}$, $$\begin{aligned}
  c &= \frac{{{\left<k\right>}}}{\sum_i\widetilde{k}_i w_i},& \gamma &= 0,  &&\text{ if } \widetilde{k}_m \ge 1 \\
  c &= \frac{{{\left<k\right>}} - 1}{\sum_i\widetilde{k}_i w_i - \widetilde{k}_m},& \gamma &= 1 - c\widetilde{k}_m,  &&\text{ otherwise, }  
\end{aligned}$$ where $\widetilde{k}_m = \min(\{\widetilde{k}_i\})$. Thus we have obtained an unconstrained minimization problem of function $\mathcal{F}$ with respect to the variables $\{\widetilde{w}_i\}$, $\{\widetilde{w}_{j\to i}\}$, $\{\widetilde{\lambda}_i\}$.

In order to find the minimum of Eq. \[eq:fe\], we employed the L-BFGS quasi-Newton algorithm [@byrd_limited_1995], with the gradient computed by finite differences.

Structural phase transition {#sec:phase}
===========================

The minimization of the free energy leads to one of two possible structures (see Fig. \[fig:trans\]): 1. For low values of $\beta$ the topology is a fully random graph; 2. For larger values of $\beta$ there is the emergence of a *segregated core* structure, where one of the blocks has a larger in-degree density and is more responsible for the regulation of the whole network.

![\[fig:trans\](Color online) Structural phase transition observed when varying the selective pressure $\beta$, as described in the text.](net-a0_5-alt.pdf "fig:"){width="\columnwidth"}\
Random topology

$\xrightarrow{\text{increasing } \beta}$

![\[fig:trans\](Color online) Structural phase transition observed when varying the selective pressure $\beta$, as described in the text.](net-a1-alt.pdf "fig:"){width="\columnwidth"}\
Segregated core

In order to precisely characterize the phase transition, we define the following order parameter, $$\phi = \frac{b^* - b_r}{b_{\text{min}} - b_r},$$ where $b_r$ is the value of $b^*$ for a fully random network, and $b_{\text{min}}$ is the smallest possible value of $b^*$ for a given ${{\left<k\right>}}$ [@peixoto_behavior_2012], given by $$b_{\text{min}} = \sum_k p_k m_k(P),$$ where $p_k$ is given by Eq. \[eq:pk\] with $k_i = {{\left<k\right>}}$. We have therefore that $\phi \in [0, 1]$ and if $\phi = 0$ the network ensemble must be fully random, and if $\phi = 1$ it has the largest possible value of fitness.

As can be seen in Fig. \[fig:diag-op\], there is a second-order phase transition at a critical value $\beta^*$, which depends on the noise level $P$. The dependence of $\beta^*$ on $P$ divides the $\beta\times
P$ phase diagram into an upper and lower branch, as can be seen in Fig. \[fig:diag-phi\]. The branches are divided at a value of $P=P^*_r$, which is the critical value of noise of a fully random network (see [@peixoto_behavior_2012] for an exact calculation of $P^*_r$). At this value of noise, a random network undergoes a dynamic phase transition, where the steady state error fraction reaches the maximum level $b^*=1/2$, and the dynamics become ergodic, as was described previously. For $ P < P^*_r$, random networks are intrinsically robust, since $b^* < 1/2$, and the critical value $\beta^*$ becomes larger with smaller $P$. In other words, the smaller is the value of noise $P$, the better is the behavior of fully random networks, such that the entropic cost of providing further improvement by creating a segregated core becomes larger, which therefore occurs only at larger values of selective pressure. The situation changes for $P>P^*_r$. Since random networks are no longer resilient, and have collapsed onto $b^*=1/2$ (see Fig. \[fig:diag-pc\]), a segregated core provides a significant improvement, for a relatively low entropic cost. This cost increases with $P$, since the core needs to be either denser, smaller or more isolated to provide the same benefit under larger noise. Thus the value of $\beta^*$ also increases with $P$. Interestingly, in the vicinity of $P=P^*_r$, the value of $\beta^*$ tends to zero. For this value of noise, the dynamics of the fully random topology lies exactly at the critical point where $b^*=1/2$, and even the smallest (structural) perturbation can move the fixed point appreciably. Since such small structural perturbations have negligible entropic cost, the value of $\beta^*$ vanishes to zero. Thus, networks with ${{\left<k\right>}}$ such that $P^*_r=P$ are the most easily evolvable.

![\[fig:diag-op\] (Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$, for different noise levels $P$, and for ${{\left<k\right>}} = 5$. The left panel shows curves for $P < P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panel shows curves for $P > P^*_r$. The curves on the left panel are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-op-k5-o2.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:diag-op\] (Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$, for different noise levels $P$, and for ${{\left<k\right>}} = 5$. The left panel shows curves for $P < P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panel shows curves for $P > P^*_r$. The curves on the left panel are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-op-k5-o1.pdf "fig:"){width="0.49\columnwidth"}

![\[fig:diag-phi\](Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$ and noise $P$, for different values of ${{\left<k\right>}}$.](diag-phi-k3.pdf "fig:"){width="1\columnwidth"}\
[${{\left<k\right>}} = 3$]{}

![\[fig:diag-phi\](Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$ and noise $P$, for different values of ${{\left<k\right>}}$.](diag-phi-k5.pdf "fig:"){width="1\columnwidth"}\
[${{\left<k\right>}} = 5$]{}

\

![\[fig:diag-phi\](Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$ and noise $P$, for different values of ${{\left<k\right>}}$.](diag-phi-k7.pdf "fig:"){width="1\columnwidth"}\
[${{\left<k\right>}} = 7$]{}

![\[fig:diag-phi\](Color online) The order parameter $\phi$ as a function of the selective pressure $\beta$ and noise $P$, for different values of ${{\left<k\right>}}$.](diag-phi-k10.pdf "fig:"){width="1\columnwidth"}\
[${{\left<k\right>}} = 10$]{}

The value of $b^*$ itself can be seen in Fig. \[fig:diag-pc\]. The upper branch $P>P^*_r$ corresponds to transitions from $b^*=1/2$ (ergodic dynamics) to $b^*<1/2$ (nonergodic dynamics), whereas the lower branch $P<P^*_r$ shows $b^*<1/2$ for both phases.

![\[fig:diag-pc\](Color online) *Left:* Value of the steady-state average error $b^*$ as a function of the selective pressure $\beta$ and noise $P$, for ${{\left<k\right>}}=5$. *Right:* Maximum tolerable noise $P^*$, as a function of the selective pressure $\beta$ and noise $P$, for ${{\left<k\right>}}=5$.](diag-e-k5.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:diag-pc\](Color online) *Left:* Value of the steady-state average error $b^*$ as a function of the selective pressure $\beta$ and noise $P$, for ${{\left<k\right>}}=5$. *Right:* Maximum tolerable noise $P^*$, as a function of the selective pressure $\beta$ and noise $P$, for ${{\left<k\right>}}=5$.](diag-pc-k5.pdf "fig:"){width="0.49\columnwidth"}

The topological properties of each phase can be seen in detail in Fig. \[fig:diag-top\], where are shown the average in-degrees $\{k_i\}$, block sizes $\{w_i\}$ and the total fraction of inputs which lead to the segregated core, $E_c=\sum_j w_{c\to j} w_j k_j / {{\left<k\right>}}$, where $c$ is the core block. The core block emerges at $\beta =
\beta^*$, with an infinitesimal size, but with a value of $k_i$ which is usually significantly above average. For $P>P^*_r$ the segregated core usually has a significantly larger average in-degree than for $P<P^*_r$. The dominance and segregation of the core block increases continuously with $\beta$, reaching values close to $E_c=1$, for larger values of $\beta$, especially for values of $P>P^*_r$.

![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-ki-k5-o2.pdf "fig:"){width="49.00000%"} ![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-ki-k5-o1.pdf "fig:"){width="49.00000%"}\
![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-wi-k5-o2.pdf "fig:"){width="49.00000%"} ![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-wi-k5-o1.pdf "fig:"){width="49.00000%"}\
![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-wij-core-k5-o2.pdf "fig:"){width="49.00000%"} ![\[fig:diag-top\](Color online) Block average in-degrees $\{k_i\}$, block sizes $\{w_i\}$, and total fraction $E_c$ of inputs originating from the core block, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$. The left panels show curves for $P <
    P^*_r$, where $P^*_r$ is the critical value of noise for a fully random network, and the right panels show curves for $P >
    P^*_r$. All curves on the left panels are shown in order of increasing $P$ from right to left, and on the right panel, from left to right.](diag-wij-core-k5-o1.pdf "fig:"){width="49.00000%"}

Each network on the evolved ensemble has a critical value of noise $P^*$ (different from the value of $P$ for which it was evolved), for which its dynamics undergoes the aforementioned ergodicity transition and which represents the maximum tolerable noise (see [@peixoto_behavior_2012] for an exact calculation of $P^*$ for arbitrary blockmodels). Interestingly, the evolution of $b^*$ does not automatically result in larger values of $P^*$, as is shown in Fig. \[fig:diag-pc\]: Some ensembles evolved under larger selective pressure possess a lower value of $P^*$ than others evolved under lower selective pressure (for the same value of $P$ under evolution). This means the evolution is reasonably specialized for the level of noise it is under, and the behavior of the networks under larger values of noise for which they were evolved is not automatically better than that of other networks with smaller fitness. However, despite these deviations, the general tendency is that, for larger values of $\beta$, $b^*$ and $P^*$ are decreased and increased, respectively.

Monte-Carlo simulations {#sec:mc}
=======================

We have also performed Monte Carlo simulations to observe the phase-transition obtained in the previous section. We employed the Metropolis-Hastings [@metropolis_equation_1953; @hastings_monte_1970] algorithm, starting from a random network with $N$ vertices, with a given average in-degree ${{\left<k\right>}}$ and a partition into $n$ blocks, represented by assigning block labels to each vertex (which is initially randomly chosen). At each iteration, one of the following moves is attempted with equal probability:

1.  *Block label move*: A random vertex $v$ is chosen, and its block label is randomly chosen among all $n$ possible values. \[it:label\]

2.  *Input move*: A vertex $v$ is chosen with probability $p\propto k(k-1)$, where $k$ is the in-degree of $v$. Another vertex $u$ is randomly chosen with uniform probability. Two random inputs from $v$ are deleted and moved to $u$. \[it:input\]

3.  *Source move*: A random vertex $v$ is chosen. A random input from $v$ is deleted and replaced by a randomly chosen one.

A move is rejected if it generates parallel edges or self-loops. The difference $\Delta b^*$ of the value of $b^*$ after and before the move is computed. The move is then accepted with a probability $p_a$ given by $$p_a =
  \begin{cases}
    1 &\text{ if } \Delta b^* \le 0,\\
    e^{-\beta N\Delta b^*} &\text{ otherwise.}
  \end{cases}$$ The probability $p\propto k(k-1)$ in move (\[it:input\]) is chosen to correspond to *two* independent single-edge moves affecting the same vertices $v$ and $u$, where in each move a random edge is chosen, and its target is moved to a randomly chosen vertex. This guarantees that there is no topological bias, and that the in-degrees are always odd.

![\[fig:mc-trans\](Color online) Average order parameter ${{\left<\phi\right>}}$ and steady-state error level $b^*$ as functions of the selective pressure $\beta$, for different values of noise $P$ and ${{\left<k\right>}}=5$, obtained with Monte Carlo simulations, for network sizes shown in the legend. On the right plot, the red star symbols ($\star$) correspond to empirical values of $b^*$ as obtained with Eq. \[eq:bn\_dyn\]. The solid gray lines are theoretical values obtained by minimizing the free energy.](monte-carlo-phi-k5-ns2.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:mc-trans\](Color online) Average order parameter ${{\left<\phi\right>}}$ and steady-state error level $b^*$ as functions of the selective pressure $\beta$, for different values of noise $P$ and ${{\left<k\right>}}=5$, obtained with Monte Carlo simulations, for network sizes shown in the legend. On the right plot, the red star symbols ($\star$) correspond to empirical values of $b^*$ as obtained with Eq. \[eq:bn\_dyn\]. The solid gray lines are theoretical values obtained by minimizing the free energy.](monte-carlo-es-k5-ns2.pdf "fig:"){width="0.49\columnwidth"}

![\[fig:mc-top\](Color online) Block average in-degrees $\{k_i\}$ and sizes $\{w_i\}$, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$, obtained with Monte-Carlo simulations, for network sizes shown in the legend. The solid gray lines are theoretical values obtained by minimizing the free energy.](monte-carlo-ks-k5-ns2.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:mc-top\](Color online) Block average in-degrees $\{k_i\}$ and sizes $\{w_i\}$, as functions of selective pressure $\beta$, for ${{\left<k\right>}} = 5$, obtained with Monte-Carlo simulations, for network sizes shown in the legend. The solid gray lines are theoretical values obtained by minimizing the free energy.](monte-carlo-ws-k5-ns2.pdf "fig:"){width="0.49\columnwidth"}

![\[fig:mc-multiblock\](Color online) *Left:* Block average in-degrees $\{k_i\}$, as a function of selective pressure $\beta$, for ${{\left<k\right>}} = 5$, and obtained with Monte-Carlo simulations, for $N=10000$ and different number of blocks. The solid gray lines are theoretical values obtained by minimizing the free energy. *Right:* Evolution in time of the average in-degrees $\{k_i\}$ in a Monte-Carlo simulation with $n=20$ blocks, $N=10000$ and $\beta=10^3$, showing the eventual merging into only two blocks.](monte-carlo-ks-k5-ns2-3.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:mc-multiblock\](Color online) *Left:* Block average in-degrees $\{k_i\}$, as a function of selective pressure $\beta$, for ${{\left<k\right>}} = 5$, and obtained with Monte-Carlo simulations, for $N=10000$ and different number of blocks. The solid gray lines are theoretical values obtained by minimizing the free energy. *Right:* Evolution in time of the average in-degrees $\{k_i\}$ in a Monte-Carlo simulation with $n=20$ blocks, $N=10000$ and $\beta=10^3$, showing the eventual merging into only two blocks.](monte-carlo-ks-k5-ns20-merge.pdf "fig:"){width="0.49\columnwidth"}

The value of $b^*$ is computed by obtaining the values of $\{w_i\}$, $\{w_{j\to i}\}$ and $\{k_i\}$, and iterating Eq. \[eq:dynblock\]. This is much faster than actually measuring the error level on the network, and produces deterministic values [^6].

Since we have employed the block label move (\[it:label\]), which tends to partition the network evenly into $n$ blocks of equal sizes, we have included an entropic cost associated with the size of a block, which did not exist in the original blockmodel above. In the original model, the partitions themselves are not relevant, and only the resulting graph topology contributes to the entropy. However, move (\[it:label\]) makes the algorithm very efficient and easy to implement, and it should not fundamentally change the results. But in order to compare with the theory, we need to include the following correction in the number of possible networks: $$\Omega' = \Omega \times \frac{N!}{\prod_i(N w_i)!}$$ which leads to the slightly modified entropy, $$\label{eq:part}
  \mathcal{S}'/N =   \mathcal{S}/N - \sum_iw_i\ln w_i.$$

In Fig. \[fig:mc-trans\] we can see the same phase transition observed previously, which matches very well the theoretical predictions. In Fig. \[fig:mc-top\] the topology can be assessed more closely, and the emergence of the segregated core is clear. Due to the partition entropy introduced in Eq. \[eq:part\], the core does not vanish at the transition; it merges continuously with the other block instead. However, the critical value $\beta^*$ is identical with the non-modified model.

The inclusion of the partition entropy also introduces the fact that different solutions are obtained for a different number of blocks, since this has a direct effect on the preferred sizes of the blocks (see Fig. \[fig:mc-multiblock\], left). However, this *does not* change the fact that for any number of blocks the preferred topology will always be an effective two-block structure. This follows from the argumentation presented previously based on the reduction of entropy resulting from block splits, and can be observed in simulations with many blocks, as shown in Fig. \[fig:mc-top\], which shows a comparison between the topologies obtained with two and three blocks, as well as the outcome of a typical simulation with $20$ blocks, which shows the eventual collapse of into an effective two-block structure.

Gene regulatory networks {#sec:real}
========================

Here we make a comparison with some features observed in actual gene regulatory networks. We consider the networks for *Saccharomyces cerevisiae* (yeast) and *Escherichia coli*, extracted from the YEASTRACT [@abdulrehman_yeastract:_2010] and RegulonDB [@gama-castro_regulondb_2010] databases, respectively. We are interested in extracting the “functional core” of the network, i.e. those nodes which are solely responsible for global regulation, like those belonging to the *segregated core* which emerges in the phase transition observed in the evolutionary process above. We will characterize the core nodes in two ways: 1. Nodes which have an out-degree larger than zero; 2. Nodes which belong to a *strongly connected component* (SCC) of the graph (i.e. the maximal set of nodes which can directly or indirectly regulate each other). The first criterion is a necessary condition, since if the out-degree is zero, then a node is not a regulator. The second criterion is stronger, since even if a node is a regulator, it can have its dynamics completely enslaved by other nodes. The nodes in the SCC are exactly those which are not necessarily enslaved, and can mutually regulate each other. Without a least one SCC in the network, an autonomous behavior with dynamical attractors other than simple fixed-points is not possible.

The yeast network is composed of $N=6402$ nodes, with an average in-degree of ${{\left<k\right>}} \approx 7.51$. The *E. coli* network is smaller and sparser, with $N=1658$ and ${{\left<k\right>}} \approx 2.32$. In both networks the number of transcription factor (TF) genes is much smaller than the total number: $N_{\text{TF}} = 182$ in yeast, and $N_{\text{TF}} = 154$ in *E. coli*. These are core genes according to the first criterion, since they have an out-degree larger than zero, as can be seen in Fig. \[fig:gene\].

\[tb\] ![\[fig:gene\](Color online) Gene regulatory networks for *Saccharomyces cerevisiae* (left) and *Escherichia coli* (right), extracted from the YEASTRACT [@abdulrehman_yeastract:_2010] and RegulonDB [@gama-castro_regulondb_2010] databases, respectively. The nodes in purple (towards the middle) are transcription factor genes, and are the only ones with out-degree larger than zero.](yeast.pdf "fig:"){width="0.49\columnwidth"} ![\[fig:gene\](Color online) Gene regulatory networks for *Saccharomyces cerevisiae* (left) and *Escherichia coli* (right), extracted from the YEASTRACT [@abdulrehman_yeastract:_2010] and RegulonDB [@gama-castro_regulondb_2010] databases, respectively. The nodes in purple (towards the middle) are transcription factor genes, and are the only ones with out-degree larger than zero.](ecoli.pdf "fig:"){width="0.49\columnwidth"}

In yeast, the average in-degree of the core nodes is higher than average, ${{\left<k\right>}}_c \approx 10.03$, as observed in the segregated core phase of the evolved networks obtained. For the SCC, the number of nodes decreases slightly to $N_{\text{SCC}}=146$, and the average in-degree changes negligibly ${{\left<k\right>}}_{CC} \approx 10.48$ (if one counts only edges between vertices of the SCC, this value is virtually identical, ${{\left<k\right>}}_{CC} \approx 10.42$). This is similar to what was found previously in [@maslov_computational_2005] for the yeast network (using an older, and less complete dataset with only $837$ genes). They have also found that the TF genes have different connection patterns, and those with the largest out-degree tend to regulate genes with lower than average in-degree. However they did not find that the TF genes form a *denser* subgraph, with an larger than average in-degree, which is most likely due to the incompleteness of the dataset used. Very similar numbers to those presented here were obtained more recently in [@balaji_comprehensive_2006], using a more complete dataset (which is not identical to the one used in this work).

For *E. coli* the situation changes somewhat: The average in-degree of the transcription factor nodes is ${{\left<k\right>}}_c \approx 1.97$, which is in fact *lower* than the global average. However, if one extracts the largest SCC, the number of nodes drops dramatically to $N_{\text{SCC}}=8$. These nodes are responsible for the regulation of $411$ genes. A majority of $1093$ genes are instead enslaved to the dynamics of SCC with only two mutually regulating nodes. Although the largest SCC *does* have an average in-degree ${{\left<k\right>}}_{CC} = 6$, the core topology seems significantly more sparse than for yeast and the evolved networks, at least with the data currently available [@cosentino_lagomarsino_hierarchy_2007]. Arguably, such a sparse regulating core is suspect from the point of view of data set completeness, since it would mean that the range of dynamical behavior for the regulatory network is very restricted. As previously mentioned, and older and less complete data set for yeast also did not reveal a denser regulating core [@maslov_computational_2005]. Nevertheless, one should also consider that such real networks are simultaneously under different, possibly competing selective pressures which also influence the resulting topology, robustness against noise being only one of them. These other factors, which are neglected in the model, could be one reason for such a discrepancy. We emphasize, however, that although apparently it is not denser, a regulating core certainly exists in the measured *E. coli* network, which is at least in partial qualitative agreement with what is observed in the model.

There are other factors that may contribute to this observed segregation which are not in principle related to noise resilience. For instance, non-regulating genes exist mostly to transcribe proteins which have some specific metabolic or structural function in the cell, and it may be difficult for these proteins to have a dual role as transcription factors, and therefore become specialized (although non-specialized proteins are not impossible, since a protein can in principle bind both to DNA and to other proteins). Nevertheless, there are good reasons to consider robustness to noise as a very plausible driving force toward this type of topology. This is corroborated, for instance, by evidence that core TF genes tend to be less noisy [@fraser_noise_2004; @jothi_genomic_2009], and that the vast majority of TF genes in yeast are not vital for the survival of the cell if repressed in isolation [@balaji_uncovering_2006]. This is fully compatible with the idea of a highly redundant functional core, which provides robustness for the rest of the network.

Another feature which is commonly investigated in empirical networks is the in- and out-degree distributions. The in-degree distribution is often narrow, while the out-degree distribution is broader, and as some suggest [@maslov_computational_2005], compatible with a power law. The model considered in this work is parametrized as a stochastic blockmodel, where each block has in- and out-degrees that are Poisson distributed. When the segregated core emerges, the system is composed of only two blocks, thus both the in- and out-degree distributions are bimodal. The in-degree distribution is indeed narrower, since the difference between the average in-degree of the two blocks is not very large for most networks obtained. The out-degree distribution is also much broader, since the average out-degree of the non-core block tends to zero, while for the core block it tends very rapidly to infinity, when the selective pressure is increased. However, the homogeneous and seemingly scale-free properties of the empirical distributions are not reproduced by the model. This implies that these features are not a direct result of evolved robustness against noise, and may, for instance, be due simply to mutational bias caused by gene duplication, which is known to qualitatively reproduce these types of in- and out-degree distributions [@ispolatov_duplication-divergence_2005; @enemark_gene_2007].

Conclusion {#sec:conclusion}
==========

We have investigated the effect of selective pressure favoring robustness against noise on the structural evolution of Boolean networks with optimal majority functions, functioning as a conceptual model for gene regulation. We have mapped the evolutionary process onto a Gibbs ensemble, and obtained its outcome by minimizing the associated free energy. We showed that the structural properties of the system undergo a phase transition at a critical value of selective pressure, from a random topology to a segregated-core structure, where a smaller fraction of the nodes form an isolated core, which is denser than the rest of the network and is responsible for most of the regulation. Since the core is denser, its nodes can profit from more regulatory redundancy, which greatly diminishes the effect of noise. This robustness is propagated to the rest of the network, which relies on the core for most of the regulation. The segregated core becomes denser, smaller, and more isolated as the selective pressure increases. We have compared the theoretical predictions with Monte-Carlo simulations of actual networks, and found perfect agreement.

We have also shown that this segregated-core structure is present in the gene regulatory network of yeast and *E. coli*. Both networks are composed of a much smaller fraction of transcription factor genes which are responsible for all regulation. In yeast, the existing core structure is very similar qualitatively to the outcome of the evolutionary process considered, with transcription factor genes forming a denser subgraph, with an average in-degree above the average for the whole network. In *E. coli* the isolated transcription factor core is composed of few very small regulating cores (strongly connected components), the largest of which has only eight nodes. We conjecture that such a sparse regulating core is possibly due to data set incompleteness, since it would severely restrict the range of possible dynamical behavior for the network. A less complete data set for yeast also did not show a denser regulating core structure [@maslov_computational_2005], although it is clearly seen with more up-to-date datasets including more genes and interactions [@balaji_comprehensive_2006]. However, one should not rule out other selection criteria which are not incorporated in the model.

It should also be noted that regulating cores of transcription factors are a common feature of other organisms, such as *Mycobacterium tuberculosis* [@balazsi_temporal_2008]. Additionally, a similar (but not identical) “bow-tie” structure was also observed in the mammalian signal transduction network [@maayan_formation_2005; @maayan_insights_2009], where most pathways are funneled through a central core.

It is possible to formulate other reasons for the existence of such a core structure, such as a forced specialization of genes into either transcription factors or target genes. Furthermore one should mention that the effects of noise are not always detrimental, and can in some circumstances even be beneficial [@blake_phenotypic_2006; @eldar_functional_2010]. Nevertheless there is compelling evidence that the core genes provide a degree of robustness to the cell. Not only are the best-connected TF nodes less noisy [@fraser_noise_2004; @jothi_genomic_2009], they are usually found — if removed individually — not to be vital for cell survival [@balaji_uncovering_2006]. This corroborates the idea that one of the major functions of the regulating core is to provide robustness via redundancy.

Furthermore, aside from the direct applicability to gene regulation, we have identified a fundamental mechanism of robustness against noise, which emerges naturally when networks are randomly selected for that purpose [^7]. Although most interesting systems require more than just robustness for their functioning, it is reasonable to conclude that the emergence of regulating cores is to be expected when there is enough selective pressure favoring noise resilience.

[^1]: Another realistic source of noise are perturbations in the update sequence of nodes, since gene regulation lacks a global synchronizing clock [@klemm_topology_2005]. However, it can be shown that absolute robustness against this type of noise can be achieved in an independent manner, and with a very small effect to the global topological characteristics of the system [@peixoto_boolean_2009].

[^2]: The definition above will lead to a bias if $k_i$ is an even number, since if the sum happens to be exactly $k_i / 2$ the output will be $0$, arbitrarily. Alternative definitions could be used, which would remove the bias [@szejka_phase_2008]. However, for the analysis presented here, this is not an issue since $k_i$ is always odd.

[^3]: Blockmodels are essentially equivalent to the hidden-variable model [@boguna_class_2003], when the hidden variables are discrete, and their multiplicity is smaller than the number of nodes.

[^4]: This is not the only two-block structure which can generate arbitrary values of $b^*$. In [@peixoto_behavior_2012] it is shown how a bipartite “restoration” structure also achieves this, albeit less efficiently.

[^5]: We have empirically verified this by minimizing the free energy with up to 10 blocks, and the solutions were always *identical* to that of the two-block case presented in the following section. The general character of the two-block topology was also verified by Monte Carlo simulations with up to $20$ blocks, as discussed later in the text (see also Fig. \[fig:mc-multiblock\]).

[^6]: One could argue that this may overlook the buildup of correlations in the network, since it assumes that the blocks are homogeneous and have an random in-degree distribution as given by Eq. \[eq:pk\]. However, we are *only* interested in networks which have this property, which, as discussed in the text, correspond to a partial maximization of entropy when the remaining constraints are in place, so this should not be an issue. To be sure, we have compared the value of $b^*$ computed this way with the actual empirical value and found a very good agreement (see Fig. \[fig:mc-trans\]).

[^7]: Other evolutionary, non-equilibrium pathways are also possible, see e.g. [@perotti_emergent_2009].
