Fundamental Properties of the Evolution of Mutational Robustness

# Fundamental Properties of the Evolution of Mutational Robustness

Lee Altenberg 111The Konrad Lorenz Institute for Evolution and Cognition Research, Martinstrasse 12, Klosterneuburg, A3400 Austria, Lee.Altenberg@kli.ac.at
###### Abstract

Evolution on neutral networks of genotypes has been found in models to concentrate on genotypes with high mutational robustness, to a degree determined by the topology of the network. Here analysis is generalized beyond neutral networks to arbitrary selection and parent-offspring transmission. In this larger realm, geometric features determine mutational robustness: the alignment of fitness with the orthogonalized eigenvectors of the mutation matrix weighted by their eigenvalues. “House of cards” mutation is found to preclude the evolution of mutational robustness. Genetic load is shown to increase with increasing mutation in arbitrary single and multiple locus fitness landscapes. The rate of decrease in population fitness can never grow as mutation rates get higher, showing that “error catastrophes” for genotype frequencies never cause precipitous losses of population fitness. The “inclusive inheritance” approach taken here naturally extends these results to a new concept of dispersal robustness.

Keywords: genetic load — spectral gap — lethal mutagenesis — epigenetic mutation — dispersal load

Based on theoretical considerations, Kimura [1] predicted that the majority of evolutionary changes in the genome in mammals should consist of neutral mutations. Since this time, research has been focused on understanding the extent and nature of neutral genetic variation in organisms. One approach is to attempt to derive them from molecular first principles [2, 3]. The idea that neutrality may not merely be a derived consequence of molecular interactions, but actually an evolved property shaped by evolutionary dynamics, had its first manifestation in Fisher’s theory for the evolution of dominance [4], followed by Waddington [5] who proposed it as a consequence of stabilizing selection, and Conrad (cf. [6, 7]) who proposed it could result from higher-order epistatic mutations which smooth the adaptive landscape at other sites and consequently enhance evolvability. Another mechanism proposed is that natural selection for adaptations robust to environmental variation entails robustness to mutation as a generic side-effect [8].

An altogether different mechanism for the evolution of mutational robustness was proposed by Bornberg-Bauer and Chan [9] and van Nimwegen et al. [10, 11], which is that evolution along neutral networks—sets of mutationally connected genotypes with equivalent fitnesses—would “concentrate at highly connected parts of the network, resulting in phenotypes that are relatively robust against mutations”[11], and “Independent of functional fitness, topology per se can lead to concentration of evolutionary population at some sequences.”[9].

Numerous citations of these papers repeat the finding that “populations will evolve toward highly connected regions of the genome space” (e.g. [12], [13]). However, neither the original papers nor subsequent studies (to my knowledge) provide analytic proofs of this observation. With few exceptions [14] progress has been limited in identifying exactly which properties of neutral networks determine their evolved mutational robustness. As is shown by the following example, more is going on than simply a “tendency to evolve toward highly connected parts of the network”.

Figure 1 compares the equilibrium population distribution for two neutral networks that have 63 equally fit genotypes, and one inviable genotype off the network. The only difference between the networks is their mutational topology. One network is the set of all 64 nucleotide triplets, while the other is the set of copy-number variants, from to copies. Two properties of the population equilibrium distributions on the neutral networks stand out: first, the equilibrium frequency of a genotype is determined not solely by how many neutral neighbors it has, but also by its position within the network. In the copy-number network, 62 of the genotypes are identical in having only neutral neighbors. Yet the stationary distribution increases 20-fold over these 62 genotypes as they get more mutationally distant from the lethal genotype [15]. The same phenomenon is seen in the trinucleotide network but to a lesser degree. Secondly, dramatic differences are seen between the two neutral networks in the size of the genetic loads they maintain. The genetic load of the trinucleotide network is 44 times the genetic load of the copy-number network.

The differences between the evolutionary outcomes on these two mutational graphs can only be the result of their different topologies, but what properties of their topologies? Numerous results from the field of spectral graph theory are applicable to this question; here no review is attempted. While graph theory has been widely used in models of mutation, the results typically impose the assumptions that mutations occur only once during replication, mutation is symmetric and occurs at a single rate, and fitness is either zero or a single other value (all assumptions in [9, 11]).

