The distribution of fitness effects during adaptive walks using a simple genetic network

The tempo and mode of adaptation depends on the availability of beneficial alleles. Genetic interactions arising from gene networks can restrict this availability. However, the extent to which networks affect adaptation remains largely unknown. Current models of evolution consider additive genotype-phenotype relationships while often ignoring the contribution of gene interactions to phenotypic variance. In this study, we model a quantitative trait as the product of a simple gene regulatory network, the negative autoregulation motif. Using forward-time genetic simulations, we measure adaptive walks towards a phenotypic optimum in both additive and network models. A key expectation from adaptive walk theory is that the distribution of fitness effects of new beneficial mutations is exponential. We found that both models instead harbored distributions with fewer large-effect beneficial alleles than expected. The network model also had a complex and bimodal distribution of fitness effects among all mutations, with a considerable density at deleterious selection coefficients. This behavior is reminiscent of the cost of complexity, where correlations among traits constrain adaptation. Our results suggest that the interactions emerging from genetic networks can generate complex and multimodal distributions of fitness effects.


Introduction
A lingering question in the study of adaptation concerns the distribution of effect sizes of adaptive alleles that affect trait and fitness variation in natural populations.Research across various species and traits has yielded mixed results.In some cases, natural selection favors alleles with small phenotypic effects (e.g.Visscher et al. 2012;Boyle et al. 2017;Kosheleva and Desai 2018), such as the alleles contributing to the evolution of body size in mammals (Kemper et al. 2012).Other times, however, alleles with larger effects are preferred, such as those found in flower color changes that affect pollinator preferences (e.g.Dorweiler et al. 1993;Bradshaw et al. 1995;Bomblies and Peichel 2022).Understanding when allelic effects of different sizes are favored is crucial to comprehend the "adaptive walk"the metaphorical path a population takes through genetic space as it adapts to its environment (Walsh and Lynch 2018b).Each "step" of the walk represents a beneficial allele fixing in the population (i.e.reaching 100% frequency), moving it closer to phenotypic states that maximize fitness (Maynard Smith 1970; Kauffman and Levin 1987).The distribution of fitness effects 20 (DFE) among beneficial alleles is a key feature of adaptation, 21 as it describes the tempo and mode of the adaptive walk in phe-22 notypic space: how many steps there are, how large each step 23 is, and in which direction each step drives the population.

24
Adaptive walks and the distribution of fitness ef-25 fects

26
The leading theory of adaptive walks was developed by Gille-27 spie and Orr (Gillespie 1983, 1984;Orr 1998).Under this the-28 ory, adaptive walks are characterized by a sequence of muta-29 tions with diminishing effect sizes.A key prediction of the 30 Gillespie-Orr model is that the distribution of beneficial fitness 31 effects (DBFE) during adaptation is negative exponential (we 32 will refer to this DBFE as exponential from here on) (Orr 1998(Orr , 33 2003a)).This shape emerges as populations approach the opti-34 mum phenotype: large effect mutations become progressively 35 disadvantageous due to the increased risk of "overshooting" the 36 optimum.Gillespie recognized that extreme-value theory -a 37 statistical branch focused on sampling from extreme tails of distributions -could be integrated into the genetic theory of adaptation to study the DBFE.To understand this connection, consider the distribution of fitnesses among all genotypes.This distribution is analogous to the distribution of fitness effects (DFE) among all new mutations.Under normal circumstances, the wild-type genotype should have relatively high fitness, and hence it belongs to the right tail of this distribution.The wildtype is a threshold which separates the deleterious genotypes to the left of the wild-type in the distribution from the beneficial right tail.This tail is analogous to the DBFE.Extreme value theory predicts that the DBFE can be described by a family of distributions termed the extreme value distribution (EVD).
Gillespie posited that adaptive steps are samples from the EVD (Gillespie 1984).Studying the shape of the EVD highlights the availability of beneficial mutations to a population and which effect sizes are likely to contribute to adaptation.
The EVD can be broadly categorized into one of three "domains of attraction", or shapes: the Gumbel, Weibull, or Fréchet family of distributions (Fig. 1; (Rokyta et al. 2008;Joyce et al. 2008;Martin and Lenormand 2008)).Gillespie argued that the Gumbel is most likely to represent empirical DFEs as it captures common distributions such as the normal, gamma, and log-normal (Gillespie 1983(Gillespie , 1984)).The Gumbel domain describes the DBFE's exponential expectation outlined above (Orr 1998).Mutagenesis studies and spontaneous mutation experiments largely support this theoretical expectation (e.g.Orr 1998;Imhof and Schlötterer 2001;MacLean and Buckling 2009;Kassen and Bataillon 2006;Sanjuán et al. 2004;Vale et al. 2012;Eyre-Walker and Keightley 2007).However, evidence in favor of alternative domains does exist, suggesting the field is still open to exploration.Adaptive walks with Weibulldistributed DBFEs are characterized by fewer large-effect beneficial mutations compared to Gumbel EVDs (Charras-Garrido and Lezaud 2013, Fig. 1).Rokyta et al. (2008) observed a Weibull distribution in the ID11 ssDNA phage, hinting at an upper limit on the size of beneficial fitness effects.Further, stabilizing selection can generate a Weibull-distributed EVD by limiting the size of beneficial mutations as populations approach the optimum (Martin and Lenormand 2008).On the other hand, Fréchet EVDs have more frequent beneficial mutations (Charras-Garrido and Lezaud 2013, Fig. 1).Schenk et al. (2012) observed a Fréchet EVD in Escherichia coli adapting to antibiotics.Although it is unclear as to which conditions cause Fréchet EVDs to arise, environments which invoke strong selective pressures on populations are associated with their appearance (Neidhart et al. 2014;Schenk et al. 2012;Bank et al. 2014;Foll et al. 2014).
There are a large number of factors which can influence the DBFE and these can produce non-Gumbel EVDs.Some of these factors are environmental, such as stabilizing selection driving Weibull EVDs (Martin and Lenormand 2008).However, developmental and selective constraints also contribute to the shape of the DBFE, implicating the structure of the genotype-phenotype-fitness (GP W ) map in the DBFE's shape (Connallon and

