Stabilization of Microbial Communities by Responsive Phenotypic Switching

Clonal microbes can switch between different phenotypes and recent theoretical work has shown that stochastic switching between these subpopulations can stabilize microbial communities. This phenotypic switching need not be stochastic, however, but could also be in response to environmental factors, both biotic and abiotic. Here, motivated by the bacterial persistence phenotype, we explore the ecological effects of such responsive switching by analyzing phenotypic switching in response to competing species. We show that the stability of microbial communities with responsive switching differs generically from that of communities with stochastic switching only. To understand the mechanisms by which responsive switching stabilizes coexistence, we go on to analyze simple two-species models. Combining exact results and numerical simulations, we extend the classical stability results for the competition of two species without phenotypic variation to the case in which one species switches, stochastically and responsively, between two phenotypes. In particular, we show that responsive switching can stabilize coexistence even when stochastic switching on its own does not affect the stability of the community.


I. INTRODUCTION
One of the classical results of theoretical ecology is that large random ecological communities are very likely to be unstable [1].This is a statistical result, based on the analysis of random matrices representing the Jacobians of otherwise unspecified population dynamics [1].It may appear to contradict Nature's large biodiversity, but actual biological interactions are not random: Rather, they are the product of a long history of evolution, during which population dynamics pruned a possibly much larger set of species [2].Nevertheless, the mathematical constraints from random matrix theory that restrict the evolution of this biodiversity can be revealed by this statistical take on population stability, because it enables some analysis of generic large ecological communities and their huge parameter space, and therefore complements the exact results for small systems that are not available for these larger systems.The power of this statistical approach has been demonstrated in a large body of work which, for example, analyzed the effect of the interaction type and structure on the stability of the community [3][4][5][6], revealed the stabilizing effect of higher-order interactions and explicit resource dynamics [7,8], or explored yet other related problems [9][10][11][12][13].
In this context, recent theoretical work has revealed that subpopulation structure such as phenotypic variation can stabilize microbial communities [14,15].In particular, we have argued that abundant phenotypic variation is generically destabilizing [15], essentially because introducing phenotypic variation increases the effective number of species in the system, which is known to be destabilizing [1,3].More subtly however, stochastic switching to a rare phenotype such as the bacterial persister phenotype [16][17][18] can stabilize communities [15].
Although such stochastic phenotypic switching is optimal in infrequently changing environments [19], frequent environmental cues, be they biotic or abiotic, favor responsive phenotypic switching associated to some kind of sensing mechanism [19].Indeed, recent experimental evidence suggests that such sensing is implicated in the bacterial stress response and, in particular, in switching to stress-resilient phenotypes akin to bacterial persisters [16][17][18].For example, formation of persisters under stress has recently been associated with production of the "alarmone" ppGpp, suggesting a stress-dependent persistence response [16,18].Data showing that sublethal antibiotic concentrations increase persister concentration [20] are also consistent with responsive contributions to the persistence switching rates.(While early experiments exposing Escherichia coli to antibiotics had suggested that direct sensing is absent from the persistence response [21], these could not, as already noted in Ref. [19], exclude effects such as a dependence of switching rates on antibiotic concentration.)The cues for this responsive switching may be toxins produced by other bacterial species, such as the colicinogens produced by certain strains of E. coli that have garnered attention in the context of bacterial "rock-paper-scissors" games [22][23][24][25].Indeed, very recent work [26] provides experimental evidence linking such toxins to a stress response: a "suicidal" subpopulation of Pseudomonas aeruginosa actively migrates up antibiotic gradients while upregulating the release of its own bacteriocins, suggesting an attack response against toxin-producing competitors [26].Responsive mechanisms are thus starting to be appreciated as a feature of microbial populations, but their ecological role in microbial communities remains unclear.
Here, we address these ecological implications of responsive switching theoretically by analyzing its effects on the stability of competitive microbial communities with a rare, slowly growing and weakly competing, persister-like phenotype.In the first part of this paper, we extend the two-phenotype model of Ref. [15] to include responsive phenotypic switching, and show how the statistical stability properties of this model dif-fer from those of a model with stochastic switching only, even if the second phenotype is rare.These statistical results emphasize the importance of the type of phenotypic switching for stability, but leave unanswered the question: Under which conditions does phenotypic switching, be it stochastic or responsive, stabilize or destabilize this community?We address this question in the second part of this paper: We extend the classical results for the stability of Lotka-Volterra systems [27] by analyzing a minimal model of two competing species in which one species switches, both stochastically and in response to the other species, between two phenotypes.Using numerical simulations and by deriving exact results for still simpler models, we show in particular that responsive switching can promote coexistence even in cases in which stochastic switching on its own does not affect the stability of the community.

II. STATISTICAL STABILITY OF RESPONSIVE PHENOTYPIC SWITCHING
In this Section, we introduce a model for the competition of  species that switch, both responsively and stochastically, between two phenotypes.This model extends that of Ref. [15].We show that its statistical stability properties are generically different from those of the corresponding model with stochastic switching only.

A. Model
We consider the competition of  2 well-mixed species that have two phenotypes, B and P, each, between which they switch stochastically.In addition, the B phenotypes of each species respond to the other species by switching to the corresponding P phenotype [Fig.1(a)].We denote by   and   the respective abundances of the B and P phenotypes of species , for 1  .With Lotka-Volterra competition terms [27], the dynamics of the vectors  = ( 1 ,  2 , . . .,   ) and  = ( 1 ,  2 , . . .,   ) are thus [28] where ,  are growth rates, the nonnegative entries of the matrices C, D, E, F are competition strengths, , ℓ are nonnegative rates of stochastic switching, and R, S are nonnegative rates of responsive switching.The diagonal entries of R, S vanish, so that microbes do not switch phenotype in response to the presence of other microbes of their own species.We stress that Eqs.(1) are deterministic, as are all subsequent models in this paper: in particular, they do not resolve the individual, stochastic switching events, but only their mean behavior expressed by the deterministic rates of stochastic switching.The functional form of the responsive switching rates in Eqs.(1), with a direct dependence on the competitor abundances, implies a neglect of the dynamics of the chemical cues of the responsive switching.This is justified since our model aims to elucidate the effect of this responsive switching on stability, rather than the effect of these chemical dynamics that was analyzed, e.g., in Ref. [29].The same argument justifies, for instance, neglecting resource dynamics in spite of their known effect on stability [8].The linear dependence of responsive switching on the abundances of other species is, however, a simplifying assumption, just as the logistic Lotka-Volterra interaction terms [27] in Eqs.(1) are the simplest choice of interaction terms.While the full dynamics of the system (and in particular, the existence of non-steady-state attractors like limit cycles) are expected to depend on the details of the functional forms of the interaction terms and switching rates, any such functional forms linearize to the logistic interaction terms and linear switching rates in Eqs.(1) close to an equilibrium.These details do not therefore affect the stability properties of equilibria.Moreover, as shown in Appendix A, the logistic nonlinearities in Eqs.(1) are sufficient for their dynamics to be bounded.For these reasons, we believe these simplifying assumptions to be appropriate for the qualitative analysis of phenotypic dynamics and in particular of the stability of the equilibria of Eqs.(1) in this paper.The parameter  in Eqs.(1) scales the growth and competition of the P phenotypes, and the stochastic and responsive switching rates into the P phenotypes.In the limit  1, the P phenotype is therefore a slowly growing and weakly competing phenotype such as bacterial persisters [16][17][18].

B. Reduced and Averaged Models
How does the stability of a microbial community with responsive phenotypic switching differ from that of a community with stochastic switching only?To answer this biological question and thus understand the ecological effects of the type of phenotypic switching, it is tempting to compare the mathematical stability of equilibria of Eqs.(1) to that of equilibria of a corresponding reduced model with stochastic switching only, in which R = S = O.Such a comparison does not however answer our biological question, because the populations at equilibrium of Eqs.(1) and this reduced model are in general different.This direct comparison cannot therefore decide whether stability differences result from this difference of populations or from the differences in phenotypic switching.
To understand the two types of phenotypic switching, we do not therefore compare Eqs.(1) to this reduced model, but instead to an averaged model that has the same population at equilibrium as Eqs.(1).While the reduced model takes values of interaction parameters or switching rates from Eqs. ( 1), the averaged model replaces these values with effective values that are determined by the condition of equality of populations at equilibrium.We introduced such averaging of models in Ref. [15].There, we pointed out that equality of populations at equilibrium does not lead to simple relations between the eigenvalues of the corresponding Jacobians, so there is no reason to expect their stability properties to be the same.The present discussion develops these ideas.
While stability differences between Eqs. ( 1) and the corresponding averaged model thus reveal the ecological roles of responsive and stochastic phenotypic switching, stability differences between the reduced and averaged models stem from the differences of the reduced and the averaged, effective model parameters.In an actual biological community, these parameters can evolve independently from the evolution of a phenotypic substructure, and so stability differences need not result from this phenotypic substructure.This raises interesting questions akin to those asked in Ref. [19]: For example, does evolving a more complex phenotypic substructure require more evolutionary adaptations than evolving the effective parameters directly?Answering such questions requires, however, coupling Eqs. ( 1) to an evolutionary model.This is beyond the scope of this paper, in which we will therefore focus on the stability differences between Eqs. (1) and the corresponding averaged model, which can be imputed directly to responsive phenotypic switching.
After this rather abstract discussion of reduced and averaged models, we now write down the averaged model with stochastic switching only: This is the model that we have analyzed in Ref. [15].To establish the correspondence between Eqs. (1)