An alternative approach is taken here, which has proven valuable in previous work [16, 17, 18, 19, 20]. The specific problem—evolution on neutral networks in nucleotide sequence space—is embedded into the larger space of problems, which includes arbitrary mutation patterns and arbitrary selection values. The generality forces one to seek the appropriate fundamental mathematical properties.

It also expands the applicability of any results. Mutation is treated generally enough to apply to non-genetic information transmission, the principal example being organismal location, where the analog of mutation is dispersal. The results here thus automatically apply to dispersal load, and in the process define a new concept of dispersal robustness. Moreover, the generality of the treatment has the potential to apply widely to diverse mechanisms of “inclusive inheritance” [21].

With arbitrary fitnesses, one can no longer characterize mutation neighborhoods simply by the fraction of mutations that are neutral, since now the distribution of fitness effects of mutation (DFE) includes advantageous or deleterious mutations, as well as neutral and lethal. The more general statistic is the expected fitness of offspring. Averaged over the population, one obtains an aggregate population mutational robustness — the degree to which the population maintains its fitness in the face of mutation pressure. A complete absence of mutational robustness occurs when the genetic load is the mutation rate, corresponding to the Haldane-Muller principle [22, 23]. Complete mutational robustness, on the other hand, would mean that the population suffers no genetic load as the mutation rate increases. A given adaptive landscape will fall somewhere between these two extrema.

The main results found here are that the population mutational robustness at a mutation-selection balance is determined by abstract spectral properties: the alignment between the fitnesses and the eigenvectors of the mutation matrix, weighted by their eigenvalues.

This spectral analysis provide a lens through which to examine models of mutation, dispersal, an inheritance generally. Models of mutation and selection are usually constructed without appreciation of how the assumptions manifest in the eigenvalues and eigenvectors of the mutation matrix (or variation operator in the case of continuous variation). We will see, for example, that the widely encountered [24] “house-of-cards” mutation model [25] is incapable of supporting mutational robustness. The results here provide direction for analyzing a wide variety of models for their capacity to support the evolution of mutational robustness, and enable comparisons to be made between different kinds of mutation processes, including nucleotide base mutation, epigenetic mark mutation, and gene copy number variation. The comparisons extend to the spectral properties of dispersal matrices, thereby unifying the results for genetic robustness and genetic load with those for “dispersal robustness” and dispersal load.

## 1 The Setting

The object of interest is the state of the population when it has converged in frequencies to an equilibrium under the following asexual, haploid evolutionary dynamics, where the population is large enough to be treated as infinite and has discrete non-overlapping generations (semelparity). The only event that changes the genotypes during reproduction is mutation; there is no recombination.

 zi(t+1) =1¯¯¯¯w(t) n∑j=1Mij(μ)wjzj(t), (1)

where

is the number of possible haploid genotypes,

is the frequency of haploid genotype in the population at time , , ,

is the fitness of genotype , ,

is the probability, when the mutation rate is , that parent of genotype has offspring of genotype , , , and

is the mean fitness of the population at time , its Malthusian rate of increase in size.

In the case of a multilocus genotype where mutation occurs independently at each of loci, the mutation probabilities may be decomposed as products of the mutation probabilities at each locus, , where index the alleles at locus .

Mutation at each locus may be parameterized by a mutation rate and a mutation distribution, , given that mutation occurs: , where and when . In matrix form this is

 M(ξ)=(1−μ)I(ξ)+μP(ξ), (2)

where each is an column-stochastic matrix, is the number of possible alleles at locus , and is the identity matrix.

The multilocus mutation matrix can be represented using the Kronecker product, , as

 M(μ)=L⨂ξ=1[(1−μ)I(ξ)+μP(ξ)]. (3)

Matrices of the form represent the situation where only a single transforming event occurs during reproduction and will be referred to as single-event mutation matrices, while those of the form (3) represent the situation where multiple independent transforming events occur during reproduction and will be referred to as multiple-event mutation matrices. In the limit of small , multiple-event matrices converge to single-event matrices, even when they model multiple loci.

From (1), equilibrium frequencies must satisfy , or in vector form,

 M(μ)D^z =(w⊤^z)^z=ˆ¯¯¯¯w ^z, (4)

where

is the vector of ones, its transpose,