Genotype-phenotype maps in adaptation 58
The genotype-phenotype-fitness (GP W ) map is composed of 59 two distinct components: the genotype-phenotype (GP ) and 60 the phenotype-fitness (P W ) relationships.The GP map de-61 scribes how developmental and physiological processes trans-62 late genotypes into phenotypes, while the P W map describes 63 how ecological and selection regimes create fitness differences 64 between phenotypes.For example, stabilizing selection favors 65 intermediate phenotypes, while directional selection favors ex-66 treme phenotypes in one direction (Simpson 1944).The P W 67 map, along with the related genotype-fitness (GW ) map, has 68 been extensively studied in the fields of quantitative and popu-69 lation genetics (e.g.Wright 1931Wright , 1932;;Thornton 2019;May-70 nard Smith et al. 1985;Lande and Shannon 1996).Because the 71 GP map is notoriously challenging to estimate, many quan-72 titative genetics models assume an additive relationship be-73 tween genotype and phenotype (Walsh and Lynch 2018a;Lande 74 1976).These models derive from Fisher's (1918) infinitesimal 75 model, which supposes that continuous trait distributions can 76 be produced by loci under Mendelian segregation as long as 77 those loci are a) many, and b) have small, additive effects on 78 the phenotype (Fisher 1918).This model has become the basis 79 of modern quantitative genetics, providing a simple GP map 80 that requires no information about the underlying genetic sys-81 tems that describe the developmental and physiological under-82 pinnings of traits (Falconer 1996;Walsh and Lynch 2018a).
Overwhelming empirical evidence suggests that non-additive gene interactions (epistasis) are ubiquitous in nature (Álvarez-Castro and Carlborg 2007;Hansen 2013;Ang et al. 2023) and play a crucial role in adaptation (Draghi et al. 2011;Breen et al. 2012;Bendixsen et al. 2017;Weinreich et al. 2006).Additive models, including the infinitesimal, often fail to capture these gene interactions, and when identified, it is unclear how they represent biological gene action (Huang and Mackay 2016).
Hence, it is unknown how functional epistasis might contribute to the developmental constraints that limit evolution (Hansen 2013).As we delve deeper into the nature of the genotype-phenotypefitness landscape, it becomes important to re-evaluate our foundational quantitative genetic models in light of our current understanding of the molecular basis of traits.For instance, how does the distribution of fitness effects change when non-additive GP maps are considered?Does this affect a population's chance to adapt to a changing environment?To address these questions, we consider a simple gene regulatory network motif to motivate our approach and implement a nonlinear GP map into adaptive walk theory.
At the foundation of any trait lies a complex network of interacting genes and regulatory elements that respond to environmental cues by producing specific proteins (Doane and Elemento 2017).Systems biologists employ mathematical networks known as gene regulatory networks (GRNs) to model such systems (Alon 2019).Empirical networks harbor startling complexity: for example, consider the circadian clock network in Arabidopsis thaliana.This network is driven by a number of interacting network motifs: small, common subnetworks with particular effects on expression behavior.For instance, a feedforward loop motif in this network generates pulses of expression in PRR9/7 (Joanito et al. 2018).In addition to this behavior, a negative feedback loop (another motif) between CCA1/LHY and PRR5/TOC1 generates a bistable switch (Joanito et al. 2018).This switch is toggled at daybreak by the aforementioned PRR9/7 feedforward loop and again at dusk by yet another circuit (Joanito et al. 2018).Given the complexity of many GRNs, it is common to study motifs as separate units to elucidate the reasons for their recruitment into so many systems (Alon 2007;Hallinan and Jackway 2005;Seo et al. 2009).This allows systems biologists to probe the general effects of motifs on a broad range of networks.In this study, we focus on the simplest network motif, negative autoregulation (NAR; Fig. 2).This motif consists of two genes: gene X activates gene Z (indicated by the pointed arrow in Fig. 2A), while the presence of Z's product limits its further expression (indicated by the flat arrow in Fig. 2A).
NAR motifs are highly prevalent in biological networks.For The NAR motif consists of two genes, X and Z.The expression of X activates Z, and the expression of Z inhibits further Z production (A).This results in a characteristic expression curve (B).Z production begins when X is activated (blue shaded area) and stops when it is inactivated.Gene Z product approaches a steady state concentration and then quickly falls off in the absence of X.
example, close to half of all transcription factors in E. coli are 58 negatively autoregulated (Stewart et al. 2013;Shen-Orr et al. 59 2002).In plants, expression of the oil biosynthesis transcrip-60 tion factor WRINKLE1 (WRI1) is driven via a NAR network 61 that has been evolutionarily conserved from at least the split 62 between monocotyledons and dicotyledons (Snell et al. 2019).63 Snell et al. (2019) found that WRI1 expression was limited by 64 the concentration of WRI1 protein in planta.The ubiquity of 65 the NAR motif in nature comes as a result of its self-balancing 66 property: when a gene is overexpressed, the NAR mechanism 67 triggers to reduce that gene's production, leading to a steady 68 state of gene expression (Alon 2019).This reduces variabil-69 ity in gene expression between cells and accelerates transcrip-70 tional responses to environmental cues (Rosenfeld et al. 2002;71 Nevozhay et al. 2009).Owing to its simplicity and ubiquity, we 72 consider the NAR motif a reasonable toy model for beginning 73 to explore the evolution of complex traits mediated by genetic 74 networks.

75
How might we expect networks to affect adaptation?Genetic 76 networks impose functional epistasis, biological interactions 77 between genes, which can create nonlinear GP maps (Het-78 her and Hohenlohe 2014).This results in rugged fitness land-79 scapes where the presence of multiple local optima might ob-80 struct adaptation to a global fitness peak (Romero and Arnold 81 2009).In other words, populations can get trapped at local op-82 tima because the path towards the global peak (which maxi-83 mizes fitness) is paved by low-fitness intermediate genotypes 84 (Wright 1931).Adaptation should then be limited by the struc-85 ture of the fitness landscape: the more local peaks there are, 86 the less likely it is for a population to be able to find the global 87 optimum.In turn, the structure of the fitness landscape will 88 depend on the nature of the underlying network.For instance, 89 the NAR network provides a relatively simple fitness landscape.90 Kozuch et al. (2020) investigated the fitness landscape of E. 91 coli's lexA NAR network.The authors found that the fitness 92 landscape was ridge-like, with fitness maximized along a pa-93 rameter space that balanced lexA production with the strength 1 of autoregulation (Kozuch et al. 2020) The coefficients in this equation have biological relevance.X 52 and Z represent the cellular concentrations of the two genes.
53 n XZ and n Z are steepness coefficients (Hill coefficients) reflecting the sensitivity of the system to the presence of X and Z products.Higher values indicate a more rapid activation (n XZ ) or repression (n Z ) response as X or Z concentration increases (Alon 2019).α Z is the rate at which Z is removed from the cell.This might reflect, for example, the activity of ubiquitinase in tagging Z for removal and the 26S proteasome in breaking down the protein (Morimoto et al. 2016;Schmidtke et al. 2006).
β Z represents the production rate of Z, which is influenced by factors such as transcription factor binding affinity, enhancers, silencers, and trans-acting regulatory elements (Halfon 2020).These biological interpretations provide some realism to modeling quantitative traits.It also gives interpretability to the action of mutations that affect these coefficients, which will be explained in the sections following.
X follows a step function and is only expressed within the interval t ∈ [t start , t stop ]: where t start and t stop are the time points at which X is activated and deactivated.Values for these parameters are given below.
The solution of the ODE is characterized by a nonlinear approach to maximum Z expression, followed by a decline after X expression ceases (Figure 2B).The area underneath the expression curve (Figure 2B) is the total amount of Z produced, which we can take as a quantitative trait value.To create variation in this trait value, the coefficients of Eq. 1 can be varied.We refer to these coefficients as "molecular components".In this study, we focus on modeling mutations in two of the NAR's molecular components: the Z degradation rate, α Z , and the Z production rate, β Z .The constant values of the other molecular components (K Z , K XZ , n Z , and n XZ ) can be found in Supp.Table 1.Mutations in our model have direct effects on either α Z or β Z .This leads to a differently-shaped expression curve when the ODE is solved, and hence a (potentially) different trait value.
We adopt a model where loci contribute multiplicatively to the values of molecular components.We refer to these loci as molecular quantitative trait loci (mQTLs).The multiplicative transformation ensures that the molecular components are always positive, which is essential as these values represent rates which cannot be negative.The evolution of these mQTLs can be studied using the Wright-Fisher model, a foundational model in population genetics that describes the stochastic process of allele frequency change in a finite population (Wright 1931;Fisher 1930).In our implementation, individuals in the Wright-Fisher population are diploid.
For a more abstract model that applies to any network and any number of traits with any number of mQTLs, refer to S1 Appendix.

Phenotype calculation
To calculate an individual's phenotype from their set of mQTLs, we follow the below algorithm: 1. Take the exponent of the sum of all alleles for α Z across 1 all loci and both homologous chromosomes (this is anal-2 ogous to an additive model summing effects across all 3 loci, but with a multiplicative transformation as described 4 above).4. Take the area under the expression curve to get the total 10 amount of Z expression (the phenotype).
5. Repeat steps 1-4 for each individual in the population.
12 This algorithm is described mathematically below. 13 Let C be a vector of the molecular components: An allele at locus i and chromosome j has an effect a ij on one 16 of the molecular components. 17 Given an individual's alleles, the value for each molecular com- 21 where 35 where F is the function that defines the differential equation.

