Evolution favors protein mutational robustness in sufficiently
large populations
Jesse D Bloom∗1 , Zhongyi Lu1 , David Chen1 , Alpan Raval2 , Ophelia S Venturelli1 and
Frances H Arnold∗1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA
Keck Graduate Institute of Applied Life Sciences and School of Mathematical Sciences, Claremont Graduate University,
Claremont, CA 91711, USA
Email: Jesse D Bloom∗- jesse.bloom@gmail.com; Zhongyi Lu - lu07@caltech.edu; David Chen - davidc@caltech.edu;
Alpan Raval - alpan raval@kgi.edu; Ophelia S Venturelli - opheliav@stanford.edu; Frances H Arnold∗-
frances@cheme.caltech.edu;
Corresponding author
Abstract
Background: An important question is whether evolution favors properties such as mutational robustness
or evolvability that do not directly benefit any individual, but can influence the course of future evolution.
Functionally similar proteins can differ substantially in their robustness to mutations and capacity to
evolve new functions, but it has remained unclear whether any of these differences might be due to
evolutionary selection for these properties.
Results: Here we use laboratory experiments to demonstrate that evolution favors protein mutational
robustness if the evolving population is sufficiently large. We neutrally evolve cytochrome P450 proteins
under identical selection pressures and mutation rates in populations of different sizes, and show that
proteins from the larger and thus more polymorphic population tend towards higher mutational robust-
ness. Proteins from the larger population also evolve greater stability, a biophysical property that is
known to enhance both mutational robustness and evolvability. The excess mutational robustness and
stability is well described by existing mathematical theories, and can be quantitatively related to the
way that the proteins occupy their neutral network.
Conclusions: Our work is the first experimental demonstration of the general tendency of evolution to
favor mutational robustness and protein stability in highly polymorphic populations. We suggest that
this phenomenon may contribute to the mutational robustness and evolvability of viruses and bacteria
that exist in large populations.
http://arxiv.org/abs/0704.1885v1
Background
Proteins are quite tolerant of mutations, al-
lowing evolution to produce highly diverged
sequences that fold to similar structures and
perform conserved biochemical functions [1, 2].
However, proteins with nearly identical struc-
tures and functions may differ in their robust-
ness to mutation [3–5], as well as in their ca-
pacity to acquire new functions [5]. The fact
that mutational robustness and evolvability can
vary among the functionally equivalent proteins
produced by natural sequence divergence makes
these properties important hidden dimensions
in evolution — direct selection for protein func-
tion is blind to them, yet they can play a cru-
cial role in enabling future evolution. Whether
the evolutionary process somehow promotes the
acquisition of mutational robustness and evolv-
ability therefore remains a major question [6–8].
Previous experiments have identified several
specific evolutionary conditions that can affect
mutational robustness. For example, genetic
complementation decreases the mutational ro-
bustness of viruses [9], while high mutation rates
favor mutational robustness in simulated digi-
tal organisms [10]. However, theory [11] makes
the much broader — and heretofore experimen-
tally untested — prediction that extra muta-
tional robustness will arise quite generally in
sufficiently large populations. This prediction
cannot be understood in the standard frame-
work of Kimura’s neutral theory [12], since one
of the usual assumptions of the neutral theory
is that mutational robustness is constant. (Al-
though Takahata [13] treated the consequences
of stochastically fluctuating neutrality on the
molecular clock, he did not describe how muta-
tional robustness might change systematically
during evolution.) However, changes in mu-
tational robustness can be described by envi-
sioning evolution as occurring on neutral net-
works, or sets of functionally equivalent pro-
teins that are connected by single mutational
steps [14–17]. In a seminal theoretical analysis
of evolution on neutral networks, van Nimwe-
gen and coworkers [11] predicted that the extent
of mutational robustness should depend on the
degree of population polymorphism. Here we
briefly summarize their reasoning, since it moti-
vates our experimental work. We also refer the
reader to chapter 16 of [8], which contains an
excellent explanation of the densely mathemat-
ical work of van Nimwegen and coworkers [11].
If an evolving population is mostly
monomorphic, then each mutation is either
lost or goes to fixation before another muta-
tion occurs. The population is therefore usually
clustered at a single genotype and rarely expe-
riences mutations, meaning that selection does
not distinguish between genotypes of different
mutational robustness. All nodes of the neutral
network are thus equivalent and will be occupied
by the population with equal probability [11].
On the other hand, a highly polymorphic pop-
ulation is always spread across many nodes of
the neutral network. When mutations occur,
the members of the population at highly con-
nected nodes have a better chance of surviving,
causing them to be favored by evolution and
increasing the average mutational robustness
[11, 17–20]. Specifically, a highly polymorphic
population occupies each node with a proba-
bility proportional to its eigenvector central-
ity [11, 17], a measure of how connected it is
to other connected nodes (a variant of eigen-
vector centrality is used by Google’s PageRank
algorithm to rank a webpage’s importance in
the network of internet links [21]). Figure 1A
illustrates how mostly monomorphic and highly
polymorphic populations are predicted to oc-
cupy a neutral network. For proteins, changes
in neutral network occupancy should be man-
ifested by changes in thermodynamic stabil-
ity [22], with proteins from highly polymorphic
populations predicted to be more stable than
their counterparts from mostly monomorphic
populations (Figure 1B). Note that the extent
of polymorphism depends on the product of
the mutation rate and population size, meaning
that protein populations of different sizes are
predicted to evolve to different levels of mu-
tational robustness and stability even if they
experience the same mutation rate.
Results and Discussion
Design of neutral evolution experiment
To test whether high population polymorphism
drives an increase in mutational robustness and
protein stability, we performed laboratory evo-
lution experiments on cytochrome P450 pro-
teins. The basic idea was to neutrally evolve
P450s under a constant selection pressure in
populations that were either monomorphic or
highly polymorphic, and observe whether the
proteins evolved to different levels of mutational
robustness and stability. The evolution experi-
ments started with a P450 BM3 heme domain
that had been engineered to hydroxylate 12-p-
nitrophenoxydodecanoic acid (12-pNCA) [23].
We imposed the selection criterion that Es-
cherichia coli cells expressing the P450 had to
yield lysate with enough active enzyme to hy-
droxylate a specified amount of 12-pNCA in 40
minutes. This criterion roughly corresponds to
the case in which an enzyme must catalyze a
biochemically relevant reaction at some mini-
mal level in order for its host to survive. Note
that other properties such as stability and ex-
pression level can vary freely, provided that the
criterion for total activity is met.
The properties of a neutrally evolving pro-
tein eventually “equilibrate,” much as the prop-
erties of an isolated physical system under some
macroscopic constraint tend towards the values
that maximize the system’s internal entropy.
For proteins, this usually means that stability,
expression, and activity drift towards their low-
est tolerable values, since the vast majority of
random sequences do not encode stable, well-
expressed enzymes (that is, natural selection
must work against sequence entropy to maintain
a functional protein) [22, 24]. The initial P450
had been engineered for maximal activity [23],
meaning that it was not equilibrated to the more
mild selection criterion of the experiments. We
therefore neutrally evolved this initial P450 for
16 generations, introducing random mutations
with error-prone PCR and retaining all mu-
tants that met the selection criterion for total
activity on 12-pNCA. The procedure used for
this equilibration evolution was similar to that
for the polymorphic neutral evolution described
below. As expected, expression, stability, and
activity all dropped during the equilibration
evolution. At the end of the equilibration evo-
lution, we chose a single sequence as the parent
for the neutral evolution experiments. The gene
encoding this parent sequence contained 29 nu-
cleotide mutations and 13 amino acid mutations
relative to the initial P450 (Additional File 1).
We used this parent gene to begin three
parallel sets of neutral evolution experiments,
which we named “monomorphic,” “polymor-
phic,” and “unselected” (Figure 2). The
monomorphic experiments capture the case
where the population moves as a single entity,
the polymorphic experiment captures the case
where the population spreads across many se-
quences, and the unselected experiments show
how the gene evolves in the absence of selection
for protein function. In all experiments, at each
generation we used error-prone PCR to intro-
duce an average of 1.4 nucleotide mutations per
P450 gene (Table 1). The mutant genes were
ligated into a plasmid and transformed into E.
coli [25], and transformants were selected using
the plasmid’s antibiotic resistance marker. For
the unselected case, we randomly picked one of
the mutants, recovered the mutant gene with
a plasmid mini-prep, and used this mutant as
the template for the next generation of error-
prone PCR. We performed four independent
replicates of unselected evolution, evolving each
for 12 generations.
For the monomorphic and polymorphic pop-
ulations, we imposed the selection criterion that
the P450s hydroxylate 12-pNCA with at least
75% of the total activity of the original par-
ent gene. We expressed the P450s in E. coli,
and then assayed the cell lysates for activity in
a high-throughput 96-well plate format. The
total amount of product produced by 80 µl of
clarified lysate in 40 minutes was compared to
the median of four control wells containing the
original parent P450 to determine if the mutant
met the selection criterion. The only differ-
ence between the monomorphic and polymor-
phic experiments was the size of the evolving
populations. In the monomorphic limit, each
mutation is either lost or goes to fixation before
the next occurs. We enforced this evolutionary
dynamic by holding the population size to a
single protein sequence, similar to the “blind
ant” random walk of [11]. At each generation,
we assayed a single mutant. If this mutant met
the selection criterion, then it was carried over
to the next generation, corresponding to a neu-
tral mutation going to fixation. If the mutant
failed the selection criterion, then the popu-
lation stayed at the previous sequence for the
next generation, corresponding to a mutation
lost to selection. If all of the mutants assayed
had zero or one mutations, then this proto-
col would correspond exactly to the equations
of [11, 22]. However, in order to achieve appre-
ciable sequence evolution on a laboratory time
scale, we used a mutation rate that sometimes
produced multiple mutations in a generation.
We mathematically describe this situation in
the Mathematical Appendix; here we simply
note that it is possible to think of each gener-
ation as introducing a single mutational event
rather than a single mutation. We performed
22 independent replicates of monomorphic evo-
lution, evolving each for 25 generations.
In the polymorphic limit, the population
spreads across many sequences. To implement
this experimentally, we assayed 435 mutants at
each generation. The selection criterion was
used to classify each mutant as functional or
nonfunctional. In neutral evolution, all func-
tional mutants reproduce with equal proba-
bility. We therefore pooled equal volumes of
stationary-phase cultures of each functional mu-
tant and recovered the pooled genes with a mini-
prep. The polymorphic evolution experiment
therefore approaches the equations of [11, 22],
again with the exception that a sequence may
undergo multiple mutations at a single gener-
ation. We give the equations describing this
situation in the Mathematical Appendix. Since
the population evolves deterministically in the
polymorphic limit [11, 22], a single replicate
was performed. Because mutations accumu-
late more rapidly in the polymorphic experi-
ments than the monomorphic ones, we evolved
the polymorphic population for 15 generations
rather than 25.
Mutations and mutational robustness
Figure 3 shows how mutations accumulated dur-
ing the course of the neutral evolution experi-
ments (full data are in Table 2 and Additional
File 2). Since the unselected protein popula-
tions evolve without constraint, mutations ac-
cumulate at the same rate at which they are
introduced by error-prone PCR, 1.4 nucleotide
mutations per generation. Because selection
eliminates mutations that disrupt P450 activ-
ity, mutations accumulate more slowly in the
monomorphic and polymorphic populations.
Mutations accumulate more rapidly in the poly-
morphic population than in the monomorphic
populations. This difference in rates is pre-
dicted by the equations in the Mathematical
Appendix to be a consequence of the fact that
the polymorphic population is more mutation-
ally robust, and so can tolerate more of the
possible mutations.
To test directly whether the polymorphic
population evolves higher average mutational
robustness, we measured the fraction of 435
random mutants that met the selection crite-
rion. Figure 4 shows that the polymorphic pop-
ulation neutrally evolved to a markedly higher
mutational robustness than the monomorphic
populations, with 50± 2% of the final polymor-
phic population mutants continuing to function
versus 39 ± 2% for the final monomorphic pop-
ulations (Chi-square P -value of 10−3 that these
values are significantly different). The only dif-
ference between the two types of populations
was their size, so evolution has clearly favored
mutational robustness in the larger and thus
more polymorphic population. This finding
represents the first experimental support for
the prediction that highly polymorphic popu-
lations evolve excess mutational robustness [11].
Theory predicts that the excess mutational
robustness of a highly polymorphic protein
population comes from increased protein sta-
bility [22]. Because the P450 variants unfold
irreversibly, an equilibrium thermodynamic sta-
bility ∆Gf cannot be measured. We therefore
determined stability to irreversible thermal and
chemical denaturation, two highly correlated
measures of P450 stability that have previously
been shown to contribute to mutational ro-
bustness [5] (see Additional Files 3, 4, and 5).
Figure 5 shows that proteins from the polymor-
phic population were in fact more stable than
their counterparts from the monomorphic pop-
ulation. We also observed that proteins in the
polymorphic population tended to accumulate
to higher levels in E. coli (Figure 5). Elevated
expression could be a byproduct of increased
stability, or it could independently increase mu-
tational robustness by allowing the proteins to
better tolerate mutations that decrease codon
adaptation or reduce folding efficiency. It is pos-
sible that additional unrecognized biophysical
factors also contributed to the excess mutational
robustness of the polymorphic population, but
no such factors were immediately obvious.
Interpretation in terms of the P450 neutral
network
The higher mutational robustness of the poly-
morphic population is due to the fact that it
occupies the P450 gene neutral network dif-
ferently than the monomorphic populations.
Measurements from the evolution experiments
can therefore be used to infer basic proper-
ties of the underlying neutral network of P450
genes, as originally noted by van Nimwegen and
coworkers [11]. In the Mathematical Appendix,
we derive approximations for the normalized
principal eigenvalue 〈ν〉∞ and the normalized
average connectivity 〈ν〉o of the neutral net-
work, where in both cases the normalization
is obtained by dividing by the network coordi-
nation number. We obtain 〈ν〉∞ = 0.51 and
〈ν〉o = 0.35 for the P450 gene neutral net-
work. Our ability to consistently estimate these
two parameters from four different experimental
measurements supports the idea that the the-
ory that we elaborate in the Mathematical Ap-
pendix appropriately describes the experiments.
The difference between 〈ν〉∞ and 〈ν〉o is a mea-
sure of the extent to which some P450 neutral
network nodes have more connections than oth-
ers. We note that 〈ν〉∞ is approximately equal
to the exponential decline parameter for the
asymptotic decline in the fraction of functional
mutants with increasing numbers of random nu-
cleotide mutations [3, 26,27] (see Mathematical
Appendix). Previous studies looking at this ex-
ponential decline have reported 〈ν〉∞ = 0.7 for
subtilisin [26], 〈ν〉∞ = 0.7 for 3-methyladenine
DNA glycosylase [27], and 〈ν〉∞ = 0.7 - 0.8
for TEM1 β-lactamase [3]. These comparisons
suggest that P450 has a sparser neutral net-
work (smaller 〈ν〉∞) than these other proteins.
We suspect, however, that these earlier studies
(one of which is our own) overestimate 〈ν〉∞
due to insufficient equilibration of the starting
sequence. We believe that the approach of the
current work is more accurate for determining
〈ν〉∞ because the measurements are made after
many mutations have equilibrated the initial se-
quence. This approach could be used in future
experiments to compare the neutral network
connectivities of proteins from different fami-
lies.
Conclusions
We have demonstrated that neutral evolution
favors more mutationally robust proteins when
the evolving population is highly polymorphic.
Strikingly, the excess mutational robustness is
due only to population polymorphism, and so
will arise in any population of sufficiently large
size. Our work is the first experimental demon-
stration of this phenomenon, which is predicted
to occur quite generally in neutrally evolving
proteins and nucleic acids [11]. Furthermore,
we were able to identify one of the biophysical
factors underlying the increase in mutational
robustness by showing that proteins from the
highly polymorphic population are more sta-
ble. We recognize that evolution in a biological
context will be more complex. In our exper-
iments, fitness was the P450’s ability to be
expressed in active form by bacteria grown to
saturation in an environment with plentiful nu-
trients. Biological fitness, however, depends on
numerous additional and subtle effects such as
the metabolic costs of synthesis or the burdens
imposed by misfolded molecules. Some muta-
tions that are neutral in the experiments may
therefore have deleterious effects in a biologi-
cal setting [28]. The experiments nonetheless
capture the overriding constraint that proteins
retain their biochemical functions. Our success
in quantitatively explaining the results supports
the notion that important aspects of protein
evolution can be described simply in terms of
mutational effects on stability [22,28].
An obvious question is whether evolution
in nature favors mutational robustness by the
process we have demonstrated. Whether natu-
ral populations will neutrally evolve mutational
robustness depends on whether they are suffi-
ciently polymorphic, which will be the case if
the product of their effective population size N
and per protein per generation mutation rate µ
is much greater than one [11, 12]. Accurately
estimating Nµ, which is closely related to the
widely used parameter θ in population genet-
ics, for natural populations is difficult [29, 30]
(note that since mutational robustness is a
protein-wide property, the relevant mutation
rate is per protein, which is ≈ 102 to 103 larger
than the per codon mutation rate). For hu-
mans and other multicellular organisms, Nµ
is probably too small [31] for their proteins to
neutrally evolve mutational robustness. But
estimates [31, 32] place Nµ ≈ 10 to 100 for
typical-length proteins in bacteria, and it is
probably much higher for many viruses [33,34].
It is therefore likely that many viral and some
bacterial proteins have neutrally evolved extra
mutational robustness.
The neutral evolution of protein mutational
robustness may also contribute to adaptive evo-
lution. Experiments have shown that extra sta-
bility increases a protein’s evolvability by allow-
ing it to tolerate a wider range of functionally
beneficial but destabilizing mutations [5]. A
similar phenomenon seems to occur in natural
evolution, where functionally neutral but stabi-
lizing mutations can play a key role in adaptive
evolution by counterbalancing the destabilizing
effects of other functionally beneficial muta-
tions [35]. Viruses and perhaps bacteria may
thus benefit from large population sizes and
high mutation rates that drive an increase in
the mutational robustness and stability of their
proteins, which in turn enhances the capacity of
these proteins to rapidly change their sequences
and evolve new functions.
Methods
Equilibration evolution of the P450 protein
We began with a 21B3 P450 peroxygenase that
had been engineered for highly efficient hydrox-
ylation of 12-pNCA [23] (sequence shown in Ad-
ditional File 6). This P450 was not well equi-
librated to the constant selection criterion that
we planned to impose, since it had substantially
higher total activity. We therefore neutrally
evolved it for 16 generations in order to cre-
ate P450s that were better equilibrated to the
selection criterion. We evolved two parallel pop-
ulations, which we named R1 and R2. The pro-
cedure was exactly identical to that described
below for the polymorphic evolution with the
following exceptions:
• Starting sequence: the starting sequence
for the equilibration evolution was the
21B3 sequence.
• Population size: each of the two equili-
bration evolution populations had a size
of 174 sequences rather than the 435 used
for the polymorphic evolution.
• Selection criterion: the sequences were re-
quired to have at least 75% of the total
activity of the 21B3 P450.
• Mutation rate: the mutation rate for the
equilibration evolution was much higher
than for the polymorphic evolution. The
error-prone PCR protocol used 200 µM
manganese chloride (MnCl2), rather than
the 25 µM used for the polymorphic evo-
lution. We estimate that this error-prone
PCR protocol introduced ≈ 4 nucleotide
mutations per P450 gene at each genera-
tion during the equilibration evolution.
We performed 16 generations of equilibration
evolution, and then randomly selected 23 func-
tional mutants from each of the R1 and R2
populations (Additional File 7). We picked one
of these mutants, R1-11, for use as the parent
for the neutral evolution experiments.
Detailed protocol for evolution experiments
We began with the R1-11 P450 BM3 heme
domain variant (sequence in Additional File
1) cloned into the pCWori [25] plasmid with
a 5’ BamH1 and 3’ EcoR1 site as described
in [5]. The cloning primers were pCWori for (5’-
GAAACAGGATCCATCGATGCTTAGGAGGTCAT-
3’ and pCWori rev clone (5’-GCTCATGTTTGACAGCTTATCATCG-
3’). We used error-prone PCR to generate
mutants, taking great care to make the error-
prone PCR protocol repeatable by using a rela-
tively small number of thermal cycles. This was
both to control the mutation rate by ensuring
that the reaction did not saturate the reagents
(which would cause the number of doublings to
become sensitive to the initial template concen-
tration), and to avoid the PCR-based recombi-
nation events which can occur during with the
last few thermal cycles of PCR reactions [36,37].
The PCR reactions were 100 µl in volume, and
contained ≈ 13 ng of plasmid template (cor-
responding to ≈ 3 ng of template gene), 7
mM magnesium chloride MgCl2, 1 × Applied
BioSystems PCR Buffer II without MgCl2, 25
µM MnCl2, 0.5 µM pCWori for primer, 0.5
µM pCWori rev primer, 200 µM of dATP and
dGTP, 500 µm of dTTP and dCTP, and 5 units
of Applied Biosystems AmpliTaq polymerase.
The reactions were run on the BLOCK setting
of a MJ Research PCR machine with a program
of 95oC for 2 minutes, then 15 cycles of (95oC
for 30 seconds, 57oC for 30 seconds, 72oC for 90
seconds), and then cooling to 4oC. This proto-
col yielded roughly 1-1.5 µg of product gene (as
quantified by gel electrophoresis versus a known
standard), for a PCR efficiency of ≈ 0.5. Se-
quencing the unselected populations at the end
of the experiment indicated that this protocol
introduced an average of 1.4 ± 0.2 nucleotide
mutations, with the nucleotide error-spectrum
shown in Table 1. Because the number of PCR
doublings is large compared the average muta-
tion rate, the distribution of mutations among
sequences should be well-described by the Pois-
son distribution [38,39].
The mutant genes from the error-prone PCR
were purified over a ZymoResearch DNA clean
and concentrator column, and digested at 37oC
with EcoR1 and BamH1. The digested genes
were then purified from an agarose gel with
ZymoResearch DNA gel extraction columns,
and ligated into pCWori plasmid that had been
digested with BamH1 and EcoR1 and dephos-
phorylated. The ligations were transformed
into electro-competent catalase-free Escherichia
coli [25] (the catalase is removed because it
breaks down the hydrogen peroxide utilized by
the P450 peroxygenase), plated on Luria Broth
(LB) plates containing 100 µg/ml of ampicillin
to select for the plasmid’s antibiotic resistance
marker, and grown at 37oC. Transformation of
a control ligation reaction without any digested
gene yielded at least 100-fold fewer colonies,
indicating that the rate of plasmid self-ligation
was less than one percent.
Individual mutant colonies from the plates
were picked into 96-well 2 ml deep-well plates
containing 400 µl of LB supplemented with 100
µg/ml ampicillin. Each plate contained four
parental control wells with cells carrying the
parent R1-11 gene, four null control wells with
cells carrying the pCWori plasmid without a
P450 gene, and a non-inoculated well to check
for contamination. For the polymorphic pop-
ulation, we picked five such plates with all 87
other wells containing different mutants for a
total population size of 5 × 87 = 435 mutants.
For the 22 monomorphic populations (we be-
gan with 24 populations but two had to be
discarded due to contamination), we picked a
single colony for growth and screening. For the
unselected populations we picked a single colony
for growth without screening. The LB deep-well
plates were grown for 16-20 hours at 30oC, 210
revolutions per minute (rpm), and 80% relative
humidity in a Kuhner humidified shaker. To ex-
press the P450 mutants, we prepared 2 ml deep
well plates containing 400 µl per well of terrific
broth (TB) supplemented with 200 µM iso-
propyl β-D-thiogalactoside (IPTG), 100 µg/ml
ampicillin, and 500 µM of δ-aminolevulinic acid.
We used a pipetting robot inoculated these TB
plates with 100 µl from the LB plates. We
stored the LB deep-well plates at 4oC, and
grew the TB deep-well plates in the humidi-
fied shaker at 30oC, 210 rpm, and 80% relative
humidity for 22-24 hours. After this growth,
the cells were harvested by centrifuging the
TB deep-well plates at 4000×g for 5 minutes
and discarding the liquid. The cell pellets were
flash-frozen in liquid nitrogen to aid in cell lysis.
To lyse the cells for the assays, we resus-
pended the cell pellets in 300 µl of 100 mM [4-
(2-hydroxyethyl)-1-piperazinepropanesulfonic
acid] (EPPS) (pH 8.2) with 0.5 mg/ml lysozyme
and 4 units/ml of deoxyribonuclease by pipet-
ting 40 times with the pipetting robot. We then
incubated the plates at 37oC for 30 minutes,
again resuspended with the pipetting robot, and
put back at 37oC for an addition 30 minutes.
We then pelleted the cell debris by centrifu-
gation at 6000×g for 5 minutes at 4oC. The
pipetting robot was used to dispense 80 µl of
the clarified lysate into 96-well microtiter plates
(Rainin). We prepared a 6× stock of 1.5 mM
12-pNCA in 36% dimethyl sulfoxide (DMSO)
and the EPPS buffer (the 12-pNCA was stored
in the DMSO solution and combined with the
buffer immediately before use). We used a mul-
tichannel pipette to add 20 µl of this substrate
stock to each well of the microtiter plate. We
briefly mixed the plates with “shake” setting
of a 96-well plate spectrophotometer, and read
an absorbance baseline at 398 nm. We then
immediately added 20 µl of a freshly prepared
solution of 24 mM hydrogen peroxide in the
EPPS buffer to initiate the reaction, and mixed
again. The final reaction conditions were there-
fore the EPPS buffer with 6% DMSO, 4 mM
hydrogen peroxide, and 250 µM 12-pNCA. After
40 minutes we quantified the amount of enzy-
matic product by the increase in absorbance
at 398 nm. This absorbance increase is due
to the 4-nitrophenolate molecule released after
the P450 hydroxylates the twelfth carbon of the
12-pNCA molecule [23]. To score the mutants
as functional or nonfunctional, we compared
their gain in absorbance minus the median null
control reading to that of the median parental
control reading minus the median null control
reading. All mutants that had at least 75% of
the parental gain were scored as functional, all
other mutants were scored as nonfunctional.
We used the information from these assays
to select the parents for the next generation. For
the unselected population we did not require the
mutants to be functional, so the selected mutant
was used to start a 4 ml culture of LB with 100
µg/ml ampicillin, and the plasmid DNA was
harvested with a mini-prep. This plasmid DNA
was used as the template for the next round of
error-prone PCR. Therefore, after the first gen-
eration the four unselected replicates diverged
into four separate error-prone PCR reactions.
These unselected replicates were evolved for a
total of twelve generations, and were sequenced
at every third generation.
For the polymorphic population, all mu-
tants that were functional contributed an equal
amount of plasmid DNA as template for the
next generation. In order to do this, we col-
lected 50 µl of the culture from the LB deep-well
plate for each mutant that was scored as func-
tional. All of these LB aliquots were pooled,
and then the plasmid DNA was collected with
a mini-prep. The pool of plasmid DNA was
used as template for the next generation’s error-
prone PCR reactions. We performed 15 gener-
ations of evolution for this polymorphic pop-
ulation. Note that at each generation we are
assaying 435 mutants as part of the evolution-
ary procedure, so this provides information on
mutational robustness. At every third gen-
eration, we also selected a random sample of
functional mutants for sequencing. After 15
generations, we randomly selected 22 mutants
for stability measurements and sequencing anal-
ysis. The random selections were made from all
functional mutants with the Python computer
language random number generator.
For the monomorphic populations, at each
generation we assayed just a single mutant. If
that mutant was nonfunctional, then at that
generation the population stayed at its original
sequence. In that case, for the next generation
we simply picked a new mutant from the previ-
ous generation’s plate of transformed mutants.
If the mutant we screened was functional, then
that mutant represented the new population.
We therefore grew a 4 ml LB culture with 100
µg/ml of ampicillin, and collected the plasmid
DNA with a miniprep. That plasmid DNA was
then used as the template for the next genera-
tion’s error-prone PCR reaction. We thus had
22 (actually 24 before 2 were contaminated) in-
dependent monomorphic populations that were
being evolved in parallel. Each was evolved for
25 generations, and at the end of these 25 gen-
erations we measured the stability of the final
sequence of each population. Each time an as-
sayed mutant was functional, we sequenced the
new P450 gene. We also measured the average
mutational robustness of the monomorphic pop-
ulations at every fifth generation. To do this,
we did a pooled mini-prep of equal volumes
of LB cultures of all 22 replicates to obtain a
equal mix of plasmid DNA. We then performed
error-prone PCR on this mix, and assayed 435
mutants to measure the fraction functional. Full
neutral evolution data are in Additional File 2.
Test for recombination during error-prone
During the polymorphic population evolution,
we performed error-prone PCR on a mix of
different plasmids. It is common for PCR
on mixed templates to lead to recombination
events during the reaction [36, 37]. We at-
tempted to reduce this recombination by using
a small number of thermal cycles. However, in
order to test for recombination, we analyzed
the sequences of the final 22 selected members
of the polymorphic population. There are a
variety of statistical tests to detect recombi-
nation in a set of sequences. A comparison
of these tests by Posada [40] found that the
Max-Chi2 method developed by John Maynard
Smith [41] performs well. A publicly avail-
able implementation of this method [42] is at
http://www.lifesci.sussex.ac.uk/CSE/test/maxchi.php.
We used this implementation to analyze the 22
final polymorphic sequences, and the resulting
P -value was 0.29 after 100 random permuta-
tions, indicating that there is not significant
recombination.
Measurement of P450 stabilities
We measured the stabilities to both irreversible
thermal and irreversible urea denaturation
of the final (generation 25) member of each
monomorphic population, as well as of the 22
randomly selected members of the polymorphic
population. As discussed in the Supplemen-
tary Information of [5], cytochrome P450 BM3
heme domains (and indeed most P450s) dena-
ture irreversibly, forcing us to use resistance
to irreversible denaturation to quantify protein
stability. The first stability measure is the T50,
defined as the temperature at which half of the
protein is denatured after a 10 minute incu-
bation. The second stability measure is the
[urea]50, defined as the urea concentration at
which half of the protein denatures after a 4
hour room-temperature incubation. Each set
of measurements (those of T50 and [urea]50)
was performed on all of the mutants in the
same day, and each mutant was treated identi-
cally. Therefore, it is possible to make accurate
comparisons of the relative values of the mea-
surements within the data set. However, the
absolute values of the T50 and [urea]50 values
may be less accurate. Therefore, care should be
taken in comparing the absolute value of these
measurements to those of other studies (such
as [5]).
Both the T50 and [urea]50 measurements
were performed in clarified cell lysate. The
protein was expressed using catalase-free E.
coli [25] containing the encoding gene on the
IPTG inducible pCWori [25] plasmid. We used
freshly streaked cells to inoculate 2 ml cultures
of LB supplemented with 100 µg/ml of ampi-
cillin, and grew these starter cultures overnight
with shaking at 37oC. We then used 0.5 ml
http://www.lifesci.sussex.ac.uk/CSE/test/maxchi.php
from these starter cultures to inoculate 1 L
flasks containing 200 ml of TB supplemented
with 100 µg/ml of ampicillin. The TB cultures
were grown at 30oC and 210 rpm until they
reached an optical density at 600 nm of ≈0.9,
at which point IPTG and δ-aminolevulinic acid
were added to a final concentration of 0.5 mM
each. The cultures were grown for an additional
19 hours, then the cells were harvested by pel-
letting 50 ml aliquots at 5,500 g and 4oC for 10
min, and stored at -20oC. To obtain clarified
lysate, each pellet was resuspended in 8 ml of
100 mM EPPS (pH 8.2) and lysed by sonica-
tion, while being kept on ice. The cell debris
was pelleted by centrifugation at 8,000 g and
4oC for 10 minutes, and the clarified lysate was
decanted and kept on ice.
For the T50 measurements, 125 µl of clar-
ified lysate from a single mutant was added
to all 12 wells in a row of a 96-well hard-shell
thin-wall microplate (MJ Research). The plate
was heated for 10 minutes using the gradient
method of an Eppendorf Mastercycler gradient
PCR machine, with the gradient set at either
33oC-45oC or 46oC-58oC (each mutant was ex-
posed to both of these gradients), the machine
on the BLOCK setting, and the heated lid set to
75oC with the lid WAIT option. The plate was
then cooled to 4oC, removed from the PCR ma-
chine, and centrifuged at 5,500 g and 4oC for 5
minutes to pellet any debris. A pipetting robot
was used to dispense 80 µl of the supernatent
into a 96-well microtiter plate (Rainin), and the
amount of remaining properly folded P450 was
quantified from the carbon monoxide difference
spectrum as described below. The T50 values
were determined by fitting sigmoidal curves the
percent of remaining protein as shown in Addi-
tional File 3. Our ability to accurately compare
T50 values within the data set requires that each
well in a given column of the gradient PCR ma-
chine be at the same temperature. We used
a thermocouple to measure the temperature of
the wells with the machine lid open, and con-
firmed that the wells were within a few tenths
of a degree of the same temperature. Further
evidence for the consistency of our T50 val-
ues comes from the fact that two independent
measurements of the T50 for our R1-11 par-
ent yielded values that differed by only 0.1oC.
However, the absolute values of the measured
temperatures are less accurate. Thermocouple
measurements indicated that, with the machine
lid open, the wells were ≈ 1oC cooler than the
indicated temperature. We were unable to as-
certain the temperatures with the heated lid
closed, but based on comparisons water bath
measurements, the temperatures with the lid
closed slightly exceeded the indicated tempera-
tures.
For the [urea]50 measurements, 125 µl of the
clarified lysate from a single mutant was added
to all 12 wells in a row of a 96-well microtiter
plate. A pipetting robot was then used to add
and mix 125 µl of a 2X solution of urea in 100
mM EPPS (pH 8.2) so that each subsequent
column had a higher concentration of urea, and
so that the final urea concentrations were those
shown in Additional File 4. The plates were
left on the bench at room temperature for 4
hours, and the amount of remaining properly
folded P450 was quantified from the carbon
monoxide difference spectrum as described be-
low. The [urea]50 values were determined by
fitting sigmoidal curves to the percent of re-
maining protein. Evidence for the consistency
of the [urea]50 measurements comes from the
fact that two independent measurements of the
[urea]50 for our R1-11 parent yielded values that
differed by only 0.01 M. In addition, the [urea]50
and T50 values are highly correlated (Additional
File 5), indicating that they provide consistent
measures of stability.
For both the T50 and [urea]50 measurements,
the folded P450 was quantified from the car-
bon monoxide difference spectrum [43]. The
microtiter plates containing the P450 samples
were first used to read blank spectra at 450 and
490 nm using a Tecan Safire 2 plate reader. The
plates were then incubated for 10 minutes in an
airtight oven with carbon monoxide. The plates
were removed form the oven and 10 µl of 0.1 M
sodium hydrosulfite in 1.3 M potassium phos-
phate (pH 8.0) was immediately added to each
well. After 5-10 minutes, spectra were again
read at 450 and 490 nm. The amount of P450
is proportional to the increase in the signal at
450 nm after this procedure minus the change
in the signal at 490 nm.
Mathematical Appendix
Contents
A.1 Mathematical background 12
A.2 Monomorphic limit 14
A.3 Polymorphic limit 15
A.4 Approximations for polymorphic limit 16
A.5 Approximations for monomorphic limit 18
A.6 Interpretation in terms of neutral networks 20
A.7 Detailed justification for approximating pM by po 22
A.1 Mathematical background
The first purpose of this appendix is to provide mathematical equations that describe the experi-
ments. The second is to show how four measurements from the experiments can be used to calculate
two quantities that describe the topology of the underlying protein neutral network. We will derive
two equations for both quantitites, each in terms of a different measurement. The fact that the four
equations will be seen to yield consistent results provides evidence for the accuracy of the following
calculations.
Our calculations are based on a view of neutral protein evolution as a process constrained by a
stability threshold, a view that we originally introduced to explain experimental protein mutagen-
esis results [3]. The calculations closely parallel our earlier work [22], which is in turn based on a
general theoretical treatment of evolution on neutral networks by van Nimwegen and coworkers [11].
These calculations will probably be most thoroughly understood by first reading those works. The
primary difference between the current calculations and [22] is that previously we assumed that
the per generation per protein mutation rate µ was ≪ 1, so that at each generation a protein was
either unmutated (with probability 1−µ) or experienced a single mutation (with probability µ). In
contrast, here we allow the mutation rate to be arbitrarily large, so that a protein may experience
multiple mutations in a single generation (in this sense the calculations resemble the generalization
by Wilke [18] of [11]). Specifically, let fm be the probability that a protein experiences m mutations
in a single generation. Here we derive results for arbitrary fm, and then approximations relevant to
the form of fm in the experiments. In the limiting case of small mutation rate (where f0 = 1− µ,
f1 = µ, and fm = 0 for m > 1), the calculations here reduce to those in [22]. Proteins evolving
in nature typically experience very low mutation rates, so [22] probably offers the best description
of natural protein evolution. The calculations presented here are designed to specifically treat the
evolutionary dynamics of the experiments.
A protein’s thermodynamic stability is described by its free energy of folding, ∆Gf , with more
negative values indicating more stable proteins. As described in several previous papers [3, 5, 22],
we assume that selection requires a protein to fold with some minimal stability ∆Gminf , so that a
protein adequately folds if and only if ∆Gf ≤ ∆G
f . The amount of extra stability a protein
possesses relative to the stability threshold is given by ∆Gextraf = ∆Gf − ∆G
f ; note that all
folded proteins will have ∆Gextraf ≤ 0. We further assume that as long as ∆G
extra
f ≤ 0, selection is
indifferent to the exact amount of extra stability that a protein possesses (see [22] for a discussion
of the limitations of this assumption). We conceptually divide the continuous variable of protein
stability into small discrete bins of width b. Specifically, a protein is in bin i if it has ∆Gextraf be-
tween (1− i) b and −ib, where i = 1, 2, . . .. Mutating a protein changes its stability by an amount
∆∆G (defined as the stability of the mutant protein minus the stability of the initial protein), and
so may move it to a new stability bin. In [22], we defined a matrix W with elements Wij giving
the transition probabilities that a single mutation changes a protein’s stability from bin j to bin i.
We noted that W could be computed from the distribution of ∆∆G values for all single mutations,
and argued that W remains fairly constant during neutral evolution since the distribution of ∆∆G
values remains relatively unchanged. However, we emphasize that (as discussed in detail in [22])
the constancy of the ∆∆G distribution remains an assumption, albeit one that has now been shown
to be quite accurate for lattice proteins [3, 22,44] and provide a consistent theoretical explanation
for a growing body of experimental results (the current work as well as [3]).
Since we are allowing for larger mutation rates, and we must consider the possibility that a
protein’s stability might change due to multiple mutations at a single generation. Therefore, we
make a more general definition of Wij,m as the probability that m random mutations to a protein
in stability bin j change its stability to bin i, and let Wm be the matrix with elements Wij,m. Note
that Wm only describes mutations that cause transitions from one folded protein to another, since
the stability bins i = 1, 2, . . . all correspond to folded proteins. As before [22], we assume that Wm
is roughly constant during evolution, meaning that the distribution of ∆∆G values for multiple
mutations is roughly constant during neutral evolution. Note that if m = 1, then Wm is just the
matrix W that can be computed from the distribution of single-mutant ∆∆G values [22]. We will
now use the matrices Wm to calculate the following characteristics of a population that has evolved
to equilibrium: the distribution of stabilities, the average number of mutations 〈m〉T accumulated
after T generations, and the average fraction 〈F〉 of stably folded proteins in the population. We
then introduce a few approximations (that should be quite accurate for the experimental work in
this paper) that greatly simplify these calculations. Finally, we relate the calculations to properties
of the underlying protein neutral network.
As described generally by van Nimwegen and coworkers [11], the evolutionary dynamics depend
on whether the evolving population tends to be monomorphic or highly polymorphic. When the per
sequence per generation mutation rate µ is ≪ 1, whether the population is mostly monomorphic or
highly polymorphic is determined by the product of the population size N and µ: when Nµ ≪ 1
the population is mostly monomorphic, and when Nµ ≫ 1 the population is highly polymorphic.
However, with multiple mutations per generation, Nµ is no longer an appropriate parameter to
distinguish between mono- and polymorphism, since if the population size is sufficiently small the
population can still be monomorphic even if there are multiple mutations per generation. Specifi-
cally, in one set of experiments we constrained the population to be monomorphic (by maintaining
a population size of one), but still allowed the single protein in this population to experience more
than one mutation at a generation. So we instead denote the populations as either monomorphic or
polymorphic. We indicate quantities calculated for the monomorphic population by the subscript
M (i.e. 〈F〉M ) and those calculated for the polymorphic population by the subscript P (i.e. 〈F〉P ).
A.2 Monomorphic limit
In the limit of a completely monomorphic population, all of the proteins are in a single stability
bin. Let pi (t) be the probability that the population is in stability bin i at time t, and let p (t)
be the column vector with elements pi (t). At each generation there is a probability f0 that there
is no mutation that becomes fixed in the population, a probability of
Wij,mpj that the
population experiences a mutational event (which could be a single mutation or several simultaneous
mutations) that moves it into bin i, and a probability
Wji,m that the population is in
bin i and experiences one or more mutations that move it to another bin of stably folded proteins.
Define νi,m =
Wji,m to be the fraction of m-mutants of a protein in bin i that still fold, and let
Vm be the matrix with diagonal elements given by Vii,m = νi,m and all other elements zero. The
time evolution of p is
p (t+ 1) =
fm (Wm −Vm)
p (t) (1)
where I is the identity matrix. Note that mutations that destabilize a protein beyond the stabil-
ity threshold are immediately lost to natural selection, and so leave the population in its original
stability bin. This describes the experiments for the monomorphic populations, where we retain
the parental sequence if the single mutant we generate is nonfunctional. Equation 1 corresponds to
Equation (1) of [22], and the blind ant random walk described by van Nimwegen and coworkers [11].
Equation 1 describes a Markov process with a non-negative, irreducible, and acyclic transition
matrix, and so p approaches a unique stationary distribution (equilibrium value) of pM given by
the eigenvector equation
fm (Wm −Vm)
pM. (2)
Once p has reached equilibrium, the average fraction of proteins that still stably fold at each
generation is
〈F〉M = e
pM (3)
where e = (1, . . . , 1) is the unit row vector.
To calculate 〈m〉T,M , the average number of mutations accumulated after T generations once the
population has equilibrated, we note that at each generation there is a probability of fmpj
Wij,m
that a randomly chosen protein is in bin j, experiences m mutations, and still stably folds. The
average number of mutations accumulated in a single generation is simply the average ofm weighted
over this probability. So summing over all values of m and j, we see that
〈m〉T,M = Te
mfmWmpM. (4)
This equation corresponds to Equation (6) of [22], which was derived using an embedded Markov
process formalism. Here we have foregone this formalism for the more intuitive argument presented
above, since we do not attempt to calculate higher moments of the number of mutations.
A.3 Polymorphic limit
In the limit when the population is highly polymorphic, at each generation there are sequences in
many different stability bins. In this case, we describe the distribution of stabilities by the column
vector x (t), with element xi (t) giving the fraction of proteins in stability bin i at time t. At
generation t, the fraction of mutants that continue to fold is
〈F〉t = e
x (t) . (5)
Therefore, in order to maintain a constant population size, each remaining protein must produce
an average of αt = 〈F〉t
offspring. The population therefore evolves according to
x (t+ 1) = αt
x (t) . (6)
After the population evolves for a sufficiently long period of time, x will approach an equilibrium
value of xP. At this equilibrium, the average fraction of mutants that fold at each generation is
〈F〉P = e
xP, (7)
and the equilibrium reproduction rate is α = 〈F〉P
. Therefore,
xP = α
xP. (8)
Equations 7 and 8 can be combined to show that xP and 〈F〉P can be calculated from the eigenvector
equation
(〈F〉P − f0)xP =
fmWmxP, (9)
with (〈F〉P − f0) the principal eigenvalue of the nonnegative and irreducible matrix
fmWm.
Equation 9 corresponds to Equation (14) of [22], Equation (6) of the work by van Nimwegen and
coworkers [11], and Equation (13) of the work by Wilke [18].
We now calculate 〈m〉T,P , the average number of mutations accumulated in T generations after
the population has equilibrated. At equilibrium, there is a probability of fmxj
Wij,m that a
protein is in bin j, experiences m mutations, and still stably folds. Subsequently, all of these folded
proteins produce an average of α offspring. The average number of mutations accumulated in a
single generation is simply the average of m weighted over this probability, and then multiplied by
the average reproduction rate. So summing over all values of m and j, we obtain
〈m〉T,P = αTe
mfmWmxP =
mfmWmxP. (10)
This equation is the counterpart of Equation (18) of [22], where we have again foregone the em-
bedded Markov process formalism for a more intuitive derivation.
A.4 Approximations for polymorphic limit
We can dramatically simplify the results from the previous sections with several reasonable ap-
proximations. The first approximation is that the ∆∆G values for random mutations are roughly
additive, and is supported by a number of experimental studies of the thermodynamic effects of
mutations [45–47]. We have previously shown that this approximation can be used to accurately
describe experimental protein mutagenesis data with a simple stability threshold model [3]. Under
this approximation, the distribution of net ∆∆G values for multiple mutations can be computed
from the distribution of ∆∆G values for single mutations by performing convolutions of the single-
mutation ∆∆G distribution [3], meaning that Wm for arbitrary m can be computed solely from the
distribution of ∆∆G values for single mutations. However, to simplify the equations from previous
sections, we need to express Wm for arbitrary m only in terms of W (recall that W = W1). Since
W only contains information about stability transitions from folded proteins to other folded pro-
teins, if we make the second approximation that a protein that is destabilized beyond the minimal
stability threshold by one mutation is not re-stabilized to a folded protein by a subsequent muta-
tion, then Wm = W
m. This approximation that unfolded proteins are not re-stabilized should be
quite accurate since stabilizing mutations tend to be relatively rare and small in magnitude [48–51]
(this is the underlying idea behind the Markov chain approximation that was shown to be highly
accurate for lattice proteins [44]). To summarize, if ∆∆G values are roughly additive and stabilizing
mutations are rare, we have the approximation
Wm ≈ W
m. (11)
Simplifying the equations of the previous sections also requires assigning a specific functional
form to fm, the probability that a sequence undergoesmmutations. Here we assume that mutations
are Poisson distributed among sequences, so that
e−µµm
where µ =
mfm is the average number of mutations per protein per generation. When the mu-
tations are introduced by error-prone PCR, the Poisson distribution is an excellent approximation
to the true theoretical distribution of mutations created by error-prone PCR [38,39] provided that
µ is much less than the number of PCR doublings, as is the case in all of the experiments in the
current work.
We now use the approximations of Equations 11 and 12 to simplify the results given above for
the highly polymorphic limit. We begin by using these approximations to rewrite Equation 9 as
〈F〉P − e
xP = e
WmxP. (13)
This equation makes clear that xP is the principal eigenvector of the matrix
Wm, therefore
xP must also be the principal eigenvector of W. Now in our earlier work [22], we defined the
principal eigenvector of W as x∞, called the corresponding eigenvalue 〈ν〉∞, and showed that this
eigenvalue is shown the average fraction of single mutations that are neutral in a population that
is evolving with Nµ ≫ 1 and µ ≪ 1. Therefore, with the approximation of Equation 11, xP and
x∞ are equal, and are both defined by the same eigenvector equation,
〈ν〉∞xP = WxP = Wx∞ = 〈ν〉∞x∞. (14)
Combining Equations 13 and 14 we have,
〈F〉PxP = e
(µ〈ν〉∞)
1−〈ν〉∞
xP (15)
Equation 15 can be solved to yield
〈ν〉∞ = 1 +
ln 〈F〉P
. (16)
Similarly, we can simplify Equation 10,
〈m〉T,P =
mfmWmxP
1−〈ν〉∞
eWmxP
= Te−µ〈ν〉∞
(µ〈ν〉∞)
= Tµ〈ν〉∞e
−µ〈ν〉∞
(µ〈ν〉∞)
= Tµ〈ν〉∞. (17)
Solving this equation for 〈ν〉∞ yields
〈ν〉∞ =
〈m〉T,P
. (18)
A.5 Approximations for monomorphic limit
We now simplify the equations for the monomorphic limit. This requires several further approxima-
tions. We begin by approximating that the stability probability distribution pM given by Equation
2 by the distribution po defined in [22] as satisfying
0 = (W −V)po. (19)
The basic rationale behind approximating pM with po is that Equation 2 can be viewed as a per-
turbation to Equation 19 [52]. Essentially, po is an eigenvector of the matrix W −V while pM is
the corresponding eigenvector of the matrix W − V +
(Wm −Vm). The latter matrix
can be viewed as a perturbation to the first, since the sum
(Wm −Vm) is small. This
smallness is due to the fact that Wm tends to zero with large m, causing Vm to tend towards the
identity matrix. In addition, the µm/m! terms tend to zero with large m. Therefore, the terms in
the summation are all simply either a perturbation to W − V or involve subtracting terms that
are fractions of the identity matrix. The perturbations lead to bounded changes in the eigenvec-
tors [52], while the identity matrix terms do not change the eigenvectors. Below we give a more
rigorous justification of the assumption that pM is approximately equal to po.
We need one additional approximation to make further progress. Both Equations 3 and 4 con-
tain terms of the form Wmpo, and even if we use Equation 11 to rewrite these terms as W
there are no further clear simplifications. However, any probability vector that is multiplied re-
peatedly by W and normalized will eventually converge to x∞ = xP (since this is the principal
eigenvector of W). We make the approximation that this convergence is sufficiently rapid to be
essentially complete after a single multiplication. This approximation is supported by both protein
mutagenesis studies [3,26,27] that indicate that proteins rapidly converge to an exponential decline
in the fraction folded (indicating the stability distribution has equilibrated, as discussed below, and
by lattice protein studies showing the same [3, 44]. Therefore, we make the approximation that
eWmpo = 〈ν〉oeW
m−1x∞ = 〈ν〉o〈ν〉∞
where 〈ν〉o = eWpo has the same definition as in [22],
where it was defined as the average fraction of functional single mutants of a population evolving
with µ ≪ 1 and Nµ ≪ 1.
We use these approximations to simplify Equation 3 as
〈F〉M = e
= e−µ
eWmpo
= e−µ
1 + µ〈ν〉o
(µ〈ν〉∞)
= e−µ
(µ〈ν〉∞)
= e−µ
eµ〈ν〉∞ − 1
. (20)
Solving this equation for 〈ν〉o, we find
〈ν〉o =
〈ν〉∞ (〈F〉Me
µ − 1)
eµ〈ν〉∞ − 1
. (21)
We now use the approximations to simplify Equation 4 as
〈m〉T,M = Te
mfmWmpM
= Te−µ
eWmpo
= Te−µ〈ν〉o
= µTe−µ〈ν〉o
(µ〈ν〉∞)
= µT 〈ν〉oe
〈ν〉∞−1
. (22)
Solving this equation for 〈ν〉o yields
〈ν〉o =
〈m〉T,Me
1−〈ν〉∞
. (23)
To recap, we now have equations to calculate 〈ν〉∞ and 〈ν〉o from experimentally measurable
quantities. Equations 16 and 18 allow us to calculate 〈ν〉∞ from 〈F〉P and 〈m〉T,P , respectively.
Given this calculated value of 〈ν〉∞, Equations 21 and 23 then allow us to calculate 〈ν〉o from 〈F〉M
and 〈m〉T,M , respectively. The fact that we have two equations each for 〈ν〉∞ and 〈ν〉o allows us
to assess the self-consistency of the approach.
A.6 Interpretation in terms of neutral networks
Throughout the preceding calculations, we have referred to 〈ν〉∞ and 〈ν〉o as we defined them
in [22]: namely, as the average neutrality of protein populations evolving with µ ≪ 1 and Nµ
either ≫ 1 or ≪ 1, respectively. However, van Nimwegen and coworkers [11] have shown that
they can also be interpreted in terms of the underlying neutral network. In the experiments we
make mutations at the nucleotide (rather than amino acid) level, so each point in our sequence
space corresponds to a different gene. Every gene that yields an amount of protein sufficient to
hydroxylate the twelfth carbon of 12-p-nitrophenoxydodecanoic acid with at least 75% of the total
activity conferred by the original R1-11 parent gene represents a node on this neutral network.
We note that in the experiments (and also usually in natural evolution), the edges on the neutral
network are not all completely equivalent or fully undirected, since some mutations are more likely
to occur than others (for example, error-prone PCR with Taq polymerase is more likely to cause an
A→G mutation than an A→C mutation). In the analysis that follows, we ignore this complication
and assume all neutral network edges are equivalent.
In an extremely insightful analysis, van Nimwegen and coworkers [11] have shown that impor-
tant characteristics of a neutral network can be inferred from evolutionary quantities. Specifically,
they have shown that if a population is evolving with µ ≪ 1 and Nµ ≫ 1, then the average
neutrality (which we have denoted by 〈ν〉∞) is equal to the principal eigenvalue of the adjacency
matrix of the neutral network, normalized by the network coordination number (number of possi-
ble connections per node). In addition, they pointed out that a population evolving with µ ≪ 1
and Nµ ≪ 1 moves like a blind ant random walk, meaning that the average neutrality (which we
have denoted by 〈ν〉o) is equal to the average connectivity of a neutral network node divided by
the network coordination number. In our P450 experiments, we have measured the values needed
to estimate 〈ν〉∞ and 〈ν〉o using Equations 16, 18, 21, and 23. Using the final values listed in
Table 2, 〈F〉P = 0.50 and 〈F〉M = 0.39. Taking the final nucleotide mutation values from Table
2, 〈m〉T,P /T = 0.69 and 〈m〉T,M/T = 0.31. The average mutation rate, computed from the un-
selected population, is µ = 1.40. So using Equation 16, 〈ν〉∞ = 0.53, while using Equation 18,
〈ν〉∞ = 0.49. The consistency of these two values supports the idea that the calculations above
accurately describe the evolutionary process. Taking the average value of these two measurement
as 〈ν〉∞ = 0.51, we can then use Equations 21 and 23 to calculate 〈ν〉o. We calculate values of
0.28 and 0.43, respectively. These estimates differ by more than those for 〈ν〉∞, perhaps because
additional approximations have gone into the derivation of the relevant equations (in addition, we
have made no attempt to carry out the rather complex propagation of the sampling errors of Table
2). However, the values are still in a similar range. Taking the average of these two values, we
estimate that 〈ν〉o = 0.35. So overall, we predict that each functional P450 gene should have an
average fraction of 0.35 of its sequence nearest neighbors also encoding a functional gene, for an
average of about 1,500 neighbor genes. We predict that the principal eigenvalue of the neutral
network adjacency matrix is 0.51 ×3L. The fact that principal eigenvalue exceeds the average
connectivity indicates that the neutral network is not a regular graph, but instead has some nodes
more highly connected than others.
The value for 〈ν〉∞ calculated above can also be related to measurements from protein mu-
tagenesis experiments. A number of studies [3, 26, 27] have observed that the probability that a
protein remains functional after m mutations falls off exponentially with the number of mutations.
In fact, the decline is not always exponential for the first few mutations if the starting protein has
especially high or low stability [3] or activity [53], but will still converge to this exponential form
after a few mutations [3,44,54]. The stability threshold model can be used to relate this decline to
〈ν〉∞, as is done indirectly in the Markov chain approximation of [44]. Here we make that connec-
tion explicit. The initial protein has a stability that falls into some stability bin i. Therefore, its
stability can be described by the column vector y0, which has element i equal to one and all other
elements equal to zero. Now imagine constructing all single mutants of this protein. The fraction
of these single mutants that still fold is just eWy0, and the distribution of stabilities among the
single mutants is y1 = Wy0 (note that the elements of y1 no longer sum to one). Similarly, after
m mutations, the fraction of mutants that still fold is eWmy0, and the distribution of stabilities
among the m-mutants is ym = Wmy0. With the approximation of Equation 11, ym = W
This makes it clear that ym will converge to a vector proportional to x∞, the principal eigenvector
of W. Once this convergence is complete, each new mutation simply reduces the fraction of mu-
tants that fold by a factor of 〈ν〉∞, the principal eigenvalue of W (and the spectral radius of the
neutral network normalized by the coordination number). Therefore, what we have called 〈ν〉∞ in
the present work and [22] is equal to what is called x in [27], q in [26], and 〈ν〉 in [3]. The major
difficulty that is faced in extracting 〈ν〉∞ by the method of those three studies [3, 26,27] is that it
is not possible to directly assay mutants with m mutations, but instead only possible to assay a set
of mutants with a distribution of m. All three studies use different (and valid) methods to account
for this distribution, but this accounting is still difficult because most of the functional mutants
come from the low m end of the distribution. This makes it hard to get accurate values for the
fraction functional after large numbers of mutations, since most of the functional mutants in the
set come from sequences with few mutations. For this reason, we believe the current method of
measuring 〈ν〉∞ is more accurate. A second caution about comparing values of 〈ν〉∞ from different
studies is that its value depends on the nucleotide error-spectrum of the experiment, since differ-
ent mutagenesis methods create different distributions of nucleotide and amino acid mutation types.
We also briefly mention how we arrived at an estimate of 〈ν〉∞ for 3-methyladenine DNA
glycosylase from the data of [27]. This paper reports that a fraction x = 0.34 of amino acid
mutations inactivate the protein. We would like to determine the fraction 〈ν〉∞ of nucleotide mu-
tations that do not inactivate the protein. Roughly 75% of random mutations to a gene will be
synonymous. Therefore, m amino acid mutations should cause about 4m/3 nucleotide mutations.
The study of [27] measures that after m mutations, a fraction (1− x)
of the mutants are func-
tional. That means that 〈ν〉∞
fraction should be functional. Equating these expressions yields
〈ν〉∞ = exp
log (1− x)
. So using x = 0.34, we arrive at 〈ν〉∞ = 0.73.
A.7 Detailed justification for approximating pM by po
Here we provide a detailed justification for the approximation that pM is about equal to po. In the
monomorphic limit, the time evolution of p is given by Equation 1, and the stationary distribution
pM is given by Equation 2. We assume the approximations of Equations 11 and 12 and show that
we can approximate pM by po, where po is given by Equation 19. To justify this approximation,
we insert po into the right hand side of Equation 1 and ask to what extent po is left unaltered by
the dynamics. If po is found to be stationary to good approximation then, by uniqueness of the
stationary distribution of an ergodic process, po would be a good approximation to pM.
We therefore suppose that at some time t the distribution is given by po and compute, using
Equation 1, the change in po after one generation. The new distribution at time t+ 1 is given by
p (t+ 1) =
fm (W
m −Vm)
po. (24)
Using (V −W)po = 0, and taking components of the above equation, we obtain
pi (t+ 1) = p0i +
fm [(W
m −Vm)po]i . (25)
Thus po would be an approximately stationary distribution of the dynamics if |
fm [(W
m −Vm)po]i| ≪
p0i. We now proceed to show that this will be the case in most situations of interest by deriving
upper and lower bounds on the second term of the right hand side of Equation 25.
Consider first the term (Wmpo)i, which can be written as
(Wmpo)i =
k1,...,km
Wik1Wk1k2 · · ·Wkm−1kmp0km
k1,...,km−1
Wik1Wk1k2 · · ·Wkm−2km−1νkm−1p0km−1 , (26)
where we have used Wpo = Vpo in the second equality. We now note that νk ≤ νmax for all
k, where νmax is the maximum neutrality, maximized over all bins. This leads to the successive
inequalities
(Wmpo)i ≤ νmax
k1,...,km−1
Wik1Wk1k2 · · ·Wkm−2km−1p0km−1
= νmax
k1,...,km−2
Wik1Wk1k2 · · ·Wkm−3km−2νkm−2p0km−2
k1,...,km−2
Wik1Wk1k2 · · ·Wkm−3km−2p0km−2
≤ νm−1
Wik1p0k1 , (27)
yielding the upper bound
(Wmpo)i ≤ ν
νip0i. (28)
In an identical manner, we obtain the lower bound
(Wmpo)i ≥ ν
νip0i, (29)
where νmin is the smallest neutrality, minimized over all bins. Note that both inequalities above
become exact equalities when all bins have the same neutrality ν, which could be interpreted as
either νmin or νmax.
Having obtained inequality constraints on (Wmpo)i, we now consider the term (Vmpo)i, which
can be written as
(Vmpo)i = p0iνi,m
= p0i
(Wm)ji
= p0i
j,k1,...,km−1
Wjk1Wk1k2 · · ·Wkm−1i
= p0i
k1,...,km−1
νk1Wk1k2 · · ·Wkm−1i
≤ p0iνmax
k1,...,km−1
Wk1k2 · · ·Wkm−1i
≤ p0iν
Wkm−1i, (30)
which yields an identical upper bound to that on (Wmpo)i, namely,
(Vmpo)i ≤ ν
νip0i, (31)
and similarly
(Vmpo)i ≥ ν
νip0i. (32)
It should again be noted that both the above inequalities become exact equalities when all bins
have a common neutrality ν.
We are now in a position to estimate bounds on the magnitude of the second term of Equation
25. Using the four inequalities of Equations 28, 29, 31, and 32 above, we have
− νm−1
νip0i ≤ [(W
m −Vm)po]i ≤
− νm−1
νip0i, (33)
or equivalently,
|[(Wm −Vm)po]i| ≤
− νm−1
νip0i, (34)
where the inequality above becomes an exact equality when all bins have the same neutrality. How-
ever, in this limit, the right hand side of the above equation vanishes, and therefore the second
term of Equation 25 is identically zero in this case, giving the result that pM is exactly equal to
po when all bins have the same neutrality, even if µ is arbitrarily large.
We now carry out the sum over m to obtain an upper bound on the second term of Equation 25
in the more general and realistic case of unequal neutrality bins. Using Equation 34 and the specific
Poisson form of fm, we obtain an upper bound on the fractional change in p0i in one generation:
pi(t+ 1)− p0i
≤ νie
− νm−1
= νie
eµνmax − 1
eµνmin − 1
. (35)
The above bound vanishes for small µ, is an increasing function of νmax − νmin, and is typically
much smaller than 1. An extreme estimate of the size of the fractional change can be made when
νmax = 1 and νmin = 0. In this case, using µ = 1.4 (the value in our experiments), the above
inequality simplifies to
pi(t+ 1)− p0i
1− e−µ − µe−µ
≃ 0.41νi. (36)
Noting that νi < 1, the fractional change in p0i is therefore reasonably controlled even in the most
extreme case. For realistic situations, the fractional change in p0i is expected to be much lower,
thus justifying the use of po as the stationary distribution of the dynamics of Equation 1.
Author Contributions
JDB and FHA designed the project and wrote the paper. JDB and ZL performed the bulk of the
experiments; OSV assisted with the experiments. JDB and DC analyzed the data. JDB and AR
performed the theoretical work.
Acknowledgments
We thank Claus O Wilke for helpful advice and comments. JDB is supported by a HHMI predoctoral
fellowship. ZL and DC were supported by Summer Undergraduate Research Fellowships from the California
Institute of Technology. AR is supported by NSF grants CCF 0523643 and FIBR 0527023.
References
1. Zuckerkandl E, Pauling L: Evolutionary divergence and convergence in proteins. In Evolving genes and
proteins, New York, NY: Academic Press 1965:97–166.
2. Lesk AM, Chothia C: How different amino acid sequences determine similar protein structures: The
structure and evolutionary dynamics of the globins. J. Mol. Biol. 1980, 136:225–270.
3. Bloom JD, Silberg JJ, Wilke CO, Drummond DA, Adami C, Arnold FH: Thermodynamic prediction of
protein neutrality. Proc. Natl. Acad. Sci. USA 2005, 102:606–611.
4. Besenmatter W, Kast P, Hilvert D:Relative tolerance of mesostable and thermostable protein homologs
to extensive mutation. Proteins 2007, 66:500–506.
5. Bloom JD, Labthavikul ST, Otey CR, Arnold FH: Protein stability promotes evolvability. Proc. Natl. Acad.
Sci. USA 2006, 103:5869–5874.
6. Lenski RE, Barrick JE, Ofria C: Balancing robustness and evolvability. PLoS Biology 2006, 4:e428.
7. Sniegowski PD, Murphy HA: Evolvability. Current Biology 2006, 16:R831–R834.
8. Wagner A: Robustness and Evolvability in Living Systems. Princeton, New Jersey: Princeton University Press
2005.
9. Montville R, Froissart R, Remold SK, Tenaillon O, Turner PE: Evolution of mutational robustness in an
RNA virus. PLoS Biol. 2005, 3:e381.
10. Wilke CO, Wang JL, Ofria C, Lenski RE, Adami C: Evolution of digital organisms at high mutation
rates leads to survival of the flattest. Nature 2001, 412:331–333.
11. van Nimwegen E, Crutchfield JP, Huynen M: Neutral evolution of mutational robustness. Proc. Natl. Acad.
Sci. USA 1999, 96:9716–9720.
12. Kimura M: The Neutral Theory of Molecular Evolution. Cambridge, U.K.: Cambridge Univ. Press 1983.
13. Takahata N: On the overdispersed molecular clock. Genetics 1987, 116:169–179.
14. Smith JM: Natural selection and the concept of a protein space. Nature 1970, 225:563–564.
15. Lipman DJ, Wilbur WJ: Modeling neutrality and selective evolution of protein folding. Proc. R. Soc.
London Ser. B 1991, 245:7–11.
16. Huynen MA, Stadler PF, Fontana W: Smoothness within ruggedness: the role of neutrality in adapta-
tion. Proc. Natl. Acad. Sci. USA 1996, 93:397–401.
17. Bornberg-Bauer E, Chan HS: Modeling evolutionary landscapes: mutational stability, topology, and
superfunnels in sequence space. Proc. Natl. Acad. Sci. USA 1999, 96:10689–10694.
18. Wilke CO: Adaptive evolution on neutral networks. Bull. Math. Biol. 2001, 63:715–730.
19. Taverna DM, Goldstein RA: Why are proteins so robust to site mutations? J. Mol. Biol. 2002, 315:479–
20. Xia Y, Levitt M: Funnel–like organization in sequence space determines the distributions of protein
stability and folding rate preferred by evolution. Proteins 2004, 55:107–114.
21. Brin S, Page L: The anatomy of a large-scale hypertextual Web search engine. Computer Networks and
ISDN Systems 1998, 30(1–7):107–117.
22. Bloom JD, Raval A, Wilke CO: Thermodynamics of neutral protein evolution. Genetics 2007, 175:255–
23. Cirino PC, Arnold FH: A self-sufficient peroxide-driven hydroxylation biocatalyst. Angew. Chem. Int.
Ed. 2003, 42:3299–3301.
24. Taverna DM, Goldstein RA: Why are proteins marginally stable? Proteins 2002, 46:105–109.
25. Barnes HJ, Arlotto MP, Waterman MR: Expression and enzymatic activity of recombinant cytochrome
P450 17 α-hydroxylase in Escherichia coli . Proc. Natl. Acad. Sci. USA 1991, 88:5597–5601.
26. Shafikhani S, Siegel RA, Ferrari E, Schellenberger V: Generation of large libraries of random mutants in
Bacillus subtilis by PCR-based plasmid multimerization. BioTechniques 1997, 23:304–310.
27. Guo HH, Choe J, Loeb LA: Protein tolerance to random amino acid change. Proc. Natl. Acad. Sci. USA
2004, 101:9205–9210.
28. DePristo MA, Weinreich DM, Hartl DL: Missense meanderings in sequence space: a biophysical view
of protein evolution. Nat. Rev. Genetics 2005, 6:678–687.
29. Plotkin JB, Dushoff J, Desai MM, Fraser HB: Codon usage and selection in proteins. J. Mol. Evol. 2006,
63:635–653.
30. Berg OG: Selection intensity for codon bias and the effective population size of Escherichia coli .
Genetics 1996, 142:1379–1382.
31. Lynch M, Conery JS: The origins of genome complexity. Science 2003, 302:1401–1404.
32. Hartl DL, Moriyama EN, Sawyer SA: Selection intensity for codon bias. Genetics 1994, 138:227–234.
33. Moya A, Holmes EC, Gonzalez-Candelas F: The population genetics and evolutionary epidemiology of
RNA viruses. Nature Reviews Microbiology 2004, 2:279–288.
34. Vignuzzi M, Stone JK, Arnold JJ, Cameron CE, Andino R:Quasispecies diversity determines pathogenesis
through cooperative interactions in a viral population. Nature 2006, 439:344–348.
35. Wang X, Minasov G, Shoichet BK: Evolution of an antibiotic resistance enzyme constrained by stability
and activity trade–offs. J. Mol. Biol. 2002, 320:85–95.
36. Meyerhans A, Vartanian JP, Wain-Hobson S: DNA recombination during PCR. Nucleic Acids Res. 1990,
18:1687–1691.
37. Kanagawa T: Bias and artifacts in multitemplate polymerase chain reactions (PCR). Journal of
Bioscience and Bioengineering 2003, 96:317–323.
38. Sun F: The polymerase chain reaction and branching processes. J. Comput. Biol. 1995, 2:63–86.
39. Drummond DA, Iverson BL, Georgiou F, Arnold FH: Why high–error–rate random mutagenesis libraries
are highly enriched in functional and improved proteins. J. Mol. Biol. 2005, 350:806–816.
40. Posada D: Evaluation of methods for detecting recombination from DNA sequences: empirical data.
Mol. Biol. Evol. 2002, 19:708–717.
41. Smith JM: Analyzing the mosaic structure of genes. J. Mol. Evol. 1992, 34:126–129.
42. Piganeau G, Gardner M, Eyre-Walker A: A broad survey of recombination in animal mitochondria. Mol.
Biol. Evol. 2004, 21:2319–2325.
43. Otey CR: High-throughput carbon monoxide binding assay for cytochromes P450. In Directed enzyme
evolution: screening and selection methods, Volume 230 of Methods in Molecular Biology. Edited by Arnold FH,
Georgiou G, Humana press 2003:137–139.
44. Wilke CO, Bloom JD, Drummond DA, Raval A: Predicting the tolerance of proteins to random amino
acid substitution. Biophysical J. 2005, 89:3714–3720.
45. Wells JA: Additivity of mutational effects in proteins. Biochemistry 1990, 29:8509–8517.
46. Zhang XJ, Baase WA, Shoichet BK, Wilson KP, Matthews BW: Enhancement of protein stability by the
combination of point mutations in T4 lysozyme is additive. Protein Eng. 1995, 8:1017–1022.
47. Serrano L, Day AG, Fersht AR: Step-wise mutation of barnase to binase: a procedure for engineering
increased stability of proteins and an experimental analysis of the evolution of protein stability.
J. Mol. Biol. 1993, 233:305–312.
48. Godoy-Ruiz R, Perez-Jimenez R, Ibarra-Molero B, Sanchez-Ruiz JM: Relation between protein stability,
evolution and structure as probed by carboxylic acid mutations. J. Mol. Biol. 2004, 336:313–318.
49. Pakula AA, Sauer RT: Genetic analysis of protein stability and function. Annu. Rev. Genet. 1989, 23:289–
50. Matthews BW: Structural and genetic analysis of protein stability. Annu. Rev. Biochem. 1993, 62:139–
51. Bava KA, Gromiha MM, Uedaira H, Kitajimi K, Sarai A: ProTherm, version 4.0: thermodynamic
database for proteins and mutants. Nucleic Acids Res. 2004, 32:D120–D121.
52. Franklin JN: Matrix Theory. Englewood Cliffs, N.J.: Prentice-Hall, Inc. 1968.
53. Bershtein S, Segal M, Bekerman R, Tokuriki N, Tawfik DS: Robustness-epistasis link shapes the fitness
landscape of a randomly drifting protein. Nature 2006, 444:929–932.
54. Bloom JD, Arnold FH, Wilke CO: Breaking proteins with mutations: threads and thresholds in evo-
lution. Molecular Systems Biology 2007, 3:76.
Figures
monomorphic polymorphic
BPSfrag replacements
0.2probability
threshold
more stable← ∆Gf → less stable
PSfrag replacements
threshold
more stable← ∆Gf → less stable
monomorphic
polymorphic
Figure 1 - Theoretical views of the evolution of protein mutational robustness and stability.
(A) Theory predicts that a mostly monomorphic population is equally likely to occupy any node of its neutral
network, while a highly polymorphic population will prefer more connected nodes [11]. Node sizes are drawn
proportional to the occupation probabilities. (B) Proteins evolving in a highly polymorphic population
are predicted to be more stable than their counterparts in a mostly monomorphic population [22]. The
histograms illustrate the distributions of stabilities for the two cases. The increased stability is a biophysical
manifestation of excess mutational robustness, since more stable proteins are more mutationally robust [3–5].
polymorphic monomorphic unselected
Figure 2 - Outline of the neutral evolution experimental procedure.
For the polymorphic population, error-prone PCR was used to generate mutant P450 genes. These genes
were ligated into a plasmid and transformed into E. coli. Individual mutants (435) were picked, expressed
in E. coli, and assayed for enzymatic activity. All mutants that met the selection criterion contributed an
equal amount of plasmid DNA as template for the next generation of error-prone PCR. The monomorphic
populations were treated similarly, except only a single mutant was assayed at each generation. If this mutant
met the selection criterion then it became the template for the next generation of error-prone PCR; otherwise
at the next generation another colony was picked from the same plate. In the unselected populations a
single mutant was picked and used as the template for the next generation of error-prone PCR.
PSfrag replacements
generation
〈mnt〉
unselected
polymorphic
monomorphic
PSfrag replacements
0 10 20
generation
〈maa〉 〈
unselected
polymorphic
monomorphic
Figure 3 - Accumulation of nucleotide (〈mnt〉) and nonsynonymous (〈maa〉) mutations in the
experimentally evolved P450 populations.
For the unselected and monomorphic populations, numbers are the average over all replicates at the indicated
generation; for the polymorphic population they are from a random sample, with sampling standard error
shown.
PSfrag replacements
0 10 20
generation
polymorphic
monomorphic
Figure 4 - The polymorphic population neutrally evolved a higher average mutational robustness
than the monomorphic populations.
The fraction functional was determined by assaying 435 mutants (average of 1.5 nucleotide mutations per
gene). Error bars show binomial standard error. For the monomorphic population, numbers are the average
over all replicates.
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
38 40 42 44
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
number
0.4 0.6 0.8 1.0
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
number
20 60 100 140
Figure 5 - The more mutationally robust proteins are more stable.
The P450s from the polymorphic population neutrally evolved higher stability and expression levels than
their counterparts from the monomorphic populations. The histograms show the distributions for the final
protein from all monomorphic replicates and for the same number of randomly chosen proteins from the final
polymorphic population. The plots show (left to right) the temperature at which half the protein irreversibly
denatured after 10 minutes (T50), the urea concentration at which half the protein denatured after 4 hours
([urea]50), and the expression level relative to that of the original parental P450. The means are significantly
different, with unequal variance t-test P -values of 0.02, 0.005, and 0.04, respectively.
Tables
Table 1 - Error-prone PCR nucleotide mutation spectrum.
Spectrum of nucleotide mutations introduced by the error-prone PCR procedure used in the neutral evolution
experiments. The spectrum was determined by sequencing the four final (generation 12) sequences from the
unselected population, since in these sequences the mutations accumulate without constraint. As has been
previously noted for error-prone PCR with Taq polymerase [3,5,26], the nucleotide error spectrum is biased
towards certain types of mutations.
Total nucleotide mutations 67
% synonymous mutations 25
Mutation types (%)
A →T, T →A 19.4
A →C, T →G 1.5
A →G, T →C 64.2
G →A, C →T 4.5
G →C, C →G 0.0
G →T, C →A 1.5
frameshift 9.0
Table 2 - Neutral evolution robustness and mutation data.
Each row is for a different generation, T . Entries of NA indicate that no measurement was made. The 〈mnt〉
and 〈maa〉 are the average number of nucleotide mutations and nonsynonymous mutations, respectively.
Numbers in parentheses are total counts over the total samples. Subscripts indicate the population type:
U for unselected, P for polymorphic, and M for monomorphic. For the unselected and monomorphic
populations, numbers represent averages of all replicates. For the polymorphic population, numbers are for
a random sample of functional mutants. 〈F〉P and 〈F〉M are the fraction of functional mutants out of 435
assayed.
T 〈mnt〉U 〈maa〉U 〈mnt〉P 〈maa〉P 〈mnt〉M 〈maa〉M 〈F〉P 〈F〉M
0 0 0 0 0 0 0 0.48 (210 / 435) 0.48 (210 / 435)
1 NA NA NA NA 0.1 (3 / 22) 0.3 (6 / 22) 0.48 (208 / 435) NA
2 NA NA NA NA 0.4 (9 / 22) 0.8 (17 / 22) 0.49 (215 / 435) NA
3 5.0 (20 / 4) 3.5 (14 / 4) 2.7 (27 / 10) 1.4 (14 / 10) 1.0 (23 / 22) 0.4 (9 / 22) 0.49 (215 / 435) NA
4 NA NA NA NA 1.5 (32 / 22) 0.7 (15 / 22) 0.48 (208 / 435) NA
5 NA NA NA NA 2.2 (48 / 22) 1.1 (25 / 22) 0.45 (197 / 435) 0.43 (185 / 435)
6 9.8 (39 / 4) 7.5 (30 / 4) 5.5 (55 / 10) 2.1 (21 / 10) 2.6 (58 / 22) 1.4 (31 / 22) 0.46 (198 / 435) NA
7 NA NA NA NA 3.1 (69 / 22) 1.8 (39 / 22) 0.52 (227 / 435) NA
8 NA NA NA NA 3.4 (74 / 22) 1.8 (40 / 22) 0.46 (200 / 435) NA
9 13.0 (52 / 4) 10.3 (41 / 4) 6.7 (61 / 9) 3.1 (28 / 9) 3.7 (82 / 22) 2.1 (46 / 22) 0.47 (203 / 435) NA
10 NA NA NA NA 4.2 (92 / 22) 2.4 (52 / 22) 0.46 (199 / 435) 0.40 (175 / 435)
11 NA NA NA NA 4.6 (102 / 22) 2.5 (56 / 22) 0.48 (207 / 435) NA
12 16.8 (67 / 4) 12.5 (50 / 4) 7.8 (70 / 9) 3.3 (30 / 9) 4.9 (107 / 22) 2.6 (58 / 22) 0.52 (228 / 435) NA
13 NA NA NA NA 5.0 (110 / 22) 2.7 (60 / 22) 0.52 (227 / 435) NA
14 NA NA NA NA 5.3 (116 / 22) 2.9 (64 / 22) 0.50 (216 / 435) NA
15 NA NA 10.3 (227 / 22) 3.8 (83 / 22) 5.6 (123 / 22) 3.0 (67 / 22) 0.50 (219 / 435) 0.39 (171 / 435)
16 NA NA NA NA 5.8 (127 / 22) 3.0 (67 / 22) NA NA
17 NA NA NA NA 6.0 (133 / 22) 3.1 (69 / 22) NA NA
18 NA NA NA NA 6.3 (137 / 22) 3.2 (71 / 22) NA NA
19 NA NA NA NA 6.3 (138 / 22) 3.3 (72 / 22) NA NA
20 NA NA NA NA 6.6 (145 / 22) 3.4 (75 / 22) NA 0.37 (160 / 435)
21 NA NA NA NA 6.9 (152 / 22) 3.6 (79 / 22) NA NA
22 NA NA NA NA 7.1 (156 / 22) 3.7 (81 / 22) NA NA
23 NA NA NA NA 7.2 (158 / 22) 3.7 (81 / 22) NA NA
24 NA NA NA NA 7.3 (161 / 22) 3.8 (83 / 22) NA NA
25 NA NA NA NA 7.7 (169 / 22) 4.0 (87 / 22) NA 0.39 (169 / 435)
Additional material
Additional file 1 - Sequence of the parent P450 used start neutral evolution.
FASTA file with sequence of the R1-11 P450 BME used as the neutral evolution parent. This sequence was
isolated after the equilibration evolution.
Additional file 2 - Information about sequences from neutral evolution experiments.
The entries give the name of the mutant, the number of nonsynonymous and nucleotide mutations relative
to the R1-11 parent, the [urea]50 value if measured, the T50 value if measured, the percent of the parental
expression level if measured, and then a list of all of the mutations. Amino acid mutations are numbered in
the standard P450 numbering scheme. The names of the mutants indicate their origin. Names beginning
with “P-G3” are randomly chosen functional mutants from generation 3 of the polymorphic population,
etc. Names of the form “P1,” “P2,”, etc. are the 22 functional mutants that were randomly chosen from
the final (generation 15) polymorphic population. Numbers P5 and P12 are missing because two of the
original 24 randomly selected polymorphic population members were randomly chosen to be discarded after
it was discovered that two of the 24 monomorphic replicates were contaminated. Names beginning with
“U1” indicate that sequences are from the first unselected replicate, etc. Names beginning “M1” indicate
sequences are from the first monomorphic replicate, etc. Replicates “M9” and “M10” were discarded due to
contamination during the experiment. For each replicate, we sequenced each new functional mutant. The
last functional mutant after 25 generations represents the final sequence for that replicate, and is given an
abbreviated name without the generation suffix.
Additional file 3 - Thermostability measurements.
Raw data from the T50 thermostability measurements.
Additional file 4 - Urea stability measurements.
Raw data from the [urea]50 thermostability measurements.
Additional file 5 - Correlation of thermal and urea stabilities.
The T50 and [urea]50 values are highly correlated.
Additional file 6 - Sequence of initial P450 used to start equilibration evolution.
FASTA file with sequence of the 21B3 P450 BM3 heme domain described in [23]. This P450 was used as
the initial parent to start the equilibration evolution.
Additional file 7 - Mutations accumulated during equilibration evolution.
The file lists the mutations in the 46 P450 variants selected at the end of the equilibration evolution. Each line
gives the name of the variant, with the prefix indicating whether it came from the R1 or R2 population. The
next entries give the number of nucleotide and nonsynonymous mutations. All of the individual mutations
relative to 21B3 are then listed. Amino acid mutations are numbered in the standard P450 numbering
scheme, with the threonine after the N-terminal methionine given the number one.
	A.1 Mathematical background
	A.2 Monomorphic limit
	A.3 Polymorphic limit
	A.4 Approximations for polymorphic limit
	A.5 Approximations for monomorphic limit
	A.6 Interpretation in terms of neutral networks
	A.7 Detailed justification for approximating pM by po

  BACKGROUND: An important question is whether evolution favors properties such
as mutational robustness or evolvability that do not directly benefit any
individual, but can influence the course of future evolution. Functionally
similar proteins can differ substantially in their robustness to mutations and
capacity to evolve new functions, but it has remained unclear whether any of
these differences might be due to evolutionary selection for these properties.
  RESULTS: Here we use laboratory experiments to demonstrate that evolution
favors protein mutational robustness if the evolving population is sufficiently
large. We neutrally evolve cytochrome P450 proteins under identical selection
pressures and mutation rates in populations of different sizes, and show that
proteins from the larger and thus more polymorphic population tend towards
higher mutational robustness. Proteins from the larger population also evolve
greater stability, a biophysical property that is known to enhance both
mutational robustness and evolvability. The excess mutational robustness and
stability is well described by existing mathematical theories, and can be
quantitatively related to the way that the proteins occupy their neutral
network.
  CONCLUSIONS: Our work is the first experimental demonstration of the general
tendency of evolution to favor mutational robustness and protein stability in
highly polymorphic populations. We suggest that this phenomenon may contribute
to the mutational robustness and evolvability of viruses and bacteria that
exist in large populations.

Evolution favors protein mutational robustness in sufficiently
large populations
Jesse D Bloom∗1 , Zhongyi Lu1 , David Chen1 , Alpan Raval2 , Ophelia S Venturelli1 and
Frances H Arnold∗1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA
Keck Graduate Institute of Applied Life Sciences and School of Mathematical Sciences, Claremont Graduate University,
Claremont, CA 91711, USA
Email: Jesse D Bloom∗- jesse.bloom@gmail.com; Zhongyi Lu - lu07@caltech.edu; David Chen - davidc@caltech.edu;
Alpan Raval - alpan raval@kgi.edu; Ophelia S Venturelli - opheliav@stanford.edu; Frances H Arnold∗-
frances@cheme.caltech.edu;
Corresponding author
Abstract
Background: An important question is whether evolution favors properties such as mutational robustness
or evolvability that do not directly benefit any individual, but can influence the course of future evolution.
Functionally similar proteins can differ substantially in their robustness to mutations and capacity to
evolve new functions, but it has remained unclear whether any of these differences might be due to
evolutionary selection for these properties.
Results: Here we use laboratory experiments to demonstrate that evolution favors protein mutational
robustness if the evolving population is sufficiently large. We neutrally evolve cytochrome P450 proteins
under identical selection pressures and mutation rates in populations of different sizes, and show that
proteins from the larger and thus more polymorphic population tend towards higher mutational robust-
ness. Proteins from the larger population also evolve greater stability, a biophysical property that is
known to enhance both mutational robustness and evolvability. The excess mutational robustness and
stability is well described by existing mathematical theories, and can be quantitatively related to the
way that the proteins occupy their neutral network.
Conclusions: Our work is the first experimental demonstration of the general tendency of evolution to
favor mutational robustness and protein stability in highly polymorphic populations. We suggest that
this phenomenon may contribute to the mutational robustness and evolvability of viruses and bacteria
that exist in large populations.
http://arxiv.org/abs/0704.1885v1
Background
Proteins are quite tolerant of mutations, al-
lowing evolution to produce highly diverged
sequences that fold to similar structures and
perform conserved biochemical functions [1, 2].
However, proteins with nearly identical struc-
tures and functions may differ in their robust-
ness to mutation [3–5], as well as in their ca-
pacity to acquire new functions [5]. The fact
that mutational robustness and evolvability can
vary among the functionally equivalent proteins
produced by natural sequence divergence makes
these properties important hidden dimensions
in evolution — direct selection for protein func-
tion is blind to them, yet they can play a cru-
cial role in enabling future evolution. Whether
the evolutionary process somehow promotes the
acquisition of mutational robustness and evolv-
ability therefore remains a major question [6–8].
Previous experiments have identified several
specific evolutionary conditions that can affect
mutational robustness. For example, genetic
complementation decreases the mutational ro-
bustness of viruses [9], while high mutation rates
favor mutational robustness in simulated digi-
tal organisms [10]. However, theory [11] makes
the much broader — and heretofore experimen-
tally untested — prediction that extra muta-
tional robustness will arise quite generally in
sufficiently large populations. This prediction
cannot be understood in the standard frame-
work of Kimura’s neutral theory [12], since one
of the usual assumptions of the neutral theory
is that mutational robustness is constant. (Al-
though Takahata [13] treated the consequences
of stochastically fluctuating neutrality on the
molecular clock, he did not describe how muta-
tional robustness might change systematically
during evolution.) However, changes in mu-
tational robustness can be described by envi-
sioning evolution as occurring on neutral net-
works, or sets of functionally equivalent pro-
teins that are connected by single mutational
steps [14–17]. In a seminal theoretical analysis
of evolution on neutral networks, van Nimwe-
gen and coworkers [11] predicted that the extent
of mutational robustness should depend on the
degree of population polymorphism. Here we
briefly summarize their reasoning, since it moti-
vates our experimental work. We also refer the
reader to chapter 16 of [8], which contains an
excellent explanation of the densely mathemat-
ical work of van Nimwegen and coworkers [11].
If an evolving population is mostly
monomorphic, then each mutation is either
lost or goes to fixation before another muta-
tion occurs. The population is therefore usually
clustered at a single genotype and rarely expe-
riences mutations, meaning that selection does
not distinguish between genotypes of different
mutational robustness. All nodes of the neutral
network are thus equivalent and will be occupied
by the population with equal probability [11].
On the other hand, a highly polymorphic pop-
ulation is always spread across many nodes of
the neutral network. When mutations occur,
the members of the population at highly con-
nected nodes have a better chance of surviving,
causing them to be favored by evolution and
increasing the average mutational robustness
[11, 17–20]. Specifically, a highly polymorphic
population occupies each node with a proba-
bility proportional to its eigenvector central-
ity [11, 17], a measure of how connected it is
to other connected nodes (a variant of eigen-
vector centrality is used by Google’s PageRank
algorithm to rank a webpage’s importance in
the network of internet links [21]). Figure 1A
illustrates how mostly monomorphic and highly
polymorphic populations are predicted to oc-
cupy a neutral network. For proteins, changes
in neutral network occupancy should be man-
ifested by changes in thermodynamic stabil-
ity [22], with proteins from highly polymorphic
populations predicted to be more stable than
their counterparts from mostly monomorphic
populations (Figure 1B). Note that the extent
of polymorphism depends on the product of
the mutation rate and population size, meaning
that protein populations of different sizes are
predicted to evolve to different levels of mu-
tational robustness and stability even if they
experience the same mutation rate.
Results and Discussion
Design of neutral evolution experiment
To test whether high population polymorphism
drives an increase in mutational robustness and
protein stability, we performed laboratory evo-
lution experiments on cytochrome P450 pro-
teins. The basic idea was to neutrally evolve
P450s under a constant selection pressure in
populations that were either monomorphic or
highly polymorphic, and observe whether the
proteins evolved to different levels of mutational
robustness and stability. The evolution experi-
ments started with a P450 BM3 heme domain
that had been engineered to hydroxylate 12-p-
nitrophenoxydodecanoic acid (12-pNCA) [23].
We imposed the selection criterion that Es-
cherichia coli cells expressing the P450 had to
yield lysate with enough active enzyme to hy-
droxylate a specified amount of 12-pNCA in 40
minutes. This criterion roughly corresponds to
the case in which an enzyme must catalyze a
biochemically relevant reaction at some mini-
mal level in order for its host to survive. Note
that other properties such as stability and ex-
pression level can vary freely, provided that the
criterion for total activity is met.
The properties of a neutrally evolving pro-
tein eventually “equilibrate,” much as the prop-
erties of an isolated physical system under some
macroscopic constraint tend towards the values
that maximize the system’s internal entropy.
For proteins, this usually means that stability,
expression, and activity drift towards their low-
est tolerable values, since the vast majority of
random sequences do not encode stable, well-
expressed enzymes (that is, natural selection
must work against sequence entropy to maintain
a functional protein) [22, 24]. The initial P450
had been engineered for maximal activity [23],
meaning that it was not equilibrated to the more
mild selection criterion of the experiments. We
therefore neutrally evolved this initial P450 for
16 generations, introducing random mutations
with error-prone PCR and retaining all mu-
tants that met the selection criterion for total
activity on 12-pNCA. The procedure used for
this equilibration evolution was similar to that
for the polymorphic neutral evolution described
below. As expected, expression, stability, and
activity all dropped during the equilibration
evolution. At the end of the equilibration evo-
lution, we chose a single sequence as the parent
for the neutral evolution experiments. The gene
encoding this parent sequence contained 29 nu-
cleotide mutations and 13 amino acid mutations
relative to the initial P450 (Additional File 1).
We used this parent gene to begin three
parallel sets of neutral evolution experiments,
which we named “monomorphic,” “polymor-
phic,” and “unselected” (Figure 2). The
monomorphic experiments capture the case
where the population moves as a single entity,
the polymorphic experiment captures the case
where the population spreads across many se-
quences, and the unselected experiments show
how the gene evolves in the absence of selection
for protein function. In all experiments, at each
generation we used error-prone PCR to intro-
duce an average of 1.4 nucleotide mutations per
P450 gene (Table 1). The mutant genes were
ligated into a plasmid and transformed into E.
coli [25], and transformants were selected using
the plasmid’s antibiotic resistance marker. For
the unselected case, we randomly picked one of
the mutants, recovered the mutant gene with
a plasmid mini-prep, and used this mutant as
the template for the next generation of error-
prone PCR. We performed four independent
replicates of unselected evolution, evolving each
for 12 generations.
For the monomorphic and polymorphic pop-
ulations, we imposed the selection criterion that
the P450s hydroxylate 12-pNCA with at least
75% of the total activity of the original par-
ent gene. We expressed the P450s in E. coli,
and then assayed the cell lysates for activity in
a high-throughput 96-well plate format. The
total amount of product produced by 80 µl of
clarified lysate in 40 minutes was compared to
the median of four control wells containing the
original parent P450 to determine if the mutant
met the selection criterion. The only differ-
ence between the monomorphic and polymor-
phic experiments was the size of the evolving
populations. In the monomorphic limit, each
mutation is either lost or goes to fixation before
the next occurs. We enforced this evolutionary
dynamic by holding the population size to a
single protein sequence, similar to the “blind
ant” random walk of [11]. At each generation,
we assayed a single mutant. If this mutant met
the selection criterion, then it was carried over
to the next generation, corresponding to a neu-
tral mutation going to fixation. If the mutant
failed the selection criterion, then the popu-
lation stayed at the previous sequence for the
next generation, corresponding to a mutation
lost to selection. If all of the mutants assayed
had zero or one mutations, then this proto-
col would correspond exactly to the equations
of [11, 22]. However, in order to achieve appre-
ciable sequence evolution on a laboratory time
scale, we used a mutation rate that sometimes
produced multiple mutations in a generation.
We mathematically describe this situation in
the Mathematical Appendix; here we simply
note that it is possible to think of each gener-
ation as introducing a single mutational event
rather than a single mutation. We performed
22 independent replicates of monomorphic evo-
lution, evolving each for 25 generations.
In the polymorphic limit, the population
spreads across many sequences. To implement
this experimentally, we assayed 435 mutants at
each generation. The selection criterion was
used to classify each mutant as functional or
nonfunctional. In neutral evolution, all func-
tional mutants reproduce with equal proba-
bility. We therefore pooled equal volumes of
stationary-phase cultures of each functional mu-
tant and recovered the pooled genes with a mini-
prep. The polymorphic evolution experiment
therefore approaches the equations of [11, 22],
again with the exception that a sequence may
undergo multiple mutations at a single gener-
ation. We give the equations describing this
situation in the Mathematical Appendix. Since
the population evolves deterministically in the
polymorphic limit [11, 22], a single replicate
was performed. Because mutations accumu-
late more rapidly in the polymorphic experi-
ments than the monomorphic ones, we evolved
the polymorphic population for 15 generations
rather than 25.
Mutations and mutational robustness
Figure 3 shows how mutations accumulated dur-
ing the course of the neutral evolution experi-
ments (full data are in Table 2 and Additional
File 2). Since the unselected protein popula-
tions evolve without constraint, mutations ac-
cumulate at the same rate at which they are
introduced by error-prone PCR, 1.4 nucleotide
mutations per generation. Because selection
eliminates mutations that disrupt P450 activ-
ity, mutations accumulate more slowly in the
monomorphic and polymorphic populations.
Mutations accumulate more rapidly in the poly-
morphic population than in the monomorphic
populations. This difference in rates is pre-
dicted by the equations in the Mathematical
Appendix to be a consequence of the fact that
the polymorphic population is more mutation-
ally robust, and so can tolerate more of the
possible mutations.
To test directly whether the polymorphic
population evolves higher average mutational
robustness, we measured the fraction of 435
random mutants that met the selection crite-
rion. Figure 4 shows that the polymorphic pop-
ulation neutrally evolved to a markedly higher
mutational robustness than the monomorphic
populations, with 50± 2% of the final polymor-
phic population mutants continuing to function
versus 39 ± 2% for the final monomorphic pop-
ulations (Chi-square P -value of 10−3 that these
values are significantly different). The only dif-
ference between the two types of populations
was their size, so evolution has clearly favored
mutational robustness in the larger and thus
more polymorphic population. This finding
represents the first experimental support for
the prediction that highly polymorphic popu-
lations evolve excess mutational robustness [11].
Theory predicts that the excess mutational
robustness of a highly polymorphic protein
population comes from increased protein sta-
bility [22]. Because the P450 variants unfold
irreversibly, an equilibrium thermodynamic sta-
bility ∆Gf cannot be measured. We therefore
determined stability to irreversible thermal and
chemical denaturation, two highly correlated
measures of P450 stability that have previously
been shown to contribute to mutational ro-
bustness [5] (see Additional Files 3, 4, and 5).
Figure 5 shows that proteins from the polymor-
phic population were in fact more stable than
their counterparts from the monomorphic pop-
ulation. We also observed that proteins in the
polymorphic population tended to accumulate
to higher levels in E. coli (Figure 5). Elevated
expression could be a byproduct of increased
stability, or it could independently increase mu-
tational robustness by allowing the proteins to
better tolerate mutations that decrease codon
adaptation or reduce folding efficiency. It is pos-
sible that additional unrecognized biophysical
factors also contributed to the excess mutational
robustness of the polymorphic population, but
no such factors were immediately obvious.
Interpretation in terms of the P450 neutral
network
The higher mutational robustness of the poly-
morphic population is due to the fact that it
occupies the P450 gene neutral network dif-
ferently than the monomorphic populations.
Measurements from the evolution experiments
can therefore be used to infer basic proper-
ties of the underlying neutral network of P450
genes, as originally noted by van Nimwegen and
coworkers [11]. In the Mathematical Appendix,
we derive approximations for the normalized
principal eigenvalue 〈ν〉∞ and the normalized
average connectivity 〈ν〉o of the neutral net-
work, where in both cases the normalization
is obtained by dividing by the network coordi-
nation number. We obtain 〈ν〉∞ = 0.51 and
〈ν〉o = 0.35 for the P450 gene neutral net-
work. Our ability to consistently estimate these
two parameters from four different experimental
measurements supports the idea that the the-
ory that we elaborate in the Mathematical Ap-
pendix appropriately describes the experiments.
The difference between 〈ν〉∞ and 〈ν〉o is a mea-
sure of the extent to which some P450 neutral
network nodes have more connections than oth-
ers. We note that 〈ν〉∞ is approximately equal
to the exponential decline parameter for the
asymptotic decline in the fraction of functional
mutants with increasing numbers of random nu-
cleotide mutations [3, 26,27] (see Mathematical
Appendix). Previous studies looking at this ex-
ponential decline have reported 〈ν〉∞ = 0.7 for
subtilisin [26], 〈ν〉∞ = 0.7 for 3-methyladenine
DNA glycosylase [27], and 〈ν〉∞ = 0.7 - 0.8
for TEM1 β-lactamase [3]. These comparisons
suggest that P450 has a sparser neutral net-
work (smaller 〈ν〉∞) than these other proteins.
We suspect, however, that these earlier studies
(one of which is our own) overestimate 〈ν〉∞
due to insufficient equilibration of the starting
sequence. We believe that the approach of the
current work is more accurate for determining
〈ν〉∞ because the measurements are made after
many mutations have equilibrated the initial se-
quence. This approach could be used in future
experiments to compare the neutral network
connectivities of proteins from different fami-
lies.
Conclusions
We have demonstrated that neutral evolution
favors more mutationally robust proteins when
the evolving population is highly polymorphic.
Strikingly, the excess mutational robustness is
due only to population polymorphism, and so
will arise in any population of sufficiently large
size. Our work is the first experimental demon-
stration of this phenomenon, which is predicted
to occur quite generally in neutrally evolving
proteins and nucleic acids [11]. Furthermore,
we were able to identify one of the biophysical
factors underlying the increase in mutational
robustness by showing that proteins from the
highly polymorphic population are more sta-
ble. We recognize that evolution in a biological
context will be more complex. In our exper-
iments, fitness was the P450’s ability to be
expressed in active form by bacteria grown to
saturation in an environment with plentiful nu-
trients. Biological fitness, however, depends on
numerous additional and subtle effects such as
the metabolic costs of synthesis or the burdens
imposed by misfolded molecules. Some muta-
tions that are neutral in the experiments may
therefore have deleterious effects in a biologi-
cal setting [28]. The experiments nonetheless
capture the overriding constraint that proteins
retain their biochemical functions. Our success
in quantitatively explaining the results supports
the notion that important aspects of protein
evolution can be described simply in terms of
mutational effects on stability [22,28].
An obvious question is whether evolution
in nature favors mutational robustness by the
process we have demonstrated. Whether natu-
ral populations will neutrally evolve mutational
robustness depends on whether they are suffi-
ciently polymorphic, which will be the case if
the product of their effective population size N
and per protein per generation mutation rate µ
is much greater than one [11, 12]. Accurately
estimating Nµ, which is closely related to the
widely used parameter θ in population genet-
ics, for natural populations is difficult [29, 30]
(note that since mutational robustness is a
protein-wide property, the relevant mutation
rate is per protein, which is ≈ 102 to 103 larger
than the per codon mutation rate). For hu-
mans and other multicellular organisms, Nµ
is probably too small [31] for their proteins to
neutrally evolve mutational robustness. But
estimates [31, 32] place Nµ ≈ 10 to 100 for
typical-length proteins in bacteria, and it is
probably much higher for many viruses [33,34].
It is therefore likely that many viral and some
bacterial proteins have neutrally evolved extra
mutational robustness.
The neutral evolution of protein mutational
robustness may also contribute to adaptive evo-
lution. Experiments have shown that extra sta-
bility increases a protein’s evolvability by allow-
ing it to tolerate a wider range of functionally
beneficial but destabilizing mutations [5]. A
similar phenomenon seems to occur in natural
evolution, where functionally neutral but stabi-
lizing mutations can play a key role in adaptive
evolution by counterbalancing the destabilizing
effects of other functionally beneficial muta-
tions [35]. Viruses and perhaps bacteria may
thus benefit from large population sizes and
high mutation rates that drive an increase in
the mutational robustness and stability of their
proteins, which in turn enhances the capacity of
these proteins to rapidly change their sequences
and evolve new functions.
Methods
Equilibration evolution of the P450 protein
We began with a 21B3 P450 peroxygenase that
had been engineered for highly efficient hydrox-
ylation of 12-pNCA [23] (sequence shown in Ad-
ditional File 6). This P450 was not well equi-
librated to the constant selection criterion that
we planned to impose, since it had substantially
higher total activity. We therefore neutrally
evolved it for 16 generations in order to cre-
ate P450s that were better equilibrated to the
selection criterion. We evolved two parallel pop-
ulations, which we named R1 and R2. The pro-
cedure was exactly identical to that described
below for the polymorphic evolution with the
following exceptions:
• Starting sequence: the starting sequence
for the equilibration evolution was the
21B3 sequence.
• Population size: each of the two equili-
bration evolution populations had a size
of 174 sequences rather than the 435 used
for the polymorphic evolution.
• Selection criterion: the sequences were re-
quired to have at least 75% of the total
activity of the 21B3 P450.
• Mutation rate: the mutation rate for the
equilibration evolution was much higher
than for the polymorphic evolution. The
error-prone PCR protocol used 200 µM
manganese chloride (MnCl2), rather than
the 25 µM used for the polymorphic evo-
lution. We estimate that this error-prone
PCR protocol introduced ≈ 4 nucleotide
mutations per P450 gene at each genera-
tion during the equilibration evolution.
We performed 16 generations of equilibration
evolution, and then randomly selected 23 func-
tional mutants from each of the R1 and R2
populations (Additional File 7). We picked one
of these mutants, R1-11, for use as the parent
for the neutral evolution experiments.
Detailed protocol for evolution experiments
We began with the R1-11 P450 BM3 heme
domain variant (sequence in Additional File
1) cloned into the pCWori [25] plasmid with
a 5’ BamH1 and 3’ EcoR1 site as described
in [5]. The cloning primers were pCWori for (5’-
GAAACAGGATCCATCGATGCTTAGGAGGTCAT-
3’ and pCWori rev clone (5’-GCTCATGTTTGACAGCTTATCATCG-
3’). We used error-prone PCR to generate
mutants, taking great care to make the error-
prone PCR protocol repeatable by using a rela-
tively small number of thermal cycles. This was
both to control the mutation rate by ensuring
that the reaction did not saturate the reagents
(which would cause the number of doublings to
become sensitive to the initial template concen-
tration), and to avoid the PCR-based recombi-
nation events which can occur during with the
last few thermal cycles of PCR reactions [36,37].
The PCR reactions were 100 µl in volume, and
contained ≈ 13 ng of plasmid template (cor-
responding to ≈ 3 ng of template gene), 7
mM magnesium chloride MgCl2, 1 × Applied
BioSystems PCR Buffer II without MgCl2, 25
µM MnCl2, 0.5 µM pCWori for primer, 0.5
µM pCWori rev primer, 200 µM of dATP and
dGTP, 500 µm of dTTP and dCTP, and 5 units
of Applied Biosystems AmpliTaq polymerase.
The reactions were run on the BLOCK setting
of a MJ Research PCR machine with a program
of 95oC for 2 minutes, then 15 cycles of (95oC
for 30 seconds, 57oC for 30 seconds, 72oC for 90
seconds), and then cooling to 4oC. This proto-
col yielded roughly 1-1.5 µg of product gene (as
quantified by gel electrophoresis versus a known
standard), for a PCR efficiency of ≈ 0.5. Se-
quencing the unselected populations at the end
of the experiment indicated that this protocol
introduced an average of 1.4 ± 0.2 nucleotide
mutations, with the nucleotide error-spectrum
shown in Table 1. Because the number of PCR
doublings is large compared the average muta-
tion rate, the distribution of mutations among
sequences should be well-described by the Pois-
son distribution [38,39].
The mutant genes from the error-prone PCR
were purified over a ZymoResearch DNA clean
and concentrator column, and digested at 37oC
with EcoR1 and BamH1. The digested genes
were then purified from an agarose gel with
ZymoResearch DNA gel extraction columns,
and ligated into pCWori plasmid that had been
digested with BamH1 and EcoR1 and dephos-
phorylated. The ligations were transformed
into electro-competent catalase-free Escherichia
coli [25] (the catalase is removed because it
breaks down the hydrogen peroxide utilized by
the P450 peroxygenase), plated on Luria Broth
(LB) plates containing 100 µg/ml of ampicillin
to select for the plasmid’s antibiotic resistance
marker, and grown at 37oC. Transformation of
a control ligation reaction without any digested
gene yielded at least 100-fold fewer colonies,
indicating that the rate of plasmid self-ligation
was less than one percent.
Individual mutant colonies from the plates
were picked into 96-well 2 ml deep-well plates
containing 400 µl of LB supplemented with 100
µg/ml ampicillin. Each plate contained four
parental control wells with cells carrying the
parent R1-11 gene, four null control wells with
cells carrying the pCWori plasmid without a
P450 gene, and a non-inoculated well to check
for contamination. For the polymorphic pop-
ulation, we picked five such plates with all 87
other wells containing different mutants for a
total population size of 5 × 87 = 435 mutants.
For the 22 monomorphic populations (we be-
gan with 24 populations but two had to be
discarded due to contamination), we picked a
single colony for growth and screening. For the
unselected populations we picked a single colony
for growth without screening. The LB deep-well
plates were grown for 16-20 hours at 30oC, 210
revolutions per minute (rpm), and 80% relative
humidity in a Kuhner humidified shaker. To ex-
press the P450 mutants, we prepared 2 ml deep
well plates containing 400 µl per well of terrific
broth (TB) supplemented with 200 µM iso-
propyl β-D-thiogalactoside (IPTG), 100 µg/ml
ampicillin, and 500 µM of δ-aminolevulinic acid.
We used a pipetting robot inoculated these TB
plates with 100 µl from the LB plates. We
stored the LB deep-well plates at 4oC, and
grew the TB deep-well plates in the humidi-
fied shaker at 30oC, 210 rpm, and 80% relative
humidity for 22-24 hours. After this growth,
the cells were harvested by centrifuging the
TB deep-well plates at 4000×g for 5 minutes
and discarding the liquid. The cell pellets were
flash-frozen in liquid nitrogen to aid in cell lysis.
To lyse the cells for the assays, we resus-
pended the cell pellets in 300 µl of 100 mM [4-
(2-hydroxyethyl)-1-piperazinepropanesulfonic
acid] (EPPS) (pH 8.2) with 0.5 mg/ml lysozyme
and 4 units/ml of deoxyribonuclease by pipet-
ting 40 times with the pipetting robot. We then
incubated the plates at 37oC for 30 minutes,
again resuspended with the pipetting robot, and
put back at 37oC for an addition 30 minutes.
We then pelleted the cell debris by centrifu-
gation at 6000×g for 5 minutes at 4oC. The
pipetting robot was used to dispense 80 µl of
the clarified lysate into 96-well microtiter plates
(Rainin). We prepared a 6× stock of 1.5 mM
12-pNCA in 36% dimethyl sulfoxide (DMSO)
and the EPPS buffer (the 12-pNCA was stored
in the DMSO solution and combined with the
buffer immediately before use). We used a mul-
tichannel pipette to add 20 µl of this substrate
stock to each well of the microtiter plate. We
briefly mixed the plates with “shake” setting
of a 96-well plate spectrophotometer, and read
an absorbance baseline at 398 nm. We then
immediately added 20 µl of a freshly prepared
solution of 24 mM hydrogen peroxide in the
EPPS buffer to initiate the reaction, and mixed
again. The final reaction conditions were there-
fore the EPPS buffer with 6% DMSO, 4 mM
hydrogen peroxide, and 250 µM 12-pNCA. After
40 minutes we quantified the amount of enzy-
matic product by the increase in absorbance
at 398 nm. This absorbance increase is due
to the 4-nitrophenolate molecule released after
the P450 hydroxylates the twelfth carbon of the
12-pNCA molecule [23]. To score the mutants
as functional or nonfunctional, we compared
their gain in absorbance minus the median null
control reading to that of the median parental
control reading minus the median null control
reading. All mutants that had at least 75% of
the parental gain were scored as functional, all
other mutants were scored as nonfunctional.
We used the information from these assays
to select the parents for the next generation. For
the unselected population we did not require the
mutants to be functional, so the selected mutant
was used to start a 4 ml culture of LB with 100
µg/ml ampicillin, and the plasmid DNA was
harvested with a mini-prep. This plasmid DNA
was used as the template for the next round of
error-prone PCR. Therefore, after the first gen-
eration the four unselected replicates diverged
into four separate error-prone PCR reactions.
These unselected replicates were evolved for a
total of twelve generations, and were sequenced
at every third generation.
For the polymorphic population, all mu-
tants that were functional contributed an equal
amount of plasmid DNA as template for the
next generation. In order to do this, we col-
lected 50 µl of the culture from the LB deep-well
plate for each mutant that was scored as func-
tional. All of these LB aliquots were pooled,
and then the plasmid DNA was collected with
a mini-prep. The pool of plasmid DNA was
used as template for the next generation’s error-
prone PCR reactions. We performed 15 gener-
ations of evolution for this polymorphic pop-
ulation. Note that at each generation we are
assaying 435 mutants as part of the evolution-
ary procedure, so this provides information on
mutational robustness. At every third gen-
eration, we also selected a random sample of
functional mutants for sequencing. After 15
generations, we randomly selected 22 mutants
for stability measurements and sequencing anal-
ysis. The random selections were made from all
functional mutants with the Python computer
language random number generator.
For the monomorphic populations, at each
generation we assayed just a single mutant. If
that mutant was nonfunctional, then at that
generation the population stayed at its original
sequence. In that case, for the next generation
we simply picked a new mutant from the previ-
ous generation’s plate of transformed mutants.
If the mutant we screened was functional, then
that mutant represented the new population.
We therefore grew a 4 ml LB culture with 100
µg/ml of ampicillin, and collected the plasmid
DNA with a miniprep. That plasmid DNA was
then used as the template for the next genera-
tion’s error-prone PCR reaction. We thus had
22 (actually 24 before 2 were contaminated) in-
dependent monomorphic populations that were
being evolved in parallel. Each was evolved for
25 generations, and at the end of these 25 gen-
erations we measured the stability of the final
sequence of each population. Each time an as-
sayed mutant was functional, we sequenced the
new P450 gene. We also measured the average
mutational robustness of the monomorphic pop-
ulations at every fifth generation. To do this,
we did a pooled mini-prep of equal volumes
of LB cultures of all 22 replicates to obtain a
equal mix of plasmid DNA. We then performed
error-prone PCR on this mix, and assayed 435
mutants to measure the fraction functional. Full
neutral evolution data are in Additional File 2.
Test for recombination during error-prone
During the polymorphic population evolution,
we performed error-prone PCR on a mix of
different plasmids. It is common for PCR
on mixed templates to lead to recombination
events during the reaction [36, 37]. We at-
tempted to reduce this recombination by using
a small number of thermal cycles. However, in
order to test for recombination, we analyzed
the sequences of the final 22 selected members
of the polymorphic population. There are a
variety of statistical tests to detect recombi-
nation in a set of sequences. A comparison
of these tests by Posada [40] found that the
Max-Chi2 method developed by John Maynard
Smith [41] performs well. A publicly avail-
able implementation of this method [42] is at
http://www.lifesci.sussex.ac.uk/CSE/test/maxchi.php.
We used this implementation to analyze the 22
final polymorphic sequences, and the resulting
P -value was 0.29 after 100 random permuta-
tions, indicating that there is not significant
recombination.
Measurement of P450 stabilities
We measured the stabilities to both irreversible
thermal and irreversible urea denaturation
of the final (generation 25) member of each
monomorphic population, as well as of the 22
randomly selected members of the polymorphic
population. As discussed in the Supplemen-
tary Information of [5], cytochrome P450 BM3
heme domains (and indeed most P450s) dena-
ture irreversibly, forcing us to use resistance
to irreversible denaturation to quantify protein
stability. The first stability measure is the T50,
defined as the temperature at which half of the
protein is denatured after a 10 minute incu-
bation. The second stability measure is the
[urea]50, defined as the urea concentration at
which half of the protein denatures after a 4
hour room-temperature incubation. Each set
of measurements (those of T50 and [urea]50)
was performed on all of the mutants in the
same day, and each mutant was treated identi-
cally. Therefore, it is possible to make accurate
comparisons of the relative values of the mea-
surements within the data set. However, the
absolute values of the T50 and [urea]50 values
may be less accurate. Therefore, care should be
taken in comparing the absolute value of these
measurements to those of other studies (such
as [5]).
Both the T50 and [urea]50 measurements
were performed in clarified cell lysate. The
protein was expressed using catalase-free E.
coli [25] containing the encoding gene on the
IPTG inducible pCWori [25] plasmid. We used
freshly streaked cells to inoculate 2 ml cultures
of LB supplemented with 100 µg/ml of ampi-
cillin, and grew these starter cultures overnight
with shaking at 37oC. We then used 0.5 ml
http://www.lifesci.sussex.ac.uk/CSE/test/maxchi.php
from these starter cultures to inoculate 1 L
flasks containing 200 ml of TB supplemented
with 100 µg/ml of ampicillin. The TB cultures
were grown at 30oC and 210 rpm until they
reached an optical density at 600 nm of ≈0.9,
at which point IPTG and δ-aminolevulinic acid
were added to a final concentration of 0.5 mM
each. The cultures were grown for an additional
19 hours, then the cells were harvested by pel-
letting 50 ml aliquots at 5,500 g and 4oC for 10
min, and stored at -20oC. To obtain clarified
lysate, each pellet was resuspended in 8 ml of
100 mM EPPS (pH 8.2) and lysed by sonica-
tion, while being kept on ice. The cell debris
was pelleted by centrifugation at 8,000 g and
4oC for 10 minutes, and the clarified lysate was
decanted and kept on ice.
For the T50 measurements, 125 µl of clar-
ified lysate from a single mutant was added
to all 12 wells in a row of a 96-well hard-shell
thin-wall microplate (MJ Research). The plate
was heated for 10 minutes using the gradient
method of an Eppendorf Mastercycler gradient
PCR machine, with the gradient set at either
33oC-45oC or 46oC-58oC (each mutant was ex-
posed to both of these gradients), the machine
on the BLOCK setting, and the heated lid set to
75oC with the lid WAIT option. The plate was
then cooled to 4oC, removed from the PCR ma-
chine, and centrifuged at 5,500 g and 4oC for 5
minutes to pellet any debris. A pipetting robot
was used to dispense 80 µl of the supernatent
into a 96-well microtiter plate (Rainin), and the
amount of remaining properly folded P450 was
quantified from the carbon monoxide difference
spectrum as described below. The T50 values
were determined by fitting sigmoidal curves the
percent of remaining protein as shown in Addi-
tional File 3. Our ability to accurately compare
T50 values within the data set requires that each
well in a given column of the gradient PCR ma-
chine be at the same temperature. We used
a thermocouple to measure the temperature of
the wells with the machine lid open, and con-
firmed that the wells were within a few tenths
of a degree of the same temperature. Further
evidence for the consistency of our T50 val-
ues comes from the fact that two independent
measurements of the T50 for our R1-11 par-
ent yielded values that differed by only 0.1oC.
However, the absolute values of the measured
temperatures are less accurate. Thermocouple
measurements indicated that, with the machine
lid open, the wells were ≈ 1oC cooler than the
indicated temperature. We were unable to as-
certain the temperatures with the heated lid
closed, but based on comparisons water bath
measurements, the temperatures with the lid
closed slightly exceeded the indicated tempera-
tures.
For the [urea]50 measurements, 125 µl of the
clarified lysate from a single mutant was added
to all 12 wells in a row of a 96-well microtiter
plate. A pipetting robot was then used to add
and mix 125 µl of a 2X solution of urea in 100
mM EPPS (pH 8.2) so that each subsequent
column had a higher concentration of urea, and
so that the final urea concentrations were those
shown in Additional File 4. The plates were
left on the bench at room temperature for 4
hours, and the amount of remaining properly
folded P450 was quantified from the carbon
monoxide difference spectrum as described be-
low. The [urea]50 values were determined by
fitting sigmoidal curves to the percent of re-
maining protein. Evidence for the consistency
of the [urea]50 measurements comes from the
fact that two independent measurements of the
[urea]50 for our R1-11 parent yielded values that
differed by only 0.01 M. In addition, the [urea]50
and T50 values are highly correlated (Additional
File 5), indicating that they provide consistent
measures of stability.
For both the T50 and [urea]50 measurements,
the folded P450 was quantified from the car-
bon monoxide difference spectrum [43]. The
microtiter plates containing the P450 samples
were first used to read blank spectra at 450 and
490 nm using a Tecan Safire 2 plate reader. The
plates were then incubated for 10 minutes in an
airtight oven with carbon monoxide. The plates
were removed form the oven and 10 µl of 0.1 M
sodium hydrosulfite in 1.3 M potassium phos-
phate (pH 8.0) was immediately added to each
well. After 5-10 minutes, spectra were again
read at 450 and 490 nm. The amount of P450
is proportional to the increase in the signal at
450 nm after this procedure minus the change
in the signal at 490 nm.
Mathematical Appendix
Contents
A.1 Mathematical background 12
A.2 Monomorphic limit 14
A.3 Polymorphic limit 15
A.4 Approximations for polymorphic limit 16
A.5 Approximations for monomorphic limit 18
A.6 Interpretation in terms of neutral networks 20
A.7 Detailed justification for approximating pM by po 22
A.1 Mathematical background
The first purpose of this appendix is to provide mathematical equations that describe the experi-
ments. The second is to show how four measurements from the experiments can be used to calculate
two quantities that describe the topology of the underlying protein neutral network. We will derive
two equations for both quantitites, each in terms of a different measurement. The fact that the four
equations will be seen to yield consistent results provides evidence for the accuracy of the following
calculations.
Our calculations are based on a view of neutral protein evolution as a process constrained by a
stability threshold, a view that we originally introduced to explain experimental protein mutagen-
esis results [3]. The calculations closely parallel our earlier work [22], which is in turn based on a
general theoretical treatment of evolution on neutral networks by van Nimwegen and coworkers [11].
These calculations will probably be most thoroughly understood by first reading those works. The
primary difference between the current calculations and [22] is that previously we assumed that
the per generation per protein mutation rate µ was ≪ 1, so that at each generation a protein was
either unmutated (with probability 1−µ) or experienced a single mutation (with probability µ). In
contrast, here we allow the mutation rate to be arbitrarily large, so that a protein may experience
multiple mutations in a single generation (in this sense the calculations resemble the generalization
by Wilke [18] of [11]). Specifically, let fm be the probability that a protein experiences m mutations
in a single generation. Here we derive results for arbitrary fm, and then approximations relevant to
the form of fm in the experiments. In the limiting case of small mutation rate (where f0 = 1− µ,
f1 = µ, and fm = 0 for m > 1), the calculations here reduce to those in [22]. Proteins evolving
in nature typically experience very low mutation rates, so [22] probably offers the best description
of natural protein evolution. The calculations presented here are designed to specifically treat the
evolutionary dynamics of the experiments.
A protein’s thermodynamic stability is described by its free energy of folding, ∆Gf , with more
negative values indicating more stable proteins. As described in several previous papers [3, 5, 22],
we assume that selection requires a protein to fold with some minimal stability ∆Gminf , so that a
protein adequately folds if and only if ∆Gf ≤ ∆G
f . The amount of extra stability a protein
possesses relative to the stability threshold is given by ∆Gextraf = ∆Gf − ∆G
f ; note that all
folded proteins will have ∆Gextraf ≤ 0. We further assume that as long as ∆G
extra
f ≤ 0, selection is
indifferent to the exact amount of extra stability that a protein possesses (see [22] for a discussion
of the limitations of this assumption). We conceptually divide the continuous variable of protein
stability into small discrete bins of width b. Specifically, a protein is in bin i if it has ∆Gextraf be-
tween (1− i) b and −ib, where i = 1, 2, . . .. Mutating a protein changes its stability by an amount
∆∆G (defined as the stability of the mutant protein minus the stability of the initial protein), and
so may move it to a new stability bin. In [22], we defined a matrix W with elements Wij giving
the transition probabilities that a single mutation changes a protein’s stability from bin j to bin i.
We noted that W could be computed from the distribution of ∆∆G values for all single mutations,
and argued that W remains fairly constant during neutral evolution since the distribution of ∆∆G
values remains relatively unchanged. However, we emphasize that (as discussed in detail in [22])
the constancy of the ∆∆G distribution remains an assumption, albeit one that has now been shown
to be quite accurate for lattice proteins [3, 22,44] and provide a consistent theoretical explanation
for a growing body of experimental results (the current work as well as [3]).
Since we are allowing for larger mutation rates, and we must consider the possibility that a
protein’s stability might change due to multiple mutations at a single generation. Therefore, we
make a more general definition of Wij,m as the probability that m random mutations to a protein
in stability bin j change its stability to bin i, and let Wm be the matrix with elements Wij,m. Note
that Wm only describes mutations that cause transitions from one folded protein to another, since
the stability bins i = 1, 2, . . . all correspond to folded proteins. As before [22], we assume that Wm
is roughly constant during evolution, meaning that the distribution of ∆∆G values for multiple
mutations is roughly constant during neutral evolution. Note that if m = 1, then Wm is just the
matrix W that can be computed from the distribution of single-mutant ∆∆G values [22]. We will
now use the matrices Wm to calculate the following characteristics of a population that has evolved
to equilibrium: the distribution of stabilities, the average number of mutations 〈m〉T accumulated
after T generations, and the average fraction 〈F〉 of stably folded proteins in the population. We
then introduce a few approximations (that should be quite accurate for the experimental work in
this paper) that greatly simplify these calculations. Finally, we relate the calculations to properties
of the underlying protein neutral network.
As described generally by van Nimwegen and coworkers [11], the evolutionary dynamics depend
on whether the evolving population tends to be monomorphic or highly polymorphic. When the per
sequence per generation mutation rate µ is ≪ 1, whether the population is mostly monomorphic or
highly polymorphic is determined by the product of the population size N and µ: when Nµ ≪ 1
the population is mostly monomorphic, and when Nµ ≫ 1 the population is highly polymorphic.
However, with multiple mutations per generation, Nµ is no longer an appropriate parameter to
distinguish between mono- and polymorphism, since if the population size is sufficiently small the
population can still be monomorphic even if there are multiple mutations per generation. Specifi-
cally, in one set of experiments we constrained the population to be monomorphic (by maintaining
a population size of one), but still allowed the single protein in this population to experience more
than one mutation at a generation. So we instead denote the populations as either monomorphic or
polymorphic. We indicate quantities calculated for the monomorphic population by the subscript
M (i.e. 〈F〉M ) and those calculated for the polymorphic population by the subscript P (i.e. 〈F〉P ).
A.2 Monomorphic limit
In the limit of a completely monomorphic population, all of the proteins are in a single stability
bin. Let pi (t) be the probability that the population is in stability bin i at time t, and let p (t)
be the column vector with elements pi (t). At each generation there is a probability f0 that there
is no mutation that becomes fixed in the population, a probability of
Wij,mpj that the
population experiences a mutational event (which could be a single mutation or several simultaneous
mutations) that moves it into bin i, and a probability
Wji,m that the population is in
bin i and experiences one or more mutations that move it to another bin of stably folded proteins.
Define νi,m =
Wji,m to be the fraction of m-mutants of a protein in bin i that still fold, and let
Vm be the matrix with diagonal elements given by Vii,m = νi,m and all other elements zero. The
time evolution of p is
p (t+ 1) =
fm (Wm −Vm)
p (t) (1)
where I is the identity matrix. Note that mutations that destabilize a protein beyond the stabil-
ity threshold are immediately lost to natural selection, and so leave the population in its original
stability bin. This describes the experiments for the monomorphic populations, where we retain
the parental sequence if the single mutant we generate is nonfunctional. Equation 1 corresponds to
Equation (1) of [22], and the blind ant random walk described by van Nimwegen and coworkers [11].
Equation 1 describes a Markov process with a non-negative, irreducible, and acyclic transition
matrix, and so p approaches a unique stationary distribution (equilibrium value) of pM given by
the eigenvector equation
fm (Wm −Vm)
pM. (2)
Once p has reached equilibrium, the average fraction of proteins that still stably fold at each
generation is
〈F〉M = e
pM (3)
where e = (1, . . . , 1) is the unit row vector.
To calculate 〈m〉T,M , the average number of mutations accumulated after T generations once the
population has equilibrated, we note that at each generation there is a probability of fmpj
Wij,m
that a randomly chosen protein is in bin j, experiences m mutations, and still stably folds. The
average number of mutations accumulated in a single generation is simply the average ofm weighted
over this probability. So summing over all values of m and j, we see that
〈m〉T,M = Te
mfmWmpM. (4)
This equation corresponds to Equation (6) of [22], which was derived using an embedded Markov
process formalism. Here we have foregone this formalism for the more intuitive argument presented
above, since we do not attempt to calculate higher moments of the number of mutations.
A.3 Polymorphic limit
In the limit when the population is highly polymorphic, at each generation there are sequences in
many different stability bins. In this case, we describe the distribution of stabilities by the column
vector x (t), with element xi (t) giving the fraction of proteins in stability bin i at time t. At
generation t, the fraction of mutants that continue to fold is
〈F〉t = e
x (t) . (5)
Therefore, in order to maintain a constant population size, each remaining protein must produce
an average of αt = 〈F〉t
offspring. The population therefore evolves according to
x (t+ 1) = αt
x (t) . (6)
After the population evolves for a sufficiently long period of time, x will approach an equilibrium
value of xP. At this equilibrium, the average fraction of mutants that fold at each generation is
〈F〉P = e
xP, (7)
and the equilibrium reproduction rate is α = 〈F〉P
. Therefore,
xP = α
xP. (8)
Equations 7 and 8 can be combined to show that xP and 〈F〉P can be calculated from the eigenvector
equation
(〈F〉P − f0)xP =
fmWmxP, (9)
with (〈F〉P − f0) the principal eigenvalue of the nonnegative and irreducible matrix
fmWm.
Equation 9 corresponds to Equation (14) of [22], Equation (6) of the work by van Nimwegen and
coworkers [11], and Equation (13) of the work by Wilke [18].
We now calculate 〈m〉T,P , the average number of mutations accumulated in T generations after
the population has equilibrated. At equilibrium, there is a probability of fmxj
Wij,m that a
protein is in bin j, experiences m mutations, and still stably folds. Subsequently, all of these folded
proteins produce an average of α offspring. The average number of mutations accumulated in a
single generation is simply the average of m weighted over this probability, and then multiplied by
the average reproduction rate. So summing over all values of m and j, we obtain
〈m〉T,P = αTe
mfmWmxP =
mfmWmxP. (10)
This equation is the counterpart of Equation (18) of [22], where we have again foregone the em-
bedded Markov process formalism for a more intuitive derivation.
A.4 Approximations for polymorphic limit
We can dramatically simplify the results from the previous sections with several reasonable ap-
proximations. The first approximation is that the ∆∆G values for random mutations are roughly
additive, and is supported by a number of experimental studies of the thermodynamic effects of
mutations [45–47]. We have previously shown that this approximation can be used to accurately
describe experimental protein mutagenesis data with a simple stability threshold model [3]. Under
this approximation, the distribution of net ∆∆G values for multiple mutations can be computed
from the distribution of ∆∆G values for single mutations by performing convolutions of the single-
mutation ∆∆G distribution [3], meaning that Wm for arbitrary m can be computed solely from the
distribution of ∆∆G values for single mutations. However, to simplify the equations from previous
sections, we need to express Wm for arbitrary m only in terms of W (recall that W = W1). Since
W only contains information about stability transitions from folded proteins to other folded pro-
teins, if we make the second approximation that a protein that is destabilized beyond the minimal
stability threshold by one mutation is not re-stabilized to a folded protein by a subsequent muta-
tion, then Wm = W
m. This approximation that unfolded proteins are not re-stabilized should be
quite accurate since stabilizing mutations tend to be relatively rare and small in magnitude [48–51]
(this is the underlying idea behind the Markov chain approximation that was shown to be highly
accurate for lattice proteins [44]). To summarize, if ∆∆G values are roughly additive and stabilizing
mutations are rare, we have the approximation
Wm ≈ W
m. (11)
Simplifying the equations of the previous sections also requires assigning a specific functional
form to fm, the probability that a sequence undergoesmmutations. Here we assume that mutations
are Poisson distributed among sequences, so that
e−µµm
where µ =
mfm is the average number of mutations per protein per generation. When the mu-
tations are introduced by error-prone PCR, the Poisson distribution is an excellent approximation
to the true theoretical distribution of mutations created by error-prone PCR [38,39] provided that
µ is much less than the number of PCR doublings, as is the case in all of the experiments in the
current work.
We now use the approximations of Equations 11 and 12 to simplify the results given above for
the highly polymorphic limit. We begin by using these approximations to rewrite Equation 9 as
〈F〉P − e
xP = e
WmxP. (13)
This equation makes clear that xP is the principal eigenvector of the matrix
Wm, therefore
xP must also be the principal eigenvector of W. Now in our earlier work [22], we defined the
principal eigenvector of W as x∞, called the corresponding eigenvalue 〈ν〉∞, and showed that this
eigenvalue is shown the average fraction of single mutations that are neutral in a population that
is evolving with Nµ ≫ 1 and µ ≪ 1. Therefore, with the approximation of Equation 11, xP and
x∞ are equal, and are both defined by the same eigenvector equation,
〈ν〉∞xP = WxP = Wx∞ = 〈ν〉∞x∞. (14)
Combining Equations 13 and 14 we have,
〈F〉PxP = e
(µ〈ν〉∞)
1−〈ν〉∞
xP (15)
Equation 15 can be solved to yield
〈ν〉∞ = 1 +
ln 〈F〉P
. (16)
Similarly, we can simplify Equation 10,
〈m〉T,P =
mfmWmxP
1−〈ν〉∞
eWmxP
= Te−µ〈ν〉∞
(µ〈ν〉∞)
= Tµ〈ν〉∞e
−µ〈ν〉∞
(µ〈ν〉∞)
= Tµ〈ν〉∞. (17)
Solving this equation for 〈ν〉∞ yields
〈ν〉∞ =
〈m〉T,P
. (18)
A.5 Approximations for monomorphic limit
We now simplify the equations for the monomorphic limit. This requires several further approxima-
tions. We begin by approximating that the stability probability distribution pM given by Equation
2 by the distribution po defined in [22] as satisfying
0 = (W −V)po. (19)
The basic rationale behind approximating pM with po is that Equation 2 can be viewed as a per-
turbation to Equation 19 [52]. Essentially, po is an eigenvector of the matrix W −V while pM is
the corresponding eigenvector of the matrix W − V +
(Wm −Vm). The latter matrix
can be viewed as a perturbation to the first, since the sum
(Wm −Vm) is small. This
smallness is due to the fact that Wm tends to zero with large m, causing Vm to tend towards the
identity matrix. In addition, the µm/m! terms tend to zero with large m. Therefore, the terms in
the summation are all simply either a perturbation to W − V or involve subtracting terms that
are fractions of the identity matrix. The perturbations lead to bounded changes in the eigenvec-
tors [52], while the identity matrix terms do not change the eigenvectors. Below we give a more
rigorous justification of the assumption that pM is approximately equal to po.
We need one additional approximation to make further progress. Both Equations 3 and 4 con-
tain terms of the form Wmpo, and even if we use Equation 11 to rewrite these terms as W
there are no further clear simplifications. However, any probability vector that is multiplied re-
peatedly by W and normalized will eventually converge to x∞ = xP (since this is the principal
eigenvector of W). We make the approximation that this convergence is sufficiently rapid to be
essentially complete after a single multiplication. This approximation is supported by both protein
mutagenesis studies [3,26,27] that indicate that proteins rapidly converge to an exponential decline
in the fraction folded (indicating the stability distribution has equilibrated, as discussed below, and
by lattice protein studies showing the same [3, 44]. Therefore, we make the approximation that
eWmpo = 〈ν〉oeW
m−1x∞ = 〈ν〉o〈ν〉∞
where 〈ν〉o = eWpo has the same definition as in [22],
where it was defined as the average fraction of functional single mutants of a population evolving
with µ ≪ 1 and Nµ ≪ 1.
We use these approximations to simplify Equation 3 as
〈F〉M = e
= e−µ
eWmpo
= e−µ
1 + µ〈ν〉o
(µ〈ν〉∞)
= e−µ
(µ〈ν〉∞)
= e−µ
eµ〈ν〉∞ − 1
. (20)
Solving this equation for 〈ν〉o, we find
〈ν〉o =
〈ν〉∞ (〈F〉Me
µ − 1)
eµ〈ν〉∞ − 1
. (21)
We now use the approximations to simplify Equation 4 as
〈m〉T,M = Te
mfmWmpM
= Te−µ
eWmpo
= Te−µ〈ν〉o
= µTe−µ〈ν〉o
(µ〈ν〉∞)
= µT 〈ν〉oe
〈ν〉∞−1
. (22)
Solving this equation for 〈ν〉o yields
〈ν〉o =
〈m〉T,Me
1−〈ν〉∞
. (23)
To recap, we now have equations to calculate 〈ν〉∞ and 〈ν〉o from experimentally measurable
quantities. Equations 16 and 18 allow us to calculate 〈ν〉∞ from 〈F〉P and 〈m〉T,P , respectively.
Given this calculated value of 〈ν〉∞, Equations 21 and 23 then allow us to calculate 〈ν〉o from 〈F〉M
and 〈m〉T,M , respectively. The fact that we have two equations each for 〈ν〉∞ and 〈ν〉o allows us
to assess the self-consistency of the approach.
A.6 Interpretation in terms of neutral networks
Throughout the preceding calculations, we have referred to 〈ν〉∞ and 〈ν〉o as we defined them
in [22]: namely, as the average neutrality of protein populations evolving with µ ≪ 1 and Nµ
either ≫ 1 or ≪ 1, respectively. However, van Nimwegen and coworkers [11] have shown that
they can also be interpreted in terms of the underlying neutral network. In the experiments we
make mutations at the nucleotide (rather than amino acid) level, so each point in our sequence
space corresponds to a different gene. Every gene that yields an amount of protein sufficient to
hydroxylate the twelfth carbon of 12-p-nitrophenoxydodecanoic acid with at least 75% of the total
activity conferred by the original R1-11 parent gene represents a node on this neutral network.
We note that in the experiments (and also usually in natural evolution), the edges on the neutral
network are not all completely equivalent or fully undirected, since some mutations are more likely
to occur than others (for example, error-prone PCR with Taq polymerase is more likely to cause an
A→G mutation than an A→C mutation). In the analysis that follows, we ignore this complication
and assume all neutral network edges are equivalent.
In an extremely insightful analysis, van Nimwegen and coworkers [11] have shown that impor-
tant characteristics of a neutral network can be inferred from evolutionary quantities. Specifically,
they have shown that if a population is evolving with µ ≪ 1 and Nµ ≫ 1, then the average
neutrality (which we have denoted by 〈ν〉∞) is equal to the principal eigenvalue of the adjacency
matrix of the neutral network, normalized by the network coordination number (number of possi-
ble connections per node). In addition, they pointed out that a population evolving with µ ≪ 1
and Nµ ≪ 1 moves like a blind ant random walk, meaning that the average neutrality (which we
have denoted by 〈ν〉o) is equal to the average connectivity of a neutral network node divided by
the network coordination number. In our P450 experiments, we have measured the values needed
to estimate 〈ν〉∞ and 〈ν〉o using Equations 16, 18, 21, and 23. Using the final values listed in
Table 2, 〈F〉P = 0.50 and 〈F〉M = 0.39. Taking the final nucleotide mutation values from Table
2, 〈m〉T,P /T = 0.69 and 〈m〉T,M/T = 0.31. The average mutation rate, computed from the un-
selected population, is µ = 1.40. So using Equation 16, 〈ν〉∞ = 0.53, while using Equation 18,
〈ν〉∞ = 0.49. The consistency of these two values supports the idea that the calculations above
accurately describe the evolutionary process. Taking the average value of these two measurement
as 〈ν〉∞ = 0.51, we can then use Equations 21 and 23 to calculate 〈ν〉o. We calculate values of
0.28 and 0.43, respectively. These estimates differ by more than those for 〈ν〉∞, perhaps because
additional approximations have gone into the derivation of the relevant equations (in addition, we
have made no attempt to carry out the rather complex propagation of the sampling errors of Table
2). However, the values are still in a similar range. Taking the average of these two values, we
estimate that 〈ν〉o = 0.35. So overall, we predict that each functional P450 gene should have an
average fraction of 0.35 of its sequence nearest neighbors also encoding a functional gene, for an
average of about 1,500 neighbor genes. We predict that the principal eigenvalue of the neutral
network adjacency matrix is 0.51 ×3L. The fact that principal eigenvalue exceeds the average
connectivity indicates that the neutral network is not a regular graph, but instead has some nodes
more highly connected than others.
The value for 〈ν〉∞ calculated above can also be related to measurements from protein mu-
tagenesis experiments. A number of studies [3, 26, 27] have observed that the probability that a
protein remains functional after m mutations falls off exponentially with the number of mutations.
In fact, the decline is not always exponential for the first few mutations if the starting protein has
especially high or low stability [3] or activity [53], but will still converge to this exponential form
after a few mutations [3,44,54]. The stability threshold model can be used to relate this decline to
〈ν〉∞, as is done indirectly in the Markov chain approximation of [44]. Here we make that connec-
tion explicit. The initial protein has a stability that falls into some stability bin i. Therefore, its
stability can be described by the column vector y0, which has element i equal to one and all other
elements equal to zero. Now imagine constructing all single mutants of this protein. The fraction
of these single mutants that still fold is just eWy0, and the distribution of stabilities among the
single mutants is y1 = Wy0 (note that the elements of y1 no longer sum to one). Similarly, after
m mutations, the fraction of mutants that still fold is eWmy0, and the distribution of stabilities
among the m-mutants is ym = Wmy0. With the approximation of Equation 11, ym = W
This makes it clear that ym will converge to a vector proportional to x∞, the principal eigenvector
of W. Once this convergence is complete, each new mutation simply reduces the fraction of mu-
tants that fold by a factor of 〈ν〉∞, the principal eigenvalue of W (and the spectral radius of the
neutral network normalized by the coordination number). Therefore, what we have called 〈ν〉∞ in
the present work and [22] is equal to what is called x in [27], q in [26], and 〈ν〉 in [3]. The major
difficulty that is faced in extracting 〈ν〉∞ by the method of those three studies [3, 26,27] is that it
is not possible to directly assay mutants with m mutations, but instead only possible to assay a set
of mutants with a distribution of m. All three studies use different (and valid) methods to account
for this distribution, but this accounting is still difficult because most of the functional mutants
come from the low m end of the distribution. This makes it hard to get accurate values for the
fraction functional after large numbers of mutations, since most of the functional mutants in the
set come from sequences with few mutations. For this reason, we believe the current method of
measuring 〈ν〉∞ is more accurate. A second caution about comparing values of 〈ν〉∞ from different
studies is that its value depends on the nucleotide error-spectrum of the experiment, since differ-
ent mutagenesis methods create different distributions of nucleotide and amino acid mutation types.
We also briefly mention how we arrived at an estimate of 〈ν〉∞ for 3-methyladenine DNA
glycosylase from the data of [27]. This paper reports that a fraction x = 0.34 of amino acid
mutations inactivate the protein. We would like to determine the fraction 〈ν〉∞ of nucleotide mu-
tations that do not inactivate the protein. Roughly 75% of random mutations to a gene will be
synonymous. Therefore, m amino acid mutations should cause about 4m/3 nucleotide mutations.
The study of [27] measures that after m mutations, a fraction (1− x)
of the mutants are func-
tional. That means that 〈ν〉∞
fraction should be functional. Equating these expressions yields
〈ν〉∞ = exp
log (1− x)
. So using x = 0.34, we arrive at 〈ν〉∞ = 0.73.
A.7 Detailed justification for approximating pM by po
Here we provide a detailed justification for the approximation that pM is about equal to po. In the
monomorphic limit, the time evolution of p is given by Equation 1, and the stationary distribution
pM is given by Equation 2. We assume the approximations of Equations 11 and 12 and show that
we can approximate pM by po, where po is given by Equation 19. To justify this approximation,
we insert po into the right hand side of Equation 1 and ask to what extent po is left unaltered by
the dynamics. If po is found to be stationary to good approximation then, by uniqueness of the
stationary distribution of an ergodic process, po would be a good approximation to pM.
We therefore suppose that at some time t the distribution is given by po and compute, using
Equation 1, the change in po after one generation. The new distribution at time t+ 1 is given by
p (t+ 1) =
fm (W
m −Vm)
po. (24)
Using (V −W)po = 0, and taking components of the above equation, we obtain
pi (t+ 1) = p0i +
fm [(W
m −Vm)po]i . (25)
Thus po would be an approximately stationary distribution of the dynamics if |
fm [(W
m −Vm)po]i| ≪
p0i. We now proceed to show that this will be the case in most situations of interest by deriving
upper and lower bounds on the second term of the right hand side of Equation 25.
Consider first the term (Wmpo)i, which can be written as
(Wmpo)i =
k1,...,km
Wik1Wk1k2 · · ·Wkm−1kmp0km
k1,...,km−1
Wik1Wk1k2 · · ·Wkm−2km−1νkm−1p0km−1 , (26)
where we have used Wpo = Vpo in the second equality. We now note that νk ≤ νmax for all
k, where νmax is the maximum neutrality, maximized over all bins. This leads to the successive
inequalities
(Wmpo)i ≤ νmax
k1,...,km−1
Wik1Wk1k2 · · ·Wkm−2km−1p0km−1
= νmax
k1,...,km−2
Wik1Wk1k2 · · ·Wkm−3km−2νkm−2p0km−2
k1,...,km−2
Wik1Wk1k2 · · ·Wkm−3km−2p0km−2
≤ νm−1
Wik1p0k1 , (27)
yielding the upper bound
(Wmpo)i ≤ ν
νip0i. (28)
In an identical manner, we obtain the lower bound
(Wmpo)i ≥ ν
νip0i, (29)
where νmin is the smallest neutrality, minimized over all bins. Note that both inequalities above
become exact equalities when all bins have the same neutrality ν, which could be interpreted as
either νmin or νmax.
Having obtained inequality constraints on (Wmpo)i, we now consider the term (Vmpo)i, which
can be written as
(Vmpo)i = p0iνi,m
= p0i
(Wm)ji
= p0i
j,k1,...,km−1
Wjk1Wk1k2 · · ·Wkm−1i
= p0i
k1,...,km−1
νk1Wk1k2 · · ·Wkm−1i
≤ p0iνmax
k1,...,km−1
Wk1k2 · · ·Wkm−1i
≤ p0iν
Wkm−1i, (30)
which yields an identical upper bound to that on (Wmpo)i, namely,
(Vmpo)i ≤ ν
νip0i, (31)
and similarly
(Vmpo)i ≥ ν
νip0i. (32)
It should again be noted that both the above inequalities become exact equalities when all bins
have a common neutrality ν.
We are now in a position to estimate bounds on the magnitude of the second term of Equation
25. Using the four inequalities of Equations 28, 29, 31, and 32 above, we have
− νm−1
νip0i ≤ [(W
m −Vm)po]i ≤
− νm−1
νip0i, (33)
or equivalently,
|[(Wm −Vm)po]i| ≤
− νm−1
νip0i, (34)
where the inequality above becomes an exact equality when all bins have the same neutrality. How-
ever, in this limit, the right hand side of the above equation vanishes, and therefore the second
term of Equation 25 is identically zero in this case, giving the result that pM is exactly equal to
po when all bins have the same neutrality, even if µ is arbitrarily large.
We now carry out the sum over m to obtain an upper bound on the second term of Equation 25
in the more general and realistic case of unequal neutrality bins. Using Equation 34 and the specific
Poisson form of fm, we obtain an upper bound on the fractional change in p0i in one generation:
pi(t+ 1)− p0i
≤ νie
− νm−1
= νie
eµνmax − 1
eµνmin − 1
. (35)
The above bound vanishes for small µ, is an increasing function of νmax − νmin, and is typically
much smaller than 1. An extreme estimate of the size of the fractional change can be made when
νmax = 1 and νmin = 0. In this case, using µ = 1.4 (the value in our experiments), the above
inequality simplifies to
pi(t+ 1)− p0i
1− e−µ − µe−µ
≃ 0.41νi. (36)
Noting that νi < 1, the fractional change in p0i is therefore reasonably controlled even in the most
extreme case. For realistic situations, the fractional change in p0i is expected to be much lower,
thus justifying the use of po as the stationary distribution of the dynamics of Equation 1.
Author Contributions
JDB and FHA designed the project and wrote the paper. JDB and ZL performed the bulk of the
experiments; OSV assisted with the experiments. JDB and DC analyzed the data. JDB and AR
performed the theoretical work.
Acknowledgments
We thank Claus O Wilke for helpful advice and comments. JDB is supported by a HHMI predoctoral
fellowship. ZL and DC were supported by Summer Undergraduate Research Fellowships from the California
Institute of Technology. AR is supported by NSF grants CCF 0523643 and FIBR 0527023.
References
1. Zuckerkandl E, Pauling L: Evolutionary divergence and convergence in proteins. In Evolving genes and
proteins, New York, NY: Academic Press 1965:97–166.
2. Lesk AM, Chothia C: How different amino acid sequences determine similar protein structures: The
structure and evolutionary dynamics of the globins. J. Mol. Biol. 1980, 136:225–270.
3. Bloom JD, Silberg JJ, Wilke CO, Drummond DA, Adami C, Arnold FH: Thermodynamic prediction of
protein neutrality. Proc. Natl. Acad. Sci. USA 2005, 102:606–611.
4. Besenmatter W, Kast P, Hilvert D:Relative tolerance of mesostable and thermostable protein homologs
to extensive mutation. Proteins 2007, 66:500–506.
5. Bloom JD, Labthavikul ST, Otey CR, Arnold FH: Protein stability promotes evolvability. Proc. Natl. Acad.
Sci. USA 2006, 103:5869–5874.
6. Lenski RE, Barrick JE, Ofria C: Balancing robustness and evolvability. PLoS Biology 2006, 4:e428.
7. Sniegowski PD, Murphy HA: Evolvability. Current Biology 2006, 16:R831–R834.
8. Wagner A: Robustness and Evolvability in Living Systems. Princeton, New Jersey: Princeton University Press
2005.
9. Montville R, Froissart R, Remold SK, Tenaillon O, Turner PE: Evolution of mutational robustness in an
RNA virus. PLoS Biol. 2005, 3:e381.
10. Wilke CO, Wang JL, Ofria C, Lenski RE, Adami C: Evolution of digital organisms at high mutation
rates leads to survival of the flattest. Nature 2001, 412:331–333.
11. van Nimwegen E, Crutchfield JP, Huynen M: Neutral evolution of mutational robustness. Proc. Natl. Acad.
Sci. USA 1999, 96:9716–9720.
12. Kimura M: The Neutral Theory of Molecular Evolution. Cambridge, U.K.: Cambridge Univ. Press 1983.
13. Takahata N: On the overdispersed molecular clock. Genetics 1987, 116:169–179.
14. Smith JM: Natural selection and the concept of a protein space. Nature 1970, 225:563–564.
15. Lipman DJ, Wilbur WJ: Modeling neutrality and selective evolution of protein folding. Proc. R. Soc.
London Ser. B 1991, 245:7–11.
16. Huynen MA, Stadler PF, Fontana W: Smoothness within ruggedness: the role of neutrality in adapta-
tion. Proc. Natl. Acad. Sci. USA 1996, 93:397–401.
17. Bornberg-Bauer E, Chan HS: Modeling evolutionary landscapes: mutational stability, topology, and
superfunnels in sequence space. Proc. Natl. Acad. Sci. USA 1999, 96:10689–10694.
18. Wilke CO: Adaptive evolution on neutral networks. Bull. Math. Biol. 2001, 63:715–730.
19. Taverna DM, Goldstein RA: Why are proteins so robust to site mutations? J. Mol. Biol. 2002, 315:479–
20. Xia Y, Levitt M: Funnel–like organization in sequence space determines the distributions of protein
stability and folding rate preferred by evolution. Proteins 2004, 55:107–114.
21. Brin S, Page L: The anatomy of a large-scale hypertextual Web search engine. Computer Networks and
ISDN Systems 1998, 30(1–7):107–117.
22. Bloom JD, Raval A, Wilke CO: Thermodynamics of neutral protein evolution. Genetics 2007, 175:255–
23. Cirino PC, Arnold FH: A self-sufficient peroxide-driven hydroxylation biocatalyst. Angew. Chem. Int.
Ed. 2003, 42:3299–3301.
24. Taverna DM, Goldstein RA: Why are proteins marginally stable? Proteins 2002, 46:105–109.
25. Barnes HJ, Arlotto MP, Waterman MR: Expression and enzymatic activity of recombinant cytochrome
P450 17 α-hydroxylase in Escherichia coli . Proc. Natl. Acad. Sci. USA 1991, 88:5597–5601.
26. Shafikhani S, Siegel RA, Ferrari E, Schellenberger V: Generation of large libraries of random mutants in
Bacillus subtilis by PCR-based plasmid multimerization. BioTechniques 1997, 23:304–310.
27. Guo HH, Choe J, Loeb LA: Protein tolerance to random amino acid change. Proc. Natl. Acad. Sci. USA
2004, 101:9205–9210.
28. DePristo MA, Weinreich DM, Hartl DL: Missense meanderings in sequence space: a biophysical view
of protein evolution. Nat. Rev. Genetics 2005, 6:678–687.
29. Plotkin JB, Dushoff J, Desai MM, Fraser HB: Codon usage and selection in proteins. J. Mol. Evol. 2006,
63:635–653.
30. Berg OG: Selection intensity for codon bias and the effective population size of Escherichia coli .
Genetics 1996, 142:1379–1382.
31. Lynch M, Conery JS: The origins of genome complexity. Science 2003, 302:1401–1404.
32. Hartl DL, Moriyama EN, Sawyer SA: Selection intensity for codon bias. Genetics 1994, 138:227–234.
33. Moya A, Holmes EC, Gonzalez-Candelas F: The population genetics and evolutionary epidemiology of
RNA viruses. Nature Reviews Microbiology 2004, 2:279–288.
34. Vignuzzi M, Stone JK, Arnold JJ, Cameron CE, Andino R:Quasispecies diversity determines pathogenesis
through cooperative interactions in a viral population. Nature 2006, 439:344–348.
35. Wang X, Minasov G, Shoichet BK: Evolution of an antibiotic resistance enzyme constrained by stability
and activity trade–offs. J. Mol. Biol. 2002, 320:85–95.
36. Meyerhans A, Vartanian JP, Wain-Hobson S: DNA recombination during PCR. Nucleic Acids Res. 1990,
18:1687–1691.
37. Kanagawa T: Bias and artifacts in multitemplate polymerase chain reactions (PCR). Journal of
Bioscience and Bioengineering 2003, 96:317–323.
38. Sun F: The polymerase chain reaction and branching processes. J. Comput. Biol. 1995, 2:63–86.
39. Drummond DA, Iverson BL, Georgiou F, Arnold FH: Why high–error–rate random mutagenesis libraries
are highly enriched in functional and improved proteins. J. Mol. Biol. 2005, 350:806–816.
40. Posada D: Evaluation of methods for detecting recombination from DNA sequences: empirical data.
Mol. Biol. Evol. 2002, 19:708–717.
41. Smith JM: Analyzing the mosaic structure of genes. J. Mol. Evol. 1992, 34:126–129.
42. Piganeau G, Gardner M, Eyre-Walker A: A broad survey of recombination in animal mitochondria. Mol.
Biol. Evol. 2004, 21:2319–2325.
43. Otey CR: High-throughput carbon monoxide binding assay for cytochromes P450. In Directed enzyme
evolution: screening and selection methods, Volume 230 of Methods in Molecular Biology. Edited by Arnold FH,
Georgiou G, Humana press 2003:137–139.
44. Wilke CO, Bloom JD, Drummond DA, Raval A: Predicting the tolerance of proteins to random amino
acid substitution. Biophysical J. 2005, 89:3714–3720.
45. Wells JA: Additivity of mutational effects in proteins. Biochemistry 1990, 29:8509–8517.
46. Zhang XJ, Baase WA, Shoichet BK, Wilson KP, Matthews BW: Enhancement of protein stability by the
combination of point mutations in T4 lysozyme is additive. Protein Eng. 1995, 8:1017–1022.
47. Serrano L, Day AG, Fersht AR: Step-wise mutation of barnase to binase: a procedure for engineering
increased stability of proteins and an experimental analysis of the evolution of protein stability.
J. Mol. Biol. 1993, 233:305–312.
48. Godoy-Ruiz R, Perez-Jimenez R, Ibarra-Molero B, Sanchez-Ruiz JM: Relation between protein stability,
evolution and structure as probed by carboxylic acid mutations. J. Mol. Biol. 2004, 336:313–318.
49. Pakula AA, Sauer RT: Genetic analysis of protein stability and function. Annu. Rev. Genet. 1989, 23:289–
50. Matthews BW: Structural and genetic analysis of protein stability. Annu. Rev. Biochem. 1993, 62:139–
51. Bava KA, Gromiha MM, Uedaira H, Kitajimi K, Sarai A: ProTherm, version 4.0: thermodynamic
database for proteins and mutants. Nucleic Acids Res. 2004, 32:D120–D121.
52. Franklin JN: Matrix Theory. Englewood Cliffs, N.J.: Prentice-Hall, Inc. 1968.
53. Bershtein S, Segal M, Bekerman R, Tokuriki N, Tawfik DS: Robustness-epistasis link shapes the fitness
landscape of a randomly drifting protein. Nature 2006, 444:929–932.
54. Bloom JD, Arnold FH, Wilke CO: Breaking proteins with mutations: threads and thresholds in evo-
lution. Molecular Systems Biology 2007, 3:76.
Figures
monomorphic polymorphic
BPSfrag replacements
0.2probability
threshold
more stable← ∆Gf → less stable
PSfrag replacements
threshold
more stable← ∆Gf → less stable
monomorphic
polymorphic
Figure 1 - Theoretical views of the evolution of protein mutational robustness and stability.
(A) Theory predicts that a mostly monomorphic population is equally likely to occupy any node of its neutral
network, while a highly polymorphic population will prefer more connected nodes [11]. Node sizes are drawn
proportional to the occupation probabilities. (B) Proteins evolving in a highly polymorphic population
are predicted to be more stable than their counterparts in a mostly monomorphic population [22]. The
histograms illustrate the distributions of stabilities for the two cases. The increased stability is a biophysical
manifestation of excess mutational robustness, since more stable proteins are more mutationally robust [3–5].
polymorphic monomorphic unselected
Figure 2 - Outline of the neutral evolution experimental procedure.
For the polymorphic population, error-prone PCR was used to generate mutant P450 genes. These genes
were ligated into a plasmid and transformed into E. coli. Individual mutants (435) were picked, expressed
in E. coli, and assayed for enzymatic activity. All mutants that met the selection criterion contributed an
equal amount of plasmid DNA as template for the next generation of error-prone PCR. The monomorphic
populations were treated similarly, except only a single mutant was assayed at each generation. If this mutant
met the selection criterion then it became the template for the next generation of error-prone PCR; otherwise
at the next generation another colony was picked from the same plate. In the unselected populations a
single mutant was picked and used as the template for the next generation of error-prone PCR.
PSfrag replacements
generation
〈mnt〉
unselected
polymorphic
monomorphic
PSfrag replacements
0 10 20
generation
〈maa〉 〈
unselected
polymorphic
monomorphic
Figure 3 - Accumulation of nucleotide (〈mnt〉) and nonsynonymous (〈maa〉) mutations in the
experimentally evolved P450 populations.
For the unselected and monomorphic populations, numbers are the average over all replicates at the indicated
generation; for the polymorphic population they are from a random sample, with sampling standard error
shown.
PSfrag replacements
0 10 20
generation
polymorphic
monomorphic
Figure 4 - The polymorphic population neutrally evolved a higher average mutational robustness
than the monomorphic populations.
The fraction functional was determined by assaying 435 mutants (average of 1.5 nucleotide mutations per
gene). Error bars show binomial standard error. For the monomorphic population, numbers are the average
over all replicates.
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
38 40 42 44
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
number
0.4 0.6 0.8 1.0
PSfrag replacements
polymorphic
monomorphic
T50 (
[urea]50 (M)
percent parental expression
number
20 60 100 140
Figure 5 - The more mutationally robust proteins are more stable.
The P450s from the polymorphic population neutrally evolved higher stability and expression levels than
their counterparts from the monomorphic populations. The histograms show the distributions for the final
protein from all monomorphic replicates and for the same number of randomly chosen proteins from the final
polymorphic population. The plots show (left to right) the temperature at which half the protein irreversibly
denatured after 10 minutes (T50), the urea concentration at which half the protein denatured after 4 hours
([urea]50), and the expression level relative to that of the original parental P450. The means are significantly
different, with unequal variance t-test P -values of 0.02, 0.005, and 0.04, respectively.
Tables
Table 1 - Error-prone PCR nucleotide mutation spectrum.
Spectrum of nucleotide mutations introduced by the error-prone PCR procedure used in the neutral evolution
experiments. The spectrum was determined by sequencing the four final (generation 12) sequences from the
unselected population, since in these sequences the mutations accumulate without constraint. As has been
previously noted for error-prone PCR with Taq polymerase [3,5,26], the nucleotide error spectrum is biased
towards certain types of mutations.
Total nucleotide mutations 67
% synonymous mutations 25
Mutation types (%)
A →T, T →A 19.4
A →C, T →G 1.5
A →G, T →C 64.2
G →A, C →T 4.5
G →C, C →G 0.0
G →T, C →A 1.5
frameshift 9.0
Table 2 - Neutral evolution robustness and mutation data.
Each row is for a different generation, T . Entries of NA indicate that no measurement was made. The 〈mnt〉
and 〈maa〉 are the average number of nucleotide mutations and nonsynonymous mutations, respectively.
Numbers in parentheses are total counts over the total samples. Subscripts indicate the population type:
U for unselected, P for polymorphic, and M for monomorphic. For the unselected and monomorphic
populations, numbers represent averages of all replicates. For the polymorphic population, numbers are for
a random sample of functional mutants. 〈F〉P and 〈F〉M are the fraction of functional mutants out of 435
assayed.
T 〈mnt〉U 〈maa〉U 〈mnt〉P 〈maa〉P 〈mnt〉M 〈maa〉M 〈F〉P 〈F〉M
0 0 0 0 0 0 0 0.48 (210 / 435) 0.48 (210 / 435)
1 NA NA NA NA 0.1 (3 / 22) 0.3 (6 / 22) 0.48 (208 / 435) NA
2 NA NA NA NA 0.4 (9 / 22) 0.8 (17 / 22) 0.49 (215 / 435) NA
3 5.0 (20 / 4) 3.5 (14 / 4) 2.7 (27 / 10) 1.4 (14 / 10) 1.0 (23 / 22) 0.4 (9 / 22) 0.49 (215 / 435) NA
4 NA NA NA NA 1.5 (32 / 22) 0.7 (15 / 22) 0.48 (208 / 435) NA
5 NA NA NA NA 2.2 (48 / 22) 1.1 (25 / 22) 0.45 (197 / 435) 0.43 (185 / 435)
6 9.8 (39 / 4) 7.5 (30 / 4) 5.5 (55 / 10) 2.1 (21 / 10) 2.6 (58 / 22) 1.4 (31 / 22) 0.46 (198 / 435) NA
7 NA NA NA NA 3.1 (69 / 22) 1.8 (39 / 22) 0.52 (227 / 435) NA
8 NA NA NA NA 3.4 (74 / 22) 1.8 (40 / 22) 0.46 (200 / 435) NA
9 13.0 (52 / 4) 10.3 (41 / 4) 6.7 (61 / 9) 3.1 (28 / 9) 3.7 (82 / 22) 2.1 (46 / 22) 0.47 (203 / 435) NA
10 NA NA NA NA 4.2 (92 / 22) 2.4 (52 / 22) 0.46 (199 / 435) 0.40 (175 / 435)
11 NA NA NA NA 4.6 (102 / 22) 2.5 (56 / 22) 0.48 (207 / 435) NA
12 16.8 (67 / 4) 12.5 (50 / 4) 7.8 (70 / 9) 3.3 (30 / 9) 4.9 (107 / 22) 2.6 (58 / 22) 0.52 (228 / 435) NA
13 NA NA NA NA 5.0 (110 / 22) 2.7 (60 / 22) 0.52 (227 / 435) NA
14 NA NA NA NA 5.3 (116 / 22) 2.9 (64 / 22) 0.50 (216 / 435) NA
15 NA NA 10.3 (227 / 22) 3.8 (83 / 22) 5.6 (123 / 22) 3.0 (67 / 22) 0.50 (219 / 435) 0.39 (171 / 435)
16 NA NA NA NA 5.8 (127 / 22) 3.0 (67 / 22) NA NA
17 NA NA NA NA 6.0 (133 / 22) 3.1 (69 / 22) NA NA
18 NA NA NA NA 6.3 (137 / 22) 3.2 (71 / 22) NA NA
19 NA NA NA NA 6.3 (138 / 22) 3.3 (72 / 22) NA NA
20 NA NA NA NA 6.6 (145 / 22) 3.4 (75 / 22) NA 0.37 (160 / 435)
21 NA NA NA NA 6.9 (152 / 22) 3.6 (79 / 22) NA NA
22 NA NA NA NA 7.1 (156 / 22) 3.7 (81 / 22) NA NA
23 NA NA NA NA 7.2 (158 / 22) 3.7 (81 / 22) NA NA
24 NA NA NA NA 7.3 (161 / 22) 3.8 (83 / 22) NA NA
25 NA NA NA NA 7.7 (169 / 22) 4.0 (87 / 22) NA 0.39 (169 / 435)
Additional material
Additional file 1 - Sequence of the parent P450 used start neutral evolution.
FASTA file with sequence of the R1-11 P450 BME used as the neutral evolution parent. This sequence was
isolated after the equilibration evolution.
Additional file 2 - Information about sequences from neutral evolution experiments.
The entries give the name of the mutant, the number of nonsynonymous and nucleotide mutations relative
to the R1-11 parent, the [urea]50 value if measured, the T50 value if measured, the percent of the parental
expression level if measured, and then a list of all of the mutations. Amino acid mutations are numbered in
the standard P450 numbering scheme. The names of the mutants indicate their origin. Names beginning
with “P-G3” are randomly chosen functional mutants from generation 3 of the polymorphic population,
etc. Names of the form “P1,” “P2,”, etc. are the 22 functional mutants that were randomly chosen from
the final (generation 15) polymorphic population. Numbers P5 and P12 are missing because two of the
original 24 randomly selected polymorphic population members were randomly chosen to be discarded after
it was discovered that two of the 24 monomorphic replicates were contaminated. Names beginning with
“U1” indicate that sequences are from the first unselected replicate, etc. Names beginning “M1” indicate
sequences are from the first monomorphic replicate, etc. Replicates “M9” and “M10” were discarded due to
contamination during the experiment. For each replicate, we sequenced each new functional mutant. The
last functional mutant after 25 generations represents the final sequence for that replicate, and is given an
abbreviated name without the generation suffix.
Additional file 3 - Thermostability measurements.
Raw data from the T50 thermostability measurements.
Additional file 4 - Urea stability measurements.
Raw data from the [urea]50 thermostability measurements.
Additional file 5 - Correlation of thermal and urea stabilities.
The T50 and [urea]50 values are highly correlated.
Additional file 6 - Sequence of initial P450 used to start equilibration evolution.
FASTA file with sequence of the 21B3 P450 BM3 heme domain described in [23]. This P450 was used as
the initial parent to start the equilibration evolution.
Additional file 7 - Mutations accumulated during equilibration evolution.
The file lists the mutations in the 46 P450 variants selected at the end of the equilibration evolution. Each line
gives the name of the variant, with the prefix indicating whether it came from the R1 or R2 population. The
next entries give the number of nucleotide and nonsynonymous mutations. All of the individual mutations
relative to 21B3 are then listed. Amino acid mutations are numbered in the standard P450 numbering
scheme, with the threonine after the N-terminal methionine given the number one.
	A.1 Mathematical background
	A.2 Monomorphic limit
	A.3 Polymorphic limit
	A.4 Approximations for polymorphic limit
	A.5 Approximations for monomorphic limit
	A.6 Interpretation in terms of neutral networks
	A.7 Detailed justification for approximating pM by po
