Sam F Greenbury

Iain G Johnston

Matthew A Smith

Jonathan P K Doye

Ard A Louis

The effect of scale-free topology on the robustness and evolvability of genetic regulatory networks

Keywords: gene regulatory networks, Boolean threshold networks, robustness, evolvability, scale-free, oscillatory gene expression

, is compared to the ER topology in further detail through an extension to multiple target outputs, with either an oscillatory or a non-oscillatory nature. For multiple oscillatory targets of the same length, the differences between SF out and ER networks are enhanced, but for non-oscillatory targets both types of networks show fairly similar evolvability. We find that SF networks generate oscillations much more easily than ER networks do, and this may explain why SF networks are more evolvable than ER networks are for oscillatory phenotypes. In spite of their greater evolvability, we find that networks with SF out topologies are also more robust to mutations than ER networks. Furthermore, the SF out topologies are more robust to changes in initial conditions (environmental robustness). For both topologies, we find that once a population of networks has reached the target state, further neutral evolution can lead to an increase in both the mutational robustness and the environmental robustness to changes in initial conditions.

Introduction

Robustness & evolvability in genetic regulatory networks

Biological systems require their phenotype to be robust to a variety of perturbations. They must be mutationally robust to minimize the possibility of a permanent deleterious mutation in the genome. They must be environmentally robust to mitigate the effects of temporary macroscopic environmental variation, or stochastic noise in microscopic biochemical processes. By possessing a robust phenotype, functions required for survival can be performed consistently and descendants can expect a similarly reliable genotype [START_REF] Wagner | Robustness and Evolvability in Living Systems[END_REF].

At the same time, environments can change permanently and populations must adapt by producing heritable phenotypic variation. This ability to generate phenotypic innovation from genetic changes under natural selection is often called evolvability, and is observed throughout the biological world [START_REF] Dawkins | The evolution of evolvability[END_REF][START_REF] Kirschner | Evolvability[END_REF][START_REF] Wagner | Robustness and Evolvability in Living Systems[END_REF]. It may itself be an evolvable trait [START_REF] Earl | Evolvability is a selectable trait[END_REF]. The relationship between environmental robustness and evolvability is very complex [START_REF] West-Eberhard | Developmental plasticity and evolution[END_REF][START_REF] Kirschner | The plausibility of life: Resolving Darwin's dilemma[END_REF]. On the other hand, it would seem at first sight that evolvability and mutational robustness are straightforwardly antagonistic: the more easily mutations change the properties of an organism, the less mutationally robust it should be.