C. Results
In the spirit of the random matrix approach to ecological stability, we compare the stability of Eqs.(1) and Eqs.(2) by sampling their parameters randomly and computing stability statistics.Since the coexistence equilibria of Eqs.(1) and Eqs.(2) cannot be found in closed form, we cannot sample the model parameters directly; rather, as discussed in more detail in Appendix A, we sample the coexistence state itself and some model parameters directly, leaving linear equations to be solved for the remaining parameters to ensure that the chosen coexistence state is steady [15].
For these random systems, we find that, as the number  of species increases, stable coexistence states of Eqs.(1) are increasingly unlikely to be stable with the dynamics of Eqs.(2), and vice versa: Stable coexistence with responsive switching is increasingly unlikely to be stable with stochastic switching only, and vice versa [Fig.2(a)].This trend persists for  1, although the probabilities are reduced in magnitude [Fig.2(b)].Responsive switching can thus stabilize and destabilize equilibria of Eqs.(1).However, coexistence is less likely to be stable with responsive switching than with stochastic switching only [Fig.2(c)], although the probabilities are nearly equal for small  [Fig.2(c), inset].
In the limit  1, coexistence states can be determined in closed form by asymptotic expansion in the small parameter .
as derived in more detail in Appendix A. We thence obtain asymptotic expressions for the Jacobians (Appendix A), similarly to calculations in Ref. [15].Sampling asymptotic coexistence equilibria using this solution, we confirm our findings above, although results differ at a quantitative level because of the different sampling methods [Fig.2(b)].
The models with and without responsive switching thus generically lead to different stability results.The Jacobians at E, J * with responsive switching and K * with stochastic switching only, are related by which follows from the calculations in Appendix A. We stress that this is an exact result and not an asymptotic approximation.Even in view of this linear relation, the fact that the two models give different stability results is not fundamentally surprising [15]: this simply reflects the fact that linear relations between matrices do not imply linear relations between their determinants.What is perhaps more unexpected is that, although J * = K * +  (), the stability results differ even in the limit  1.We have previously related behavior of this ilk to the possibility of eigenvalues with small real parts being stabilized or destabilized by higher-order terms in the expansion [15].same plot, focused on small probability differences.Probabilities were estimated from up to 5 • 10 8 random systems each.Parameter values:  = 1 and  = 0.01 [30]; for the latter value, both exact and asymptotic equilibria were computed.Error bars are 95% confidence intervals [31] larger than the plot markers.
We had to restrict to small values of  when computing the numerical results in Fig. 2, because it becomes increasingly difficult to sample increasingly rare [1,3] stable systems as  increases.The strength of the classical random matrix approach [1,3] to stability is that it can often circumvent this difficulty by analyzing random Jacobians without further specification of the underlying population dynamics and hence of their equilibria.This simplification is not however available for the questions addressed here, because we need to relate "full" and averaged models (in this case, with and without responsive switching).Nonetheless, the trends in Fig. 2 suggest that the stability differences found numerically there are amplified for larger values of .
There is, however, a more important limitation of this stability analysis: one may argue that instability of an equilibrium is not biologically or ecologically significant, because it does not imply extinction of a species, since the dynamics, perturbed away from the unstable equilibrium, may converge to a different equilibrium, a limit cycle, or a more complex attractor, and so the species may still coexist.A more relevant question is therefore: does stability of one model (with responsive switching or with stochastic switching only) not only fail to predict stability, but also fail to predict coexistence in the other model?Unlike predicting stability, predicting coexistence requires the full dynamics of the system and its non-steady-state attractors, and so the predictions depend, as noted above, on the details of the interactions and switching rates.Analyzing the effect of these details on coexistence is beyond the scope of this paper.Here, we restrict to the logistic interactions and linear switching rates in Eqs.(1), evolve unstable equilibria numerically and, as described in Appendix A, determine whether the species coexist permanently [32].
Figure 3 shows the distributions of the long-time dynamics of equilibria that are unstable with responsive switching or with stochastic switching only, respectively, but are stable in the other model.The possible long-time dynamics are permanent coexistence or extinction of some species; in both cases, we distinguish between convergence of the remaining species to an  (2)], of extinction of some species of a random system perturbed away from an equilibrium that is stable in the other model, as a function of the number of species in the system, .(b) Same plot, but focused on low probabilities.Exact equilibria were obtained for  = 1 and  = 0.01 [30].Each probability was computed from up to 5 • 10 8 random systems and up to 5 • 10 3 random systems having an equilibrium with different stability in models R and S. Thick error bars correct the estimated probabilities for the systems in which the numerical solution of the long-time behavior did not converge (Fig. 3); thin error bars add 95% confidence intervals [31].Only those error bars larger than the plot markers are shown.equilibrium or to a limit cycle.The small proportion of systems for which the numerical solution does not converge (Fig. 3) may include systems in which more complex attractors arise.We note in particular that extinction of some species is not a rare outcome (Fig. 3): it is actually the most likely outcome if both phenotypes have similar abundances [ = 1, Fig. 3(a)], while convergence to a limit cycle of all species is more likely if one phenotype is rare [ 1, Fig. 3(b)].From these distributions, we estimate, for both models, the probabilities of extinction of some species for systems perturbed away from an equilibrium that is stable in the other model (Fig. 4).The increase of these probabilities with the number of species  qualitatively matches the increase of the probabilities of an equilibrium that is stable in one model being unstable in the other model (Fig. 2).We can therefore extend our previous conclusion: stability of an equilibrium in one model does not even in general imply permanent coexistence in the other model.
Finally, we could similarly compare the stability of Eqs.(1) to that of an averaged model without phenotypic switching (Appendix A).However, in view of the correspondence of Eqs.(1) and Eqs.(2), the qualitative conclusions of such a comparison must parallel those of comparing models with stochastic switching only and averaged models without phenotypic variation.These two models we compared in Ref. [15], where we concluded that stochastic switching to an abundant phenotype is destabilizing, but stochastic switching to a rare phenotype, corresponding to  1, is stabilizing.Here, we can therefore conclude similarly that responsive switching to an abundant phenotype generically destabilizes the community (compared to the case in which there is no phenotypic variation), but that responsive switching to a rare phenotype has a stabilizing effect.
However, all of these results are fundamentally statistical in nature.While they show how responsive phenotypic switching affects stability on average, they do not yield any insight into the conditions under which phenotypic switching stabilizes or destabilizes the community.To understand the mechanisms underlying these stability differences, we therefore complement this statistical analysis with an analysis of reduced twospecies models in the next Section.

III. STABILIZATION OF TWO-SPECIES COEXISTENCE BY RESPONSIVE SWITCHING
In this Section, we reduce the -species model (1) to a minimal two-species model of responsive phenotypic switching.In the context of this model, we analyze the mechanisms by which stochastic and responsive switching affect stability numerically.We confirm some of these numerical results by deriving exact results for simplified models in Appendix B. These exact results extend the classical results [27], which we rederive in Appendix C, for the stability of the two-species Lotka-Volterra competition model.

A. Minimal two-species model
The simplest setting for responsive phenotypic switching is the competition of two well-mixed species [Fig.1(b)].The first species has a single phenotype, A, while the second one has two phenotypes, B and P, between which it switches stochastically.Moreover, phenotype B responds to the competitor phenotype A by switching to P [Fig.1(b)].To interpret this model in terms of actual biological systems, it is useful to think of P as a phenotype resilient to stress conditions, akin to the bacterial persister phenotype [16][17][18], and to which the common, "normal" phenotype B switches in response to toxins produced by A. With this interpretation and for simplicity, we shall sometimes refer to B, P, A and their respective abundances , ,  as bacteria, persisters, and competitors [33], and we shall sometimes assume the model parameters to scale accordingly.
Thus, again using Lotka-Volterra competition terms [27], the nondimensionalized dynamics of the system are wherein , , , , , , , , , , , ,  0 are dimensionless nonnegative parameters.To obtain this form of the equa- tions, we have scaled time and the population sizes so as to remove the parameters that would otherwise appear in the logistic growth term (1 − ) of the bacteria in the absence of persisters and competitors [34].The interpretation of the model parameters in Eqs. ( 4) is given in Table I.
We have not explicitly introduced a positive parameter  1 scaling the competition dynamics of and switching rates to persisters, but we expect , , , ,  =  (1), , , , , , , ,  =  () from the interpretation of the model in terms of bacteria, persisters, and competitors and by comparison of Eqs. ( 1) and (4).We will not assume  1, but we will sometimes invoke  1 below and in the calculations in Appendices B, C, D to restrict parameter ranges or impose inequalities between parameters.

