O R I G I NA L A RT I C L E
doi:10.1111/j.1558-5646.2011.01526.x

ENVIRONMENTAL ROBUSTNESS AND
THE ADAPTABILITY OF POPULATIONS
Alexander J. Stewart,1,2, Todd L. Parsons,1 and Joshua B. Plotkin1,3
1

Department of Biology, University of Pennsylvania, Philadelphia, 19104, Pennsylvania

2

CoMPLEX, University College London, Physics Building, Gower Street, London WC1E 6BT, United Kingdom
3E-mail: jplotkin@sas.upenn.edu

Received March 30, 2011
Accepted November 15, 2011
Recent work has shown that genetic robustness can either facilitate or impede adaptation. But the impact of environmentalrobustness on adaptation remains unclear. Environmental robustness helps ensure that organisms consistently develop the same
phenotype in the face of “environmental noise” during development. Under purifying selection, those genotypes that express the
optimal phenotype most reliably will be selectively favored. The resulting reduction in genetic variation tends to slow adaptation
when the population is faced with a novel target phenotype. However, environmental noise sometimes induces the expressionof an alternative advantageous phenotype, which may speed adaptation by genetic assimilation. Here, we use a populationgenetic model to explore how these two opposing effects of environmental noise influence the capacity of a population toadapt. We analyze how the rate of adaptation depends on the frequency of environmental noise, the degree of environmentalrobustness in the population, the distribution of environmental robustness across genotypes, the population size, and the strengthof selection for a newly adaptive phenotype. Over a broad regime, we find that environmental noise can either facilitate or impede
adaptation. Our analysis uncovers several surprising insights about the relationship between environmental noise and adaptation,and it provides a general framework for interpreting empirical studies of both genetic and environmental robustness.

KEY WORDS:

Adaptation, mutations, population genetics.

Robustness and adaptability are fundamental and seemingly conflicting properties of biological systems (de Visser et al. 2003;
Lenski et al. 2006; Wagner 2008; Draghi et al. 2010). An organism is robust if it can reliably produce a given phenotype acrossa range of environmental and genetic perturbations. But adaptability requires that alternative phenotypes be expressed so thata population can adapt to new conditions. Naively, it seems thatthe more robust an organism is, the less frequently it expresses
alternative phenotypes, and so the less adaptable (or evolvable)
the population will be. However, recent work demonstrates thatgenetic robustness can either facilitate or impede a population’s
capacity to adapt to a novel environment (Wagner 2008; Draghi
et al. 2010). The relationship between adaptation and environmental robustness—that is, the robustness of an individual’s expressedphenotype against nonheritable perturbations—remains unclear.

There are reasons to believe that environmental noise may facilitate adaptation. It is well known that genetic and environmental
perturbations often have similar phenotypic effects on an organism. As a result, occasionally expressing alternative phenotypes
by environmental noise might help speed adaptation via “genetic
assimilation”—that is, by promoting those individuals who willsubsequently evolve the adaptive phenotype by mutation (Ancel
and Fontana 2000; Rutherford 2003; Earl and Deem 2004; Gibson
and Dworkin 2004; Masel 2005; Meyers et al. 2005). But there
are also reasons to believe that environmental noise may impede
adaptation: noise will select for the most environmentally robust
genotypes, thus reducing neutral standing genetic variation in agiven environment and thereby slowing adaptation to a new target phenotype (Hermisson and Wagner 2004; Wagner et al. 1997;
van Nimwegen et al. 1999; Masel and Trotter 2010). Therefore,C 2012 The Author(s). Evolution 
C 2012 The Society for the Study of Evolution.

1598

Evolution 66-5: 1598–1612

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

we might expect a complex relationship between environmental
noise, environmental robustness, and the rate at which a population adapts to a novel environment. It is this relationship that is
the focus of this article.

of phenocopy (e.g., a chaperone-mediated capacitor), but rather
we will explore in greater generality how environmental robustness influences a population’s capacity to adapt, assuming thatgenetic and environmental perturbations tend to produce similar
phenotypes.

Phenocopy and Genetic Assimilation

Genetic and environmental perturbations often have similar phenotypic effects—a phenomenon known as “phenocopy.” Phenocopy was first discovered in the classic studies of Goldschmidt(1935) and Waddington (1953a), and it has since been observed
at many different levels of biological organization, including
protein-coding sequence, mRNA and protein structures, individual gene expression patterns, regulatory or metabolic network
dynamics, and entire developmental pathways (Table 1).
The classic example of the phenocopy concept comes fromstudies of heat shock in Drosophila (Waddington 1942, 1953a,b,
1959). Waddington found that the type of phenotypic variation induced by heat shock could subsequently be produced by genetic
changes accrued over several generations of artificial selection(Waddington 1953a). In other words, a phenotype that was initially produced as the result of an environmental stimulus was
subsequently produced in the absence of that stimulus as the result of genetic changes—a process called “genetic assimilation.”
The mechanistic basis for genetic assimilation in these heat shock
experiments has been traced, in part, to the behavior of the molecular chaperone Hsp90 (Rutherford and Lindquist 1998; Chow andChan 1999; Rutherford 2003; Sangster et al. 2004), which buffers
the effects of standing genetic variation by aiding protein folding. When the function of Hsp90 is impaired, either through heatshock or through mutation, previously cryptic genetic variation isrevealed as phenotypic variation. Any such phenotypic variationthat is adaptive can quickly become fixed through genetic assimilation (Ancel and Fontana 2000; Gibson and Dworkin 2004;
Jarosz and Lindquist 2010; Masel and Trotter 2010).
Environmental perturbations can also lead to phenocopy in aconstant genetic background (i.e., even without standing cryptic
genetic variation). One of the most important examples of thisoccurs in the secondary structures of RNA and protein molecules(Ancel and Fontana 2000; Bloom et al. 2006). Although the majority of genetically identical sequences fold to the same, minimum
energy structure, some adopt higher energy structures. The alternative structures that a molecule can assume generally correspondto the minimum energy structures associated with the genotypes
available by mutation to that molecule (Ancel and Fontana 2000;
Bloom et al. 2006). The degree to which these alternative phenotypes are expressed depends on environmental factors, such astemperature. In this case, phenocopy occurs on a constant genetic
background, without the need for cryptic genetic variation.
In this study, we will not specify a specific phenotype ofinterest (e.g., RNA secondary structure) or a specific mechanism

Variation in the Degree of Environmental
Robustness among Genotypes

The degree to which an individual is robust to environmental noise
may depend upon the individual’s genotype. When a populationis under stabilizing selection for a given phenotype, and whengenotypes vary in their environmental robustness, selection willtend to concentrate the population at the highly robust genotypes.
The strength of this effect will increase as the frequency of environmental perturbations increases. As a result of this effect, thepopulation will harbor less genetic diversity when environmental
perturbations occur more frequently, which in turn will reduce
the speed with which the population adapts to a new target phenotype. Thus, genotypic variation in environmental robustness
can influence the adaptability of a population (Masel and Trotter
2010). (We use the term adaptability, which is more accurate inthis context than the term evolvability). Standing genetic diversity
increases the speed of adaptation both by increasing the probability that the population harbors adaptable individuals at the timeof the environmental shift and by increasing the rate at whichsuch individuals are produced following the environmental shift.
Therefore, standing genetic diversity is important for adaptation
when adaptable genotypes are rare.
Presentation Outline

We will construct a general, population-genetic model to explore
the effects of environmental robustness on adaptability. Our model
assumes the concept of phenocopy and it also allows for variationin environmental robustness across genotypes. On the one hand,the phenocopy phenomenon has been shown to aid the adaptationof a population following an environmental shift in the presence ofenvironmental noise (Ancel and Fontana 2000; Rutherford 2003;
Earl and Deem 2004; Gibson and Dworkin 2004; Masel 2005;
Meyers et al. 2005). But on the other hand, as discussed above,
variation in environmental robustness across a neutral network
tends to reduce genetic variation prior to an environmental shift(i.e., under stabilizing selection) when noise is present, and thus
slows adaptation (Wagner et al. 1997; Hermisson and Wagner
2004; Masel and Trotter 2010). These two forces therefore have
opposite influences on the capacity of a population to adapt toa new environment. Here, we resolve how these opposing forcesinteract to determine the adaptability of populations subject toenvironmental noise.
The article is structured as follows. We begin by introducing a general model for the evolution of a population on an

EVOLUTION MAY 2012

1599

A L E X A N D E R J. S T E WA RT E T A L .

Table 1.

Examples of phenocopy observed in different biological systems.

Phenotype

Genetic perturbation

Environmental perturbation

References

Protein misfolding

Mutations at protein-coding
regions
Mutations at stop codons

Errors in transcription ortranslation
Accidental read through ofstop codons
Temperature induced noise

Goldsmith and Tawkik (2009),
Meyerovich et al. (2010)
Masel and Bergman (2003)

Incorrectly transcribedmRNA
Protein/mRNA misfolding
Protein/mRNA misfolding
Incorrect protein expression
level
Incorrect gene expression
pattern

Incorrectly executed
developmental/ metabolic
pathway

Mutations to coding regions
Knockout of molecularchaperone Hsp90
Mutations to regulatory
regions
Gene knockout resulting inloss of buffering in gene
networkHeat shock

Mutation resulting in loss ofredundancy in pathways

Large environmental
perturbation

arbitrary neutral network, with an arbitrary distribution of environmental robustness across the network. We consider a population initially adapted to one environment (one target phenotype)
and then subsequently exposed to a new environment (a different
target phenotype). We treat our model analytically to recover thefull probability distribution of the time required to adapt to the
new target phenotype. We next discuss a simplified neutral network with constant mutational robustness in which we separate
genotypes into classes of high and low environmental robustness.
This simple model is used to develop intuition for how environmental robustness affects adaptability. Finally, we summarize ourresults and discuss future directions, including how our framework might be used to interpret systematic studies in quantitative
genetics.

A Population-Genetic Model
Our general model is based on the Moran process from population genetics, in which each haploid individual in a populationof constant size reproduces at a rate determined by its relative
fitness, with its offspring replacing a randomly chosen individual.
During reproduction, the offspring may acquire a mutation, whichchanges its genotype from that of its parent. Different genotypes
may encode for different phenotypes; fitness is determined bycomparing an individual’s phenotype to the phenotype selectedby the current environment, as described below.
In a given environment that selects for a specific phenotype,
we will focus our analysis on the set of all genotypes that encodethat phenotype. The set of such genotypes is called a “neutral
network” because mutations among them do not change an individual’s encoded phenotype. The structure of mutations among

1600

EVOLUTION MAY 2012

Intrinsic/extrinsic noise
Large environmental
perturbation

Ancel and Fontana (2000),
Bloom et al. (2006)
Rutherford and Lindquist(1998)
Raj and van Oudenaarden(2008)
Levy and Siegal (2008), Siegal
and Bergman (2002),
Lopez-Maury et al. (2008),
Wagner and Wright (2007)
McLoughlin and Copley(2008), Tawkik (2010)

the genotypes on the neutral network is defined by the adjacency
matrix M, which has entries m i j = 1, if genotype i can mutate to
genotype j without changing the phenotype, and m i j = 0 otherwise. We define the mutational robustness, qi , of genotype i as
the proportion of all mutations that do not change its phenotype.
More specifically, we defineqi =1
m ij ,
L j(1)where L denotes the total number of different mutations that could
arise in an individual (e.g., the length of its binary genome). The
mutational robustness, qi , therefore quantifies the chance that amutant offspring of an individual with genotype i will remainon the neutral network (i.e., will still express the same phenotype), or not. Each genotype also has an associated environmentalrobustness, φi , discussed in more detail below.
When an individual of genotype i is chosen to reproduce,
a mutation occurs with probability μ. The mutation rate μ istherefore the genome wide mutation. We have typically chosen amutation rate of μ = 10−3 in our simulations, corresponding to a
genome of about 105 nucleotides (see Supporting Information). If
a mutation occurs, the resulting offspring remains on the neutral
network (i.e., encodes the same phenotype) with probability qi .
In this case, the offspring nonetheless has a new genotype, which
is drawn uniformly among i’s neutral neighbors (i.e., from theset of genotypes j with m i j = 1). With probability 1 − qi , on theother hand, the resulting offspring encodes a different phenotype
and thus lies outside of the neutral network.
Following the work of Draghi et al. (2010), we assume
that each genotype has a set of K alternative phenotypes that

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

constitute its phenotypic neighborhood—that is, the set of possible phenotypes that can be produced by a nonneutral mutation.
These K phenotypes are drawn uniformly from a total of P possible alternative phenotypes. We make the simplifying assumptionthat phenotypic neighborhoods are independent, such that theK phenotypes accessible to a genotype are redrawn whenever amutation occurs. When a mutation produces an individual lying
outside the neutral network, its resulting phenotype is drawn fromthe K alternatives that belong to its phenotypic neighborhood. The
form of genotype–phenotype map described above corresponds toevolution in a space of infinite genotypes, such that each mutation
results in an entirely new genotype entering the population. Ourmodel specifies the effects of such mutations on the phenotype,using summary statistics such as K , q, and P. Such a statistical description of the genotype–phenotype map provides a very
good approximation to evolution on an explicit, finite genotype–phenotype map, such as the RNA folding landscape (Draghi
et al. 2010). Furthermore, we demonstrate in the Supporting Information that our results hold even when we relax the statisticalassumptions described here and allow for correlations between
the phenotypic neighborhoods of neighboring genotypes.
As described so far, our model is identical to the one exploredby Draghi et al. (2010), and it can be used to study how mutational robustness influences the process of adaptation. To study
environmental robustness, we will extend the model by assumingthat when an individual produces an offspring, the offspring experiences an environmental perturbation with probability . Such anenvironmental perturbation may result in the offspring developing a different phenotype from the “intended” phenotype. In otherwords, when a perturbation occurs, an offspring may express a
different phenotype than the one encoded by its genotype. The perturbation is specific to the development of that offspring alone,
and it does not affect other individuals in the population. The
parameter  therefore measures the amount of noise that occurs
during the development of an individual, which we call environmental noise. If such an environmental perturbation occurs during
a reproduction event, the resulting individual adopts its intended
phenotype with probability φi . With probability 1 − φi , however,
the offspring is not robust to the environmental perturbation andthe offspring adopts an alternative phenotype. The parameter φitherefore describes the robustness of genotype i against environmental noise.
We assume that the alternative phenotypes that result frommutation in an individual and those that result from environmental perturbation are the same—so that when an environmental
perturbation produces an individual with an alternative phenotype, that alternative phenotype is drawn from the K phenotypes
that constitute its phenotypic neighborhood (i.e., the same setof phenotypes the individual can produce by mutation). This assumption captures the concept of phenocopy—namely that ge-