However, this simple picture can be challenged on numerous fronts (Bloom et al., 2007a,b;Wagner, 2008b,a;[START_REF] Daniels | Sloppiness, robustness, and evolvability in systems biology[END_REF][START_REF] Draghi | Mutational robustness can facilitate adaptation[END_REF]. It is important to take into account the fact that evolution acts on populations and not on individuals. Furthermore, the structure of the evolutionary landscape, the population size N pop , the mutation rate μ and the type of phenotypic variation required are all factors that can influence the relationship between mutational robustness and evolvability.

For example Wagner, using RNA secondary structure as a model (Wagner, 2008b), demonstrated the importance of distinguishing between genotype (sequence) robustness and evolvability, which do share an antagonistic relationship, and phenotype (structure) robustness and evolvability which, by contrast, share a beneficial relationship. For a single sequence (genotype), mutational robustness is simply the fraction of mutations that alter the phenotype. However, a single RNA phenotype (a particular secondary structure) may be the minimal free energy structure for a large number of different genotypes. Such a set of genotypes is called a neutral space, and during evolution, a population can spread on the neutral space to generate a large amount of genetic diversity. The larger the neutral space, the more likely that random mutations will generate another member of the neutral space, so the larger the mutational robustness. It can be demonstrated (Wagner, 2008b) that for RNA, larger neutral spaces have a larger number of different phenotypes accessible with a single mutation. Thus a population that can generate greater genetic diversity through non-deleterious mutations across a larger neutral space, also has more evolutionary innovation available, making such phenotypes more evolvable. Furthermore, for RNA, the utility of the neutral space with respect to robustness and evolvability is shown to be independent of the evolutionary dynamical regime (being of importance in both the N pop μ 1 and N pop μ 1 cases, where N pop is population size and μ is mutation rate).

These arguments depend on the structure of the neutral space (Wagner, 2008a). If it is made up of many disconnected pockets, then a population may not be able to explore the entire neutral space. If on the other hand all the neutral neighbours are clustered together tightly in one part of the space, then they may not have a large diversity of other phenotypes within easy reach. Models of RNA secondary structure have the advantage that they can be relatively easily analyzed by computational methods [START_REF] Wagner | Robustness and Evolvability in Living Systems[END_REF]. Thus detailed questions about the topology of neutral spaces and the effects of changing mutation rates and population sizes can be analyzed in some detail. For other systems, however, the relationship between mutational robustness, neutral spaces and evolvability may be more complex [START_REF] Draghi | Mutational robustness can facilitate adaptation[END_REF].

In this paper we focus on gene regulatory networks (GRNs), another set of biological systems that can be analyzed computationally. There has been an enormous amount of interest in modeling GRNs because they play a central role in biological functionality [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. Moreover, it is increasingly being recognized that altering gene expression patterns is a highly effective evolutionary mechanism. Important phenotypic changes are often achievable through very small numbers of mutations to regulatory regions [START_REF] Carroll | Evolution at two levels: On genes and form[END_REF]. For these reasons, the mutational robustness, environmental robustness and evolvability of GRNs have received much theoretical attention. For example, the reliability under perturbed system configurations (environmental robustness) has been investigated [START_REF] Aldana | Boolean dynamics of networks with scale-free topology[END_REF][START_REF] Klemm | Topology of biological networks and reliability of information processing[END_REF] whilst the evolution of this type of robustness is considered in [START_REF] Braunewell | Reliability of regulatory networks and its evolution[END_REF]. Robustness to mutation has also been shown to be an evolvable property, correlating with robustness to initial gene expression states (Ciliberti et al., 2007b). Furthermore, there is evidence suggesting that neutral spaces can enhance evolvability of expression patterns (Ciliberti et al., 2007a). The effect of varying environments has been demonstrated to increase evolvability [START_REF] Crombach | Evolution of evolvability in gene regulatory networks[END_REF], whilst the role of overall network topology has also been considered [START_REF] Aldana | Boolean dynamics of networks with scale-free topology[END_REF][START_REF] Oikonomou | Effects of topology on network evolution[END_REF][START_REF] Aldana | Robustness and evolvability in genetic regulatory networks[END_REF].

Modelling GRNs

Realistically modelling large complex interacting GRNs is a difficult task. Even a modest sized network can have an enormous number of parameters depending on which details of the fundamental processes involved in transcription and translation are taken into account [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. In practice, therefore, many different approaches to simplify the dynamics of GRNs can be found in the literature, each with their strengths and limitations [START_REF] Polynikis | Comparing different ode modelling approaches for gene regulatory networks[END_REF].

A particularly simplified model, the Boolean network, was proposed by [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF]. The generic Boolean model makes two important assumptions. First, the gene states are the variables under time evolution, as opposed to the concentrations of gene products. These states are further assumed to be either active or inactive, i.e. either "on" or "off". This approximation is based upon evidence that the changes in gene state are cooperative transitions, allowing the approximation of sigmoidal functions as step-wise ones [START_REF] Kauffman | The Origins of Order: Self-Organization and Selection in Evolution[END_REF]. The second assumption is synchronous updating of gene states throughout the network in discrete time steps through physically reasonable rules. In the original model, these rules were random Boolean functions of all incoming nodes, whose output determined a given node's state in the next time step. Other update rules, based upon neural threshold network models [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF], have also been applied in Boolean threshold networks [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF]. The Boolean model has been extensively used in theoretical biology and, despite its coarse-grained nature, has been successfully applied, for example, to the yeast cell cycle network [START_REF] Li | The yeast cell-cycle network is robustly designed[END_REF], the segment polarity gene network of D. melanogaster [START_REF] Albert | The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster[END_REF] and the floral cell fate of A. thaliana [START_REF] Espinosa-Soto | A gene regulatory network model for cell-fate determination during arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles[END_REF].

An important feature of Boolean network models is the manifestation of attractors. Following a certain amount of time evolution of the system, all initial configurations converge on a very small subset of possible configurations -the attractors of the network, also referred to as the network's attractor landscape [START_REF] Kauffman | The Origins of Order: Self-Organization and Selection in Evolution[END_REF]. Boolean network models can exhibit ordered, chaotic and critical regimes [START_REF] Kauffman | The Origins of Order: Self-Organization and Selection in Evolution[END_REF]. In the ordered regime, small perturbations in an initial configuration tend not to spread, with the same attractor reached as in the unperturbed case [START_REF] Aldana | Boolean dynamics of networks with scale-free topology[END_REF]. However, in the chaotic regime these small perturbations spread, possibly resulting in the network arriving at a different attractor. The transition from the ordered to the chaotic regime occurs as the average connectivity increases. The critical regime lies on the boundary between the ordered and chaotic regimes.

It has been hypothesised that the attractor landscape represents the possible "cell fates" given a cell's regulatory interactions. Different sets of initial configurations induce the production of different attractors, regarded as different cell types [START_REF] Kauffman | The Origins of Order: Self-Organization and Selection in Evolution[END_REF]. Recently, some evidence has been produced supporting this theory [START_REF] Espinosa-Soto | A gene regulatory network model for cell-fate determination during arabidopsis thaliana flower development that is robust and recovers experimental gene expression profiles[END_REF][START_REF] Huang | Cell fates as highdimensional attractor states of a complex gene regulatory network[END_REF]. Within this interpretation, there is a diverse range of possibilities, from "frozen" attractors of unit length, as required in developmental gene expression patterns, to ones including longer oscillatory periods capable of more complex function, for example, the changes in gene expression over the course of the cell cycle [START_REF] Li | The yeast cell-cycle network is robustly designed[END_REF].

The role of topology

Within the Boolean approximation that coarse-grains out many biochemical parameters, a GRN can be viewed as a directed graph, with genes as nodes and gene interactions as edges. The number of incoming and outgoing connections for Power-law Power-law a given node -its incoming/outgoing degree -can be described statistically by distributions that characterize the overall structures or topologies [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF]. If a set of nodes is randomly assigned connections, the resulting degrees are binomially distributed, forming Erdős-Rényi graphs [START_REF] Erdős | On the evolution of random graphs[END_REF]. Another class is the scale-free (SF) topology, where the node degrees are power-law distributed. Given that the incoming and outgoing degrees are distinct, these two distributions can be combined to produce four different topologies, which are defined in Table 1. The probability of any SF topology occuring in a randomly constructed network decreases rapidly as the number of nodes becomes larger than the average connectivity [START_REF] Oikonomou | Effects of topology on network evolution[END_REF].

Since the observation of the Internet's SF topology by [START_REF] Albert | Internet: Diameter of the worldwide web[END_REF], SF networks have been identified in many other fields, ranging from citation networks [START_REF] Barabási | Emergence of scaling in random networks[END_REF] to biological networks, including metabolic networks [START_REF] Jeong | The large-scale organization of metabolic networks[END_REF] and GRNs [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF]. Given the improbability of large SF topologies occurring randomly, this ubiquity has been suggested to be a consequence of common generating mechanisms [START_REF] Barabási | Emergence of scaling in random networks[END_REF]. One possible mechanism is preferential attachment, whereby highly connected nodes have a greater probability of receiving a new connection -"the rich get richer". This could certainly explain the power law distribution of networks such as the Internet [START_REF] Barabási | Emergence of scaling in random networks[END_REF], but may not be the only responsible mechanism [START_REF] Keller | Revisiting "scale-free" networks[END_REF]. A modified view of preferential attachment has been suggested to be responsible for the power laws observed in large GRNs [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF]. A fundamental process in biological evolution and GRN growth is gene duplication [START_REF] Carroll | Evolution at two levels: On genes and form[END_REF], where a gene is simply copied twice into the daughter system. In the corresponding GRN, all connections from the duplicated gene are doubled. As a highly connected gene is more likely to be connected to the randomly duplicated gene, new genes preferentially attach to highly connected genes.

The process outlined above would give rise to an SF degree distribution. However, in GRNs, it is only the outgoing degree that has been observed to possess a power-law distribution [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF][START_REF] Aldana | Robustness and evolvability in genetic regulatory networks[END_REF]. Either evolutionary rewiring has undone SF incoming degrees (due to the vast number of random network configurations relative to scale-free ones), or it is not gene duplication that is responsible for the generation of power-laws within GRNs. In either case, there must be another force either maintaining or generating an outgoing SF distrbution.

The effect of the SF topology on Boolean networks has already been investigated by other authors. For example, [START_REF] Aldana | Boolean dynamics of networks with scale-free topology[END_REF] and [START_REF] Aldana | A natural class of robust networks[END_REF] investigated the effect of dynamical perturbations in Boolean networks with an SF topology. In random networks, the ordered regimes is achievable through all nodes having low connectivities, with fine tuning of mean degree connectivity. For SF topologies, however, the ordered regime can be attained whilst allowing the presence of highly connected nodes, or "hubs". In more recent work, [START_REF] Aldana | Robustness and evolvability in genetic regulatory networks[END_REF] have also demonstrated increased robustness of the attractor landscape in SF topologies under the process of gene duplication. In an important paper, [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] studied large (N = 500 node) networks with Boolean threshold dynamics [START_REF] Kurten | Correspondence between neural threshold networks and Kauffman Boolean cellular automata[END_REF]. They used standard genetic algorithms (GAs) to model evolution and found that when the selection pressure is for a single node in the network to produce a desired oscillatory target, an SF in topology will on average reach its target in significantly less generations than an ER network will, suggesting that the former is more evolvable than the latter. The enhanced evolvability for oscillatory targets holds for a number of different average connectivities. They argue that, in contrast to ER networks where the parameters must be fine-tuned to achieve critical behaviour (and associated enhanced evolvability), the SF network will promote evolution for a wide set of different parameter values.

In this study we consider the work of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] further, with a particular focus on how the enhanced evolvability they find depends on the specific target (phenotype) and how it relates to both mutational and environmental robustness. Besides the SF in topology used in [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] we also investigate the other topologies in Table 1 with a particular focus on the SF out topology that better resembles biological GRNs. We note that in our study, as in [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], the emphasis is on general results for very schematic models of GRNs rather than on networks directly extracted from particular experimental systems. The advantage of this approach is that general trends are easier to extract, while the disadvantage is of course that any direct connections to biological systems will need to be worked out further down the line.

We proceed as follows: In Section 2 we describe in some detail the methods that we use to simulate the evolutionary dynamics of the networks. In Section 3 we study the evolvability of the four different topologies described in Table 1 for single oscillatory targets of length L = 10, confirming that all types of SF topology are more evolvable than an ER topology. We also demonstrate that SF topologies are much more likely to generate oscillatory outputs even without any evolutionary pressure. Section 4 introduces two new kinds of targets. For multiple oscillatory targets of the same length L, the difference between ER and SF out topologies is enhanced for increasing number of independent targets. However, these differences are much smaller for multiple targets of different length. The difference in evolvability between the ER and the SF out topologies is also quite small when evolving towards a target that fixes 100 of the 500 output nodes, suggesting that the enhanced evolvability of the SF out topology depends in part on the kind of phenotype that is being selected for. In Section 5, the mutational and environmental robustness of the ER and SF out topologies are studied. Even though the SF out topologies are much more evolvable than the ER topologies for oscillatory outputs, they are simultaneously more mutationally and environmentally robust. Under neutral evolution both the mutational robustness and the environmental robustness increase. We discuss our main results in Section 6. To show that these results are not limited to the model used in [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], we demonstrate in Appendix A that another GRN model [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF] also exhibits a similar interplay between robustness, evolvability and topology. Appendices B -E are shorter and mainly focus on technical points left out of the main paper.

Methods

We model the evolutionary dynamics of a population of GRNs in several stages. Firstly, sample networks are created, generating a population of genotypes. This is followed by the simulation of network dynamics determining the attractor, with the phenotype obtained from measurements at a node or set of nodes randomly selected at the beginning of the run. This population of networks then undergoes evolution through recurrent mutation of the genotype, followed by selection for the fittest individuals based on their phenotype, which is re-measured every generation.

Structure

Networks were created using the "Configuration model", a standard method of generating arbitrary directed networks [START_REF] Milo | On the uniform generation of random graphs with prescribed degree sequences[END_REF][START_REF] Milo | On the uniform generation of random graphs with prescribed degree sequences[END_REF]. The number of nodes in the networks, N, was set to 500, whilst the average connectivity of network topologies was chosen to be 1.88. These values provided a direct comparison with the work of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], as well as being biologically reasonable [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF].

To allocate node degrees for SF topologies, numbers were drawn randomly from a power-law distribution, with the probability of degree k given by

P S F (k) = k -γ N m=1 m -γ 1 ≤ k ≤ N (1)
where γ parameterises the distribution.

For the ER topology, in the case where the number of nodes is much larger than the average connectivity, the binomial distribution may be approximated as a Poisson distribution [START_REF] Milo | On the uniform generation of random graphs with prescribed degree sequences[END_REF]. In contrast to the SF distribution, where all nodes have at least one connection (P S F (0) = 0), the Poisson distribution permits some nodes to have either an in or out degree of 0, or even both. Such nodes cannot both affect and be affected by the network and thus, their presence introduces a reduction in the effective size of the networks. To combat this effect, we imposed a condition on the Poisson distribution, whereby degrees of 0 are disallowed and the distribution is re-normalised accordingly. We therefore restrict ourselves to considering networks consisting of a single, large connected component. This adapted Poisson distribution for a degree k is given by

P ER (k) = K k e -K (1 -e -K )k! 1 ≤ k ≤ N (2)
where K parameterises the distribution.

The average connectivity of each distribution can be determined simply by k = N k=1 kP(k). The corresponding values of γ and K for an average connectivity of 1.88 are γ = 2.5 and K = 1.431.

Dynamics

A Boolean threshold model was used to model the dynamics. At the beginning of a dynamical run, each node σ i is initialized by a random process to to be either in an "on" state (σ i = 1) or an "off" state (σ i = 0). These are the initial conditions (ICs) of the network. In an evolutionary run, we either impose the same set of initial node states at each generation (constant ICs) or choose initial node states randomly at each generation (random ICs). We investigate both situations in this work (in [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] only random ICs were used). From an initial configuration, the network dynamics are iterated over a set of discrete time steps, with the state of each node updated synchronously between time steps (see Fig. 1a).

The rules for updating the state of a particular node depend on the state of the incoming nodes, and the strength of these signals. An N × N matrix of weights, w, defines the interactions between all nodes. Connections between nodes are assigned weights in matrix w randomly from the set of real numbers on the interval (-1, 1).

The state of node i at the following time step is denoted σ i (t+ 1) and is determined by a sum over all incoming nodes

S i = N j=1 w i j σ j (t) (3)
combined with the following threshold rule:

σ i (t + 1) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 ifS i > 0 σ i (t) if S i = 0 0 ifS i < 0 (4)
To physically interpret the matrix of weights, consider a connection from node i to node j. If w i j > 0, node j promotes node i whilst for w i j < 0, node j inhibits node i. If w i j = 0 there is no connection from j to i.

To prevent frozen dynamics, there is a modified rule for nodes with only a single incoming connection [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. If node j is incoming to node i the the dynamics depend purely on the sign of the weight with Figure 1: A summary of the dynamical and evolutionary processes acting on a network. In (a) an example of the network dynamics is demonstrated. The node coloured blue is selected as the output node. When one time step is performed, the state of the nodes will have been updated synchronously by the rules associated with Eqns. 3 and 4. Further times steps are then performed to produce an output bit string as shown above. This is then compared to the target and its fitness calculated with Eqn. 6. In this particular case maximal fitness is achieved. In (b) the two types of mutation are demonstrated. The red arrows represent a swapping mutation where the grey edges have been switched but preserve each node's degree distribution. The green edge, with a new weight of 0.79, represents a weight mutation from its previous value of 0.23. When a mutation occurs during an evolutionary run, one mutation type is chosen, with each type being equally probable.

σ i (t + 1) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ σ j (t) ifw ji > 0 ¬σ j (t) if w ji < 0, (5)
where ¬ is the NOT operator. The behaviour of nodes connected by this modified rule is simple. Node states are either copied or inverted for positive and negative weights respectively.

The critical parameter values, γ c and K c , separating the ordered and chaotic regimes for each topology, have been previously calculated by [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] based upon the annealed approximation introduced by [START_REF] Derrida | Random networks of automata: A simple annealed approximation[END_REF]. We repeat the calculation (see Appendix B) finding the critical value for the SF topology to be γ c = 2.42 for networks of size N = 500. We also calculate the critical connectivity for the adapted Poisson distribution, which we determine to be K c = 3.538. The primary parameter values used in this study of γ = 2.5 and K = 1.431, both produce average connectivities below that needed for criticality, so that the networks are within the ordered regime.

After a run typically lasting 350 time steps, the steady state dynamics of the output node(s) in each network is determined and this sequence defines the individual's phenotype. An example of a single dynamical time step is performed on model networks in Fig. 1a.

Evolution

To model evolution of the population towards a desired output (phenotype) we used a simple genetic algorithm. First a population of N pop = 50 independent networks is generated. Then, for each generation, the following cutoff selection regime is applied. Each network first produces three daughter networks, enlarging the population from N pop = 50 to N pop = 200. Each node in these daughter networks has a mutation probability μ = 0.02 of undergoing a mutation, a value previously used by [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. There are two types of mutation: either a random change in the weight of any of a node's connections or an incoming connection is swapped with another incoming connection in the network, provided the new connections do not already exist (Fig. 1). Both mutation types preserve the node degree distribution throughout an evolutionary run.

In single target experiments, the output of a node, randomly selected at the start of the evolutionary run, is used to calculate the network's fitness. The fitness is calculated from the minimum Hamming distance between the node's periodic output and the target output over all cyclic permutations. We define this mathematically as

F = max j∈{0,...,δ} ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 - 1 L o L tgt LoLtgt i=1 |σ o (i + j) -σ tgt (i)| ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (6)
where L o is the output period, L tgt is the target period, σ o (i) is the ith bit in the output sequence and σ tgt (i) is the ith bit in the target sequence. The sequences are, for comparison, extended by repeating them periodically to construct sequences of equal length L o L tgt . The sum is then a Hamming distance of these two strings. Taking the maximum over an offset j between 0 and δ = max{L o , L tgt } performs the Hamming distance calculation over all cyclic permutations of the strings. The Hamming distance measures the fraction of the sequence that matches the target, providing a linear measure for the fitness. Once all networks in the enlarged population have been assigned a fitness, the networks are ordered by fitness and the top 50 kept to form the population for the next generation.

Evolution towards single L=10 targets

All types of SF topology are more evolvable

Expanding on the work of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], we tested evolutionary performance or evolvability of each of the four different topologies (ER, SF out , SF in and SF). Evolvability is measured through adaptation speed towards the target output, i.e. the number of generations before every member of a population has reached maximum fitness. Evolvability has been frequently used in different contexts, sometimes creating confusion [START_REF] Poole | Prokaryote and eukaryote evolvability[END_REF][START_REF] Wagner | Robustness and Evolvability in Living Systems[END_REF]. In this case, however, the adaptation speed would seem the most natural measure of evolvability.

A randomly chosen node from each network was evolved towards a randomly chosen L = 10 target for a maximum of 10 000 generations. This was repeated 50 times for each topology, with both constant and random ICs. When evolution is performed with random ICs, robustness to this stochasticity must evolve for an individual to maintain itself in the population. For constant ICs there is no such constraint.

We observed that occasionally an evolutionary run continues for a very large number of generations before converging to a solution, or sometimes not converging at all within the limits of our calculations. The mean of the adaptation time distribution can be dominated by these rare events. For that reason, we use the median adaptation time, arguing that this provides a more accurate measure of the typical adaptation time than the mean does. We define τ to be the number of generations for all members of a population to reach maximum fitness in a single evolutionary run, whilst τ is the median over a batch of runs. The values of τ for the different topologies and ICs, are presented in Table 2.

Out of the three SF topologies, the SF both evolves more rapidly than the SF in topology does, which in turn is followed by SF out and ER. This result supports the findings of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], where the SF in topology is tested against the ER topology and shown to evolve much more rapidly. The order of the topologies, with respect to adaptation speed, is unchanged between constant and random ICs, although the decrease in speed from constant to random ICs varies amongst the topologies, with SF in and SF both topologies showing greater relative changes than ER and SF out do. This latter property gives an indication of the natural robustness to ICs of each topology, a property we explore in later sections.

The SF topology exhibits more oscillations and greater synchronicity than the ER topology

To investigate whether the SF topology more naturally generates oscillations than the ER topology does, we examined the dynamics of networks without evolution towards a target. For each topology 10 000 networks were generated and two measurements made. For the first, the period of every node in a network was found and the largest period recorded. A histogram of these periods is presented in Fig. 2. For the second measurement, the period of a randomly chosen individual nodes was recorded. A histogram of periods from these 10 000 networks is presented in Fig. 3.

In Fig. 2, we observe a clear divide between the SF out topology and the other two SF topologies. The SF out topology produces a probability distribution similar to that of the ER topology. The SF both and SF in topologies have similar distributions to each other. They are significantly more likely to have maximum periods at larger values of L. Fig. 3 examines the probability of finding oscillations of period L at a randomly chosen node. This result differs to the previous one as SF both now shows more oscillations than SF in , whilst SF out is clearly superior to ER at generating oscillations. Again the SF in and SF both topologies are considerably more likely than the SF out topology to have an oscillating node. In both Fig. 2 and 3 even oscillations are more easy to generate than odd ones.

We next measured the distribution of oscillations in the network during evolutionary runs towards L = 10 targets. The period of each node was measured in the first generation where all individuals in the population have a fitness of unity (defined to be the first maximally fit generation), giving an indication of the extent to which the target node's period is spread throughout the network. This measures the synchronicity within evolved networks.

The mean fraction of nodes oscillating with each period are presented in Fig. 4 for each of the four network topologies, averaged over 50 independent evolutionary runs. This figure demonstrates that a greater fraction of ER networks are frozen (L = 1). All the SF topologies have an increased number of oscillatory nodes and particularly L = 10 nodes. Periods that are factors of 10 (2 and 5) also appear to be much more prevalent than other periods.

In general we observe that the scale-free topologies are more likely to exhibit oscillations at a randomly chosen node, have larger maximum periods and show more synchronicity during an evolutionary run. We note that this increased ability of SF networks to produce oscillatory outputs is not universally observed: studies of random Boolean networks [START_REF] Drossel | Critical Boolean networks with scale-free in-degree distribution[END_REF] showed that critical SF networks possess more frozen nodes than networks with a fixed number of inputs. The increased oscillatory nature of SF networks in our model may be due to the modified update rule for single incoming nodes. Through this rule, any oscillatory node connected to a node with a single incoming connection will always propagate the oscillation to that node. Given that P S F (1) = 0.745 compared The probability is calculated by taking the mean fraction of nodes with each period, taken over 50 independent evolutionary runs. The "oth" column refers to other periods, i.e. all periods not labelled. This figure measures the degree of synchronicity in the networks. Note the difference in scale on the y-axes of the right and left graphs.

with P ER (1) = 0.450, guaranteed propagators are much more likely in SF in networks. The adaptive advantage generated by the SF out topology must have different origins, however, as the corresponding P( 1) is the same as for ER networks, and hence the modified rule affects both cases similarly.

In the rest of this paper we focus mainly on comparisons of the ER with the SF out distribution because the latter most closely resembles the topology of GRNs found in nature [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF][START_REF] Aldana | Robustness and evolvability in genetic regulatory networks[END_REF].

Evolution towards multiple oscillatory and frozen targets

In nature GRNs may employ phenotypes where multiple nodes require oscillating periods. This would be the case, for example, if a set of genes played a part in an oscillatory process within an organism. The segmentation clock [START_REF] Pourquie | The Segmentation Clock: Converting Embryonic Time into Spatial Pattern[END_REF] and circadian rhythms [START_REF] Hirata | Oscillatory Expression of the bHLH Factor Hes1 Regulated by a Negative Feedback Loop[END_REF] rely on multiple genes working together in this way.

To study this scenario, evolution was performed towards a target phenotype defined by requiring the output of several randomly chosen nodes' to have specific periods. We define the target period set, L, to comprise the period of each target node. For example, for N tgts = 3 with two L = 5 targets and one L = 10 target, L = {5, 5, 10}.

To determine the fitness of a network, each target is assumed to make an equal contribution. A random set of nodes are assessed against each of the randomly chosen target outputs of length specified in L and each node's fractional fitness is added. As such, we define fitness of the networks with multiple targets as s. the number of L = 7 targets for the ER (red, × markers) and SF out (blue, + markers) topologies with corresponding fits (the lines are guides to the eye). Both topologies experience an increase in adaptation time with an increase in the number of targets. This increase is more rapid for the ER topology, typically requiring more than 20 000 generations for N tgts > 6. The slower rate of increase of τ seen in the SF out topology is likely to be due to the greater synchronicity of this topology.

where F i is the fitness of the i th node in the phenotype and N tgts is the number of target nodes defining the phenotype. Each F i is given by Eq. ( 6).

The adaptation time required to achieve multiple targets of the same length increases rapidly with the number of targets

As the ER topology evolves rather slowly towards L = 10 targets, a smaller target length was selected. This allows a greater range of target numbers to be tested within a computationally reasonable number of generations. If the target length is too small, however, only a few different targets will be possible. For example, there are only 6 possible L = 5 targets (see Appendix C). Given these considerations, multiple L = 7 targets were chosen as this length is easily attainable by the ER topology whilst there are 18 possible L = 7 targets.

To measure speed of adaptation, we used the median number of generations until a population is maximally fit, τ, as in Section 3.1. In Fig. 5, τ is plotted against N tgts . As expected, for both ER and SF out topologies there is an increase in the adaptation time as N tgts increases. The change in adaptation time is non-linear, increasing rapidly with the number of targets. The rate of increase is much more rapid in ER networks than SF out ones. For example, SF out networks are still capable of adapting to N tgts = 20 in just over half the time ER takes to adapt for N tgts = 8.

As observed earlier in Section 3.2, the SF out topology shows greater synchronicity for oscillatory signals than the ER topology does. This difference may explain why the SF out topology performs so much better for multiple targets of the same length.

The adaptation time required to acheive multiple targets of different lengths depends on the number of targets and their relative periods

To test whether the enhanced synchronicity of the SF out topology is a key reason it performs so much better for multiple targets of the same length, we also examined different sets of multiple targets with different lengths. All the evolutionary runs were with ICs, and the results are shown in Table 3. Clearly having targets with multiple lengths greatly slows down the median adaptation time compared to having multiple targets of the same length for both topologies. It is also interesting to compare the performance of the networks for the L = {5, 10} targets to the performance for the L = {5, 8} targets. Although the former targets are slightly longer than the latter, they are easier to reproduce. Presumably this is because L = 10 is a multiple of L = 5 whereas L = 8 and L = 5 are coprime. We also note that the difference between the SF out and ER topologies is smaller for the coprime targets, suggesting that the enhanced synchronicity of the SF out topology, observed in Fig. 4, plays an important role in its enhanced ability to adapt to certain types of targets (phenotypes). Finally, we also tested the target L = {5, 10, 20} but were not able to find solutions within the 40 000 generation cutoff we used in our simulations, showing that multiple longer targets are significantly more difficult to achieve than a larger set of shorter targets (as in Fig. 5).

The evolvabilities of SF out and ER networks for multiple L = 1 targets are much more similar than for oscillatory targets

Thus far only targets of an oscillatory nature have been considered. This is because single L = 1 targets would be trivial to evolve; nodes of L = 1 are by far the most common in all attractors. With the extension to multiple target nodes, however, evolving to a large number of specific L = 1 targets is a non-trivial task. Biologically, a set of "frozen" targets could correspond to a particular constant gene expression pattern over some length of time, for example in an organisms' development.

As was done earlier, networks of each topology were evolved over 10 000 generations in 50 independent runs. Out of the N = 500 nodes, 100 were randomly chosen and required to possess a specific L = 1 target. We will refer to this target set as L = {100 × 1}. The values of τ with random and constant ICs are displayed in Table 4. The SF out topology maintains its advantage over the ER topology but the difference in evolutionary is the first maximally fit generation) for L = {10}. The errorbars represent the standard error on the mean. Robustness is still increasing after 10 000 generations for both topologies. The SF out topology is more robust from the outset and remains so throughout. Both topologies are able to evolve greater robustness. Noting the y-axis range, the change in robustness is small but statistically significant for both topologies.

