Compensatory mutation can drive gene regulatory network evolution

Gene regulatory networks underlie every aspect of life; better understanding their assembly would better our understanding of evolution more generally. For example, evolutionary theory typically assumed that low-fitness intermediary pathways are not a significant factor in evolution, yet there is substantial empirical evidence of compensatory mutation. Here we revise theoretical assumptions to explore the possibility that compensatory mutation may drive rapid evolutionary recovery. Using a well-established in silico model of gene regulatory networks, we show that assuming only that deleterious mutations are not fatal, compensatory mutation is surprisingly frequent. Further, we find that it entails biases that drive the evolution of regulatory pathways. In our simulations, we find compensatory mutation to be common during periods of relaxed selection, with 8-15% of degraded networks having regulatory function restored by a single randomly-generated additional mutation. Though this process reduces average robustness, proportionally higher robustness is found in networks where compensatory mutations occur close to the deleterious mutation site, or where the compensatory mutation results in a large regulatory effect size. This location- and size-specific robustness systematically biases which networks are purged by selection for network stability, producing emergent changes to the population of regulatory networks. We show that over time, large-effect and co-located mutations accumulate, assuming only that episodes of relaxed selection occur, even very rarely. This accumulation results in an increase in regulatory complexity. Our findings help explain a process by which large-effect mutations structure complex regulatory networks, and may account for the speed and pervasiveness of observed occurrence of compensatory mutation, for example in the context of antibiotic resistance, which we discuss. If sustained by in vitro experiments, these results promise a significant breakthrough in the understanding of evolutionary and regulatory processes.


