## Abstract

The ability of bacteria to spontaneously switch their expressed phenotype from an identical underlying genotype is now widely acknowledged. Mechanisms behind these switches have been shown to be evolvable. Important questions thus arise: In a fluctuating environment, under what conditions can stochastic switching evolve and how is the evolutionarily optimal switching rate related to the environmental changes? Here we derive exact analytical results for the long-term exponential population growth rate in a two-state periodically changing environment, where the environmental states vary in both their duration and in their impact on the fitness of each phenotype. Using methods from statistical physics we derive conditions under which nonswitching is evolutionarily optimal, and we furthermore demonstrate that the transition between the nonswitching and switching regimes is discontinuous (a first-order phase transition). Our general analytical method allows the evolutionary effects of asymmetries in selection pressures and environmental growth rates to be quantified. The evolutionary implications of our findings are discussed in relation to their to real-world applications in the light of recent experimental evidence.

GENETICALLY homogenous populations under uniform environmental conditions have been observed to exhibit phenotypic heterogeneity due to the inherent stochastic nature of cellular processes, including gene expression (Elowitz *et al.* 2002; Maheshri and O'Shea 2007). While many regulatory circuits and their design principles are aimed both at minimizing noise (Faser *et al.* 2004; Lehner 2008) and making cellular functions robust to the inherent noise (Becskei and Serrano 2000; Maheshri and O'Shea 2007), it has also come to our attention that cellular noise can be exploited to generate phenotypically distinct subpopulations in isogenic populations (McAdams and Arkin 1997; Rao *et al.* 2002) as a possible survival strategy in fluctuating environments (van der Woude and Baumler 2004). Of particular interest have been stochastic epigenetic switches that exploit molecular noise to stochastically switch between states of gene expression (Krishna *et al.* 2005; Smits *et al.* 2006; Choi *et al.* 2008; Leisner *et al.* 2008). Feedback-based multistability (Smits *et al.* 2006) and methylation-based switches (Lim and van Oudenaarden 2007; van der Woude and Henderson 2008) have also been shown to generate phenotypic diversity.

Maintenance of phenotypic diversity has mainly been considered a bet-hedging strategy that can maximize the fitness of the genotype in uncertain environments (Stumpf *et al.* 2002; Evans and Dennehy 2005; Beaumont *et al*. 2009), but alternative uses have also been proposed (Ackermann *et al.* 2008). Maintaining phenotypic diversity at all times means that some fraction of the population is always maladapted, but when the environment changes, the previously maladapted subpopulation becomes the fitter. Theoretical studies have demonstrated that if the different environments select for different phenotypes, stochastic interphenotype switching as a bet-hedging strategy can be advantageous: if the environmental change cannot be sensed at all or there is a delay in signal transduction (Wolf *et al.* 2005), if the rate of response is too low (Thattai and van Oudenaarden, 2004), or if the cost of sensing the state of the environment is higher than the cost of maintaining a maladapted subpopulation (Kussell and Leibler 2005).

Gene expression noise has been shown to be genetically tunable at both the single-gene and network levels (Ozbudak *et al.* 2002; Landry *et al.* 2007). This suggests that stochastic interphenotype switching rate is an evolvable trait, as directly confirmed in experiments combining repeated strong selection events with sequencing analysis (Beaumont *et al.* 2009). It has been experimentally demonstrated that different rates of stochastic switching between phenotypes confer different long-term fitness in any particular fluctuating environment (Acar *et al.* 2008). It is also established that when the environmental durations and growth rate differences between environments are large enough, so that population structure closely approaches steady state before the next environmental change, in two-state systems the optimal rate of switching is approximately equal to the rate of environmental change (Lachmann and Jablonka 1996; Thattai and van Oudenaarden 2004; Kussell *et al.* 2005). In general it has been established to be proportional to the rate of environmental change (Kussell and Leibler 2005). Simulations by Salathé *et al.* (2009) extend these results by allowing the selection pressures in the different environments to differ and to take values that do not necessarily imply that steady state is reached. A threshold phenomenon is observed: below a certain sharp threshold selection pressure difference between the environments, the evolutionarily stable interphenotype switching rates approximate the rate of environmental change; above the threshold value, phenotypic homogeneity is favorable.