adaptation times are vastly reduced. The effect of random ICs does not slow evolution down that much, and there is a similar relative decrease for both topologies.

Robustness in single & multiple target experiments

In the previous section we have investigated the evolvability of ER and various SF topologies towards oscillatory and frozen targets. In this section we study the mutational and environmental robustness of these networks under similar evolutionary situations. We particularly focus on the change in robustness under neutral evolution, after the population has achieved maximal fitness.

Mutational robustness is higher in SF topologies and increases under neutral mutation

To examine the mutational robustness of large networks with a large range of weights, a sampling approach must be taken as the number of possible network genotypes in a 1mutation neighbourhood is astronomical. To sample these genotypes, each node in a network was tested with 6 sample mutations of each mutation type in Fig. 1. We employ constant ICs throughout, ensuring that IC robustness effects do not play a role. An individual's mutational robustness, R ind m , is defined as the proportion of sample mutations that do not result in a network's phenotype changing. The robustness testing was performed on every member of a population during an evolutionary run at the generation where maximum fitness is attained (T = 1), and then at intervals (T = 2, 5, 10, 20, 50, 100, 200, 500, 1 000, 2 000, 5 000, 10 000) afterwards. At each interval, R ind m is averaged over the whole population producing an average population mutational robustness, R m . By averaging over the population, the statistical error is reduced. We show in Appendix E that the diversity of genotypes is maintained during neutral evolution.