is the diagonal matrix of fitness coefficients,

is the vector of fitness coefficients, and

is the mean fitness of the population.

Equation (4) shows that the mean fitness is an eigenvalue of .

It is further assumed that is irreducible when , meaning any genotype can mutate in some number of steps to any other genotype . is irreducible as well if no genotype is lethal, so from Perron-Frobenius theory we know there is only one possible equilibrium, the strictly positive Perron vector , and eigenvalue is the Perron root and spectral radius of , referred to as .

If some mutations are lethal, then is reducible, and more interesting equilibrium behavior becomes possible. The Frobenius normal form of a reducible matrix ([26, pp. 74-77]) partitions the set of genotypes into blocks, such that the restriction of to any block produces an irreducible matrix . The blocks constitute quasispecies [27]. The block with the largest spectral radius among the genotypes present at equilibrium sets the mean fitness of the population. The analysis will therefore focus on such irreducible block matrices.

### 1.1 Diploidy, frequency-dependent selection, recombination, assortative mating, and dispersal

The mutational robustness of population at equilibrium can be analyzed without knowing whether this equilibrium is dynamically stable. In the case of haploid selection on a quasispecies without recombination, the system is essentially linear and the polymorphic equilibrium is always stable. Relaxation of these assumptions to accommodate greater biological variety may alter the stability and existence conditions of equilibria but does not change the robustness analysis. The results here therefore apply to equilibria in systems with diploid and frequency-dependent selection, assortative mating, and also recombination. In this case, letting be the probability that parents and produce offspring , one gets a a frequency-dependent version of (1):

 zi(t+1) =1¯¯¯¯w(t)n∑j=1zj(t)n∑k=1TijkwjkXjk(z), (5) =1¯¯¯¯w(t)n∑j=1Mij(z)wj(z)zj(t), (6)

where is the fitness of diploid genotype , is the probability that the mate of genotype is genotype (frequency dependent), , and Under random mating, . Selfing cannot be accommodated within the product structure in (5). In models of dispersal, typically represents the density of organisms, without normalization by . Carrying capacities then makes the growth rates density dependent.

Because frequency dependence is relevant only to stability, which is not considered, the analysis to follow applies to both (1) and (6). However, frequency dependence of both and makes them composite formal quantities, rather than biological essential quantities as in the case of mutation and haploid selection.

### 1.2 Mean Fitness, Genetic Load, and Mutational Robustness

Let us retrace the calculation of population mutational robustness in the neutral network model of van Nimwegen et al. [11]. Genotypes have only two possible fitnesses, , and . The mutational robustness of a genotype with fitness is the probability that its offspring have fitness . The set of genotypes in the neutral quasispecies is . Then the mutational robustness of a genotype—it’s neutrality— is . The population neutrality at equilibrium is the expectation of censused before reproduction, with parent frequencies for . Using (4), since or for ,

 ^zi =n∑j=1Mijwj^zjˆ¯¯¯¯w=∑j∈NMijw^zjˆ¯¯¯¯w, ˆ¯¯¯¯w=n∑i=1wi^zi=∑i∈Nw^zi.

From these we see that the equilibrium average neutrality is

 ^¯¯¯ν =∑j∈Nνjw^zjˆ¯w=∑i,j∈NMijw^zjˆ¯¯¯¯w=∑i∈N^zi=ˆ¯¯¯¯ww=r(MD)w. (7)

The concept of mutational robustness can be extended beyond neutral networks by simply generalizing to .

###### Definition 1.

At a mutation-selection balance, the population mutational robustness is defined to be

 Ro(M,D):=ˆ¯¯¯¯wmaxi[wi]=r(MD)maxi[wi]. (8)

The quantity is the complement of the classical genetic load, [22, 23]. We would like to know where falls within the possible range of values given the mutation rate in . The maximum possible value is , meaning no loss of fitness from mutation. The minimum possible value will be called the Haldane limit, which is for loci.

###### Result 2 (Haldane Limits on Mutational Robustness).
 Ro(M(μ),D) :=r(M(μ)D)maxi[wi]≥(1−μ)L.

See also [28, eq. (2.3)], [29, pp. 149–150, eq. (5.31)–(5.34)]. The proof for this result and those following are provided in the Supporting Information (SI).