B. Results
In this subsection, we analyze model (4) numerically and ask: under which conditions and to what extent does responsive switching stabilize or destabilize coexistence?

Numerical setup
Although the full model ( 4) is too complicated for meaningful analytical progress to be made, its equilibria can be found numerically and efficiently by precomputing, using M -, the exact polynomial equations satisfied by the equilibria from Gröbner bases [35].For each computed equilibrium, we check its accuracy using the values of the right-hand sides of Eqs.(4) evaluated there, and verify that all solutions have been found using a test based on Sturm's theorem [36].Finally, we determine the stability of the computed equilibria using the Routh-Hurwitz conditions [27].Similarly, we determine the stability of equilibria of the averaged models (with stochastic switching only and without phenotypic variation) corresponding to Eqs. (4).Unstable equilibria are integrated numerically, similarly to the calculations in Sec.II, to determine whether the species coexist permanently notwithstanding a particular coexistence equilibrium being unstable.We reduce the number of systems that need to be integrated in this way by showing, in Appendix D, that coexistence is permanent if all trivial steady states are unstable provided that the persister scalings (5) are satisfied.The setup for our numerical calculations is as follows: we fix a population of bacteria and competitors by fixing the competition parameters , , , and compute stability diagrams in the (, ) plane for random choices of the remaining model parameters , , , , , , , , which we constrain to satisfy the persister scalings (5).We report these numerical results by plotting, for each point in the (, ) diagram, stability statistics (such as the proportion of random systems in which coexistence is stable or permanent at this point), but it is important to note that these statistics, while providing a convenient way of visualizing the generic properties of the stability diagrams, have no real biological meaning because the averaged models vary between these random instantiations of the full model, since the bacteria-persister competition parameters are held fixed while the switching parameters and the remaining competition parameters are varied.