Mutational robustness was measured for L = {10} and L = {100 × 1} targets, although we only show results for the former target in Fig. 6. The main and perhaps unexpected result for both types of target is that throughout these evolutionary periods the SF out topology maintains a greater robustness to random mutation than the ER topology, even though it is simultaneously more evolvable. We also note that R m is fairly high for both topologies, suggesting that the vast majority of mutations don't affect the phenotype. The most likely reason for this is that a large number of nodes simply do not participate in the core network that drives the output of the target node. Although it is not shown, the value of R m is also fairly high for the L = {100 × 1} targets (R m > 0.8) and moreover the SF out topology outperforms the ER topology here as well. Again, since 20% of the nodes are fixed by the target output, this suggests that a large number of the nodes do not participate in the output that determines the phenotype of these systems.

The second effect this figure clearly shows is that the mutational robustness increases with evolutionary time. It was shown in [START_REF] Van Nimwegen | Neutral evolution of mutational robustness[END_REF] that under neutral mutation, the mutational robustness of a population can increase. One mechanism for this phenomenon is that those genomes that are more mutationally robust are more likely to have offspring that share the same phenotype. Under continued selection then, such genomes are more likely to survive than genomes with low mutational robustness. These authors also showed that for selection to optimize mutational robustness in neutrally evolving systems, the condition μN pop 1 must be satisfied. Here the mutation rate per node is set to μ node = 0.02. For an N = 500 node network, this gives an average of μ = μ node N = 10 mutations per mutated offspring. With a fraction m = 0.75 mutated offspring in the extended population before the cut off is imposed, in each generation there are μ × m × N pop = μ N pop = 10 × 0.75 × 50 = 375 mutations in a future generation's genepool. Given that the condition μN pop 1 is satisfied, along with the genetic diversity observed in the populations (see Appendix E), further optimisation of mutational robustness is certainly expected. However, selection may not be the only cause of the observed effects. It may be possible for the increase to be obtained purely through entropic effects -for example, if many more genotypes are available at higher mutational robustness, the population will discover such possibilities simply through diffusing over the neutral space. Working out the exact reasons of this optimisation is thus a subtle question, and we don't think that our results unambiguously fix the cause of the increase in R m .