netic and environmental perturbations tend to produce a similar set of phenotypes—whose longstanding empirical support wehave described above.
We consider a haploid population of N asexuals that reproduce according to the Moran process and occupy a neutral
network of the type described above. We assume that initially thephenotype associated with the focal neutral network is optimalfor the environment, and we assign it fitness one without lossof generality. All alternative phenotypes are strongly deleterious,and assigned fitness zero (or, equivalently, any selective deficit
 1/N [van Nimwegen et al. 1999; Draghi et al. 2010]). After along period of stabilizing selection, we assume that an environmental shift occurs, such that one of the P alternative phenotypes
now has a selective advantage, s, over the initial phenotype. The
newly adaptive phenotype is chosen uniformly from among the P
alternative phenotypes. Once the environmental shift occurs, westudy the adaptation time—defined as the time required for the
population to produce an individual whose genotype encodes the
newly adaptive phenotype. (Once one such adaptive genotype is
produced, the subsequent probability of and time to its fixation
are well understood for asexuals, according to classical population genetics [Ewens 2004]; however, explicit calculation of thisquantity may be difficult or impossible in the presence of clonalinterference or in sexual populations). The adaption time is a random variable whose properties depend on how the population is
distributed across the neutral network, and on the extent to which
phenocopy speeds adaptation by occasionally producing individuals with the adaptive phenotype. As the rate of environmental
perturbation, , increases, the first of these two effects, increased
selection for genotypes with high environmental robustness, will
reduce genetic diversity and make it more and more difficult forthe population to explore the neutral network, thus slowing adaptation. However, the second of these two effects, phenocopy, will
result in more frequent expression of the adaptive phenotype,which will tend to speed adaptation by promoting those individuals who are a single mutation away from the adaptive phenotype.
How strongly these two opposing effects influence adaptation
time, and under what conditions, is the central problem we seek
to resolve.
Following Draghi et al. (2010), we separate the neutral
network into adaptable genotypes and nonadaptable genotypes.
Adaptable genotypes are one mutation away from the newly adaptive phenotype, whereas nonadaptable genotypes do not containthe newly adaptive phenotype in their phenotypic neighborhoods(see Fig. 1 ). In the presence of environmental perturbations  > 0(and provided environmental robustness is not complete, i.e., provided φ = 1), adaptable genotypes gain a direct selective advantage in our model, compared to nonadaptable individuals withthe same environmental robustness, because they occasionallyexpress the newly adaptive phenotype.

EVOLUTION MAY 2012

1601

A L E X A N D E R J. S T E WA RT E T A L .

Figure 1. Depiction of genotypes (dots) and phenotypes (colors).

In this example, the initially adaptive phenotype is indicated byred, and the newly adaptive phenotype following the environmental shift is indicated by blue. Dark red indicates genotypes
that are adaptable (i.e., those that can express the adaptive phenotype either by environmental noise or genetic mutation); light
red indicates genotypes that are nonadaptable. All other phenotypes, shown here as black or green, are deleterious. Althoughsome genotypes are shown as only neighboring other genotypesof the same color, in reality, all genotypes will always have neighbors that code for a different phenotype.

Analysis
We focus on the situation in which adaptable individuals are rareat the time of the environmental shift, meaning that the number√of adaptable individuals is O( N ) or smaller. We further assume
that the probability of environmental perturbation in an individual,
, is sufficiently small that the selective advantage of adaptablegenotypes after the environmental shift is O( √1N ) or smaller. The
less realistic case in which adaptable individuals are common isanalyzed in Supporting Information.
In our regime of interest, following the environmental shift,
the frequencies of nonadaptable genotypes behave nearly deterministically, and so their dynamics can be accurately approximated by ordinary differential equations (ODEs). In particular,
if we let X i denote the frequency of nonadaptable genotype iin the population, then we can approximate the dynamics by thefollowing system of ODEs:
⎞
⎛

d Xi
σ j X j ⎠ Xi .
= ⎝σi −
dt
j ∈A
/

1602

EVOLUTION MAY 2012

(2)

Here, A denotes the set of adaptable genotypes, and σi denotes
the selective coefficient of genotype i, which we define below.
The sumj ∈A/ σ j X j gives the mean fitness of genotypes thatare not adaptable. Loosely speaking, the nonadaptable types behave in a deterministic manner because we have assumed they arecommon at the time of the environmental shift. We formally justify this deterministic approximation for nonadaptable genotypes
in Supporting Information, and we demonstrate below that thisapproach provides a highly accurate approximation of our modelby comparing it to Monte Carlo simulations of the exact Moran
process.
Unlike the nonadaptable types, the adaptable types are assumed rare at the time of the environmental shift, and so we cannot use a deterministic approximation. Instead, we approximate
the behavior of the adaptable genotypes by deriving a diffusion
approximation of the Moran process (see Supporting Information). The resulting diffusion approximation can be described bya stochastic differential equation (SDE), which is similar to an
ODE except that the instantaneous dynamics involve Brownian
motion. In particular, we have derived the following SDE to approximate the (stochastic) dynamics of the frequencies of each
adaptable genotype, j:



σj −
σi X i Y j + θ
m k j X k dt + 2Y j dB j (t).
dY j =
i ∈A
/
k ∈A
/
(3)
The notation used to express equation (3) is the standard notationfor SDEs (Gardiner 2004), in which B j (t) denotes independent
Brownian motions. Intuitively, the SDE above indicates that thechange in frequency of genotype j, Y j , in a small interval oftime, t, is a Gaussian random variable with meanσj −

σi X i Y j + θ

i ∈A
/



m kj X kk ∈A/and standard deviation
2Y j .
Here, Y j denotes the number of individuals, rescaled by √1N , atadaptable genotype j (see Supporting Information). A denotes
the set of adaptable genotypes, θ = N μ denotes the populationscaled mutation rate, and σ j denotes selection coefficients defined below. The various terms in this SDE each have intuitive
interpretations. As is typical in population-genetic models, theinstantaneous mean change in the frequency of genotype j isdetermined by two terms, the first due to selection and the second due to mutation; whereas the instantaneous variance in thechange of genotype j’s frequency is determined by its current
frequency, which encodes the effect of genetic drift. The variance term arising from genetic drift has a somewhat unusual form
and it involves a 2Y j because we have scaled the population of

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

adaptable genotypes by √1N , instead of the more typical N1 (seeSupporting Information). We formally justify this diffusion approximation for adaptable genotypes in Supporting Information,
and we demonstrate below that this approach provides a highly
accurate approximation of our model by comparing it to Monte
Carlo simulations of the exact Moran process.
The selection terms in equations (2) and (3) arise from environmental perturbations, which cause selective differences among
genotypes on the neutral network. For a nonadaptable genotypei, such a perturbation will either have no effect (with probability φi ), or it will produce a strongly deleterious phenotype withprobability (1 − φi ). Thus, the expected fitness of a nonadaptablegenotype i is 1 − (1 − φi ), and the selection coefficient on suchindividuals is:σi√ = −(1 − φi ),
N

i∈
/ A.

(4)

Likewise, adaptable genotypes express a deleterious phenotype
with probability KK−1 (1 − φi ), and they express the newly advantageous target phenotype, which has fitness 1 + s, with probability
1(1 − φi ). Therefore, the selection coefficient for an adaptable
K
genotype is1+sσi(1 − φi ) − (1 − φi ), i ∈ A.
(5)
√ =
K
N
√
The factor 1/ N appears the two equations above because wehave assumed the selection coefficients are O √1N (see Supporting Information).
Finally, we note that, following the environmental shift,
adaptable genotypes
mutate to the newly adaptive phenotype at
1−qirate μ K . This flux into the newly adaptive phenotype, alongwith equations (2) and (3), can be used to solve analytically forthe distribution of waiting times for a newly adaptive genotype toarise in the population. In particular, we can compute the meantime it takes for a population to adapt to a novel environment, for
an arbitrary neutral network with arbitrary variation in environmental and genetic robustness across the neutral genotypes. The
full analytical solution to this problem is given in the Supporting
Information.

A Simplified Model
We now consider a simplified version of the model above, which
will allow us to develop intuition about the effects of environmental robustness on adaptability. We assume that all genotypes havethe same mutational robustness, q. Each genotype on the neutral
network is assigned either high or low environmental robustness,denoted by φ H and φ L . Following a neutral mutation, genotypes
with high environmental robustness produce other genotypes withhigh environmental robustness with probability π H (and there-

fore produce genotypes with low environmental robustness withprobability 1 − π H ), whereas genotypes with low environmentalrobustness produce other genotypes with low environmental robustness with probability π L (and therefore produce genotypes
with low environmental robustness with probability 1 − π L ). The
neutral network now consists of four distinct classes of genotypes:
high environmental robustness and adaptable, high environmentalrobustness and nonadaptable, low environmental robustness andadaptable, and low environmental robustness and nonadaptable.
The mutation rates among these four genotypes, and between
the adaptable and the newly adaptive genotypes, are summarized
in Figure 2. Henceforth, we analyze our model in terms of the
population density at these four classes of genotypes, rather thanin terms of the frequencies of each individual genotype. Eachclass represents a large number of distinct genotypes and therefore harbors much genetic diversity. As the rates of mutations
into a particular class increase, this corresponds to increasing thenumber of distinct genotypes within that class, and thus increases
the genetic diversity within the class. Figure 2 shows the rates ofmutations between and within classes.
In our simplified model, the mutational robustness, q, isthe same for all genotypes. Therefore, we henceforth refer togenotypes of high and low robustness, with the understanding that this means genotypes of high and low environmentalrobustness.
The parameters π H and π L provide a measure of how clustered the genotypes of high and low robustness are on the neutral network. In particular, an individual with a high-robustness
1
neutral mutations, on average, begenotype will undergo 1−πHfore it encounters a genotype of low robustness. Similarly, an1
individual with a low-robustness genotype will undergo 1−πLmutations, on average, before it encounters a genotype of highrobustness. Moreover, the overall fraction of genotypes in the
1−π L
neutral network that have high robustness is given by 2−(πL +π H )(Fig. 3).
In the case of this simplified model, we can use equations(3)–(5) to write down the SDEs for the diffusion process associated with the two adaptable types (those of high and low
environmental robustness):KdY H = σ H Y H + qθ (π H X H + (1 − π L )X L ) dt
P(6)√
+ 2Y H d B H (t)KdY L = σ L Y L + qθ ((1 − π H )X H + π L X L ) dt
P√
+ 2Y L d B L (t),(7)

where X H and X L denote the fraction of the population at highand low robustness nonadaptable genotypes, and Y H and Y L

EVOLUTION MAY 20121603

A L E X A N D E R J. S T E WA RT E T A L .

μq 1 − KP

μq(1 − πH)

μqπH 1 − K
PμqπH KP

High robustness

High robustness

Non-adaptableμq(1 − πL)

q
μ1−
K

Newly
adaptive

Adaptable

Low robustness

Low robustness

μqπL 1 − K
PKμqπL P

μ 1−q
K

μq K
P
Figure 2. Rates of mutation among the four classes of genotypes, in our simplified model. Genotypes are either adaptable withhigh environmental robustness (top right), adaptable with low environmental robustness (bottom right), nonadaptable with highenvironmental robustness (top left), or nonadaptable with low environmental robustness (bottom left). Mutations change genotypes ofhigh robustness (filled circles) to genotypes of low robustness (open circles) with probability q(1 − π H ), and from low (pink) to high (red)robustness with probability q(1 − π L ). Mutations change nonadaptable genotypes (left) to adaptable genotypes (right) with probability
q 1 − KP and adaptable to nonadaptable with probability q KP . Adaptable genotypes mutate to a genotype encoding the newly adaptive
phenotype with probability 1−q
K . Mutations within classes are also indicated.

denote the number of individuals, rescaled by √1N (see Supporting
Information), at high- and low-robustness adaptable genotypes.
As before, the deterministic part of the system consists of a selection term, and an influx term due to mutation from nonadaptable
to adaptable types. We note that when the adaptable types are
rare, the probability of back mutation is rare, and thus it does notappear in our limiting diffusions, as justified in the Supporting Information. Following the environmental shift, the selection terms
σ H and σ L , for the high- and low-robustness adaptable genotypes,
respectively, are given byσH
1+s(1 − φ H ) +  (φ H − φ L ) X L
√ =
K
N

(8)

1+s
σL
(1 − φ L ) −  (φ H − φ L ) X H .
√ =
K
N

(9)

These selection coefficients have an intuitive interpretation. The
first term in equation (8) and equation (9) represents the selective
advantage of adaptable genotypes over the rest of the population that results from their occasionally expressing the newly
adaptive phenotype. In particular, such individuals express an

1604

EVOLUTION MAY 2012

alternative phenotype with probability  (1 − φ), and that alternative phenotype will be adaptive, with selective benefit s, withprobability K1 . The second terms in equations (8) and (9) represent the strength of selection on genotypes of high and low
robustness in the population. These can be understood as follows:
high-robustness genotypes gain a selective advantage  (φ H − φ L )over low-robustness genotypes (because they express deleterious alternative phenotypes less frequently). This results in a selective advantage for high-robustness genotypes in equation (8)and a selective disadvantage for low-robustness genotypes inm
equation (9).
Since both terms in equation (8) favor high-robustness
adaptable genotypes, σ H is always positive. However, for lowrobustness genotypes, one term gives a selective advantage andthe other term a selective disadvantage; therefore, σ L can be either
positive or negative.
Finally, the influx terms in equation (6) and equation (7)depend on π H and π L —that is, on how genotypes of high andlow robustness are clustered across the network.
In the Results section, we will use the analytical equations
above to study how the various factors in our model influence theE N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