39
To calculate Z, we integrate the expression curve: where t max = 10.
In Eq. 2, t start = 1 and t stop = 6.Z activity is evaluated for

Fitness calculation
To map the trait value (P ) to fitness (w), we utilize a Gaussian fitness function, which allows us to model directional and stabilizing selection.The Gaussian fitness function, originally described by Lande (Lande 1976), is defined as follows: In this function, ∆z = P − P O , where P O is the optimal phenotype, and σ w is the width of the fitness function.Wider functions represent weaker selection.

Evolutionary simulation
To model the evolution of the network over multiple generations, we employ a Wright-Fisher (WF) model (Fisher 1923;Wright 1931).We consider a diploid population of N = 5000 individuals with random mating and non-overlapping generations.This satisfies the requirement that 2N µ 1 (where 2N is the total number of genomes in the population; 2N µ = 0.0915), which is necessary to meet the Gillespie-Orr model's assumption of strong selection relative to mutation (see section Model validation for more information on this) (Gillespie 1983; Orr 2010).Offspring are generated by sampling with replacement from parents, with the sampling probability weighted by their fitness.Individuals possess a pair of homologous chromosomes with two mQTLs: one for each molecular component.Random mutations occur at a rate of µ = 9.1528 × 10 −6 per locus per generation.This mutation rate is based on the average mutation rates observed in A. thaliana and adjusted to per-locus rates by multiplying it with the average length of a eukaryotic gene (1346bp) (Aston et al. 2017;Xu et al. 2006).A. thaliana was chosen as it is a model plant species with good estimates of mutation rates.The mutational effects of mQTLs on molecular components were drawn from a standard normal distribution and the effect added to the previous allele at that locus such that a new ∼ N (a old , 1), where a new is the new allele at a given locus and a old is the previous allelic effect at that locus.After being sampled, the allelic effects were exponentiated as per Eq. 3. We assume free recombination between loci.
We also considered an additive model.In this scenario, genomes also consisted of two causal loci to match the mutational target size of the NAR model.However, the trait value was given by summing allelic effects across QTLs instead of solving an ODE: where a ij k is the allelic effect of locus i on chromosome j.Both models underwent 50,000 generations of burn-in.The burn-in ensured that populations had high fitness just prior to the optimum shift, an assumption of the Gillespie-Orr model (Gillespie 1983(Gillespie , 1984;;Orr 1998).During burn-in, populations adapted to a phenotypic optimum at P O = 1.The phenotypic optimum was then instantaneously shifted to P O = 2, and the population was monitored for 10,000 generations of adaptation.

Computational implementation
We measured the phenotypic means and allelic effects of segregating and fixed mutations every 50 generations.We replicated each model 2,880 times with 32-bit integer seeds sampled from a uniform distribution in R 4.3.1 (R Core Team 2023).Detailed simulation parameters can be found in Table 1.Simulations were executed on the National Computational Infrastructure's Gadi HPC system.Of the simulation parameters, the most important were the fitness function width, σ w , and the optimum shift amount, O shif t .
σ w and O shif t are important because the Gillespie-Orr model assumes that the wild-type has relatively high fitness, so that in the total space of genotypes, it exists on the right-tail of the distribution (Orr 1998).Hence, the fitness function and optimum shift needs to be chosen so that at the optimum shift (generation 50,000) fitness is still relatively high.We chose σ w and O shif t so that this assumption is met: at the optimum shift, individuals perfectly adapted to the burn-in optimum suffered a ∼ 5% drop in fitness.

Model validation
The Gillespie-Orr model assumes that there is strong per-locus selection relative to mutation (Gillespie 1983(Gillespie , 1984)).This regime is often referred to as the strong selection weak mutation (SSWM) paradigm, as opposed to the weak selection strong mutation regimes found in polygenic models (e.g. the infinitesimal model) (Fisher 1918;Barghi et al. 2020).
To ensure we were in the SSWM domain assumed by adaptive walk theory, we measured population heterozygosity and the effects of segregating alleles on trait variation.To measure heterozygosity, H O , we used the equation where N het is the number of individuals heterozygous at locus i, and N is the total population size.Since we had two genes in our simulations (one for each molecular component, and two for the additive model), we took the average heterozygosity across both.To measure the effects of segregating alleles on the trait, we measured the ratio between the phenotype cre-56 ated by only fixations and the population mean phenotype, 57 r f s = P F ixed / P (8) 58 where P F ixed is the phenotype due to only fixed alleles and P 59 is the mean population phenotype.P F ixed was calculated by 60 removing any segregating effects from the α Z and β Z values 61 and recalculating the phenotype via Eq. 4. 62

Fitness effect calculation 63
An important part of adaptive walk theory is the shape of the 64 DFE.To assess this, we needed to measure the fitness effects of 65 mutations that arose during the adaptive walk.We performed 66 single gene knockouts and measured the difference between the 67 "wild-type" and knockout phenotypes and fitness.A flowchart 68 is given in Supp.Fig. 2A.To measure a given mutation's 69 effect on fitness, we first subtracted that mutation's effect on 70 the molecular component (α Z or β Z ) and recalculated the trait 71 value via Eq. 4. We then calculated the fitness of this knocked-72 out individual via Eq. 5. To get our final fitness effect, we sub-73 tracted the fitness of the individual with the mutation knocked 74 out (w k ) from the fitness of a "wild-type" individual (w): We built a custom tool to recalculate the phenotype and fitness 77 recalculations, implementing the Ascent C++ library to solve 78 ODEs (Berry et al. 2021).This library was also used in the 79 SLiM simulations, keeping the phenotype calculation consis-80 tent.The source code for our tool can be found at https: 81 //github.com/nobrien97/odeLandscape.The fit-82 ness effects were calculated in a similar manner for additive 83 simulations.However, the phenotype was computed as the sum 84 of additive effects instead of using the network ODE solution.

85
In silico mutant screen 86 In addition to the DFE of mutations that arose during the adap-87 tive walk, we were also curious about the total DFE across 88 the space of possible mutations (including deleterious muta-89 tions that would be lost quickly).We ran an in silico mutant 90 screen experiment to determine whether the fitness distribution 91 of possible mutations differed between additive and NAR pop-92 ulations.A flowchart is given in Supp.Fig. 2B.We sampled 93 1,000 mutations from a standard normal distribution and added 94 them to the molecular component values/quantitative trait val-95 ues of individuals at each step of the adaptive walk.This was 96 done for both network and additive models.We chose a stan-97 dard normal distribution as it matched the distribution of new 98 mutations that populations faced in the simulations.For the 99 network models, mutations were applied separately to α Z and 100 β Z , while in the additive models, the mutation effect was added 101 to the quantitative trait value.We then recalculated the phe-102 notypic and fitness effects of the sampled mutations using the 103 fitness effect calculation methods described above.The source 104 code for this experiment can be found at https://github.105 com/nobrien97/NARAdaptiveWalk2023.

Expected waiting times to beneficial mutation
To tie the mutant screen findings to a quantitative measure of adaptation, we used the mutant screen fitness effect distribution to estimate the expected waiting time to a beneficial mutation between models.We first calculated waiting times, t wait for a particular simulation replicate at a given adaptive step using the equation where N is the population size, µ is the per-locus, pergeneration mutation rate, and p s>0 is the probability that a new mutation is beneficial.N is multiplied by 4 because the population is diploid, and because there were two loci contributing to the trait.p s>0 is given by where n ISM S is the total number of mutations generated in the in silico mutation screen experiment and n ISM S s>0 is the number of mutation screen alleles with beneficial effects on fitness.
Waiting times were calculated for each adaptive step and model.
Comparisons between models were done using a bootstrap analysis.We sampled 100,000 different NAR-Additive model pairs across all adaptive steps, calculating the difference between their waiting times.We then calculated the means and 95% confidence intervals (CIs) for each model's waiting times and the difference between models (∆t wait ).The source code for this analysis can be found at https://github.com/nobrien97/NARAdaptiveWalk2023.A flowchart describing this methodology is shown in Supp.Fig. 3.
We also measured the difference in p s>0 between models and across adaptive steps using a linear model: Where α is the intercept, β i is the slope coefficient for depen-31 dent variable i, and is the residual error.To estimate marginal 32 means and contrasts, we used the R package emmeans 1.8.5 in 33 R 4.3.1 (Lenth 2023; R Core Team 2023).

34
Fitness landscapes

35
With the mutant DFE evaluated, we turned to how the underly-36 ing network might drive the complex distribution that we dis-37 covered.This involved constructing a fitness landscape to show 38 how fitness changed with different combinations of α Z and β Z , 39 introducing epistasis to the system.We generated 160,000 com-40 binations of α Z and β Z , sampling α Z and β Z values between 41 0 and 3. We then calculated the phenotype and fitness for each 42 combination using the ODE landscaper tool we used to calcu-43 late fitness effects.We plotted the resulting fitnesses against 44 α Z and β Z using ggplot2 3.4.2(Wickham 2016) in R 4.3.1 45 (R Core Team 2023).A similar method was used to estimate 46 the fitness landscape of the ratio β Z /α Z as α Z increased.For 47 this analysis, we generated another 160,000 combinations of 48 α Z and β Z /α Z ratios.Both α Z and β Z /α Z values were sam-49 pled between 0 and 3.

51
The final part of our analysis involved measuring the shape 52 of the distribution of beneficial fitness effects (DBFE), which 53 should be Gumbel-distributed under Gillespie-Orr predictions 54 (Orr 2003a;Gillespie 1984).We used a method developed 55 by Beisel et al. (2007) to fit a generalized Pareto distribution (GPD) to the DBFE obtained from the mutant screen.A flowchart describing this approach is provided (Supp.Fig. 4).
The GPD characterizes a distribution of extreme events that exceed a threshold.The threshold is the high-fitness wild-type and the extreme events are beneficial mutations.The shape parameter of the GPD, κ, (shape) indicates which domain of attraction the underlying distribution belongs to.Beisel et al.'s (2007) method uses a likelihood ratio test to evaluate if κ = 0, corresponding to the Gumbel case.If κ > 0, the DFE belongs to the Fréchet domain, and if κ < 0, it belongs to the Weibull domain (Fisher and Tippett 1928;Beisel et al. 2007).
We used two different sampling approaches to fit our data to GPDs (Supp.Fig. 4).In the first method, we pooled mutant screen alleles across all replicates to create a pooled DBFE (pDBFE) for additive and network models.We sampled 1,000 mutations from the pDBFE and fit the GPD to the resulting sampling distribution.We repeated this process 10,000 times, yielding a distribution of estimates of κ.The pooled approach shows the average shape of the DBFE across many adaptive walks.In the second method, we sampled up to 100 mutations from each simulation's mutant screen DBFE and fit the GPD to that sample.This was done per adaptive step.Simulation-step pairs with fewer than 10 beneficial alleles were discarded as they could not be reliably used to estimate the GPD shape (Beisel et al. 2007).We obtained a distribution of κ estimates, one for each simulation-step pair.This approach shows the variability in the shape of the DFE during different adaptive walks.In both methods, p-values from the likelihood ratio tests were combined across replicates using Fisher's method (Fisher 1936).Both approaches fit the GPD using the R package GenSA

Assumption validation
We first checked if we were in the strong selection weak mutation (SSWM) regime assumed by the Gillespie-Orr model.
Under SSWM, each beneficial mutation should fix before another arises.Hence, heterozygosity at QTLs/mQTLs should be close to zero most of the time.We found observed heterozygosity (H O ) at the QTLs/mQTLs was low on average, peaking just after the optimum shift at 11.39% ± 0.423% (95% CI) and 9.637% ± 0.400% (95% CI) in additive and NAR models, respectively (Supp.Fig. 8).The grand mean H O across all time points was 8.368% ± 0.028% (95% CI) and 7.841% ± 0.027% (95% CI) for additive and NAR models.Although H O was not zero, the majority of variation arising was not beneficial (Fig. 5B), meaning that the chance that two beneficial mutations could co-segregate was lower than H O suggests.
To further test the SSWM assumption, we also measured the ratio between phenotypes generated by only fixations (i.e. the phenotype in absence of all segregating variants) and the mean population phenotype (Eq.8.This ratio describes how much trait variation was contributed by segregating variants compared to fixations.When r F S = 1, segregating variants contribute no trait variation.When r F S > 1, segregating variants decrease the trait value on average compared to the phenotype due to 57 fixations, and vice versa when r F S < 1.We found that most 58 populations had no variation contributed by segregating vari-59 ants (Supp.Fig. 9).The average r F S was 1.019 ± 0.007 (95% 60 CI) and 1.024 ± 0.011 (95% CI) for additive and NAR popu-61 lations respectively.From these figures, it seems that the pop-62 ulations are in the SSWM domain, enabling us to explore the 63 extreme value theory expectations set out by Gillespie and Orr 64 (Gillespie 1983(Gillespie , 1984;;Orr 1998).

65
NAR and additive populations differ in their abil-66 ity to adapt 67 We first examined the phenotypic response to selection during 68 adaptive walks in populations with either an additive or a net-69 work (NAR motif) GP map.We defined an adapted popula-70 tion as one with a mean phenotype within 10% of the optimum 71 (1.9 < P < 2.1).We chose 10% as the threshold as it only 72 decreased fitness by 0.04% relative to a 5% threshold, but in-73 creased the sample size of adapted populations by 42% (2636 74 adapted populations vs. 1529).By the end of the simulation, 75 50.5% of the additive simulations and 41% of the network sim-76 ulations had adapted.For the remainder of the results, we con-77 sider only the adapted populations.

78
Among the adapted populations, there were differences in the 79 rate of trait evolution between the models (Fig. 3).Network 80 model populations took longer to reach the optimum on average 81 (Fig. 3A), although the changes in phenotype across adaptive 82 steps were similar between the models (Fig. 3B).Note that 83 this difference cannot be solely attributed to the NAR ODE, as 84 the multiplicative scaling of the allelic effects also contributes 85 to differences between the models.See the discussion for more 86 on this distinction.There was no difference between the models 87 in the number of adaptive steps taken.Of the populations that 88 adapted, 66.3% did so in a single step, 28% in two steps, and the 89 remaining 5.7% adapted in three or more steps.On average, the 90 number of steps taken was 1.4 ± 0.026 (95% CI).This agrees 91 with results from Orr (Orr 2003b), who found a mean walk of ≈ 92 1.72 steps, assuming that the genotype was comprised of many 93 loci and that the starting genotype of the walk was randomly 94 sampled.

95
Among populations that had not yet adapted after one step (the 96 left tail in Fig. 3), the median step in a network model occurred 97 900 (95% CI [150,1600]) generations later than in the additive 98 models (Supp.Fig. 5).To investigate the underlying causes 99 of the differences in adaptation success between models, we 100 examined the DBFE among fixations during the adaptive walk.101 The distribution of beneficial fitness effects among 102 fixations is similar between NAR and additive 103 populations 104 Fixations had similar effect sizes in both models and became 105 smaller over the course of the adaptive walk (Fig. 4).We found 106 that the mode of the distribution of beneficial fixations was off-107 set from 0 (Supp.Fig. 6), supporting a gamma distribution of 108 fixations rather than an exponential distribution.We fit gamma 109 distributions to the DFE among fixations in both the additive 110 and network models using fitdistrplus 1.1.8in R 4.3.1 111 The grey area shows the overlap between additive and network distributions.
(R Core Team 2023; Delignette-Muller and Dutang 2015).The deleterious fixations shown in Fig. 4 were excluded from the dataset when fitting the gamma distribution.These deleterious fixations were neutral or adaptive prior to the shift in the optimum and typically reached high frequencies before the adaptive walk commenced (Supp.Fig. 7).Therefore, these fixations can be attributed to the action of genetic drift.
We compared the shape and rate parameters of the gamma fits between the models.The network model populations exhibited slightly fewer small effect fixations and had a shorter tail compared to the additive models (NAR: Shape = 2.032 ± 0.147 95% CI, Rate = 44.392± 3.651 95% CI; Additive: Shape = 1.722 ± 0.108 95% CI, Rate = 37.588 ± 2.730 95% CI).
However, the DFE among beneficial fixations did not differ strongly between the models.To investigate this further, we needed to determine the availability of beneficial mutations to the populations.We did this by conducting an in silico mutant screen experiment, examining the DFE among new mutations at each adaptive step.
The distribution of fitness effects among possible 20 mutants is complex and deleterious in NAR pop-21 ulations

22
The in silico mutation screen experiment revealed that the DFE 23 among new mutations differed dramatically between models 24 (Fig. 5A).Both were strongly negatively skewed; however, mu-25 tations in NAR populations were overall more deleterious than 26 those in additive models, with a longer and more pronounced 27 tail of strongly deleterious mutations.In addition, the DFE in 28 NAR populations was bimodal (Fig. 5A).The first mode was at 29 s ≈ 0, similar to the additive DFE mode (NAR: s = −0.0003;30 Additive: s = −0.0004).The second mode was at s = −0.113.31 Among beneficial mutations, the distributions were more simi-32 lar.

33
Data pooling affects generalized Pareto distribu-34 tion fitting results

35
We fit a generalized Pareto distribution (GPD) to the DBFEs 36 generated by the mutant screen experiment using methods de-37 sampling each replicate at each of its adaptive steps separately and estimating the GPD for each estimate (for more information, see the Methods section).Under the pooled approach, we found evidence for a minor deviation from an exponential distribution towards the Weibull domain in both the additive and network models (Network: κ = −0.130± 0.0007; Additive: κ = −0.118± 0.0005).However, this deviation is small enough that the adaptive walk should behave as if its DFE belonged to the Gumbel domain (Joyce et al. 2008).
Under the per-simulation method, we saw a different result.
The DBFE in both NAR and additive populations were strongly Weibull distributed.In the NAR model, κ tended to decrease over the adaptive walk.However, the additive model κ estimate was relatively stable at around κ ≈ −2, slightly increasing over the adaptive walk (Supp.Table 2; Supp.Fig. 11).In ad-22 dition, these non-pooled fits were considerably more variable 23 than under the pooled approach (Supp.Fig. 10), particularly 24 at early adaptive steps.Under the per-simulation/adaptive-step 25 approach, across all adaptive steps, at least 99.8% of additive 26 replicates and at least 99.7% of NAR replicates had DFEs of 27 beneficial mutations that fit the Weibull domain.Under the 28 pooling method, all replicates in both models were approxi-29 mately Gumbel-distributed (i.e.κ ≈ 0).3.630 ± 1.107% (95% CI) less likely to produce a beneficial mutation.At further adaptive steps, there was no difference between models.We measured the effect of the lower beneficial mutation rate on adaptation by estimating the waiting time for a new beneficial mutation to arise.We found that the mean expected waiting time for a beneficial mutation in NAR populations was 98.555 ± 5.816 (95% CI) generations longer than in additive populations (Additive mean waiting time: 136.382 ± 2.131 (95% CI) generations; NAR mean waiting time: 234.937 ± 5.414 generations; Fig. 5C).Given that the DFE differed between models, we wanted to explore the causes underpinning this.To do so, we investigated the fitness landscape of α Z and β Z .
The NAR fitness landscape is ridge-like across molecular component space We found a pronounced fitness ridge on the NAR fitness landscape, surrounded by valleys of low fitness (Fig. 6A).This ridge lay diagonally across α Z and β Z space, suggesting that fitness depended on the ratio of the molecular components.This relates to the underlying ODE: β Z /α Z is the level of steady state Z expression under simple gene regulation (i.e.without any negative feedback or induction) (Alon 2019).Plotting this ratio over α Z confirmed this: β Z /α Z was the arbiter of fitness as long as α Z 0.5 (Fig. 6B).We investigated how the β Z /α Z 24 ratio contributed to the shape of the GP map, finding a non-25 linear relationship between β Z /α Z and phenotype (Fig. 6C).26 Around the optimum, increasing β Z /α Z quickly increases the 27 trait value, however with larger β Z /α Z increases, the trait value 28 changes less (Fig. 6C).We approximated the optimum β Z /α Z 29 ratio (θ β Z /α Z ), by substituting our fixed molecular components 30 We considered 31 the case where X = 1, simplifying to The equilibrium fulfills the equation If we assume a model of simple gene regulation, where θ is the steady state of gene expression under simple gene 40 regulation (Alon 2019).For small values of Z (where Z < 1), 41 this is well approximated by Using this, we approximate the solution of the ODE by the equi-3 librium whenever X = 1 and by 0 when X = 0. Hence, the 4 solution can be given by When the optimum is at P = M = 2, the optimum steady state is This assumes that the increase in Z after X is activated is rapid, and likewise for the decrease in Z when X is deactivated (i.e. in Fig. 2B, the green curve should be close to rectangular around X activation).This behavior is driven by n Z and n XZ : for our chosen values, the above approximation predicts the observed To explore how the nonlinearity of the GP map might have 3 influenced how populations navigated molecular component 4 space, we examined the ratio of absolute additive effect sizes 5 on molecular components among populations that adapted by 6 two steps or more, denoted by φ (Fig. 6D).Mathematically, where a is the allelic effect on α Z or β Z without the multiplica- tive transformation (i.e. the value composed by Eq. 6 instead 2007).We suggest that researchers should be cautious of the homogenization effect of the pooling approach when estimating the DBFE.Given these considerations, we next discuss the variation in estimates of the shape parameter κ when using pooled and non-pooled approaches.
When we fit GPDs to each replicate with their own set of mutations, almost all the fitted GPDs belonged to the Weibull domain (Supp.Fig. 10).In addition, κ was quite variable between simulations, models, and adaptive steps.κ affects the GPD by modulating an upper bound on the effect size of beneficial mutations.With decreasing κ, this upper bound approaches zero.
We observed that κ tended to decrease over longer walks under the network model.In the additive model, κ was relatively stable, however it slightly increased (particularly when comparing κ before the optimum shift and at adaptive steps ≥ 3 (Supp.Fig. 11; Supp.Table 2).Regardless of this increase, both models were largely Weibull-distributed (over 99% of simulations/adaptive steps belonged to the Weibull domain), which rejects Gillespie's Gumbel expectation.
Weibull-distributed DFEs of beneficial effects have been observed in a number of empirical studies in viral and bacterial populations (e.g.Rokyta et al. 2008;MacLean and Buckling 2009;Foll et al. 2014).Furthermore, a Weibull distribution of effects is predicted when populations are close to an optimum (Martin and Lenormand 2008).Hence, our results match theoretical expectations.Evidence for limits on adaptation seems common and both models seem to be affected by such limitations.While both models were starved of large-effect beneficial mutations relative to the exponential expectation, this effect was more limiting for network models than additive as the walk went on.This was further highlighted by a major difference between network and additive responses to selection: the rate of trait evolution.

Trait adaptation and network evolution
In our simulations, we found that network simulations approached the optimum more slowly than additive models on average (Fig. 3A).This difference in rate might not have only been due to the NAR structure, but also because of the difference in allelic effect scaling between additive and network models, as shown in Eq. 3 and Eq. 6.This scaling difference can influence rates of adaptation, mandating caution in attributing observed differences solely to the network's complexity.Hence, we do not make conclusions on how much the observed differences are due to the NAR network's structure compared to the difference in allelic effect scaling.Nonetheless, there was a difference in the rate of adaptation between the models and the NAR DFE does not completely reflect the expectation under a multiplicative model, so the distribution of effects is at least partly mediated by the NAR structure.This will be the subject of a future study.
We explored the distribution of new mutations and discovered a complex DFE in the network models, contrasting with the simple distributions observed in additive models (Fig. 5).The additive model distribution largely met Gillespie's expectations: the DFE was comprised of mostly neutral mutations with a small, exponential tail of beneficial mutations, and a longer tail of deleterious mutations (Fig. 5, (Gillespie 1991) pg.267, Fig. 1 6.5). 2 In the network model, we found a bimodal distribution of fit-3 ness effects.The bimodality is more pronounced than expected 4 compared to a multiplicative model, suggesting that the NAR 5 contributes to this shape.This observation aligns with previous 6 studies that have reported complex and multimodal distribu-  netic combinations have greater effects on fitness than the sum of their parts (positive epistasis), or the combinations can reduce fitness compared to the sum expectation (negative epistasis).Drift and mutation generate positive and negative epistasis in similar proportions (Ortiz-Barrientos et al. 2016).However, when recombination is rare, selection will quickly fix synergistic combinations, leading to an excess of deleterious gene combinations which persist through drift and limit the efficiency of selection (Hill and Robertson 1966).When a gene network introduces further epistasis on traits, this might compound the imbalance between negative and positive epistasis.Thus, network complexity may trade off with adaptive potential.
Trade-offs between complexity and adaptability echo the cost of complexity described by Orr (Orr 2000).The cost of complexity postulates that an increasing number of traits under selection leads to an increasingly deleterious mutation space, impeding adaptation.This cost could also apply to network models of traits, where molecular components under selection limit adaptation.This notion finds support in studies that have explored network topology in different species.For instance, Costanzo et al. (2010) found that increasing the degree of genetic interaction in Saccharomyces cerevisiae genes was negatively correlated with single mutant fitness, suggesting an adaptive cost to maintaining a highly connected gene network.This cost might lead to the long-term stability of gene networks or complex traits: comparisons between S. cerevisiae and Schizosaccharomyces pombe genetic interaction networks revealed remarkable similarity in network structure (Koch et al. 2012).
Another study by Barua and Mikheyev (Barua and Mikheyev 2021) found features of housekeeping gene networks involved in reptile venom production were conserved between amniote clades, with those genes contributing to saliva production in mammals (Barua and Mikheyev 2021).The cost of complexity suggests some rather strong limitations on the evolution of complex genetic networks: in Orr's (Orr 2000) model, the rate of adaptation declines according to the inverse square root of the number of traits, 1/ √ n.Empirical evidence suggests that biological networks can be extremely complicated, connecting hundreds of nodes.Adaptation involving such a highly connected network would be extremely slow if Orr's (Orr 2000) theory is correct.However, aspects of network structure, including local connectivity and modularity, might alleviate such a cost of molecular complexity.
One empirical example of network structure modulating the rate of evolution is the biofilm production network in Candida species (Mancera et al. 2021).In this network, seven master regulators control the expression of about one-sixth of the C. albicans genome and are required to drive biofilm production (Mancera et al. 2021).Genes that are connected to one or more of these master regulators are much less connected than the master regulators themselves and more "free" to evolve without affecting the expression of hundreds of other genes.Mancera et al. (2021) found that the master regulators evolve slowly compared to the target genes, possibly due to the structure of the network.Since the master regulators affect the expression of hundreds of genes, the space of beneficial mutations at master regulator loci might scale according to the cost of complex-ity.On the other hand, the authors found that target genes are considerably more divergent between Candida species and populations, perhaps due to target genes harboring fewer interactions within the network than master regulators (Mancera et al. 2021).
The NAR motif serves as a simplified regulatory network, effectively sidestepping the full cost of complexity outlined by Orr (Orr 2000).Despite this, it seems to impose greater adaptive constraints than the additive model, a point underlined by the density of strongly deleterious alleles in the DFE of new mutations (Fig. 5).Fixing multiple NAR molecular components (e.g.K Z , K XZ , n Z , n XZ ) renders the model's behavior more predictable.For sufficiently large α Z , individual fitness correlates solely with the ratio β Z /α Z (Fig. 6B), reflecting the steady-state cellular concentration under simple regulation (Alon 2019).This outcome may be attributed to our parameter choices; specifically, the rapid onset of Z activa- Should K Z be allowed to evolve, the fitness landscape would likely exhibit increased ruggedness.In this context, the optimal β Z /α Z ratio would be contingent upon the K Z value.
A similar effect is anticipated for K XZ , the regulator of Z's production rate in response to X. Consequently, both steadystate and fitness would depend on K Z and K XZ .Introducing more evolvable molecular components is expected to amplify the landscape's ruggedness due to interdependencies that influence fitness, exacerbating the previously noted cost of complexity.To substantiate these projections, future research should explore more intricate networks, as exemplified by empirical studies such as Bertheloot et al. (2020).In addition, exploring the behavior of other network motifs will elucidate the generality of our findings.
There are other motifs that are common in nature that could be reasonable choices for a toy model of network-mediated traits.
Other options include feedback loops, feedforward loops, single input modules, and cascades.Each of these motifs have different behaviors that could lead to different evolutionary dynamics.For instance, feedback loops can create oscillatory expression patterns (Pigolotti et al. 2007).Autoregulation is a simple example of a feedback loop, however multi-locus feedback loops also exist.The p53-Mdm2 feedback loop is an example of a two-gene feedback loop in vertebrates.p53 is an important tumor suppressor transcription factor.Its expression is kept at an equilibrium cell concentration by the presence of Mdm2 (Wu et al. 1993).p53 is a positive regulator of Mdm2, whilst Mdm2 inhibits p53 expression.When this balance is interrupted, increased Mdm2 production results in cell proliferation, whilst increased p53 production leads to cessation of growth (Wu et al. 1993).Another common motif, the feed-forward loop, has eight different configurations, each with dif-58 ferent behaviors (Alon 2019).Some of these configurations are 59 more common in transcription networks than others, suggesting 60 a selective pressure favoring these ubiquitous forms (Mangan 61 and Alon 2003).

62
Ramifications of strong selection weak mutation 63 (SSWM) and alternative approaches 64 This study focused on adaptation by employing a Gillespie-Orr 65 model, which operates under the assumption of SSWM, char-66 acterized by N s 1 and N µ 1.We ensured adherence 67 to the SSWM assumption by choosing an appropriate popula-68 tion size and mutation rate (see Table 1).In our simulations, 69 N µ = 0.092.However, N s varies due to selection coefficients 70 being dependent on individuals' genetic backgrounds and phe-71 notypes.Nonetheless, N s > 1 holds at s > 0.0002 for this 72 model, so most beneficial mutations should meet the N s 1 73 criterion.In the instance they do not, these mutations are driven 74 by drift (Lynch 2010).This also applies to slightly deleteri-75 ous mutations.N |s| = 1 represents a "drift barrier", where the 76 fitness effects of alleles become effectively too small to drive 77 allele frequency changes, leading to drift-dominated dynamics.78 The drift barrier explains our finding that the DFE among fixa-79 tions followed a gamma distribution as opposed to an exponen-80 tial in both network and additive models (6).As populations 81 approach the optimum, selection's efficiency is weakened by 82 the drift barrier (Lynch 2010).At the barrier, the probability of 83 fixation (π) is no longer predicated by π ≈ 2s when the pop-84 ulation is close to the optimum.Instead, π = 1/2N becomes 85 more applicable (Charlesworth and Charlesworth 2010).The 86 switch from selection-dominated to drift-dominated dynamics 87 occurs when N |s| = 1 (Koonin 2016).In our simulations, this 88 means that when |s| < 0.0002, drift is expected to dominate.89 Owing to the relatively large population size in our simula-90 tions (N = 5000), drift is unlikely to drive allele fixations early 91 in the adaptive walk (i.e. when most beneficial alleles have 92 s > 0.0002), however later steps might be more affected as the 93 possible fitness improvement of any given allele faces dimin-94 ishing returns.

95
Another assumption of the SSWM model is that only one muta-96 tion segregates at any given time point, i.e. before a new benefi-97 cial mutation arises, the previous beneficial mutation must have 98 fixed.This was met most of the time (see Supp.Fig. 8, 9), 99 however, some simulations revealed instances where segregat-100 ing variation contributed to adaptation, including the presence 101 of alleles under balancing selection (Supp.Fig. 12).Future re-102 search should therefore explore the model using a polygenic ap-103 proach, where many loci contribute small amounts to the molec-104 ular components underlying complex traits.Adaptation would 105 then occur by small frequency shifts at these loci (Walsh and 106 Lynch 2018c).The genetic architecture of polygenic traits, that 107 is, the number of loci, frequencies, and the distributions of their 108 additive and non-additive effects on phenotypes, can greatly af-109 fect the response to selection (Barghi et al. 2020).There are 110 also expectations for how genetic networks might determine 111 the genetic architectures of polygenic traits: large-effect loci 112 are likely linked to cis-regulatory features and coding region 113 mutations, whilst small-effect loci are more likely to represent trans-regulatory elements (Boyle et al. 2017;Liu et al. 2019).
However, the non-additive effects of genetic networks on phenotypes remain understudied.These interactions have ramifications for theories of the evolution of recombination (Ritz et al. 2017), and the persistence of linkage disequilibrium in natural populations (Ortiz-Barrientos et al. 2016).

Opportunities and challenges in modeling complex genetic networks
Our study serves as an introductory venture into the evolutionary modeling of genetic networks, elucidating how they guide evolutionary trajectories.The strong selection and weak mutation assumptions grant computational tractability but sideline rich polygenic dynamics (Barghi et al. 2020;Thornton 2019).
The NAR motif offers initial insights but lacks empirical complexity, making future work with more complex architectures important.Moving forward, our focus will shift towards a spectrum of genetic architectures to better capture non-additive effects.This will include investigating other network motifs, including several feedforward loops.By doing so, we will enrich our understanding of how topology influences evolutionary pathways (Siegal et al. 2007).This transition will require further statistical and computational innovations to maintain feasibility.Collaborations with molecular biologists will enrich our models with further mechanistic details.Rigorous validation against empirical data and existing theories will be a cornerstone in constructing more realistic models of adaptation (e.g.Bertheloot et al. 2020).
As well as this simple approach, we also constructed a general form of our network model in S1 Appendix.In the full model, we envision a phenotype as a combination of "molecular traits": intermediate traits like methylation or gene expression that represent cellular or developmental processes (Claringbould et al. 2017).Our NAR model represents the simplest case, where there is one molecular trait, Z expression, which we treat as the phenotype.In more complex implementations, phenotypes can emerge from a composite of molecular traits.Each molecular trait has its own differential equation and its own set of molecular components driving gene expression, creating a hierarchical model of gene expression driving quantitative trait variation.
As we delve into more complex systems, computational hurdles become prominent.The complexity of the ODE system significantly impacts in silico performance (Agocs 2020;Städter et al. 2021).We are considering optimization strategies such as smarter caching of common genotypes, data-oriented design, and adaptive step-size algorithms.Pre-computed ODE solutions stored on disk are also under exploration as a means to tackle the computational load.In summary, our study, while preliminary, opens up avenues for understanding the complex interplay between genetic networks and evolutionary dynamics, leaving us optimistic about the future of this research area.

Conclusion
Our network approach contributes methodologically to the field and holds broad applicability to many open questions in evolutionary biology, from quantitative genetics to the mechanisms underpinning adaptability.It offers a robust framework to in-56 vestigate issues such as the maintenance of genetic variation, 57 the evolution of recombination, and the dynamics of adaptation 58 from both standing variation and de novo mutations.We have 59 introduced how networks can generate epistasis, with implica-60 tions for the evolution of recombination, maintenance of link-61 age, and the trade-offs between network complexity and adapt-62 ability.Our study here shows that networks can create complex 63 distributions of fitness effects that differ from additive expec-64 tations.While promising, it is essential to acknowledge the 65 model's limitations, particularly concerning its computational 66 demands and simplifications, which future work should aim to 67 address.Empirical validation through evolutionary experiments 68 will be crucial in assessing our model's real-world applicabil-69 ity.Expanding into adaptation by standing variation will fur-70 ther illuminate the role of genetic interactions during polygenic 71 adaptation.As detailed systems models of complex traits be-72 come increasingly available, our framework can be employed to 73 simulate the evolution of traits based on empirically-described 74 networks.The utility of our model also extends beyond evolu-75 tionary biology, offering valuable insights for interdisciplinary 76 collaborations that aim to integrate molecular biology, compu-77 tational science, and statistical genetics.Future work should fo-78 cus on leveraging this interdisciplinary potential and on explor-79 ing more complex genetic networks to better understand how 80 gene interactions constrain evolution.Supplementary Table 2 Mean generalized Pareto distribution (GPD) parameters fit to mutant screen distributions of fitness effects over adaptive walks in additive and network populations.Brackets indicate 95% confidence intervals.Parameters were estimated by fitting a GPD to a random sample of mutations from the mutant screen experiments conducted on each replicate at each adaptive step.κ is the shape parameter of the GPD.κ is the mean κ value across n replicate simulations.The log-likelihood ratio tests the null hypothesis that a sample of beneficial mutations belongs to an exponential distribution (which meets κ = 0).P-values across the n replicates were combined using Fisher's method.

Figure 1
Figure 1 Examples of how the shape parameter, κ, can change the generalized Pareto distribution (GPD).The GPD, an extreme value distribution, characterizes the behavior of extreme values: random samples from the right tail of a continuous distribution.Extreme values are used to represent beneficial mutations during an adaptive walk.Depending on κ, the GPD can belong to one of three domains: Gumbel, Weibull, or Fréchet.When κ = 0, the extreme values align with the Gumbel distribution, which is the anticipated distribution for beneficial mutations during an adaptive walk (solid line).If κ < 0, the extreme values take on a Weibull distribution, exhibiting a truncated right tail (dotted line).Decreasing κ shifts the maximum extreme value (the truncation point) towards zero.Conversely, for κ > 0, the extreme values follow the Fréchet distribution, characterized by an extended right tail (dashed line).The specific κ values used in this plot are 0, -0.75, and 1 for Gumbel, Weibull, and Fréchet, respectively.
Integrating information about the underlying systems describing trait development enables us to capture genetic interactions affecting fitness.Modeling the molecular networks underpinning trait development and expression can help us understand how gene interactions and regulatory processes shape the survival and reproductive success of organisms, and provide a mechanistic view of how variation arises in natural populations.

Figure 2 A
Figure2A negative autoregulation (NAR) motif and its expression curve.The NAR motif consists of two genes, X and Z.The expression of X activates Z, and the expression of Z inhibits further Z production (A).This results in a characteristic expression curve (B).Z production begins when X is activated (blue shaded area) and stops when it is inactivated.Gene Z product approaches a steady state concentration and then quickly falls off in the absence of X.

10
In this paper, we use a novel approach to model a quantitative 11 trait as the product of a NAR motif.Previous attempts to rec-12 oncile quantitative and population genetics with systems biol-13 ogy have used a variety of approaches.Some have focused on 14 the explicit modeling of network structures (e.g.Wagner 1996; expression patterns of the NAR motif, we first 43 translate its network diagram (Figure 2A) to a system of ordi-44 nary differential equations (ODE).ODEs are commonly used 45 to model gene networks due to their balance between efficiency 46 and realism (Schlitt and Brazma 2007; Karlebach and Shamir 47 2008).The solution to the ODE predicts gene expression over

48 a
time period, such as during cell development.The NAR ODE 49 is given by: 50

5 2 .
Repeat step 1 for β Z .6 3. Substitute the α Z and β Z values into the ODE, Eq. 1, 7 and solve the ODE between time points 0 and 10 to get 8 an expression curve. 9

18 ponent
is calculated by exponentiation of the sum of allelic ef-19 fects across all mQTLs and chromosomes.For molecular com-20 ponent k: number of causal loci 23 along the genome.This treatment of allelic effects assumes no 24 explicit epistasis or dominance deviations, however these might 25 arise as consequences of the nonlinear ODE solution.26 After calculating the molecular component values (C), we can 27 determine the amount of Z expression, the phenotype.We 28 achieve this by constructing a system of ODEs (Eq. 1) and 29 solving for the area under the expression curve (i.e. the total Z 30 produced during a time period).31 We can express the ODE as an initial value problem where 32 given the "initial" or starting value of the function we can deter-33 mine the function's behavior for subsequent time points using 34 the ODE: 1.1.8(Xiang et al. 2013) in R 4.3.1 (R Core Team 2023), with code modified from Lebeuf-Taylor et al. (2019).

Figure 3
Figure 3 Phenotypic evolution over 10,000 generations of adaptation across 2,880 replicates of additive and network models.(A) The distribution of mean phenotypes among replicate populations during the adaptive period.The -1000 timepoint represents 1,000 generations before the optimum shifted.(B) The mean phenotype at each adaptive step (each fixation along the adaptive walk).The last step is pooled across all steps greater than or equal to step three in the walk.The dotted line in A and B represents the post-shift optimum.The grey area shows the overlap between additive and network distributions.

Figure 4
Figure 4 Fitness effects among fixations during adaptive walks to a phenotypic optimum across 2,880 replicates of additive and network models.(A) The overall distribution of all fixations across the entire walk.(B) The distribution of fixations at each adaptive step.The grey area shows overlap between the additive and network distributions of effects.

Figure 5
Figure 5 Fitness effects among possible mutations and the consequences for adaptation in 2,880 replicates of additive and network models.(A) The distribution of fitness effects among 1,000 randomly sampled mutations in additive and network models.The grey area represents overlap between additive and network models.(B) The proportion of mutations that are beneficial over the adaptive walk.(C) The change in beneficial mutation probability corresponds to the expected waiting time for a beneficial mutation.

Figure 6
Figure 6 The fitness and phenotype landscapes of the network model.(A) The fitness landscape with regards to the molecular components of the network, α Z and β Z .The landscape consists of a high-fitness ridge (yellow) surrounded by low-fitness valleys (purple).(B) The ratio of β Z /α Z , the steady state of the network, is what matters for adaptation.As α Z increases relative to β Z /α Z , fitness becomes constant over α Z .(C) Fitness is a nonlinear function of the network steady state, β Z /α Z , although largely linear around the phenotypic optimum.An optimum ratio at θ β Z /α Z = 0.8 maximizes fitness.The dashed black line shows the phenotypic optimum.(D) The contributions of mutations in each molecular component to the adaptive walk.This is taken for populations with at least 2 steps in the walk.φ shows the difference in sums of absolute allelic effects among fixations in α Z and β Z across an adaptive walk.Values greater than 0 represent larger changes in α Z than β Z , and vice versa for values smaller than 0. At 0, both α Z and β Z changed by the same amount during the adaptive walk.

7 12 Williamson
tions of deleterious fitness effects (e.g.Kousathanas and Keight-8 ley 2013; Zeyl and DeVisser 2001), potentially attributed to 9 different classes of mutations with distinct DFEs (Eyre-Walker 10 and Keightley 2007).Mutations in gene regulatory regions can 11 have widely varying effects on fitness (e.g.Halligan et al. 2013; et al. 2014; Ohta 2002), suggesting that network 13 structure can strongly affect the fitness distribution of mutants 14 on a given molecular component (Erwin and Davidson 2009).15 The ridge-like shape of the α Z /β Z fitness landscape reflects 16 this, showing that dependencies between the molecular compo-17 nents drive the fitness of organisms in the NAR model and that 18 both the phenotypic and fitness effect of a mutation depends on 19 its genetic background (Fig. 6).20 Unlike the additive model, the fitness landscape between 21 the molecular components of the NAR revealed non-22 interchangeable per-locus effects (Fig. 6A,C).The ridge-like 23 shape of the landscape corroborates findings by Kozuch et al.

24 ( 28 (
2020), who investigated a similar NAR motif in the E. coli lexA 25 transcription factor.NAR populations also faced lower benefi-26 cial mutation rates and a slight bias favoring α Z mutations over 27 β Z , which might have been caused by the ridge-like landscape Fig. 5, 6D).The bias suggests α Z mutations should be less 29 deleterious on average than β Z mutations.Supporting this, α Z 30 mutations were slightly less likely to be strongly deleterious, 31 instead falling between the distribution's modes (Supp.Fig. 32 13).However, both mutation types remained largely deleteri-33 ous, and had similar DBFEs.A potential explanation for the 34 bias towards α Z mutations involves the optimal β Z /α Z ratio 35 of 0.8 post-shift, versus 0.4 at burn-in.To reach this optimum 36 from burn-in, populations could increase β Z or decrease α Z .37 However, the required α Z decrease exceeded the β Z increase.38 In essence, larger α Z mutations were needed to reach the new 39 optimum.Crucially, α Z and β Z mutations did not affect the 40 phenotype equally. of the deleterious DFE for NAR populations (Fig. 5A).The 49 NAR model's increased complexity through multiple molecu-50 lar components and dynamic biochemical interactions renders 51 it more susceptible to deleterious mutations.This vulnerability 52 could amplify Hill-Robertson interference, with beneficial alle-53 les overshadowed by deleterious backgrounds (Hill and Robert-54 son 1966).For instance, consider the persistence of epistasis 55 over time.Gene interactions can lead to situations where ge-56 tion/suppression due to n Z = n XZ = 8.Varying these Hill coefficients or treating them as evolvable components in future work could offer deeper insights into how more gradual responses to X activation/deactivation and Z production might modulate total Z expression (Alon 2019).Moreover, evolving parameters like K Z and K XZ could elucidate the full NAR fitness landscape.In essence, our simplified NAR model balances complexity and tractability, avoiding certain complexities while retaining meaningful constraints.
. However, other motifs

Table 1 .
Simulation parameters.Table of symbols, names, descriptions, and values for relevant parameters used in the forward-time Wright-Fisher simulations.

Table 1
Simulation parameters