Environmental robustness increases markedly with time after maximal fitness is reached for random initial conditions, and is higher in SF topologies

GRN oscillators in vivo are often required to function under a variety of environmental conditions. By possessing environmental robustness, function can be maintained under such stochasticity and is thus a desirable trait. Here we examine the change in environmental robustness over time during neutral evolution and consider runs with both random and constant ICs.

Given the size of the networks (there are 2 500 ≈ 3 × 10 150 different initial configurations) a sampling approach was necessary. To estimate the environmental robustness, a random set of initial node states was taken, and the output phenotype was measured. This process was repeated for 1 000 different initial conditions (chosen randomly for each dynamics simulation) and the average fraction of initial conditions where the target phenotype was produced for an individual is called R ind d . The average of R ind d was taken over all the networks in the population to give R d , a measure of the environmental robustness (the subscript d refers to the dynamic nature of this measure of environmental robustness). This quantity was then further averaged over 50 evolutionary runs to reduce statistical error, and measured at the same intervals as in Section 5.1 after maximum fitness is attained.

Results for the mean environmental robustness Rd for L = {10} targets are shown for random ICs in Fig. 7 and for runs with constant ICs in Fig. 8. Firstly we observe that the SF out topologies are more environmentally robust than the ER topologies are, an effect that mirrors the mutational robustness. Secondly, the runs with random ICs show a marked increase of R d during neutral evolution. The runs with fixed ICs show a modest increase of R d during neutral evolution, especially for the ER networks.