Figure 3. Illustration of mutations within and between classes
for different clusterings of individual genotypes. Left—genotypes
with high environmental robustness (filled circles) are tightly clus-

tered, and tend to neighbor one another rather than genotypes
with low environmental robustness (empty circles). This corresponds to π H ∼ 1. Right—genotypes with high environmentalrobustness tend not to neighbor other genotypes of the same
type. This corresponds to π H ∼ 0. When genotypes of the highrobustness type tend to neighbor each other, the population canspread to many different genotypes and the degree of genetic
diversity in the population increases.

capacity of a population to adapt. We will also compare our analytic approximation for the mean adaptation time to exact Monte
Carlo simulations of the Moran process. Simulations are carriedout with populations of N = 10,000 individuals. An individual ischosen to reproduce with a probability determined by its fitness
relative the population as a whole. The chosen individual produces
a viable offspring provided the offspring does not suffer a deleterious mutation or a deleterious environmental perturbation. Mutations occur at rate μ = 10−3 , across the entire genome. A mutation
will change the offspring’s genotype, phenotypic neighborhood,and environmental robustness as described above (see Fig. 2). The
offspring produced in such a reproduction event replaces an individual selected randomly from the population. Simulations are
stopped whenever an individual arises whose genotype encodes
the newly adaptive phenotype. The waiting time before such anadaptation arises is averaged over 10,000 replicate simulations, toproduce the mean waiting time plotted in the figures below.

Results
We focus on the case in which the population has evolved for along time in response to selection for a given phenotype, and hasreached its equilibrium distribution on the neutral network prior
to the environmental shift. Since we are interested in the influence of environmental robustness on adaptation, we keep fixed
the population-scaled mutation rate, θ = N μ, mutational robustness, q, number of phenotypic neighbors, K , and total numberof alternative phenotypes, P; and we vary the other parameters:
the degree of environmental robustness, φ H , φ L , the clustering of
high- and low-robustness genotypes, π H , π L , the selective advantage of the newly adaptive phenotype, s, and the probability thatan individual experiences environmental noise .

Prior to the environmental shift, approximately all nonadaptable individuals will have high robustness (i.e., X H = 1 and
X L = 0), since we have assumed high-robustness genotypes have√a selective advantage over low-robustness genotypes O( N ) (seeSupporting Information). Following the environmental shift, thedynamics of the nonadaptable individuals are approximated by
the ODE in equation (2). Substituting these values in equations (6)–(9), we see that the selection terms vary with the environmental robustness of the different genotypes, φ H and φ L ,with the rate of environmental perturbations, , and with thestrength of selection for the newly adaptive type, s. However,
the mutational influx of adaptable genotypes depends only on
how genotypes of high robustness are distributed in the neutral
network, π H .
Our central question—how do opposing influences of environmental robustness resolve to affect adaptation time?—is mostnaturally addressed by considering how adaptation time varies
with the probability that an individual experiences an environmental perturbation during development, . However, this relationship is complex, and its form depends on the distribution ofhigh-robustness genotypes, π H , the degree of environmental robustness of different genotypes, φ H and φ L , and the strength of
selection for the newly adaptive phenotype, s. Therefore, we begin by considering the relationship between these parameters andadaptation time, before considering the frequency of environmental perturbations.
DEPENDENCE OF ADAPTATION TIME ON THE
CLUSTERING OF ROBUST GENOTYPES

We begin by considering how the clustering of high-robustness
genotypes influences adaptation time. Varying π H from 0 to 1
will increase the clustering of high-robustness genotypes. Whenπ H ∼ 1, the neutral neighbors of high-robustness genotypes arealmost always other high-robustness genotypes; whereas whenπ H ∼ 0, the neutral neighbors of high-robustness genotypes arealmost always low-robustness genotypes (Fig. 3).
There are two distinct selective regimes that correspond toqualitatively different ways in which adaptation time varies with
< 1, adaptable genotypes with high robustness
π H . When 1+s
Kare selected for more strongly than adaptable genotypes with low
robustness (i.e., σ H > σ L ), and adaptation time decreases as π H
increases (Fig. 4 ). In contrast, when 1+s> 1, adaptable genotypes
K
with low robustness are selected for more strongly than adaptablegenotypes with high robustness (σ H < σ L ), and adaptation timeincreases as π H increases (Fig. 4).
< 1, is more realistic beThe first of these two regimes, 1+sK
cause we expect the size of the phenotypic neighborhood, K , tobe much greater than 1 and the selective advantage, s, of the
newly adaptive phenotype to be less than 1, typically. In themost realistic scenario, therefore, adaptation time will decrease

EVOLUTION MAY 20121605

A L E X A N D E R J. S T E WA RT E T A L .

namics in these differing regimes can be illustrated through timeseries plots of the frequency of high- and low-robustness adaptable genotypes in the different selective regimes (see Supporting
Information).
Here, we have focused on the distribution of high-robustness
genotypes, π H , assuming that the clustering of low-robustness
genotypes, π L , remains fixed. Varying the clustering of lowrobustness genotypes does not affect adaptation time if, as wehave assumed, the population is at equilibrium prior to the environmental shift (eqs. 6 and 7).
DEPENDENCE OF ADAPTATION TIME ON THE
STRENGTH OF ENVIRONMENTAL ROBUSTNESS
Figure 4. Mean adaptation time as a function of the clusteringof high-robustness genotypes, π H . As π H increases, the probability that a high-robustness genotype neighbors another highrobustness genotype increases (i.e., more clustering). When 1+s
K <1 (red line), this causes a decrease in the mean adaptation time ofthe population. When 1+s
K > 1 (blue line), this causes a increase inmean adaptation time. Plots show populations of N = 10,000 individuals, with μ = 0.001,  = 0.1, P = 100, K = 5, q = 0.5, φ H = 0.9,
and φ L = 0.1 with s = 1 (red line) and s = 5 (blue line). Lines indicate the analytic solution to our model, whereas dots indicate themeans of 10,000 replicate Monte Carlo simulations.

with increasing π H —that is, it takes longer for the population
to adapt when high-robustness genotypes are not clustered. The
intuition underlying this result is simple. Prior to the environmental shift, high-robustness genotypes are always selectively
favored (because environmental perturbations are always deleterious). At the same time, after the environmental shift, highrobustness adaptable genotypes are more selectively favored than
low-robustness adaptable genotypes, that is, σ H > σ L in equations 8–9. The more tightly clustered high-robustness genotypes
are (i.e., the larger π H ), the more often neutral mutations to genotypes of high robustness will produce other genotypes of highrobustness. As a result, there is more standing neutral genetic
variation in the population prior to the environmental shift, because there are a greater number of distinct genotypes within
the class of high environmental robustness. This in turn causes
adaptable genotypes to be produced at a higher rate after the envi> 1, lowronmental shift, speeding adaptation. In contrast, if 1+s
Krobustness adaptable genotypes are selectively favored over highrobustness adaptable genotypes, that is, σ H < σ L in equations(8)–(9). Therefore, increasing the clustering of high-robustness
genotypes increases the neutral genetic variation in the population prior to the environmental shift, but decreases the rate at
which advantageous adaptable genotypes are produced after theenvironmental shift. These two effects oppose each other, and wefind that the overall effect is to slow adaptation. The adaptive dy-

1606

EVOLUTION MAY 2012

We now consider how variation in the strength of environmental robustness (parameters φ H and φ L ) affects adaptation time.
As Figure 5 shows adaptation time always increases with increasing φ H (i.e., when the robustness of the high-robust type isincreased, adaptation is slower). There are two distinct reasons for
this, associated with the selection terms in equation (8) and equation (9). First, as φ H increases, the (postshift) selection term forhigh-robustness adaptable genotypes decreases (this can be seen
directly from eq. 8, setting X H = 1 and X L = 0, but it is also
intuitive: as φ H increases, the chance that a high-robust genotype
expresses the adaptive phenotype decreases), so that the population is less likely to flow toward these adaptable genotypes,
which slows adaptation. Second, as φ H increases, the selective
term for low-robustness adaptable genotypes also decreases (thiscan be seen directly from eq. 9, with X H = 1 and X L = 0 but
it is also intuitive: as φ H increases, the selective disadvantage ofa low-robustness genotype compared to the rest of the population becomes greater), which again retards flow toward adaptablegenotypes and slows adaptation.
Varying the environmental robustness of low-robustness
genotypes, φ L , has a more complex effect on adaptation time.
Increasing φ L has no effect on the selective advantage of highrobustness adaptable genotypes. However, it has conflicting
effects on the selective advantage of low-robustness adaptablegenotypes (as seen in eq. 9). On the one hand, increasing φ L decreases the probability that an environmental perturbation on a
low-robustness adaptable genotype will produce the newly adaptive phenotype. However, increasing φ L also reduces the probability that an environmental perturbation on a low-robustness adaptable genotype will produce a deleterious phenotype. Whether
the overall selection term σ L is increasing or decreasing with
φ L therefore depends on how advantageous it is to produce the. Whennewly adaptive phenotype. This depends on the ratio 1+s
K
1+s
<
1,
the
selection
term
σ
for
the
low-robustness
adaptableLK
genotypes increases as φ L increases (eq. 9); as a result, adapta> 1, the selection
tion time decreases as φ L increases. When 1+s
Kterm for the low-robustness adaptable genotypes decreases as φ L

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

mentally robust genotypes would always prefer to reduce theirrobustness, as we saw above. Thus, populations with less variation in environmental robustness will generally adapt faster thanthose with greater variation in environmental robustness.
DEPENDENCE OF ADAPTATION TIME ON THE
FREQUENCY OF ENVIRONMENTAL PERTURBATION, 

Figure 5. Top—mean adaptation time as a function of the amountof environmental robustness, φ L , among low-robust genotypes.
When selection for the newly adaptive phenotype is relatively
weak ( 1+sK < 1, red line), increasing φ L reduce the frequencyof adaptable individuals, which retards adaptation. When selection is stronger ( 1+s
K > 1, blue line), increasing φ L promotes
adaptable individuals and thus facilitates adaptation. Plots show
populations of N = 10, 000 individuals, with μ = 0.001,  = 0.1,
P = 100, K = 5, q = 0.5, φ H = 1, and π H = 0.5 with s = 2 (redline) and s = 6 (blue line). Bottom—mean adaptation time with
φ H . As φ H increases, the selection term for both high- and lowrobustness adaptable genotypes decreases. Plots show populations of N = 10,000 individuals, with μ = 0.001,  = 0.1, P = 100,
K = 5, q = 0.5, φ L = 0, and π H = 0.5 with s = 2 (red line) ands = 5 (blue line). Lines indicate the analytic solution to our model,
whereas dots indicate the means of 10,000 replicate Monte Carlosimulations.

increases (eq. 9); as a result, adaptation time increases as φ L increases (Fig. 5). As in the previous section, the more realistic caseoccurs when 1+s< 1. Therefore, we expect that increasing φ L
Kwill generally decrease adaptation time.
In summary, from the perspective of adaptation, the least environmentally robust genotypes in a population would generally
prefer to increase their robustness; whereas the most environ-

We now turn to the central question of this study—how does
adaptation time vary with the frequency of environmental perturbations during development? Increasing the frequency of suchperturbations, , will increase the magnitude of the selection terms
associated with both high- and low-robustness adaptable genotypes (eqs. 8 and 9). However, increasing  does not change therelative strength of these selection terms, the direction of selection, or the total rate of influx into the adaptable genotypes. These
depend on the factors discussed in the previous two sections—the distribution of high-robustness genotypes, and the strengthof environmental robustness associated with both high- and lowrobustness genotypes.
There are two selective regimes of interest when considering
how adaptation time depends on the frequency of environmental perturbations. The first regime occurs when low-robustness
adaptable genotypes have lower fitness than the population as awhole, that is, when σ L is negative. According to equation (9),H −φ L
this occurs when 1+s< φ1−φ. The second regime occurs whenKL
low-robustness adaptable genotypes are selectively favored compared to the population as a whole, that is, when σ L is positive:
1+s
H −φ L> φ1−φ. High-robustness adaptable genotypes are always
KL
selectively favored compared to the population as a whole, thatis, σ H > 0 (eq. 8).
In the first selective regime, increasing  increases thestrength of selection for high-robustness adaptable genotypes,
while increasing the strength of selection against low-robustness
adaptable genotypes. As a result increasing  has conflicting effects on the adaptation time—it increases the adaptation time fromlow-robustness adaptable genotypes but decreases the adaptation
time from high-robustness adaptable genotypes. Initially, as  increases from 0, adaptation time increases, as selection against
low-robustness adaptable genotypes reduces the density of individuals at all adaptable genotypes. However, as  continues toincrease, selection for high-robustness adaptable genotypes concentrates the population at these genotypes further, and causes theoverall adaptation time to decrease. This results in a nonmonotonic relationship between adaptation time and : increasing thefrequency of environmental perturbation can either facilitate or
impede adaptability (Fig. 6 ). In this regime, the population adapts
most quickly when environmental noise is either very frequent orvery infrequent.
In the second selective regime, increasing  increases thestrength of selection for both high- and low-robustness adaptable

EVOLUTION MAY 20121607

A L E X A N D E R J. S T E WA RT E T A L .

the relationship between mutational robustness and adaptation(see Supporting Information).

Discussion

Figure 6. Mean adaptation time as a function of the probabilityof environmental perturbation, . When selection for the newly
φ H −φ L
adaptive phenotype is relatively strong ( 1+s
K > 1−φ L , blue line),
increasing amounts of environmental perturbation promotes all
types of adaptable genotypes, facilitating adaptation. When se-

φ H −φ L
lection is weaker ( 1+sK < 1−φ L , red line) increasing  has contrary effects on the low-robust and high-robust adaptable individuals; in this regime, the relationship between mean adapta-

