Phenotypic switching can speed up biological evolution of microbes

0.63 MB·

by  Andrew C. Tadrowski

Phenotypic switching can speed up biological evolution of microbes
Andrew C. Tadrowski1, Martin R. Evans1, and Bartl omiej Waclaw1,2
1SUPA, School of Physics and Astronomy, The University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom
2Centre for Synthetic and Systems Biology, The University of Edinburgh
Stochastic phenotype switching has been suggested to play a beneficial role in microbial pop-
ulations by leading to the division of labour among cells, or ensuring that at least some of the
population survives an unexpected change in environmental conditions. Here we use a computa-
tional model to investigate an alternative possible function of stochastic phenotype switching - as
a way to adapt more quickly even in a static environment. We show that when a genetic mutation
causes a population to become less fit, switching to an alternative phenotype with higher fitness
6
(growth rate) may give the population enough time to develop compensatory mutations that in-
1
crease thefitnessagain. Thepossibility of switching phenotypescan reducethetimetoadaptation
0
by orders of magnitude if the “fitness valley” caused by the deleterious mutation is deep enough.
2
Ourworkhasimportantimplicationsfortheemergenceofantibiotic-resistantbacteria. Inlinewith

b recent experimental findings we hypothesise that switching to a slower growing but less sensitive
e phenotypehelpsbacteriatodevelopresistancebyexploringalargersetofbeneficialmutationswhile
F
avoiding deleterious ones.