Introduction
Compensatory mutations are common and relatively scale-invariant 124 We first test whether compensatory mutation is common in the context of the synthetic GRNs. We 125 find that, unlike deleterious mutation, the frequency of compensatory mutation is almost scale-invariant. 126 From Fig. S1A and B, we can see that the stability and robustness in initial networks are quite different 127 among varying sizes and levels of connectivity of gene regulatory networks. Which type of network, once Figure 2: The frequency of compensatory mutation is relatively insensitive to network size and network connectivity in the context of gene regulatory networks. For each network size (N = 5, 10, 15, 20, 30 and 40 genes) for each value of network connectivity (proportion of regulatory relationships) given from a range of values in continuous intervals ([0.2, 1], step size 0.02), the frequency of experiencing first deleterious followed by a compensatory mutation on just two rounds of mutation was tested based on an initial 10, 000 randomly generated stable gene networks. The shaded areas represent 95% confidence intervals based on 100 independent runs.
increases, but very gradually. In contrast, for the larger networks N = 30 and 40, with the rise in connectivity, 132 the compensatory mutation rates decrease slightly. However, overall the results indicate that the frequency 133 of individuals that can be fixed by compensatory mutation is more sensitive to network size than to network 134 connectivity, and not particularly sensitive to either. The implied probability of compensatory mutation 135 from the relative frequencies observed ranges from 8% to 15% of compromised networks recovering, with the 136 larger rates associated with larger networks. This marked scale-invariance (see Fig. S1C, which is identical  Next, we investigate the occurrence of compensatory mutations in populations that have been exposed 140 to bouts of generations of relaxed selection and selection for network stability. Again as a slight modification 141 of what is standard for this type of model, we assume selection favours network stability. We find that seen in a larger network with size N = 20 genes (see Fig. S4B), where the frequency of mutations being compensatory, if they occur on the original deleterious site, is 85%. The percentages beside each edge in these figures indicate the proportion of mutations that occur on that edge that are compensatory, out of the 181 1, 000 simulated second rounds of mutation we ran on each edge for each network after it had previously 182 suffered a single deleterious mutation. In general, as these representative figures indicate, the compensatory 183 effect could happen in many positions in a broken network, but it is more likely to be observed at sites that 184 are close to a deleterious mutation's site.  Fig. S4. It illustrates the 186 frequency among 10, 000 initially-stable gene networks of compensatory mutation against different spatial 187 distances from the single deleterious mutation suffered by each network. As can be seen, compensatory 188 mutations generally occur in edges between genes close to the deleterious mutation site. We restrict the 189 analysis to these five categories because there is only a narrow range of distribution distances for randomly 190 sampled mutations (see Appendix for more details). We further conducted similar experiments for networks 191 with neutral mutations to investigate whether compensatory mutations have any special property in terms of 192 location. We found that, compared with the results of compensatory mutations, neutral mutations are more 193 evenly distributed. Specifically, instead of measuring the frequency of a second, compensatory mutation (that 194 restores network stability for a compromised network with a single deleterious mutation), we measure the 195 frequency of a second, neutral mutation with different distance effects that retains the stability for a network 196 that has already had one neutral mutation. From Fig. 3A (dashed line) we can see that the distance effect 197 has a much less profound role in networks with two consecutive neutral mutations than in networks with 198 one deleterious mutation and one compensatory mutation. In fact, neutral mutations tend to be enriched if 199 they are far apart in larger networks (see Appendix). 200 The point of compensation is of course to recover the network's fitness, or here, stability. We therefore 201 next investigate whether there is an impact on the robustness of networks after deleterious followed by com-202 pensatory mutation varies by the location of the compensatory mutation, and contrast this with differences 203 in robustness after two neutral mutations. We again find that patterns are very dependent on location. 204 Specifically, we compare robustness of stable networks following one round of deleterious and compensatory 205 mutation with that of stable networks with two consecutive neutral mutations, as shown in Fig. 3C. In gen-206 eral, robustness is far higher when compensatory mutation occurs closer to the original deleterious mutation 207 site (see the solid line in Fig. 3C), whereas after two neutral mutations, closer distances are not better asso-208 ciated with higher robustness (see the dashed line in Fig. 3C) 1 . Even though networks with compensatory 209 mutations occurring near to the site of the deleterious mutation exhibit profoundly more robustness than 210 1 Note that robustness is higher overall for networks having experienced two neutral mutations-this is unlikely to be caused by the mutations, and more likely to be a characteristic of the network likely to contain neutral mutations. Figure 3: Both frequency and robustness of networks with compensatory mutations exhibit biases in location and mutation size. For N = 5 and c = 0.4, we first collected a pool of compromised networks with deleterious mutations after a single mutation round. We then forced second mutations, classifying these as being 0 (on the same site), 1, 2, 3 and 4 steps away from the original deleterious mutations (A) or adding a weight from [−5, +5] (step size 0.5) to the original regulatory impact (B). For each of these distance or weight categories, we measured the probability that the mutation was compensatory (that it returned the network to stability, see solid lines in (A) and (B)), based on 10, 000 sample networks collected for each category. The sample networks for control groups (see dashed lines in (A) and (B)) were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Similarly, for N = 5 and c = 0.4, we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation that was 0, 1, 2 and 3 steps away from the previous deleterious mutation (C) or with different shifts in gene regulation from [−5, +5] (step size 1 and with four additional regulation shifts: −0.5, −0.1, 0.1 and 0.5) (D). Then, we assessed the robustness of the sample networks at each category (see solid lines in (C) and (D)). The sample networks for control groups (see dashed lines in (C) and (D)) were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Error bars represent 95% confidence intervals based on 100 independent runs. those at other locations, their actual robustness is much lower than that of networks with neutral mutations 211 (see Fig. 3C and Appendix). Nevertheless, these theoretical results indicate that these co-localised compen-212 satory mutations are more likely to be accumulated, whereas compensatory mutations that are far apart 213 from the previous deleterious mutations are more likely to be lost during subsequent selection for network 214 stability.
Independent of location, we also investigated how different mutation size influences the probability of com-216 pensation in compromised networks. We found that compensation is more likely to be driven by large-effect 217 mutations. Fig. 3B (solid line) presents the frequency of compensatory mutation against various intensities 218 of up or down regulation among 10, 000 randomly-generated stable gene networks that had experienced a 219 single deleterious mutation. For a randomly-chosen site in each network, we experimented with mutations 220 across a range of regulatory strengths. As can be seen, larger regulation changes, both positive and negative, 221 are up to a point associated with an increased frequency of compensatory mutation. However, the shape 222 of the curve for compensatory mutations across all edges is not a symmetrical 'V'. Rather, compensatory 223 mutations occur more by positive changes to gene regulation than by negative changes. The explanation for 224 this phenomenon is rooted in the fact that there are two edge types that can be affected by compensatory 225 mutation: inter-gene regulation connecting two different genes and self-regulating edges. In the simulations, 226 almost no compensatory mutations are both negative and self-regulating (see Appendix). The 'V' shape 227 for only inter-gene regulation is almost symmetrical (see Appendix), suggesting that for these, negative and 228 positive regulations are equally likely to be useful. It is true for both the negative and positive cases that 229 compensatory mutation is increasingly likely with greater regulatory strength up to a certain extent.

230
Although we found that compensatory mutation tends to positive, this is not a property special only 231 to compensatory mutations. From Fig. S5A, we can see that there is more positive regulation in both 232 initially-stable networks and networks with compensatory mutations, whereas deleterious mutations in com-233 promised networks tend to be more negative. By separating self-and non-self-regulatory edges, we find that 234 compensatory mutations have a larger effect (in terms of shifting gene regulation) on self-regulatory edges 235 than non-self-regulatory edges (see Fig. S5B and C). We then conduct similar experiments for networks with 236 neutral mutations to investigate whether compensatory mutations have any special property in terms of 237 mutation size. We find that, compared with the results of compensatory mutations, small-size mutations are 238 more likely to be observed in networks with neutral mutations. Specifically, similar to the location experi-239 ments, we measured the frequency of a second mutation (neutral mutation) with different mutation effects 240 that can retain the stability for a network that has already had one neutral mutation. From Fig. 3B (dashed 241 line) we can see that that neutral mutations are more likely to be of small magnitude than large, and where 242 they are large they are more likely to be positive than negative. Compensatory mutations are more likely to 243 be of large magnitude than small, but are also more likely to be positive than negative.

244
As with location, we also investigated the robustness of networks subject to mutations of different sizes. 245 We found that patterns of shifting regulation-generating robustness are also quite different. Specifically, 246 we compared robustness of stable networks having one deleterious mutation and compensatory mutation 247 with that of stable networks having two consecutive neutral mutations, as shown in Fig. 3D. In general, the 248 robustness is significantly higher when compensatory mutation has a larger shift in gene regulation (see the 249 solid line in Fig. 3D). Although networks with neutral mutations tend to have a similar pattern (see the 250 dashed line in Fig. 3D), by measuring the percentage change in robustness (see Appendix), we can clearly 251 see that size has a greater impact on robustness in compensatory mutations. Here again, it should also 252 be noted that although networks with compensatory mutations exhibit a more profound biased change in 253 robustness with respect to mutation size, their actual robustness is lower than that of networks with neutral 254 mutations (see Fig. 3D). These theoretical results indicate that these large-effect compensatory mutations 255 are more likely to be accumulated, whereas small-effect compensatory mutations are more likely to be lost 256 by subsequent selection for network stability.

257
Note that similar patterns to those described here are also observed in networks with different sizes and 258 connectivity. See more supporting figures in the Appendix.

259
Compensatory mutation generates regulatory complexity 260 We now explore the long-term evolutionary consequences of compensatory mutations. Given the biases 261 identified in the previous section concerning location and magnitude, we might predict that the effects 262 of these two fundamental network properties would facilitate an altered neutral evolution, at least during 263 periods of relaxed selection. Recall that we assume such relaxed periods will be interspersed between bouts 264 of selection for network stability, and also that selection in our model favours network stability. We in fact do 265 observe in our simulations an increase in the complexity of gene regulatory networks, but only in a context 266 where they have been withdrawn from the selection for network stability for at least some proportion of 267 generations. Specifically, we first generate a pool of 10, 000 stable networks (N = 10) with a simple 'Star' 268 topology (see Fig. 4A), then evolve the population under different evolutionary scenarios. Fig. 4B shows four 269 evolutionary scenarios where the population is exposed to selection for network stability in every generation 270 such that there is no opportunity for compensatory mutation. From the typical results (networks with a 271 median connectivity), we find that: 272 1. the median connectivity is the same as the initial population's if it is evolved without mutation or : Compensatory mutation generates regulatory complexity in stable networks without an initial variation in network structure. The initial population pool was composed of 10, 000 sample stable networks with N = 10 genes. These networks had a similar "Star" topology (one hub node and nine non-hub nodes), varying network connectivity in [0.10, 0.26]. The detailed description of generating initial population can be found in Appendix. A representative network from the initial population is shown in (A). The initial population was evolved for 5000 generations with the selection for network stability (f RS= 0) under no mutation and no recombination regime, mutation but no recombination regime, recombination but no mutation regime, mutation and recombination regime (B). The initial population was also evolved for 5000 generations under relaxed selection regime (unstable networks will not be eliminated in generations under relaxed selection) with frequency f RS = 1/10, 1/25, 1/50 (C). Note that compensatory mutation cannot happen when the population is subject to selection for network stability, since no deleterious mutations survive to be compensated. The representative networks were selected randomly with a median connectivity in the evolved populations.  selection, invoking mutation (including compensatory mutation) and recombination. From these typical and increases and is higher than in the case when the population is subjected exclusively to the selection for 283 network stability so that no compensatory mutation can occur. 284 We would like to quantify the impact of relaxed selection on regulatory complexity. In another experiment, 285 we collect 10, 000 stable networks and then evolve them for 5, 000 generations, allowing both mutation and 286 recombination. From Fig. S6, we can see that if there is no relaxed selection at all, the mean connectivity 287 of the population is highly preserved during evolution, whereas the network connectivity increases if we 288 allow compensatory mutations to occur in periods of relaxed selection. It should be noted that in the first 289 experiment, as shown in Fig. 4, we fix the network structure but vary the network connectivity in the 290 initial population, whereas we fix the network connectivity but vary the network structure in this second 291 experiment. These results demonstrate that selection for network stability where it impedes deleterious and 292 compensatory mutations constricts complexity, whereas compensatory mutations contribute to regulatory 293 complexity as a part of neutral process. Compensatory mutations have long been considered the primary means by which low-fitness lineages 296 might be able to be restored to high fitness [Levin et al., 2000, Crawford et al., 2007, Meer et al., 2010.

297
More recently, Dunai et al. [2019] suggest compensatory mutation may account for robust adoption of costly 298 traits that are of critical importance to an organism in certain circumstances, such as antibiotic resistance.

299
Given that most mutations are believed to be deleterious at least initially, some such process would be 300 essential for mutation to contribute to genetic innovation and evolution more generally. However, the extent 301 of the role of compensatory mutations has often been considered to be negligible because they were considered 302 to be highly improbable and therefore rare. As such, they have not been studied extensively, and many of 303 their general properties have been unknown.

304
If the results presented here in simulation hold for in vivo regulatory networks, then compensation 305 may be far more probable and frequent than had previously been anticipated. Our results indicate that gene 306 networks may by their nature be surprisingly robust, such that a wide variety of alterations to a compromised 307 network may effect its recovery. Significantly, we find that the frequency of compensatory mutation -unlike 308 deleterious mutations -is relatively invariant to the size of the network. This may mean that iterations 309 of deleterious and compensatory mutation play a far larger role in evolution than previously thought. Our  Dunai et al. [2019] find that compensatory mutation rather than reversion are actually fairly outcomes for 382 bacteria developing resistance to multiple antibiotics at once. In our present paper, we have demonstrated 383 theoretical support and explanation for these empirical findings, by showing that compensatory mutations 384 can be greatly increased if the population evolves during a phase of relaxed selection regime (Fig. S2). Here,

385
'relaxed' is more obviously a relative term. Organisms challenged by antibiotics are already stressed, and 386 trade-off the costs of initial mutations with the benefits in surviving the antibiotic assault. In this climate, 387 'compensatory' mutations are mutations selected because they mitigate these additional costs, reducing the 388 chance that the mutations leading to antibiotic resistance are swept through reversion from the population.

389
Although antibiotic resistance is obviously a problem in human contexts, more generally this sort of dynamic 390 illustrates a context in which 'relaxed' selection might come into play -when a mutation produces both 391 costs and benefits, or these vary with the ecological context of the organism.

392
Many studies have shown that conventional de novo mutations are widely distributed throughout the 393 genome and have a wide distribution of phenotypic effects, from complete lethality to weak benefit with 394 respect to fitness [Sanjuán et al., 2004 Mezmouk and Ross-Ibarra, 2014]. Although there have been no predictive tests of the location of compen-396 satory mutations, empirical studies show that compensatory mutations are often found in proteins that are 397 in or interact with proteins that exhibit a deleterious mutation , 398 Davis et al., 2009, Comas et al., 2012, Bhattacherjee et al., 2015. Our findings concur with this. In this 399 paper, we have showed that there is a bias with respect to where compensatory mutations happen such 400 that compensatory mutations tend to generate regulatory circuits that closely interact with each other (solid  Previous work has indicated that the origin of mutational robustness may come from the non-adaptive 407 results of biophysical principles or non-adaptive evolutionary forces [Ruths andNakhleh, 2013, Payne and408 Wagner, 2015]. During periods of relaxed selection, regulatory networks with otherwise-lethal mutations 409 have the potential to be compensated by additional mutations. If compensatory mutation occurs frequently 410 enough and generates different patterns of gene regulation than networks with neutral mutations, then the 411 processes observed here could alter which types of network are lost when selection for network stability 412 does occur. Systematic biases in the loss of particular network configurations could allow network features 413 associated with compensatory mutation to accumulate in the population, even when the features do not 414 at least initially confer differential reproductive success. In addition-as we have shown-the combination 415 of recombination, deleterious mutation and compensatory mutation under moderately effective population 416 sizes could then permit the evolution of increased regulatory complexity. In this paper, we have shown 417 that stable networks with compensatory mutations generating a profound change in robustness compared 418 to the impact on stable networks of neutral mutations ( Fig. 3C and D). These results indicate that over 419 time, compensatory mutations that occur during generations of relaxed selection could be biased such that 420 regulatory circuits that closely interact and have larger interactive impacts are more likely to be maintained. 421 We have also shown that, at least in our system, over time these can have profound impact on the complexity 422 of the networks.

423
Taken together, we believe these findings demonstrate that the nature of compensatory mutation has 424 been misunderstood theoretically. Periods of relaxed selection (as per our model) or indeed of increased 425 differential selection (as per ecological challenges e.g. of antibiotics) produce a context in which natural 426 innovations may be tolerated long enough to be combined. Combining mutations allows for a larger range of 427 genomic innovation. Both our models and the empirical data of others show a surprising level of resilience in 428 complex biological systems. Our models indicate that resilience may scales well with increasing complexity.

429
Overall, the biases that emerge in this process as innovations accumulate may be an important new factor 430 in understanding the evolution of gene regulatory networks, and evolution more broadly.

432
We employed a well-established synthetic model of gene regulatory networks to simulate compensatory where f (·) is a sigmoidal function, and i is a constant which reflects either a basal transcription rate of gene In all the simulations here, we define network developmental stability as the progression from an arbitrary 450 initial expression state, S(0), to an equilibrium expression state (reaching a fixed phenotypic pattern), S EQ , 451 by iterating Equation (1)  for network stability is also referred to as selection for network stability in which unstable networks will be 454 eliminated. The equilibrium expression state can be reached when the following equation is met: where ξ is a small positive integer and set to be 10 −4 in all simulations, and D(S, of transcription or the time necessary to export mRNA into the cytoplasm for translation [Wagner, 1994].

462
Initialisation 463 Each individual network in the population was generated with a gene regulatory matrix W associated 464 with an expression state vector S(0). Specifically, the matrix was generated by randomly filling W with 465 c × N 2 non-zero elements w i,j that was drawn from a standard normal distribution N (0, 1). The associated 466 initial expression state S(0) was also set by randomly choosing each s i (0) = +1 or −1.

467
In the mutation operation, exactly one element w i,j picked at random in each regulatory matrix W would 469 be replaced by w i,j ∼ N (0, 1). Note that the mutation only occurs among non-zero elements. In other words, 470 the mutation process will not change the topology of the original network W in terms of forming new edges 471 or deleting existing edges between two genes.

472
Recombination 473 In some simulations presented in this paper, we allowed individual networks to recombine with each 474 other. A recombinant was produced by picking two individuals and selecting rows of the W matrices from 475 each parent with an equal probability. This process is similar to free recombination between units formed 476 by each gene and its cis-regulatory elements, but with no recombination within regulatory regions.

477
Strong and relaxed selection for network stability

478
In the selection for network stability regime, only individuals which were able to attain developmental 479 stability after the mutation process were selected. In contrast, all individuals can survive regardless of they 480 were capable or incapable of reaching equilibrium when the selection for network stability was relaxed.  Figure S1: The influence of the size and connectivity of a gene regulatory network on its initial stability, robustness and frequency of compensatory mutation. For each network size (N = 5, 10, 15, 20, 30 and 40) with each connectivity given from a range of values in continuous intervals ([0.2, 1], step size 0.02), we tested the proportion of gene networks that are stable based on an initial 10, 000 randomly generated networks (A), the robustness of stable networks after exposure to a single round of mutation based on an initial 10, 000 randomly generated stable networks (B), and the frequency of compensatory mutation based on an initial 10, 000 randomly generated stable networks (C) (rescaled from Fig. 2). The shaded areas represent 95% confidence intervals based on 100 independent runs. Figure S2: Compensatory mutations can occur in networks regardless of different patterns of selection. For each network size (N = 5, 10, 15, 20, 30 and 40) with network connectivity c = 0.76, we collected 10, 000 only stable, both stable and unstable, and only unstable networks with one to fifteen rounds of mutation. For each round of mutation, each network was subjected to one single mutation (for unstable networks) or two single mutations (for stable networks). Then, we measured the frequency of compensatory mutation in networks that have been subjected to bouts of selection for network stability (A), networks that have been subjected to bouts of relaxed selection (B), and networks with cumulative deleterious mutations (C). The error bars represent 95% confidence intervals based on 100 independent runs.  , for a particular compromised network that was stable initially, we executed one additional mutation round 1, 000 times on each edge. Then, the percentage of each mutation on that edge that restored GRN stability after mutation was measured. Unmarked edges had a CM 0% of the time. Note the solid line with width also indicates the probability an edge's mutation was compensatory, and the dashed line to represent the edges for which mutation never compensated this particular deleterious mutation. The original deleterious mutation occurred on the edge marked in red. Note: The directed edge represents the interaction between two connected genes. But we do not distinguish negative or positive regulation in the provided examples. Figure S5: The distribution of regulation in initially-stable, compromised and restored networks. For randomly generated stable networks with N = 5 and c = 0.4, we collected 10, 000 sample regulations. We also collected 10, 000 sample regulation weights from deleterious mutations that compromised initially-stable networks as well as from compensatory mutations that restored the stability of previously-broken networks. We then measured the distributions in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). Given that the regulations are continuous values, we grouped them into 19 bins from [−4.5, +4.5] (step size 0.5). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S6: Compensatory mutation generates regulatory complexity in stable networks without an initial variation in connectivity. For the network size N = 40 and the connectivity c = 0.15, we collected 10, 000 stable networks, then evolved them for 5000 generations, allowing recombination at each generation. In every 200 generations, we measured the network connectivity of the population (stable) in which the relaxed selection occurs in every 2, 10, 25 50 and 200 generations. We also measured the network connectivity of the population when there was no relaxed selection as a control group. Shaded areas represent 95% confidence intervals based on 10 independent runs.

Supplementary Text
Estimating the relative frequency of compensatory mutation To gain an impression of the properties of the initial gene regulatory networks, we first tested the probability of network stability in randomly-generated networks. As illustrated in Fig. S1A, smaller networks are more likely to be stable. Moreover, the relative frequency of stability in networks with low levels of connectivity is higher than that of networks with high levels of connectivity. This is in general accordance with previous work, typically done at connectivity c = 0.75, e.g. Azevedo et al. [2006], which indicates that larger networks with complex topology tend to be unstable.
In the second experiment, we explored the robustness of initially-stable networks; that is, we investigated the probability that stable networks remain stable after a single round of mutation. Here, a single mutation means exactly one non-zero entry in an individual's genotype would be mutated. Given that the initiallystable networks were collected from the original randomly-generated ones, it would seem reasonable to predict that the small stable networks are more likely to break after one mutation round, since they contain fewer pathways and a single mutation, therefore, has a greater proportional effect. However, the results in Fig. S1B show the opposite effect: the stability of the small networks is still high. The mutation operation is effectively an alternative way of generating new networks; thus, the mutated networks have the same properties as the initial ones.
In our third experiment, we measured the compensatory mutation frequency in previously-stable networks. Specifically, we started from a population pool where each stable network was randomly generated.
Then, we exposed these initially-stable networks to a single round of mutation. We focused on those unstable networks where each network contained a single deleterious mutation. Next, we exposed these compromised networks to an additional round of mutation. Finally, we tested the stability of the resulting networks. The stable networks at this point had experienced compensatory mutation. We then measured the frequency of individuals that experienced compensatory mutation. As shown in Fig. S1C (also see Fig. 2, main text), the frequency of compensatory mutation is largely scale invariant both to the network size and the network connectivity.

Exploring strong and relaxed selection for network stability on compensatory mutation frequency
In this set of experiments, we investigated the frequency of compensatory mutation after many generations of both strong and relaxed selection for network stability to test whether compensatory mutation continues to occur even after lengthy evolution (see Fig. S2A and Fig. S2B). Specifically, under the selection for network stability regime, we collected 10, 000 stable networks at each generation where each network in the population was subjected to one single mutation. Then, we performed another round of mutation, focusing on the unstable networks that resulted from the previous round, and measured the probability of a second mutation that can restore the network stability of those compromised networks. Similarly, under the relaxed selection regime, we collected 10, 000 networks at each generation where each network in the population was subjected to one single mutation. However, for each relaxed selection generation, there were both stable and unstable networks after the population was subjected to the single mutation, since we did not perform selection restricting networks to being stable. The overall frequency of compensatory mutation for the population during each relaxed selection generation was averaged over the results of stable networks and unstable networks that were calculated separately.
In addition, we also measured the frequency of compensatory mutation among unstable networks during each relaxed selection event to further confirm that compensatory mutation can occur even in seriously damaged networks (see Fig. S2C). Specifically, we collected 10, 000 unstable networks at each generation where each network in the population was subjected to one single mutation, so really in this case we had selected against network stability. Then, we performed another round of mutations and measured the probability of a second mutation that could restore network stability. Note that this set of experiments is similar to those experiments described above, but here we only focus on unstable networks, whereas we consider both stable and unstable networks in the relaxed selection regime.
Exploring the frequency of relaxed selection in simulating compensatory mutations In this set of experiments, we tested whether frequent relaxed selection can generate more compensatory mutations (see Fig. S3A). Specifically, we collected a population pool of 10, 000 stable networks that were generated randomly. The initial population was then evolved under a relaxed selection regime with a frequency of 1/2, 1/5, 1/10, 1/25, 1/100, 1/200 and 1/500 for a total of 1, 000 generations. Note that during relaxed selection event, both stable and unstable networks can survive when the population is subjected to one single round of mutation. The number of compensatory mutations was recorded immediately after each relaxed selection event when the population was subjected to another single round of mutation. The reported results are the total number of compensatory mutations (see Fig. S3A) and frequency of compensatory mutation (per network per relaxed selection event, see Fig. S3B) arising over 1, 000 generations.

Exploring population diversity for highly stable networks
In this set of experiments, we investigated how the population diversity is impacted in networks that have been exposed to many generations of strong selection for network stability (see Fig. S7). Specifically, we tested whether the increased compensatory mutation frequency shown in Fig. S2A was due to the property of particular networks that had been selected for, or whether it was the property of a diverse population.
Following the measurement used in Azevedo et al. [2006], the genetic diversity is defined as: where n is the total number of alleles, i.e., the unique values contained in the same site crossing all individual networks, and p i is the frequency of allele i. The genetic variation in a population is calculated as the mean gene diversity over non-zero sites of the interaction matrix for a given genotype W . Figure S7: Population diversity of highly stable networks. For each network size (N = 5, 10, 15, 20, 30 and 40) with network connectivity c = 0.76, we tested population diversity for 10, 000 networks that had been exposed to selection for network stability following up to fifteen rounds of mutation as described in Fig. S2A. The error bars represent 95% confidence intervals based on 100 independent runs.
From Fig. S7, we can see that networks that have been though many generations of selection for network stability can still maintain a high network diversity.

Location of compensatory mutations
In this set of experiments, we first sought to visualise locations at which the compensatory mutations are more likely to occur (see Fig. S4). To this end, in a set of compromised networks (those stable networks that proved fragile to a single round of mutation), we marked the site of the deleterious mutation, then measured the relative frequency of compensatory mutation that occurred at each possible site, including the site of the deleterious mutation, within this compromised network. For each possible site, we measured the outcomes over 1, 000 simulated mutations on that site (so that only the extent of regulation was mutated randomly, not the location).
To quantify the distance between deleterious and (potentially) compensatory mutation, we first define 'distance' as used in this paper. Suppose a given gene regulatory network, denoted as W , has two marked edges denoted as − − → AB (deleterious mutation) and − − → CD (compensatory mutation), where A, B, C and D represent different genes in W and − → · marks the edge direction. The distance between − − → AB and − − → CD can be calculated as where dis(A, C) is the fewest edges possible from A to C and path(A, C) includes the vertices on the shortest path between A and C in network W .
An example process of compensatory mutation in a gene regulatory network can be seen in Fig. S8. This stable network can be compromised by a single deleterious mutation (marked in red) and compensated by an additional mutation (marked in blue). According to Equation (4), the distance from deleterious mutation site −→ CA to compensatory mutation site − − → CE can be calculated as: Next, we compared the relative frequencies of compensatory mutation among gene networks whose marked edges (caused by additional mutation) were 0, 1, 2, 3, and 4 steps away from the deleterious mutation (see Fig. 3, main text, and also see Fig. S9). We also performed similar experiments for medium (N = 20) and large networks (N = 40), as shown in Fig. S10 and Fig. S11. Figure S8: An example process of compensatory mutation in a gene regulatory network. The initially stable gene network contains five genes: A, B, C, D and E. In the initial network (on the left side), each directional edge represents the strength (weight) of interaction between the linked two genes. The initial gene expression pattern is s(0) = (−1, −1, +1, +1, +1). In the compromised network (in the middle), a mutation occurs on −→ CA (indicated in red), which leads to the failure of stabilising the gene expression patterns (marked by dashed circles). In the compensated network (on the right side), the compromised network is fixed by an additional mutation that occurs on − − → CE (indicated in blue), reaching an equilibrium expression s EQ = (−1, −1, +1, +1, +1). Figure S9: The compensatory mutation location and distance distribution of all mutations relative to the original deleterious mutation sites (Small Networks). For initially stable networks with size N = 5 and connectivity c = 0.4, we first collected a pool of compromised networks with deleterious mutations after a single mutation round. We then forced second mutations, classifying these as being 0 (on the same site), 1, 2, 3 and 4 steps away from the original deleterious mutations. For each of these mutation-site-distance categories, we measured the probability that the mutation was compensatory (that it returned the network to stability), based on 10, 000 sample networks collected for each distance category as shown in the solid line. We also recorded the spatial distribution of second mutations (10, 000 sample networks) occurring randomly in those compromised networks with respect to their original deleterious mutation sites, shown in the dashed line. The error bars represent 95% confidence intervals based on 100 independent runs. Figure S10: The compensatory mutation location and distance distribution of all mutations relative to the original deleterious mutation sites (Medium Networks). For initially stable networks with size N = 20 and connectivity c = 0.2, we first collected a pool of compromised networks with deleterious mutations after a single mutation round. We then forced second mutations, classifying these as being 0 (on the same site), 1, 2, 3, 4 and 5 steps away from the original deleterious mutations. For each of these mutation-site-distance categories, we measured the probability that the mutation was compensatory (that it returned the network to stability), based on 10, 000 sample networks collected for each distance category as shown in the solid line. We also recorded the spatial distribution of second mutations (10, 000 sample networks) occurring randomly in those compromised networks with respect to their original deleterious mutation sites, shown in the dashed line. The error bars represent 95% confidence intervals based on 100 independent runs. Figure S11: The compensatory mutation location and distance distribution of all mutations relative to the original deleterious mutation sites (Large Networks). For initially stable networks with size N = 40 and connectivity c = 0.15, we first collected a pool of compromised networks with deleterious mutations after a single mutation round. We then forced second mutations, classifying these as being 0 (on the same site), 1, 2, 3, 4 and 5 steps away from the original deleterious mutations. For each of these mutation-site-distance categories, we measured the probability that the mutation was compensatory (that it returned the network to stability), based on 10, 000 sample networks collected for each distance category as shown in the solid line. We also recorded the spatial distribution of second mutations (10, 000 sample networks) occurring randomly in those compromised networks with respect to their original deleterious mutation sites, shown in the dashed line. The error bars represent 95% confidence intervals based on 100 independent runs.

Exploring the size of gene regulation on compensatory mutation frequency
In this set of experiments, we investigated effective changes in gene regulation associated with these mutations (see Fig. 3, main text, and also see Fig. S12). Specifically, we conducted experiments to measure the frequency of compensatory mutation when the second mutation had an additional weight added to it.
We studied a range of weight changes from (w = [−5, 5]) with a step size of 0.05. For each step size, we first performed one mutation round as usual on the initial population of stable networks, creating a subpopulation of 10, 000 compromised networks. Then, for these mutated networks we performed a second mutation round; however, this time instead of replacing one entry in the interaction matrix with N (0, 1), we added a fixed value w drawn from [−5, 5] to the original value of the randomly picked site. Then, we measured the frequency of second mutations restoring the network stability. We also performed similar experiments for medium (N = 20) and large networks (N = 40), as shown in Fig. S13 and Fig. S14. Figure S12: The influence of different intensities of gene regulations on the frequency of compensatory mutation (Small Networks). We first collected 10, 000 sample networks that had been made unstable by a single mutation from a pool of initially stable networks with N = 5 and c = 0.4. Then, we experimented with how a new mutation of varying intensities of gene regulation altered the chances of restoring gene stability. Specifically, we performed new mutations to those compromised networks with deleterious mutations by adding a weight from [−5, +5] (step size 0.5) to the original regulatory impact, then assessed the resulting patterns in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). The shaded areas represent 95% confidence intervals based on 100 independent runs. Figure S13: The influence of different intensities of gene regulation on frequency of compensatory mutation (Medium Networks). We first collected 10, 000 sample networks that had been made unstable by a single mutation from a pool of initially stable networks with N = 20 and c = 0.2. Then, we experimented with how a new mutation of varying intensities of gene regulation altered the chances of restoring gene stability. Specifically, we performed new mutations to those compromised networks with deleterious mutations by adding a weight from [−5, +5] (step size 0.5) to the original regulatory impact, then assessed the resulting patterns in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). The shaded areas represent 95% confidence intervals based on 100 independent runs. Figure S14: The influence of different intensities of gene regulation on frequency of compensatory mutation (Large Networks). We first collected 10, 000 sample networks that had been made unstable by a single mutation from a pool of initially stable networks with N = 40 and c = 0.15. Then, we experimented with how a new mutation of varying intensities of gene regulation altered the chances of restoring gene stability. Specifically, we performed new mutations to those compromised networks with deleterious mutations by adding a weight from [−5, +5] (step size 0.5) to the original regulatory impact, then assessed the resulting patterns in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). The shaded areas represent 95% confidence intervals based on 100 independent runs.
Exploring the distribution of regulation in initially-stable, compromised and restored networks In this set of experiments, we investigated the distribution of regulation in initially stable, compromised and restored networks (see Fig. S5). Specifically, we collected 10, 000 sample regulatory values each from edges of randomly generated stable networks, edges where deleterious mutations occurred (compromising network stability), and edges where compensatory mutations occurred (restoring previously compromised networks). We then measured their corresponding distributions, discriminating between self-and non-selfregulatory edges. We also performed similar experiments for medium (N = 20) and large networks (N = 40), as shown in Fig. S15 and Fig. S16. Figure S15: The distribution of regulation in initially stable, compromised and restored networks (Medium Networks). For randomly generated stable networks with N = 20 and c = 0.2, we collected 10, 000 sample regulations. We also collected 10, 000 sample regulation weights from deleterious mutations that compromised initially stable networks as well as from compensatory mutations that restored the stability of previously broken networks. We then measured the distributions in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). Given that the regulations are continuous values, we grouped them into 19 bins from [−4.5, +4.5] (step size 0.5). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S16: The distribution of regulation in initially stable, compromised and restored networks (Large Networks). For randomly generated stable networks with N = 40 and c = 0.15, we collected 10, 000 sample regulations. We also collected 10, 000 sample regulation weights from deleterious mutations that compromised initially stable networks as well as from compensatory mutations that restored the stability of previously broken networks. We then measured the distributions in all regulatory edges (A), in self-regulatory edges (B) and ignoring self-regulatory edges (C). Given that the regulations are continuous values, we grouped them into 19 bins from [−4.5, +4.5] (step size 0.5). The error bars represent 95% confidence intervals based on 100 independent runs.