tion time and probability of environmental perturbation is nonmonotonic (increasing and then decreasing). Plots show populations of N = 10,000 individuals, with μ = 0.001, P = 100, K = 5,
q = 0.5, φ H = 0.9, φ L = 0.1, and π H = 0.1 with s = 2 (red line) ands = 5 (blue line). Lines indicate the analytic solution to our model,
whereas dots indicate the means of 10,000 replicate Monte Carlosimulations.

genotypes, and it therefore always decreases adaptation time. In
this regime, more environmental noise always speeds adaptation.
In summary, when environmental fluctuations are rare( ∼ 0), increasing their frequency can either impede or facilitate
adaptation depending on the strength of selection for the newly
adaptive phenotype. However, when environmental fluctuations
are common ( ∼ 1), increasing their frequency will always speed
adaptation.
We have also studied how the clustering of high-robustness
genotypes influences the relationship between environmental perturbations and mean adaptation time (see Supporting Information).
DEPENDENCE OF ADAPTATION TIME ON
MUTATIONAL ROBUSTNESSWe have also studied how environmental noise can influence therelationship between mutational robustness and adaptability. Previous work has shown that, in the absence of environmental noise,intermediate levels of mutational robustness produce the fastest
rates of adaptation (Wagner 2008; Draghi et al. 2010). We find thatthe addition of environmental noise does not qualitatively change

1608

EVOLUTION MAY 2012

The expression of a phenotype is inherently noisy. As a result,when a population experiences purifying selection, those genotypes that express the optimal phenotype more reliably will beselectively favored. This is the key insight behind Waddington’sconcept of environmental canalization, developed over half a century ago (Waddington 1942, 1953a,b, 1959). Nevertheless, theextent to which environmental robustness evolves in response tonoise, and how this phenomenon impacts evolution over longer
time scales, as the optimal phenotype itself changes, is not fully
understood.
When canalization fails, environmental noise elicits the expression of an alternative phenotype in an individual. The resulting phenotype is often then a phenocopy—that is, a phenotype
that mimics the effect of a mutation, but is not heritable. When apopulation is faced with a new environment, the phenocopy phenomenon may speed adaptation, because it allows some individuals to express adaptive phenotypes without waiting for mutations.
However, when a population is adapted to the current environment, environmental noise tends to concentrate the population on
the most robust genotypes, which will reduce genetic diversity
and the capacity for adaptation. We have sought to elucidate howthese two conflicting effects of environmental robustness resolve
in determining the adaptability of populations.
We have employed a general population-genetic model thatwas previously used to study the interplay between mutational
robustness and adaptability (Draghi et al. 2010). Our model specifies a fitness landscape in terms of the properties of a neutral
network of genotypes. Whenever an individual develops an alternative phenotype as a result of environmental perturbation, it
develops into one of the K phenotypes accessible by mutation;
this assumption encodes the concept of phenocopy in our modeling framework.
We have studied the relationship between environmental
noise and adaptation time—that is, the time required for a population, initially in equilibrium, to acquire a new target phenotype following an environmental shift. In our simplest model, theadaptation time depends on several parameters: the frequency ofenvironmental perturbations (), the selective advantage of thenovel target phenotype (s), the levels of environmental robustness
among genotypes on the neutral network (φ H and φ L ), and theclustering of genotypes with high environmental robustness (π H ).
The interplay of these parameters gives rise to a surprisingly rich
set of behaviors, which are summarized in Table 2.
The relationship between environmental and mutational robustness may be further complicated if both are allowed to vary

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

Table 2.

Summary of results.

Parameter varied

Selective regime

Mean adaptation timeHigh-robustness genotype clustering, π H
High-robustness genotype clustering, π H
Low-robustness genotype clustering, π L
High environmental robustness level, φ H
Low environmental robustness level, φ L
Low environmental robustness level, φ L
Rate of environmental perturbation, Rate of environmental perturbation, Mutational robustness, q

1+s
<1
K
1+s
≥1
K

Decreasing
Increasing
Constant
Increasing
Decreasing
Increasing
Nonmonotonic (increasing then decreasing)
Decreasing
Nonmonotonic (decreasing then increasing)

All
All
1+s
<1
K
1+s
≥1
K
1+s
H −φ L< φ1−φKL
φ H −φ L
1+s
≥
K
1−φ L

All

across the neutral network simultaneously. Our analysis has assumed the simplest case—when environmental and mutational
robustness vary independently. Even in this simple case, we havefound a rich diversity of possible behaviors (Table 2). The assumption of independence is partly supported by the few systematic empirical studies of mutational and environmental robustness. Such studies have found that polymorphisms affecting
mutational robustness map to different loci than those influencing
environmental robustness (Proulx et al. 2007; Fraser and Schadt
2010; Lehner 2010), suggesting that these traits may be largely
independent. Nonetheless, the impact of covariation in these two
types of robustness (Hermisson and Wagner 2004; Lehner 2010;
Meiklejohn and Hartl 2002; Proulx et al. 2007; Wagner et al.
1997) remains an important topic for future research.
It is interesting to contrast our results on environmental noisewith quantitative studies of the Baldwin effect, which also describe the genetic reinforcement of initially nonhereditary traits(Ancel 1999, 2000). Whereas we have focused here on the tension between costs and benefits of randomly expressing an alternative phenotype during development, studies of the Baldwin
effect have been concerned with the tension between the adaptive
advantage and inherent cost of phenotypic plasticity—that is, theability of an individual to express a range of phenotypes and to
deterministically choose among these the most fit phenotype in
the current environment. In other words, in our study, the phenotype expressed by an individual is subject to random perturbation,
whereas in studies of the Baldwin effect, the expressed phenotype is chosen deterministically: if an individual has the ability
to express the optimal phenotype with its range of alternatives,
it will. Models of the Baldwin effect have uncovered a complex
interplay between the rate at which the environment changes, thedegree of phenotypic plasticity in a population, and the rate ofadaptation (Ancel 1999, 2000). We have studied a different phenomenon, but likewise found a complex relationship between theamount of development noise and the rate of adaptation. Bothphenomena (phenotypic plasticity and environmental robustness)
highlight the complex relationship between environmentally in-duced phenotypic variation and the capacity for a population toadapt.
A phenomenon analogous to the Baldwin effect may occur
in our system if the newly adapted individual that first fixes in
the population has low robustness (analogous to high plasticity
in the models of [Ancel 1999, 2000]); whereas subsequent stabilizing selection will then select for a high-robustness (adapted)
individual. This situation is more likely to occur in the regime
σL > σ H .
In deriving our analytical results, we have relied on threestrong assumptions: a neutral mutation completely redraws an individual’s phenotypic neighborhood; the number of phenotypes,
K , in a genotype’s neighborhood does not vary across the neutral
network; and alternative phenotypes are generally lethal. Relaxing
each of these assumptions does not change our qualitative results(see Supporting Information). In particular, we relax the first assumption by introducing a parameter, f , which is the fraction of K
neighbors that are redrawn following a neutral mutation. Allowing
such correlations between the phenotypic neighborhoods of neutral neighbors affects the rate at which individuals arrive at adaptable genotypes through mutation, but it preserves the qualitative
relationships between adaptability and environmental robustnessreported above. Our qualitative results are also unchanged if thesize of the phenotypic neighborhood, K , varies with genotype.
In this case, however, the selective regimes described in Table 2. Finally, our qualitative
now depend on the average value of 1+s
K
results are unaltered if alternative phenotypes are moderately deleterious rather than lethal; however, the different selective regimes
described in Table 2 would then depend on s, K , and the sizeof the fitness penalty suffered by nonadaptive alternative phenotypes (see Supporting Information). We have limited our analysis
to the regime in which the newly adaptive type is always favored
over the initially fit phenotype, regardless of its environmentalrobustness. However, it is possible to imagine scenarios in whichs is sufficiently small that newly adaptive genotypes with low
robustness are deleterious compared to high-robustness individuals of the initially fit phenotype. This situation will occur when s

EVOLUTION MAY 20121609

A L E X A N D E R J. S T E WA RT E T A L .

is O() or smaller (see Supporting Information). To study adaptability in this case, we must determine the waiting time until
a newly adaptive individual with high robustness is produced.
However, the qualitative relationships we observe in this caseare the same as those described in the main text (see Supporting
Information).
In this article, we have considered adaptation to a single environmental shift. It would also be interesting to study repeatedenvironmental shifts. To do so will require analyzing the adaptation time of a population which is not at equilibrium on its neutral
network when the next environmental shift occurs. Accounting
for such repeated transient dynamics should be possible in ourmodeling framework, but it will certainly be nontrivial. The relationship between the rate of environmental shifts and the evolutionof the environmental robustness thus remains an interesting areafor future research.
Our modeling framework has allowed for substantial variation in robustness among a large set of genotypes. This framework
is appropriate if many sites across an organism’s genome can independently influence its robustness, such as might be expected
if networks of interacting genes jointly determine an organism’s
phenotype. In fact, the widespread contribution of many genomic
sites to robustness has been directly observed in organisms ranging from yeast to plants and vertebrates (Fraser and Schadt 2010).
The framework for our study is therefore rather different frommodels restricted to a specific locus (such as the HSP90 locus or a
DNA repair locus) that alone determines variation in robustness.
Despite longstanding interest in phenotypic robustness, therole of robustness in adaptation remains a topic of active investigation. In contrast to most theoretical work, which has focused on
how mutational robustness influences adaptation (de Visser et al.
2003; Lenski et al. 2006; Wagner 2008; Draghi et al. 2010), thisstudy provides a framework for analyzing the effects of environmental robustness. Although we have uncovered a complex set of
possible behaviors (summarized in Table 2), one important general observation emerges: environmental perturbations can either
facilitate or impede adaptation. These results have been derived
in a general theoretical framework, and so they should apply toabroad range of biological situations and systems.
Theoretical studies of environmental robustness must eventually be reconciled with empirical data, either from surveys ofstanding variation or, ideally, from evolution experiments. Fortunately, there is an extensive literature based on quantitative genetics that describes empirical patterns of environmental robustness in organisms, including yeast, plants, and mammals (Mackay
and Lyman 2005; Hall et al. 2007; Nogami et al. 2007; Sangster
et al. 2008; Fraser and Schadt 2010). Such studies typically finda broad array of traits (morphological, life-history, and gene expression levels) that are buffered against environmental noise and,furthermore, extensive standing genetic variation in the degree of1610

EVOLUTION MAY 2012

such buffering among individuals. Thus, environmental robustness is apparently selected for in many natural populations, inaccordance with theoretical expectations for a population understabilizing selection for a target phenotype (Wagner et al. 1997;
Hermisson and Wagner 2004); but a natural population typicallyalso contains variation in the degree of environmental robustness
among genotypes. Our analysis suggests that the existence of robust variants in a population, and also the variation in robustness
among individuals, will have a significant impact on the population’s capacity to adapt to a novel environment. Thus, the conceptof canalization, first considered by Waddington for populations
under stabilizing selection, plays an important role in shaping theevolution of populations as they adapt to new environments overlonger timescales.
ACKNOWLEDGMENTS
We thank members of the Plotkin lab for constructive comments and input.
JBP acknowledges support from the James S. McDonnell Foundation, theAlfred P. Sloan Foundation, the David and Lucille Packard Foundation,
the Burroughs Wellcome Fund, Defense Advanced Research Projects
Agency (HR0011-05-1-0057), and the US National Institute of Allergy
and Infectious Diseases (2U54AI057168).

LITERATURE CITED
Ancel, L. W. 1999. A quantitative model of the Simpson-Baldwin effect. J.
Theor. Biol. 196:197–209.
———. 2000. Undermining the Baldwin expediting effect: does phenotypic
plasticity accelerate evolution? Theor. Popul. Biol. 58: 307–319.
Ancel, L. W., and W. Fontana. 2000. Plasticity, evolvability, and modularity
in rna. J. Exp. Zool. 288: 242–283.
Bloom, J., S. Labthavikul, C. Otey, and F. Arnold. 2006. Protein stability
promotes evolvability. PNAS 103: 5869–5874.
Chow, K., and K. Chan. 1999. Stress-induced phenocopy of C. elegans defines functional steps of sensory organ differentiation. Develop. Growth.
Differ.41: 629–637.
Draghi, J., T. Parsons, G. Wagner, and J. Plotkin. 2010. Mutational robustness
can facilitate adaptation. Nature 436: 353–355.
Earl, D., and M. Deem. 2004. Evolvability is a selectable trait. PNAS 101:
11531–11536.
Ewens, W. J. 2004. Mathematical population genetics, vol. v. 27.
2nd ed. Springer, New York. Available at http://www.loc.gov/
catdir/enhancements/fy0818/2003065728-d.html.
Fraser, H., and E. Schadt. 2010. The quantitative genetics of phenotypicrobustness. PLoS ONE 5: e8635.
Gardiner, C. W. 2004. Handbook of stochastic methods for physics, chemistry,
and the natural sciences. 3rd ed. Springer, Berlin. Available at http://
www.loc.gov/catdir/enhancements/fy0818/2004043676-d.html.
Gibson, G., and I. Dworkin. 2004. Uncovering cryptic genetic variation. Nat.
Rev. Genet. 5: 681–690.
Goldschmidt, R. 1935. Gen und aussencharakter. Biologisches Zentralblatt
55: 535–554.
Goldsmith, M., and D. Tawkik. 2009. Potential role of phenotypic mutationsin the evolution of protein expression and stability. PNAS 106: 6197–
6202.
Hall, M. C., I. Dworkin, M. C. Ungerer, and M. Purugganan. 2007. Genetics
of microenvironmental canalization in Arabidopsis thaliana. Proc. Natl.
Acad. Sci. USA 104: 13717–13722.

E N V I RO N M E N TA L RO B U S T N E S S A N D T H E A DA P TA B I L I T Y O F P O P U L AT I O N S

