Vol 463 | 21 January 2010 | doi:10.1038/nature08694

LETTERS
Mutational robustness can facilitate adaptation
Jeremy A. Draghi1, Todd L. Parsons1, Günter P. Wagner3 & Joshua B. Plotkin1,2

Robustness seems to be the opposite of evolvability. If phenotypes
are robust against mutation, we might expect that a population
will have difficulty adapting to an environmental change, as
several studies have suggested1–4. However, other studies contend
that robust organisms are more adaptable5–8. A quantitative
understanding of the relationship between robustness and evolvability will help resolve these conflicting reports and will clarify
outstanding problems in molecular and experimental evolution,
evolutionary developmental biology and protein engineering.
Here we demonstrate, using a general population genetics model,that mutational robustness can either impede or facilitate adaptation, depending on the population size, the mutation rate and thestructure of the fitness landscape. In particular, neutral diversity
in a robust population can accelerate adaptation as long as thenumber of phenotypes accessible to an individual by mutation issmaller than the total number of phenotypes in the fitness landscape. These results provide a quantitative resolution to a significant ambiguity in evolutionary theory.
The relationship between robustness and evolvability is complex
because robust populations harbour a large diversity of neutral genotypes
that may be important in adaptation9–11. Although neutral mutations do
not change an organism’s phenotype, they may nevertheless have epistatic consequences for the phenotypic effects of subsequent mutations12–18. In particular, a neutral mutation can alter an individual’s‘phenotypic neighbourhood’, that is, the set of distinct phenotypes thatthe individual can access through a further mutation. Pioneering studies
based on RNA folding and network dynamics suggest that genotypes
expressing a particular phenotype are often linked by neutral mutations
into a large neutral network, and that members of a neutral network differ
widely in their phenotypic neighbourhoods1,19–21. Numerous studies
have documented the importance of neutral variation in allowing a
population to access adaptive phenotypes5,17,18,22–24, and neutral networks
have consequently been proposed to facilitate adaptation9–11.
Here we analyse the relationship between robustness and evolvability
using a population genetics model that specifies statistical properties ofthe fitness landscape. Our approach bypasses the tremendous complexity of explicit neutral networks1,11,17,19,21,22 to focus instead on theessential evolutionary consequences of epistatic mutations. We consider a population of N individuals reproducing according to thediscrete-time, infinite-sites Moran model. In each time step, a randomly
chosen individual produces one offspring. Upon replication, a mutation occurs with probability m, producing a novel genotype. With
probability q, the mutation is neutral. The parameter q therefore quantifies robustness, which is assumed to be the same for all genotypes on a
network. With probability 1 2 q, the mutation is non-neutral andchanges the offspring’s phenotype to one of K phenotypes accessible
from a given genotype. Each genotype has a specific set of K accessible
phenotypes that constitute its phenotypic neighbourhood; these K
phenotypes are drawn uniformly from P possible alternatives. Phenotypic neighbourhoods are assumed to be independent, such that the Kaccessible phenotypes are redrawn whenever a mutation occurs (werelax this and other assumptions below). When the number of phenotypes accessible to an individual, K, is significantly smaller the total
number of alternative phenotypes in the landscape, P, neutral mutations can profoundly alter an individual’s phenotypic neighbourhood.
This genotype–phenotype map is illustrated in Fig. 1.
Our model implicitly represents a space of adjacent neutral networks.
Neutral mutations produce other genotypes on the focal network,whereas non-neutral mutations produce genotypes on adjacent networks, each expressing one of P alternative phenotypes. To study evolution on the focal network, we assume that initially all of the P alternative
phenotypes are lethal (our results hold more generally; see Supplementary Information, section 5). We analyse the relationship between
robustness, q, and the time required to adapt to a novel environment;
this analysis is outlined in Box 1 and detailed in Supplementary
Information, section 1.
We find that a robust population may adapt either more slowly or
more quickly than one that is less robust (Fig. 2). Starting from asteady-state population with robustness q, we consider an environmental shift that assigns one of the P alternative phenotypes thegreatest fitness. We have derived an analytic expression for the mean
waiting time before this fittest phenotype subsequently arises in the
population (Supplementary Information, section 1.4). When all
phenotypes are accessible from any genotype (K 5 P), neutral mutations have no epistatic consequences and we observe what is naively
expected: more robust populations always adapt more slowly (Fig. 2).
However, when the phenotypic neighbourhood size, K, is smaller
than the total number of phenotypes, P, we find an unexpected
pattern: the relationship between robustness and evolvability is
non-monotonic. In particular, populations with an intermediate

Figure 1 | The genotype–phenotype model. Schematic representation of the
genotype–phenotype map used in our analysis. Each circle corresponds to a
genotype; colours denote phenotypes. The model parameter q quantifies
robustness: a proportion q of mutations are neutral (solid lines) and theremaining mutations are non-neutral (dashed lines). A non-neutral
mutation changes an individual’s phenotype to one of the K accessible
alternatives that form the individual’s phenotypic neighbourhood. When K
is smaller than the total number of alternative phenotypes in the landscape,P, individuals may have different phenotypic neighbourhoods. The central
pair of adjacent genotypes shown here express the same phenotype, but they
have different phenotypic neighbourhoods.

1
Department of Biology, 2Program in Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA. 3Department of Ecology
and Evolutionary Biology, Yale University, New Haven, Connecticut 06520, USA.

353
©2010 Macmillan Publishers Limited. All rights reserved

LETTERS

NATURE | Vol 463 | 21 January 2010

Box 1 | Analysis of adaptation time
We study the time, following an environmental change, until the newly
beneficial phenotype arises in a population with robustness q. A
genotype is said to be ‘adaptable’ if its phenotypic neighbourhood
contains the beneficial phenotype; our analysis links the stochastic
evolution of these adaptable types to the adaptation time. Let p(t, y)denote the probability density of there being p
y adaptable
individuals at
ﬃﬃﬃﬃ
time t, scaling space and time by the factor N. Then p(t, y) is wellapproximated by the solution to