Exploring properties of location and size effects in neutral mutations
In this set of experiments, we investigated properties of location and size effects in neutral mutations which served as control groups for solid lines in Fig. 3A and B, main text, Specifically, to test the location effect, we collected a population pool of stable networks that had been subjected to one round of mutation (neutral). Then, we measured the probability of stable networks after performing a second mutation that was 0, 1, 2, 3, and 4 steps away from the previous neutral mutation site based on 10, 000 sample networks for each distance category (see dashed line in Fig. 3A). Similarly, to test the mutation size effect, we collected a population pool of stable networks that had been subjected to one round of mutation (neutral). Then, we measured the probability of stable networks after performing a second mutation that had a particular shift in gene regulation from [−5, +5] based on 10, 000 sample networks for each shifted-weight category (see  For medium networks (N = 20, c = 0.2) and large networks (N = 40, c = 0.15), we first collected a pool of stable networks with neutral mutations after a single mutation round. We then forced second mutations, classifying these as being 0 (on the same site), 1, 2, 3 and 4 steps away from the previous neutral mutations. For each of these mutation-site-distance categories, we measured the probability that the mutation was neutral (did not impair network stability) based on 10, 000 sample networks collected for each distance category. The error bars represent 95% confidence intervals based on 100 independent runs. Figure S18: Mutation size effect in networks with neutral mutations (Medium and Large Networks). We first collected 10, 000 stable networks with neutral mutations after a single mutation round from a pool of initially stable medium networks (N = 20, c = 0.2) and large networks (N = 40, c = 0.15). Then, we experimented with how new mutations of varying intensities of gene regulation altered the chance of retaining network stability. Specifically, we performed new mutations to those networks with neutral mutations by adding a weight from [−5, +5] (step size 1 and with four additional regulation shifts: −0.5, −0.1, 0.1 and 0.5) to the original regulatory impact, then assessed the resulting patterns. The error bars represent 95% confidence intervals based on 100 independent runs.

Exploring the impact of distance and size effects on network robustness
In this set of experiments, we explored the effects of location and mutation size on robustness in networks with one deleterious mutation and one compensatory mutation and in networks with two consecutive neutral mutations to investigate whether networks with compensatory mutations have a different evolutionary consequence compare with networks with neutral mutations (see Fig. 3C and D and also see Fig. S19 and Fig. S22).
Specifically, to test the distance effect, we collected 10, 000 sample networks at each distance (between deleterious mutation and compensatory mutation). Then, for each category of distance, we measured the proportion of stable networks after one additional round of single mutation. The reported results are both actual robustness (see the solid line in Fig. S19A) and percentage change in robustness (see the solid line in Fig. S19B). Similarly, for the control group, instead of collecting networks that were subjected to one deleterious mutation and one subsequent compensatory mutation, we collected 10, 000 sample networks that were subjected to two consecutive neutral mutations at each distance (between two neutral mutations), and then assessed the actual robustness (see the dashed line in Fig. S19A) as well as the percentage of robustness change (see the dashed line in Fig. S19B). We also performed similar experiments for medium (N = 20) and large networks (N = 40), as shown in Fig. S20 and Fig. S21.
Likewise, to test size effect, we collected 10, 000 sample networks that were compensated by mutations with different shifts in gene regulation. Then, for each category of mutation size, we measured the proportion of stable networks after one additional round of single mutation. The reported results are both actual robustness (see the solid line in Fig. S22A) and percentage change in robustness (see the solid line in Fig. S22B). Similarly, for the control group, instead of collecting networks that were subjected to one normal deleterious mutation and one subsequent compensatory mutation with different shifts in gene regulation, we collected 10, 000 sample networks that were subjected to two consecutive neutral mutations, one normal neutral mutation and the other neutral mutation with different shifts in gene regulation, and then assessed the actual robustness (see the dashed line in Fig. S22A) as well as the percentage of robustness change (see the dashed line in Fig. S22B). We also performed similar experiments for medium (N = 20) and large networks (N = 40), as shown in Fig. S23 and Fig. S24. Figure S19: The impact of distance effect on network robustness (Small Networks). For small networks (N = 5, c = 0.4), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation that was 0, 1, 2 and 3 steps away from the previous deleterious mutation. The sample networks for control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Then, we assessed robustness of sample networks at each distance step. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then divided by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S20: The impact of distance effect on network robustness (Medium Networks). For medium networks (N = 20, c = 0.2), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation that was 0, 1, 2, 3 and 4 steps away from the previous deleterious mutation. The sample networks for the control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Then, we assessed the robustness of the sample networks at each distance step. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then dividing by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S21: The impact of distance effect on network robustness (Large Networks). For large networks (N = 40, c = 0.15), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation that was 0, 1, 2, 3 and 4 steps away from the previous deleterious mutation. The sample networks for the control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Then, we assessed the robustness of the sample networks at each distance step. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then dividing by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S22: The impact of mutation size effect on network robustness (Small Networks). For small networks (N = 5, c = 0.4), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation with different shifts in gene regulation from [−5, +5] (step size 1 and with four additional regulation shifts: −0.5, −0.1, 0.1 and 0.5). The sample networks for control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Note that the second neutral mutation has different shifts in gene regulation as the compensatory mutation. Then, we assessed robustness of sample networks at each category. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then divided by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S23: The impact of mutation size effect on network robustness (Medium Networks). For medium networks (N = 20, c = 0.2), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation with different shifts in gene regulation from [−5, +5] (step size 1 and with four additional regulation shifts: −0.5, −0.1, 0.1 and 0.5). The sample networks for the control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Note that the second neutral mutation has different shifts in gene regulation to the compensatory mutation. Then, we assessed the robustness of the sample networks at each category. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then dividing by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs. Figure S24: The impact of mutation size effect on network robustness (Large Networks). For large networks (N = 40, c = 0.15), we collected 10, 000 sample stable networks that were subjected one deleterious mutation and then restored by one subsequent compensatory mutation with different shifts in gene regulation from [−5, +5] (step size 1 and with four additional regulation shifts: −0.5, −0.1, 0.1 and 0.5). The sample networks for the control group were collected in a similar way, except that the networks were subjected to two consecutive neutral mutations. Note that the second neutral mutation has different shifts in gene regulation to the compensatory mutation. Then, we assessed the robustness of the sample networks at each category. The reported results are actual robustness (A), and change in robustness (B) (the actual robustness was normalised by subtracting the minimal value among all categories, and then dividing by the minimal value). The error bars represent 95% confidence intervals based on 100 independent runs.

Exploring how network connectivity evolves under a relaxed selection regime
In this set of experiments, we investigated whether regulatory complexity (increased network connectivity) could arise under a relaxed selection regime where compensatory mutations could occur and accumulate (see Fig. 4 and Fig. S6).
In the first set of experiments, we tested whether we could observe greater complexity arising using a population pool of 10, 000 stable networks of N = 10 genes with a simple 'Star' topology (see fig 4, main text). Specifically, the initial population pool was generated using the following rules: • Randomly select a gene to be the hub node.
• There is at least one edge between the hub node and non-hub nodes (either inward or outward); there is a possibility (0.5) of having both inward and outward edges.
• Each node has a possibility (0.5) of having a self-regulatory edge (including the hub node).
• The value (interaction strength) of each edge is drawn from the standard normal distribution N (0, 1).
In theory, for network size N = 10, the minimum connectivity is c min = 0.09 (9 edges) and the maximum connectivity is c max = 0.28 (28 edges). In the randomly generated initial population pool used in this paper, the minimum connectivity was c min = 0.10 (10 edges), the maximum connectivity was c max = 0.26 (26 edges), the median connectivity wasc = 0.17 (17 edges) and the average connectivity wasc ≈ 0.17. Then, the initial population was evolved for 5, 000 generations under strong and relaxed selection regimes: In four scenarios with selection for network stability, the initial population was evolved under: a no mutation and no recombination regime, a mutation but no recombination regime, a recombination but no mutation regime, a mutation and recombination; in three other scenarios, the initial population was evolved under a relaxed selection regime with a frequency of 1/10, 1/25, and 1/50. The statistical details for connectivity in initial and evolved populations can be found in Table S1. Note that compensatory mutation could only occur during periods of relaxed selection.
In order to further test the hypothesis that relaxed selection can facilitate regulatory complexity, in the second set of experiments, we further investigated how network connectivity evolves under a relaxed selection regime using randomly generated networks (see Fig. S6). Specifically, for a network size N = 40 with connectivity c = 0.15, we collected 10, 000 stable networks, each of which had the same initial gene expression pattern, all activation, i.e., s(0) = (+1, +1, . . . , +1). This population was then evolved for 5, 000 generations, in this case allowing for recombination with other individuals from the same generation. Note that in the previously-described experiments in this paper, a mutation could not change the topology of an individual network; that is, it could not change zero elements into non-zero or vice versa. In contrast,