Hermisson, J., and G. P. Wagner. 2004. The population genetic theory of
hidden variation and genetic robustness. Genetics 168: 2271–2284.
Jarosz, D. F., and S. Lindquist. 2010. Hsp90 and environmental stress transform the adaptive value of natural genetic variation. Science 330: 1820–
1824.
Lehner, B. 2010. Genes confer similar robustness to environmental, stochastic,
and genetic perturbations in yeast. PLoS ONE 5: e9035.
Lenski, R., J. Barrick, and C. Ofria. 2006. Balancing robustness and evolvability. PLoS Biol. 4: e428.
Levy, S. F., and M. L. Siegal. 2008. Network hubs buffer environmental
variation in Saccharomyces cerevisiae. PLoS Biol. 6: e264.
Lopez-Maury, L., S. Marguerat, and J. Bahler. 2008. Tuning gene expression tochanging environments: from rapid responses to evolutionary adaptation.
Nat. Rev. Genet. 9: 583–593.
Mackay, T. F., and R. F. Lyman. 2005. Drosophila bristles and the natureof quantitative genetic variation. Philos. Trans. R. Soc. Lond. B 360:
1513–1527.
Masel, J. 2005. Evolutionary capacitance may be favored by natural selection.
Genetics 170: 1359–1371.
Masel, J., and A. Bergman. 2003. The evolution of the evolvability propertiesof the yeast prion [psi+]. Evolution 57: 1498–1512.
Masel, J., and M. Trotter. 2010. Robustness and evolvability. Trends Genet.
26: 406–414.
McLoughlin, S., and S. Copley. 2008. A compromise required by gene sharing
enables survival: implications for evolution of new enzyme activities.
PNAS 105: 13497–13502.
Meiklejohn, C. D., and D. L. Hartl. 2002. A single mode of canalization.
Trends Ecol. Evol. 17: 468–473.
Meyerovich, M., G. Mamou, and S. Ben-Yehuda. 2010. Visualizing higherror levels during gene expression in living bacterial cells. PNAS 107:
11543–11548.
Meyers, L., F. Ance, and M. Lachmann. 2005. Evolution of genetic potential.
PLoS Comput. Biol. 1: e32.
van Nimwegen, E., J. P. Crutchfield, and M. Huynen. 1999. Neutral evolutionof mutational robustness. PNAS 96: 9716–9720.
Nogami, S., Y. Ohya, and G. Yvert. 2007. Genetic complexity and quantitative
trait loci mapping of yeast morphological traits. PLoS Genet. 3: e31.
Proulx, S., S. Nuzhdin, and D. Promislow. 2007. Direct selection on genetic
robustness revealed in the yeast transcriptome. PLoS ONE 2: e911.

Raj, A., and A. van Oudenaarden. 2008. Nature, nurture, or chance: stochastic
gene expression and its consequences. Cell 135: 216–226.
Rutherford, S. 2003. Between genotype and phenotype: protein chaperones
and evolvability. Nat. Rev. Genet. 4: 263–274.
Rutherford, S., and S. Lindquist. 1998. Hsp90 as a capacitor for morphologicalevolution. Nature 396: 336–342.
Sangster, T. A., S. Lindquist, and C. Queitsch. 2004. Under cover: causes,
effects and implications of hsp90-mediated genetic capacitance. BioEssays 26: 348–362.
Sangster, T. A., N. Salathia, S. Undurraga, R. Milo, K. Schellenberg,
S. Lindquist, and C. Queitsch. 2008. Hsp90 affects the expression ofgenetic variation and developmental stability in quantitative traits. Proc.
Natl. Acad. Sci. USA 105: 2963–2968.
Siegal, M. L., and A. Bergman. 2002. Waddington’s canalization revisited: developmental stability and evolution. PNAS 99: 10528–
10532.
Tawkik, D. 2010. Messy biology and the origins of evolutionary innovations.
Nat. Chem. Biol. 6: 692–696.
de Visser, J. A., J. Hermisson, G. P. Wagner, L. Ancel Meyers, H. BagheriChaichian, J. L. Blanchard, L. Chao, J. M. Cheverud, S. F. Elena, W.
Fontana, et al. 2003. Perspective: Evolution and detection of genetic
robustness. Evolution 57: 1959–1972.
Waddington, C. 1942. Canalization of development and the inheritance ofacquired characteristics. Nature 150: 563–565.
Waddington, C. 1953a. The Baldwin effect, genetic assimilation and homeostasis. Evolution 7: 386–387.
Waddington, C. 1953b. Genetic assimilation of an acquired character. Evolution 7: 118–126.
———. 1959 Evolutionary systems—animal and human. Nature 183:1634–
1638.
Wagner, A. 2008. Robustness and evolvability: a paradox resolved. Proc. R.
Soc. Lond. B 275: 91–100.
Wagner, A., and J. Wright. 2007. Alternative routes and mutational
robustness in complex regulatory networks. Biosystems 88: 163–
172.
Wagner, G. P., G. Booth, and H. Bagheri-Chaichian. 1997. A population
genetic theory of canalization. Evolution 51:329–347.

Associate Editor: Max Reuter

EVOLUTION MAY 2012

1611

A L E X A N D E R J. S T E WA RT E T A L .

Supporting Information
The following supporting information is available for this article:
Figure S1. Correlated phenotypic neighbourhood—Mean adaptation time as a function of the probability of environmental
perturbation, .
Figure S2. Variation in size of phenotypic neighbourhood—plots show populations of N = 10,000 individuals, with μ = 0.001,
P = 100, K̄ = 5, q = 0.5, φH = 0.9, φL = 0.1 and πH = 0.5 with s = 2 (red line) and s = 8 (blue line).
Figure S3. Nonlethal deleterious mutations—plots show populations of N = 10,000 individuals, with μ = 0.001, P = 100, K =10, q = 0.2, φH = 0.9, φL = 0.1 and πH = 0.5, d = 0.2 with s = 0.5 (red line) and s = 3 (blue line).
Figure S4. Deleterious low robustness adaptive genotype—plots show populations of N = 10,000 individuals, with μ = 0.001,
P = 100, K = 5, q = 0.5 and s = 0.05.
H −φ L< φ1−φ.
Figure S5. Mean adaptation time as a function of the probability of environmental perturbation, , in the regime 1+s
K
L
Figure S6.The mean adaptation time as a function of mutational robustness, q.
Figure S7. Most realistic selective regime with s/K = 0.2, πH = 0.5 and πL = 0.5.
Figure S8. Less realistic selective regime with s/K = 1.2, πH = 0.5 and πL = 0.5.
Figure S9. Most realistic selective regime with s/K = 0.2, πH = 0.5 and πL = 0.1.
Table S1. Model parameters.
Supporting Information may be found in the online version of this article.
Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting information supplied by theauthors. Any queries (other than missing material) should be directed to the corresponding author for the article.

1612

EVOLUTION MAY 2012

Environmental robustness and the adaptability of
populations: Supporting Information

Alexander J. Stewart12 , Todd L. Parsons1 , Joshua B. Plotkin1

November 9, 2011

1 Department of Biology, University of Pennsylvania, Philadelphia, PA, USA
2 CoMPLEX, University College London, Physics Building, Gower Street, London WC1E 6BT,

UK

1Analytical results

In this section we discuss in detail how adaptation time is related to environmental robustness
by considering the full population dynamics and distribution of adaptation times for thegeneral model presented in the main text.

1.1

Model Description

We focus on the general model, described in the main text, consisting of a population of N
individuals evolving under the infinite alleles Moran model. Genotypes which belong to theneutral network, are indexed by a set I of which a subset A ⊂ I are adaptable, meaning thatthey are a single mutation away from a genotype which has a selective advantage once anenvironmental shift occurs. We assume that at the time of the environmental shift individuals
with an adaptable genotype are rare, meaning that the number of adaptable individuals is�√ �ON or smaller. The less realistic case in which adaptable individuals are common isanalyzed in Supplementary Information.

Following Draghi et al. (2010) we adopt the concept of phenotypic neighbourhoods. A
mutation at genotype i is neutral with probability qi and results in a genotype which codes fora new phenotype with probability 1−qi . When a mutation is non-neutral, the new phenotype
which results is drawn from K alternative phenotypes, which constitute the phenotypicneighbourhood of the original genotype. The K phenotypes of the phenotypic neighbourhood
are drawn uniformly form a total P possible alternatives. All alternative phenotypes are
assumed to be lethal prior to the environmental shift. Following the environmental shift
we assume that an alternative phenotype gains a selective advantage, and we calculate theexpected time for the population to produce an individual with such a genotype.
1

We also extend the concept of phenotypic neighbourhoods to describe the eﬀects ofenvironmental perturbations on phenotypes. Environmental perturbations occur with rate �.
With probability φi an organism which experiences an environmental perturbation develops
the correct (target) phenotype, i. With probability 1 − φi it develops the phenotype of oneof it’s K phenotypic neighbours. The parameter φi therefore measures the environmentalrobustness of genotype i.
The per-genome mutation rate is µ and the genome length is L. We assume that all sites
mutate at the same rate so that mutations within the neutral network are described by theadjacency matrix M which has entries mij = 1 if genotype i can mutate to genotype j, andmij = 0 otherwise, with the mutational robustness qi of a genotype given by Eq. 1 in themain text. In this section for notational convenience we define m̃ij = L1 mij .
In order to determine the adaptability of a population of the type described above westudy the continuous time Moran model in which all genotypes have the same death rate
and birth rates vary according to the fitness of the genotype. Therefore the birth rate of
an individual with genotype i is 1 + Si . If the number of individuals with genotype i is Zi ,the population dynamics of the Moran model are described by a continuous-time Markov
process with transitionsrij (Z(t))

rij (Z(t))

Zi (t) −−−−−−→ Zi (t + ∆t) = Zi (t) + 1,

Zj (t) −−−−−−→ Zj (t + ∆t) = Zj (t) − 1.

that occur at rate rij (Z(t)), fordef

rij (z) =�(1 − µ)zi (1 + Si ) + µ�k∈Im̃ki zk (1 + Sk )

�

�z �
j

N

.

The first term describes the rate at which individuals with genotype i are born. The second2term describes the rate at which individuals of type j die. Individuals with genotype i areborn at a rate that depends on zi as well as on the strength of selection for that genotype, Si ,and the rate of mutation into and out of that genotype. Once an individual is born it replaces
an existing member of the population. The individual that is replaced is selected uniformly
at random from the existing individuals, and as a result the death rate of individuals withz

genotype j is simply Nj .
We note that alternative formulations of the Moran model, for example in which death
rates vary, or in the discrete time version in which the birth rates are normalised, will resultin the same limiting diﬀusions as those described below up to a rescaling of time (see e.g.
(Etheridge, 2011)). For the sake of improved clarity in what follows we will work withtransition rates which are described in terms of θ = N µ and a rescaled selective coeﬃcient√σi = N Si , such that we haverij (z) =

��

θ
1−
N

�

zi

�

σi
1+ √
N

�

�
�� � �θ �σkzj+m̃ki zk 1 + √
N k∈I
N
N

The notations adopted in what follows are summarised in Table 3 below:[Table 1 about here.]

In order to recover the full population dynamics of the model we consider the infinitesimal
generator A of the continuous time stochastic process Z(t), which is given by (Ethier andKurtz, 2005)(Af )(z) =�i,jrij (z) (f (z + ei − ej ) − f (z)) .

where ei is the unit vector in the direction of i.
3

1.2

Limit When Adaptable Types are RarePrevious work (Draghi et al., 2010) showed that a closely related model to the one discussed�√ �here has a diﬀusive limit when the number of adaptable individuals is O
N , while the√number of non-adaptable individuals is O (N ), and time is rescaled by N . To that end, we
consider the rescaled population Z̃i at each genotype i, where adaptable types are rescaled by√a factor N and non-adaptable types are rescaled by a factor N . Since the rescaling diﬀers
between adaptable and non-adaptable genotypes it is convenient to introduce the parameter
hi which tells us how a given genotype i is rescaled, i.e

Z̃i (t) = hi Zi (t),
where
hi =




 √1

if i ∈ A

N



1

otherwise.

N

We can now find the generator for the limiting process. Following (Ethier and Kurtz, 2005),
this is given bylim ANN →∞for(AN f )(z̃) =�√i,jN r̃ij (z̃) (f (z̃ + hi ei − hj ej ) − f (z̃)) ,√��where r̃ij (z̃) = rij (h−1
N increases the rate of eventsi z̃i )i∈I . Intuitively, rescaling time byby the same factor, while increasing or decreasing Zi corresponds to changing Z̃i by ±hi .

Therefore, since hi is small in the large N limit, to find the limiting generator we Taylor

4expand f (z̃ + hi ei − hj ej ) about z̃ to obtain(AN f )(z̃) =

√

N

�

�

hi

i∈I

�
j∈I

(r̃ij (z̃) − r̃ji (z̃))

∂f
1� 2�
∂ 2f
+
hi
(r̃ij (z̃) + r̃ji (z̃)) 2∂ z̃i 2 i∈I
∂ z̃i
j∈I

�
�
�
3
∂ 2f
−
hi hj r̃ij (z̃)+ O N−2 .
∂z̃
i ∂ z̃ji�=j�

In what follows, we derive the advective terms; the other terms are derived in a similar
fashion.
To begin, we haverij (z) − rji (z) =and, recalling
�

�

3

N2

�

θ
1−
N

�

�
�
θ �
σk
(σi − σj )zi zj + 2(m̃ki zj − m̃kj zi )zk 1 + √N k∈INj∈I zj = N ,(rij (z)−rji (z)) =

j∈I

1

�
1 �
3

N 2 j∈Iθ
1−N�

(σi −σj )zi zj +θ �

N 2 k∈I�N m̃ki − zi

�
j∈I

m̃kj

�

zk

�

σk
1+ √
N

Thus,
√

�
�
1 � −1θN hi(r̃ij (z̃) − r̃ji (z̃)) =h1−(σi − σj )z̃i z̃j
N j∈I j
N
j∈I
�
� �
�
�
θ � −1σk−1+ hi 3hk
N m̃ki − hi z̃im̃kj z̃k 1 + √N
N 2 k∈I
j∈I
�
�
�
�
�
1 �
θ
θ
=√
1−
(σi − σj )z̃i z̃j +1−(σi − σj )z̃i z̃j
N
N
N j∈A
j�∈A
�
� �
�
�
1 �
σk
+θ
hi m̃ki − z̃im̃kj z̃k 1 + √
N
N
j∈I
k∈A
�
� �
�
� √
�
1
σk
+θ
N hi m̃ki − √ z̃im̃kj z̃k 1 + √
.
N
N
j∈I
k�∈A
�

5

�

If we now take the large population limit we find,

��




�
�
(σ
−
σ
)z̃
√
j j z̃i + θj�∈A ik�∈A m̃ki z̃klim N hi(r̃ij (z̃) − r̃ji (z̃)) = �
�
N →∞

�

j∈I

(σ
−
σ
)z̃
j j z̃ij�∈A i�

if i ∈ A
otherwise.

Proceeding similarly,

lim

N →∞

√

N

�

�

h2i

i∈I

(r̃ij (z̃) + r̃ji (z̃)) =

j∈I

while
lim

√

N →∞

N

�

hi hj r̃ij (z̃)

i�=j




2z̃

i



0

if i ∈ A
otherwise,

∂ 2f
= 0.
∂ z̃i ∂ z̃jThus, the limiting process has generator(Af )(z̃) =

�
i�∈A

��

�
j�∈A

(σi − σj )z̃j

� �
z̃i

+

�
i∈A

∂f
∂ z̃i
��

�
j�∈A

(σi − σj )z̃j�z̃i + θ�k�∈Am̃ki z̃k

�

∂f
∂ 2f
+ z̃i 2 ,
∂ z̃i
∂ z̃i

i.e. the limiting system may be decomposed into two components since the first sum describes
the dynamics of the non-adaptable genotypes (i ∈/ A) and the second sum describes the
dynamics of the adaptable genotypes (i ∈ A). For readability, in what follows we distinguishbetween the adaptable and non-adaptable types by defining

Z̃i =




X

i



 Yi

if i ∈/A
if i ∈ A

Decomposing the limiting system into equations for the non-adaptable and adaptable types,
6

we have that the non-adaptable types correspond to the terms
�
i�∈A

��

�
j�∈A

(σi − σj )xj

�

xi

�

∂f
∂xi

This describes a system which evolves deterministically (i.e which has no stochastic component), according to the dynamical system (see e.g. Gardiner (2004))
dXi
=
dt

�

�
j�∈A

(σi − σj )Xj

�

Xi ,

which is Eq. 2 of the main text. However, the adaptable types correspond to
�
i∈A

��

�
j�∈A

(σi − σj )xj�yi + θ

�k�∈Am̃ki xk�

∂f∂ 2f+ yi 2 .
∂yi∂yi

Therefore, given X(t), each adaptable type evolves according to an independent diﬀusion
processes, with non-autonomous rates dependent upon the non-adaptable types:

dYi =

��

�
j�∈A

(σi − σj )Xj�

Yi + θ

�k�∈A�m̃ki Xk dt +�
2Yi dBi (t),for independent Brownian motions Bi (t), i ∈ A (Gardiner, 2004), which is Eq. ?? of themain text.

1.3

The Waiting Time Distribution When Adaptable Types are
Rare

In the previous section we showed that the population at each adaptable genotype is described by an independent diﬀusion process. We are interested in the time it takes for the7population to produce an individual with the newly adaptive genotype following an environmental shift. Because adaptable genotypes evolve independently it suﬃces to consider thewaiting time problem for a single diﬀusion,dY = (σY + M ) dt +√2Y dB(t),�where for readability we omit the subscript i, and write σ = j�∈A (σi −σj )Xj for the selection�term and M = θ k�∈A m̃ki Xk for the influx to each adaptable type. The first newly adaptive
genotype occurs when an individual at an adaptable genotype mutates to the newly adaptive

genotype. For a given adaptable genotype the waiting time T for the arrival of the first newly
adaptive mutant to arise is then given by the distribution��
P {T > t} = E e−ζV (t) ,where ζ is the rate at which the newly adaptive type is produced from an adaptable typethrough mutation and V (t) is a stochastically varying rate function which depends on the
population at the adaptable genotype

V (t) =

� t

Y (s) ds

(1)

0

We now make two observations:

1. While the process V (t) is not by itself Markovian, the pair (Y (t), V (t)) does satisfy theMarkov property. Now (1), is equivalent to the SDE dV (t) = Y (t) dt, so (Y (t), V (t))

8satisfies the system of SDEs
dY = (σY + M ) dt +√2Y dB(t),dV = Y dt.
Thus, their joint density satisfies the Forward Kolmogorov Equation (Gardiner, 2004)
∂u
∂
∂
∂2
= − [(σy + M )u] −[yu] + 2 [yu]
∂t
∂y
∂v
∂y

(2)

with initial condition u(y, v, t) = δ(y0 ,0) , the Dirac point mass at the single point (y0 , 0).
2. If we set��
φ(p, q, t) = E e−pY (t)−qV (t) =def

� ∞� ∞
0e−py−qv u(y, v, t) dy dv,0which is the joint Laplace transform of Y (t) and V (t), we have
P {T > t} = φ(0, ζ, t).
Taking the Laplace transform of Eq. 10 yields
�
�
∂φ
∂φ
∂
2
= (−p + σp + q)
− M pφ + lim (σy + M )u −[yu] − yu + lim vu.
y→0
v→0
∂t
∂p
∂y
Now,
lim (σy + M )u −

y→0

9

∂
[yu]
∂y

is the flux, J(v, t) across the y = 0 boundary, which we assume to be 0 on biological grounds(see Draghi et al. (2010) §1.2.2 for a more detailed discussion on boundary conditions), while
lim yu = 0.

y→0

and
lim vu = 0.

v→0We are thus left with solving∂φ∂φ+ (p2 − σp − q)
= −M pφ,∂t∂pwhich may be done via the method of characteristics. To that end, we solvedpp2 − σp − q

= dt,

�
�
�2 �
2
which, upon completing the square, p2 − σp − q = p − σ2 − q + σ4 , has solutionp(t) =

2α(q)
σ
−
α(q)
+
1 − Ce2α(q)t2(3)

where C is a constant of integration that identifies the characteristic curve and we have set
��
�
2�
�
σ
α(q) = ��q + ��.
4

defChoosing the form for p(t) given in Eq. 11, we now set φ̃(t) = φ(p(t), q, t), so thatdφ̃
= −M p(t)φ̃.
dt

10This linear equation has solution�tφ̃(t) = e−M 0 p(τ ) dτ F (C),for an arbitrary function F . Integrating Eq.11 yields���� (1 − C)e2α(q)t � � σ�+p(τ ) dτ = ln ��−α(q)t1 − Ce2α(q)t �
2
0

� t

We can now recover φ(p, q, t) by noting from Eq. 11 thatC=p − α(q) − σ2 −2α(q)tep + α(q) − σ2so that, after some simplification��
� e2α(q)t − p−α(q)− σ2 �−M
�
�
σ
�
p
−
α(q)
−
σ
p+α(q)− σ2 �
−M
−α(q)
t
−2α(q)t2()
�2
φ(p, q, t) = ��
e
F
e
�
2α(q)
p + α(q) − σ2σ��p+α(q)−
2

Lastly, to determine F , we observe thatφ(p, q, 0) ==

� ∞� ∞
�0 ∞ �0 ∞
0

0

e

−py0

e−py−qv u(y, z, 0) dy dv
e−py−qv δ(y0 ,0) (y, v) dy dv = e−py0 ,so that=F�p − α(q) − σ2p + α(q) − σ2

�

i.e.
2α(q)

F (x) = e−( 1−x + 2 −α(q))y011

σThus,




�
�


2α(q)
� e2α(q)t − p−α(q)− σ2 �−M− p−α(q)− σ+ σ2 −α(q)y0
σ ��
2 e−2α(q)tσ1−p+α(q)− 2
σp+α(q)− 2�φ(p, q, t) = ��e−M ( 2 −α(q))t e
�
2α(q)
�
�
p+α(q)− σ2In order to recover the desired waiting time distribution we must set p = 0 and q = ζ, whichyields




��
�


2α(ζ)
− α(ζ)+ σ+ σ2 −α(ζ)y0
� α(ζ) − σ � eα(ζ)t + �α(ζ) + σ � e−α(ζ)t �−M
2 e−2α(ζ)t
σ
1+
�
�
2
2
α(ζ)− σ2
P {T > t} = �
e−M 2 t e
�
�
�
2α(ζ)
(4)

Notice that taking σ = 0 and y0 = 0, we get

P {T > t} = sechM

�� �
ζt ,

which is, up to factors of two owing to the diﬀerence between Wright-Fisher and Moranmodels, the result given in (Draghi et al., 2010).
Eq. 12 gives the waiting time distribution for a given adaptable genotype i to giverise to an individual with the newly adaptive genotype. Since the population at each of
the diﬀerent adaptable genotypes evolves according to an independent diﬀusion process, we
can calculate the waiting time distribution for each adaptable genotype independently, andrecover the expected waiting time for a full network of genotypes to give rise to the newly
adaptive genotype. This is calculated as follows: take P {Ti > t} to be the probability that
we wait for longer than t for a newly adaptive genotype to arise from adaptable genotype i.
Because these distributions are independent, the probability that we wait for longer than t
for a newly adaptive genotype to arise from any adaptable genotype is just the product over�
the distributions, i.e.
i P {Ti > t}. The mean waiting time, E [T ], for a newly adaptive
12

genotype to arise anywhere in the network is then given by
E [T ] =

� ∞�
0

i

P {Ti > t} dt(5)

In the simplified model with genotypes of only high and low robustness, this can be reduced
to the product over two distributions – the waiting time for a newly adaptive type to arise
from a high robustness adaptable genotype, P {TH > t}, and the waiting time for a newly
adaptive type to arise from a low robustness adaptable genotype, P {TL > t}, where wemust take the influx term M in each distribution to be the total influx to all high and low
robustness adaptable types respectively. These are given in Eqs. 6-7 as

MH = qθK(πH XH + (1 − πL )XL )P

ML = qθK((1 − πH )XH + πL XL )P

andThe selection terms, σH and σL used for each distribution are given in Eqs. ??–??. The
expected waiting time, E [T ], for a newly adaptive genotype to arise in the simplified modelis then given by
E [T ] =� ∞0

P {TH > t} P {TL > t} dt

It is this quantity, calculated for the simplified model assuming that adaptable types are sufficiently rare that prior to the environmental shift no adaptable types exist in the population,
which is used in our results (Figures ??–6).

13

1.4

Limit when adaptable individuals are common

In this section we determine the adaptation time when the numer of adaptable individuals
at the time of the environmental shift is O(N ). In this case we say that the adaptable√
individuals are “common”, compared to the case in which they are rare (i.e. O( N )) as
described in the main text. We assume that a subset A� ⊆ A of the adaptable genotypes√
are represned by O(N ) individuals, so that Zi (t) = N Z̃i (t) for i ∈ A� , Zj (t) = N Z̃j (t) forj ∈ A − A� .
To investigate the first arrival of the advantageous type, we will consider the process stopped
when the first advantageous type arrives: we introduce a “graveyard” state ∆, and assumethatri,∆ (z) = µmi,∆ zi

for i ∈ A,r∆,i = 0 for all i ∈ I
i.e. we assume the process jumps to the graveyard state whenever some individual produces
an oﬀspring of the selected type, whilst the process, once in the graveyard state, remains
there. This extended process has generator�i,jrij (z) (f (z + ei − ej ) − f (z)) −�i∈AProceeding as before, we takeZ̃i (t) = hi Zi (t),14ri,∆ (z) (f (∆) − f (z)) .

where now
hi =




 √1

if i ∈ A − A�

N



1

otherwise,N

and time is left unscaled. Taking the limit as N → ∞, we find that the limiting generator is�

i∈A�θmi,∆ z̃i (f (∆) − f (z̃)) .

Thus, the process reduces to a simple Markov jump process, in which after an exponentiallydistributed time with rate�θmi,∆ z̃i ,i∈A�the process jumps to the graveyard state, i.e. an individual of the advantageous type is
produced, while the frequencies of each type remain fixed at their initial values over thistimescale. Further, the rate at which advantageous oﬀspring are produced from an individualof type i ∈ A� isθmi,∆ z̃i ,from which we conclude that the probability that the first advantageous individual had aparent of type i ∈ A is

P {type i} =

while




 � θmi,∆ z̃i0j∈A� θmj,∆ z̃j

P {T ≥ t} = e−(�

if i ∈ A�
otherwise,i∈A� θmi,∆ z̃i)t .

This expression gives the waiting time before arrival of an individual with the newly adaptive
genotype, in the case that (some) adaptable types are common at the time of the environ15

mental shift. In contrast to the case in which all adaptable types are rare, we see that thewaiting time in this case is exponentially distributed, and that the individual who gives riseto the newly adaptive genotype aways comes from an adaptable genotype i ∈ A� , i.e. thecommon subset of adaptable types at the time of the environmental shift. As a result, there
is a diﬀerent relationship between between environmental robustness, rate of environmental
perturbations, etc. compared to the case case in which adaptable types are rare.

2

Relaxation of model assumptions

In this section we relax each of the main model assumptions used in the main text, and werevise our analysis accordingly. We demonstrate that our qualitative results are unchanged
by relaxing these assumptions.

2.1

Mutation rates:

In our simulations we have chose a per-genome mutation rate of µ = 10−3 . The reason forthis is as follows: for a non-synonymous substitution rate of 10−8 , µ = 10−3 is equivalent toa genome of 105 nucleotides. The rate of mutations that result in a newly adaptive genotype
is not µ = 10−3 but is µ(1 − q)/K (see Fig 2 of the main text). In our simulations we havetypically picked K = 5 and q = 0.5 so that µ(1 − q)/K = 10−4 . It corresponds to around
104 possible adaptive mutations (given our assumed per-nucleotide substitution rate). It
therefore corresponds to the case in which many genes are initially involved in adaptation.
Mutation rates between other genotype classes can be calculated in a similar way. Smaller
mutation rates between genotype classes can be achieved by choosing diﬀerent values forother parameters,. For example, taking K = 100 results in a rate of adaptive mutations16

