# Neutral Evolution of Mutational Robustness

## I Introduction

Kimura’s contention that a majority of genotypic change in evolution is selectively neutral  has gained renewed attention with the recent analysis of evolutionary optimization methods  and the discovery of neutral networks in genotype-phenotype models for RNA secondary structure  and protein structure . It was found that collections of mutually neutral genotypes, which are connected via single mutational steps, form extended networks that permeate large regions of genotype space. Intuitively, a large degeneracy in genotype-phenotype maps, when combined with the high connectivity of (high-dimensional) genotype spaces, readily leads to such extended neutral networks. This intuition is now supported by recent theoretical results .

In in vitro evolution of ribozymes, mutations responsible for an increase in fitness are only a small minority of the total number of accepted mutations . This indicates that, even in adaptive evolution, the majority of point mutations is neutral. The fact that only a minority of loci is conserved in sequences evolved from a single ancestor similarly indicates a high degeneracy in ribozymal genotype-phenotype maps . Neutrality is also implicated in experiments where RNA sequences evolve a given structure starting from a range of different initial genotypes . More generally, neutrality in RNA and protein genotype-phenotype maps is indicated by the observation that their structures are much better conserved during evolution than their sequences .

Given the presence of neutral networks that preserve structure or function in sequence space, one asks, How does an evolving population distribute itself over a neutral network? Can we detect and analyze structural properties of neutral networks from data on biological or in vitro populations? To what extent does a population evolve toward highly connected parts of the network, resulting in sequences that are relatively insensitive to mutations? Such mutational robustness has been observed in biological RNA structures  and in simulations of the evolution of RNA secondary structure . However, an analytical understanding of the phenomenon, the underlying mechanisms, and their dependence on evolutionary parameters—such as, mutation rate, population size, selection advantage, and the topology of the neutral network—has up to now not been available.

Here we develop a dynamical model for the evolution of populations on neutral networks and show analytically that, for biologically relevant population sizes and mutation rates, a population’s distribution over a neutral network is determined solely by the network’s topology. Consequently, one can infer important structural information about neutral networks from data on evolving populations, even without specific knowledge of the evolutionary parameters. Simulations of the evolution of a population of RNA sequences, evolving on a neutral network defined with respect to secondary structure, confirm our theoretical predictions and illustrate their application to inferring network topology.

## II Modeling Neutrality

We assume that genotype space contains a neutral network of high, but equal fitness genotypes on which the majority of a population is concentrated and that the neighboring parts of genotype space consist of genotypes with markedly lower fitness. The genotype space consists of all sequences of length $`L`$ over a finite alphabet $`𝒜`$ of $`A`$ symbols. The neutral network on which the population moves can be most naturally regarded as a graph $`G`$ embedded in this genotype space. The vertex set of $`G`$ consists of all genotypes that are on the neutral network; denote its size by $`|G|`$. Two vertices are connected by an edge if and only if they differ by a single point mutation.

We will investigate the dynamics of a population evolving on this neutral network and analyze the dependence of several population statistics on the topology of the graph $`G`$. With these results, we will then show how measuring various population statistics enables one to infer $`G`$’s structural properties.

For the evolutionary process, we assume a discrete-generation selection-mutation dynamics with constant population size $`M`$. Individuals on the neutral network $`G`$ have a fitness $`\sigma `$. Individuals outside the neutral network have fitnesses that are considerably smaller than $`\sigma `$. Under the approximations we use, the exact fitness values for genotypes off $`G`$ turn out to be immaterial. Each generation, $`M`$ individuals are selected with replacement and with probability proportional to fitness and then mutated with probability $`\mu `$. These individuals form the next generation.

This dynamical system is a discrete-time version of Eigen’s molecular evolution model . Our analysis can be directly translated to the continuous-time equations for the Eigen model. The results remain essentially unchanged.

Although our analysis can be extended to more complicated mutation schemes, we will assume that only single point mutations can occur at each reproduction of an individual. With probability $`\mu `$ one of the $`L`$ symbols is chosen with uniform probability and is mutated to one of the $`A1`$ other symbols. Thus, under a mutation, a genotype $`s`$ moves with uniform probability to one of the $`L(A1)`$ neighboring points in genotype space.

### A Infinite-Population Solution

The first step is to solve for the asymptotic distribution of the population over the neutral network $`G`$ in the limit of very large population sizes.