Lp
L2
L bKq
b(1{q)~ 2 ½yp(t, y){p(t, y) {yp(t, y)
Lt Ly
Ly P
K
where b 5 Nm. The first term in this expression quantifies genetic drift, thesecond term quantifies the increase in adaptable individuals throughmutation and the third term describes the rate of mutations that produce
the beneficial phenotype. The conflicting effects of robustness on
adaptation are evident in this expression: an increase in robustness (q)increases the supply of adaptable individuals, but it also reduces the rate
at which beneficial mutations arise in such individuals. Solving aboundary-value problem related to this equation produces an analytic
expression for the expected arrival time of the beneficial phenotype(Supplementary Information, section 1), a graph of which is shown in Fig. 2.

10
2,000
1,000

K = 5, P = 100

500

100

K = 30, P = 100

50
K = P = 100
10
0.0

0.2

0.4
0.6
Robustness, q

0.8

1.0

Figure 2 | Robustness and adaptation time. The relationship between
robustness, q, and the average waiting time before the arrival of a specificbeneficial mutation, for three fitness landscapes. Points show the means of
10,000 replicate Monte Carlo simulations, and lines show our analytic
predictions (Box 1 and Supplementary Information, section 1). When all
possible phenotypes in the landscape are directly accessible by a mutationfrom any genotype (that is, when K 5 P), robustness always inhibits
adaptation (red curve). However, when phenotypic neighbourhoods are
small (that is, when K , P), neutral mutations have epistatic consequences
and the resulting relationship between robustness and adaptation time isnon-monotonic: adaptation is most rapid at intermediate levels ofrobustness. N 5 10,000, m 5 0.001.

Number of mutant phenotypes

Mean adaptation time (generations)amount of robustness adapt more quickly than populations withlittle or no robustness (Fig. 2).
There is a simple explanation for this counterintuitive result. In a
population with little robustness (small q), most mutations are lethal
and little genetic variation accumulates. As a result, the population
may not contain any adaptable individuals, that is, those that are asingle mutation away from the beneficial phenotype. Thus, when q issmall the population may need to wait a long time before an adaptable
individual arises, and then wait further for the adaptive phenotype toarise. However, slightly more robust populations contain a greater
diversity of neutral genotypes, each of which has an independent
chance (probability K/P) of being adaptable; thus, more robust populations may adapt more quickly.
Adaptation is most rapid when a population has an intermediatelevel of robustness. Moreover, this optimal level of robustness
increases as the ratio K/P decreases (Fig. 2). This trend confirms
the primary intuition behind our result: when phenotypic neighbourhoods are small, less robust populations contain few individuals
who are ‘prepared to adapt’. In this range (q is small and K , P),
increasing robustness results in a larger repertoire of phenotypes
accessible to the population, thereby accelerating adaptive evolution.

In addition to adaptation time, we have also studied another measureof evolvability, namely the diversity of phenotypes produced by mutations in a population in a steady state. Again, the naive expectation isthat as robustness increases, fewer non-neutral mutants are produced
each generation and, as a result, the diversity of mutant phenotypes
should decrease. However, an increase in robustness also increases theneutral genetic diversity within a population, and when K is less than P,
each additional neutral type may increase the number of phenotypes
accessible to the population through mutation. Thus, as with adaptation time, an unexpected, non-monotonic, relationship is apparent
when K , P: more robust populations can produce greater phenotypic
diversity than their less robust equivalents (Fig. 3). We have derived
an analytic expression to quantify the range of parameters for which
this relationship is non-monotonic (Supplementary Information,section 2). Our analysis shows that when K is smaller than a thresholddetermined by P, N and m, the diversity of mutant phenotypes is maximized at an intermediate level of robustness.
There is an interesting difference between adaptation times andphenotypic diversity: increasing the population size or mutation rate
makes the relationship between robustness and adaptation time more
like the naive monotonic prediction, whereas it makes the relationshipbetween robustness and phenotypic diversity less like the naive monotonic prediction (Supplementary Information, section 4). Althoughthese influences of population size and mutation rate have some
intuitive basis, they demonstrate that even qualitative predictions
about the robustness–adaptability relationship require an explicit
population genetics model.
Our analysis relies on four strong assumptions: a neutral mutation
completely redraws the phenotypic neighbourhood; the number ofphenotypes, K, in a genotype’s neighbourhood is independent of its
robustness, q; the values of K and q do not vary across the neutral
network; and alternative phenotypes are generally lethal. Relaxing
each of these assumptions does not change our qualitative results(Supplementary Information, section 5). Briefly, we relax the first
assumption by introducing a parameter, f, which is the fraction of
K neighbours that are redrawn following a neutral mutation.
Allowing correlations between the phenotypic neighbourhoods ofneutral neighbours (that is, allowing f , 1) still preserves the nonmonotonic relationship between robustness and evolvability.
Furthermore, a strong linear correlation between K and q, or variation in either quantity across the network, does not change ourresults. When q varies across the network, the population evolves

K = P = 100

86
K = 5, P = 100
4

2

0
0.0

0.2

0.4
0.6
Robustness, q

0.8

1.0

Figure 3 | Robustness and diversity. The relationship between robustness, q,
and the diversity of phenotypes produced by mutation in each generation, for
two fitness landscapes. Points show the means of 100,000 replicate simulations;arrows depict slopes calculated analytically (Supplementary Information,section 2). As these results demonstrate, an increase in robustness can increase
phenotypic diversity, but only when the level of robustness, q, is small and thenumber of phenotypes accessible from a single genotype, K, is less than the total
number of phenotypes in the landscape, P. N 5 10,000, m 5 0.001.

354©2010 Macmillan Publishers Limited. All rights reserved

LETTERS

NATURE | Vol 463 | 21 January 2010

towards greater robustness as predicted by previous studies25,26.
Nonetheless, the time required to acquire a new adaptive phenotype
is still accurately described by our analytic formula, replacing the
fixed value of q by the average q in the population. The same relationship between robustness and adaptation also holds when alternative phenotypes are moderately deleterious, as opposed to lethal.
Therefore, our conclusions are not sensitive to any of the strong
assumptions used to derive our analytical results.
Our results reveal a complex relationship between robustness andevolvability. In some situations, increasing robustness will decrease
evolvability, whereas in other situations it will accelerate adaptation.
The latter phenomenon can occur only when the number of phenotypes accessible to an individual, K, is smaller than the total numberof alternative phenotypes in the landscape, P. To assess the plausibilityof this condition, and to test the assumptions and predictions of ourabstract model using an empirical, mechanistic genotype–phenotype
map, we examined the folding and evolution of simulated RNA
molecules, using the Vienna RNA Package (version 1.6.1) to estimate
reasonable values of K, P, and f for RNA. Because these parameters
vary among genotypes in an RNA neutral network, we determined
appropriate averages of K, P, and f (Supplementary Information,section 6.1). For sequences of length 40 nucleotides, we estimated thatK < 19 and that P . 60,000, confirming that K , P for RNA.
Furthermore, we found that f < 0.3, indicating that neutral mutations
substantially alter phenotypic neighbourhoods. Finally, we evolved
RNA populations in silico with varying levels of robustness, andobserved a non-monotonic relation between evolvability and robustness, which was predicted accurately by our abstract model (Supplementary Information, section 6.2).
Recent studies have used theoretical27,28 or biological5,8 examples toargue that robustness increases evolvability. Another study hasargued that robustness can either increase or decrease evolvability,
depending upon the level at which robustness is described11.
Although that study provided important intuition, it did not quantify
the effects of robustness on adaptation in an evolving population. By
contrast, our analysis describes the population genetics connecting
these important properties. This perspective allows a quantitative
resolution to opposing informal arguments, and highlights the complex interplay of influences shaping mutational robustness29,30.
Our analysis also reveals general patterns that may guide future
experimental studies. First, the relationship between robustness andevolvability can be non-monotonic. In light of this complexity,
empirical studies must go beyond pairwise comparisons of highand low-robustness strains8, to measure evolvability over a broad
range of robustness values. Second, the population size and mutation
rate in part determine whether robustness increases or decreases
adaptation time. This insight was not apparent from informal arguments linking robustness and evolvability9–11, and has not yet been
considered in any empirical work. Finally, the parameters K, P, and fprovide a new way to quantify epistasis beyond the conventionalframework of synergistic and antagonistic interactions among
selected sites.
Even though most standing genetic variation is neutral, the epistatic consequences of neutral mutations have received little experimental study. Our results demonstrate that conditionally neutral
mutations strongly influence a population’s capacity to adapt; thisform of ‘neutral epistasis’ therefore deserves direct experimental
interrogation.
Received 4 November; accepted 19 November 2009.
1.
2.

Ancel, L. W. & Fontana, W. Plasticity, evolvability, and modularity in RNA. J. Exp.
Zool. 288, 242–283 (2000).
Carter, A. J. R., Hermisson, J. & Hansen, T. F. The role of epistatic gene interactionsin the response to selection and the evolution of evolvability. Theor. Popul. Biol. 68,
179–196 (2005).

3.

Cowperthwaite, M. C., Economo, E. P., Harcombe, W. R., Miller, E. L. & Meyers, L.
A. The ascent of the abundant: how mutational networks constrain evolution.
PLoS Comput. Biol. 4, e1000110 (2008).
4. Parter, M., Kashtan, N. & Alon, U. Facilitated variation: how evolution learns frompast environments to generalize to new environments. PLoS Comput. Biol. 4,
e1000206 (2008).
5. Bloom, J. D., Labthavikul, S. T., Otey, C. R. & Arnold, F. H. Protein stability
promotes evolvability. Proc. Natl Acad. Sci. USA 103, 5869–5874 (2006).
6. Aldana, M., Balleza, E., Kauffman, S. & Resendiz, O. Robustness and evolvability ingenetic regulatory networks. J. Theor. Biol. 245, 433–448 (2007).
7. Elena, S. F. & Sanjuan, R. The effect of genetic robustness on evolvability in digitalorganisms. BMC Evol. Biol. 8, 284 (2008).
8. McBride, R. C., Ogbunugafor, C. B. & Turner, P. E. Robustness promotes
evolvability of thermotolerance in an RNA virus. BMC Evol. Biol. 8, 231 (2008).
9. de Visser, J. A. G. M. et al. Evolution and detection of genetic robustness. Evolution
57, 1959–1972 (2003).
10. Lenski, R. E., Barrick, J. E. & Ofria, C. Balancing robustness and evolvability. PLoS
Biol. 4, e428 (2006).
11. Wagner, A. Robustness and evolvability: a paradox resolved. Proc. R. Soc. B 275,
91–100 (2008).
12. Rutherford, S. L. & Lindquist, S. Hsp90 as a capacitor for morphological evolution.
Nature 396, 336–342 (1998).
13. Bergman, A. & Siegal, M. L. Evolutionary capacitance as a general feature ofcomplex gene networks. Nature 424, 549–552 (2003).
14. Kirschner, M. & Gerhart, J. The Plausibility of Life: Resolving Darwin’s Dilemma (Yale
Univ. Press, 2005).
15. Wagner, A. Robustness and Evolvability in Living Systems (Princeton Univ. Press, 2005).
16. Meyers, L., Ancel, F. & Lachmann, M. Evolution of genetic potential. PLoS Comput.
Biol. 1, e32 (2005).
17. van Nimwegen, E. Influenza escapes immunity along neutral networks. Science
314, 1884–1886 (2006).
18. Blount, Z. D., Borland, C. Z. & Lenski, R. E. Historical contingency and the evolutionof a key innovation in an experimental population of Escherichia coli. Proc. Natl
Acad. Sci. USA 105, 7899–7906 (2008).
19. Fontana, W. & Schuster, P. Continuity in evolution: on the nature of transitions.
Science 280, 1451–1455 (1998).
20. Ciliberti, S., Martin, O. C. & Wagner, A. Innovation and robustness in complex
regulatory gene networks. Proc. Natl Acad. Sci. USA 104, 13591–13596 (2007).
21. Sumedha, Martin O. C., &. Wagner, A. New structural variation in evolutionary
searches of RNA neutral networks. Biosystems 90, 475–485 (2007).
22. Koelle, K., Cobey, S., Grenfell, B. & Pascual, M. Epochal evolution shapes thephylodynamics of interpandemic influenza a (H3N2) in humans. Science 314,
1898–1903 (2006).
23. Cambray, G. & Mazel, D. Synonymous genes explore different evolutionary
landscapes. PLoS Genet. 4, e1000256 (2008).
24. Isalan, M. et al. Evolvability and hierarchy in rewired bacterial gene networks.
Nature 452, 840–845 (2008).
25. van Nimwegen, E., Crutchfield, J. & Huynen, M. Neutral evolution of mutational
robustness. Proc. Natl Acad. Sci. USA 96, 9716–9720 (1999).
26. Forster, R., Adami, C. & Wilke, C. O. Selection for mutational robustness in finite
populations. J. Theor. Biol. 243, 181–190 (2006).
27. Daniels, B. C., Chen, Y.-J., Sethna, J. P., Gutenkunst, R. N. & Myers, C. R.
Sloppiness, robustness, and evolvability in systems biology. Curr. Opin. Biotechnol.
19, 389–395 (2008).
28. Wilds, R., Kauffman, S. A. & Glass, L. Evolution of complex dynamics. Chaos 18,
033109 (2008).
29. Wilke, C. O., Wang, J. L., Ofria, C., Lenski, R. E. & Adami, C. Evolution of digitalorganisms at high mutation rates leads to survival of the flattest. Nature 412,
331–333 (2001).
30. Krakauer, D. C. & Plotkin, J. B. Redundancy, antiredundancy, and the robustness ofgenomes. Proc. Natl Acad. Sci. USA 99, 1405–1409 (2002).

Supplementary Information is linked to the online version of the paper at
www.nature.com/nature.
Acknowledgements We thank P. Turner and members of the Plotkin laboratory foradvice and feedback. J.B.P. acknowledges funding from the Burroughs Wellcome
Fund, the David and Lucile Packard Foundation, the James S. McDonnell
Foundation, the Alfred P. Sloan Foundation, the Defense Advanced Research
Projects Agency (HR0011-05-1-0057) and the US National Institute of Allergy and
Infectious Diseases (2U54AI057168). G.P.W. acknowledges funding from theJohn Templeton Foundation and the Perinatology Research Branch of the US
National Institutes of Health.
Author Contributions J.A.D., J.B.P. and G.P.W. designed the project. J.A.D. and
J.B.P. wrote the paper; G.P.W. and T.L.P. edited the paper. J.A.D. performed thesimulations. T.L.P. performed most of the analysis, with contributions from J.B.P.
and J.A.D. J.A.D., T.L.P. and J.B.P. wrote the Supplementary Information.
Author Information Reprints and permissions information is available atwww.nature.com/reprints. The authors declare no competing financial interests.
Correspondence and requests for materials should be addressed to J.B.P.
(jplotkin@sas.upenn.edu).

355©2010 Macmillan Publishers Limited. All rights reservedMutational robustness can facilitate adaptation
Supplementary Information

doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Jeremy A. Draghi1 , Todd L. Parsons1 , Günter P. Wagner3 , Joshua B. Plotkin1,2

1

Department of Biology, University of Pennsylvania, Philadelphia, PA, 19104, USA
Program in Applied Mathematics and Computational Science, University of Pennsylvania, Philadelphia, PA, 19104, USA
3
Department of Ecology & Evolutionary Biology, Yale University, New Haven, CT,
06520, USA
2

1

Robustness and adaptation time

In this section we study how robustness is related to adaptation time.

1.1

Model Description

As in the main text, we consider a population of N individuals under the inﬁnite
alleles Moran model. In each discrete time step, a randomly-chosen individual produces one oﬀspring, which replaces a random individual. A mutation occurs withprobability µ and produces a unique genotype. With probability q a mutation isneutral; q therefore quantiﬁes robustness. Otherwise, the mutation is non-neutral
and changes the phenotype to one of K phenotypes accessible from a given genotype. Each genotype has a speciﬁc set of K accessible phenotypes which constituteits phenotypic neighborhood; these K phenotypes are drawn uniformly from P possible alternatives. Genotypes have independent phenotypic neighborhoods, so theK accessible phenotypes are redrawn whenever a mutation occurs. The model is
described by the ﬁve parameters K, P , µ, q, and N .
Starting from a population in steady state, we consider an environmental shiftthat assigns one of the P alternative phenotypes the highest ﬁtness. Before the environmental shift, we assume that all P of the alternative phenotypes are inviable,
such that only genotypes expressing the wild-type phenotype survive and reproduce.
We consider a population evolving in this regime of stabilizing selection until itreaches steady-state, as described in Section
1.4 below. Sometime after the popu1
lation reaches steady-state, the selective environment shifts such that one of the P
alternative phenotypes is no longer inviable, but is instead more ﬁt than the wildtype. In this section we derive an analytic expression for the mean adaptation time
– i.e. the average amount of time elapsed before the newly beneﬁcial phenotype ﬁrst
arises in the population.
Let t = 0 denote the time at which the alternative phenotype becomes more ﬁtthan the current phenotype – i.e. the time of the environmental shift. At this time
we classify all genotypes in the population into three distinct groups. Genotypes that
express the newly beneﬁcial phenotype are called class C. Genotypes that express thewild-time phenotype are divided into two classes: those that can reach the optimal
phenotype by a single point mutation, called class B or ’adaptable’, and those thatcannot reach the optimal phenotype by a single mutation, called class A. This
simpliﬁcation into three classes is equivalent to the inﬁnite-allele model described in
the main text because all genotypes within a class have identical mutations rates toother classes.
A mutation that arises in a genotype of class A produces a genotype of class B
with probability q K, and it produces another genotype in class A with probabilityP
q(1 − K).
The
sameholds for mutations from class B to class A. However, aP
www.nature.com/nature
mutation arising in a genotype of class B might alternatively produce a genotype ofclass C; this probability of this occurrence is given by the chance that the mutation is1doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

We consider a population evolving in this regime of stabilizing selection until itreaches steady-state, as described in Section 1.4 below. Sometime after the population reaches steady-state, the selective environment shifts such that one of the P
alternative phenotypes is no longer inviable, but is instead more ﬁt than the wildtype. In this section we derive an analytic expression for the mean adaptation time
– i.e. the average amount of time elapsed before the newly beneﬁcial phenotype ﬁrst
arises in the population.
Let t = 0 denote the time at which the alternative phenotype becomes more ﬁtthan the current phenotype – i.e. the time of the environmental shift. At this time
we classify all genotypes in the population into three distinct groups. Genotypes that
express the newly beneﬁcial phenotype are called class C. Genotypes that express thewild-time phenotype are divided into two classes: those that can reach the optimal
phenotype by a single point mutation, called class B or ’adaptable’, and those thatcannot reach the optimal phenotype by a single mutation, called class A. This
simpliﬁcation into three classes is equivalent to the inﬁnite-allele model described in
the main text because all genotypes within a class have identical mutations rates toother classes.
A mutation that arises in a genotype of class A produces a genotype of class B
with probability q K, and it produces another genotype in class A with probabilityPKq(1 − P ). The same holds for mutations from class B to class A. However, amutation arising in a genotype of class B might alternatively produce a genotype ofclass C; this probability of this occurrence is given by the chance that the mutation is
non-neutral, (1 − q), times the probability that a non-neutral mutation will produce
the distinguished beneﬁcial phenotype out of K phenotypic neighbors, K1 . Thesemutation rates are summarized in Supplementary Figure 1.
We assume that the population dynamics follow a discrete-time Moran modelwith a total population size N : at each time step, a randomly-chosen individual
produces one oﬀspring that replaces any individual, including the parent, with equal
probability. We assume that all individuals are equally likely to give birth, though
our results remain unchanged if we assume one of class A or B is weakly selected.
µ( K
q
P )
µ(1− K
q
P )

 A 

µ(

1− K
P

)q



B

µ
(1−q)
K

 C

µ( K
q
P )

Supplementary Figure 1 – Probability of mutation within and between genotype
classes2

www.nature.com/nature2doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

We are interested in the ﬁrst time to arrival of an individual of type C. Prior to
the ﬁrst arrival, the total number of individuals is held ﬁxed at N , and so it suﬃces
to keep track of the number of individuals in class B until the ﬁrst arrival. When theﬁrst individual of class C arises, we assume that the process jumps to an absorbing“graveyard” state ∆, at which time we stop the process. We denote this Markov
stopping time by τ∆ . We wish to compute the expectation of τ∆ in terms of theparameters K, P , µ, N , and q.
Let XN (t) denote the number of individuals in class B in the tth time step fort < τ∆ . For t ≥ τ∆ we let XN (t) = ∆. Then XN (t) is a discrete-time Markov chaindef
on the set I = {0, . . . , N } ∪ {∆}. The transition probabilities of this chain,QNi,j = P {XN (t + 1) = j|XN (t) = i} ,are given by 
2   

N
−
i
N
−
i
i
+µ
q
,
+
QN
i,i+1 = (1 − µ)
N
N
N


  2   

 
i
N −i
N −i
i
i
K
N
Qi,i−1 = (1 − µ)
+µ 1−
q
,
+
N
N
P
N
N
N
 
i
µ
N
Qi,∆ = (1 − q)
,
K
N
N
N
N
QN
i,i = 1 − Qi,i+1 − Qi,i−1 − Qi,∆ ,


i
N



N −i
N





K
P

N
for 0 ≤ i ≤ N . Furthermore, QN∆,∆ = 1, and Qi,j = 0 for all other pairs i, j ∈ I.
We use the notation QN to denote the matrix with entries QNi,j , the generator of theMarkov chain XN (t).

1.2

Two continuous-time approximations

The discrete process XN (t) is too complicated to consider directly, so we will introduce two diﬀerent limiting processes that asymptotically capture the essential
behavior. We will use these continuous-time approximations to derive an analytic
expression for the mean adaptation time. This approximation will be asymptoticallyaccurate for large N .
We analyze two diﬀerent time- and mass-rescaled processes:−1
Y1,N (t) = η1,N
XN (α1,N t) ,−1
XN (α2,N t)
Y2,N (t) = η2,N

3

www.nature.com/nature3doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

where we deﬁne the choices of α1,N  1, η1,N  1, α2,N  1, and η2,N  1 below.
For each n = 1, 2 note that Yn,N is a continuous-time processes, whereas XN isdiscrete. One unit of time for Yn,N (t) corresponds to αn,N time units in the discrete
process XN (t). By deﬁnition, Yn,N (t) is O(1) when XN (αn,N t) ∝ ηn,N .
1
We show below that for small initial values XN (0), by choosing η1,N = N 2 and3α1,N = N 2 the process Y1,N (t) converges in distribution to a diﬀusion process withdef
killing, Y1 (t), on R∆+ = [0, ∞) ∪ {∆}, with transition density functionp(t, z, y) dz = P{Y1 (t) ∈ [z, z + dz)|Y1 (0) = y},satisfying the forward Kolmogorov equation∂2β(1 − q)∂ βKq∂p
= 2 [zp(t, z, y)] −p(t, z, y) −zp(t, z, y).
∂t
∂z
∂z
P
K
where β = N µ. Thus, when the number of adaptable individuals is small, backmutations have negligible eﬀect and the rate of mutation to the adaptable type is
eﬀectively constant. Changes in number are a result of three eﬀects, encapsulated
in the three terms: genetic drift, mutations from non-adaptable to adaptable types
and ﬁnally mutations from adaptable types to the target phenotype. Even when
the speciﬁc form of the mutation rates is varied, as we will consider below, thisasymptotic form, and the interpretation of its three terms, remains unchanged.
We will also show below that by choosing η2,N = N and α2,N = N the process
Y2,N (t) converges to a pure jump process that jumps into state ∆ with probability
one, after a random time which is exponentially distributed with rateβ(1 − q)Y2 (0).
K
In the following sections, we give a proof of these two limits. Readers who areless interested in the formal details may skip ahead to Section 1.3 for the derivationof the expected ﬁrst arrival time.
1.2.1

Generators of Stochastic Processes

Let Y (t) be some continuous-time stochastic process, and let Py and Ey denoteprobability and expectation conditioned on Y (0) = y, respectively. We recall from
(Karlin & Taylor 1981, Ethier & Kurtz 2005) that the inﬁnitesimal generator A of
Y (t) isEy [f (Y (t))] − f (y)def
Af (y) = lim
,
t↓0
t
4

www.nature.com/nature4doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

for all f for which the limit on the right exists. We refer to all such f as the domainof A, D(A). For continuous time and continuous state processes, the generator
plays a role analogous to the transition matrix, and may be used, in conjunction
with restrictions on the domain, to uniquely characterize the process. Inﬁnitesimal
generators provide a convenient unifying framework within which to study Markov
processes.
For example, for a diﬀusion process with killing Z(t), with probability transition
density p(t, z, y), i.e.
 b
Py {Z(t) ∈ (a, b)} =p(t, z, y) dz,asatisfying the Komolgorov forward equation
∂p
1 ∂2
∂
=
[b(z)p(t, z, y)] − k(z)p(t, z, y),[a(z)p(t, z, y)] −
2
∂t
2 ∂z
∂z
the corresponding generator is1
Af (y) = a(y)f  (y) + b(y)f  (y) − k(y)f (y).
2
The domain of A consists of all twice diﬀerentiable functions satisfying appropriate
boundary conditions (absorbing, reﬂecting, etc.) and vanishing at the graveyard
point ∆.
A continuous time Markov jump process, where the time to the next jump, starting from y, is exponentially distributed with rate r(y), and µ(y, dz) is the probability
that the process jumps from y to z, has generator
Af (y) = r(y) f (z) − f (y) µ(y, dz).
For the two limits described above, we have,
A1 f (y) = yf  (y) +βKq β(1 − q)f (y) −yf (y).
P
K

(1.1)

1

corresponding to ηn,N = N 2 , for ηn,N = N , we get generatorβ(1 − q)y (f (∆) − f (y)) ,Kβ= − (1 − q)yf (y)K

A2 f (y) =

(1.2)

5

www.nature.com/nature5doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

for f continuous in [0, 1] (and vanishing at ∆). The latter generator describes a
Markov jump process with jump rate r(y) = Kβ (1 − q)y and jump distribution given
by a Dirac point mass, µ(y, dz) = δ∆ (dz), i.e. all jumps are to the graveyard state.
A2 uniquely characterizes the corresponding process. For A1 however, we also needto determine appropriate boundary conditions, which we discuss below.
1.2.2

Domain and Boundary Conditions for A1

In general, the inﬁnitesimal generator does not uniquely specify a diﬀusion process.
We must also characterize the domain to which we apply the generator by determining or imposing appropriate boundary conditions. We do so by following Feller’sboundary classiﬁcation (Feller 1954a,b) (see also (Karlin & Taylor 1981, Ethier &Kurtz 2005)).
For A1 , ∞ is always a natural boundary, and cannot be reached in ﬁnite�time. Ac
cording to Feller’s classiﬁcation scheme we consider the parameter ν = 12 1 − βqK.
PIf ν < 0, then 0 is a regular boundary, and it is an entrance boundary if ν > 0. In theformer case, the process can enter or leave at 0, and we may specify any boundary
condition from absorption to reﬂection. In the latter, the process can enter the interval (0, ∞) in ﬁnite time if started from 0, but can never reach 0 when started froman interior point. Since our original ﬁnite Markov chain model allows for mutation
away from 0, we will impose reﬂecting boundary conditions at 0 for ν < 0.
Whether the boundary behavior at 0 is entrance or reﬂecting, the corresponding
condition for a function f belonging to D(A) is the same,f  (y)= 0.
y↓0 s(y)

lim
where

s(y) = y 2ν−1 .
The natural boundary at ∞ leads to the requirement
lim f (y) = 0.

y→∞We also require that ∆ be an absorbing state, which implies thatf (∆) = 0
for all f ∈ D(A1 ). Lastly, we require that all functions in D(A1 ) be continuous on[0, ∞) and be twice continuously diﬀerentiable on (0, ∞). Together, these conditions
specify the domain of A1 and uniquely determine the process Y1 (t).
6www.nature.com/nature6doi: 10.1038/nature08694

1.2.3

SUPPLEMENTARY INFORMATION

Deriving the Limiting Generators

To prove convergence, Theorem 6.5 in Chapter 1 of (Ethier & Kurtz 2005) provesthat it suﬃces to show that, for appropriate choices of αn,N and ηn,N ,sup |αn,N (QN − I)f (y) − An f (y)| = 0

lim

N →∞ y∈Gn,N

for n = 1, 2 and f ∈ D(Ai ) ∩ C 3 (0, ∞), where−1
j|j ∈ N, 1 ≤ j ≤ N },Gn,N = {ηn,N

is the set of all possible values of Yn,N (t).
−1Now, for y = ηn,Ni ∈ Gn,N ,−1i) = αn,N
αn,N (QN − I)f (ηn,N= αn,Nj=i,∆

QN
i,j


j∈I

−1
−1
QN
i,j f (ηn,N j) − αn,N f (ηn,N i)



�
−1
−1
−1
j) − f (ηn,Ni) − αn,N QN
f (ηn,Ni,∆ f (∆) − f (ηn,N i)�−1i, and recalling that f (∆) = 0, yields
Taylor expanding f (y) about y = ηn,N

= αn,N



1 −2
1 −3
−1
 −1
2  −1
3 
η
η
(j
−i)f
(η
i)+
QN
(j
−i)
f
(η
i)+
(j
−i)
f
(ξ)
η
i,j
n,N
n,N
n,N
2 n,N
6 n,N
j=i,∆
−1− αn,N QNi,∆ f (ηn,N i)� −1  −1
� N
 −2  −1
1
N
N
η
Q
α
= αn,N QN
−
Q
f
(η
i)
+
+
Q
n,N
i,i+1
i,i−1
i,i+1
i,i−1 ηn,N f (ηn,N i)
n,N
n,N
2
�
 −3 
1−1NN+ αn,N QN−
Qi,i+1
i,i−1 ηn,N f (ξ) − αn,N Qi,∆ f (ηn,N i)
6

−1
−1
−1
for some ξ between ηn,Ni and ηn,Nj. Substituting y = ηn,Ni and β = N µ into QNi,jand simplifying yields


−1
ηn,N
1
K
= αn,N βq− 2 y f  (y)
N P
N

 
 

−1
−1
−2
−2
ηn,N
ηn,N
ηn,N
ηn,N
αn,N
K
1 2
K
+
y − 2 2 y − 2 2 βy +βq2
+ 4 βq 1 − 2
y f  (y)
2
N
N
N
N
P
N
P


−3
−1
αn,N ηn,N K ηn,N
αn,N β+− 2 y f  (ξ) − −1 2 (1 − q)yf (y)6N PNηn,N N K

7

www.nature.com/nature7doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

We look for scalings for which all terms are ﬁnite, so that the limit is wellposed, and for which the killing term is non-zero. For the former, we must have1−1
ηn,N = O(N 2 ) or smaller, while for the latter, we must have αn,N = ηn,N
N 2.
3
1
If we take α1,N = N 2 and η1,N = N 2 , thenα1,N (QN − I)f (y) → A1 f (y)as N → ∞, giving the diﬀusion limit, while if we take η2,N = N and α2,N = N , thenα2,N (QN − I)f (y) → A2 f (y)1

Scalings with N  ηN  N 2 give rise to a generator identical to A2 ,Af (y) = −β(1 − q)yf (y),Konly the domain now consists of functions continuous on [0, ∞) and vanishing at ∆.
All other choices of αN and ηN result in a generator that becomes unbounded or
tends to 0 as N → ∞.

1.3

Adaptation time from ﬁxed initial condition

In this section we calculate the expected value of τ∆ , the ﬁrst time at which the
newly beneﬁcial phenotype arises, for each of our two limiting processes Y1 (t) andY2 (t), conditioned on Y1 (0) = y or Y2 (0) = y:def

T1 (y) = EYy 1 [τ∆ ]defand T2 (y) = EYy 2 [τ∆ ] .

Once we account for the rescaling of time and mass, each of these may be used
as an asymptotic approximation to the expected ﬁrst arrival time starting from a
population with i individuals of class B, in our original discrete model:defN
T N (i) = EX
[τ∆ ] .
i1

In particular, for i ∝ N 2 , we will use Y1 (t) to approximate XN (t), so
31

T N (i) ∼ N 2 T1 (N − 2 i),(1.3)

while for i ∝ N , we will use Y2 (t) to approximate XN (t), so
T N (i) ∼ N T2 (N −1 i).

(1.4)

8

www.nature.com/nature8doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

In fact, as we show in Section 1.4, Eq.(1.3) is valid for all 1 ≤ i ≤ N .
When we consider generator A2 (Eq.(1.2)), the mean adaptation time is simply
the mean of an exponential distribution:
1
.
(1.5)
T2 (y) = β(1 − q)yK

On this timescale, the numbers of individuals in class B remains ﬁxed before C arrives. Unfortunately, this simple expression becomes unbounded as y → 0, because
individuals in class C can only arise from B individuals. Moreover, the expression above cannot be integrated against the steady-state distribution (see below).
Therefore, we must use a diﬀerent time scaling in order to determine the expectedarrival time starting from very small numbers in class B. We will use the diﬀusion
approximation with generator (1.1).
Since ∆ is the only absorbing state for our process, T1 (y) is the expected ﬁrst
time to absorption for Y1 (t) conditioned on starting from y, which can be written asa Dirichlet problem for the generator (Karlin & Taylor 1981),
A1 T1 (y) = −1
T ∈ D(A1 )i.e.
yT1 (y) +βKq β(1 − q)T1 (y) −yT1 (y) = −1PKlim y 1−2ν T1 (y) = 0y↓0

lim T1 (y) = 0.

y→∞This can be readily solved via the method of Green’s functions (Karlin & Taylor
1981).
We ﬁrst ﬁnd solutions u0 (y), u∞ (y) to the homogeneous problem A1 u(y) = 0,uwith u0 (y) and u∞ (y) satisfying the boundary conditions at zero ( s0 (0+) = 0) andinﬁnity (u∞ (∞) = 0), respectively:
u0 (y) = y ν I−ν (ay)andu∞ (y) = y ν Kν (ay).

�

βqK
1
1
−
. Iν (z) and Kν (z) are the modiﬁed
and,asbefore,ν=
where a = β(1−q)
K
2P
Bessel functions (Abramowitz & Stegun 1965). The Green’s function G(y, ξ) is thengiven bym(ξ)u0 (y)u∞ (ξ) if ξ ≤ yG(y, ξ) =m(ξ)u0 (ξ)u∞ (y) if y ≤ ξ.
9www.nature.com/nature9doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

wherem(y) = y −2ν .
Intuitively, G(y, ξ) dξ represents the expected time Y1 (t) spends in [ξ, ξ + dξ),given that Y1 (0) = y. Thus the expected adaptation time is given by
 ∞
 ∞
 y
T1 (y) =G(y, ξ) dξ = u0 (y)m(ξ)u∞ (ξ) dξ + u∞ (y)m(ξ)u0 (ξ) dξ
0
y
0
√

π  a ν � 1νΓ 2 −ν
= y I−ν (ay)2a 2 � νa� ay 2 
y
1 π
2
1 3
+
1 F2 2 ; 2 , 1 + ν; 2
2 sin(νπ) Γ(1 + ν)
� a −ν 1−2ν



�
y
2
ay
2
1
3
−
1 F2 2 − ν; 2 − ν, 1 − ν; 2
(1 − 2ν)Γ(1 − ν)
� a −ν 1−ν� ay 2 y2
1
3+ Kν (ay), (1.6)
1 F2 2 − ν; 2 − ν, 1 − ν; 2
(1 − 2ν)Γ(1 − ν)

where 1 F2 (a1 ; b1 , b2 ; z) is a generalized hypergeometric function (Slater 1966). Thisequation, while unwieldy, gives an analytic expression for the mean adaptation time,
starting from a ﬁxed frequency of B-type individuals.
Using asymptotic properties of the Bessel functions, it is possible to show that 1 
T1 (y) = β(1−q) 1 + O y1 ,
(1.7)
y
K
in agreement with (1.5), the expression we previously obtained for large initial frequencies. Thus (1.6) is valid not only for small initial frequencies of XN (0), but infact gives a uniform asymptotic estimate for all starting frequencies XN (0).
1.3.1

Approximation near y = 0

From Abramowitz & Stegun (1965), we have ν11Iν (z) ∼z ,Γ(1 + ν) 2and
1 πKν (z) =2 sin(νπ)

1Γ(1 − ν)



1
z
2

−ν

0 ≤ z  1.

1−
Γ(1 + ν)



1
z
2

ν 

(1.8)

�
� 
1 + O z1

(1.9)

10

www.nature.com/nature10doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Combining these
expressions
 ∞ asymptotics
 y −νwith hypergeometric series expressions for−ν
the integrals y ξ Kν (aξ) dξ and 0 ξ I−ν (aξ) dξ, and simplifying using the reﬂection formula for the gamma function (Abramowitz & Stegun 1965),
Γ(ν)Γ(1 − ν) =π,sin(νπ)we obtain an asymptotic expression for T1 (y) for y very small,
�



√
y
1
π Γ 12 − ν
T1 (y) =+1−+ o(y)
2a Γ(1 − ν)2ν
1 − 2ν
In particular, we obtain a simple expression for the expected ﬁrst arrival time froma population consisting only of A-type individuals,�√π Γ 12 − ν
T1 (0+) =.
2a Γ(1 − ν)We may further simplify this using Gauss’ duplication formula for the gamma function (Abramowitz & Stegun 1965),

to obtain


�
22z−1
Γ(2z) = √ Γ(z)Γ z + 12 ,π�
22−2ν−1 Γ 12 − ν
T1 (0+) =a Γ(1 − 2ν)� βqK 2βqK
2 P −2 Γ 12 P� βqK  .
=β(1−q) ΓP(1.10)

K

1.4

Adaptation time from steady state

Although (1.6) gives an analytic expression for the mean adaptation time startingfrom a ﬁxed initial number of class-B genotypes, XN (0), we are actually interested
in the adaptation time starting from a population in steady state. Therefore, wemust integrate (1.6) over the probability distribution for the frequency of class-B
genotypes in steady state.
Prior to the environmental shift, we have assumed that all genotypes not expressing the wild-type, including those in class C, are inviable. Therefore, the relative
frequencies of class A and class B genotypes follow a standard, neutral Moran model
11

www.nature.com/nature11doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

with asymmetric mutations between two types. Genotypes of� class A mutate to classK
K
B at rate βq P , and class B mutates to class A at rate βq 1 − P . Therefore thesteady state frequency of class B individuals follows the well-known beta distribution(Ewens 2004) with probability density function
P {XN (0) = x} = g(x) = Γ β 1− KΓ(βq)
xβ (1− P )q−1 (1 − x)β ( P )q−1 .
( ( P )q)Γ(β ( KP )q)K

K

Thus, the expected adaptation time, in generations, starting from a population
at steady state is11

T = N2

1

g(x)T1 (xN 2 ) dx.

(1.11)

0

The integrand above is diﬃcult to compute numerically for large arguments, x. Fortunately, for x large we have the asymptotic expansion given by Eq. 1.7. Therefore,
in practice we numerically evaluate the expression above by dividing the integral into
two regimes: 3/√N
 1
1
1
2
2
T ≈N
g(x)T1 (xN ) dx + √ g(x)T2 (x) dx
3/ N0When K  P , almost all of the mass of the Dirichlet distribution is concentrated
near x = 0; i.e. with very high probability there are few or no adaptable types attime t = 0. Thus, some degree of robustness is necessary for the arrival of suﬃcient
adaptable types to survive drift and mutate to the target phenotype. As q increases
from 0, the rate of arrival of adapted phenotypes increases, accelerating the arrivalof the target type. However, when q approaches 1, the rate of mutation to thetarget type becomes exceedingly small, creating the observed non-monotonicity inthe expected arrival time.

2

The relationship between robustness and phenotypic diversity

In this section we analyze how robustness, q, is related to the diversity of phenotypes produced by mutation, each generation, in a population at steady state. This
calculation avoids the temporal complexity of waiting times and addresses a simple,
general question: can populations that are more robust to mutation ever produce
greater phenotypic diversity?
As before, we assume that one phenotype is ﬁt and that the P alternative phenotypes are lethal. Mutations to a given genotype can produce only K phenotypes,
12

www.nature.com/nature12doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

q is the proportion of mutations that are neutral, and a neutral mutation produces anovel genotype with a new set of K mutational neighbors. Diversity is measured as
the expected number of distinct lethal phenotypes produced in a single generation,
which we refer to as D(q).
We wish to understand under what parameters a population with some levelof robustness can produce more phenotypic diversity than a comparable population
with no robustness in the genotype-phenotype
map. In other words, we wish todD(q) understand the derivative dq  . In this section we derive an analytic expression
q=0

for this derivative in terms of the model parameters K, P , µ, and N .
As before we consider a haploid population of N asexuals which mutate at a
rate µ and reproduce according to the Moran process. Let A denote the number ofdistinct neutral genotypes present in the population. We deﬁne
βq = N µ(1 − q)andβn = N µq

1 − µN +1 q N +1.
1 − µqβq denotes the expected number of non-neutral mutations that arise in the population, each generations. βn is the neutral mutation rate appropriate to the inﬁnitealleles Moran model (Ewens 2004), which we use to determine the limiting number
and distribution of neutral alleles.
Let αA,m denote the expected number of novel phenotypes among m mutations
arising in a population with A neutral types, given that each mutation gives rise to
one of K equally-probable neighboring phenotypes. Then, the number of non-neutral
mutations is Poisson distributed with mean βq , andD(q) =∞∞
βqm e−βq 

m=0

m!

αk,m P{A = k}.

(2.1)

k=0

If q = 0, there there is only one neutral genotype present, i.e. A = 1 withprobability one. This genotype can produce up to K alternative phenotypes by
mutation. The expected number of non-neutral mutations in a single generation isthen β0 = N µ. We thus haveD(0) =∞
β m e−β0
0

m=0

m!

α1,m

(2.2)

13

www.nature.com/nature13doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

We now turn to computing α1,m . Let ξi = 1 if the ith mutation gives rise to a novel
phenotype, and 0 otherwise. Then,
ξ=

m


ξi

i=1

is the number of distinct phenotypes, and
i−1

1
.
P{ξi = 1} = 1 −
K
Thus,



m




1
E [ξi ] = K 1 − 1 −α1,m = E [ξ] =
K
i=1

m 

We
 now turn to D(q) for q > 0. Since we are only interested in the derivative, we need only determine D(q) up to o(q). To this end, we observe thatdD(q) dq q=0

βn = N µq

 N

µn q nn=0

= β0 q + O(q 2 ).

Next, we consider P{A = k}. The distribution of neutral genotypes in the inﬁnitealleles Moran model may be obtained via Hoppe’s urn (Ewens 2004). Brieﬂy, we starti−1with a population of size 1. At each time step, with probability i−1+βa randomnβnindividual produces a clone, while with probability i−1+βn , we add an individual witha novel genotype. After N iterations, we have a population with N individuals, for
which
k



βn
i−1
P {A = k} =
.
j
i
−
1
+
β
i − 1 + βn
n
2≤i≤N
i=1
2≤j <j <...≤N12i∈{j1 ,...,jk }Thus,
P {A = 1} =

N


i−1
i − 1 + βn
i=2

(2.3)

= 1 − β0 qHN −1 + O(q 2 ),where Hn is the nth partial Harmonic sum,
Hn ≡

n

1

i
i=1

.

14

www.nature.com/nature14doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Similarly,
P {A = 2} =

N


 i−1
βn
j − 1 + βn i=j i − 1 + βn
j=2

(2.4)

= β0 qHN −1 + O(q 2 ),while P {A = k} = O(q k−1 ) for k > 2. Thus, in the limit q → 0, we need only
consider populations composed of one or two neutral genotypes.
We next turn to α2,m , the expected number of unique phenotypes given m mutations and A = 2. If pi is the probability that there are N − i individuals of the ﬁrst
genotype and i of the second, the expected number of unique phenotypes is:
α2,m =

N−1


pi γm,i ,

(2.5)

i=1

where γm,i is the expected number of phenotypes, resulting from m mutations in a
population with (N − i) type 1 individuals and i type 2 individuals. Note that types
1 and 2 are equivalent (they both encode the wild-type phenotype) except that eachhas an independent set of K phenotypic neighbors.
We ﬁnd γm,i as before: let ξj,k (j = 1, 2, k = 1, . . . , i) be 1 if the k th mutation
to one of the type j individuals gives rise to a new phenotype. Without loss ofgenerality, we ﬁrst consider mutations in type 2. As above, the total number ofmutations to type 1 is
i

ξ1 =
ξ1,k
k=1

and E[ξ1 ] = α1,i . Now ξ2,k is 1 if and only if the k th mutation gives rise to diﬀerent
phenotype from all k − 1 previous mutations to type 2 individuals, and moreover,that phenotype was not produced by a mutation to a type 1 individual. For thelatter to be true, either the new phenotype is inaccessible to type 1, with probability
1− K, or the phenotype is adjacent to type 1, with probability K, but was not onePPof the ξ1 phenotypes that arose. Thus,k−1 
KK K − ξ1
1
+1−
P {ξ2,k = 1} = 1 −
K
P K
P
m−i
Thus, if ξ2 = k=1 ξ2,k ,K
K K − α1,i+1−
E [ξ1 + ξ2 ] = α1,i + α1,m−i
PKP
α1,i α1,m−i= α1,i + α1,m−i −,P15

www.nature.com/nature15

doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

and, since the m mutations are binomially distributed among the two genotypes,
j 
m−j m  
mi
α1,j α1,m−j i
γm,i =.
(2.6)
α1,j + α1,m−j −
1−
N
N
P
j
j=0

Finally, the probability that there are i individuals with the ﬁrst genotype, giventhat A = 2, can be determined via Hoppe’s urn. At each step, a new individualis added. The ﬁrst individual of the second genotype arrives at step k + 1 withprobability
βnβn +kN −i 1 ,k=1 βn +kjoining the k individuals of type 1 already present. Treating the k type 1’s and thesingle type 2 as separate lineages, the probability of i individuals of type 2 in the
total population of N is simply the proportion of partitions of N with one partitionof size i (Joyce & Tavaré 1987),
�Summing over all k gives
N −i

pi =

N −i−1
�Nk−1

−1
k

1
βn + k
k=1

=

1

.

−1 N −i�N −i−1
1
�Nk−1

−1
β
+
k
n
k
k=1

N −i

(N − i − 1)!(N − k − 1)!

HN −1 k=1

(N − i − k)!(N − 1)!

+ O(q)

(2.7)
(2.8)

Combining (2.1), (2.3) and (2.4), we have,∞βqm e−βq
D(q) =(P {A = 1} α1,m + P {A = 2} α2,m ) + O(q 2 )m!
m=0∞  m −β0β0 e=((1 − β0 qHN −1 )α1,m + β0 qHN −1 α2,m ) + O(q 2 )m!
m=0

while subtracting (2.2), dividing by q, and taking q → 0 yields∞
dD(q) β0m e−β0[β0 HN −1 (α2,m − α1,m ) − (m − β0 )α1,m ] .
=
dq q=0 m=0 m!

(2.9)

Supplementary Figure 2 plots Eq. (2.9) for P = 90, N = 300, and µ = 0.04 fora range of K. Also plotted are simulation results for ∆D
calculated at ∆q = 0.0001,∆q
for 100 million replicates for each point.
16

www.nature.com/nature16doi: 10.1038/nature08694

0dD
dq q=020

4060

SUPPLEMENTARY INFORMATION

0

20

40

K

60

80

Supplementary Figure 2 – Predictions using Eq. 2.9 (line) compared with themeans of 100 million replicate simulations at q = 0.0001 (circles).

3

Simulation Methods

Monte Carlo simulations were performed in two main ways. For the data in Figure 2& 3 in the main text and Supplementary Figures 4, matrices of transition probabilities
were pre-calculated and used to simulate individual replicates of a Markov chain.
Starting distributions for a two-allele Moran model were calculated using Eq. 3.58 in(Ewens 2004). These simulations were validated by comparison to individual-based
simulations of populations of N genotypes. Similar individual-based simulations werealso used to generate the data shown in Supplementary Figure 9.
Code for all simulations was written in C or C++ and compiled using gcc 4.0.1.
The GSL libraries (v.1.9) were used for pseudorandom number generation, special
functions, and probability distributions.

4

Quantitative eﬀects of N and µ on the optimallevel of robustness

We used numerical methods to explore how the adaptation time from steady state,
given by Eq. 1.11, and the mutational diversity, described above in Section 2, dependon the population parameters N and µ. The tables below show the approximate q
17

www.nature.com/nature17doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

which minimizes the adaptation time, or maximizes the mutational diversity, for agiven N and µ. Note that the eﬀects of N and µ can be approximated by considering
only their product, β, but that β aﬀects our measurements of evolvability in qualitatively diﬀerent ways: increasing β decreases the q which minimizes adaptation time,
but increases the q which maximizes mutational diversity.

N
5000
10000
10000
10000
15000
10000
30000

µ
0.001
0.00050.001
0.0015
0.001
0.003
0.001

β q of minimum time
5
0.67
5
0.68
10
0.66
15
0.64
15
0.65
30
0.58
30
0.6

q of maximum diversity
0.1
0.13
0.17
0.19
0.19
0.2
0.2

Table 1 – Eﬀects of population size and mutation rate on the q that minimizes adaptation time, and the q that maximizes mutational diversity. P = 100 and K = 5.

5

Relaxation of model assumptions

We relax each of the four main assumptions underlying the model presented in themain text: (1) a single neutral mutation redraws completely the phenotypes neighboring a genotype; (2) the number of phenotypes, K, in a genotype’s mutational
neighborhood is independent of its robustness, q; (3) K and q do not vary among
the genotypes on a neutral network and; (4) alternative phenotypes are lethal (oradaptive).
To relax the ﬁrst assumption we introduced a new parameter, f , deﬁned as thefraction of the K phenotypes neighboring a genotype which are redrawn after amutation. Our original model is therefore equivalent to this generalization whenf = 1. Considering this new parameter in light of the analysis above, it is clear thatthe only eﬀect of f is to scale the mutation rates between neutral genotypes that areor are not adaptable – i.e. class B �and class A genotypes, respectively. The mutation
rate from A to B is therefore µqf K. Since decreasing f is equivalent to increasingP
P , values of f < 1 actually broaden the range of parameters for which robustness andevolvability are positively correlated. This is illustrated in Supplementary Figure 3,
which shows analytic predictions of adaptation times for a range of f .
18

www.nature.com/nature18doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Supplementary Figure 3 – The predicted mean time before the arrival of a beneﬁcial
mutation with varying f . N = 10, 000, µ = 0.001, K = 5; a) P = 15; b) P = 100.
f = 1, bottom curves; f = 0.5, middle curves; and f = 0.1, top curves. Decreasing f
increases adaptation time, but the shape of the relationship with robustness remains
consistent.

The second assumption is an implicit description of the topology of neutral networks. In our original model, the fraction of mutations that are neutral (q) was
independent of the phenotypic diversity of non-neutral mutations (K). We relax thisassumption by allowing these two quantities to be correlated. As q increases, thenumber of mutational neighbors with distinct phenotypes might be expected to decrease, and so the richness of those phenotypic mutants, K, might also be expected todecrease. Therefore, in relaxing the second assumption we explored a linear, negative
relationship between q and an eﬀective value of K, denoted K̂(q):K̂(q) = K(1 − q)

(5.1)

This gives rise to an alternate limiting process


∂p
∂2
β
∂ βKq(1 − q)
= 2 [zp(t, z, y)] −p(t, z, y) − zp(t, z, y),
∂t
∂z
∂z
P
K
for which the expected arrival time, T1 (y) still takes the form of Equation (1.6),
given an appropriate choice of a and ν (see Section 1). This process also results
in a non-monotonic relationship between the arrival time of the adapted phenotype
and the robustness, q, but the mechanism is slightly diﬀerent. Here, the rate of
mutation from the adaptable types to the target phenotype is constant, while therate of mutation to the adaptable class approaches 0 for q near 0 or 1, and therefore
leads to a non-monotonic relationship between robustness and evolvability.
19

www.nature.com/nature19doi: 10.1038/nature08694●

K decreases with increasing q
K increases with increasing q
Constant K

●

●
●

●
●

●

●
●

●
●

200

500 1000

5000

Arrival Time for Beneficial Mutant

SUPPLEMENTARY INFORMATION

0.2

●

●
●

0.4

● ● ● ● ●

q

0.6

●

0.8

Supplementary Figure 4 – The mean time before arrival of a beneﬁcial mutation,with either K constant or K̂ correlated with q according to Eq. 5.1 or Eq. 5.2.
N = 10000, µ = 0.001, P = 200, and K = 20. Points are means of 10,000 replicate
simulations.

Supplementary Figure 4 shows the results of simulations when K is negatively correlated with q according to Eq. 5.1. Comparing the open circles to the closed circles,
which show numerical results for the same parameters and K independent of q, wenote a quantitative diﬀerence for larger values of q. However, the non-monotonicity
is preserved, suggesting that our qualitative conclusions are not sensitive even to astrong negative correlation between K and q.
For the sake of completeness, we also explored a positive correlation between q
and K:
K̂(q) = Kq + 1(5.2)
Simulations in this regime produced qualitatively similar results; see Supplementary
Figure 4.
To relax the third assumption, we performed simulations in which either K or q
20

www.nature.com/nature20doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

●

2500

Arrival Time for Beneficial Mutant

3000

were allowed to vary by genotype. (P is a property of an entire ﬁtness landscape andso cannot vary with genotype). To vary K, we assigned each genotype in our neutral network a value of K drawn independently from a distribution. K is therefore
eﬀectively re-drawn following a neutral mutation. Supplementary Figure 5 shows
the relationship between robustness, q, and mean adaptation time in three diﬀerent
situations: constant K = 5, K drawn from a Poisson distribution with mean 5,
and K drawn from a shifted geometric distribution with mean 5 (distributions weretruncated at K = P ). Populations were allowed to evolve before the environmental
shift until the number of adaptable individuals, and the distribution of K, had equilibrated. As the ﬁgure demonstrates, our results are remarkably robust even when
we allow K to vary across the neutral network. In fact, our analytic prediction forthe mean adaptation time is highly accurate, after replacing the ﬁxed value of K, in
the original model, with the mean of the distribution from which K is drawn, in theextended model.

1500

2000

●

Fixed K
Poisson−distributed K
Geometric−distributed K

●
●
●

●
●

●

●

●

●

0

500

1000

●

0.2

0.4

q

0.6

0.8

Supplementary Figure 5 – The mean time before arrival of a beneﬁcial mutation,
for constant K and variable K. N = 10000, µ = 0.001, P = 100, and mean K = 5.
Points show the means of at least 4,000 replicate simulations.

21

www.nature.com/nature21

doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATIONAs before, this can be explained by considering the rate of mutation into thefamily of adaptable genotypes and the rate of mutation from the adapted types tothe target phenotype. We omit the calculations here, but it can be shown thatthe rate of mutation into the adaptable class takes the form µq E[K], while thePrate of mutation into the target phenotype is obtained by averaging β(1−q)
over allK
individuals in the population. As before, these two rates determine the expectedarrival time of the target phenotype.
We also relaxed the third assumption by allowing q to vary among genotypes ina neutral network. To do so, we performed simulations in which an oﬀspring inherits
its parent’s value of q, plus a random perturbation whenever a neutral mutation
occurs. These perturbations are drawn from a Gaussian with standard deviation
0.1; if a perturbation would produce a q outside of [0, 1], then q is left unaltered.
This method of mutating q allows robustness to evolve. Prior theory suggests thata population should evolve towards elevated robustness (higher q) in a ﬁxed environment, because more robust individuals have an eﬀective selective advantage oforder µ (van Nimwegen et al. (1999), Forster et al. (2006)). This is indeed what
we observe: the mean q in a population, denoted q, tends to increase over time ina ﬁxed environment. Figure 6 shows this behavior by plotting the distribution of q
across the ensemble of replicate simulations. At the beginning of a simulation, eachindividual is assigned a q from the uniform distribution, and q is therefore approximately normally distributed with a mean of 0.5. After each population evolves for
10,000 generations in a ﬁxed environment, the ensemble distribution of q has shifted
signiﬁcantly towards larger values (Supplementary Figure 6), reﬂecting the evolutionof robustness in a ﬁxed environment.
After 10,000 generations, we assign one of the P alternative phenotypes the highest ﬁtness, and we record the subsequent arrival time of the ﬁrst such adaptive
mutation. The relationship between q in a population at the time of this environmental shift and the subsequent time before the arrival of a beneﬁcial mutation issummarized in Supplementary Figure 7. In this ﬁgure, the results of almost 300,000
replicates are binned according to q and plotted as grey points. The line illustrates
analytical results results with ﬁxed q, while black dots show comparable simulation
results for ﬁxed q. Although q may continue to evolve between the 10,000th generation and the eventual arrival of a beneﬁcial mutant, we nonetheless ﬁnd that q is
an excellent predictor of adaptation time. In fact, our analytic formula for the mean
adaptation time is highly accurate, after replacing the ﬁxed value of q in the original model with observed mean q in the population. We also performed simulations
in which we waited 100,000 generations before allowing beneﬁcial mutations; these22

www.nature.com/nature22doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Count

30000
20000

Count

0

10000
0

10000 20000 30000 40000 50000

q at 10000 generations

40000

q at 0 generations

0.46

0.48

0.50

q

0.52

0.54

0.0

0.2

0.4

q

0.6

0.8

1.0

Supplementary Figure 6 – (a) Histograms of q among nearly 300,000 replicatesat generation zero (left) and generation 10,000 (right). Parameters: N = 10000,µ = 0.001, P = 100, and K = 5.

simulations produced very similar results. Thus, our results are robust with respect
to signiﬁcant variation in q across a neutral network and the associated evolution ofq within a population.
For the sake of completeness, we also allowed f to vary across the neutral network.
We assigned each genotype in our neutral network a value of f drawn independently
from a distribution. f is therefore eﬀectively re-drawn following a neutral mutation.
Supplementary Figure 8 shows the relationship between robustness, q, and mean
adaptation time in two diﬀerent situations: constant f = 0.5, f drawn from a normal
distribution with mean 0.5 and variance 0.1 (truncated at zero and one). Populations
were allowed to evolve before the environmental shift until the number of adaptable
individuals, and the distribution of f , had equilibrated. As the ﬁgure demonstrates,
our results are robust even when we allow f to vary across the neutral network. In
fact, our analytic prediction for the mean adaptation time is highly accurate, after
replacing the ﬁxed value of f , in the original model, with the mean of the distribution
from which f is drawn, in the extended model.
Finally, we considered the eﬀects of relaxing the fourth assumption. In our original model all mutations were either neutral or lethal, prior to the environmental
shift. We simulated an alternative model in which the previously-lethal phenotypes

23

www.nature.com/nature23doi: 10.1038/nature08694

Arrival Time to Beneficial Mutation

SUPPLEMENTARY INFORMATIONq
Supplementary Figure 7 – Mean arrival time of the ﬁrst beneﬁcial mutant. The grey
crosses depict means binned according to the value of q at the time of the environmental
shift. The solid line shows comparable analytical results for ﬁxed q, while the black
circles show simulations for ﬁxed q. N = 10000, µ = 0.001, P = 100, and K = 5.

were assigned ﬁtness 1 − s, reﬂecting a disadvantage of size s compared the ﬁtnessof the focal phenotype. As before, individuals are chosen to reproduce in proportion
to their ﬁtnesses.
To describe mutation among the phenotypic classes in this more general model
it is helpful to extend the notation used above. Genotypes of type C express thephenotype that will be beneﬁcial after the environment changes, and they have ﬁtness
1 − s prior to that shift. Genotypes of type B can produce type C mutants, andthey can have ﬁtness 1, which we will call Bf it , or 1 − s, which we denote by Bdel .
Similarly, type A genotypes cannot produce type C, and they may either be Af it orAdel . Any mutation, whether neutral, deleterious, or beneﬁcial, causes the set of K
neighbors to be redrawn. Therefore, Adel -types may or may not be able to mutate
back to Af it , depending on the phenotypic neighborhood of that type. Finally,
we impose a necessary constraint on phenotypic neighbors: a mutant’s phenotypicneighborhood must contain its parent’s phenotype. So, if a type C genotype mutates24

www.nature.com/nature24doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Supplementary Figure 8 – The mean time before arrival of a beneﬁcial mutation,
for constant f and variable f . N = 10000, µ = 0.001, P = 100, and mean f = 0.5.
Points show the means of at least 4000 replicate simulations.

to a deleterious phenotype, it must become Bdel .
We used individual-based Moran simulations to explore the time elapsed until
the ﬁrst type-C individual arises. Each simulation began with an average proportionK
individuals of type Bf it , with the remaining of type Af it . Populations evolved forP
30,000 generations before the environmental shift, at which time type C individuals
become beneﬁcial. 30,000 generations was determined to be suﬃcient for equilibriumby observing the distributions of Bf it , Bdel , and type C individuals across an ensembleof evolving populations. We recorded the arrival time of the beneﬁcial type in termsof the number of generations after the environmental shift; if type C individuals were
present at the environmental shift, this time was recorded as zero.
When alternative phenotypes are strongly deleterious (s ≈ 1) we expect to see
little diﬀerence between the results of the original, lethal-mutation model and the
new model. However, when s ≈ 0 we anticipate a monotonic, negative relationshipbetween robustness and evolvability. When s = 0 non-neutral mutations have nocost, contribute to diversity in phenotypic neighborhoods, and generate type C individuals. Therefore, we ask how large must s be to match qualitatively the resultsof our original, lethal-mutation model.
Supplementary Figure 9(a) conﬁrms that the relationship between robustness andevolvability is still non-monotonic when alternative phenotypes are moderately dele25

www.nature.com/nature

25

SUPPLEMENTARY INFORMATION

(b)

●

●

●

●
●

●

●

500

10

1000

●
● ●

15

s = 0.02
s = 0.05
s = 0.3●(a)Distinct Mutant Phenotypes

2000

Arrival Time for Beneficial Mutant

doi: 10.1038/nature08694

200

●

●

5

●

●

●
●

●

●
●

0.2

●

●

0.4

●

q

●

0.6

0.8

●

s = 0.02
s = 0.05
s = 0.3

0.2

●

0.4

q

0.6

0.8

Supplementary Figure 9 – Adaptation time and phenotypic diversity when alternative phenotypes are deleterious, but not lethal. N = 10, 000, µ = 0.001, P = 100, andK = 5 for all simulations shown. The dotted lines displays the means for the originalmodel with lethal alternative phenotypes. (a) Mean adaptation time for three valuesof s, the ﬁtness penalty of alternative phenotypes. Each point is the mean of 10,000
replicate simulations. (b) Mean number of unique mutant phenotypes. Each point is
the mean of 5000 simulations.

terious. Similarly, we also used Moran simulations to measure phenotypic diversity
in steady-state populations when the alternative phenotypes were deleterious but
not lethal. These results, shown in Supplementary Figure 9(b), demonstrate thatrobustness and phenotypic diversity can also exhibit a non-monotonic relationshipwhen alternative phenotypes are moderately disadvantageous. While the assumptionthat other phenotypes are lethal is mathematically convenient in our analyses, it canbe relaxed without changing our qualitative results.

6

RNA Simulations

6.1

Measuring Epistasis in the RNA Landscape

We used the Vienna RNA package to calculate the minimum-free-energy secondary
structures of RNA nucleotide sequences. The resulting genotype-phenotype map

26

www.nature.com/nature26

doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

was analyzed to estimate P, K, and f for a range of sequence lengths, L. While there
are signiﬁcant diﬀerences between the RNA landscapes and the type of abstract
landscapes considered in our model, these estimated values provide some intuition
about reasonable parameter values.
Our model assumes that out of P equally-likely phenotypes, K are equallyaccessible from a given genotype. Since phenotypes are not evenly distributed in
the RNA landscape, either globally or in the neighborhood of a single genotype, wemust calculate some eﬀective number of types for RNA. We note that, in our model,
if we introduce two mutations into individuals of the same genotype, then the probability these mutation produce the same phenotype is K1 . Similarly, if the mutations
occur in individuals with diﬀerent genotypes, the probability of producing the same
phenotype is P1 . Since these probabilities have an intuitive connection to the rolesof K and P in our model, we deﬁne the eﬀective K and P in the RNA landscapeto yield the same probabilities. Additionally, some RNA sequences have the trivial,
unfolded shape as their minimum free energy structure. We choose to regard thisunfolded phenotype as inviable. We therefore deﬁne Ke as the inverse of the probability, pclone , that two viable phenotypes produced by non-neutral mutations in clonesof the same genotype, are the same. Similarly, we deﬁne Pe as the inverse of theprobability, prandom , that two viable phenotypes produced by non-neutral mutationsin individuals with random, viable genotypes, are the same. Both pclone and prandom
are measured by sampling random genotypes of length L with viable structures, andrecording the phenotypes of random, non-neutral mutants of those genotypes; one
hundred thousand genotypes were sampled for each L. Ke is therefore an average K
across all genotypes in the landscape, while Pe is a lower bound on the number ofaccessible phenotypes; both are plotted in Supplementary Figure 10 for a range ofL.
Estimating f is complicated by variation in K across the network. In keepingwith the approach for Ke and Pe above, we relate f to the probability that two
mutants have identical phenotypes. In our model, f determines the expected fractionof K neighbors that diﬀer between two immediate neighbors on a neutral network,for P large. We therefore measure the probability, pneighbor , that a viable, nonneutral mutant of genotype x has the same phenotype as a viable, non-neutral mutant
produced from a neutral neighbor of x. pneighbor can be related to f , pclone , andprandom as follows. First, note that we expect an adjacent pair of genotypes tohave (1 − f )K phenotypic neighbors in common, and f K neighbors drawn from theremaining P − (1 − f )K possible phenotypes. A given phenotype can appear only
once in a set of K neighbors, so if a phenotype falls among the shared (1 − f ) portionof the neighborhoods, then it cannot fall among the unshared, redrawn f portion
27

www.nature.com/nature27

50000SUPPLEMENTARY INFORMATION

(a)

(b)

Pe
20

40

60

80

Sequence Length L

100

100

10

500

20

30

2000

Ke

40

10000

50

60

doi: 10.1038/nature08694

20

40

60

80

Sequence Length L

100

Supplementary Figure 10 – Ke and Pe for Simulated RNA Molecules of Length Lof either genotype’s neighborhood. Therefore, we only need to consider the cases
where both mutants produce phenotypes in the shared part of their neighborhoods,
with probability (1 − f )2 , or both produce phenotypes in the unshared portion of
their neighborhoods, with probability f 2 . In the former case, the probability thatcloneboth mutants will have the same phenotype is (1−f1)Ke , or p1−f
. In the latter case,each mutant can be one of P − (1 − f )K possible phenotypes; if K  P then theprobability that they are the same phenotype is approximately P1e , or prandom . These
considerations yield the following relationship among pneighbor , pclone , prandom , and f :pneighbor ≈ (1 − f ) pclone + f 2 prandom(6.1)

We use the equation above to deﬁne the eﬀective value of f for the RNA landscape, after ﬁrst measuring pneighbor , pclone , prandom . Supplementary Figure 11 shows
values of fe estimated in this way, for several sequence lengths L. These estimatescomplement the ﬁndings of (Sumedha et al. 2007), who examined the overlap inphenotypic distance of related, but not neighboring, RNA genotypes.
Together these estimates for the RNA folding landscape suggest that the spaceof possible phenotypes is very rich compared to the phenotypic neighborhood of asingle genotype – i.e. K  P . Furthermore, the phenotypic neighborhoods of neutral neighbors diﬀer signiﬁcantly. These measurements in RNA conﬁrm the central
28

www.nature.com/nature28

doi: 10.1038/nature08694

0.26

0.28

fe

0.30

0.32

0.34

SUPPLEMENTARY INFORMATION

20

25

30

35

Sequence Length L

40

Supplementary Figure 11 – fe for simulated RNA molecules of length L.

premise of our model: a neutral mutation can have profound, epistatic consequences
for the phenotypes of subsequent mutations.

6.2

Robustness and Evolvability in Evolving RNA Populations

To more directly test our abstract model using the RNA genotype-phenotype landscape, we performed evolutionary simulations with populations of RNA sequences.
In these simulations, genotypes were 24-base-pair sequences and phenotypes werethe minimum-free-energy structures of these sequences. For a set of simulations, weﬁrst chose M structures to designate as ‘high-ﬁtness;’ these correspond to the single,
high-ﬁtness type in our abstract model. Adjusting M allows us to tune the easeof adapting to high ﬁtness, much as we set the ratio K
in our general model. For
Peach replicate, we chose genotypes at random until one was found that folds into a
non-trivial structure which is not a high-ﬁtness type. The chosen phenotype represents the ﬁt wild-type; all other structures, except the high-ﬁtness types, are thenconsidered to be inviable. The chosen genotype is then used to found a population of
N clones, which then evolve according to a Wright-Fisher, discrete-generation model
until the ﬁrst high-ﬁtness phenotype arises. The ’adaptation time’ is deﬁned as the29

www.nature.com/nature29doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

number of generations which pass before a high-ﬁtness type arises. Note that each
replicate begins from a potentially unique genotype, but the high-ﬁtness phenotype
remain consistent throughout a set of replicates.
We measured the robustness, q, of each initial genotype as the frequency of point
mutations that do not alter the phenotype, as we have done throughout our analysis.
We used this initial q to predict adaptation times according to the analytical formulaederived above. For these simulations, the initial population was clonal, and so either
all individuals were pre-adapted or none were pre-adapted. Thus, we obtained the
expected adaptation time by suming expected adaptation times from each of these
initial conditions, weighted by their probability:
P{none pre-adapted}T N (0) + P{all pre-adapted}T N (N )∼ P{none pre-adapted}T1 (0) + P{all pre-adapted}T2 (1),for large N , where T1 (0) is given by (1.10) and T2 (1) is obtained from (1.5). To calculate predicted waiting times as a function of q, we measured three properties of these
initial genotypes corresponding to the three parameters of our general model: theprobability that a genotype is pre-adapted, the average mutation rate of pre-adapted
genotypes to the high-ﬁtness types, and the probability that a genotype which is not
pre-adapted has a pre-adapted neighbor. These measurements correspond to K
, 1,
P K
and f in our abstract model. Although we estimated f above for RNA landscapesin general, we found that speciﬁc landscapes, determined by the particular set of M
high-ﬁtness types, have values of f which are more relevant to evolution of thoselandscape. Also, the mutation rate to high-ﬁtness types declined monotonically withq in the RNA data, and so we used a linear regression of this relationship in our
analytic predictions.
Supplementary Figure 12 shows the mean adaptation time as a function of robustness, q, for two sets of RNA simulations. In both cases, N = 500 and µ = 0.0002;
M = 500 in the lower curve, and M = 200 in the upper curve. Note that the empirical relationship between robustness and evolvability is somewhat more complex
and diﬀers from our analytical prediction; however, the primary trend matches the
U-shaped prediction of our general, population-genetic model. Moroever, our general model accurately predicts the scale of the waiting times in the RNA simulations,over several orders of magnitude.
Our analysis of the RNA genotype-phenotype map suggests that robustness canfacilitate adaptation under some circumstances. However, Ancel & Fontana (2000)
found that evolving RNA populations became less able to adapt as their robustness
increased. While several diﬀerences between their study and ours may contributeto this discrepancy, the most fundamental one is that Ancel and Fontana did not30

www.nature.com/nature30

Mean adaptation time

doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

10000
5000

1000
500

100
0.0

0.2

0.4

0.6

0.8

1.0

Robustness, q
Supplementary Figure 12 – The mean time before arrival of a beneﬁcial mutation in
RNA simulations. Each point is the mean of replicate simulations, and lines represent
analytical predictions as described in the text. N = 500 and µ = 0.0002; M , thenumber of high-ﬁtness phenotypes, is 500 in the lower curve, and 200 in the upper.

manipulate robustness to test its eﬀect on adaptation. In fact, populations typicallybecame more robust as they evolved, such that robustness and distance from theoptimum were confounded. By specifying a set of high-ﬁtness types, rather thanevolving population toward a distant optimum, our simulations reduce the eﬀectof confounding variables. We also note that Cowperthwaite et al. (2008) found no
relation between the size of the neutral network on which a population began, andits success at evolving to a distant phenotypic optimum; similarly, Wagner (2008)
found no relationship between the robustness of a sequence, and a measure of its
evolvability. The major diﬀerence between our results and those negative ﬁndings isthat we designate a set of structures as high-ﬁtness, as opposed to picking a single,
optimal structure. As a result, the populations evolved much more rapidly in our
simulations, compared with the slow, step-wise process of adaptation to a distant
optimum. The relative speed of adaptation in our simulations suggests that therobustness of the initial genotype has a much greater inﬂuence on adaptation thanin previous studies.
31

www.nature.com/nature31doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

For all simulations we used version 1.6.1 of the Vienna package. Sequences were
folded at 37◦ C with the option ‘dangles’ set to ‘2.’ These settings produced the samelandscapes as those described in (Cowperthwaite et al. 2008).

References
Abramowitz, M. & Stegun, I. A. (1965), Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables, Dover Publications, New York.
Ancel, L. W. & Fontana, W. (2000), ‘Plasticity, evolvability, and modularity in rna.’,
J Exp Zool 288(3), 242–283.
Cowperthwaite, M. C., Economo, E. P., Harcombe, W. R., Miller, E. L. & Meyers,
L. A. (2008), ‘The ascent of the abundant: How mutational networks constrain
evolution’, Plos Computational Biology 4(7).
Ethier, S. N. & Kurtz, T. G. (2005), Markov Processes: Characterization and Convergence, Wiley Interscience, Hoboken.
Ewens, W. J. (2004), Mathematical population genetics, Vol. v. 27, 2nd ed edn,
Springer, New York.
URL: http://www.loc.gov/catdir/enhancements/fy0818/2003065728-d.html
Feller, W. (1954a), ‘Diﬀusion processes in one dimension’, Trans. Amer. Math. Soc.
77(1), 1–31.
Feller, W. (1954b), ‘The general diﬀusion operator and positivity preserving semigroups in one dimension’, Ann. of Math. 60(3), 417–436.
Forster, R., Adami, C. & Wilke, C. O. (2006), ‘Selection for mutational robustness
in ﬁnite populations’, Journal of Theoretical Biology 243(2), 181–190.
Joyce, P. & Tavaré, S. (1987), ‘Cycles, permutations and the structure of the Yuleprocess with immigration’, Stoch. Proc. Appl, 25, 309–314.
Karlin, S. & Taylor, H. M. (1981), A Second Course in Stochastic Processes, Academic Press, San Diego.
Slater, L. J. (1966), Generalized hypergeometric functions, University Press, Cambridge.

32

www.nature.com/nature32doi: 10.1038/nature08694

SUPPLEMENTARY INFORMATION

Sumedha, Martin, O. C. & Wagner, A. (2007), ‘New structural variation in evolutionary searches of rna neutral networks’, Biosystems 90(2), 475–485.
van Nimwegen, E., Crutchﬁeld, J. & Huynen, M. (1999), ‘Neutral evolution of mutational robustness’, Proceedings of the National Academy of Sciences of the UnitedStates of America 96(17), 9716–9720.
Wagner, A. (2008), ‘Robustness and evolvability: a paradox resolved.’, Proc Biol Sci
275(1630), 91–100.

33

www.nature.com/nature33