When random ICs are used, it may be possible for selection to increase the environmental robustness of the population, for similar reasons to the arguments for the increase in mutational robustness. During the evolutionary runs those genotypes that The SF out topology naturally possesses greater environmental robustness. Over time the robustness of SF out changes insignificantly. The ER topology, on the other hand, almost doubles its robustness over the period measured. However, the ER topology still remains much less robust than the SF out topology.

are more environmentally robust, with larger attractor basins associated with the desired target, are more likely to have progeny that generate the required phenotype. Thus they are more likely to survive than less environmentally robust genotypes. Furthermore, we can make an estimate of the minimal environmental robustness (R min d ) for individuals to simply be sustained in the population. Let us assume individuals possess maximal mutational robustness and mutations affect their environmental robustness negligibly. Under these conditions, at least one of the four genotypes must produce the target phenotype again. This gives a value of R min d = 0.25. As the solution must propagate through the population, and is not fully mutationally robust, this is the lower bound and the actual value of R d required may be significantly larger.

For constant ICs there is no minimum environmental robustness which explains in part why the values of Rd are lower than for random ICs. Moreover, there is no direct selection pressure for environmental robustness that acts on runs with constant ICs. We observe that Rd does increase for ER networks. The reasons for this are not clear. A possible cause could be a correlation between mutational robustness, which is affected by selection, and environmental robustness. Such a correlation was also found by [START_REF] Mihaljev | Evolution of a population of random Boolean networks[END_REF], in a study of random Boolean networks. Another possibility would be an entropic effect where diffusion on a neutral space is more likely to find phenotypes with higher robustness.

Although they are not shown, the L = {100 × 1} targets show similar trends to the L = 10 targets: the SF out topology is more environmentally robust than the ER topology and furthermore both topologies show a modest increase in R d under neutral mutation for both random and fixed ICs.

Discussion

In this paper we have studied the interplay between evolvability and robustness of Boolean threshold dynamics models with ER and SF topologies. Our main conclusions are:

1) The enhanced evolvability of SF networks over ER depends on the phenotype under selection. We confirm the results of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] that networks with SF in topology are more evolvable towards single oscillatory targets than ER networks are, and find the same enhanced evolvability for the SF both topology and also for the SF out topology that most closely resembles the global topology of biological GRNs. The evolvability advantage of the SF out topology over the ER topology increases significantly for multiple targets with the same length. By contrast, for multiple targets of a different length, the difference in evolvability between SF out and ER topologies is reduced, and if the targets are coprime, the differences are even smaller. We also studied targets where 100 of the 500 nodes were given a fixed output. Here again the differences in evolvabilities between the topologies was greatly reduced compared to the L = 10 oscillatory outputs studied in [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. These results show that the evolvability advantage of a particular network topology depends strongly on the phenotype under selection.

2) SF networks show more oscillatory phenotypes and greater synchronicity than ER networks do. We compared the likelihood of finding an oscillatory output of period L 1 at a randomly chosen output node in networks without any evolutionary selection. The SF networks show a greater probability of finding oscillatory outputs than the ER networks do. Similarly, once a population has reached maximum fitness under evolution towards an L = 10 target for a single output node, a different node in the SF networks is much more likely than an ER network to show L = 10 oscillatory output: SF networks show much greater synchronicity than ER networks do. We argue that the increased evolvability towards oscillatory targets observed for SF networks is due to the fact that such phenotypes are more naturally generated by these topologies.

3) SF out networks are more mutationally robust than ER networks are. For both network topologies, the mutational robustness increases under neutral mutations. One might naively expect that the enhanced evolvability would lead to a reduced mutational robustness. However, as has been pointed out by numerous authors (Bloom et al., 2007a,b;[START_REF] Aldana | Robustness and evolvability in genetic regulatory networks[END_REF]Wagner, 2008b,a;[START_REF] Daniels | Sloppiness, robustness, and evolvability in systems biology[END_REF][START_REF] Draghi | Mutational robustness can facilitate adaptation[END_REF], this simple expectation does not necessarily hold for populations. Even without selection, SF topologies show a greater propensity than ER networks to exhibit oscillatory phenotypes. This difference suggests that the number of different genotypes that generate the same oscillatory phenotype is larger in SF networks, i.e. that the neutral spaces are larger. Larger neutral spaces may explain how greater mutational robustness can correlate with greater evolvability (Wagner, 2008b). The size of dynamic network neutral spaces has been investigated in studies of the evolution of random Boolean networks [START_REF] Mihaljev | Evolution of a population of random Boolean networks[END_REF][START_REF] Szejka | Evolution of Boolean networks under selection for a robust response to external inputs yields an extensive neutral space[END_REF], which have shown that these systems possess large neutral spaces which span search space. However, to firmly establish that the evolvability advantage of the SF networks can be explained by properties of the neutral space requires further investigation. For example, one must show that larger neutral spaces have evolutionary access to a greater phenotypic diversity.

Another possible factor is the heterogeneity introduced by the SF out topology. The measure of mutational robustness treats each node in the network equally (a biologically reasonable assumption), whereas, in an SF out topology, nodes are not equal: most nodes trivially have a single outgoing connection, whilst a few have a large number, the so-called 'hubs' of the network. The hubs in the SF out topology may play an important role, and mutations to these may be most effective, but rare as hubs only comprise a small part of the network. However, when mutations do occur to hubs (or to nodes closely connected to a hub), there is the potential for large scale influence of the mutation.

We also find that the mutational robustness increases under neutral mutations, an effect we attribute to the fact that more mutationally robust genotypes are more likely to have progeny that survive under future mutations [START_REF] Van Nimwegen | Neutral evolution of mutational robustness[END_REF].

4) SF out networks are more environmentally robust to changes in initial conditions than ER networks are, and both networks types exhibit an increase in environmental robustness under neutral mutations. The increased environmental robustness observed in SF topologies may be due to the comparatively large number of nodes acting in cohort to produce the required oscillatory signal -a feature shown by the observed high synchronicity in SF topologies. This large group of synchronised nodes may act to increase the number of initial states that converge to the required attractor, as perturbations to node states may be "damped out" by the overall tendency of the group to produce a particular oscillation. By contrast, the smaller groups of synchronised nodes in ER networks may be more sensitive to the initial states of nodes.

For random ICs, the environmental robustness increases under neutral mutations for reasons similar to those that explain how mutational robustness increases: more environmentally robust genotypes are more likely to survive in future generations when the environment changes, due to the increased attractor basin size associated with the desired target sequence. We find a similar, but more modest increase in environmental robustness in ER networks for fixed ICs, where there should be no selection for more robust phenotypes. At present it is not clear what causes this increase, but it may be correlated to the increase in mutational robustness.

An important question is whether the conclusions listed above hold only for the particular models we studied, or whether they are valid more generally. We show in Appendix A that similar conclusions also hold for a different threshold model [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF] that is frequently applied to study GRNs. Another consideration is the effect of varying network connectivity. For example, higher network connectivity can lead to the chaotic phase, and the network properties we discuss may be significantly altered in that regime [START_REF] Szejka | Evolution of Boolean networks under selection for a robust response to external inputs yields an extensive neutral space[END_REF].