Motivated by the findings of Salathé *et al.* (2009) that the optimal strategy appears to change discontinuously with selection pressure strength and asymmetry, in this article we use a quantitative analytical framework to ask more generally how the factors determining the selection pressures in the different environments, and the environmental durations, affect whether stochastic switching or phenotypic homogeneity is the evolutionarily stable strategy. In our study we model exponential growth in an infinite population. Therefore there is no feedback between the environment and the populations or between competing genotypes and it is the evolutionarily stable switching rate that maximizes the long-term exponential growth rate of a population.

In this work we derive exactly the expression for the maximal eigenvalue of the time evolution matrix, from which we can compute the long-term growth rate for given parameter values. We introduce our methodology for this procedure in the next section. We then demonstrate the implications of our findings with reference to exact numerical computations. We give our results primarily for the case where there is asymmetry between the fitness of the two phenotypes but the environmental durations are identical, thereby providing a rigorous theoretical basis that confirms previous observations from simulations, and which demonstrates that any change from nonswitching to switching must be discontinuous (a first-order phase transition). We go on to consider the relationship between the environmental durations and the strength of selection pressure in determining the evolutionarily optimal switching strategy. We then extend our analysis to show how a more complex phase diagram emerges when the environmental durations are also asymmetric. Finally we discuss the evolutionary implications of our findings, the possibility of hysteresis effects, and how growth-limiting conditions will perturb our exact results. Our analytical approach therefore confirms existing results, allows these to be extended to more realistic environmental scenarios, and provides deeper insight by capturing the discontinuous nature of the predicted evolutionary outcomes.

## ANALYSIS

We consider an asexually reproducing population in an environment that periodically fluctuates between two states (Ishii *et al.* 1989). Each individual in the population has the capacity to exhibit either one of two phenotypes, A or B, each of which confers higher exponential growth rate in environments *E*_{A} and *E*_{B}, respectively. The growth rate of phenotype *i* in environment *E* is γ. In environment *E*_{A} phenotype A grows faster so that γ > γ, and similarly γ < γ so that phenotype A grows relatively slowly in environment *E*_{B}. The interphenotype switching rates are constant and independent of the state of the environment. The rate of switching from phenotype A to B is denoted *k* and is assumed to equal the rate of switching from phenotype B to A. We define fitness of a genotype in terms of its evolutionary dynamics. The selection pressure in each environment is defined as the growth rate difference between the two phenotypes, and selection pressures may therefore differ between the two environmental states. We take *E*_{A} to be the environment in which the selection pressure (Δ_{A}, defined by Δ_{A} = γ − γ) is smaller. We can then write Δ_{B} = γ − γ as Δ_{B} = Δ_{A} + δ, where δ ≥ 0 is the selection pressure difference between the environments (we henceforth suppress the subscript on the Δ). Optimal switching rate is defined as the evolutionarily stable strategy (Maynard-Smith 1984), where that strategy, over a long epoch of periodic environmental change, is resistant to invasion by differing strategies.

The time evolution of the two phenotypes is given by the following differential equations, where *n _{i}* denotes the number of individuals exhibiting phenotype

*i*, where

*i*can be A or B. The growth rate of phenotype

*i*in environment

*E*is γ, where

*E*can be

*E*

_{A}or

*E*

_{B}:(1)(2)

The population vectorwhose entries are the number of individuals in each phenotype, changes at the rate given by Equation 3, where *M*_{E} is the projection matrix describing the population behavior in environment *E*:(3)

