Home

Search

Collections

Journals

About

Contact us

My IOPscience

Emergence of phenotype switching through continuous and discontinuous evolutionary transitions

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2015 Phys. Biol. 12 046004 (http://iopscience.iop.org/1478-3975/12/4/046004) View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address: 150.135.239.97 This content was downloaded on 02/10/2015 at 22:00

Please note that terms and conditions apply.

Phys. Biol. 12 (2015) 046004

doi:10.1088/1478-3975/12/4/046004

PAPER

RECEIVED

27 November 2014

Emergence of phenotype switching through continuous and discontinuous evolutionary transitions

REVISED

18 March 2015 ACCEPTED FOR PUBLICATION

Pintu Patra and Stefan Klumpp

8 April 2015

Max Planck Institute of Colloids and Interfaces, Science Park Golm, D-14424 Potsdam, Germany

PUBLISHED

E-mail: [email protected]

28 May 2015

Keywords: bacterial persistence, phenotype switching, population heterogeneity, bet-hedging, antibiotic tolerance

Abstract Bacterial persistence (phenotypic tolerance to antibiotics) provides a prime example of bet-hedging, where normally growing cells generate slow-growing but antibiotic-tolerant persister cells to survive through periods of exposure to antibiotics. The population dynamics of persistence is explained by a phenotype switching mechanism that allows individual cells to switch between these different cellular states with different environmental sensitivities. Here, we perform a theoretical study based on an exact solution for the case of a periodic variation of the environment to address how phenotype switching emerges and under what conditions switching is or is not beneficial for long-time growth. Specifically we report a bifurcation through which a fitness maximum and minimum emerge above a threshold in the duration of exposure to the antibiotic. Only above this threshold, the optimal phenotype switching rates are adjusted to the time scales of the environment, as emphasized by previous theoretical studies, while below the threshold a non-switching population is fitter than a switching one. The bifurcation can be of different type, depending on how the phenotype switching rates are allowed to vary. If the switching rates for both directions of the switch are coupled, the transition is discontinuous and results in evolutionary hysteresis, which we confirm with a stochastic simulation. If the switching rates vary individually, a continuous transition is obtained and no hysteresis is found. We discuss how both scenarios can be linked to changes in the underlying molecular networks.

1. Introduction In recent years it has become increasingly clear that heterogeneity in microbial populations is the rule rather than the exception. Even isogenic populations exhibit heterogeneous phenotypes, either via broad distributions of traits due to stochastic effects or, most importantly, via multiple distinct phenotypes that arise from bistability or multistability in the underlying gene networks [1, 2]. Hallmarks of such phenotypic heterogeneity are bimodal or multimodal gene expression patterns and stochastic switching between these phenotypes [3, 4]. Examples for this type of heterogeneity have been reported in metabolic systems such as central metabolism in yeast [5, 6] and Escherichia coli (E. coli) [7, 8] and the diauxic shift in Lactococcus lactis [9], in strategies for stress such as competence [10] and sporulation [11], in virulence [12] and in tolerance to antibiotics [3, 13]. © 2015 IOP Publishing Ltd

Several different roles of heterogeneity have been discussed, the most important ones being specialization and division of labor [12] on the one hand and bet-hedging on the other hand [5, 9, 14, 15]. In particular, the case of phenotypic tolerance to antibiotic (persistence) has been studied quite extensively in the last years as a prime example of phenotypic heterogeneity and bet-hedging, as one of the most well defined cases where stochasticity in gene expression networks has a clear beneficial role and due to its relevance to antibiotic treatment failures [3, 15–20]. The central idea of bet-hedging is to phenotypically diversify a population in order to be prepared for changes in environmental conditions. Thus, cells generate different phenotypes, some of which are maladapted to the present growth conditions, but provide an advantage when the conditions change, because they are adapted to the new condition. The phenotypes that are selected against by the current conditions are re-generated by

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

phenotype switching and provide a benefit to the whole population in the long run. Thus a bet-hedging strategy allocates some individuals in a population to a protective state as an ‘insurance policy’ to survive adverse environmental conditions [14, 15, 21]. In the case of bacterial persistence, the population benefits from ‘persister’ cells that are more tolerant when exposed to antibiotics [13, 22, 23], but grow slowly in the absence of antibiotics [19], thereby causing a growth reduction and thus a cost to the population in growth supporting conditions. In a changing environment, the interplay between cost and benefit determines the circumstances under which phenotype switching is advantageous. Several theoretical studies have been performed to understand the basic feature of phenotype switching mechanisms [15, 17, 18, 20, 24–26]. Using models for stochastically or periodically varying environments, these studies have shown that phenotypic heterogeneity is indeed advantageous in a varying environment. An aspect that received particular attention is the fact that optimal phenotype switching rates can be found that maximize the long-term growth rate [15, 17, 20, 25, 27]. These optimal switching rates are adapted to the characteristic time scales of environmental changes and indicate that heterogeneity can be evolved and tuned to a changing environment. However, two recent studies have provided evidence that this may not always be the case and rather than tuning the switching rates, a change in the environmental time scales may abolish heterogeneity altogether, possibly in an abrupt fashion [24, 27, 28]. Specifically, a threshold phenomenon was observed, where a marginal difference in selection pressure between the environmental phases selects very different evolutionary stable switching rates, i.e., either relatively frequent switching or extremely rare or even no switching [27]. Further, an analytical study explained this threshold phenomenon for the case of equal switching rates between the phenotypes and found that a non-switching population can be selected over a switching one when there is an asymmetry between two environmental phases in the selection pressure or in the duration [24]. Here we use a combination of an analytical model and evolutionary simulations towards a deeper understanding of the conditions under which phenotype switching is beneficial. We derive an exact expression for the long-term fitness of a heterogeneous population consisting of normal cells and persisters in a periodically changing environment. For this result we do not make use of assumption such as long periods of constant environment [15, 20] or symmetric switching rates [24]. Thus our result can be used to study the population dynamics and the evolutionary dynamics of bacterial persistence over the full range of cellular and environmental parameters in order to identify conditions which undermine the selection for phenotype switching in changing environments. Here, we 2

use these results to identify conditions for the existence of heterogeneity, focusing on the duration of the stress phases. Repeated exposure to stress such as antibiotic treatment is a primary facilitator for the selection of persistent phenotypes. For example, high persistence mutants of E. coli have been isolated through alternating exposure to normal growth conditions and antibiotic treatments [29, 30], and periodic application of high doses of antibiotics led to increased persister levels of Pseudomonas aeruginosa in cystic fibrosis patients [31]. Specifically, we show that a minimal duration is needed for the long-time benefit of phenotype switching. Surprisingly, we find that the transitions between different evolutionary stable states or strategies can be continuous or discontinuous, depending on how the parameters are allowed to vary. Correspondingly, the evolution of phenotype switching may exhibit hysteresis in some, but not all modes of evolution of cellular parameters.

2. Model and methods 2.1. Phenotype switching model We consider a bacterial population where an individual cell can have two distinct phenotypic states, normal and persister, which are characterized by different sensitivities to given environmental conditions. The environmental sensitivity is reflected in growth and death rates of the subpopulation in the given environment. The normal cells generate persister cells through a reversible phenotype switch [19]. A cell in the normal state can switch to the persister state with rate a and a cell in the persister state can switch back to the normal state with rate b. We denote the growth rate of normal cells (n) and persister cells (p) by μnm and μpm respectively, where m indicates the growth medium or, more generally the growth conditions. Below we will use indices ‘g’ and ‘s’ to denote unstressed growth and stress conditions, respectively (e.g., growth medium not containing or containing an antibiotic). The spontaneous switching between phenotypic states leads to distinct subpopulation of normal and persister cells which compose the total population. The subpopulation dynamics can be described by the following set of linear differential equations [19] dn = μnm n − an + bp , dt dp = μ pm p + an − bp . dt

(1)

The time evolution equation for the total population N (t ) = n + p can be written as ⎡ ⎤ d 1 N = ⎢ μnm + μ pm − μnm ⎥N. dt 1 + f (t ) ⎦ ⎣

(

)

(2)

The term in brackets is the time-dependent growth rate of the total populations, which depends on the

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

subpopulation ratio f (t ) = n (t ) p (t ). The average growth rate of the total population in a growth period tg can be computed as

(

μ Ng = μng + μ pg − μng

) t1 ∫ g

t 0+tg t0

1 dτ . (3) 1 + f (τ )

Note that this average depends on the initial conditions at t0 (via f (t 0 ) = f0 ). Below, where we consider periodically varying environments, this dependence will lead to a dependence on the duration of the previous environment. The expression for f(t) can be obtained by solving the differential equation for the time evolution of f(t) (equation (5) in [25]), and has the following form ⎛ ⎞ f0 − f ′ ⎟ −t Δ 2 + 4ab g e f ′ − f *⎜ ⎜ ⎟ ⎝ f0 − f * ⎠ , f (t ) = ⎛ ⎞ f f − ′ 0 ⎟ e −t Δg2 + 4ab 1−⎜ ⎜ ⎟ ⎝ f0 − f * ⎠

subpopulation ratio. If we start in the steady state ( f0 = f ′, then F = 0), this contribution is zero. An analogous result can be obtained for stress conditions considering h (t ) = p (t ) n (t ) instead of f(t). The average growth rate of the total population during a stress period ts is given by μ Ns

=

h′μ ps

μns

(1 + h′) (1 + h′) 2 ⎡ ⎤ 1 1 − H e −t s Δs + 4ab ⎥ , + ln ⎢ ⎥ t s ⎢⎣ 1−H ⎦

(6)

where

( )

H h0 = (4)

Δs = h′ =

where

+

( 1 + h*) ( h

) , (1 + h′) ( h 0 − h*) ( μ ps − b) − ( μns − a),

(Δ +

− h′

0

Δs2 + 4ab

s

) (2b),

and

Δg = f′ =

(

) (

μng

−a −

(Δ

Δ g2

+

g

μ pg

)

−b ,

+ 4ab

h* =

) (2a)

and f* =

(Δ

g

−

Δ g2 + 4ab

) (2a).

Thus equation (4) describes the relaxation of the subpopulation ratio to its steady state f′. This relaxation is characterized by a relaxation time scale τg = (Δg2 + 4ab)−1 2 ≈ Δg−1, which for small switching rates is approximately equal to the inverse of the growth rate difference (the selection pressure) between the phenotypes. 2.2. Growth rate during growth and stress period Using equation (4) for f(t), the average growth rate of the total population during a growth period tg is found to be μ Ng

=

f ′μng

μ pg

+ (1 + f ′) (1 + f ′) 2 ⎡ ⎤ 1 1 − F e −t g Δg + 4ab ⎥ , + ln ⎢ ⎥ t g ⎢⎣ 1−F ⎦

(1 + f *)(f0 − f ′)

3

s

Δs2 + 4ab

) (2b).