1
Keywords: evolution|populationgenetics |stochasticphenotype switching

]
E
INTRODUCTION in a fluctuating and unpredictable environment it pays
P
. to have a fraction of the population in a different state,
o
which is perhaps maladapted to the present environ-
i Biological evolution relies on two mechanisms which
b ment but better suited to possible future environments.
are instrumental in natural selection: preferential sur-
- Since only a small fraction of the population expresses
q vival of better adapted individuals (selection) and vari-
[ ations among individuals (phenotypic variation). One the maladapted phenotype at any one time, this strat-
of the sources of phenotypic variability is genetic alter- egy conserves resourceswhile allowing the population to
2
stay prepared for an unexpected change. Examples in-
v ation due to mutations and recombination. However,
clude bacterial persisters [4, 25], flu(Ag43)/fim switch
0 evengeneticallyidenticalorganismswilloftenbehavedif-
0 ferently because the same genotype may lead to many [3, 4, 19, 44] and competence to non-competence switch-
6 ing in the bacterium B. subtilis [29]. Animal cells also
different phenotypes - observable traits of an organism -
4 exhibit this behaviour. Epithelial-to-mesenhymaltransi-
as a result of environmental factors and the organism’s
0
tion–alandmarkincancerprogression–isthoughttobe
. history. Although ubiquitous and easily observedin ani-
1 a phenotypic, epigenetic change [33]. There is also some
malsandplants,phenotypicdiversitycanalreadybeam-
0
evidence that phenotypic switching may be involved in
6 ply demonstrated in microorganisms. Examples range
resistance to chemotherapy [21], although this remains
1 from different cell sizes depending on growth medium
: [49], through bistability in utilization of different food controversial[50].
v
i sources[35],todiversificationbetweenmotile/non-motile Here weinvestigateanalternativepossible rolefor the
X
cells [20]. Microorganisms are often able to switch be- evolution of phenotype switching in microbes: the facili-
r tweenthese phenotypesinresponseto achangeinexter- tation of genetic evolution in a static environment. Pre-
a
nal conditions such as the arrival of a new food source vioustheoreticalwork[28]hasdemonstratedthatswitch-
or depletion of the currently used one. A typical exam- ingtoapersistentphenotypecanprovidealargerpoolof
ple is diauxic shift - a switch to another food source, for bacteria to evolve genetic resistance to antibiotics. That
example from glucose to cellobiose in L. lactis when glu- work only addressed the situation in which resistance is
cose becomes depleted [41], which involves altering gene caused by a single, beneficial mutation. However, a sig-
expression levels without changing the genetic code. nificant change of fitness, defined here as a measure of a
Some microorganisms switch seemingly randomly be- microorganism’sgrowthrate,oftenrequiresmultiplemu-
tween two or more phenotypes even in the absence of tations. Although many fitness landscapes have at least
external stimuli. This causes the population to become one accessible pathway [14, 46], along which the fitness
phenotypicallyheterogeneous. Severalexplanationshave increases monotonically, in some cases a fitness valley -
been proposed as to why stochastic phenotype switching a genotype with a lower fitness than all its neighbours
has evolved [1]. One of them is the division of labour [2] - is unavoidable. For example, developing resistance to
in which different microbial cells perform different func- the antibiotic streptomycin involves a fitness cost which
tions and this maximizes the benefit to the population. must be counteracted by compensatory mutations, i.e.
Another theory,calledbet hedging [6,26], proposesthat subsequent mutations that increase the microorganism’s
2
A flat. We showthat this switching enablesthe population
μ μ
1B 2B 3B to avoid the costly valley, even if the alternative pheno-
Phenotype typeislessfitthantheinitialstate. Wealsodemonstrate
the existence of an optimal range for the switching rate.
α α α
In this range the time to evolve the best-adapted geno-
type can be reduced by many orders of magnitude com-
1A μ 2A μ 3A Genotype pared to the case without switching. Finally, we show
that if the switching rate is allowedto evolve it will con-
Initial Target
verge to values within this optimal range.
B
1-μ
Replication MODEL
r
μ
Mutation We consider a population of haploid cells. Each cell
α can be in one of two phenotypic states, A and B, and
Phenotype can assume one of three genotypes (“genetic states”), as
switch shown in Fig. 1A,B. Cells replicate stochastically with
Genotype d
rate r (1−N/K), where i labels one of the six possible
space i
Death states, N is the totalpopulationsize andK is the carry-
ing capacity of the environment. The logistic-like factor
C
(1−N/K) causes growth to cease when N becomes as
S largeasK,whichlimitsthetotalpopulationsize. During
1(cid:0)
r phenotype A replication a cell will produce a mutated offspring of an
e phenotype B adjacentgenotypewithprobabilityµ. Mutationdoesnot
t
a
r change the phenotypic state. All cells switch randomly
wth (cid:2)c between the states A and B at the symmetric switching
o 1(cid:1) rate α. Cells are randomly removed from the popula-
r
G
δ tion at a constant death rate d. If the carrying capacity
1(cid:1) K is small (K = 10···104) the model is appropriate to
1 2 3 describe a small microbial population growing in a mi-
Genotype crofluidic chemostat with constant dilution rate [5]. For
larger K (K =104···109) the model is relevant to pop-
Figure 1: The model. (A) Diagram showing the six pos- ulations cultured in mesoscopic (cm-size) chemostats.
sible states of a cell and the available transitions between
The population initially consists of all individuals of
them. The genotypes are labelled 1, 2 and 3, the pheno-
type 1A, i.e. of genotype 1 in phenotypic state A, which
typic states are labelled A and B. Transitions between geno-
types/phenotypes occur at rates µ and α, respectively. All hasthegrowthrater1A =1(arbitraryunits). Theinitial
cellsareinitiallyinstate1A.Evolutioncontinuesuntilasingle population size, unless otherwise stated, is equal to the
cell reaches the target state 3A. (B) Each cell can replicate, equilibrium size (1−d)K. State 3A is the global max-
switch phenotype, or die with rates b,α and d respectively. imum of both fitness landscapes with the growth rate
Upon replication a cell has the probability µ of producing a r3A =1+S, where S >0 is its selective advantage over
mutantofeachneighbouringgenotype. (C) Thefitnessland-
the ‘wild-type’ 1A (Fig. 1C). Such a beneficial mutant
scapesforbothphenotypes. PhenotypeAhasafitnessvalley
has the probability of fixation, from a single individual
at 2A while phenotypeB has uniform fitnessacross all geno-
in a population of cells in state 1A, approximately given
types.
by 1−e−S ≈ S [34] when 1/K ≪ S ≪ 1. We are in-
1−e−KS
terested in the adaptation time T that it takes for the
fitness toatleastthatofthe ancestralstrain[30,38,39]. populationtoevolvethefirstindividualinstate3A,con-
Here we study evolution in the presence of such a fit- ditioned on the population not going extinct (otherwise
ness valley, in which a population must acquire two mu- the time would be infinite).
tations to reach the best-adapted genotype. The first We shall begin by considering the case in which the
mutation is deleterious and corresponds to a valley in growth rate r2A = 1−δ, while the growth rates r1B =
the fitness landscape that is crossed by a second, com- r2B = r3B = 1, which leads to a “fitness valley” in phe-
pensatory mutation. Theory and computer modelling notype A and a flat fitness landscape in phenotype B
suggeststhat fitness valleys significantly affect biological (Fig.1C).Intheabsenceofphenotypeswitchingtheonly
evolution[9,42,47,48], often,butnotalways,slowingit waya population in state 1A canevolveanindividual in
down. However, in our model cells can also switch to an state 3A is to go through the low-fitness state 2A. How-
alternative phenotype in which the fitness landscape is ever, phenotype B allows the population to traverse an
3
alternative route that avoids state 2A. A
109 α = 0
-5
α = 10
108
RESULTS
T
107
Optimal range of the switching rate exists in which
evolution is fastest 106
105
Figure 2A shows simulation results examining the ef-
fect of phenotype switching on the adaptation time T
μ=10-5 =10-4 δ=0.6 K=106
for a range of the parameters of the model. In all cases
the presence of low frequency switching, i.e. small α, B109
=1x10-5
decreases the adaptation time by orders of magnitude. 108
Holdingallotherparametersfixedandchangingonlythe =3x10-5
107
switchingraterevealsanoptimalrangeforα(Fig.2B)in T =9x10-5
whichT isminimizedprovidedµissufficientlysmall. For 106
the smallest mutation probability (µ=10−5) considered =2.7x10-4
105
inFig.2Bthe minimaladaptationtimeisapproximately =8.1x10-4
two orders of magnitude lower than in the absence of 104 =0.00243
switching (α=0), or when switching is frequent (α very 103 =0.00729
large). Only if the mutation rate is unrealistically large
10-9 10-7 10-5 10-3 10-1 10
α
does switching not speed up evolution (SI and Fig. S1).
Intheabsenceoftransitionsbetweenstates2Aand2B
(i.e. removing this step from the model in Fig. 1A) the Figure 2: Adaptation time T for different scenarios – note
logarithmicscale. (A)AbarchartcomparingpairsofT values
adaptation time decreases monotonically with α (SI and
withandwithoutswitchingphenotypes(α=10−5 andα=0
Fig. S2).
respectively) for different parameter values. Left-most pair
of bars: K = 100,µ = 10−5,δ = 0.4 and d = 0.1. A label
underneatheachpairofbarsindicateswhichvariablehasbeen
Fastest evolutionary trajectories avoid the valley changed compared to the left-most pair. (B) T as a function
of the switching rate α for a range of mutation probabilities
µ. Parameters as in (A). For small enough µ an optimal
To understand the evolutionary trajectory selected in
(minimizing T) switching rate can beseen.
the optimal range of switching frequency we examined
the histories of successful cells, i.e. the states from Fig.
1A visited during evolution from state 1A to the final
state 3A. We then represented each cell’s trajectory in
tiontime forlargeαinFig. 2B.Asimilardecomposition
the state space as a sequence of links connecting the vis-
into three regions can be obtained from the probability
ited states (Fig. 3A). These were then classified as one
of individual transitions occurring between the states of
of 21 classes by grouping sequences with multiple (back-
the system (Fig. S3).
and-forth) transitions between the states together with
those without (Fig. 3A). Fig. 3B shows the most proba- We next examined which trajectories would be pre-
bletrajectoryasafunctionofµandα. Asexpected,sim- ferredif wevariedK insteadofµ, asthe carryingcapac-
ulated evolution favours different trajectories depending ity is more likely to vary significantly in real microbial
on the values of µ,α. The (µ,α)-space can be approxi- populationsandis mucheasierto controlexperimentally
mately separated into three regions, each corresponding than the mutation rate. Since we ascertained earlier (cf.
to one or more trajectory classes. Region 1 corresponds the previous paragraph) that the time to adaptation is
to trajectories that go through the fitness valley at 2A. mostsignificantlyaffectedby whetherthe trajectoryvis-
Inthisregionmutationsarefrequentenoughtooffsetthe its state 2A or not, we further reduced the number of
loss of fitness incurred when passing through the valley. trajectory classes to just two (Fig. 3C): the “u” (upper)
In region 2 the dominant trajectory avoids the valley by trajectory type that avoids state 2A or the “v” (valley)
switchingtophenotypeB.Thisregioncorrespondstothe trajectory type that visits it. We found a large region
fastest adaptation times from Fig. 2B. Region 3 is char- in the space (K,α) that favours trajectories that avoid
acterisedbyamixtureoftrajectorytypesusingbothphe- the state 2A (Fig. 3D). This region extends to large
notype B and visiting state 2A. This region corresponds values of K ∼ 109. The adaptation time depends non-
to large α which not only enables transitions to pheno- monotonicallyonthepopulationsizeK(Fig. S4);similar
typeBbutalsofromstate2Btothedeleteriousstate2A. non-monotonic behaviour has been seen in models with-
The latter is responsible for the increase in the adapta- out phenotype switching [13, 32].
4
A C
Example trajectories
u (upper) trajectory: v (valley) trajectory:
avoids state 2A visits state 2A
D Probability state 2A is visited
Trajectory
class
α
B α
1.0
10-3
10-2
10-5
10-4
10-7
10-6
10-5 10-4 10-3 10-2 10-1 μ 10 103 105 107 109K
Figure 3: Trajectories of successful cells in genotype/phenotype space. (A) Diagram showing how different trajectories that
go through thesame set of states are grouped together intothesame trajectory class. This class is represented asa symbolin
which bluelinescorrespond totransitions madebythesuccessful cell. (B) Themost probabletrajectory class as afunction of
µ and α, for K = 100. Three regions labelled 1,2,3 can be distinguished. (C) Trajectories can be further grouped regarding
whether they avoid (trajectory “u”) or visit (trajectory “v”) the fitness valley at state 2A. (D) The most common trajectory
group (symbols as in panel C) and the probability of a trajectory visiting state 2A (colours, see the colour bar) as a function
of switching rate α and carrying capacity K, for µ=10−6. In all simulations δ=0.4 and d=0.1.
Results are robust to small fitness costs of The same optimal switching rate range as in Fig. 2
phenotype switching can be seen in Fig. 4 for costs c ≤ 0.05 of state 1A’s
fitness, although it reduces in significance as c increases.
This is due to the increased time taken by trajectories
So far we have assumed that the alternative pheno-
that avoid state 2A as a result of the reduced growth
type B has the same fitness as the wild-type state 1A.
rates of the phenotype B states. For large c, such as
We shall now consider the case in which phenotype B is
c=0.1inFig. 4,thereisnooptimalswitchingrangebut
less adapted to the environment than the state 1A. This
the adaptation time T decreases monotonically with α.
is a common scenario; for example, phenotypes more re-
Stochasticphenotypeswitchingthusspeedsupevolution
sistanttoantibioticsoftenhavealowergrowthratethan
even in the presence of (moderate) fitness costs.
susceptible phenotypes do in the absence of the drug
[4]. To model this, we assume that the growth rates
r1B = r2B = r3B = 1−c, where c > 0 is the fitness cost
of switching to phenotype B (see Fig. 1C). A cell that Switching rate evolves until it reaches the optimal
switches from 1A to 1B grows more slowly but switch- range
ing can still be advantageous if the alternative route by
phenotype B is faster than crossing the fitness valley at Sofarwehaveshownthatthereoftenexistsanoptimal
state 2A. Fig. 4 shows the adaptation time T as a func- range for the switching rate that results in the shortest
tion of switching rate α for different fitness costs. For adaptation time. This range coincides with successful
sufficiently small c the same qualitative behaviour as for cells avoiding the deleterious state 2A by switching to
c=0 is observed. phenotype B. To see whether a variable switching rate
5
1010 space at α values below the optimal α range. However,
the plot of T in Fig. 5C suggests this will be a slow pro-
cess. The latter observation can be explained through
109
these cells first evolving a large switching rate, allow-
ing them to switch quickly to phenotype B, before their
T
switching rate decreases again. This picture is corrobo-
108
= 0.1 ratedbyobservingthatatlowαthereislittlephenotype
c = 0.05
switching (purple links in Fig. 5B) and yet the likeli-
c= 0.02
107 c = 0.01 hoodofmutationsinphenotypeBacrossgenotypespace
c = 0 is large (green links in Fig. 5B). This provides an alter-
c
native efficient trajectory for a population to cross the
10-7 10-6 10 10-4 10 10-2 10-1 α
-(cid:3) -(cid:4) fitness valley. However, such trajectories need at least
one extra mutation compared to those that evolve α di-
Figure 4: Adaptation time in the presence of fitness costs.
rectly to within the optimal range. Therefore evolution
The adaptation time T as a function of the switching rate α
is dominated by the latter trajectories, in particular for
for different fitness costs c = {0,0.01,0.02,0.05,0.1} of phe-
notypeBandforparameters K =100,µ=10−5,δ=0.9 and low mutation probabilities µ.
d=0.1.
DISCUSSION
would evolve to be in this optimal range, thus optimis-
ing the evolutionary process, we simulated a population Phenotypic plasticity has been recently shown to af-
of cells that started at 1A with α ≈ 0 initially. The fectDarwinianevolutionofanimals[15]. Ourtheoretical
switching rate was subsequently allowed to evolve: dur- researchsuggeststhatstochastic phenotypic switchingis
ing replication, offspring were assigned a new α value also able to significantly speed up biologicalevolution of
with probability µ, the same as for mutations between microbes,perhapsbymanyordersofmagnitude. Switch-
genotypes 1 ↔ 2 ↔ 3. The new α was randomly and ing phenotypes enables cells to evade a fitness valley
uniformly selected from a fixed set of possible values whichmaybedifficulttocrossotherwise,particularlyfor
(α = 2.56×10−10 ×5i for integer i ∈ [0,15]); similar small populations. Although our model uses only three
resultswereobtainedwhenαevolvedincrementally(Fig. genotypes,webelievetheresultspresentedherearemuch
S5). In this extended model the evolutionary trajectory moregeneralandapplytorealistic,multi-dimensionalfit-
of successful cells spans three dimensions: two for the ness landscapes [36, 43, 46]. Epistasis causes such land-
genotype, one of which evolves the switching rate, and scapes to contain many local fitness maxima separated
one for the two phenotypes A and B (Fig. 5A). The end by genotypes with lower growth rate [24]. In a small
point of the trajectory remains a cell produced at state population, in the absence of phenotype switching, evo-
3A regardless of its switching rate. We collected tra- lution can spend much time in a local maximum before
jectories of successful cells in the expanded state space, “tunnelling” through one of the fitness valleys. Pheno-
including states with different α, and calculated transi- typic switching opens up a second fitness landscape in
tionprobabilitiesbetweenanytwostates. Fig. 5Bshows which minima and maxima can be located at different
these probabilities as links ofdifferent thicknessbetween genotypes compared to the first landscape. Phenotype
the states instate space. We observethat successfultra- switching provides alternative routes between otherwise
jectories are unlikely to feature genotype mutations in isolatedfitness maximain both landscapes,andthus en-
phenotypeA(veryfewredlinks). Instead,mostsuccess- ables evolution to proceed faster. This mechanism, sim-
fulcellsfirstswitchtostate1B,mutateto3B,andswitch ilar to the proposed interplay between genetic and epi-
back to state 3A. The evolved switching rate often falls genetic mutations [23], adds another dimension to the
within the optimal range found in Fig. 2B for the same complexityofevolutionarypathways[36,37,46]whichis
set of parameters K,µ and d. This is further illustrated relevant for the predictability of evolution [10].
in Fig. 5C where we overlay the adaptation time from We now discuss how plausible our computational re-
Fig. 2B with a barchartthat showsthe probabilitythat sults are for real microbial populations. Although our
a particular α was selected by evolution. theory is quite general, we focus on antibiotic resistance
Fig. 5C reveals two more interesting features. First, evolution. We believe our model represents best a sit-
we see that cells that have evolved α above the opti- uation in which a population of microbes is exposed to
mal range (α ≥ 0.02) are unlikely to be successful due sub-lethal concentrations of an antibiotic. In this sce-
to fasttransitionsbetweenstates2Band2A(fitness val- nario, state 1A corresponds to a ‘wild-type’ genotype
ley). Suchtransitionseffectivelylowerthefitnessofstate whose growth rate is slightly lowered by the presence
2B and create a new fitness valley at this state. Second, of the antibiotic, but cells are still able to reproduce.
weseethatalargenumberoftrajectoriescrossgenotype The fitness landscape (Fig. 1C) for phenotype A then
6
A
Individual trajectories collected Combined trajectories
Phenotype
....
SPS rate α Genotype
Step in phenotype A
Step in phenotype B
Initial state
Step switching phenotypes
Each α value has the 6 Step evolving α
B genotype/phenotype C
states available 3A
probability of crossing genotype space
1A
0.2 (cid:6) 109
10 y
t
10-1 bili0.15 108
a
b (cid:5)
10-3 ro 0.1 107
P
10-5
B 0.05 106
10-7
A
3
10-9 2 10-9 10-7 10-5 10-3 10-1 10
1 α
Figure 5: Evolution selects switching rates within the optimal range. (A) Left: examples of evolutionary trajectories in the
extended model in which the stochastic phenotype switching (SPS) rate α also evolves. (A) Right: trajectories are used
to calculate transition probabilities between the states of the system, which are then represented by the thickness of links
connecting the states. Red links correspond to mutations in phenotype A, green links to mutations phenotype B and purple
links to switching between phenotypes. (B) Graph of transition probabilities where line thicknesses are proportional to the
probability that a successful trajectory will take that step. The parameters are K = 100,µ= 10−5,δ = 0.4 and d = 0.1, the
samethatwereusedinFig. 2B(blueline). Thepopulationbeginsatthewild-type1Awithα=2.56×10−10 andevolvesuntil
a cell in state 3A is produced. (C) The probability that genotype space is crossed at a given α in either phenotype A or B.
TheprobabilityhasamaximumwheretheadaptationtimeT forfixedα(plotfrom Fig. 2Bsuperimposedonthesamegraph)
has its minimum.
corresponds to a situation in which an initial mutation tation is 10−10 − 10−9 per replication [11, 27] but can
further lowers the growthrate by e.g. making the target be higher (10−7) in mutator strains [27, 40] or in the
enzyme less susceptible to the drug but alsoless efficient presence of antibiotics [18]. Assuming µ = 10−6 which
[31,38,39],whileasecondmutationcompensatesforthe would correspond to ∼10 point mutations in a mutator
loss of efficiency and increases the growth rate beyond strain or any loss-of-functionmutation in a non-mutator
thatoftheancestralstrain. Inthesameexampletheflat strain (assuming a typical gene length 1kbp [22]), Fig.
landscape for phenotype B can result from switching to 3Dindicatesthatphenotypeswitchingwouldbe thepre-
a state with higher gene expression of multidrug efflux ferredmode of evolutionfor α<10−4 in a population of
pumps or antibiotic-degrading enzymes. This has been K =103 cells, for 10−7 <α<5×10−3 for intermediate
observed for a biofilm-forming bacterium P. aeruginosa size K = 106, and for α ≈ 10−3 for K = 109. For lower
whichcanswitchtoaresistantphenotypeinthepresence mutation probabilities µ we expect this behaviour to ex-
of antibiotics [8, 12]. Alternatively, phenotype B can be tendtolowerswitchingrates,consistentwiththeincreas-
a persister - although the conventional view is that per- ingoptimalrangesseeninFig. 2B,whilemaintainingap-
sistent cells do not replicate,this has been recently chal- proximatelythe same upper limit ofα≈10−3,...,10−2.
lenged [45]. Slower (but non-zero) replication rates can Switching rates encounteredin nature can be as large as
enable persisters to mutate and explore a larger set of this [3, 44] and thus we believe that stochastic pheno-
genotypes than the non-persister (phenotype A) would typic switching can be a significant feature affecting the
be able to do in the absence of switching. biological evolution of microbes.
In microbes, the probability of single nucleotide mu- The model we have proposed is a mathematical ideal-
7
isation. Although it may be possible to test the quan-
titative predictions of our model in a carefully designed
experiment, many of the assumptions we made would
[1] Martin Ackermann. A functional perspective on pheno-
have to be relaxed to model naturally occurring popu-
typic heterogeneity in microorganisms. Nature Reviews
lations of bacteria. However, recent work suggests that Microbiology, 13(8):497–508, August 2015.
mechanismssimilartowhatwedescribedheremaybeat [2] MartinAckermann,BrbelStecher,NikkiE.Freed,Pascal
play in real microbial populations. In particular, Refs. Songhet,Wolf-DietrichHardt,andMichaelDoebeli.Self-
[7, 51] show that the bacterium E. coli switches to a destructive cooperation mediated by phenotypic noise.
filamentous phenotype (“phenotype B”) as a result of Nature, 454(7207):987–990, August 2008.
[3] Aileen M Adiciptaningrum, Ian C Blomfield, and
exposure to sub-inhibitory concentration of the antibi-
Sander J Tans. Direct observation of type 1 fimbrial
oticciprofloxacin,andfilamentouscellsoccasionallypro-
switching. EMBO reports, 10(5):527–532, May 2009.
duce normal-length cells (“phenotype A”) resistant to
[4] Nathalie Q. Balaban, Jack Merrin, Remy Chait, Lukasz
ciprofloxacin. It has been hypothesized that phenotype Kowalik, and Stanislas Leibler. Bacterial Persistence
B providesa “safe niche” giving bacteriaenoughtime to as a Phenotypic Switch. Science, 305(5690):1622–1625,
evolve resistance. Although the switching between the September2004.
phenotypesinthisexperimentisnotentirelyrandombut [5] FrederickK.Balagadd´e,LingchongYou,CarlL.Hansen,
Frances H. Arnold, and Stephen R. Quake. Long-
a response to stress,the observedbehaviour is similar to
Term Monitoring of Bacteria Undergoing Programmed
that predicted by our model.
Population Control in a Microchemostat. Science,
309(5731):137–140, July 2005.
[6] HubertusJ.E. Beaumont,JennaGallie, Christian Kost,
Gayle C. Ferguson, and Paul B. Rainey. Experimen-
METHODS
tal evolution of bet hedging. Nature, 462(7269):90–93,
November2009.
Custom written Java programs were used to simulate [7] JuliaBos,QiucenZhang,SaurabhVyawahare,Elizabeth
the model. For smaller K (K ≤ 100) each cell was as- Rogers, Susan M. Rosenberg, and Robert H. Austin.
signed a current type (one of the 6 possible states from Emergence of antibiotic resistance from multinucleated
bacterialfilaments. Proceedings oftheNational Academy
Fig. 1A) and a sequence of past states beginning at
of Sciences, 112(1):178–183, January 2015.
state 1A. A variant of the Gillespie kinetic Monte Carlo
[8] Elena B. M. Breidenstein, C´esar de la Fuente-Nu´n˜ez,
algorithm[17] was used to evolve the system. In each
and Robert E. W. Hancock. Pseudomonas aeruginosa:
time step, a random cell and action (replication, death, all roads lead to resistance. Trends in Microbiology,
orphenotypeswitch)wasselectedwithaprobabilitypro- 19(8):419–426, August 2011.
portional to the rate with which that action occurs. A [9] Arthur W. Covert, Richard E. Lenski, Claus O. Wilke,
second random number drawn from an exponential dis- and Charles Ofria. Experiments on the role of delete-
rious mutations as stepping stones in adaptive evolu-
tribution, with mean equal to the inverse of the total
tion. Proceedings of the National Academy of Sciences,
rate of all possible actions in the system, was used as
110(34):E3171–E3178, August 2013.
a measure of time elapsed during that step. For larger
[10] J. Arjan G. M. de Visser and Joachim Krug. Empiri-
K > 100 (or when c > 0) the evolution of the system calfitnesslandscapesandthepredictabilityofevolution.
was approximately modelled using a type of τ−leaping Nature Reviews Genetics, 15(7):480–490, July 2014.
algorithm [16]. [11] J. W. Drake. A constant rate of spontaneous muta-
Statistical analysiswas performedin Mathematica us- tioninDNA-basedmicrobes. ProceedingsoftheNational
Academy of Sciences, 88(16):7160–7164, August 1991.
ing a custom-written code. Mean values presented in all
plots were obtained by averaging over 104 runs, except [12] Eliana Drenkard. Antimicrobial resistance of Pseu-
domonas aeruginosa biofilms. Microbes and Infection,
for Fig. 3D where points were averaged over 102 −103
5(13):1213–1219, November2003.
runs. Error bars are s.e.m. [13] S.FElena,C.OWilke,C.Ofria,andR.ELenski.Effects
of population size and mutation rate on the evolution of
mutational robustness. Evolution, 61(3):666–674, 2007.
[14] J. Franke, A. Klo¨zer, J. A.G.M de Visser, and J. Krug.
Acknowledgments
Evolutionaryaccessibility ofmutationalpathways. PLoS
computational biology, 7(8):e1002134, 2011.
We thank Rosalind Allen for helpful discussions and [15] Cameron K. Ghalambor, Kim L. Hoke, Emily W. Ru-
ell, Eva K. Fischer, David N. Reznick, and Kimberly a.
many suggestions at the beginning of this project. BW
Hughes. Non-adaptiveplasticity potentiatesrapidadap-
was supported by a RSE/Scottish Government Personal
tive evolution of gene expression in nature. Nature,
Research Fellowship. MRE and AT thank the EPSRC
September2015.
for support through grant number EP/J007404/1and a [16] D.T.Gillespie. Approximateacceleratedstochasticsim-
studentship respectively. This work has made use of the ulation of chemically reacting systems. The Journal of
resourcesprovidedbythe EdinburghCompute andData Chemical Physics, 115(4):1716–1733, 2001.
Facility (ECDF) (http://www.ecdf.ed.ac.uk/). [17] Daniel T. Gillespie. Exact stochastic simulation of cou-
8
pled chemical reactions. The journal of physical chem- pathogens, 5(8):e1000541, August 2009.
istry, 81(25):2340–2361, 1977. [32] O. C Martin, L. Peliti, et al. Population size effects in
[18] Stephen H. Gillespie, Shreya Basu, Anne L. Dickens, evolutionarydynamicsonneutralnetworksandtoyland-
Denise M. O’Sullivan, and Timothy D. McHugh. Effect scapes. JournalofStatisticalMechanics: TheoryandEx-
of subinhibitory concentrations of ciprofloxacin on My- periment, 2007:P05011, 2007.
cobacterium fortuitum mutation rates. Journal of An- [33] Oliver G. McDonald, Hao Wu, Winston Timp, Akiko
timicrobial Chemotherapy, 56(2):344–348, August 2005. Doi, and Andrew P. Feinberg. Genome-scale epigenetic
[19] H. Hasman, M.A. Schembri, and P. Klemm. Antigen 43 reprogramming duringepithelial-to-mesenchymaltransi-
andtype1fimbriaedeterminecolony morphology ofEs- tion. Nature Structural & Molecular Biology, 18(8):867–
cherichiacoliK-12. Journal of bacteriology, 182(4):1089, 874, August 2011.
2000. [34] MartinA.Nowak.EvolutionaryDynamics. HarvardUni-
[20] Daniel B. Kearns and Richard Losick. Cell population versity Press, September2006.
heterogeneity during growth of Bacillus subtilis. Genes [35] Ertugrul M. Ozbudak, Mukund Thattai, Han N. Lim,
& Development, 19(24):3083–3094, December 2005. Boris I. Shraiman, and Alexander van Oudenaarden.
[21] Kristel Kemper, Pauline L. de Goeje, Daniel S. Peeper, Multistability in the lactose utilization network of Es-
andRen´eevanAmerongen. Phenotypeswitching: tumor cherichiacoli.Nature,427(6976):737–740,February2004.
cell plasticity as a resistance mechanism and target for [36] Adam C. Palmer, Erdal Toprak, Michael Baym, Se-
therapy. Cancer Research, 74(21):5937–5941, November ungsoo Kim, Adrian Veres, Shimon Bershtein, and Roy
2014. Kishony. Delayed commitment to evolutionary fate in
[22] Michael J. Kerner, Dean J. Naylor, Yasushi Ishihama, antibiotic resistance fitness landscapes. Nature Commu-
TobiasMaier,Hung-ChunChang,AnnaP.Stines,Costa nications, 6, June2015.
Georgopoulos, Dmitrij Frishman, Manajit Hayer-Hartl, [37] Frank J. Poelwijk, Daniel J. Kiviet, and Sander J.
Matthias Mann, and F. Ulrich Hartl. Proteome-wide Tans. EvolutionaryPotentialofaDuplicatedRepressor-
analysis of chaperonin-dependent protein folding in Es- Operator Pair: Simulating Pathways Using Mutation
cherichia coli. Cell, 122(2):209–220, July 2005. Data. PLoS Computational Biology, 2(5):e58, 2006.
[23] Filippos D. Klironomos, Johannes Berg, and Sin´ead [38] SJSchragandVPerrot. Reducingantibioticresistance.
Collins. Howepigeneticmutationscanaffectgeneticevo- Nature, 381(6578):120–121, May 1996.
lution: Model and mechanism. BioEssays, 35(6):571– [39] Stephanie J Schrag, V´eronique Perrot, and Bruce R
578, 2013. Levin. Adaptation to the fitness costs of antibi-
[24] R.Korona, C.H. Nakatsu,L.J. Forney,and R.E. Lenski. otic resistance in Escherichia coli. Proceedings of the
Evidence for multiple adaptive peaks from populations Royal Society of London. Series B: Biological Sciences,
of bacteria evolving in a structured habitat. Proceedings 264(1386):1287–1291, September1997.
oftheNational AcademyofSciences oftheUnitedStates [40] PDSniegowski,PJGerrish,andRELenski. Evolution
of America, 91(19):9037–41, September1994. ofhighmutationratesinexperimentalpopulationsofE.
[25] Edo Kussell, Roy Kishony, Nathalie Q. Balaban, and coli. Nature, 387(6634):703–705, June 1997.
Stanislas Leibler. Bacterial Persistence A Model of Sur- [41] Ana Solopova, Jordi van Gestel, Franz J. Weissing,
vivalinChangingEnvironments. Genetics, 169(4):1807– Herwig Bachmann, Bas Teusink, Jan Kok, and Os-
1814, April 2005. car P. Kuipers. Bet-hedging during bacterial diauxic
[26] Edo Kussell and Stanislas Leibler. Phenotypic Diver- shift. Proceedings of the National Academy of Sciences,
sity,PopulationGrowth,andInformationinFluctuating 111(20):7427–7432, 2014.
Environments.Science,309(5743):2075–2078, September [42] IvanG.Szendro,JasperFranke,J.ArjanG.M.deVisser,
2005. and Joachim Krug. Predictability of evolution depends
[27] HeewookLee,EllenPopodi,HaixuTang,andPatriciaL. nonmonotonically on population size. Proceedings of the
Foster.Rateandmolecularspectrumofspontaneousmu- National Academy of Sciences, 110(2):571–576, January
tations in the bacterium Escherichia coli as determined 2013.
bywhole-genomesequencing.ProceedingsoftheNational [43] Ivan G. Szendro, Martijn F. Schenk, Jasper Franke,
Academy of Sciences of the United States of America, Joachim Krug, and J. Arjan G. M. de Visser. Quan-
109(41):E2774–E2783, October2012. titative analyses of empirical fitness landscapes. Jour-
[28] Bruce R. Levin and Daniel E. Rozen. Non-inherited nal of Statistical Mechanics: Theory and Experiment,
antibiotic resistance. Nature Reviews Microbiology, 2013(01):P01005, January 2013.
4(7):556–562, 2006. [44] Marjan W van der Woude and Ian R Henderson. Reg-
[29] H´edia Maamar and David Dubnau. Bistability in the ulation and function of Ag43 (flu). Annual Review of
Bacillus subtilis K-state (competence) system requires a Microbiology, 62:153–169, 2008.
positive feedback loop: Bistability in B. subtilis com- [45] Y. Wakamoto, N. Dhar, R. Chait, K. Schneider,
petence. Molecular Microbiology, 56(3):615–624, March F. Signorino-Gelo, S. Leibler, and J. D. McKinney. Dy-
2005. namic Persistence of Antibiotic-Stressed Mycobacteria.
[30] Sophie Maisnier-Patin, Otto G. Berg, Lars Liljas, and Science, 339(6115):91–95, January 2013.
DanI.Andersson. Compensatoryadaptationtothedele- [46] D.M.Weinreich. Darwinian EvolutionCan FollowOnly
terious effect of antibiotic resistance in Salmonella ty- Very Few Mutational Paths to Fitter Proteins. Science,
phimurium. Molecular Microbiology, 46(2):355–366, Oc- 312(5770):111–114, April 2006.
tober 2002. [47] D.M.WeinreichandL.Chao. Rapidevolutionaryescape
[31] Linda L Marcusson, Niels Frimodt-Moller, and Di- by large populations from local fitness peaks is likely in
armaid Hughes. Interplay in the selection of flu- nature. Evolution, 59(6):1175–1182, 2005.
oroquinolone resistance and bacterial fitness. PLoS [48] D.B. Weissman, M.M. Desai, D.S. Fisher, and M.W.
9
Feldman. The rate at which asexual populations cross Adaptation time decreases monotonously with α in
fitnessvalleys. Theoretical population biology,75(4):286– the absence of switching between states 2A and 2B
300, 2009.
[49] Zhizhong Yao, Rebecca M. Davis, Roy Kishony, Daniel
We identified that region 3 in Fig. 3B appears due
Kahne, and Natividad Ruiz. Regulation of cell size in
to the existence of the step between states 2A and 2B.
responsetonutrientavailabilitybyfattyacidbiosynthesis
inEscherichiacoli. Proceedings of the National Academy For large α values a trajectory via phenotype B is likely
of Sciences, 109(38):E2561–E2568, September2012. to involve switching to the deleterious state 2A. This re-
[50] StefanoZapperiandCaterinaA.M.LaPorta. Docancer ducestheaveragefitnessofthepopulationandcausesthe
cellsundergophenotypicswitching? Thecase forimper- adaptationtimetoincreasewithα(Fig. 2B).Todemon-
fectcancerstemcellmarkers. Scientific Reports, 2,June
stratethisconsiderthemodelwithoutswitchingbetween
2012.
states 2A and 2B (Fig. S2A). The adaptation time T as
[51] Q. Zhang, G. Lambert, D. Liao, H. Kim, K. Robin, C.-
afunctionofswitchingrateαdecreasesmonotonicallyin
k. Tung, N. Pourmand, and R. H. Austin. Accelera-
tion of Emergence of Bacterial Antibiotic Resistance in this modified model (Fig. S2B).
ConnectedMicroenvironments.Science,333(6050):1764–
1767, September2011. A
6 state system of genotype/phenotype pairs
μ μ
2B
(cid:13)(cid:12) 3(cid:12) Phenotype
SUPPLEMENTARY INFORMATION
α α
Phenotypic switching will not always speed up
evolution
μ μ
2A Genotype
(cid:13)(cid:14) 3(cid:14)
In Fig. 2A we compared the adaptation time T with Initial Tar
g(cid:10)(cid:11)
and without phenotypic switching for various parame-
B
ters. In all cases switching decreased T. However, if the
mutationprobabilityµislargeenough,thepathwaythat (cid:17)
108 =2x104
goes directly through the valley genotype 2A may take K=2x106
lesstimethananypathwayinvolvingthealternativephe- K=2x108
K
notype B. Figure S1 shows T as a function of mutation
probability µ for the case with and without switching. 106
For large enough µ, switching makes no observable dif-
ference to the adaptation time as successful cells do not
use an alternative pathway. This is the same behaviour
104
observed in region 1 of Fig. 3B.
10 10-6 10-4 10-2 α
(cid:15)(cid:16)
109 α = 0 Figure S2: Time to adaptation without phenotypeswitching
α = 10-4 between states 2A and 2B. (A) The modified model. (B)
Adaptationtime T asafunction of switching rateαfor K =
107 2×104,2×106 and 2×108. Dashed lines indicate T when
α = 0. Remaining parameters are µ = 10−6,δ = 0.4 and
(cid:9) d=0.1
105
103
Probability of individual steps for successful
10-5 10 μ 10-3 10-2 trajectories
(cid:7)(cid:8)
FigureS1: Phenotypicswitchingdoesnotspeedupevolution In Fig. 3B we showed the most common trajectory
for large mutation probabilities µ. Shown is the adaptation classesforsuccessfulcellsasafunctionof(µ,α). Amore
timeT fromcomputersimulationsasafunctionofµwithand
detailed picture can be obtained by looking at the prob-
without phenotype switching (α = 10−4 and α = 0 respec-
ability that each step occurs in a successful cell’s tra-
tively),for K =100,δ=0.4 and d=0.1.
jectory. For each pair of (µ,α) we create a diagram of
transitionsbetweenstates(i.e. 1A,1B,...,3A,3B)inwhich
10
A non-monotonic behaviour can be seen: T is maximal for
intermediateK andfallsoffforbothsmallandlargecar-
. rying capacities.
N .
Simulation . (cid:19) α=8.192x10-6
1010
runs α=2.62144x10-4
Combined α=8.388608x10-3
trajectories 108
Individual trajectories
collected 106
B α
104
1.0
102 104 106 108
(cid:18)
10-4 FigureS4: Adaptationtimeasafunctionofcarryingcapacity
K for different values of α. Remaining parameters are µ =
10−6,δ=0.4 and d=0.1
10-2
Evolving the switching rate in fixed increments
10-6
Fig. 5 showed simulation results in which newly cre-
10-5 10-4 10-3 10-2 10-1 μ
atedcellscouldmutatetheirswitchingrateαbydrawing
anewvaluefromadiscretesetofallowedswitchingrates.
Figure S3: The probability of each step being taken by suc-
We observedsuccessfulcelltrajectoriesto consistlargely
cessful cells. (A) Constructing images that display the prob-
of those that evolved switching rates to within the op-
ability thateach stepis taken. By collecting manysuccessful
timal range identified in Fig. 2B. Here we consider the
trajectory classes we can compile them into a single image
thatconsistsoflinksofdifferentthicknesses,wherethethick- same set of allowed α values but, upon mutation of α,
ness of each link is proportional to the probability that that the new value is chosen as one of the two closest values
step is taken by a successful cell. (B) The probability each to the switching rate of the ancestor cell. This corre-
stepistakenasafunctionofµandα. Remainingparameters sponds to α evolving in fixed multiplicative steps, with
are K =100,δ=0.4 and d=0.1.
increaseordecreaseofαbeingequallylikely. Therestof
the algorithm is the same as before.
the thickness of each link corresponds to the probability Weagaincollectedtrajectoriesofsuccessfulcellsinthe
that that step is taken (Fig. S3A). Figure S3B shows a expanded state space, including states with different α,
graphinwhichdiagramsfordifferentpairsof(µ,α)have and calculated transition probabilities between any two
beenputtogether. ThisgraphisverysimilartoFig. 3B. states. Fig. S5B shows these probabilities as links of
This confirms that the most-common trajectories plot- different thickness between the states in state space. We
ted in Fig. 3B are good representatives of the “typical” can clearly see some new behaviour here. First, success-
trajectories of successful cells. ful trajectories are more likely to feature genotype mu-
tations inphenotype A(more redlines)thanin Fig. 5B.
Second,theprobabilityofcrossingthegenotypespaceat
Adaptation time as function of carrying capacity large α (in either phenotype) is much smaller due to the
increased number of steps it takes a cell to evolve such
Using the same simulation data as in Fig. 3D but large α values. Otherwise, Fig. S5 is similar to Fig. 5.
plotting it differently we canexplorehow the adaptation This shows that the details of how mutations of α are
time T depends on the carrying capacity K. Figure S4 implemented do not have a significant effect upon the
shows T for different values of the switching rate α. A results of our model.