The projection matrices for environments A and B arerespectively. We consider here periodically changing environments, with environment *E*_{A} and *E*_{B} having durations *T*_{A} and *T*_{B}, respectively. We set , where *t* is the total time elapsed and gives the length of one cycle (the sum length of two consecutive environmental states). Then *q* is the number of environmental cycles elapsed that together gives an epoch of changes. The solution of Equation 3 in periodically changing environments is given by Equation 4, where the repetition of the pattern depends on the time elapsed, *t*):(4)

The matrix is the projection matrix that gives the projection of the population through exactly one environmental cycle. If *q* is large, the growth rate will be dominated by the maximal eigenvalue of this matrix. The population grows exponentially as for large *q*, enabling us to define the long-term exponential growth rate of the population, *r*, as(5)

We simplify our approach at this point by assuming that the fitness of the fitness phenotype in each environment is the same, . By exploiting the underlying symmetries of the projection matrix, in a similar manner to projection operators found in statistical mechanics (Ma 1985), we are able to derive the following expression for the eigenvalues of the projection matrix:(6)where(7)and(8)

These analytical results are numerically confirmed by direct computation using the symbolic mathematics package Mathematica.

## RESULTS

Equation 6 is a closed-form expression for the maximal eigenvalue that dictates the population's long-term exponential growth rate. Analysis of this expression therefore reveals (and quantifies) the conditions and mechanisms under which stochastic switching can evolve: treating *r* as a function of *k*, if *r* is globally maximized at *k* = 0, then we predict that stochastic switching will not evolve; if *r* is globally maximized at *k*^{opt} > 0, then stochastic switching at rate *k*^{opt} is evolutionarily optimal; and finally if there are local maxima in λ^{+} at *k* = 0, and at some other *k*^{opt} > 0, then we predict a sharp, discontinuous evolutionary transition from nonswitching to switching in experiments.

We find that in symmetric periodic environments (*T*_{A} = *T*_{B} = *T*) if selection pressure difference, δ, is small, the optimal switching rate is on the order of, but larger than, the rate of environmental change. More specifically, in symmetrical environments (δ = 0) maximizing Equation 6 with respect to *k* in the limit of large *T* reveals that(9)

We can therefore confirm that in periodic environments, in the special case where the selection pressures in the two environments are equal (δ = 0), the optimal switching rate approximates the rate of environmental change (Lachmann and Jablonka 1996; Thattai and van Oudenaarden 2004). Furthermore the next to leading order correction term decreases as selection pressure increases, implying that when selection pressure is high the optimal switching rate will closely match the rate of environmental change. In the case where the environmental durations are different, a similar maximization of (6) reveals that(10)where both *T*_{A} and *T*_{B} are assumed to be large. Therefore the optimal switching rate is approximately the mean rate of environmental change. We have not been able to compute a reliable expression for the next order term, as the presence of the transition between switching and nonswitching complicates the large *T*_{A}, *T*_{B} limit.

Discrete computer simulations by Salathé *et al.* (2009) demonstrated that there is a discontinuous transition, at some threshold selection pressure difference δ between the two environments, from the optimal strategy being nonswitching to being stochastic interphenotype switching at a rate that approximates the value of the rate of environmental change. By studying the invasion dynamics the authors showed the existence of locally and globally evolutionarily stable switching rates at *k* = 0 and *k* = *T*^{−1}, respectively, for small selection pressure difference and at *k* = *T*^{−1} and *k* = 0, respectively, for large selection pressure difference. Our analytical results allow us to plot the fitness as measured by the long-term exponential population growth rate as a function of possible switching rates (Equation 6), and the resulting figure demonstrates why the transition between the different optimal strategies is sharp (Figure 1). When the cost of being maladapted is equal in the two environments, there is a single peak of fitness near the value of the rate of environmental change (*k*^{opt} ≈ *T*^{−1}). For small selection pressures, Δ, if the selection pressure difference, δ, exceeds some threshold value, then the boundary value of the eigenvalue, at *k* = 0, is larger in magnitude than the maximum at *k* ≈ *T*^{−1} and the globally evolutionarily optimal switching rate is *k* = 0. As the selection pressure difference decreases below the threshold, however, the value of the eigenvalue at the local maximum increases until it becomes the global fitness peak. The crucial observation is that even when the *k* = 0 boundary value is the greater, a local fitness peak is still present close to the environmental change frequency. This means that switching cannot continuously evolve from zero to an optimal nonzero value. Transitions of this type are discontinuous in nature, with a sudden jump between different optimal values. They are also called first-order transitions in the language of phase transitions, and such predicted sharp transitions are potentially supported by experimental evidence (see below).