2.3. Growth rate in a periodically switching environment Now we consider an environment that switches periodically between growth (of duration tg) and stress (of duration ts) conditions. We calculate the average growth rate of the population over one environmental cycle (of duration t g + t s ) as μ =

μ Ng t g + μ Ns t s

( t g + ts )

.

Substituting the previous results, we get μ =

(f ′μ

g n

)

+ μpg tg

(

(1 + f ′) ts + tg

+

1 ts +

+

1 ts +

(5)

. The first two terms (1 + f ′)(f0 − f *) are the growth contributions in the steady state from the normal cell fraction (f ′ (1 + f ′)) and persister fraction (1 (1 + f ′)). This result shows that, in the long time limit, the two subpopulation behave as isolated populations growing with a steady state growth rate. The additional term depicts the transient dynamics of the population towards the steady state

where F (f0 ) =

(Δ −

)

+

(h′μ

s p

)

+ μns ts

(

(1 + h′) ts + tg

⎡ −t g Δg2 + 4ab ⎢ 1−Fe ln⎢ tg ⎢ (1 − F) ⎣ ⎡ −ts Δs2 + 4ab ⎢ 1−He ln⎢ tg ⎢ (1 − H ) ⎣

⎤ ⎥ ⎥ ⎥ ⎦ ⎤ ⎥ ⎥ ⎥ ⎦

(

)

(

)

)

(7)

Again, the linear terms in tg and ts are the contributions from the steady fractions of exponentially growing and dying normal and persister cells in each environments, and the logarithmic term is due to the change in subpopulation ratio when the environment switches from growth to stress and stress to growth. The expression of equation (7) depends on the initial subpopulation ratio during growth (via

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

F ≡ F (f0 )) and stress (via H ≡ H (h 0 )). Here, we are most interested in the long-time behavior, where the initial subpopulation ratio of an environmental phase is given by the subpopulation ratio at the end of the h 0 = [f (t g )]−1 and previous phase, i.e.,

f0 = [h (t s )]−1 . These conditions lead to two recursion relations that can be solved to calculate the exact expression for f0 and h0 for given values ts and tg (see appendix A.2). Therefore, equation (7) gives the average growth rate of the total population without any approximation or model simplification (previous theoretical studies calculated 〈μ〉 either for long environmental durations [15, 20] or for equal switching rate between the different phenotypes [24]).

2.4. Periodically varying environment with short stress durations The exact analytical expression for the growth rate, equation (7), can be used to study the population dynamics and evolution of bacterial persistence over all the range of cellular parameters (growth, death and switching rates) and environmental parameters (different growth and stress durations). Here, we are specifically interested in the case of long periods of growth and short or variable periods of stress. We consider the case that the growth period is sufficiently long, such that the population structure reaches its equilibrium, but allow for arbitrary, in particular short stress durations. Then the subpopulation ratio h 0 (=1 f ′) at the beginning of a stress period is given by the steady state of the subpopulation ratio during growth period, which gives H (h 0 ) = constant. By contrast, the subpopulation ratio at the end of a stress period f0 = 1 h (t s ) and hence F (f0 ) ≡ F (t s ) are dependent on the stress duration ts. Substituting these values, we get

(μ μ =

g n

)

− a + μ pg − b +

(

Δ g2 + 4ab t g

2 ts + t g +

(μ

s p

)

)

− b + μns − a +

(

2 ts + t g

Δs2 + 4ab t s

)

⎡ 4ab − Δ g Δs 1 ⎢ ln ⎢ + 2 2 ts + t g ⎣ 2 Δs + 4ab Δ g + 4ab

(

× 1 − e −t s +

Δs2 + 4ab

1 1 + e −t s 2

(

)

Δs2 + 4ab

)⎤⎥⎦.

(8)

Equation (8) gives the growth rate of the population in a periodic environment with short stress durations. For small switching rates, this expression can be further simplified (the approximation scheme is discussed in appendix A) to 4

μ ≈

( μng − a) t g + ( μ ps − b) ts ( ts + t g ) ( ts + t g ) +

1 ts + t g

2 ⎡ 2 2 −t s Δs ⎢ ab Δ g + Δs + Δ g Δs e × ln ⎢ Δ g2 Δs2 ⎢⎣

(

)

⎤ ⎥ ⎥. ⎥⎦

(9)

2.5. Evolutionary simulation To test for the influence of stochastic effects and to observe the predicted hysteresis (see below), we perform hybrid stochastic-deterministic simulation. In these simulations, we focus on the case of short stress duration, where the net growth rate of the population over an environmental cycle is positive. The population grows exponentially to high values in the long time. To restrict the population size while keeping the dynamics unchanged, the population is diluted after each environmental cycle of growth and stress. Such population dilution methods are commonly carried out in experiments for studying competitions and long term evolution [32]. In our simulation, the dilution step rescales the total population size to a lower value N0. During this rescaling, the sizes of different phenotypes and genotypes are chosen stochastically with the constraint that their fractions in the total population (ϕi (t )) remain unchanged on average. This can be understood as retaining each individual cells with probability N0 N , where N is the total population size before dilution. Practically, the stochastic effect of dilution is implemented by choosing the subpopulation size for each phenotype (i) from a Poisson distribution with mean ϕi N0. For comparison with the deterministic theoretical analysis, we compare the average growth rate of the population with stochastic dilution step (figure 3). Next, we introduce a finite mutation rate at which a genotype with a modified phenotype switching rate is generated from the existing cells. In a mixed population with several genotypes, the stochastic dilution acts as a population bottleneck whereby less fit genotypes are removed from the total population. The time required by the best genotype to be fixated in a population depends on its fitness advantage and the amount of dilution. The scheme of our evolutionary simulation is as follows: (i) population growth, death and switching between the phenotypes are taken as deterministic process described by equation (1). This follows from the assumption that we are simulating a population with positive growth rate over an environment cycle. Extinction is incorporated into the deterministic dynamics, by setting population sizes to zero, whenever they reach a value below one. (ii) The environment (growth and stress conditions) is switched in a deterministic fashion with fixed durations, i.e., after times tg and ts for the switch from growth to stress and

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

Figure 1. Fitness landscape for phenotype switching: (a) the average growth rate is plotted as a function of the switching rate (b = a) for different stress durations. The solid lines are plotted using equation (8) and the points are from numerical simulations of equation (1) (population growth rate averaged over 200 cycle of growth and stress phases). The plot shows the existence of a minimum and maximum in the growth rate for intermediate stress durations. The insets shows the growth rate for very short and long stress durations, with a boundary maximum at a = 0 and a single maximum at finite a, respectively. (b) Bifurcation diagram: positions of the growth rate maximum (solid line) and minimum (broken line) as functions of the stress duration ts. The red arrow indicates an example of an evolutionary transition from a finite switching rate to no switching rate as the stress duration is decreased, the blue arrow shows the transition from no phenotype switching to the optimal switching rate as the stress duration is increased. (c) Comparison of the maximal (solid) and minimal (dashed line) growth rate to the growth rate at a = 0 (gray line). All growth rates are normalized to the one at a = 0. The distance between the dashed curve and the gray line corresponds to the depth of the fitness valley that needs to be crossed for a transition from a = 0 to a = a max . Here tg = 10 h, μng = 2 h−1, μ pg = 0.2 h−1, μns = −2 h−1, μ ps = −0.2 h−1.

from stress to growth, respectively. (iii) After each environmental cycle (t g + t s ), the total population is diluted to N0, while keeping the fraction of subpopulation fixed on average, but allowing for fluctuations by using a Poisson distribution (Poisson-distributed random numbers are generated with the Knuth method when the mean of the distribution is < 500; for larger means, we use Gaussian-distributed random numbers as an accurate approximation [33]). (iv) After each dilution step, a cell (say with a genotype giving rise to the switching rate a0) is randomly chosen that can generate a mutant cell with lower (a 0 ϵ−1) or higher (a 0 ϵ ) 1 switching rate with equal probability PM . We note 2 that depending on the numerical value of the mutation rate, multiple mutant genotypes can coexist in our simulations, in particular when the fitness difference between successive genotypes is small. These steps are repeated for sufficiently long time to allow the population to reach an evolutionary stable switching rate. Furthermore, we use changing stress conditions to understand features of evolutionary transitions between non-switching and finite switching population. We fix the duration of growth (tg = 10 h in the simulations presented here) and modulate the stress duration with t s < t g . The stress duration ts is increased or decreased linearly with a rate much slower than the environmental variation rate ((t s + t g )−1). For example, in the results shown in figure 5, the stress duration is increased from a low value to a high value in t = 105 h, whereas t s + t g ⩽ 20 hours.

3. Results One biological role of phenotypic heterogeneity and switching between phenotypes is bet-hedging, an adaptation to a variable environment. In the special 5

case of persistence, one phenotype is adjusted to rapid growth under favorable conditions, but survives poorly under stress, while the other phenotype grows slowly, but has higher tolerance to stress. Several previous theoretical studies have shown the presence of evolutionary optimal switching rates that maximize the long-term growth rate of the total population when the durations of growth and stress periods are large or when the environment fluctuates very slowly. [15, 17, 18, 20, 24, 25]. The underlying idea is that a population can carry information about an environment through the subpopulation allocation that is determined by the current environmental condition and the phenotype switching rates. This subpopulation allocation can then be optimized for maximal growth in varying environmental conditions by adapting the phenotype switching rates to the rate of environmental variation [20]. The evolutionary optimal switching rate is found to be proportional to the frequency of the varying environment [15, 17, 20, 24]. These results can easily be recovered from equation (7) (see appendix A). The optimal switching rates between the phenotypes for a population growing in a slowly varying environment of growth period tg and stress period ts are found to be ⎛ 1 1⎞ ≈ t g − 2 ⎜⎜ − ⎟⎟ , Δs ⎠ a opt ⎝ Δg ⎛ 1 1 1⎞ ≈ t s + 2 ⎜⎜ − ⎟⎟ . Δs ⎠ bopt ⎝ Δg 1

(10)

The existence of optimal switching relies on the interplay between the cost (due to switching to a slowgrowing state during growth) and benefit (due to extended survival through the slow dying persisters). If growth and stress periods are sufficiently long, the benefits of phenotypic heterogeneity outweigh the

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

costs, because the fraction of mal-adapted cells are typically small. However, when the duration of the stress periods are short, this is not necessarily the case and the presence of mal-adapted persisters during the long growth phase may incur a cost that exceeds the benefit due to the persisters during the shorter stress phases. In such a case, a homogeneous population may have a fitness advantage over a heterogeneous one. 3.1. Long-time growth rate in environments with periodic short episodes of stress To get insight into the conditions under which phenotype switching is advantageous, we determined the long-time growth rate of a population in an environment that alternates between long growth periods and short stress periods. In figure 1(a), we plot the growth rate obtained from the analytical expression given by equation (8) together with numerical results for different stress durations. This plot can be interpreted as a fitness landscape, which shows optimal values of the phenotype switching rate (note that we set a = b in this plot) that maximize the growth rate, as well as the effects of modulating this rate on fitness (on evolutionary time scales). In case of long stress duration, the growth rate indeed exhibits a unique maximum (figure 1(a), inset (b)), which indicates that phenotype switching (with a non-zero optimal switching rate) is beneficial for the population growth under these environmental conditions. For very short stress duration, however the growth rate decreases monotonically with increasing switching rate (figure 1(a), inset (a)), which suggests that a homogeneous population without phenotype switching is the fittest strategy in these conditions. In the intermediate range (figure 1(a), square, circle and triangle symbols), the fitness landscape shows two fitness maxima, one at a finite switching rate and one at the limiting case of no switching. The two maxima are separated by a (local) fitness minimum. The global fitness minimum typically occurs in the limit of infinite switching rates, but this case is not relevant for the following analysis, where we are mainly interested in the maxima of the growth rate. Thus, when we refer to the minimum of the growth rate in the following, we mean the local minimum separating the two maxima, which is an important determinant of the dynamics of evolution toward the fitness maxima. 3.2. Analytical expressions for the optimal switching rates maximizing the growth rate To obtain an analytical expressions for the maximum and minimum of the growth rate, we use the approximate expression for the growth rate, equation (9). First, we consider the case when both switching rates are equal (b = a). Differentiating with respect to a gives the maximum (a + = a max ) and the

6

minimum (a− = a min ) as a± ≈

1 ⎡ ⎢1 ± t g + ts ⎣

⎤ 2 1 − t g + t s δ 2 e −t s Δs ⎥ , (11) ⎦

(

)

where δ = (1 Δg + 1 Δs )−1. The maximum and minimum can be interpreted as the fixed points of (d 〈μ〉 da), the derivative of the growth rate. These fixed points are destroyed through a saddle-node bifurcation at a certain threshold in the stress duration. This behavior is shown in figure 1(b), which shows the optimal switching rate as a function of the stress duration ts. The maximum is the stable node (shown by solid line) and the minimum is the unstable node (shown by broken line) in evolutionary sense. The dotted line in figure 1(b) shows the boundary maximum at a = 0 (for which the condition d 〈μ〉 da = 0 is not fulfilled): for a , b → 0, equation (9) reduces to (μng t g + μns t s ) (t s + t g ) − a , which shows that the case a = 0 indeed represents a boundary maximum for any value of the stress duration. Thus, one can see that in an intermediate range there are two locally optimal solutions, i.e. two values of the switching rate that maximize growth compared to small variations of the switching rate. The two maxima are separated by the minimum (amin), and the separation is determined by the stress duration ts. Thus such populations are stable against mutations that make small changes of the switching rate. However, a population with the lower growth rate can be invaded by the one with the higher growth rate. Figure 1(c) shows the growth rates at the maximum and minimum (solid and dashed line, respectively), normalized to the growth rate at the boundary maximum, i.e. for a = 0. One can see that the global maximum of the growth rate is given by the boundary maximum for short stress durations, but beyond a critical stress duration, the growth rate at amax is higher and thus phenotype switching becomes the globally fittest, i.e., evolutionary stable strategy. However, we also observe that the minimal growth rate at amin rate is always smaller than the growth rate for a = 0 and approaches the latter for long stress duration. Therefore, if a population is allowed to evolve from zero switching rate, it has to overcome a fitness valley of depth 〈μ〉 (a = 0) − 〈μ〉 (a min ). The duration of stress that is required for such evolutionary transition will depend on the strength of the mutational noise, as a sufficiently large fluctuation is required to cross this fitness valley (e.g., as indicated by the blue arrow in figure 1(b)). In contrast, if a population has a finite switching rate and adapts to shortening stress durations, the evolutionary reverse transition to non-zero switching will occur when the fitness maximum disappears (indicated by the red arrow in figure 1(b)). Thus, the presence of such a fitness minimum and maximum may lead to an evolutionary hysteresis where different thresholds depending on the history of

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

Figure 2. Maximum and minimum in the average growth rate as functions of the stress duration for different evolutionary scenarios. (a) Different differences in the growth rates or selection pressures, Δs = (μps − μns ), with fixed Δg = (μng − μ pg ) = 1.8 and η = 1. (b) Different values of switching rate ratio η. (c) Maximum of the growth rate for the case where the switching rate b is fixed, while a is varied. (d) Maximum of the growth rate for the case where a is fixed and b is varied. (e) Dependence of the minimal stress duration required for the existence of phenotype switching, t s* , on the ratio of switching rates η and (f) on the selection pressure during stress conditions Δs . In these plots, (μng − μ pg ) = (μ ps − μns) = 1.8, unless stated otherwise.

a population are required for transitions between the non-switching and the switching strategy. 3.3. Discontinuous transition for short stress duration and a fixed ratio of the switching rates Next, we are interested in how the presence and absence of phenotype switching depends on the different parameters of the model. We therefore determine the bifurcation diagram for different scenarios (figures 2(a)–(d)). The evolution of phenotype switching reflects changes in the underlying genetic circuits, for example, in the case of persister cells, the number of toxin–antitoxin systems in the genome, the relative expression level of toxins and the respective antitoxins, or their overall expression level [34, 35]. 7

Such changes may affect only one of the switching rates between the two phenotypes or both, and if both, the effects can be similar or opposing [36–39]. Here, we specifically consider the two cases, where both the switching rates change while their ratio is fixed (e.g., through the increase or decrease of stochasticity in toxin–antitoxin expression) and the case when one switching rate varies while the other is fixed (e.g., via a change of the relative expression of toxin and antitoxin). In the first case, we write b = ηa and take η as fixed, while a is varied. This scenario includes the case of equal switching rates (a = b), which has been studied previously [15, 24]. Then the maximal and minimal growth rates are found for

Phys. Biol. 12 (2015) 046004

1− amin ≈

(

P Patra and S Klumpp

)(

)

1 − Tg + ηTs Tg + Ts δ 2 e −tsΔs η Tg + ηTs

and 1+ amax ≈

(

)(

)

1 − Tg + ηTs Tg + Ts δ 2 e −tsΔs η Tg + ηTs

, (12)

respectively, where Tg = t g − 2 Δg + 2 Δs and Ts = t s + 2 Δg − 2 Δs . The optimal switching rate amax approaches 2 (Tg + ηTs ) for large ts. The saturation value depends on the ratio of the switching rates, η = b a , as shown in figure 2(b). η = 1 reproduces the case of equal switching, where a max ≈ 2 (t s + t g ) i.e. the optimal switching rate is approximately equal to the mean rate of environmental change [24]. The minimal value of ts for which the roots are positive (i.e., the threshold above which the maximum and minimum exist) is given by e t s* Δs =

δ2 Tg + ηTs η

(

)( Tg + Ts ).

(13)

This condition is a nonlinear equation with a ts dependence on the right-hand side, which can be solved iteratively by treating the ts dependence on the right-hand side as a small perturbation. To leading order, the approximate value of t s* for long growth phases (large tg) is obtained as t s* ≈

⎡ δ2 ⎤ 1 ln ⎢ T g2⎥ . Δs ⎣ η ⎦

(14)

This approximation agree well with the exact value for small η, i.e., b ⩽ a as shown in figures 2(e) and (f) (agreement for large η can be improved by taking the approximation to the next order, see the appendix A.1). Therefore, if the environmental stress duration is larger than this threshold value the population will benefit from persister formation and persistence will be sustained. Equation (14) gives the dependence of the threshold value on the selection pressures Δs and Δg in the two environments, the switching rate ratio η and the duration of growth tg. The dominant dependence is on the relaxation time τs = 1 Δs of the subpopulation ratio during the stress period, i.e., the stress duration has to be longer than a time proportional to that relaxation time for phenotype switching to be beneficial and thus evolutionary stable. If Δs varies, the threshold value changes while the maximum remains fixed at a ≈ 10−1 and, when η varies, it changes both the threshold and the maximum of growth rate. In evolutionary terms, transitions from a population with high switching rate to a population without phenotype switching (zero switching rate) will occur at this threshold. This evolutionary transition is similar to a first order phase transition in physical systems, the evolutionary stable switching rate changes discontinuously from a finite value (≈10−1) to zero as the stress duration is decreased [24].

8

3.4. Continuous transition when only one switching rate is adapted Next, we consider the cases, where only one switching rate is varied [15, 20]. In this case, we find a single maximum. Variying a with fixed b, the maximal growth rate is found for −1 ⎛ 2 2⎞ a max ≈ ⎜⎜ t g − + ⎟⎟ Δg Δs ⎠ ⎝ ⎛ ⎞ δ 2 t g + ts e −t s Δs ⎟ . × ⎜1 − ⎜ ⎟ b ⎝ ⎠

(

)

(15)

For large ts, this expression approaches aopt (equation (10)). In order for phenotype switching to be beneficial, this expression must by positive (when it is negative, the mathematical maximum occurs for a negative switching rate, and thus, the maximal growth rate in the physically accessible range of parameters, i.e., positive rates, is the boundary value at a = 0). For long growth durations, the first term is typically positive and the second term provides a constraint on the stress duration, which must be larger than t s ≳ Δs−1 ln[(δ 2 b) t g ] (the cutoff in blue and black curves in figure 2(c)). The evolutionary transition in this scenario is similar to a second order phase transition in physical systems, as the first derivative of the evolutionary stable switching rate is discontinuous at this threshold. For short growth durations, the first term may also become negative, but only if there is an asymmetry in the selection pressures between the growth and stress periods, a point that had been emphasized in previous work [24]. Varying b for fixed a, the maximal growth rate is −1 ⎛ 2 2⎞ − ⎟⎟ . In this case bmax found for bmax ≈ ⎜⎜ t s + Δg Δs ⎠ ⎝ diverges at a certain threshold value in ts, specifically when t s + 2 Δg ≈ 2 Δs (purple and gray curves in figure 2(d)). Thus, a maximal growth rate for an optimal switching rate bmax is seen only when the stress duration is larger than a threshold that depends on the difference in the selection pressures and requires again an asymmetry between the selection pressures during the stress and growth phases [24]. Compared to the variation of a with fixed b, only one of the two mechanisms for the transition is found. The reason for this asymmetry is that we have considered long growth periods, so that a condition corresponding to the requirement that the second term of equation (15) is positive is always satisfied. 3.5. Evolutionary transitions and hysteresis in finite populations The observation of a minimal growth rate in the fitness landscape determined above suggested that one may observe hysteresis in the evolution of phenotype switching. To observe that evolutionary hysteresis, we use a

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

Figure 3. Average growth rate from the semi-stochastic simulation: growth rate averaged over 1000 cycle of growth, stress and dilution. The population dilution occurs in every cycle after growth and stress phase. (a) For short stress phase, the surviving population is large, therefore the dilution does not change the subpopulation ratio. Hence, we get the growth rate as predicted from theoretical analysis. (b) For long stress phase, the surviving population is small, therefore dilution might cause the loss of a particular phenotype (inset: population fractions) leading to an abrupt decrease in fitness. The parameters are the same as in figure 1.

Figure 4. Simulated evolution of phenotype switching: (a) scheme of stochastic transitions showing how the switching rate is modified through mutation. PM 2 is the probability of a mutation increasing (forward) or decreasing (backward) the switching rate. ϵ is a multiplicative factor by which the switching rates are modulated. (b) Trajectory of a simulation: population size as a function of time with several cycles of growth, stress and dilution.

semi-stochastic simulation method which describes cycles of growth and stress by deterministic population dynamics, but as in an experiment, the population is diluted at fixed time intervals and during this dilution steps, stochasticity is introduced by randomly discarding or retaining individual cells. Moreover, mutations are also described as stochastic events in single cells. First we checked that in the absence of mutations, the stochastic simulation showed the same effective fitness landscape as our deterministic theory, at least for large population size. This is the case; in particular the bistable landscape is seen for short stress durations (figure 3(a)). The population bottleneck is the dilution step, which is performed after the stress phase. Therefore, after a long stress phase, the minority subpopulation in the 9

population may get extinct by dilution. Such extinctions are indeed seen in the intermediate range of switching rates for small population size N0 (figure 3(b), case N0 = 100), when the normal cells get extinct during dilution and cause an abrupt decrease in the growth rate. The extinct population can get replenished during growth through phenotype switching from the surviving population and when the switching rates are large the system recovers the deterministic fitness landscape. Next, we include mutations, which are implemented as stochastic transitions of a single cell to a different genotypic state with either lower or higher phenotype switching rate (figure 4(a)). Figure 4(b), shows a short run of such a simulation, where one can observe the growth (tg = 10 hours) and stress periods

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

Figure 5. Hysteresis in evolutionary simulations: the simulation starts with a single population with a given switching rate (10−3 here). The individual cells are allowed to mutate to modify their switching rate (following the rules of the network shown in figure 4(a)) with fixed probability PM 2 = 0.25 and population size N0 = 1000 after dilution. During an evolutionary run, a population with a higher growth rate than the other populations overgrow and replace others, and eventually a population with an optimal switching rate emerges. The upper panel shows the case when ts is increased, starting from a low value, and the lower panel shows the opposite scenario. (a) The switching rates are modulated together, with b = a. The evolutionary transitions occur at different stress durations and evolutionary hysteresis is observed. (b) Only one switching rate is varied (a), while the other is fixed (b = 0.1) hour−1 fixed. In this case, the evolutionary transitions occur at the same values of the stress durations indicating the absence of hysteresis.

(ts = 6 hours) followed by a dilution of the total population. In this example, a mutation occurs after several cycles and introduces a new genotype with a different phenotype switching rate. If this new population is fitter than the existing one (which is the case here), it proliferates to form a majority in the total population and after successive dilution it becomes fixated and the less fit population is lost. Finally, to test for hysteresis, we allow the duration of stress to increase or decrease during the simulation. This represents a case of evolution in a slowly changing fitness landscape, where the population always reaches the nearest fitness peak [40]. The expected hysteresis is indeed observed when the two switching rates are varied together (i.e., for a = b): we observe that an increase in the duration of stress can induce persistence (figure 5(a)), i.e., it causes a transition from a population with very low switching rate (red curve) to a population with high switching rate (black curve). As 10

the duration of stress is decreased, the reverse transition occurs at a much lower threshold in the stress duration as we expected (figure 5(b)). By contrast, if we fix one of the switching rates (say b) and allow only the variation of the other switching rate, i.e., the one way switching rate to persister state (rate a), we do not observe this evolutionary hysteresis, as expected from the bifurcation diagrams in figure 2(d). This observation indicates that the presence of such hysteresis strictly depends on the way the parameters are allowed to vary. In a biological sense, the case when both rates of phenotype switching are affected together corresponds to the case where the overall level of noise in the gene regulatory network that affects the switching is modulated [41–43]. On the other hand, cases where the relative expression levels of network components are varied (e.g. of toxin and antitoxins), correspond to case where only one switching rate varies while the other is fixed [42, 43].

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

4. Discussion Phenotype switching allows bacterial populations to diversify into distinct phenotypes with different growth and survival properties, which serves as a bethedging strategy in varying environments [15]. Several theoretical studies have shown that optimal switching rates between the phenotypes can be found that maximize the long-term growth rate in a stochastically or periodically varying environment [15, 17, 18, 20, 24, 25, 44]. This observation indicates that switching rates are evolvable quantities that can be adapted to the time scales of environmental changes. However, further studies provided evidence that the evolutionary transition from non-switching to switching in a population could be of discontinuous nature [24, 27] and thus that phenotype switching is not always beneficial. Here we have extended this line of analysis by deriving an explicit analytical solution for the population growth rate and analyzing the conditions under which phenotype switching becomes unfavorable, as indicated by the lack of optimal switching rates. We found several conditions that make switching unfavorable. In particular, our analysis shows that a threshold duration of the stress phase is required for the existence of a fitness maximum and hence for a benefit due to persister cells and phenotype switching. This threshold depends on the selection pressure difference between the phenotype in the different environments [24, 27, 45] as well as the ratio of switching rates and the duration of the growth period, with a dominant dependence on the relaxation time of the subpopulation ratio (which can be approximated by the difference in death rates between the two subpopulations). For shorter stress durations, the benefit due to persisters does not outweigh their cost during the longer periods of unstressed growth and a non-switching population is fitter with higher growth rate in the long time average. Surprisingly, the evolutionary transitions between absence and presence of phenotype switching was found to be either continuous or discontinuous, depending on the mode of evolutionary modulation of the switching rates. For back-and-forth switching rates between the two phenotypes evolving together, there is a fitness valley (a minimum of the growth rate) for intermediate stress durations that separates switching and non-switching phenotypes. As a consequence, a population that is not the globally most fit one, may be stable against small changes in its switching rate due to mutations, but vulnerable against invasion of a genotype with very different switching rates. Under slow modulation of the stress duration, this situation results in evolutionary hysteresis. By contrast, if the switching rates between the phenotypes vary individually, the transition is continuous without a fitness valley and without metastable populations and no hysteresis is obtained. We have checked the 11

occurrence of hysteresis with an stochastic simulation of evolving phenotype switching rates. In principle, experimental tests of evolutionary hysteresis could be used to distinguish continuous and discontinuous transitions. However, such a test requires that changes in the switching rates due to individual mutations are typically small. The mode of evolutionary modulation of the switching rates can be linked to the molecular networks underlying phenotype switching. Persistence in E. coli depends on multiple toxin–antitoxin systems. The levels of the toxin and antitoxin can modulate the switching rates in various fashions [46]. Deletion of an increasing number of toxin–antitoxin systems progressively decrease the level of persisters [34], likely by reducing the overall stochasticity in the system and thus the switching rates. On the other hand, for the example of the hipBA toxin–antitoxin system, the toxin (HipA) level determines the duration of dormancy (i.e., the rate b) and the switch to persistence (via the rate a) occurs when the toxin exceed a threshold value that is determined by the antitoxin HipB level [35, 46]. Thus in this case, the relative expression of toxin and antitoxin modulates the individual switching rates. A systematic modulation of phenotype switching rates has been experimentally achieved in a synthetic fim switch in E. coli [39]. In this particular system, the expression (ON/OFF) of the fim switch is controlled by two recombinases, FimB and FimE. FimB modulates the transition rate in both directions (ON to OFF and vice-versa), whereas FimE only modulates the transition from ON to OFF state. This example directly corresponds to the two discussed scenarios in the previous section and the simulation results shown in figure 5. Further, circuit-based evolutionary models coupling single cell dynamics to population growth and environmental selection are required [47] to capture the different modes of evolutionary transitions discussed here. Finally, we note that in our analysis, we have assumed that mutation or modification occurs only in the rates of phenotype switching during evolutionary adaptation. However, evolving population will also adapt by modifying their growth rate [32] or by resistance mutation [48]. However, a recent study noted that the first adaptive response to antibiotic stress occurs through modification of phenotype switching rates (i.e., the lag-time) rather than the growth rates or mutations towards resistance [49].

Appendix A A.1. Growth rate in a slowly changing environment In the well studied scenario of a slowly varying environment (such that the steady state in the subpopulation ratio is reached in every environmental phase), f0 = 1 h′ and h 0 = 1 f ′. Thus, the average growth rate is given by

Phys. Biol. 12 (2015) 046004

( f ′μ

g n

μ =

P Patra and S Klumpp

)

+ μ pg t g

(

(1 + f ′) t s + t g

+

)

( h′μ

s p

)

+ μns t s

(

(1 + h′) t s + t g

)

) ⎤⎥⎥. ) ⎥⎦

)(

(

μ =

((

) (

−a +

μ pg

)

−b +

(

2 ts + t g +

((

) (

)

μ ps − b + μns − a +

(

2 ts + t g +

Δ g2

)

)

(

(1 + g ′) t s + t g

)

( ))

(A.4)

( )

F ts =

)

f0 =

)

0

and

⎤ ⎥ ⎥ .(A.2) ⎦

Δs2 + 4ab

1 − Be −t s

Δs2 + 4ab

g ′ − g *Be −t s

.

(A.5)

Substituting these values and using the approximation scheme, we arrive at equation (6). For b = ηa , the growth rate is μ

μ ps − ηa) t s μng − a) t g ( ( ≈ + ( ts + t g ) ( ts + t g ) 2 ⎡ 2 2 2 −t Δ ⎤ 1 ⎢ ηa ( Δ g + Δs ) + Δ g Δs e ⎥ ln . + s

ts + t g

⎢ ⎢⎣

s

⎥ ⎥⎦

Δ g2 Δs2

Differentiating with respect to a, we get

(

Tg + ηTs d μ ≈− da t g + ts +

)

1 ⎛ 2ηa + δ 2 (η − 1) Ts e −t s Δs ⎞ ⎜ ⎟ , (A.6) t g + ts ⎝ ηa 2 + δ 2 e −t s Δs ⎠

where δ = Δg Δs (Δg + Δs )−1, Tg = (t g − 2/Δg + 2 Δs ) and Ts = (t s + 2 Δg − 2 Δg ). Solving for d 〈μ〉 da = 0, we get the maximum (+) and minimum (−) at a± 1± ≈

( μng − a) t g + ( − b) ts ≈ ( ts + t g ) ( ts + t g ) 2⎤ ⎡ 1 ⎢ ab ( Δ g + Δs ) ⎥ ln . +

( 1 + f *) ( f − f ′) (1 + f ′) ( f − f * ) 0

μ ps

Δ g2 Δs2

)

+ μns t s

)