Robustness relative to the Haldane limit can be defined.

###### Definition 3.

Define relative mutational robustness as

 RRo(M(μ),D):=r(M(μ)D)maxi[wi]−(1−μ)L1−(1−μ)L∈[0,1].

One further quantity of interest is the average fitness of mutant offspring,

###### Theorem 4.

The average mutant offspring fitness in a population at mutation-selection balance is

 ¯¯¯¯w\sf mu(M(μ),D) :=n∑j=1¯¯¯¯w\sf mujwj^zjˆ¯¯¯¯w =ˆ¯¯¯¯w[1−(1−μ)L1−(1−μ)LVar[wiˆ¯¯¯¯w]],

assuming the mutation distributions have for all loci and alleles .

This identity precisely expresses Fisher’s “deterioration of the environment” argument [30, p. 42], that at any equilibrium, the average mutant fitness loss must exactly offset the gain in fitness from the fitness variance.

## 2 Basic Results

The core finding of [11] is that when genotype neutralities vary over the neutral network, “the population neutrality is typically larger than the network neutrality. …Thus, a population will evolve a mutational robustness that is larger than if the population were to spread uniformly over the neutral network.” While numerical examples are explored, no analytical results are provided that define precisely when it holds. It is claimed in [12, eq. (12)] that when the mutation matrix is symmetric, i.e. , this result can be obtained from Perron-Frobenius theory, but no proof is provided. A proof is here provided (in SI), but it requires more than Perron-Frobenius theory, specifically Rayleigh theory.

###### Theorem 5 (Neutral Network Robustness).

The equilibrium population robustness is always greater than the average robustness over the neutral network,

 E[ν]=∑i∈N∑j∈NMij∑ni=1∑j∈NMij.

if (1) the mutation matrix is symmetric, (2) fitness off the network is , and (3) not all genotypes on the neutral network have the same mutational robustness.

When mutation is not symmetric, counterexamples can easily be constructed where the population evolves to an average robustness that is less than the average robustness over the entire neutral network. One example is where mutation is cyclic. Another example is where asymmetric mutation is symmetrizable (a reversible Markov chain), but mutation is biased toward genotypes with few neutral neighbors (in SI).

The next result adopts to quasispecies a theorem of Karlin, which has been available for over 30 years, and was independently proven in [31].

###### Theorem 6 (5.2 of Karlin [32]).

The equilibrium genetic load within any quasispecies under single-event mutation strictly increases with mutation rate when there is any variation in fitness, i.e. for each irreducible block of loci , strictly decreases with when for any , or when .

The same outcome has been shown for multiple-event mutation when mutation at each locus is reversible, by which I refer to being the transition matrix of a reversible Markov chain, which is also equivalent to being symmetrizable.

###### Theorem 7 (Corollaries 1, 3 of [33]).

The equilibrium genetic load of a quasispecies under irreducible multiple-event mutation strictly increases with mutation rate , when the quasispecies has variation in fitness within its irreducible blocks.

Theorems 6 and 7 have important implications for the theory of lethal mutagenesis. They show that the population mean fitness will keep decreasing with greater mutation over at least the range .

###### Corollary 8 (Sufficiency for Lethal Mutagenesis).

If then for some , the population will go extinct for all .

The observation in [11] that “Perhaps surprisingly, the tendency to evolve toward highly connected parts of the network is independent of evolutionary parameters—such as mutation rate,” is solely due to there being only two fitnesses in the model, , and .

###### Result 9 (Single-Event Relative Robustness Increases with Mutation Rates).

As the mutation rate increases with single-event mutation, if a quasispecies has more than one [only one] nonlethal fitness value, then the relative mutational robustness strictly increases [is constant] with mutation rate.

For multiple-event mutation, the relative mutational robustness may increase or decrease in , as found in initial numerical exploration.

## 3 Results for Reversible Markov Chain Mutation

Theorem 7 was provable due to the tractability afforded by reversible Markov chains. This same tractability carries over to the analysis of mutational robustness. For a reversible chain, can be represented ([34, p. 33], [35], [33]) as:

 M=D1/2πKΛK⊤D−1/2π

where reversibility requires , and

is the stationary distribution for , i.e. its Perron vector,

is the diagonal matrix of the square roots of ,