of ∼ 10−6 . However doing this adds no new insight since we have catalogued the way theadaptation time varies with all of our model parameters.

2.2

Correlated phenotypic neighbourhoods:

In the main text we assume that, following a mutation, the phenotypic neighbourhood of
the resulting individual is completely redrawn – that, the set of its K phenotypic neighbours
is drawn uniformally from the set of P alternatives. This results uncorrelated phenotypicneighbourhoods amoung diﬀerent genotypes. We relax this assumption by introducing a parameter f , which is the fraction of the K phenotypes in a individual’s phenotypic neighbourhood that are redrawn following mutation. When f = 1 we recover the case of completely
uncorrelated phenotypic neighbourhoods. When f = 0, all genotypes have the same phenotypic neighbourhood. Introducing f results in a modified mutation rate between adaptable
and non-adaptable types. Specifically, the rate of mutation from non-adaptable to adapt��able types is µqf K
and from adaptable to non-adaptable types is µq 1 − f K. All otherPP

mutation rates and selective coeﬃcients remain unchanged from those presented in the maintext. The eﬀect of decreasing f is therefore identical to increasing P . Decreasing the rateof influx to adaptable types tends to slow overall adaptation time, but it does not alter thequalitative behaviour of the model described in the main text, which depends on the selective
coeﬃcients. This is illustrated in Fig. S1 for the variation in adaptation time with rate ofenvironmental perturbations, �.

[Figure 1 about here.]

17

2.3

Variation in the size of phenotypic neighbourhoods with genotype:

Next we relax the assumption that all genotypes have exactly K phenotypic neighbours.
To do this we instead assign each of the P alternative phenotypes to the phenotypic neighbourhood of a genotype with probability k, following mutation. The size of the phenotypicneighbourhood, K, for a given genotype is therefore drawn from a binomial distribution
with mean K̄ = kP . Because phenotypic neighbourhood size varies with genotype, diﬀerent adaptable types have diﬀerent selection coeﬃcients. Adaptable types can therefore nolonger be simply separated into categories of high and low robustness. Instead they must be
separated into categories according to their environmental robustness φ and the size of theirphenotypic neighbourhood K. The selective coeﬃcient for a high robustness genotype with
K phenotypic neighbours is given byσH (K)1+s√=�
(1 − φH ) + � (φH − φL ) XL
K
N

(6)

and for low robustness genotypes with K phenotypic neighbours by

σL (K)1+s√=�
(1 − φL ) − � (φH − φL ) XH .
K
N

Let YHK is the number of individuals, rescaled by√(7)

N , at high robustness genotypes with K

phenotypic neighbours. When adaptable types are rare, the population at such genotypes
evolves according to18

dYHK =�K
σH (K)YHK + G(K)qθ ((1 − πH )XH + πL XL )P

Similarly, let YLK is the number of individuals, rescaled by

√

�

dt +

�
2YHK dBH (t)(8)

N , at low robustness genotypes

with K phenotypic neighbours. The population at these genotypes evolves according todYLK =�K
σL (K)YLK + G(K)qθ (πH XH + (1 − πL )XL )
P

�

dt +

�
2YLK dBL (t)(9)

where G(K) is the probability that the genotype has K neighbours,�P � Kk (1 − k)P −K
G(K) = �P K � P �K
P −K
K=1 K k (1 − k)Where we assume that all genotypes must have at least one phenotypic neighbour (i.e K ≥ 1).
Finally, the rate at which individuals with the newly adaptive genotype are produced thorough mutation at adaptable genotypes with K phenotypic neighbours is µ 1−q.
K

The full distribution of adaptation times from each adaptable type can be found (see Sup��
porting Information). Let P THK > t be the waiting time distribution from high robustness��
genotypes with K phenotypic neighbours, and P TLK > t be the waiting time distribution
from low robustness genotypes with K phenotypic neighbours. The expected waiting time
for an individual with the newly adaptive phenotype to arise is then19E [T ] =

� ∞�
P
0

K=1

�
� �
�
P THK > t P TLK > t dt.

Allowing K to vary does not change the set of qualitative behaviours described in the maintext. However it does change the parameter values for which those behaviours occur. Thisis because the observed behaviours depend on whether low robustness adaptable types are
selectively favoured compared to the population prior to the environmental shift. This inturn depends on the ratio Ks (see Table 2, main text). Because K is distributed acrossa range of values, this tends to result in some low robustness genotypes being selectively
favoured and some disfavoured. Despite this, all the qualitative behaviours described in themain text are still observed, as shown in Figure S2.

[Figure 2 about here.]

2.4

Inclusion of non-lethal deleterious mutations:

We now relax the assumption that all deleterious mutations are lethal, i.e that alternative
phenotypes have fitness 0. To do this we introduce four new classes of genotypes. Previously
we had classes• A – High robustness and adaptable• B – Low robustness and adaptable• C – High robustness and non-adaptable• D – Low robustness and non-adaptable
20

Each of these four classes are now split into two – fit and deleterious. The phenotype coded
for by the “fit” genotypes has fitness 1 and the phenotype coded for by the “deleterious” geno� �
types has fitness 1 − d, where d is O √1N or larger. As before, environmental perturbations
may result in phenocopy, such that individuals develop the phenotype of their mutational

neighbours. As a result genotypes Cf it and Df it have fitness 1−d(1−φH )� and 1−d(1−φL )�
respectively. Genotypes Cdel and Ddel are assumed to both have fitness 1 − d – i.e we neglect
the influence of phenocopy on the fitness of these genotypes. Genotypes Af it and Bf it havefitness 1 − d(1 − φH )� + s+d(1 − φH )� and 1 − d(1 − φL )� + s+d(1 − φL )� respectively. Finally,K
K
genotypes Adel and Bdel have fitness 1 − d + s+d(1 − φH )� and 1 − d + s+d(1 − φL )� respecKK
tively. Mutations between the diﬀerent classes occur as described for the previous model,
with the modification that mutations at a “fit” genotype which are non-neutral result in a“deleterious” genotype. Mutations at a “deleterious” genotype which are non-neutral resultin a “fit” genotype with probability K1 . The resulting model has four classes of adaptablegenotype, with selective coeﬃcientsσ f itd+s√H = �
(1 − φH ) + �d (φH − φL ) XL
K
N

(10)

σ del
d+s
√H = �
(1 − φH ) + �d (φH − φL ) XL − d (1 − �(1 − φH ))
K
N

(11)

for genotypes Af it ,for genotypes Adel ,21

σ f itd+s√L = �
(1 − φL ) − �d (φH − φL ) XH
K
N

(12)

σ del
d+s
√L = �
(1 − φL ) − �d (φH − φL ) XH − d (1 − �(1 − φL ))
K
N

(13)

for genotypes Bf it , andfor genotypes Bdel . We have assumed that, prior to the environmental shift, the whole population can be treated as being only at “fit” genotypes.

The full distribution of adaptation times from each adaptable type can be found (see Sup�
�
�
�
f it
porting Information) for this model. Let P TH > t and P THdel > t be the waiting time
�
�
�
�
f it
distribution for high robustness genotypes, and P TL > t and P TLdel > t be the waitingtime distribution for low robustness genotypes. The expected waiting time for an individualwith the newly adaptive phenotype to arise is then

� ∞ �
� �
� �
� �
�
E [T ] =P THf it > t P TLf it > t P THdel > t P TLdel > t dt.
0

Once again, this model produces the same set of qualitative behaviours as described in themain text. However it changes the parameter values for which those behaviours occur. In
particular, the diﬀerent selective regimes corresponding to diﬀerent behaviours, described in
H −φL
Table 2 of the main text are 1+s< 1 and 1+s< φ1−φ. In this model the corresponding
K
K
L
H −φL
selective regimes are d+s< d and d+s< d φ1−φ. This has the eﬀect that, if d < 1, the
K
K
L

22

selective regime described as “unrealistic” in the main text, requiring large values of s, cannow occur for smaller values of s, making the full range of behaviours described in the main
text attainable under more realistic sets of parameters. The full range of behaviours of thismodel are shown in Figure S3.

[Figure 3 about here.]

2.5

Adaptability when low robustness adaptive genotypes are selected against:

We may also consider the case when newly adaptive genotypes with low robustness have lower
fitness than the mean population fitness. This occurs in the following scenario: Assume that
newly adaptive genotypes inherit the environmental robustness of the individual that gave
rise to them. Therefore high robustness adaptive genotypes have fitness (1 + s)(1 − (1 −φH )�) + K1 (1 − φH )� and low robustness adaptive genotypes have fitness (1 + s)(1 − (1 −φL )�) + K1 (1 − φL )�. The mean fitness of the population at the environmental shift is simply
1 − (1 − φH )�, therefore low robustness adaptive genotypes are only selectively favoured if(1 + s)(1 − (1 − φL )�) +1(1 − φL )� > 1 − (1 − φH )�
K

i.e ifs>� ((φH − φL )K − (1 − φL ))K(1 − (1 − φL )�)otherwise they are selected against. When this is the case we consider that only high23