(

Δs2 + 4ab t s

Equation (A.2) exactly predicts the numerically computed growth rate for all values of the switching rates and sufficiently long environmental durations. Compared to previous theoretical results [15, 25] that were based on additional approximations, it agrees better with the numerical results, in particular when the phenotype switching rates are large. Approximation method. To reduce the above expression for the growth rate for small switching rates, we follow a simple approximation scheme, expanding to zeroth or first order in ab for ab ≈ 0 . Thus, we assume Δi2 + 4ab ≈ Δi (in all term containing growth or death rates, as theswitching rates are small compared to them, e.g., for the linear terms in ts and tg) or Δi2 + 4ab ≈ Δi + 2ab Δi (in terms containing contributions proportional to the switching rates, e.g., for the numerator term inside the logarithm). Using this approximation, the expression for the growth rate is reduced to

⎢ ⎢⎣

)

s p

(

+ 4ab t g

1 ts + t g

ts + t g

(

( g ′μ

+

where

⎡ 2 2 ⎢ 4ab − Δ g Δs + Δ g + 4ab Δs + 4ab × ln ⎢ 2 Δs2 + 4ab Δ g2 + 4ab ⎣

μ

+ μ pg t g

⎡ ⎤ −t Δ 2 + 4ab ⎢ 1−He s s ⎥ 1 ⎥, ln ⎢ + ts + t g ⎢ 1 − F t s (1 − H ) ⎥ ⎣ ⎦

(A.1)

Substituting the results for f′ and h′, we get μng

)

g n

(1 + f ′) t s + t g

⎡ ⎢ f ′ − f * h′ − h* 1 ln ⎢ − ts + t g ⎢ 2 − f ′h* + h′f * ⎣

(

( f ′μ

μ =

(

1 − Tg + ηTs

)( Tg + Ts ) δ 2 e−t Δ s

s

η . (A.7)

Tg + ηTs

The minimal value of ts for which the roots are positive, t s* , is given by (A.3)

⎥ ⎥⎦

The optimal switching rates presented in equation (7) are obtained by maximizing this expression. A.2. Growth rate for short stress periods For short stress periods, but long growth periods, we know the subpopulation ratio h 0 (=1 f ′) at the beginning of a stress period. This gives H (h 0 ) and the subpopulation ratio at the end of a stress period, f0 = 1 h (t s ). Then, the average growth rate is

12

e t s* Δs >

δ2 Tg + ηTs η

(

)( Tg + Ts ).

(A.8)

Solving this equation iteratively, and going to the next higher order (first order) compared to the expression given in the main text (t s*(Ts → 0) ≈ (1 Δs )ln[δ 2 Tg2 η]) results in t s* (1) ≳

⎡ δ2 1 ln ⎢ Tg + ηΔs−1 ln ⎡⎣ δ 2 T g2 η⎤⎦ Δs ⎣ η ⎤ × Tg + Δs−1 ln ⎡⎣ δ 2 T g2 η⎤⎦ ⎦ ,

(

(

)

) (A.9)

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

which improves the agreement with the numerical result compared to the zeroth-order result. A.3. Initial subpopulation ratio in changing environments In a periodic environment, the ratio n p = f (t g ) at the end of a growth period tg determines the inverse ratio h 0 = p n at the beginning of a stress period and vice versa. This gives the following two conditions

()

f tg =

=

h(ts) =

=

⎛ f − f′⎞ 2 ⎟⎟ e −t g Δg + 4ab f ′ − f * ⎜⎜ 0 ⎝ f0 − f * ⎠ ⎛ f − f′⎞ 2 ⎟⎟ e −t g Δg + 4ab 1 − ⎜⎜ 0 ⎝ f0 − f * ⎠ 1 h0 ⎛ h − h′ ⎞ −t Δ 2 + 4ab h′ − h*⎜ 0 ⎟e s s ⎝ h 0 − h* ⎠ ⎛ h − h′ ⎞ −t Δ 2 + 4ab 1−⎜ 0 ⎟e s s ⎝ h 0 − h* ⎠ 1 f0

Figure A1. Steady state subpopulation ratios f0 and h0 : comparison between numerical(symbols) and analytical (lines) results.

References

(A.10)

Substitution of the value of f0 or h0 in one of the above equations leads to a quadratic equation in f0 or h0, which has the following solution f0 =

B±

B 2 − 4AC 2A

(

B = f * − f ′e −t g

(

− f ′ − f * e −t g

(

A = 1 − e −t g

(

− 1 − e −ts

(

(

− 1 − e −t g

−ts Δs2 + 4ab

−ts Δs2 + 4ab

Δg2 + 4ab

Δs2 + 4ab

−ts

)

) ) ) ) ).

Δg2 + 4ab

Δs2 + 4ab

C = 1 − e −ts

)(h′ − h*e )(h* − h′e )(h′ − h*e )h′h*(f ′ − f *e )(f * − f ′e )f ′f *(h* − h′e Δg2 + 4ab

−t g Δg2 + 4ab

Δs2 + 4ab

Δg2 + 4ab

−t g Δg2 + 4ab

−ts Δs2 + 4ab

Now, C ⩽ 0, as f* and h* are negative, thus the physical solution is f0 =

B+

B 2 − 4AC , 2A

⎛ a ⎞ − B + B 2 − 4AC h0 = ⎜ ⎟ . ⎝b⎠ 2A

(A.11)

This result gives the exact expression for f0 and h0 for any values of t s , t g and for all parameters. We have checked this solution against the numerically calculated ratio of the subpopulations as shown in figure A1. 13

[1] Veening J-W, Smits W K and Kuipers O P 2008 Bistability, epigenetics, and bet-hedging in bacteria Annu. Rev. Microbiology 62 193–210 [2] Smits W K, Kuipers O P and Veening J-W 2006 Phenotypic variation in bacteria: the role of feedback regulation Nat. Rev. Microbiology 4 259–71 [3] Dhar N and McKinney J D 2007 Microbial phenotypic heterogeneity and antibiotic tolerance Curr. Opin. Microbiology 10 30–8 [4] Di Carlo D and Lee L P 2006 Dynamic single-cell analysis for quantitative biology Anal. Chem. 78 7918–25 [5] Levy S F, Ziv N and Siegal M L 2012 Bet hedging in yeast by heterogeneous, age-correlated expression of a stress protectant PLoS Biol. 10 e1001325 [6] van Heerden J H, Wortel M T, Bruggeman F J, Heijnen J J, Bollen Y J M, Planqué R, Hulshof J, OToole T G, Wahl S A and Teusink B 2014 Lost in transition: start-up of glycolysis yields subpopulations of nongrowing cells Science 343 1245114 [7] Kotte O, Volkmer B, Radzikowski J L and Heinemann M 2014 Phenotypic bistability in Escherichia coli’s central carbon metabolism Mol. Syst. Biol. 10 736 [8] Ozbudak E M, Thattai M, Lim H N, Shraiman B I and van Oudenaarden A 2004 Multistability in the lactose utilization network of Escherichia coli Nature 427 737–40 [9] Solopova A, van Gestel J, Weissing F J, Bachmann H, Teusink B, Kok J and Kuipers O P 2014 Bet-hedging during bacterial diauxic shift Proc. Natl Acad. Sci. 111 7427–32 [10] Wylie C S, Trout A D, Kessler D A and Levine H 2010 Optimal strategy for competence differentiation in bacteria PLoS Genetics 6 e1001108 [11] Veening J-W, Hamoen L W and Kuipers O P 2005 Phosphatases modulate the bistable sporulation gene expression pattern in bacillus subtilis Mol. Microbiology 56 1481–94 [12] Arnoldini M, Vizcarra I A, Miller R P, Stocker N, Diard M, Vogel V, Beardmore R E, Hardt W-D and Ackermann M 2014 Bistable expression of virulence genes in salmonella leads to the formation of an antibiotic-tolerant subpopulation PLoS Biol. 12 1001928 [13] Lewis K 2007 Persister cells, dormancy and infectious disease Nat. Rev. Microbiology 5 48–56 [14] Acar M, Mettetal J T and van Oudenaarden A 2008 Stochastic switching as a survival strategy in fluctuating environments Nat. Genetics 40 471–5 [15] Kussell E, Kishony R, Balaban N Q and Leibler S 2005 Bacterial persistence: a model of survival in changing environments Genetics 169 1807–14 [16] Keren I, Kaldalu N, Spoering A, Wang Y and Lewis K 2004 Persister cells and tolerance to antimicrobials FEMS Microbiology Lett. 230 13–8

Phys. Biol. 12 (2015) 046004

P Patra and S Klumpp

[17] Lachmann M and Jablonka E 1996 The inheritance of phenotypes: an adaptation to fluctuating environments J. Theor. Biol. 181 1–9 [18] Thattai M and van Oudenaarden A 2004 Stochastic gene expression in fluctuating environments Genetics 167 523–30 [19] Balaban N Q, Merrin J, Chait R, Kowalik L and Leibler S 2004 Bacterial persistence as a phenotypic switch Science 305 1622–5 [20] Kussell E and Leibler S 2005 Phenotypic diversity, population growth, and information in fluctuating environments Science 309 2075–8 [21] Beaumont H J E, Gallie J, Kost C, Ferguson G C and Rainey P B 2009 Experimental evolution of bet hedging Nature 462 90–93 [22] Bigger J W 1944 Treatment of staphylococcal infections with penicillin by intermittent sterilisation Lancet 244 497–500 [23] Stewart G R, Robertson B D and Young D B 2003 Tuberculosis: a problem with persistence Nat. Rev. Microbiology 1 97–105 [24] Gaál B, Pitchford J W and Wood A J 2010 Exact results for the evolution of stochastic switching in variable asymmetric environments Genetics 184 1113–9 [25] Patra P and Klumpp S 2013 Population dynamics of bacterial persistence PLoS One 8 e62814 [26] Patra P and Klumpp S 2014 Phenotypically heterogeneous populations in spatially heterogeneous environments Phys. Rev. E 89 030702 (R) [27] Salathé M, Van Cleve J and Feldman M W 2009 Evolution of stochastic switching rates in asymmetric fitness landscapes Genetics 182 1159–64 [28] Kussell E and Vucelja M 2014 Non-equilibrium physics and evolutionadaptation, extinction, and ecology: a key issues review Rep. Prog. Phys. 77 102602 [29] Moyed H S and Bertrand K P 1983 hipA, a newly recognized gene of Escherichia coli K-12 that affects frequency of persistence after inhibition of murein synthesis J. Bacteriol. 155 768–75 [30] Moyed H S and Broderick S H 1986 Molecular cloning and expression of hipA, a gene of Escherichia coli K-12 that affects frequency of persistence after inhibition of murein synthesis J. Bacteriol. 166 399–403 [31] Mulcahy L R, Burns J L, Lory S and Lewis K 2010 Emergence of Pseudomonas aeruginosa strains producing high levels of persister cells in patients with cystic fibrosis J. Bacteriol. 192 6191–9 [32] Wiser M J, Ribeck N and Lenski R E 2013 Long-term dynamics of adaptation in asexual populations Science 342 1364–7

14

[33] Knuth D E 1969 Seminumerical algorithms The Art of Computer Programming vol 2 (Reading, MA: Addison-Wesley) [34] Gerdes K and Maisonneuve E 2012 Bacterial persistence and toxin–antitoxin loci Annu. Rev. Microbiology 66 103–23 [35] Rotem E, Loinger A, Ronin I, Levin-Reisman I, Gabay C, Shoresh N, Biham O and Balaban N Q 2010 Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence Proc. Natl Acad. Sci. USA 107 12541–6 [36] McAdams H H and Arkin A 1997 Stochastic mechanisms in gene expression Proc. Natl Acad. Sci. 94 814–9 [37] Acar M, Becskei A and van Oudenaarden A 2005 Enhancement of cellular memory by reducing stochastic transitions Nature 435 228–32 [38] Fraser D and Kærn M 2009 A chance at survival: gene expression noise and phenotypic diversification strategies Mol. Microbiology 71 1333–40 [39] Hung M, Chang E, Hussein R, Frazier K, Shin J-E, Sagawa S and Lim H N 2014 Modulating the frequency and bias of stochastic switching to control phenotypic variation Nat. Commun. 5 4574 [40] Mustonen V and Lässig M 2007 Adaptations to fluctuating selection in Drosophila Proc. Natl Acad. Sci. 104 2277–82 [41] Assaf M, Roberts E, Luthey-Schulten Z and Goldenfeld N 2013 Extrinsic noise driven phenotype switching in a self-regulating gene Phys. Rev. Lett. 111 058102 [42] Choi P J, Cai L, Frieda K and Xie X S 2008 A stochastic singlemolecule event triggers phenotype switching of a bacterial cell Science 322 442–6 [43] Raj A and van Oudenaarden A 2008 Nature, nurture, or chance: stochastic gene expression and its consequences Cell 135 216–26 [44] Jablonka E, Oborny B, Molnár I, Kisdi E, Hofbauer J and Czárán T 1995 The adaptive advantage of phenotypic memory in changing environments Phil. Trans. R. Soc. B 350 133–41 [45] Kussell E, Leibler S and Grosberg A 2006 Polymer-population mapping and localization in the space of phenotypes Phys. Rev. Lett. 97 068101 [46] Maisonneuve E, Shakespeare L J, Jørgensen M G and Gerdes K 2011 Bacterial persistence by RNA endonucleases Proc. Natl Acad. Sci. 108 13206–11 [47] Kuwahara H and Soyer O S 2012 Bistability in feedback circuits as a byproduct of evolution of evolvability Mol. Syst. Biol. 8 564 [48] Davies J and Davies D 2010 Origins and evolution of antibiotic resistance Microbiology Mol. Biol. Rev. 74 417–33 [49] Fridman O, Goldberg A, Ronin I, Shoresh N and Balaban N Q 2014 Optimization of lag time underlies antibiotic tolerance in evolved bacterial populations Nature 513 418–21