The evolutionary dynamics of the threshold models we study here is considerably simpler than that which occurs in nature and so further work is needed to probe to what extent these conclusions hold for biological GRNs. At present such studies are very challenging because they require knowledge of the evolutionary history of GRNs, much of which is not easily accessible. The conclusions from this work that will most likely carry over are: (a) that the relationship between robustness and evolvability can depend on the class of phenotypes under selection; and (b) an organism's GRN can become both more evolvable and more robust simply though changing its overall topology. that has been investigated extensively with respect to GRNs. The essential differences between the Wagner model and the one of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF], is the use of different node states and update rules. In the original Wagner model, the gene states for node i are either σ i = 1, 0, -1. The update rules differ by using purely threshold rules, with no modified rule for single incoming degrees. The node state σ i = 0, occurs when the input sum S i = 0. Node states of 1 and -1 correspond to our "on" and "off" states considered in the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model. We investigated dynamics with the essence of the Wagner model. We note that nodes with σ i = 0 make no contribution to the dynamics in the next time step and given a continuous range of weights, such nodes will only occur when lacking any incoming connections. As a result, nodes with σ i = 0 will always remain in this null state throughout a dynamical run, rendering them ineffective within a network. Effectively, the update rules only differ from the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model because we can abandon the modified rule for single incoming connections. The effect of this rule is now incorporated into the behaviour of all nodes, as when a node is in its "off" state, a signal still affects the nodes it connects to. This is due to the "off" state now corresponding to σ i = -1, resulting in a non-zero amount being added in the threshold sum S i , even when the node is "off".

To examine the networks using the Wagner model, we firstly consider for what parameters the system enters the chaotic regime. In the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model, networks were within the ordered regime and need the same to be true here for comparison and biological realism. It was found that for the adapted ER topology (or SF in topology), a phase transition occurs for an average connectivity of k = 1. This gives a critical parameter of K c → ∞ for the adapted Poisson distribution. As a consequence, networks with this topology (or SF in ) cannot be produced in the ordered regime. This is an effect due to local noise propagation in the Wagner model being more likely than in the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model (see Appendix B). As a result we considered ER graphs with a standard Poisson distribution

P ERstd (k) = K k e -K k! (A.1)
rather than an adapted one. For the standard Poisson distribution, the phase transition occurs at K c = 2.075. This allowed us to compare non-adapted ER networks and SF out (SF out also possesses a non-adapted ER incoming degree distribution here) networks at the previous average degree connectivity, whilst re-

ER SF out

Figure A.9: The probability a random node has a given period, L, in the first generation where maximum fitness is attained, within the Wagner model. The probability is derived from the mean fraction of measured nodes taken over 50 independent evolutionary runs. The "oth" column refers to other periods, i.e. all periods not labelled. As in the model of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] the SF out topology exhibits greater synchronisity. maining in the ordered regime.

Appendix A.1. The evolvability of SF networks using the Wagner model is again greater than that of ER networks As in Section 3.1, we evolved networks towards random L = 10 targets over 50 evolutionary runs. The median generation numbers for a population to reach maximum fitness, τ, are shown in Table A.5 for each topology under random and constant ICs. Evolution with constant ICs is a much simpler task for both, with τ reduced by more than an order of magnitude for both topologies. Interestingly, the SF out topological advantage has vanished in this case. For random ICs, however, the results are dramatically different. The SF out topology now outperforms the ER topology by roughly the same factor as in the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model. Both topologies struggle to a much greater extent to adapt under random ICs. This is likely to be a consequence of local noise propagation being greatly increased in this model.

Appendix A.2. Networks with the Wagner model show more oscillatory behaviour