robustness adaptive genotypes are truly adaptive, and calculate the waiting time until an
individual with this genotype arises. This is easy to do since low robustness adaptable genoL)
types produce high robustness adaptive genotypes at rate (1−q)(1−πand high robustnessK
H
adaptable genotypes produce high robustness adaptive genotypes at rate (1−q)π. We canKneglect the eﬀect of low robustness adaptive genotypes mutating to produce high robustness
adaptive genotypes, since we assume that the population is initially at the non-adaptablegenotypes which cannot mutate directly to the newly adaptive genotypes. Using these modified mutation rates we determine the waiting time until a high robustness adaptive genotype
arises (assuming that s is suﬃciently small that low robustness adaptive genotypes are selected against). These waiting times are shown varying with the model parameters in Figure
S4 below. As can be seen the same qualitative behaviours as described for our model in theregime 1+s< 1 are recovered in this case, with the addition that increasing low robustness
Kgenotype clustering increases adaptation time (since it slows the rate of mutations from low
robustness adaptable genotypes to high robustness adaptive genotypes).

[Figure 4 about here.]

2.6

Other influences of clustering of high-robustness genotypes on
adaptation time:

We also examined how the clustering of high-robust genotypes influences the relationship
H −φL
between environmental perturbations and mean adaptation time, in the regime 1+s< φ1−φ.
KL

As described in the main text, in this regime increased clustering of high-robust genotypes
facilitates adaptation, regardless of the amount of environmental noise. In addition, Figure
5 shows that the position of the maximum adaptation time occurs at a lower value of �,
when genotypes are more clustered. In other words, clustering of robust genotypes broadens24the range over which increasing the frequency of environmental perturbation will facilitate
adaptation.

[Figure 5 about here.]

Dependence of adaptation time on mutational robustness:
We have also studied how environmental noise can influence the relationship between mutational robustness and adaptability. Previous work has shown that, in the absence ofenvironmental noise, intermediate levels of mutational robustness produce the fastest rates
of adaptation (Draghi et al., 2010; Wagner, 2008). But how does adaptation time dependon the amount of mutational robustness, q, when environmental noise is present?
We find that introducing environmental noise does not qualitatively change the way inwhich adaptation time varies with mutational robustness: intermediate levels of mutational
robustness still produce the faster rates of adaptation. Nonetheless, environmental noise
does alter the quantitative aspects of this relationship (Figure 6), in two ways.
Firstly, noise can either increase or decrease overall adaptation time compared to thecase when environmental noise is absent. The direction of this eﬀect depends on whether
adaptable individuals with low environmental robustness are selectively favored or disfavored,
H −φL
H −φL
that is whether 1+s> φ1−φor 1+s< φ1−φ. As above, there is a simple intuition forK
K
LL

this behavior. In the first of these regimes, both types of adaptable individuals gain adirect selective benefit in the presence of noise, speeding adaptation. In the second regime,
however, the population concentrates at genotypes with high environmental robustness, andadaptation is generally slower (especially when noise is frequent).

25

Aside from changing the overall rate of adaptation, environmental noise also changes thevalue of q that maximizes the rate of adaptation (Figure 6). That is, environmental noise
alters the level of genetic robustness that is optimal for rapid adaption.

[Figure 6 about here.]Time series:
We plot time series data for the 4 diﬀerent classes of genotype. Prior to the environmental
shift the population is in equilibrium (Fig. 7). After the environmental shift individuals
start to accumulate at adaptable genotypes due to phenocopy. Fig. 9 shows time series forparameters identical to those shown in Fig.8, with the exception that πL has been varied(πL = 0.5 in Fig. 7 and πL = 0.1 in Fig. 9). The fact that these two plots are the same
illustrates that piL has virtually no eﬀect either on the equilibrium prior to the environmental
shift or on the adaptive dynamics following the environmental shift. Fig. 8 shows time seriesplots for the same parameters as Fig. 7, with the exception that s has been increased (s = 1
in Fig. 8 and s = 6 in Fig. 8), such that the two cases lie in diﬀerent selective regimes. We
see that in Fig. 8, in contrast to Fig. 7, individuals accumulate rapidly at low robustness
adaptable genotypes.

[Figure 7 about here.][Figure 8 about here.][Figure 9 about here.]26

References
Draghi, J., T. Parsons, G. Wagner, and J. Plotkin, 2010. Mutational robustness can facilitate
adaptation. Nature 436:353–355.
Etheridge, A. M., 2011. Some Mathematical Models from Population Genetics, Lecture Notes
in Math, vol. 2012. Springer, Berlin Heidelberg.
Ethier, S. and T. Kurtz, 2005. Markov processes : characterization and convergence. Wiley
Interscience, Hoboken.
Gardiner,C. W.,

2004.

Handbook of stochastic methods for physics,

istry, and the natural sciences.

3rd ed ed. Springer-Verlag, Berlin.

chemURL

http://www.loc.gov/catdir/enhancements/fy0818/2004043676-d.html.
Wagner, A., 2008. Robustness and evolvability: a paradox resolved. Proceedings of theRoyal Society B: Biological Sciences 275:91–100.

27

List of Figures
1

Correlated phenotypic neighbourhood . . . . . . . . . . . . . . . . . . . . . .

29

2

Variation in size of phenotypic neighbourhood . . . . . . . . . . . . . . . . .

30

3Non-lethal deleterious mutations . . . . . . . . . . . . . . . . . . . . . . . . .

31

4Deleterious low robustness adaptive genotype . . . . . . . . . . . . . . . . .

32

5Adaptation, environmental noise, and clustering . . . . . . . . . . . . . . . .

33

6Adaptation and mutational robustness . . . . . . . . . . . . . . . . . . . . .

34

7Time series 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

8Time series 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

9Time series 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

28

Adaptation time �generations�

5000
4000
3000
2000
1000
0
0.0

0.2

0.4

0.6

0.8

1.0

Frequency of environmental perturbations, ΕFigure 1: Correlated phenotypic neighbourhood – Mean adaptation time as a function of
the probability of environmental perturbation, �. Plots show diﬀerent levels of correlation
between phenotypic neighbourhoods, with f = 1 (red, uncorrelated), f = 0.5 (blue) and f =
0.3 (black). As phenotypic neighbourhoods become more correlated, mean adaptation timeincreases. Plots show means of 10,000 replicate Monte-Carlo simulations with populationsof N = 10, 000 individuals, with µ = 0.001, P = 100, K = 5, q = 0.5, φH = 0.8, φL = 0.2
and πH = 0.9 with s = 2.

29

Adaptation time �generations�

Adaptation time �generations�

1400
1200
1000
800
600
400
200
0
0.0

0.2

0.4

0.6

0.8

1000
800
600
400
200
0
0.0

1.0

High robustness genotype clustering, ΠHL

0.2

0.4

0.6

0.8

1.0

Strength of low robustness genotypes, ΦL

1000

Adaptation time �generations�

Adaptation time �generations�(a) Mean adaptation time as a function of high ro- (b) Mean adaptation time as a function of robustbustness genotype clustering, πHness of low robust genotypes, φL

800
600
400
200
0
0.0

0.2

0.4

0.6

0.8

1.0

700
600
500
400
300
200
100
0
0.0

Strength of low robustness genotypes, ΦH

0.2

0.4

0.6

0.8

1.0

Frequency of environmental perturbations, ΕAdaptation time �generations�(c) Mean adaptation time as a function of robust- (d) Mean adaptation time as a function of the probness of high robust genotypes, φHability of environmental perturbation, �
2000
1500
1000
500
0
0.0

0.2

0.4

0.6

0.8

1.0

Mutational Robustness, q(e) Mean adaptation time as a function of mutational robustness, q

Figure 2: Variation in size of phenotypic neighbourhood – Plots show populations of N =10, 000 individuals, with µ = 0.001, P = 100, K̄ = 5, q = 0.5, φH = 0.9, φL = 0.1 andπH = 0.5 with s = 2 (red line) and s = 8 (blue line). Lines indicate the analytic solution to
our model, whereas dots indicate the means of 10,000 replicate Monte-Carlo simulations.
30

Adaptation time �generations�

Adaptation time �generations�

1500

1000

500

0
0.0

0.2

0.4

0.6

0.8

1.0

1000
800
600
400
200
0
0.0

High robustness genotype clustering, ΠHL

0.2

0.4

0.6

0.8

1.0

Strength of low robustness genotypes, ΦL

1000

Adaptation time �generations�

Adaptation time �generations�(a) Mean adaptation time as a function of high ro- (b) Mean adaptation time as a function of robustbustness genotype clustering, πHness of low robust genotypes, φL

800
600
400
200
0
0.0

0.2

0.4

0.6

0.8

1.0

1000
800
600
400
200
0
0.0

Strength of low robustness genotypes, ΦH

0.2

0.4

0.6

0.8

1.0

Frequency of environmental perturbations, ΕAdaptation time �generations�(c) Mean adaptation time as a function of robust- (d) Mean adaptation time as a function of the probness of high robust genotypes, φHability of environmental perturbation, �
1400
1200
1000
800
600
400
200
0
0.0

0.2

0.4

0.6

0.8

1.0

Mutational Robustness, q(e) Mean adaptation time as a function of mutational robustness, q

Figure 3: Non-lethal deleterious mutations – Plots show populations of N = 10, 000 individuals, with µ = 0.001, P = 100, K = 10, q = 0.2, φH = 0.9, φL = 0.1 and πH = 0.5, d = 0.2
with s = 0.5 (red line) and s = 3 (blue line). Lines indicate the analytic solution to ourmodel, whereas dots indicate the means of 10,000 replicate Monte-Carlo simulations.

31Mean adaptation time �generations�Mean adaptation time �generations�

6000
5000
4000
3000
2000
1000
0
0.0

0.2

0.4

0.6

0.8

1.0

1500
1400
1300
1200
1100
1000
0.0

High robustness genotype clustering, ΠH

0.2

0.4

0.6

0.8

1.0

Low robustness genotype clustering, ΠL

Adaptation time �generations�

Adaptation time �generations�(a) Mean adaptation time as a function of high ro- (b) Mean adaptation time as a function of robustbustness genotype clustering, πHness of low robust genotypes, πL
1400
1200
1000
800
600
400
200
0
0.0

0.2

0.4

0.6

0.8

1.0

1400
1200
1000
800
600
400
200
0
0.0

Strength of low robustness genotypes, ΦH

0.2

0.4

0.6

0.8

1.0

Strength of low robustness genotypes, ΦL

Mean adaptation time �generations�(c) Mean adaptation time as a function of robust- (d) Mean adaptation time as a function of robustness of high robust genotypes, φH
ness of high robust genotypes, φL
12 00011 500
11 00010 500
10 000
9500
9000
0.0

0.2

0.4

0.6

0.8

1.0

Frequency of environmental perturbations, Ε(e) Mean adaptation time as a function of the probability of environmental perturbation, �Figure 4: Deleterious low robustness adaptive genotype – Plots show populations of N =10, 000 individuals, with µ = 0.001, P = 100, K = 5, q = 0.5 and s = 0.05. In (a), (b),(d) and (e) φH = 0.9, (a), (b), (c) and (e) in φL = 0.1, in (c) and (d) πH = 0.5 and in (e)πH = 0.9 , in (c) and (d) πL = 0.5 and in (e) πL = 0.1, in (a), (b), (c) and (d) � = 0.1. Lines
indicate the analytic solution to our model.
32

Mean adaptation time �generations�

2000
1500
1000
500
0
0.0

0.2

0.4

0.6

0.8

1.0

Frequency of environmental perturbations, ΕFigure 5: Mean adaptation time as a function of the probability of environmental perturH −φLbation, �, in the regime 1+s< φ1−φ. In general, adaptation is faster when high-robustKLgenotypes are more clustered (blue line) as opposed to less clustered (red line). Moreover,the clustering of robust genotypes broadens the range over the rate of adaptation is positively correlated with the rate of environmental perturbation. Plots show populations of
N = 10, 000 individuals, with µ = 0.001, P = 100, K = 5, q = 0.5, φH = 0.9, φL = 0.1 ands = 2 with πH = 0.1 (red line) and πH = 0.3 (blue line). Lines indicate the analytic solution
to our model, whereas dots indicate the means of 10,000 replicate Monte-Carlo simulations.

33

Mean adaptation time �generations�

2500
2000
1500
1000
500
0
0.0

0.2

0.4

0.6

0.8

1.0

Mutational Robustness, q
Figure 6: The mean adaptation time as a function of mutational robustness, q. There is
non-monotonic relationship between mutational robustness and mean adaptation time, andthe shape of this trend is preserved even in the presence of environmental noise. However,H −φL
the addition of environmental noise tends to facilitate adaptation in the regime 1+s> φ1−φ,K
Land it retard adaptation otherwise. Moreover, the presence of environmental noise changes
the value of q that maximizes the rate of adaptation. Plots show populations of N = 10, 000
individuals, with µ = 0.001, � = 0.1 P = 100, K = 5, φH = 0.9, φL = 0.1 and πH = 0.5 withs = 2 (red line) and s = 5 (blue line). The case with no environmental noise (� = 0, blackline) is also shown. Lines indicate the analytic solution to our model, whereas dots indicate
the means of 10,000 replicate Monte-Carlo simulations.

34

0.075

0.95

0.07

0.945

0.065
Frequency

Frequency

0.955

0.94

0.06

0.935

0.055

0.93

0.05

0.925
0

200

400
600
Time (generations)

800

1000

0.045
0

(a) Non-adaptable high robustness
"'*(

+,&!

200

400
600
Time (generations)

800

1000

(b) Adaptable high robustness

ï-

)'(

*+&!

ï#

"'*
)
"'%(
"'(
93/:;/2<=;51<=14>?

"'%
"')(
"')

"

&'(

"'$(
"'$

&
"'((
"'(
!

"!!

#!!
$!!
./01,23141567/849:

%!!

&!!!

(c) Non-adaptable low robustness

!'(
!

"!!

#!!
$!!
,-./+01/2/345-6278

%!!

&!!!

(d) Adaptable low robustness

Figure 7: Most realistic selective regime with s/K = 0.2, πH = 0.5 and πL = 0.5 a) Nonadaptable high robustness genotypes. This shows that the prior to the environmental shift
most of the population is at this genotype class and that following the environmental shift
there is little change in the population at this genotype class before adaptation occurs (there
is a drop of ∼ 2%). b) Adaptable high robustness genotypes. This shows that prior to
the environmental shift the population at this genotypes class is small (due to their rarity)and that following the environmental shift and the population increases by ∼ 0.02. c) Nonadaptable low robustness genotypes. This shows that the prior to the environmental shift
the population at this genotype class is negligible and that following the environmental shift
there is no obvious change in the population at this genotype class before adaptation occurs.
d) Adaptable high robustness genotypes. This shows that prior to the environmental shift
the population at this genotypes class is small (due to their rarity) and that following theenvironmental shift and the population increases by ∼ 1.2 × 10−4 , much less than the highrobustness adaptable genotypes (b). Plots are for N = 10, 000 individuals, with µ = 0.001,
P = 100, K = 5, q = 0.5 and s = 1, φH = 0.9, φL = 0.1, πH = 0.5, πL = 0.1 and � = 0.1.
Plots show averages over 1000 sample paths.
35

0.96

0.13
0.12

0.94
0.11
0.1
Frequency

Frequency

0.92

0.9

0.88

0.09
0.08
0.07
0.06

0.86
0.05
0.84
0

200

400
600
Time (generations)

800

1000

0.04
0

(a) Non-adaptable high robustness
"'*(

+,&!

200

400
600
Time (generations)

800

1000

(b) Adaptable high robustness

ï-

0.03

"'*

0.025

"'%(
0.02
Frequency;51<=14>?

"'%
"')(
"')

0.015

0.01

"'$(
0.005

"'$
"'((
!

"!!

#!!
$!!
./01,23141567/849:

%!!

&!!!

(c) Non-adaptable low robustness

0
0

200

400
600
Time (generations)

800

1000

(d) Adaptable low robustness

Figure 8: Less realistic selective regime with s/K = 1.2, πH = 0.5 and πL = 0.5 a) Nonadaptable high robustness genotypes. This shows that the prior to the environmental shift
most of the population is at this genotype class and that following the environmental shift
there is a larger but still small change in the population at this genotype class before adaptation occurs (there is a drop of ∼ 10%). b) Adaptable high robustness genotypes. Thisshows that prior to the environmental shift the population at this genotypes class is small(due to their rarity) and that following the environmental shift and the population increasesby ∼ 0.07. c) Non-adaptable low robustness genotypes. This shows that the prior to theenvironmental shift the population at this genotype class is negligible and that following theenvironmental shift there is no obvious change in the population at this genotype class before
adaptation occurs. d) Adaptable high robustness genotypes. This shows that prior to theenvironmental shift the population at this genotypes class is small (due to their rarity) andthat following the environmental shift and the population increases by ∼ 0.03, this is comparable to the increase at the high robustness adaptable genotypes (b). The non-negligible
population at low-robustness adaptable genotypes is the eﬀect discussed in the main tex.
Plots are for N = 10, 000 individuals, with µ = 0.001, P = 100, K = 5, q = 0.5 and s = 6,
φH = 0.9, φL = 0.1, πH = 0.5, πL = 0.5 and � = 0.1. Plots show averages over 1000 sample
36
paths.

0.95

0.075

0.07

0.945

0.065
Frequency

Frequency

0.94

0.935

0.06

0.055
0.93

0.05

0.925
0

200

400
600
Time (generations)

800

0.045
0

1000

(a) Non-adaptable high robustness
+,&!

800

1000

ï#

"'%

"'*

"'$

"'%(

"'#

"'%

"'"

"')(
"')

&'%
&'$

"'$

&'#

"'((

&'"
"!!

#!!
$!!
./01,23141567/849:

%!!

()&!

"

"'$(

"'(
!

400
600
Time (generations)(b) Adaptable high robustnessï-

71-89-0:;;51<=14>?

"'*(

200

&
!

&!!!

(c) Non-adaptable low robustness

"!!

#!!
$!!
*+,-)./-0-123+4056

%!!

&!!!

(d) Adaptable low robustness

Figure 9: Most realistic selective regime with s/K = 0.2, πH = 0.5 and πL = 0.1 The
information here is the same as in the previous graph, but it illustrates that changing low
robustness genotype clustering has no impact on the dynamics of adaptation. Plots are for
N = 10, 000 individuals, with µ = 0.001, P = 100, K = 5, q = 0.5 and s = 1, φH = 0.9,
φL = 0.1, πH = 0.5, πL = 0.1 and � = 0.1. Plots show averages over 1000 sample paths.

37

List of Tables
1

Model Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

39

Parameter

Interpretation

A⊂I

The set of adaptable types

ζi

The rate mutation of the ith adaptable type to the newly adaptive typeµ = Nθ

The rate of mutationm̃ij

The probability that the mutant oﬀspring of an individual of type i has type j

Si = √σNi

The selective advantage of the ith typeZi (t)The number of individuals of type i at time t
Table 1: Model Parameters39