Once the (infinite) population has come to equilibrium, there will be a constant proportion $`P`$ of the population located on the network $`G`$ and a constant average fitness $`f`$ in the population. Under selection the proportion of individuals on the neutral network increases from $`P`$ to $`\sigma P/f`$. Under mutation a proportion $`\nu `$ of these individuals remains on the network, while a proportion $`1\nu `$ falls off the neutral network to lower fitness. At the same time, a proportion $`Q`$ of individuals located outside $`G`$ mutate onto the network so that an equal proportion $`P`$ ends up on $`G`$ in the next generation. Thus, at equilibrium, we have a balance equation:

$$P=\frac{\sigma }{f}\nu P+Q.$$
(1)   
In general, the contribution of $`Q`$ to $`P`$ is negligible. As mentioned above, we assume that the fitness $`\sigma `$ of the network genotypes is substantially larger than the fitnesses of those off the neutral network and that the mutation rate is small enough so that the bulk of the population is located on the neutral network. Moreover, since their fitnesses are smaller than the average fitness $`f`$, only a fraction of the individuals off the network $`G`$ produces offspring for the next generation. Of this fraction, only a small fraction mutates onto the neutral network $`G`$. Therefore, we neglect the term $`Q`$ in Eq. (1) and obtain:

$$\frac{\sigma }{f}\nu =1.$$
(2)   
This expresses the balance between selection expanding the population on the network and deleterious mutations reducing it by moving individuals off.

Under mutation an individual located at genotype $`s`$ of $`G`$ with vertex degree $`d_s`$ (the number of neutral mutant neighbors) has probability

$$\nu _s=1\mu \left(1\frac{d_s}{(A1)L}\right)$$
(3)   
to remain on the neutral network $`G`$. If asymptotically a fraction $`P_s`$ of the population is located at genotype $`s`$, then $`\nu `$ is simply the average of $`\nu _s`$ over the asymptotic distribution on the network: $`\nu =_{sG}\nu _sP_s/P`$. As Eq. (3) shows, the average $`\nu `$ is simply related to the population neutrality $`d=_{sG}d_sP_s/P`$. Moreover, using Eq. (2) we can directly relate the population neutrality $`d`$ to the average fitness $`f`$:

$$d=L(A1)\left[1\frac{\sigma f}{\mu \sigma }\right].$$
(4)   
Despite our not specifying the details of $`G`$’s topology, nor giving the fitness values of the genotypes lying off the neutral network, one can relate the population neutrality $`d`$ of the individuals on the neutral network directly to the average fitness $`f`$ in the population. It may seem surprising that this is possible at all. Since the population consists partly of sequences off the neutral network, one expects that the average fitness is determined in part by the fitnesses of these sequences. However, under the assumption that back mutations from low-fitness genotypes off the neutral network onto $`G`$ are negligible, the fitnesses of sequences outside $`G`$ only influence the total proportion $`P`$ of individuals on the network, but not the average fitness in the population.

Equation (4) shows that the population neutrality $`d`$ can be inferred from the average fitness and other parameters—such as, mutation rate. However, as we will now show, the population neutrality $`d`$ can also be obtained independently, from knowledge of the topology of $`G`$ alone.

The asymptotic equilibrium proportions $`\{P_s\}`$ of the population at network nodes $`s`$ are the solutions of the simultaneous equations:

$$P_s=(1\mu )\frac{\sigma }{f}P_s+\frac{\mu }{L(A1)}\underset{t[s]_G}{}\frac{\sigma }{f}P_t,$$
(5)   
where $`[s]_G`$ is the set of neighbors of $`s`$ that are also on the network $`G`$. Using Eq. (4), Eq. (5) can be rewritten as:

$$dP_s=\left(𝐆\stackrel{}{P}\right)_s,$$
(6)   
where $`𝐆`$ is the adjacency matrix of the graph $`G`$:

$$𝐆_{st}=\{\begin{array}{cc}1\hfill & t[s]_G,\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$
(7)   
Since $`𝐆`$ is nonnegative and the neutral network $`G`$ is connected, the adjacency matrix is irreducible. Therefore, the theorems of Frobenius-Perron for nonnegative irreducible matrices apply. These imply that the proportions $`P_s`$ of the limit distribution on the network are given by the principal eigenvector of the graph adjacency matrix $`𝐆`$. Moreover, the population neutrality is equal to $`𝐆`$’s spectral radius $`\rho `$: $`d=\rho `$. In this way, one concludes that asymptotically the population neutrality $`d`$ is independent of evolutionary parameters ($`\mu `$, $`L`$, $`\sigma `$) and of the fitness values of the genotypes off the neutral network. It is a function only of the neutral network topology as determined by the adjacency matrix $`𝐆`$.

This fortunate circumstance allows us to consider several practical consequences. Note that knowledge of $`\mu `$, $`\sigma `$, and $`f`$ allows one to infer a dominant feature of $`G`$’s topology, namely, the spectral radius $`\rho `$ of its adjacency matrix. In in vitro evolution experiments in which biomolecules are evolved (say) to bind a particular ligand , by measuring the proportion $`\nu `$ of functional molecules that remain functional after mutation, we can now infer the spectral radius $`\rho `$ of their neutral network. In other situations, such as in the bacterial evolution experiments of Ref. , it might be more natural to measure the average fitness $`f`$ of an evolving population and then use Eq. (4) to infer the population neutrality $`d`$ of viable genotypes in sequence space.

### B Blind and Myopic Random Neutral Walks

In the foregoing we solved for the asymptotic average neutrality $`d`$ of an infinite population under selection and mutation dynamics and showed that it was uniquely determined by the topology of the neutral network $`G`$. To put this result in perspective, we now compare the population neutrality $`d`$ with the effective neutralities observed under two different kinds of random walk over $`G`$. The results illustrate informative extremes of how network topology determines the population dynamics on neutral networks.

The first kind of random walk that we consider is generally referred to as a blind ant random walk. An ant starts out on a randomly chosen node of $`G`$. Each time step it chooses one of its $`L(A1)`$ neighbors at random. If the chosen neighbor is on $`G`$, the ant steps to this node, otherwise it remains at the current node for another time step. It is easy to show that this random walk asymptotically spends equal amounts of time at all of $`G`$’s nodes . Therefore, the network neutrality $`\overline{d}`$ of the nodes visited under this type of random walk is simply given by:

$$\overline{d}=\underset{sG}{}\frac{d_s}{|G|}.$$
(8)   
It is instructive to compare this with the effective neutrality observed under another random walk, called the myopic ant. An ant again starts at a random node $`sG`$. Each time step, the ant determines the set $`[s]_G`$ of network neighbors of $`s`$ and then steps to one at random. Under this random walk, the asymptotic proportion $`P_s`$ of time spent at node $`s`$ is proportional to the node degree $`d_s`$ . It turns out that the myopic neutrality $`\widehat{d}`$ seen by this ant can be expressed in terms of the mean $`\overline{d}`$ and variance $`\mathrm{Var}(d)`$ of node degrees over $`G`$:

$$\widehat{d}=\overline{d}+\frac{\mathrm{Var}(d)}{\overline{d}}.$$
(9)   
The network and myopic neutralities, $`\overline{d}`$ and $`\widehat{d}`$, are thus directly given in terms of local statistics of the distribution of vertex degrees, while the population neutrality $`d`$ is given by $`\rho `$, the spectral radius of $`G`$’s adjacency matrix. The latter is an essentially global property of $`G`$.

## III Mutational Robustness

With these cases in mind, we now consider how different network topologies are reflected by these neutralities. In prototype models of populations evolving on neutral networks, the networks are often assumed to be or are approximated as regular graphs . If the graph $`G`$ is, in fact, regular, each node has the same degree $`d`$ and, obviously, one has $`d=\overline{d}=\widehat{d}=d`$.

In more realistic neutral networks, one expects $`G`$’s neutralities to vary over the network. When this occurs, the population neutrality is typically larger than the network neutrality: $`d=\rho >\overline{d}`$. This difference quantifies precisely the extent to which a population seeks out the most connected areas of the neutral network. Thus, a population will evolve a mutational robustness that is larger than if the population were to spread uniformly over the neutral network. Additionally, the mutational robustness tends to increase during the transient phase in which the population relaxes towards the its asymptotic distribution.

Assume, for instance, that initially the population is located entirely off the neutral network $`G`$ at lower fitness sequences. At some time, a genotype $`sG`$ is discovered by the population. To a rough approximation, one can assume that the probability of a genotype $`s`$ being discovered first is proportional to the number of neighbors, $`L(A1)d_s`$, that $`s`$ has off the neutral network. Therefore, the population neutrality $`d_0`$ when the population first enters the neutral network $`G`$ is approximately given by:

$$d_0=\overline{d}\frac{\mathrm{Var}(d)}{L(A1)\overline{d}}.$$
(10)   
Therefore, we define the excess robustness $`r`$ to be the relative increase in neutrality between the initial neutrality and (asymptotic) population neutrality:

$$r\frac{dd_0}{d_0}.$$
(11)   
For networks that are sparse, i.e. $`\overline{d}L(A1)`$, this is well approximated by $`r(d\overline{d})/\overline{d}`$. Note that, while $`r`$ is defined in terms population statistics, the preceding results have shown that this robustness is only a function of $`G`$’s topology and should thus be considered a property of the network.

## IV Finite-Population Effects

Our analysis of the population distribution on the neutral network $`G`$ assumed an infinite population. For finite populations, it is well known that sampling fluctuations converge a population and this raises a question: To what extent does the asymptotic distribution $`P_s`$ still describe the distribution over the network for small populations? As a finite population diffuses over a neutral network , one might hope that the time average of the distribution over $`G`$ is still given by $`P_s`$. Indeed, the simulation results shown below indicate that for moderately large population sizes, this seems to be the case. However, a simple argument shows that this cannot be true for arbitrarily small populations.

Assume that the population size $`M`$ is so small that the product of mutation rate and population size is much smaller than $`1`$; i.e. $`M\mu 1`$. In this limit the population will, at any point in time, be completely converged onto a single genotype $`s`$ on the neutral network $`G`$. With probability $`M\mu `$ a single mutant will be produced at each generation. This mutant is equally likely to be one of the $`L(A1)`$ neighbors of $`s`$. If this mutant is not on $`G`$, it will quickly disappear due to selection. However, if the mutant is on the neutral network, there is a probability $`1/M`$ that it will take over the population. When this happens, the population will effectively have taken a random-walk step on the network, of exactly the kind followed by the blind ant. Therefore, for $`M\mu 1`$, the population neutrality will be equal to the network neutrality: $`d=\overline{d}`$. In this regime, $`r0`$ and excess mutational robustness will not emerge through evolution.

The extent to which the initial population neutrality approaches $`d`$ is determined by the extent to which evolution on $`G`$ is dominated by sampling fluctuations. In neutral evolution, population convergence is generally only a function of the product $`M\mu `$ . Thus, as the product $`M\mu `$ ranges from values much smaller than $`1`$ to values much larger than $`1`$, we predict that the population neutrality $`d`$ shifts from the network neutrality $`\overline{d}`$ to the infinite-population neutrality, given by $`𝐆`$’s spectral radius $`\rho `$.

## V RNA Evolution on Structurally Neutral Networks

The evolution of RNA molecules in a simulated flow reactor provides an excellent arena in which to test the theoretical predictions of evolved mutational robustness. The replication rates (fitnesses) were chosen to be a function only of the secondary structures of the RNA molecules. The secondary structure of RNA is an essential aspect of its phenotype, as documented by its conservation in evolution  and the convergent in vitro evolution toward a similar secondary structure when selecting for a specific function . RNA secondary structure prediction based on free energy minimization is a standard tool in experimental biology and has been shown to be reliable, especially when the minimum free energy structure is thermodynamically well defined . RNA secondary structures were determined with the Vienna Package , which uses the free energies from . Free energies of dangling ends were set to $`0`$.

The neutral network $`G`$ on which the population evolves consists of all RNA molecules of length $`L=18`$ that fold into a particular target structure. A target structure (Fig. 1) was selected that contains sufficient variation in its neutrality to test the theory, yet is not so large as to preclude an exhaustive analysis of its neutral network topology.

Using only single point mutations per replication, purine-pyrimidine base pairs {G-C, G-U, A-U} can mutate into each other, but not into pyrimidine-purine {C-G, U-G, U-A} base pairs. The target structure contains $`6`$ base pairs which can each be taken from one or the other of these two sets. Thus, the approximately $`2\times 10^8`$ sequences that are consistent with the target’s base pairs separate into $`2^6=64`$ disjoint sets. Of these, we analyzed the set in which all base pairs were of the purine-pyrimidine type and found that it contained two neutral networks of $`51,028`$ and $`5,169`$ sequences that fold into the target structure. Simulations were performed on the largest of the two. The exhaustive enumeration of this network showed that it had a network neutrality of $`\overline{d}=12.0`$ with standard deviation $`\sqrt{Var(d)}3.4`$, a maximum neutrality of $`d_s=24`$, and a minimum of $`d_s=1`$. The spectral radius of the network’s $`51028\times 51028`$ adjacency matrix was $`\rho 15.7`$. The theory predicts that, when $`M\mu 1`$, the population neutrality should converge to this value.

The simulated flow reactor contained a population of replicating and mutating RNA sequences . The replication rate of a molecule depends on whether its calculated minimum free energy structure equals that of the target: Sequences that fold into the target structure replicate on average once per time unit, while all other sequences replicate once per $`10^4`$ time units on average. During replication the progeny of a sequence has probability $`\mu `$ of a single point mutation. Selection pressure in the flow reactor is induced by an adaptive dilution flow that keeps the total RNA population fluctuating around a constant capacity $`M`$.

Evolution was seeded from various starting sequences with either a relatively high or a relatively low neutrality. Independent of the starting point, the population neutrality converges to the predicted value, as shown in Fig. 2.

Subsequently, we tested the finite-population effects on the population’s average neutrality at several different mutation rates. Figure 3 shows the dependence of the asymptotic average population neutrality on population size $`M`$ and mutation rate $`\mu `$. As expected, the population neutrality depends only on the product $`M\mu `$ of population size and mutation rate. For small $`M\mu `$ the population neutrality increases with increasing $`M\mu `$, until $`M\mu 500`$ where it saturates at the predicted value of $`d15.7`$. Since small populations do not form a stationary distribution over the neutral net, but diffuse over it , the average population neutrality at each generation may fluctuate considerably for small populations. Theoretically, sampling fluctuations in the proportions of individuals at different nodes of the network scale inversely proportional to the square root of the population size. We therefore expect the fluctuations in population neutrality to scale as the inverse square root of the population size as well. This was indeed observed in our simulations.

Finally, the fact that $`r0.31`$ for this neutral network shows that under selection and mutation, a population will evolve a mutational robustness that is $`31`$ percent higher than if it were to spread randomly over the network.

## VI Conclusions

We have shown that, under neutral evolution, a population does not move over a neutral network in an entirely random fashion, but tends to concentrate at highly connected parts of the network, resulting in phenotypes that are relatively robust against mutations. Moreover, the average number of point mutations that leave the phenotype unaltered is given by the spectral radius of the neutral network’s adjacency matrix. Thus, our theory provides an analytical foundation for the intuitive notion that evolution selects genotypes that are mutationally robust.

Perhaps surprisingly, the tendency to evolve toward highly connected parts of the network is independent of evolutionary parameters—such as, mutation rate, selection advantage, and population size (as long as $`M\mu 1`$)—and is solely determined by the network’s topology. One consequence is that one can infer properties of the neutral network’s topology from simple population statistics.

Simulations with neutral networks of RNA secondary structures confirm the theoretical results and show that even for moderate population sizes, the population neutrality converges to the infinite-population prediction. Typical sizes of in vitro populations are such that the data obtained from experiments are expected to accord with the infinite-population results derived here. It seems possible then to devise in vitro experiments that, using the results outlined above, would allow one to obtain information about the topological structure of neutral networks of biomolecules with similar functionality.

We will present the extension of our theory to cases with multiple-mutation events per reproduction elsewhere. We will also report on analytical results for a variety of network topologies that we have studied.

Finally, here we focused only on the asymptotic distribution of the population on the neutral network. But how did the population attain this equilibrium? The transient relaxation dynamics, such as that shown in Fig. 2, can be analyzed in terms of the nonprincipal eigenvectors and eigenvalues of the adjacency matrix $`𝐆`$. Since the topology of a graph is almost entirely determined by the eigensystem of its adjacency matrix, one should in principle be able to infer the complete structure of the neutral network from accurate measurements of the transient population dynamics.

Acknowledgments We thank the participants of the Santa Fe Institute workshop on Evolutionary Dynamics for stimulating this work, which was supported in part at SFI by NSF under grant IRI-9705830, Sandia National Laboratory, and the Keck Foundation. M.H. gratefully acknowledges support from a fellowship of the Royal Netherlands Academy of Arts and Sciences.