As in Section 3.2, the synchronicity in networks under evolution towards a single L = 10 target was examined but with the Wagner model. The period of each node in the first maximally fit population was measured, with an average taken over 50 evolutionary runs. A comparison of the probability of measuring each period is presented in Fig. A.9. As for the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model, the probability of a node having the same period as the evolved output (L = 10) is significantly higher in the SF out topology than in the ER topology. In the Wagner model, oscillations spread much more widely with both topologies having a much smaller fraction of L = 1 nodes. This is due to an increase in local noise propagation under the Wagner update rules. Appendix A.3. Mutational and environmental robustness of networks with the Wagner model also increase with time after maximal fitness is reached In Section 5, we considered the mutational and environmental robustness of networks evolving towards L = {10} targets. These experiments made use of the update rules used by [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. In this section, we confirm the overall results presented previously are not a consequence of those update rules, and are a general feature of threshold networks. The same method of measurements were used as in Section 5. out and that both topologies increase their mutational robustness over time. For the Wagner model the mutational robustness changes to a greater extent for both topologies than is observed for the model of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. This difference is likely due to the greater extent of local noise propagation in the Wagner model. Fig. A.11 and Fig. A.12 show the change in average environmental robustness over the maximally fit evolutionary period, for both random and constant ICs. As before, in the case of random ICs, environmental robustness can be optimised directly by selection. A significant increase in environmental robustness is observed in this case, although for both topologies, the typical maximum environmental robustness attainable starts levelling off around Rd ≈ 0.6. This is different to the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model, where a higher environmental robustness was generated by both topologies. Again this difference may be a consequence of local noise propagation being stronger in the Wagner model. The inability to evolve more complete environmental robustness requires further research. The analysis of the attractor properties within the Wagner model would probably be helpful. The result for constant ICs in Fig. A.12 shows that under neutral mutation there is a clear optimisation for both topologies. The networks naturally possess very low environmental robustness at T = 1, particularly for the ER topology.

Fundamentally the results are very similar to before with the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model. SF out topologies are consistently more evolvable than the ER topology, whilst both topologies increase their mutational robustness and environmental robustness under neutral mutation As can be seen there is limited difference between this and the original, demonstrating any bias does not affect the evolution of mutational robustness in this case.

Appendix E. Population Diversity

To determine whether the diversity of a population genotypes undergoing neutral mutations, we measured the genetic difference between every pair of inviduals each generation after maximal fitness is reached. The fractional genotypic distance D(A, B), between two networks A and B is given by

D(A, B) = 1 2M N i=1 N j=1 f (w A i j -w B i j ) (E.1)
where

f (x) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if x 0 0 if x = 0
and M is equal to the average number of entries in a network's weight matrix, given by M = N k . Values for D(A, B) are taken over all pairs A, B ∈ {population}, and then the values are binned in distances of 50/2M up to the maximal value of 1. These measurements were performed on several different types of evolutionary run. A heat map of these distances is displayed in Fig. E.15 for a typical evolutionary run (once maximum fitness is reached) for evolution towards a L = {10} target. This figure demonstrates several features of an evolutionary run. The gradually increasing "hot" curves show the population spreading out from each other in genotype space as independent mutations are incurred. The gradient of the "hot" line decreases as mutations start occuring on top of each other. However, the "hot" lines can still increase to large values (tending to D ≈ 1 for some lines), indicating that most of the mutations are neutral ones and lead to a large population diversity.

An interesting feature of the plot are the occasional dramatic decreases in diversity. These are indicated through the collapse of a "hot" line to a much lower value on a new curve, corresponding to the introduction of a new most recent common ancestor for the population. This ancestor may have been selected for or simply fixed through random sampling.

Figure 3 :

 3 Figure3: The probability of different periods being measured at a random node, in randomly generated networks. Note the logscale. SF both generates the most oscillations, followed by SF in , SF out and ER.

Figure 4 :

 4 Figure4: The probability of a randomly chosen node possessing period L in the first generation where maximum fitness for an L = 10 target is attained. The probability is calculated by taking the mean fraction of nodes with each period, taken over 50 independent evolutionary runs. The "oth" column refers to other periods, i.e. all periods not labelled. This figure measures the degree of synchronicity in the networks. Note the difference in scale on the y-axes of the right and left graphs.

  Figure5: τ v.s. the number of L = 7 targets for the ER (red, × markers) and SF out (blue, + markers) topologies with corresponding fits (the lines are guides to the eye). Both topologies experience an increase in adaptation time with an increase in the number of targets. This increase is more rapid for the ER topology, typically requiring more than 20 000 generations for N tgts > 6. The slower rate of increase of τ seen in the SF out topology is likely to be due to the greater synchronicity of this topology.

Figure 7 :

 7 Figure 7: Mean environmental robustness Rd , against T with random ICs for L = {10}. SF out is initially much more robust, although the gap is reduced over the course of neutral evolution. The increase in Rd is dramatic for both topologies, indicating a strong optimisation.

Figure 8 :

 8 Figure8: Mean environmental robustness Rd versus T , with constant ICs for L = {10}. The SF out topology naturally possesses greater environmental robustness. Over time the robustness of SF out changes insignificantly. The ER topology, on the other hand, almost doubles its robustness over the period measured. However, the ER topology still remains much less robust than the SF out topology.

  Figure A.12: Mean environmental robustness Rd , against T for evolution of L = 10 targets in the Wagner threshold model for constant ICs. Both topologies demonstrate a significant increase from very low initial levels. The increase generated may be consequence of the selection for mutational robustness also affecting environmental robustness.

Figure E. 15 :

 15 Figure E.15: The heat of a square represents the fraction of pairs of networks in the population separated by a genetic distance of D at time T. As can be seen, the genetic distance gradually increases as neutral mutations are incurred to each network. The increasing "hot" lines demonstrate the population diffusing through genotype space creating a genetically diverse population. The sudden drops in diversity, occassionally observed, correspond to a new most recent common ancestor taking over the population.

Table 1 :

 1 The topologies investigated are defined with respect to combinations of incoming and outgoing node degree distributions. When average connectivity is small compared with the number of nodes, the binomial distribution may be approximated as a Poisson distribution.

	Topology Incoming degree	Outgoing degree
	ER	Binomial/Poisson	Binomial/Poisson
	SF out	Binomial/Poisson	Power-law
	SF in	Power-law	Binomial/Poisson
	SF both

Table 2 :

 2 The median number of generations to reach maximal fitness, τ, for runs evolving towards a single L = 10 target. Results are shown for both random and constant ICs. All SF topologies outperform the ER topology. The SF in topology is most significantly affected by the difference between constant and random ICs, whilst the SF out is least affected. The probability of different periods being the largest in randomly generated networks. SF both and SF in perform similarly, whilst SF out and ER are similar too.

	ICs		ER			SF out			SF in	SF
	Constant 950			450			51		25
	Random 3 200		680			390	200
	Relative		3.4			1.5			7.7	8.0
		0.6
												ER
		0.5										SF out
	P (L)	0.3 0.4										SF in SF both
		0.2
		0.1
		0	1	2	3	4	5	6	L	7	8	9 10 11 12
	Figure 2:

Table 3 :

 3 Values of τ for evolutionary runs evolving under random ICs to different sets of periods. τ is smaller for L = {5, 10} compared with L = {5, 8} for both topologies. For the L = {5, 10, 20} no convergence was achieved.

	Period set, L	ER	SF out
	{5, 10}	7 000	2 400
	{5, 8}	9 900	4 700
	{5, 10, 20}	> 40 000	> 40 000

Table 4 :

 4 Median adaptation time τ for runs evolving towards L = {100 × 1} targets. The performance of the topologies is much more similar than for evolution towards oscillatory targets, although the SF out topology retains a slight advanatage.

	ICs	ER	SF out
	Constant	255	193
	Random	402	278
	Relative	1.58	1.44
		1	ER
		0.99	SF
		0.98
		0.97
	Rm	0.96
		0.95
		0.94
		0.93
		0.92
		0.91
		1	10	100	1,000	10,000
				T

out Figure 6: Mean mutational robustness Rm , against generation T (where T = 1

Table A .

 A 5: The results for evolution towards L = 10 targets with the A.Wagner model. τ is shown for both constant and random ICs. ER and SF out topologies perform similarly under constant ICs. However, intrinsic robustness to ICs may play a more important role.

	ICs	ER	SF out
	Constant 25	25
	Random 19 000	3 800
	Relative	850	150

  Mean mutational robustness Rm , against T for evolution of L = 10 targets in the Wagner threshold model. As found for the[START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model, neutral evolution can optimise muational robustness over time. The relative change in mutational robustness attained in this case is larger than that in the[START_REF] Oikonomou | Effects of topology on network evolution[END_REF] model, whilst the typical maximal mutational robustness values attained after 10 000 generations of neutral evolution are similar. Mean environmental robustness Rd , against T for evolution of L = 10 targets in the Wagner threshold model for random ICs. Selection increases environmental robustness in both topologies. The limit of this increase appears to be reached towards the end of the neutral evolutionary period, levelling off at around Rd ≈ 0.6. The overall lower level of environmental robustness is a consequence of increased local noise propagation in the model.

		1
		0.98	ER SF out
		0.96
		0.94
	Rm	0.92
		0.9
		0.88
		0.86
		0.84
		1		10	100	1,000	10,000
					T
	0.2 0.3 0.4 0.5 0.7 Figure A.10: 1 Rd 0.6	ER SF out	10	100	1,000	10,000
					T
	Figure A.11:

Appendix A. Comparison with the Wagner model

To examine the generality of the dynamical properties introduced by the SF out topology, we made a comparison with another threshold model that was proposed by [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF] and

Appendix B. Calculation of order-chaos phase transition

In order to calculate the order-chaos transition, we followed the method of [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]. Let us consider two different sets of ICs, {σ (1) i } and {σ (2) i }. The difference between these sets over time can be measured as with the Ham- ming distance as

The annealed approximation of [START_REF] Derrida | Random networks of automata: A simple annealed approximation[END_REF] measures whether noise propagation diminishes or extremifies in the limit of the Hamming distance between the two configurations tending to 0

where M measures the noise propagation in the network. When M > 1, the network is in the chaotic regime, whilst M < 1 is the ordered regime. For M = 1, critical behaviour is obtained. P(k i ) is the usual degree distribution, whilst p s (k i ) is the probability of local noise propagation. This measures the probability that if one of the inputs to a node is flipped at time t, this will change the output of the node at time t + 1. For threshold models this is a function of k i , unlike in the Boolean function models. Here we calculate numerically the values of p s (k i ) for both the update rules in the [START_REF] Oikonomou | Effects of topology on network evolution[END_REF]

Appendix C. Total number of sequences for each repeat length

There are a finite number of possible configurations a sequence of bits of a given length may take. Two sequences are regarded as identical if they can be mapped on to each other through cyclic permutation. If a sequence can be broken down into a multiple number of shorter sequences, it is invalid. As such, the number of possible sequences A(n), for a sequence of a length n, is given by the following

where F is the set of all factors of n and A(1) = 2 (i.e. 0 and 1).

Appendix D. Importance of mutation method with respect to increased robustness One possible explanation for the increase in mutational robustness of the SF out topology would be a possible bias introduced by the weight mutation rules. This bias could arise in the following way. Given that any edge, incoming or outgoing, can have its weight mutated for the node under mutation, the probability of selecting an incoming node for a hub in the SF out topology will be much smaller. Explicitly, there are many more outgoing edges that could be randomly selected for mutation. To assess whether this mutation method was introducing erroneous results, we performed runs towards L = {10} with modified mutation rules. Instead of allowing the weight of any edge touching a given node to be mutated, we only allowed weight mutations of the edges incoming to that node. As both topologies have the same incoming degree distribution, this removes any possible bias.