Discussion
Figure 5 shows numerical stability diagrams of Eqs. ( 4).The classical results for two-species Lotka-Volterra model (Appendix C) suggest distinguishing the cases / >  and / < , which indeed give rise to qualitatively different stability diagrams (Fig. 5).
In both cases however, sufficient levels of responsive switching destabilize coexistence compared to the averaged model with stochastic switching only.This effect is less pronounced for  < / [Fig.5(a)] than for  > / [Fig.5(b)].This destabilization is that which we already noted when discussing the statistical stability of responsive switching in the previous section: Figure 2(c) has already shown that coexistence is more likely to be stable in models with stochastic switching only than in those with responsive switching.
In the context of the minimal two-species model (4), we can however address the mechanisms underlying this average destabilization, and hence even identify the conditions under which responsive switching can be stabilizing.To this end, we analyze the effect of the P phenotype, i.e. the effect of the parameters , , , , ,  that describe its interactions with the A and B phenotypes (Table I), on the stability and permanence of coexistence: We plot, in Figs. 6 and 7, stability diagrams of Eqs.(4) for the cases in which one and only one among , , , , ,  is nonzero.We emphasize that the these parameters are of the same order, as expressed by the persister scalings (5), so neglecting the effect represented by one (say, persister growth) compared to another (say, persistercompetitor interactions) is not in general ecologically consistent.For this reason, these "one-parameter" models are not "more minimal" than Eqs.( 4), but the "one-parameter" diagrams in Figs. 6 and 7 can reveal the individual mathematical effects of these parameters, and our results will show how they can help unravel mechanisms of (de)stabilization.
The full stability diagram of Eqs. ( 4) is not of course a trivial superposition of the stability diagrams in Figs. 6 and 7. Nonetheless, we recognize features of these diagrams in the stability diagrams of Eqs. ( 4) plotted in Fig. 5.It is therefore all the more striking that, for the same parameter values for which Fig. 5 shows destabilization of coexistence due to responsive switching, the "one-parameter" diagrams show the possibility of stabilization of coexistence: Persister-competitor interactions [ > 0, Figs.6(  This suggests that coexistence can be stabilized by tuning the persister parameters so that the stabilizing parameters ,  exceed the others, while still remaining "small", as required by the persister scalings (5).This is confirmed by the numerical results in Fig. 8.The stabilizing effect of "elevating" ,  is less pronounced for  < / [Fig.8(a)] than for  > / [Fig.8(b)], mirroring the weaker destabilizing effect of "unelevated" ,  in the former case.
How do we interpret this stabilization by elevated ,  in the context of our analysis based on a weakly-competing, persisterlike phenotype P [Fig.1(b)]?The "elevated" growth rate  that stabilizes coexistence with responsive switching might correspond to a phenotype that relies on a different food source (and therefore competes weakly with both the normal phenotype B and the competitors A); the growth rate  may still be small because the alternative food source may be scarce, or more difficult to metabolize.Similarly, an elevated, but small persister-competitor interaction rate  may result from toxins produced by the persisters and acting specifically on the competitors.That the ecological fitness of the phenotypes we have just described should be linked to responsive switching in particular is not surprising.
The elevated growth and and interaction with competitors of this phenotype are consistent with a persister phenotype, because the corresponding parameters remain "small" in the sense of the persister scalings (5).While the bacterial persister phenotype [16][17][18] is sometimes considered to be a dormant phenotype [37], small levels of growth or interaction with competitors are mathematically different from absence thereof.Our theoretical results stress this point by suggesting that even such phenotypes can contribute crucially to stability and permanence of coexistence in a microbial community.

Analytical results
The full model (4) does not allow meaningful analytical progress to be made.However, some of the effects of the parameters , ,  and hence some of the features in their stability diagrams (Fig. 6) and in the full stability diagram of Eqs.(4) in Fig. 5 can be understood analytically.(We are not aware of similar analytical results for the other persister parameters , , .)In Appendix B, we derive and discuss these analytical results, which extend the classical results for the two-species Lotka-Volterra model [27] to the simplest mathematical models of responsive switching.Again, these "one-parameter" mathematical models are not "more minimal" ecological models than Eqs.( 4), but they are sufficiently simple to allow analytical understanding of the numerical results that have led us to identifying the stabilizing phenotype discussed above.We shall postpone the detailed discussion of these analytical results to Appendix B, but emphasize three of their features here.
First, the exact calculations show that there are regions of parameter space in which all steady states of the "one-parameter" models involving , ,  are unstable, but in which coexistence is still permanent.This emphasizes the importance of non-steady-state attractors for permanence of coexistence and stabilization of coexistence by responsive switching for these models and stresses how phenotypic switching increases the complexity of the dynamics: no limit cycles (and no more complex attractors) arise in classical two-species Lotka-Volterra competition model [27].Thus analysis of the linear stability of the steady states (Ref.[27] and Appendix C) provides a complete biological picture of the community for the average model, but does not yield a similarly complete picture of a two-species community with phenotypic variation.
Second, these results and further calculations in Appendix C show that stochastic switching on its own does not affect the stability of coexistence in these models.Responsive switching can thus be stabilizing even when stochastic switching has no effect on stability.This statement can be extended numerically to the other "one-parameter" models and permanence of coexistence.Indeed, the "one-parameter" (de)stabilization diagrams comparing the stability or permanence of coexistence in the full model ( 4) and in its average without phenotypic variation (not shown) are identical to those comparing Eqs. ( 4) and its average with stochastic switching only (Figs. 6 and 7, insets), so stochastic switching on its own does not even affect the permanence of coexistence in these "one-parameter" models.However, these results do not carry over to the full model ( 4): the (de)stabilization diagrams for the comparison of Eqs.(4) and its average without phenotypic variation (not shown) are different from those in Figs. 5 and 8. From these numerical results, we can in fact conclude that elevated ,  as in Fig. 8 also stabilize the averaged model with stochastic switching only compared to that without phenotypic variation.
Finally, the analytical results stress that even small rates of stochastic switching,  =  (), affect stability, but only if the B and A phenotypes are similar enough, as expressed by the scaling requirement  −    2 .Our choice of allowing, in our numerical calculations, "large" values  =  (1), which are not consistent with the persister scalings (5), is therefore one of numerical convenience that does not affect the biological validity of our analysis.The existence of such a requirement is not unexpected: an asymptotically rare phenotype P should not change the stability of communities in which one of the A and B phenotypes strongly dominates the other, i.e.  −  =  (1).This case of strong dominance is perhaps less relevant to our analysis, because one expects such strong dominance to lead to one species simply outcompeting the other.By contrast, weak dominance of one species may build up the evolutionary pressure that could lead to the emergence of stabilizing features such as the responsive switching analyzed here.In this context, our earlier results in Sec.II and in Ref. [15], showing that the effect of a rare phenotype is amplified as the number of species increases, might be relatable to a combinatorial and statistical increase of species interactions lacking such strong dominance.

IV. CONCLUSION
In this paper, we have analyzed the ecological implications of phenotypic variation and, in particular, responsive phenotypic switching in the context of microbial communities in which the species have a rare, slowly growing and weakly competing, persister-like phenotype.We have shown that the sta-tistical properties of stability and permanence of coexistence are different in models with responsive phenotypic switching and in corresponding averaged models with stochastic phenotypic switching only, and we have emphasized the need to define these averaged models carefully.Although this statistical analysis showed that coexistence is less likely to be stable on average with responsive switching than with stochastic switching only, numerical results for a minimal two-species model revealed those parameters (and hence the ecological conditions) in combination with which responsive switching can stabilize two-species communities.Exact results for simplified mathematical models showed further that responsive switching can stabilize coexistence even when stochastic switching on its own does not affect the stability of the community.Additionally, our numerical results emphasized the importance of non-steady-state attractors for coexistence of all species, even for the simplified two-species models, but analytical understanding of these attractors is still lacking.
All of our results thus hint at a complex relationship between ecological stability and phenotypic variation, of which this and previous studies [14,15] have only scratched the surface.Focusing, as we did in this paper, on a minimal twophenotype structure with the asymptotic separation afforded by a rare, persister-like phenotype has enabled more detailed and mechanistic understanding.However, extending these results to many-species systems and more general phenotypic structures remains an important challenge for future work.Reference [14] has already begun to address these questions in the context of stochastic phenotypic switching, but to understand which phenotypic interactions stabilize or destabilize manyspecies systems, further analytical progress, for simple models of communities of more than two species, will be crucial to guide statistical exploration of many-species systems due to the quadratic increases of the number of interaction parameters with the number of species.Such analytical progress may benefit from the information-theoretic approaches used in evolutionary population dynamics [38,39].Moreover, now that a range of stability-affecting structures of microbial communities, from phenotypic variation (this work and Refs.[14,15]) and species modularity [6] to higher-order interactions [7] or resource and signaling dynamics [8,29], has been identified, further theoretical work will need to analyze the interplay of such effects and ask: Which of these effects (if any) dominates the overturning of May's stability paradigm [1] in actual ecological communities?
The deterministic equations in this paper describe wellmixed populations, with the deterministic parameter that we call the stochastic switching rate representing a mean over individual, stochastic events of phenotypic switching.For a rare phenotype such as the persister phenotype, this approximation may break down because persister abundances may be small not only in relative terms [30], but also in absolute terms [40].Future work will need to address more generally such stochastic effects in simulations resolving individual rare stochastic events to disentangle their contributions to the stability of ecological communities and the contributions of their mean deterministic behavior, as studied in this paper.
Meanwhile, our work has predicted that the combination of responsive switching to a rare persister-like phenotype and its slow growth and weak interaction with competitors can favor coexistence in microbial communities.This ecological prediction remains to be tested experimentally and the weak growth, interaction with competitors, and responsive switching rates of persisters remain to be quantified in microbial communities.Some experimentally support for the importance of this parameter combination might already be given by the recent work [26] showing that a non-reproducing subpopulation of P. aeruginosa releases bacteriocins while migrating up antibiotic gradients.More importantly, these experiments emphasize the role of spatial dynamics in the response to antibiotic stress.Future theoretical work will therefore also need to add such spatial dynamics into the models that we have analyzed here.
More generally, identifying experimental systems in which the role of responsive phenotypic switching can be addressed experimentally is another key challenge for future work.In particular, experimental systems for these questions must allow disentangling the effect of responsive phenotypic switching, say, and the effect of the change of the equilibrium population that results from "turning off" responsive switching.This is precisely the problem that we highlighted when discussing reduced models and solved, theoretically, by defining averaged models for stability comparison.By introducing these averaged models, we have been able to make complex, albeit mathematical ecological perturbations.It is this ability that explains the power of theoretical approaches for guiding experimental exploration of problems in ecology.For this reason, the importance of these theoretical approaches mirrors the outstanding importance of microbial communities [41]-and hence that of understanding the biological and mathematical effects that stabilize or destabilize them-for the world that surrounds us.
systems, provides the calculations of the Jacobians of the equilibria of Eqs.(1) and the asymptotic calculations in the limit  1, and defines an averaged model without phenotypic variation for completeness.

Bounded Dynamics of Eqs. (1)
We begin by establishing sufficient conditions for the dynamics of Eqs.(1) to be bounded, and therefore to be realistic biologically.
We assume that, for each  ∈ {1, 2, . . .,  }, the B phenotype of species  satisfies   ≠ 0 or   <   , and that its P phenotype satisfies   ≠ 0 or   < ℓ  .Then, from Eqs. (1), It follows that   < 0 if   <   /  and   >  min  , for some  min  > 0. We will not need an explicit expression for  min  , but can assume without loss of generality that  min  >   /  ; similarly,   < 0 if   <   /  and   >  min  , for some  min  >   /  .These bounds are independent of the other species, and thus show that, irrespective of the initial conditions, the dynamics of (  ,   ) will enter the bounded region in Fig. 9, and remain in that region.
A similar argument shows that the dynamics are bounded if   = 0, but   <   or   = 0, but   < ℓ  (Fig. 9).On identifying parameters appropriately, these conditions also provide sufficient conditions for the dynamics of the twospecies model (4), or the simplified models (B1), (B2), (B3), to be bounded.9. Bounded dynamics of Eqs.(1).For any nonnegative initial conditions, (  ,   ) will eventually enter the pentagonal bounding region (thick lines), and remain in that region for all times.If   = 0 and   <   or   = 0 and   < ℓ  , the sides of the bounding region are modified (dashed lines).

Sampling of Random Systems
We follow the approach that we introduced in Ref. [15] to sample random instantiations of Eqs.(1) and their coexistence equilibria.
In more detail, we choose the competition strengths C, D, E, F, the stochastic switching rates , ℓ, the responsive switching rates R, S, and the equilibria  * ,  * independently from the uniform U [0, 1] distribution.This leaves linear equations to be solved for the remaining parameters , ; to ensure that they are nonnegative, we choose a common random scaling for the switching rates , ℓ and R, S.
Using the asymptotic solution derived below, the model parameters can be sampled directly, i.e. we can sample ,  randomly, and use that solution to compute  * ,  * .To avoid a breakdown of asymptoticity, we sample parameters in the interval [ 1/4 , 1] rather than [0, 1].We also discard those sampled systems for which any component of the right-hand sides of Eqs. ( 1) is greater than .Moreover, we ensure feasibility of  * by sampling  as a linear combination of the (normalized) columns of C.
Finally, to sample exact equilibria (again indirectly) for  1, we adapt our previous strategy by imposing , up to a random   2 correction, to ensure that  =  (1).

Asymptotic Coexistence Equilibria of Eqs. (1) for 𝜺 1
As announced in the main text, we seek an expansion of the coexistence state ( * ,  * ) in powers of  1 by writing Solving at order   0 ,  0 = 0 =⇒  0 = C −1 • , unless det C = 0, which we assume not to be the case.Then, at order   1 ,  1 =  + R •  0  0 /ℓ, and hence  0 C •  1 = 0, which implies  1 = 0, as claimed in the main text.On substituting these results into Eq.(A3), we find in which O is the zero matrix and I is again the identity.In the first two terms of Eq. (A7) and up to smaller corrections, we recognize the Jacobian K * of the corresponding model with stochastic switching only, because We use these expansions to sample coexistence equilibria of random systems with  1 directly, and to determine their stability numerically.

Permanent Coexistence in Eqs. (1)
To determine whether coexistence is permanent where a coexistence equilibrium is unstable, we perturb the system away from that unstable equilibrium, and evolve it numerically using the stiff solver ode15s of M (The MathWorks, Inc.).During the numerical solution, we repeatedly test for convergence to a stable equilibrium or a stable limit cycle.In particular, limit cycles and their stability are determined using the shooting method described in Ref. [42].
There is one particular numerical difficulty associated with this: the numerical integration generally fails, even at stringent tolerances, if the dynamics approach an attractor intersecting one of the planes   = 0 or   = 0. Some species abundances become arbitrarily small (while remaining nonzero at finite times) on such a trajectory.Since such small abundances lack biological meaning, we simply remove these species from the system.In more detail, we declare species  to go extinct at time   = argmin {max {  (),   ()} <  }, where  1 is fixed, and integrate the system constituted by the remaining species for  >   .We choose  = 10 −6 and find empirically that this reduces considerably the proportion of systems for which the numerical integration would otherwise fail.

Averaged Model without Phenotypic Variation
In the same way as we have associated to Eqs.

APPENDIX B: ANALYSIS OF SIMPLIFIED MODELS
In this Appendix, we derive analytical results for three simplified models to establish some of the features seen numerically in the stability diagrams in Fig. 6 and hence Fig. 5.
The reason for introducing these simplified models is that the full two-species model (4) does not easily lend itself to analytical progress.Indeed, on computing a Gröbner basis [35] for the steady-state version of Eqs.(4) using M (Wolfram, Inc.), we find that even just computing the coexistence equilibria of Eqs.(4) requires solving a quartic equation.To enable some analytical progress, we therefore introduce a simplified version of Eqs.(4), Compared to the full system (4), all but one of the competition terms involving the persisters have been removed in this system.We stress that Eqs.(B1) are not the asymptotic limit of Eqs.(4) for slowly growing and weakly competing persisters.Again, Eqs.(B1) are thus not an ecological model, but a mathematical model: including this one persister competition term, while leaving the system amenable to analytical progress, introduces nontrivial behavior.We will also consider two other simplified models that are obtained similarly: and The three simplified models (B1), (B2), (B3) thus correspond to allowing exactly one of , ,  to be nonzero.The equilibria of the analogous simplified models corresponding to the remaining logistic parameters involving persisters, viz., , , are determined by equations that are at least cubic, and hence do not allow much analytical progress.

a. Feasibility of the coexistence equilibria
We now ask whether the coexistence equilibria are feasible.In Table II, we analyze the possible sign combinations of the variables , , ,  defined in Eqs.(B6) and that appear in the coordinates of the equilibria in Eqs.(B5).This shows that only C + is feasible if ,  > 0, while only C − is feasible if ,  < 0. Neither coexistence state is feasible if  < 0,  > 0, but it is possible for both coexistence states to be feasible if  > 0,  < 0 provided that  > 0,  < 0 and that C ± are real (Table II).These conditions reduce to The second, quadratic condition implies that This additional region of feasibility does not arise if  = 0 or  = 0. Indeed, it is immediate from Eqs. (B7) that coexistence is feasible in that case if and only if , ,  =  +  all have the same sign, and hence if and only if ,  > 0 or ,  < 0.

b. Stability of the coexistence equilibria
We now turn to the question of stability of the coexistence equilibria.Before discussing the general case  ≠ 0, we discuss two special cases with  = 0.The Jacobian evaluated at a coexistence equilibrium (, , ) is Stability of the coexistence equilibria if  = 0.In the absence of responsive switching, i.e. if  = 0, and from Eq. (B13), the characteristic polynomial of the Jacobian at C = (, , ), defined by Eqs.(B7), is wherein with  =  − /.The Routh-Hurwitz conditions [27] imply that C is stable only if  < 0. Recalling that C is feasible if and only if , ,  =  +  all have the same sign, it follows that C is stable only if ,  < 0.Moreover, if ,  < 0 and hence  < 0, then  2 > 0 and  1  2 > − 3 ; the second inequality is easily checked by direct multiplication, noting that  <  < 0. The Routh-Hurwitz conditions then imply that C is stable if and only if ,  < 0.
Stability of the coexistence equilibria if  = 0.In the case  = 0, in which the competition dynamics do not involve  directly, we find that the characteristic polynomial of the Jacobian at C = (, , ), defined by Eqs.(B7), still has the form in Eq. (B14), with modified coefficients Similarly to the case  = 0 discussed above, C is stable if and only if ,  < 0 and ĉ1 ĉ2 > − 3 by the Routh-Hurwitz conditions [27].Noting that , ,  are independent of , the latter condition can be written as a quadratic in ,  2  2 +  1  +  0 > 0. Since  0 =  1  2 +  3 > 0 for ,  < 0, this holds for small enough .Moreover, if ,  < 0 and  −  < 0, all the terms in the definitions (B16) are positive, so  1 ,  2 > 0 and hence ĉ1 ĉ2 > − 3  The condition  2 < 0 is, if ,  < 0 and hence  < 0, equivalent with  < .This says that destabilization at large  requires the death rate of competitors at steady state due to inter-species competition to exceed that due to intraspecies competition (Table I).
We begin by noting that a region of instability must arise at large  provided that  < .Indeed, using M to simplify complicated expressions, we find that, for  1, with, in particular,  0 = (1 +  + ) > 0 and  2 = −( − ) (/) 2 .Now, from the persister scalings (5), we expect  < , so that  2 < 0. Since  0 > 0, it follows that C − is unstable at large  if and only if | | is large enough.This condition is different from the one that we obtained above when discussing  = 0, for which we showed that instability must occur at large  if | | is small enough.This emphasizes that the limit  = 0 is singular.
This asymptotic condition for instability at large  is independent of , , and thus of .Hence instability may, but need not occur in the region ,  < 0. If such a region of instability exists then, because A, B are unstable if ,  < 0, all steady states are in fact unstable in this region.This discussion also shows that instability must occur at large  under the conditions described by Eqs.(B10).
Coexistence is stable, however, for small  if ,  < 0 and hence  < 0. Indeed, a straightforward calculation shows that C − → C as  → 0 at fixed  ≠ 0, with C defined as in Eqs.(B7).Since C is stable for ,  < 0, so is C − for sufficiently small  by continuity.

Direct computation then yields 𝑐
The question whether C − is stable more generally under the conditions in Eqs.(B10) requires somewhat more effort.First, we discuss the limit in which ,  1, considering all other parameters to be  (1) quantities.Moreover, we assume that  = / =  (1).With these scalings,  > 0 and  < 0, so the feasibility conditions (B10) reduce to  > 4. We then find Assuming that  >  as discussed above,  > 1.We conclude that there exists a region of parameter space in which coexistence is stable for ,  1 under the condition in Eqs.(B10) if and only if ( − ) < .We also note that ( − ) <  is implied by the condition  < 0.
Finally, we consider stability near the feasibility boundary  =  * defined in Eqs.(B10).By continuity, stability near this boundary follows from stability at the boundary.Now, by definition Δ = 0 and hence, from Eq. (B18),  − 0 = 0 at this boundary, so stability there is equivalent with The discussion around Eqs. (B10) implies that  > | |.Hence  − 1 > 0 if  >  from Eq. (B22a).If  < , then (2 − )  +  > ( + ) > 0, so Eq.(B22b) shows that  − 1 > 0 continues to hold if ( − ) <  <  provided that  < .We expect this to be true from the persister scalings (5).If  >  however, parameter values such that  − 1 < 0 can be found numerically (not shown).All of this shows that coexistence is stable if  > ( − ), provided that  < .Moreover, rearranging Eqs.(B22a) or (B22b) yields in which the value of the coefficient  is of no consequence.Hence, if  < ( − ), then  − 1 < 0 and coexistence is unstable for sufficiently large ; we obtained the same result above.Here, we note additionally that, if  < ( − ), then  > 0, and so the linear and constant terms of the cubic (B22c) are positive, while its cubic term is negative, so it has exactly one positive real root.Since C − is stable for | | sufficiently small, this root of  − 1 corresponds to a point on the feasibility boundary Δ = 0. Thus the boundary of the unstable region intersects the feasibility boundary.

c. Coexistence beyond coexistence equilibria
We now extend these results on the coexistence of the two species described by Eqs.(B1) beyond coexistence at steady state.We begin by analyzing the stability of the trivial equilibria.We go on to discuss the alternative outcomes of extinction of one species and permanent coexistence [32] of the two species by proving that the two species coexist permanently if all trivial steady states are unstable.
Stability of the trivial steady states.First, we determine the stability of the trivial steady states O, A, B defined in Eqs.(B4a), and which are feasible for all parameter values.The Jacobian of Eqs.(B1) evaluated at O is In particular, this Jacobian has an eigenvalue  > 0, so O is always unstable, with small perturbations expelled from the plane  = 0.The stability of the other trivial steady states depends on , .The Jacobian at A is in which the entries left blank clearly do not affect stability.
Similarly, the Jacobian at B is Direct computation of eigenvalues shows that A is stable if and only if  > 0, while B is stable if and only if  > 0. In more detail, B is an attractor in the plane  = 0, but expels orbits out of that plane if and only if  < 0. It follows that, if  > 0 or  > 0, then there exist (feasible) initial conditions with which Eqs. (B1) converge to A or B, and hence lead to extinction of one species, and so permanent coexistence of the two species is not possible in general.
Permanent Coexistence.The above shows that permanent coexistence irrespective of the initial conditions is only possible if ,  < 0 and hence A and B are unstable, assumed henceforth.
From the conditions obtained in Appendix A, the dynamics of Eqs.(B1) are bounded.Hence extinction of one species requires the dynamics of Eqs.(B1) to converge to a (stable) limit set L that intersects the boundary   = 0. We claim such a limit set L cannot exist if A and B are unstable.
Next, we observe that the dynamics of Eqs.(B1) do not allow L ⊂  : if this were the case, the Poincaré-Bendixson theorem [44] would imply that L is (1) a fixed point, (2) a limit cycle, or (3) a connection of equilibria.However, ( 1) is not possible because O and B are both unstable, the latter by assumption; (2) is not possible because a limit cycle would necessarily contain the only interior equilibrium, B, which is impossible because the latter is, as we have noted Extending these arguments, L cannot in fact intersect  , for if it did, then it would contain a connection O → B, which is impossible because, as noted above, the directions transverse to  are unstable for both O and B [Fig. 10(e)].
Hence L must intersect { =  = 0;  > 0}, so must contain A and and two of its connections.Since it cannot contain the connection O → A by the above, it must contain the two other connections of A. However, direct computation of the eigenvectors of the corresponding Jacobian in Eq. (B23b) shows that one of these is not feasible [Fig.10(f)].This is the final contradiction showing that L cannot intersect   = 0.
This argument shows that both species coexist permanently if A, B are unstable.

d. Stability diagrams of Eqs. (B1)
The exact results derived above yield the stability diagrams shown in Fig. 11 for  > 0; we will not discuss the singular case  = 0.They reproduce some of the features of the numerical stability diagrams in Figs.6(a),(b).
Figure 11 shows how the combination of responsive switching and persister-competitor interactions ( > 0) leads to new behavior compared to case in which these effects are absent ( = 0): There are additional regions of feasibility and stability at large enough rates of responsive switching  >  min * , with  min * given by Eq. (B12).Given the persister scalings (5), it is important to note that  min * =  () is possible even if  =  (), provided that  −    2 .This condition expresses the requirement that the intra-species competitions of Coexistence is feasible in the region bounded by the thick solid black lines.Both C + and C − are unstable in the region marked "unstable", but C − is stable (and C + is unstable) in the region marked "stable".The exact boundary of the region of instability that arises at sufficiently large  (grey lines) must be computed numerically; a region of stability at large ,  must exist in case (a), but only exists in case (b) if ( − ) < .The boundary of this region asymptotes to the straight line / ≡  =  * , where  * is defined in Eq. (B21b) and depends only on ( − )/.In the hatched region of parameter space, the trivial steady states A and B are unstable, and the two species coexist permanently.This is also the region of stable steady-state coexistence for an average of Eqs.(B1) with respect to C − without phenotypic variation and with stochastic switching only (Appendix C).In the shaded region of parameter space (which need not exist), all steady states of Eqs.(B1) are unstable.phenotypes B and A be sufficiently close to their inter-species competitions (Table I).
To understand the stabilization of coexistence by responsive switching observed numerically [Figs.6(a),(b), insets], we compare Eqs.(4) to their averages without phenotypic variation and with stochastic switching only.The calculations in Appendix C show that, if  = 0, there is a one-to-one correspondence, both in terms of parameters and in terms of feasibility and stability, between Eqs. (B1) and this averaged model without phenotypic variation.In other words, stochastic switching on its own does not affect stability.However, this correspondence breaks down if  > 0 (Appendix C): Responsive switching stabilizes coexistence in a region in which the competitor growth rate  is sufficiently large.
The conditions for stability of the trivial steady states A and B derived above are independent of the rate of responsive switching .Responsive switching does not therefore help in driving the competitors to extinction, but the above shows that it makes stable steady-state coexistence possible in  > max { + /, /} where extinction of bacteria and persisters is the only possible steady state at  = 0.All of these observations show how the combination of responsive switching and persister-competitor interactions ( > 0) favors coexistence.
The shaded region of parameter space in Fig. 11(a), in which all steady states of Eqs. ( 4) are unstable, stresses the importance of non-steady attractors: Since A and B are thus unstable, coexistence is permanent, but is not at steady state, since C − , C + are unstable, too.In particular, the mathematical observation that responsive switching destabilizes the coexistence equilibrium in this region of parameter space does not contradict the ecological picture of responsive switching stabilizing coexistence that we have painted above: It simply implies that coexistence cannot be at steady state in this case, and hence that responsive switching induces oscillatory population dynamics.
Classifying all attractors (i.e.not only the stable equilibria) of Eqs.(B1) and (even for stable steady states) the initial conditions that lead to them is beyond the scope of this paper.In particular, our result that responsive switching is stabilizing means that, for some initial conditions, Eqs.(B1) must converge to stable steady-state coexistence, while these same initial conditions must lead to extinction of one species in the absence of phenotypic variation, because the two-species model without phenotypic variation has no non-steady-state attractors [27].We have no analytical proof of the non-existence of such attractors for the averaged model with stochastic switching only, so these initial conditions could lead to permanent coexistence (albeit not at steady state) in the averaged model with stochastic switching only.These analytical results cannot therefore exclude that it might be any phenotypic switching, rather than responsive switching specifically, that makes coexistence permanent.That it is indeed responsive switching that stabilizes coexistence must be shown numerically, as we have done in Figs. 6 and 7.

Stability and permanence of coexistence in Eqs. (B2)
The simplified model (B2) has three trivial steady states similar to those of Eqs.(B1) defined in Eqs.(B4a).They are wherein  = 1 + /( − ).Clearly, O and A are feasible, but, letting  =  − , B is feasible if and only if  > 0.Moreover, model (B2) has a single coexistence equilibrium, with wherein If ,  > 0, then  > 0 if and only if  > 0, assumed henceforth.The coexistence state C is then feasible if and only if , ,  have the same sign, which is, from Eq. (B26b), if and only if ,  have the same sign.We also note that the results of Appendix A show that  > 0 is a sufficient condition for the dynamics of model (B2) to be bounded.

a. Stability of the coexistence equilibrium
We now analyze the stability of the coexistence equilibrium.We will assume that  > 0; the case  = 0 is equivalent to the case  = 0 for Eqs.(B1) analyzed in the first part of this Appendix.We will first discuss the case  = 0 before analyzing  ≠ 0. The Jacobian evaluated at C is Stability of coexistence if  = 0.If  = 0, i.e. in the absence of responsive switching, the Jacobian (B27) has characteristic polynomial wherein The Routh-Hurwitz conditions [27] imply that coexistence is stable only if  < 0, i.e.only if ,  < 0 using the feasibility conditions.Conversely, if ,  < 0 and hence  < 0, then  2 > 0 and  1  2 > − 3 , of which the second inequality is easily checked by direct multiplication.The Routh-Hurwitz conditions thus imply that coexistence is stable if and only if ,  < 0.
Stability of coexistence if  ≠ 0. If  ≠ 0, the characteristic polynomial of the Jacobian (B27) still has the form in Eq. (B28), with modified coefficients As in the above analysis of the case  = 0, the Routh-Hurwitz conditions [27] imply that coexistence is feasible and stable if and only if ,  < 0 and ĉ1 ĉ2 > − 3 .
The latter condition is not implied by ,  < 0, although coexistence is stable if || or | | is sufficiently small away from the singular point  =  = 0. We prove this claim by expanding the final Routh-Hurwitz condition in  and , assuming all other parameters to be  (1) quantities.Using M to handle complicated expressions, we obtain wherein  =  +  +  2 , which proves our claim.A region of instability must however arise for ,  1.This follows from the expansion ĉ1 ĉ2 Within this region of instability, ,  < 0, and so A, B are unstable there, too, whence all steady states of Eqs.(B2) are unstable in that region.Moreover, Eqs.(B30) show that ĉ1 →  1 and ĉ2 →  2 as  → 0. Since  1  2 +  3 > 0, it follows that, for sufficiently small , ĉ1 ĉ2 +  3 > 0, too.In other words, for small , C is stable if and only if ,  < 0.
If  < /, then ,  < 0 is possible if  = 0, but, if  > /, this requires  > (1 + /) and  As in the analysis of Eqs.(B1) in the first part of this Appendix, we discuss the stability of the trivial steady states to obtain conditions for permanent coexistence.Again, O is clearly unstable, because the Jacobian of Eqs.(B2) at this equilibrium is with an unstable eigenvalues  > 0. The Jacobian at A is Again, entries that do not affect stability have been left blank.
Similarly, the Jacobian evaluated at B is On computing the eigenvalues of these Jacobians and since we assume that  > 0, we conclude, again, that A is stable if and only if  > 0, while B is stable if and only if  > 0. We note from Eqs. (B26a) that  > 0 is not possible and hence that A is unstable if  > (/ − 1).The geometric properties of these Jacobians are identical to those of the corresponding equilibria of Eqs.(B1) analyzed in the first part of this Appendix; this follows from direct computation of their eigenvectors.Hence the argument deployed there to establish permanence of coexistence carries over to model (B2): if A and B are both unstable, i.e. if ,  < 0, then the two species coexist permanently.These results yield the stability diagrams drawn in Fig. 12 for  > 0; again, we will not discuss the singular case  = 0.These diagrams confirm some of the features present in the numerical stability diagrams in Figs.6(c),(d), too.
The combination of responsive switching and persister growth ( > 0) leads to new behavior compared to  = 0 for  >  min * , with  min * now given by Eq. (B33).Again, 2 .This is precisely the condition that we discussed when analyzing model (B1).
To assess the effect of responsive switching on coexistence in Eqs.(B2), we consider again the corresponding averaged models without phenotypic variation and with stochastic switching only.The calculations in Appendix C determine the range of stability of these averaged models (Fig. 12).Again, they show Coexistence is feasible in the region bounded by the thick solid black lines, and is stable or unstable in the regions marked "stable" or "unstable", respectively.A region of instability (grey lines) arises at sufficiently large , .The exact boundary of this region must be computed numerically.Equations (B37) and (B38c) the necessary and sufficient conditions for the two different behaviors ("or") that are possible in panel (b).In the hatched region of parameter space, the trivial steady states A and B are unstable, and the two species coexist permanently.The solidly hatched region in panel (a) is also the region of stable steady-state coexistence for an average of Eqs.(B2) without phenotypic variation and with stochastic switching only (Appendix C); the dashed hatching in panel (b) signifies that steady-state coexistence cannot be stable in the averaged model for / <  (Appendix C).In the shaded regions of parameter space, all steady states of Eqs.(B2) are unstable.
that stochastic switching on its own does not affect the stability of the system.Moreover, these results emphasize the importance of establishing correspondences to averaged models when comparing the stability properties of different models: the range of competitor growth rates  for which coexistence is stable increases with the rate of responsive switching  initially [Fig.12(a)], but the parameters of the averaged models vary correspondingly, and so responsive switching is neither stabilizing nor destabilizing for small .However, large levels of responsive switching destabilize C if  < / [Fig.12(a)].This destabilization of the coexistence equilibrium does not mean, however, that responsive switching destabilizes coexistence: as shown above, in the region of large ,  where C is feasible, but unstable, the trivial steady states A and B are unstable, too, so the two species still coexist permanently.The fact that all steady states of Eqs.(B2) are unstable there simply means that the two species do not coexist at steady state.Similarly, Fig. 12(b) shows that the combination of responsive switching and persister growth is stabilizing if  > / and the competitors grow sufficiently fast.Again, for large , , the two species cannot coexist at steady state, but our calculations imply that they do coexist permanently.This emphasizes the importance of non-steady-state attractors for coexistence.As already discussed for model (B1), this stabilization argument cannot exclude the possibility of permanent, non-steady-state coexistence in the model with stochastic switching only (which is not possible in the averaged model without phenotypic variation); this question must be addressed numerically, as we have done in the discussion of Figs. 6 and 7.
The condition for stability of B, which we have computed above, is independent of the rate of responsive switching, .Again, responsive switching cannot therefore drive the competitors to extinction.Moreover, we have shown that A, in which the bacteria and persisters are extinct, is stable only if  < (/ − 1).Not only does responsive switching thus destabilize A and hence extinction of bacteria and persisters.It also makes permanent coexistence (albeit not necessarily at steady state) possible at large enough competitor growth rates  >  max {/, } (Fig. 12), where extinction of bacteria and persisters is the only possible steady state in the absence of responsive switching, i.e. at  = 0.All of this supports the idea that the combination of responsive switching and persister growth ( > 0) favors coexistence.

a. Feasibility of the coexistence equilibria
We now have to ask whether the coexistence equilibria C ± and C are feasible.The calculations are similar to those for model (B1) in the first part of this Appendix.In Table III, we analyze the possible sign combinations of the variables , , ,  , defined in Eqs.(B42), that appear in the coordinates TABLE III.Feasibility of the coexistence equilibria C ± of Eqs.(B3): discussion of the sixteen possible sign combinations of , , ,  , defined in Eqs.(B42).For some combinations, the resulting signs of  ± or  ± , defined in Eqs.(B41), are given, and the symbol C is used for some combinations to indicate that the resulting values of  ± or  ± may have nonzero imaginary parts.Some sign combinations, marked # in the final column, are inconsistent with definitions (B42); for other sign combinations, this column gives the corresponding feasibility results.Next, we analyze the stability of the coexistence equilibria.We can assume that  > 0, since the case  = 0 is equivalent to the case  = 0 for Eqs.(B1) discussed in the first part of this Appendix.In now familiar fashion, we discuss the cases  = 0 and  ≠ 0 separately, but note that, in both cases, the Jacobian evaluated at a coexistence equilibrium (, , ) is Stability of coexistence if  = 0.If  = 0, the characteristic polynomial of Eq. (B47) evaluated at C = (, , ) defined in Eqs.(B43) is where In particular, the Routh-Hurwitz conditions [27]  In this Appendix, we briefly rederive the classical results for the stability of two-species Lotka-Volterra competition models [27].We then obtain the stability conditions for averaged models of this Lotka-Volterra form and that correspond to the simplified models (B1), (B2), (B3).We combine these results with the stability calculations in Appendix B to compare the stability of coexistence in these simplified models and in their averages without phenotypic variation and with stochastic switching only.

Coexistence in a two-species Lotka-Volterra model without phenotypic variation
Coexistence in a two-species Lotka-Volterra competition model without phenotypic variation is a classical problem, discussed, for example, in Ref. [27].With the aim in mind of comparing this model to the models with phenotypic variation considered in the main text, it will be useful to rederive the results briefly.We consider the competition between species  ,  described by the differential equations wherein  ,  ,  ,  ,  ,  0 are parameters.As in our derivation of Eqs. ( 4) and (B1), we could have scaled time and  ,  to set some parameters equal to 1, but we have not done so in order to be able to relate this model to a model with phenotypic variation in the next subsection [34].Equations (C1) have a single coexistence state C = ( ,  ), where  =  / and  =  / , with We seek to describe the populations , ,  that evolve according to Eqs. ( 4) by an averaged model without phenotypic variation and two populations  , which corresponds to  and , and  , which corresponds to  [Fig.1(b)].We have introduced this averaging in Ref. [15], and we have motivated it again here, in Section II.

a. Derivation of the averaged model
A coexistence equilibrium C = (, , ) of Eqs. ( 4) is consistent with the equilibrium C = ( ,  ) of Eqs.(C1) if and only if the populations, the births, and the competition are equal at equilibrium, i.e. if and only if = ( +  ),    = ( + ), (C6c) These conditions are to Eqs. ( 4) what conditions (A9) are to Eqs. (1).We let  = /, so that, on introducing  = (1+) −1 , they reduce to Hence Eq. (C4), the stability condition for the averaged model, becomes respectively, which, from the calculations in Appendix B, are precisely the stability conditions for Eqs.(B1), (B2), and (B3), respectively, with  = 0.This establishes the one-to-one correspondence between these simplified models and the corresponding averaged models claimed in the main text.We note that, for model (B3), this requires an additional condition on the model parameters discussed in the analysis of that model; this condition follows from the persister scalings (5).If now  =  =  = 0, but  ≠ 0, then  = ( + )/, where, again,  =  −  > 0. We discuss the three simplified models (B1), (B2), (B3) severally, simplifying the interval (C8) using the explicit expressions for  derived in Appendix B.
For model (B1),  =  = 0 and  =  ± , defined in Eq. (B5a) and corresponding to the equilibria C ± .The stability conditions (C8) of the averaged model become as in definitions (B6a).Since  ± > 0 by feasibility, a necessary condition for stability is ,  < 0, but this is inconsistent with feasibility of C + (Appendix B), which is therefore unstable in the averaged model.Now, using the expression for  − in Eq. (B5a), we find as defined in Eqs.(B42a), too.Since  ± > 0 by feasibility, the second condition is clearly satisfied if  < 0. We may therefore suppose that  < 0 and  > 0. Now, on letting  = ( + )/ +  as in Eqs.(B42b) and using the explicit form of  ± given in Eq. (B41a), the second condition in Eqs.(C14) becomes  where the second inequality holds since ( + ) 2 >  2 +  2 for ,  > 0. This shows that inequalities (C16) hold true.Moreover, < 0, (C18) since  < 0,  > 0. It follows from inequalities (C16) and (C18) that condition (C15) holds true.On comparing with the feasibility conditions derived in Appendix B, we infer that C − is stable in the averaged model if and only if it is feasible.Moreover, C + is stable (and feasible) in its averaged model (which is in general different from that for C + ) if and only if conditions (B44) hold.Appendix B shows that C + is an unstable equilibrium of Eqs.(B3).We therefore emphasize that, while responsive switching ( > 0) thus destabilizes the coexistence equilibrium C + , coexistence in the unaveraged Eqs.(B3) may still be stable at the other coexistence equilibrium C − .These results also enable us to compare the simplified models (B1), (B2), (B3) with responsive switching to averaged models with stochastic switching only: The latter have an effective switching rate  =  + .This follows similarly to the correspondence of Eqs. ( 1) and (2) established in Sec.II and implies that  =  / = ( + )/ as above, i.e. conditions (C9), with  replaced by this  , are precisely the conditions that we have just analyzed.This means that the averaged model with stochastic switching only is stable if and only if the averaged model without phenotypic variation is stable.We stress however that there is no reason to expect this "complete" correspondence to hold for the full two-species model (4).

APPENDIX D: STABILITY OF THE TRIVIAL STEADY STATES OF EQS. (4) AND PERMANENT COEXISTENCE
While we did not obtain any analytical results bearing on the stability (or indeed the feasibility) of the coexistence equilibria of Eqs.(4), more meaningful progress can be made as far as the stability of the trivial steady states is concerned.These results extend our results for the simplified models by proving that coexistence is permanent for model (4) whenever these trivial steady states are all unstable.
We begin by noting that the trivial steady state O = (0, 0, 0) of model ( 4) is always unstable.Next, Eqs.(4) have a trivial steady state A = (0, 0, /), which is feasible for all parameter values.The Jacobian of Eqs. ( 4 wherein, again,  =  −  and entries left blank clearly do not affect stability; we shall assume that  > 0, consistently with the persister scalings (5).This Jacobian has one eigenvalue − < 0, which results from the trivial connection O → A. Classical stability results [27] imply that A is stable if and only if the sub-Jacobian where  = / and, once again, entries left blank do not affect stability.This Jacobian has one eigenvalue  −  − , associated with perturbations out of the plane  = 0, which may be of either sign.We note that the sub-Jacobian respectively.We have assumed, in the final line and consistently with the persister scalings (5), that  > .[It is because of the need for this additional assumption in this argument that we have provided separate proofs of our permanence result for the simplified models (B1), (B2), (B3), which do not require such an additional assumption.]It follows that the remaining eigenvalues at B have negative real parts, and that the stability of B is determined by the sign of  −  − ; determining this sign requires solving Eqs.(D5).For our purposes, it suffices to note that this implies that if B is unstable, then the direction transverse to  = 0 is unstable.Hence the geometric properties of B match those of the corresponding steady state of Eqs.(B1) in Appendix B. Compared to the discussion in Appendix B, there remains however one more case to be discussed before we can conclude that coexistence in Eqs. ( 4) is permanent if all trivial steady states are unstable by the results of Appendix B. Indeed, there could be multiple steady states of the form B = (, , 0), and so, using the notation introduced in Appendix B, the limit set L could also be a connection of several such states.This is however impossible because each of them is stable in the plane  = 0. We can therefore conclude that, also in the full model (4), coexistence is permanent if all trivial steady states are unstable.bacteria at the microscale in soil, PLoS ONE 9, e87217 (2014) estimate that bacteria in soil communities interact with at most  10 3 neighbors, so the low relative persister abundances [30] mean that, even if the community counts many fewer than the  10 4 species typical for soil communities (vide ibid.), interactions with persisters are rare, stochastic events in such communities.[43] The approach of proving permanence by "chasing" possible limit sets in this way has also been used for general three-species Lotka-Volterra systems without phenotypic variation: See, e.g., H. I. Freedman and P.

FIG. 1 .
FIG. 1. Models of stochastic and responsive phenotypic switching.(a) In the model of Sec.II, each species has two phenotypes, B and P, and switches stochastically between them.Moreover, the B phenotype of each species responds to other species by switching to the P phenotype.(b) In the minimal two-species model of Sec.III, the second species has a single competitor phenotype A, which causes B to switch to P. Dashed lines: stochastic switching.Solid lines: responsive switching.

FIG. 5 .
FIG. 5. Numerical stability diagrams of Eqs.(4) in the (, ) diagram in the cases (a) / >  and (b) / < .The color of each point in the stability diagrams represents the proportion of  = 1000 random systems for which coexistence is stable or permanent at that point.The insets plot the proportion of systems for which coexistence is destabilized, , compared to the averaged model with stochastic switching only.The symbol ⊕ in parentheses [panel (b), inset] indicates that a nonzero proportion of systems (too small to be visualizable by the color scheme) is stabilized (or becomes permanent) compared to the averaged model.Parameter values:  = 0.8,  = 1.2,  = 1.1 [panel (a)] or  = 1.9 [panel (b)], , , , , , ,  ∼ U [, 2], with  = 0.1, and  ∼ U [0.8, 1.6].

FIG. 10 .
FIG. 10.Permanent coexistence of the two species described by Eqs.(B1) if A, B are unstable.(a) Orbits of Eqs.(B1) are expelled from  = 0 and  = 0. (b) Non-existence of a (stable) limit cycle L in  = {  = 0; ,  > 0}.(c) Non-existence of a connection L containing B in  .(d) Non-existence of a homoclinic connection of O in  .(e) Non-existence of a connection containing O → B. (f) Non-existence of a connection containing A, because one of the connections of A is not feasible.

|
FIG. 11.Feasibility and stability of the coexistence states C ± of the simplified model (B1) in the (, ) diagram for  > 0, in the cases (a) / >  + / and (b) / <  + /, assuming scalings (5).Coexistence is feasible in the region bounded by the thick solid black lines.Both C + and C − are unstable in the region marked "unstable", but C − is stable (and C + is unstable) in the region marked "stable".The exact boundary of the region of instability that arises at sufficiently large  (grey lines) must be computed numerically; a region of stability at large ,  must exist in case (a), but only exists in case (b) if ( − ) < .The boundary of this region asymptotes to the straight line / ≡  =  * , where  * is defined in Eq. (B21b) and depends only on ( − )/.In the hatched region of parameter space, the trivial steady states A and B are unstable, and the two species coexist permanently.This is also the region of stable steady-state coexistence for an average of Eqs.(B1) with respect to . (B26a).In particular, this implies the lower bound  * >  min * =  −   ( − ).(B33) We are left to discuss the stability of C near the singular point  =  = 0 for  > /.This corresponds, in the (, ) plane, to  =  min * ≡ (1 + /) and  =  min * > 0 for  > /.Near this point, we write   min * = 1 + ζ,   min * = 1 + β, (B34) with ζ, β > 0. Inserting these definitions into Eqs.(B26a) shows that the domain ,  < 0 in which C is feasible corresponds, at leading order, to β >  ζ/, with  =  −  > 0. Similarly, from Eqs. (B26), (B29), and (B30), the stability boundary ĉ1 ĉ2 +  3 = 0 corresponds, again at leading order, to the straight lines β = 0, β =  0 ζ, where c. Stability diagrams of Eqs.(B2)

FIG. 12 .
FIG.12.Feasibility and stability of the coexistence state C of the simplified model (B2) in the (, ) diagram for  > 0, in the cases (a) / >  and (b) / < .Coexistence is feasible in the region bounded by the thick solid black lines, and is stable or unstable in the regions marked "stable" or "unstable", respectively.A region of instability (grey lines) arises at sufficiently large , .The exact boundary of this region must be computed numerically.Equations (B37) and (B38c) the necessary and sufficient conditions for the two different behaviors ("or") that are possible in panel (b).In the hatched region of parameter space, the trivial steady states A and B are unstable, and the two species coexist permanently.The solidly hatched region in panel (a) is also the region of stable steady-state coexistence for an average of Eqs.(B2) without phenotypic variation and with stochastic switching only (Appendix C); the dashed hatching in panel (b) signifies that steady-state coexistence cannot be stable in the averaged model for / <  (Appendix C).In the shaded regions of parameter space, all steady states of Eqs.(B2) are unstable.
It is for this reason that correspondence at equilibrium does not, as we have already noted above, translate to corresponding stability properties.In what follows, we therefore analyze these stability properties to understand the ecological effect of responsive switching.
(2) Eqs.(2), we notice that an equilibrium E = ( * ,  * ) of Eqs.(1) is also an equilibrium of Eqs.(2) if  =  + R •  * + S •  * reduced model could therefore also result, independently from responsive switching, from evolution of the rates of stochastic switching.Despite the correspondence at E of Eqs.(1) and Eqs.(2), it is clear that their dynamics away from E are in general different.

TABLE I .
Parameters of the two-species model (4), representing the interactions in Fig.1(b), and their interpretations.The parameters have been divided into three groups.Only one of the persister logistic parameters appears in each of the simplified models (B1), (B2), (B3).
A1)Consider first the generic case in which   ,   ≠ 0. Then  +   < 0 if   >   /  and   >   /  .Moreover, if   <   /  or   <   /  , then from Eqs. (1), (1)se conditions uniquely determine the effective parameters  and C of the averaged model and its equilibrium  * .Again, the corresponding reduced model would inherit its birth rates  and competition parameters C from Eqs.(1).Stability differences between this reduced model and Eqs.(A8) need not therefore result from the phenotypic substructure, but could also stem from these parameter differences, i.e. evolution of the B phenotypes.We emphasize that this averaging does not imply that the dynamics of the sum  +  resulting from Eqs. (1) are of the averaged form (A8), and so the dynamics of Eqs.(1) and Eqs.(A8) are in general different away from E.

TABLE II .
Feasibility of the coexistence equilibria C ± of Eqs.(B1): discussion of the sixteen possible sign combinations of , , ,  , defined in Eqs.(B6).For some combinations, the resulting signs of  ± or  ± , defined in Eqs.(B5), are given, and the symbol C is used for some combinations to indicate that the resulting values of  ± or  ± may have nonzero imaginary parts.Some sign combinations, marked # in the final column, are inconsistent with definitions (B6); for other sign combinations, this column gives the corresponding feasibility results.
* >  min * = |( − )/ − (B46) This additional region of feasibility does not arise if  = 0 or  = 0. Indeed, similarly to the analysis of model (B1), Eqs.(B43) show that coexistence is feasible in that case if and only if , ,  = ( + )/ all have the same sign, and hence if and only if ,  > 0 or ,  < 0. b.Stability of the coexistence equilibria [27]e A and B are stable if  > 0 and  > 0, respectively.Comparing these parameter ranges for stability, it follows that either exactly one of A , B , C is stable, or A , B are both stable.From arbitrary initial conditions, Eqs.(C1) converge to the stable steady state if is unique; if A , B are both stable, a separatrix through C separates initial conditions converging to A from those converging to B[27].