is the matrix of orthogonalized eigenvectors of , , i.e. for each , , if ,

is the -th column of ,

or ,

is the diagonal matrix of eigenvalues of .

Applying this representation of , we find that the spectral radius of is that of a symmetric matrix:

 r(MD) =r(D1/2πKΛK⊤D−1/2πD) =r(KΛK⊤D−1/2πDD1/2π) =r(KΛK⊤D)=r(ΛK⊤DK) =r(D1/2KΛK⊤D1/2).

This allows use of the Rayleigh-Ritz variational formula for the spectral radius ([36, pp. 172–173], [37, pp. 176–180]), to give the first main result.

###### Theorem 10 (Expression for the Spectral Radius).

When mutation has the transition matrix of a reversible Markov chain, , then

 r(MD) =maxx⊤x=1x⊤D1/2KΛK⊤D1/2x (9)

Here we see a fundamental property of the population mutational robustness:

###### Corollary 11.

Holding and , fixed, is non-decreasing in each eigenvalue of the mutation matrix, , .

### 3.1 Bounds Derived from the Mutational Eigenvectors

Weinberger [38] introduced the idea of representing fitness landscapes using the eigenvectors of the matrix representing the mutation structure, and this “Fourier expansion” approach has been elaborated by numerous further studies (cf. [39, 40]). The second main result is now presented.

###### Theorem 12.

(Lower Bound from the Fourier Representation of Fitness) Represent the fitnesses as or , where is the Fourier coefficient for eigenvector . Then

 r(MD)≥1∑nj=1wjn∑j=1λja2j.

The Fourier coefficients measure the ‘alignment’ between and each column , and is its “amplitude” [41]. From the theory of discrete nodal domains [42], we know that as increases from to , the eigenvectors exhibit more sign changes between their mutationally connected domains, which can be considered the ‘frequency’ of the eigenvector, and contribute to ruggedness in the fitness landscape. Greater ruggedness of the landscape (weight on the with larger ) corresponds to a smaller lower bound for the population mutational robustness, an effect that can appear even before there are multiple fitness peaks.

Because Theorem 12 applies to arbitrary fitnesses, it includes the mutational landscape modifiers described in [35], early envisioned in [7]. Note that the range of mutational robustness values possible from different arrangements of fitness values depends on the spread of the eigenvalues .

Prior applications of Fourier expansions of landscapes (e.g. [39], [43]) do not actually examine selection-mutation dynamics. Rather, only mutation occurs, and fitnesses do not enter into reproduction but are simply recorded in each generation to produce a time series. The main object of analysis has been the autocorrelation function of this time series when it is stationary, derived from the covariance between parent and offspring fitnesses. The next result ties the mutation-only covariance to the mutation-selection dynamics.

###### Theorem 13 (Relation to Random Mutational Walks).

For a stationary mutational random walk with reversible Markov chain transition matrix , and stationary distribution , let be the -weighted covariance between parent and offspring fitnesses. Then at a mutation-selection balance, ,

 r(MD)≥(n∑j=1wjπj)[1+1(w⊤π)2\rm Covπ[WP,WO]].

The expression states that the mean fitness of the population is increased above the -weighted average fitness of the genotype space by the covariance between the normalized parent and offspring fitnesses in the mutational random walk.

###### Corollary 14.

When the eigenvalues of are positive, then .

### 3.2 Bounds Derived from the Mutational Eigenvalues

Now we consider bounds on that derive from the eigenvalues of for any fitness landscape, to conclude the main results.

###### Theorem 15 (Upper Bound from the Spectral Gap).

Assuming , then

 n∑i=1wiπi ≤r(MD) ≤n∑i=1wiπi+λ2(M)(maxi[wi]−n∑i=1wiπi),

where .

Theorem 15 when applied to a neutral network at biological mutation rates (which give for all ) yields:

###### Corollary 16 (Spectral Gap Bound for Neutral Networks).

Let the neutral network be referred to as . Then

 ∑i∈Nπi≤r(MD)w=Ro(M,D)≤∑i∈Nπi+λ2∑i∉Nπi. (10)

Here we see the explanation for the different mutational loads in Fig. 1. For the trinucleotide network, , while for the copy-number variation network, . Eq. (10) gives lower bounds on load as . Comparison of the bounds with the actual loads give , respectively, for the two networks.