To establish in what circumstances it is preferable to adopt a stochastic switching strategy to maximize fitness in changing environments, we created a phase diagram for the optimal strategy as a function of the selection pressure in one of the environments, Δ, and the selection pressure difference between the environments, δ. We can determine the exact threshold value of the evolutionary switch point by comparing the maximal eigenvalues computed from Equation 6 when *k* = 0 and *k* = *k*^{opt}. This calculation describes a nontrivial curve δ*(Δ) in the (Δ,δ) phase diagram which separates the region where switching is preferred to that where it is not. This curve defines a line of discontinuous (or first-order) transitions between the two possible outcomes and is shown in Figure 2. (Note that the local fitness peak at *k* = 0 is given by *r* = γ − Δ/2, so at the approximate optimal value of *k*^{opt} ≈ *T*^{−1} the equation *r*(1/*T*) = γ − Δ/2 provides an approximation for the curves plotted in Figure 2.) In agreement with the observations of Salathé *et al.* (2009) we demonstrate that for fixed environmental durations, the threshold δ above which a nonswitching strategy is favored increases with the smaller selection pressure. An asymptote is observed: for any environmental duration there is a threshold selection pressure above which switching is always favored for arbitrarily large values of δ. We can compute an estimate for this by expanding Equation 6 in this limit. The maximal eigenvalue at the putative approximate optimal switching rate, *k* = *T*^{−1}, becomes independent of δ leading to the asymptote. Switching occurs whenwhich is satisfied when Δ > Δ*, where Δ* ≈ 5.1 / *T* after numerical investigation.

We found that the range of selection pressure and selection pressure differences for which switching is preferred depends on the environmental duration *T*. For larger environmental durations switching is preferred for even smaller values of selection pressure and higher selection pressure difference compared to smaller environmental periods. In the completely symmetric case (δ = 0 and *T* _{A} = *T* _{B}) while switching is always a better strategy than homogeneity, the long-term fitness of a stochastic switcher genotype relative to a homogenous population increases with increasing environmental period *T* as well as with selection pressure Δ. This is because a homogenous population is maladapted half of the time regardless of the value of *T*, but if *T* is larger then the population structure changes less frequently, a process during which the stochastic switching population is growing at much lower than the long-term average rate. Longer environmental durations also mean that because *k*^{opt} = *T*^{−1} is small, typically only a small proportion of the population will be in the unfit state. If the selection pressure symmetry assumption is relaxed as in Figure 2, then for some range of parameters homogeneity is the preferred strategy, but the long-term fitness of stochastic switchers nevertheless increases with increasing selection pressure in both environments relative to that of homogenous populations. Our analytical results in Figure 2 are reflected in the observation of Salathé *et al.* (2009) that switching evolves if the selection pressures in both environments are large.

We next looked at the effect of changing the environmental durations separately on the threshold values of δ. It is possible to also consider the case where the fitness of the fittest phenotype in each environment changes separately by performing an appropriate rescaling of the derivative relative differences Δ and δ, but this also involves a rescaling of *k* in one environment, which complicates the algebra but not the underlying structure of the problem. We find that if the relative duration of the environment in which the selection pressure is higher is larger (*T*_{B} > *T*_{A}), then the parameter range becomes shorter for which switching is the better strategy and the phase diagram remains qualitatively the same. The phase diagram in the (Δ,δ) plane is shown in Figure 3 for the case where *T*_{A} > *T*_{B}. The diagram shows multiple evolutionary possibilities for switching, with sharp transitions between them. Above a particular value of δ, two possibilities may happen according to the value of Δ. Either the population remains in phenotype B without switching, or it adopts switching between the two phenotypes. Below this value of δ, three possibilities emerge: remain in A, remain in B, or switch. Interestingly, as Δ increases from zero the sequence of possibilities is: remain in B, switch, remain in A, switch. This leads us to conclude that in the symmetric fitness difference case switching is preferred only if the selection pressure is large enough. We can mathematically confirm this by comparing the eigenvalues calculated from Equation 6. We find that in the limit of large *T*_{A}, *T*_{B,} *T*_{A} > *T*_{B,} the relevant comparison for switching at leading order is(11)

This allows us to establish, by graphical comparison, that switching is always the preferred strategy when *T*_{A} = *T*_{B} (as anticipated) and that deviation from this leads to the emergence of a transitional region where switching is not favored. Switching is subsequently recovered in the limit of large *T*_{A}Δ_{A} and *T*_{B}Δ_{B}. This analytical result is consistent with the findings of Kussell and Leibler (2005) that the optimal switching rates are proportional to the frequency of environmental change in the limit of large *T*_{A}Δ_{A} and *T*_{B}Δ_{B} assumed by the authors. Numerical calculation of the maximal eigenvalue of the projection matrix allows the relaxing of the assumption that the interphenotypic switching rates are equal in the two directions. If switching is the optimal strategy then the rate of switching from phenotype A to B is approximately 1/*T*_{A} and the rate of switching from phenotype B to A is approximately 1/*T*_{B}. A full computational investigation yields a phase diagram that is qualitatively the same as the one in Figure 3.

## DISCUSSION

The manner in which microorganisms achieve optimal long-term fitness has become a question that now transcends purely theoretical considerations. It is established that the medical threat posed by the ability of bacterial species to resist antibiotic agents by phenotypic heterogeneity is significant (Lewis 2005; Dubnau and Losick 2006). Understanding the underlying ecological and evolutionary context behind phenomena such as stochastic switching is a research target complementary to uncovering the details of particular genetic and environmental mechanisms (Davidson and Surette 2008). Closed-form analytical results such as those presented here help to reveal the fundamental evolutionary driving processes, which can then inform experimentally—and ultimately medically—motivated studies.

Most previous theoretical studies have assumed equal selection pressures in the different environments (Lachmann and Jablonka 1996; Thattai and van Oudenaarden 2004). Kussell and Leibler (2005) assumed combinations of large enough environmental durations and selection pressures in each environment that guarantees that equilibrium population structure is closely approached in each environmental epoch. Numerical experiments by Salathé *et al.* (2009) relaxed the assumption of symmetry of selection pressure between the environments, allowed the selection coefficient to vary and to take very small values, and found the evolutionarily stable switching rates for environmental durations of 20 generations. They concluded that even small asymmetry in selection pressure between the two environments renders stochastic switching a suboptimal strategy unless the selection pressures are large.

Our mathematical analysis confirms that asymmetry in the selection pressure between the environments can select against the evolution of stochastic interphenotype switching if the selection pressure is small in at least one of the environments and shows that the parameter range for which homogeneity is the favored strategy depends not only on the strength and asymmetry of selection, but also on the duration and asymmetry of environmental durations. We use our exact forms to predict when switching is favored and when it is not and demonstrate why the previously described transition between the two is sharp and discontinuous.

This work generalizes previous results by considering asymmetry in both selection pressures and environmental durations, *i.e.,* by considering a less idealized ecological scenario. We have encapsulated our results by plotting a phase diagram demonstrating the boundary between the optimality of switching and nonswitching and demonstrated the reason behind the unusual property that the change between switching and nonswitching is sharp, discontinuous, and first order (in the language of phase transitions) in nature. This raises the possibility of evolutionary hysteresis: continuous changes in environmental parameters across some well-defined threshold will induce a sudden sharp evolutionary response, and this evolutionary response will not be reversed if the environmental parameters are subsequently changed in the opposite direction across the threshold. Our results for the case where the environmental durations are also asymmetric leads to further unusual features in the phase diagram: the strategy of remaining in the phenotype that is fittest in the longer-duration environment, rather than switching stochastically, is evolutionarily preferred when the environments are symmetric, or near symmetric. We note that these findings are consistent with the experimental results of Beaumont *et al.* (2009) (which used very different methods of imposing selection pressure). The predicted existence of evolutionary hysteresis effects is similarly experimentally verifiable.

Bacterial persistence (Balaban *et al*. 2004; Dhar and McKinney 2007) is a well-known mechanism that depends on stochastic switching between normal and persister phenotypes, where persisters are able to survive better in presence of antibiotics but are virtually growth arrested in growth environments. On the basis of the observed switching rates between the two phenotypes in wild-type *Escherichia coli* and using an exponential growth model, Kussell *et al.* (2005) estimate that the duration of growth environment is on the order of a hundred years whereas the duration of antibiotic environments is on the order of days, which means that it is four orders of magnitudes lower than that of growth environments. In wild-type *E. coli,* the growth rate difference between the phenotypes in the antibiotic environment (Δ_{B} = 3.6 hr^{−1}) was measured to be higher than it was in the growth environment (Δ_{A} = 2 hr^{−1}) (Balaban *et al.* 2004). Then the environmental duration, in which selection pressure is lower, is large, while the environmental duration, in which the selection pressure is higher, is much shorter. These observations are consistent with our theoretical predictions of a stochastic switching strategy being favored.

Kussell and Leibler (2005) showed that if the environmental durations are large, stochastic interphenotype switching can be a better strategy than sensing the state of environment and changing phenotypes in response to it if the cost of sensing is taken into account. As these authors pointed out, for changes in the environment that occur on shorter time-scales (tens of generations or less) sensing systems and responsive switching tend to be employed. In this work we showed that if the environmental period is large, switching is also better than phenotypic homogeneity for a large range of selection pressures and selection pressure differences between environments. For the limit of large selection pressures and large environmental periods assumed by Kussell and Leibler (2005) we recovered their result that for two-state systems the optimal switching rates are approximately equal to the rates of environmental variation, and we conclude that in such circumstances, stochastic switching is a better strategy than homogeneity. It remains an open question as to how such strategies are altered in growth-limited environments, *e.g.*, in the presence of density dependence, so that long-term fitness is no longer measured purely by growth rate in the initial exponential growth phase. It is well established that bacteria follow approximately logistic growth in a finite environment and the optimal strategies will vary according to the available carrying capacity at a particular point in the growth phase; in the notation of our model, Δ becomes a variable parameter during growth. Preliminary numerical investigation of Equations 1 and 2, but modified by including (1 − (*n*_{A} + *n*_{B})/*C*) factors to impose logistic growth, reveal no new phenomena from those identified by analysis of Equation 6, but an exhaustive investigation is beyond the scope of this study.

## Acknowledgments

The authors are grateful to the anonymous reviewers for their thorough reviews and helpful comments and to the editor for expediting this process. The authors acknowledge financial support from Engineering and Physical Sciences Research Council grant number EP/F032749/1. A.J.W.'s research is funded by a Research Councils United Kingdom academic fellowship.

## Footnotes

Communicating editor: L. M. Wahl

- Received December 21, 2009.
- Accepted January 12, 2010.

- Copyright © 2010 by the Genetics Society of America