Next we shall see that the convexity of in seen for single-event mutation in Result 9 extends to multiple-event reversible mutation.

###### Theorem 17 (Convexity of the Spectral Radius in μ).

In the case where loci mutate independently at rate , each in a reversible Markov chain, then for ,

 r(M(μ)D)= maxx⊤x=1x⊤D1/2[L⨂ξ=1K(ξ)[(1−μ)I(ξ)+μΛ(ξ)]K(ξ)⊤]D1/2x

is convex in .

Theorem 17 has an important implication: “error catastrophes” for individual genotypes as the mutation rate increases never correspond to a precipitous decline in population mean fitness. The convexity of in means that no such declines are ever possible, that in contrast, steadily declines at rates that never increase as increases. These observation, noted by [35, 44, 45] based on specific fitness landscapes, are seen to hold for all possible fitness landscapes under reversible mutation.

### 3.3 House-of-Cards Mutation

The “house-of-cards” (HoC) mutation model introduced by [25] is a single-event mutation model where for all , which makes a rank-one matrix, and therefore for .

###### Result 18 (house-of-cards Mutation).

Let , where . Then

 r(M(μ)D) ≤(1−μ)maxi[wi]+μn∑i=1πiwi.

Under biological conditions the term is near zero, since it reflects the average fitness of genotypes under the stationary distribution of mutation — i.e. a genome made out of random nucleotides. In this case, , which means that house-of-cards mutation is incapable of supporting any mutational robustness. Note that HoC mutation at individual loci does not create HoC mutation for the entire genome. But also note that HoC mutation at each locus combined with multiplicative non-epistasis () decomposes into multiple HoC systems. Therefore, with multi-locus HoC mutation, multiplicative epistasis is required for population mutational robustness to rise above the Haldane limit.

## 4 Discussion

Priority here has been given to presenting the new mathematical results on these decades old problems rather than to their biological interpretations or to particular illustrations. But a few words are in order. The chief contribution is to show how the elevation of a population’s mean fitness above the Haldane limit at a mutation-selection depends on underlying spectral relationships between the mutation structure and the array of fitnesses. The results are obtained in almost complete generality, encompassing arbitrary fitnesses and for some results, arbitrary mutation structures, and in others, reversible multilocus mutation. The results directly extend to a novel concept of “dispersal robustness”.

Infinite population theory involving the spectral radius and Perron vector was originally developed in one or two locus theory, or dispersal models, where the population size could easily exceed the number of genotypes or patches, and large populations were well approximated by infinite population models. In quasispecies theory, however, results were aimed to encompass the entire genome, entailing genotype spaces orders of magnitude larger than any possible population, rendering questionable the biological relevance of and of .

Their relevance may possibly be retained in several circumstances: (1) when a small number of loci contribute to a multiplicative fitness component, (2) when transient dynamics produce a meta-stable mutation-selection balance [46], (3) when the dynamics can be approximated by a coarse-graining to give a “phenotypic quasispecies” [47], and (3) in a finite population multitype branching process. In the latter, noted in [35], the expected number of offspring produced from one parent defines a matrix , and is the trajectory of expected numbers from an initial population, . The long-term expected growth rate of the population is then .

Each of these conditions merits further investigation. With these considerations as caveats, let us review the results here.

The first main result (Theorem 15) is that the spectral gap of the mutation matrix sets an upper bound on how much above than the Haldane limit the equilibrium population fitness can be. The spectral gap measures how rapidly mixing the mutation operator is. Theorem Theorem:Main:Eigenvalue thus shows that the more rapidly mixing mutation is, the small the mutational robustness at a mutation-selection balance. In Kingman’s well-known house-of-cards model, so the spectral gap is nearly maximal, so almost no mutational robustness is possible. Theorem 15 provides an impetus to compare different kinds of mutation models, as well as dispersal models, for their spectral gaps and therefore how mutational robustness is limited.

The second main result (Theorem 12) is that the alignment of the fitnesses with the orthogonalized eigenvectors of the mutation matrix, weighted by the eigenvalues, sets a lower bound on the population mutational robustness. This bound is low when the landscape is rugged and fitnesses align with the higher-frequency eigenvectors, and high when the landscape is smooth and fitnesses align with the low-frequency eigenvectors.

Several additional results generalize results that have been found for special cases of fitness landscapes to arbitrary landscapes: Theorems 6 and 7 show that the aggregate population mutational robustness can only decline with increasing mutation rate. Theorem 17 shows that this decrease happens at an at ever lessening rate as the mutation rate increases. So “error catastrophes” for genotypes can never cause the population fitness to plummet [44]. But high enough mutation rates will inescapably drive the population to the fitness of a random genotype sequence, so for some rate less than this, lethal mutagenesis is assured. The issue of whether lethal mutagenesis can be evaded by the evolution of mutational robustness [45] is thus resolved: with high enough mutation rates, it cannot.

These applications are to be seen as only initial illustrations of the potential usefulness of the analysis presented here. Other population dynamical processes may similarly be clarified by viewing them through this lens of their eigenvalues an eigenvectors.

## Acknowledgements

Insight for this paper occurred while hearing Ludwig Geroldinger’s dissertation defense on stepping stone and island migration models [48].I thank Erik van Nimwegen, Reinhard Bürger, Joachim Hermisson, and Nick Barton for their insightful comments. I gratefully acknowledge support from The KLI Institute, Austria, and the Mathematical Biosciences Institute at Ohio State University, USA, through National Science Foundation Award #DMS 0931642.

Supporting Information

## Proofs of the Results

For convenience, the results are restated along with their proofs.

###### Result 2 (Haldane Limits on Mutational Robustness).
 Ro(M(μ)D) :=r(M(μ)D)maxi[wi]≥(1−μ)L.
###### Proof.

At equilibrium, we examine where . Then,

 ˆ¯¯¯¯w^z1 =M11w1^z1+n∑j=2M1jwj^zj

which implies

 (ˆ¯¯¯¯w−M11w1)^z1 =n∑j=2M1jwj^zj≥0 hence [16, p. 62] Ro(M(μ),D)=ˆ¯¯¯¯ww1 ≥M11,

with equality if and only if for all either , or . Further, so , with equality if and only if for all loci . ∎

###### Theorem 4.

Letting the average fitness of mutant offspring produced by parent be

 ¯¯¯¯w\sf muj=∑ni≠jwiMij∑ni≠jMij,

the average mutant offspring fitness in a population at mutation-selection balance is therefore

 ¯¯¯¯w\sf mu(M(μ),D) :=n∑j=1¯¯¯¯w\sf mujwj^zjˆ¯¯¯¯w =ˆ¯¯¯¯w[1−(1−μ)L1−(1−μ)LVar[wiˆ¯¯¯¯w]],

assuming the mutation distributions have for all loci and alleles .

###### Proof.

The average mutant offspring offspring fitness from parent , , averaged over the population, censused after selection, before reproduction, is

 ¯¯¯¯w\sf mu (M(μ),D)=n∑j=1¯¯¯¯w\sf mujwj^zjˆ¯¯¯¯w=n∑j=1∑ni≠jwiMij∑ni≠jMijwj^zjˆ¯¯¯¯w =n∑j=1∑ni=1wiMij−wjMjj1−Mjjwj^zjˆ¯¯¯¯w.

Recall that the multiple-event mutation rates are

 Mjj=L∏ξ=1M(ξ)jξjξ=L∏ξ=1[(1−μ)+μP(ξ)jξjξ].

Assuming for all , , then . Substitution gives

 ¯¯¯¯w\sf mu(M(μ),D) =n∑j=1n∑i=1wiMij−wj(1−μ)L1−(1−μ)Lwj^zjˆ¯¯¯¯w =n∑i,j=1wiMijwj^zjˆ¯¯¯¯w−n∑j=1w2j^zjˆ¯¯¯¯w(1−μ)L1−(1−μ)L.

Since

 n∑j=1w2j^zjˆ¯¯¯¯w=ˆ¯¯¯¯wn∑j=1w2j^zjˆ¯¯¯¯w2=ˆ¯¯¯¯w(1+Var[wjˆ¯¯¯¯w]),
 ¯¯¯¯w\sf mu(M(μ),D) =ˆ¯¯¯¯w−ˆ¯¯¯¯w[(1−μ)L1−(1−μ)L]Var[wjˆ¯¯¯¯w] =ˆ¯¯¯¯w[1−(1−μ)L1−(1−μ)LVar[wiˆ¯¯¯¯w]].\qed
###### Theorem 5 (Neutral Network Robustness).

The equilibrium population robustness is always greater than the average robustness over the neutral network,

 E[ν]=∑i∈N∑j∈NMij∑ni=1∑j∈NMij.

if

1. the mutation matrix is symmetric,

2. fitness off the network is , and

3. not all genotypes on the neutral network have the same mutational robustness.

###### Proof.

Repeated use will be made of the identity for arbitrary values . We may write the network’s average robustness as

 E[ν] =∑i∈N∑j∈NMij∑ni=1∑j∈NMij=1w∑ni=1wi∑nj=1Mijwj∑ni=1∑nj=1Mijwj =1we⊤DMDee⊤De.

We have since .

Since , , so , hence

 E[ν] =1we⊤DMDee⊤De=1w2(e⊤D)DMD(De)(e⊤D)(De) ≤1w2maxx≠0x⊤DMDxx⊤x (11) =1w2 r(DMD)=r(MD)w=Ro(M,D).

Equality holds if and only if is a left eigenvector of [49, Theorem 4.2.2 (Rayleigh)], in which case

 λe⊤D=(e⊤D)DMD=we⊤DMD,

and so is also an eigenvector of with eigenvalue . This is equivalent to the condition that all genotypes on the network have the same neutrality, i.e. for all , , or

 ce⊤D=1w2e⊤DMD

with . Therefore, if there are any differences in the neutralities on the network, the inequality (11) is strict. ∎

###### Theorem 7 (Corollaries 1, 3 of [33]).

The equilibrium genetic load of a quasispecies under irreducible multiple-event mutation strictly increases with mutation rate , when the quasispecies has variation in fitness within its irreducible blocks.

###### Proof.

Corollaries 1, 3 of [33] show that

 r(M(μ)D)=r(L⨂ξ=1[(1−μ)I(ξ)+μP(ξ)]D) (12)

strictly decreases in when for any , and when there is only one equilibrium for each . In the case where there are lethal genotypes, there may be more than one equilibrium, corresponding to multiple quasispecies. In that case, one must focus on the irreducible blocks, , from the Frobenius normal form of the reducible matrix . The irreducible restriction has a unique equilibrium for each , and Corollaries 1 and 3 apply, meaning that if there is any variation in fitness among , then decreases strictly in . The quasispecies may comprise more than one block, but its mean fitness, , will therefore also be decreasing in . ∎

###### Result 9 (Single-Event Relative Robustness Increases with Mutation Rates).

As the mutation rate increases with single-event mutation, if a quasispecies has more than one [only one] nonlethal fitness value, then the relative mutational robustness strictly increases [is constant] with mutation rate.

###### Proof.

We examine the restriction of to an irreducible block . In the neutral network case where for all , then

 r (M(μ)[κ]D[κ])=wr([(1−μ)I[κ]+μP[κ])] =w[(1−μ)+μr(P[κ])]=w[1+μ(r(P[κ])−1)],

so with single-event mutation, , and

 RRo(M(μ)[κ],D[κ]) =r(M(μ)[κ]D[κ])/w−1+μ1−1+μ =1+μ(r(P[κ])−1)−1+μμ =r(P[κ])<1.

which shows the relative robustness is a constant for all .

For brevity write and

 g(μ) :=RRo(M(μ)[κ],D[κ])=f(μ)−1+μμ =1+f(μ)−1μ.

Then

 g(μ1)−g(μ2) =f(μ1)−1μ1−f(μ2)−1μ2.

Let and , .

 g(μ)−g(hμ) =hf(μ)−h−(f(hμ)−1)μh. (13)

In the case where there are more than one nonlethal fitness in the irreducible block, is nonscalar ( for any ). In [19] it is shown, by applying “dual convexity” to theorems Cohen [50] and Friedland [51], that is strictly convex in when is nonscalar. Convexity means that for all , and ,

 (1−p)f(α1)+pf(α2)>f((1−p)α1+pα2). (14)

Let and . We know that

 f(0) =r(D[κ])maxi=1,…,nwi=maxi∈κwimaxi=1,…,nwi≤1.

So 14 becomes

 <