Subpopulations and Stability in Microbial Communities

In microbial communities, each species may have multiple, distinct phenotypes. How do these subpopulations affect the stability of the community? Here, we address this question theoretically, showing that simple models with subpopulation structure and averaged versions thereof generically give contradictory linear stability results. Specialising to the bacterial persister phenotype, we analyse stochastic switching between phenotypes in detail in an asymptotic limit. Abundant phenotypic variation tends to be linearly destabilising, but, surprisingly, a rare phenotype can have a stabilising effect. Finally, we extend these results by showing numerically that subpopulations also modify the stability of the system to large perturbations such as antibiotic treatments.

In microbial communities, each species may have multiple, distinct phenotypes. How do these subpopulations affect the stability of the community? Here, we address this question theoretically, showing that simple models with subpopulation structure and averaged versions thereof generically give contradictory linear stability results. Specialising to the bacterial persister phenotype, we analyse stochastic switching between phenotypes in detail in an asymptotic limit. Abundant phenotypic variation tends to be linearly destabilising, but, surprisingly, a rare phenotype can have a stabilising effect. Finally, we extend these results by showing numerically that subpopulations also modify the stability of the system to large perturbations such as antibiotic treatments.
In his seminal 1952 paper [1], Turing revealed the patternforming instability that now bears his name: although diffusion is stabilising, coupling a dynamical system that is stable in well-mixed conditions to diffusion can lead to an instability. Mathematically, this instability is a direct reflection of the fact that the determinant of the sum of two matrices is not the sum of their determinants.
Here, we reveal and study an instability associated with subpopulations in microbial communities that arises from the same mathematical principle. Our starting point is one of the most studied, yet simplest models for the dynamics of species competition is the Lotka-Volterra model [2], in which N species with abundances A follow where α is a vector of growth rates and D a matrix of nonnegative competition strengths [3]. If det D 0, there is a unique equilibrium point A * = D −1 · α of coexistence of all N species. This equilibrium is feasible (i.e. A * > 0) if and only if α lies in the positive span of the columns of D.
We now introduce subpopulation structure into these Lotka-Volterra systems: if the N species each have two subpopulations with abundances B, C, they obey where β, γ are growth rates, and the non-negative entries of E, F, G, H are competition strengths [4]. To a coexistence equilibrium (B * , C * ) of such a "subpopulation-resolving" system, we now associate an "averaged" system (1) and a corresponding equilibrium A * by requiring that, at equilibrium, population sizes, births, and competition be consistent, i.e. that The averaged system and corresponding equilibrium associated with Eq. (2) are defined uniquely by these conditions: they determine α and D, and hence an averaged model of the form (1). Also, A * defined in this way is an equilibrium of this averaged model, and is feasible if (B * , C * ) is. We note that this averaging is a property of the model, depending on the interpretation of its different terms.
We analyse the stability of random coexistence equilibria of these averaged and subpopulation-resolving systems [5]. As N increases, a stable equilibrium of the averaged model is increasingly likely to be unstable in the subpopulationresolving model [ Fig. 1(a)]. This reflects the fact that stable equilibria become increasingly rare as the number of species increases [14,15]. It is therefore all the more striking that stable equilibria of the subpopulation-resolving model are also increasingly unlikely to be stable in the averaged model as N increases [ Fig. 1(a)]. The subpopulation-resolving model (2) can be extended to species with M 2 subpopulations each, but increasing the number of subpopulations per species at fixed number N = 2 of species does not significantly affect the probability that a stable equilibrium of the subpopulationresolving model destabilises in the averaged model [ Fig. 1(b)].
We conclude that averaged models not resolving the subpopulation structure generically lead to incorrect stability results. These results share their mathematical origin with the Turing instability: simple algebraic relations between the Jacobians of dynamical systems evaluated at equilibrium points cannot be expected to lead to simple relations between the determinants of these Jacobians and hence the stability of these equilibria.
The toy model has thus revealed how the subpopulation structure can crucially affect the stability of its equilibria. The remainder of this Letter extends these observations to biologically relevant models, taking phenotypic variation in microbial communities as an example of subpopulation structure. In particular, we focus on the bacterial persister phenotype [12,16]: bacteria such as Escherichia coli switch between a normal growth state and a persister state, in which they significantly suppress growth, but are resilient to conditions of stress, such as competition, but also exposure to antibiotics [12,17,18]. In humans, this can lead to infections that are difficult to treat, even in the absence of genetic antiobiotic resistance; for this reason, the persister phenotype has received much attention in the biomedical community [12,17].
Incorporating this switching between bacteria and persisters into a phenotype-resolving model, we extend Eq. (2) to wherein κ, λ are switching rates. This stochastic form of switching has previously been used to study a single species of bacteria and persisters in the absence of competitive interactions [13,19]. Since these rates are balanced, given an equilibrium of Eq. (4), the consistency conditions (3) still define an equilibrium of the corresponding averaged model (1). Steady states of Eq. (4) cannot be found in closed form; we must therefore sample random systems indirectly [5]. As in our earlier analysis of the toy model, we find that random stable coexistence states of either this phenotype-resolving model or the corresponding averaged model are increasingly unlikely to be stable in the other model as the number N of species increases (Fig. 2).
This setup of the model does not however take into account the weak growth and competition of the persister phenotype or the large separation of switching rates [19]. We therefore modify Eq. (4) to become where ε 1 is a small parameter, so that B are bacteria and C are persisters. For wild-type E. coli, ε ≈ 10 −5 [19], but here, we take ε = 0.01 for numerical convenience, since the analysis below shows that qualitatively identical results are expected for smaller values of ε. A more intricate asymptotic separation arises for the hipQ mutant of E. coli [19], but we do not pursue this further here.
The separation of growth and competition terms and switching rates in Eq. (5) allows for steady states to be found by expansion in ε. Writing B * = B 0 + εB 1 + O ε 2 and C * = C 0 +εC 1 +O ε 2 , we find B 0 = E −1 · β and C 1 = κ B 0 /λ, but B 1 = C 0 = 0 [5]. This limit has thus led to the expected asymptotic separation of the sizes of the bacterial and persister populations: (very) few cells are persisters, at least under laboratory conditions [19]. This asymptotic solution will enable us to uniformly sample all model parameters directly [5].
To analyse the stability of equilibria of Eq. (5), we expand the Jacobian of the full model, J * = J 0 + εJ 1 + O ε 2 , finding with I the identity and O the zero matrix [5]. The Jacobian of the averaged model is Since J 0 is block-upper-triangular, its eigenvalues are those of −λI, which are stable, and those of −B 0 E = K 0 . Hence any unstable eigenvalues of J * and K * are equal to lowest order in the expansion. Equivalently, at this order, the full, phenotyperesolving model (5) is stable if and only if the corresponding averaged model is stable. This conclusion is not borne out however by the analysis of realisations of random systems: as N increases, the probability that a random stable equilibrium of the phenotype-resolving model is destabilised in corresponding averaged model still increases (Fig. 2), although the probability is reduced compared to the previous case. Much more strikingly, the probability that a stable equilibrium of the averaged model is destablised in the phenotype-resolving low high 3. Stability of equilibria of the phenotype-resolving model (4) and the corresponding averaged model, for ε = 0.01. (a) Distribution of eigenvalues, for N = 10, of the Jacobian K * of the averaged system. Histogram obtained from 10 5 random systems. Inset shows distribution of real eigenvalues. (b) Corresponding plot for the Jacobian J * of the full system. Boxes on the real axis and in the inset show the contribution to the distribution of real eigenvalues from −λI. (c) Spectral distribution of K 0 featuring a stable outlier eigenvalue, and a circular core of radius r ∼ √ N: the stability of eigenvalues ν 0 with small real part, |Re(ν 0 )| ε r, can be changed by higher-order terms. The average distance between eigenvalues is δ = O(1). model is vastly reduced (Fig. 2, inset). This is all the more surprising as we argued earlier that the opposite behaviour was to be expected for the toy model since larger systems are more likely to be unstable. We have also sampled exact equilibria of Eq. (5), similarly to our analysis of Eq. (4) above, yielding results in qualitative agreement with those based on the asymptotic equilibria (Fig. 2, inset). We are therefore justified in basing the detailed analysis of the destabilisation mechanism on the asymptotic results, but we shall use these exact equilibria in the final part of this Letter.
To explain these surprising results, we thus return to the asymptotic results and analyse the spectral structure of the Jacobians in more detail. The eigenvalues of K * approximately lie within a circle, as expected from the circular law [20], with the exception of single real eigenvalue that is large and negative [ Fig. 3(a)]. The spectral distribution of J * is the sum of this distribution and the (uniform) distribution of eigenvalues of −λI [ Fig. 3(b)]. The outliers of the spectral distribution can be analysed in great generality [21], but, for our purposes, a heuristic approach is useful: denoting by −µ 0 < 0 the mean of the distribution of entries of K 0 and neglecting correlations between entries, each row of K 0 has approximate sum −N µ 0 , and so K 0 has an approximate eigenvector 1 = (1, 1, . . . , 1) with eigenvalue −N µ 0 < 0, as argued in the Supplemental Material of Ref. [15]. The Perron-Frobenius theorem [9] immediately implies that this eigenvalue is indeed real. The other eigenvalues of K 0 [Fig. 3(c)] are uniformly distributed on a disk of radius r ∼ √ N for N 1 by the circular law [20]. An eigen-value ν 0 has small real part, |Re(ν 0 )| ε, with probability ∼ (εr) r 2 ∼ ε √ N, i.e. = cε √ N, for some c = O(1). The average distance δ between eigenvalues is determined by Nδ 2 ∼ r 2 , so δ = O(1) and the eigenvalues of K 0 are pairwise different at order O(1) [Fig. 3(c)]. Hence, expanding asymptotically, if ν 0 is an eigenvalue of K 0 , then K * has an eigenvalue ν * = ν 0 + εν 1 +O ε 2 , with ν 1 = O(1) a diagonal entry of K 1 in the eigenvector basis of K 0 [5]. Hence ν * is stable if either (i) Re(ν 0 ) < 0 and |Re(ν 0 )| ε or (ii) |Re(ν 0 )| ε and the small real part of ν * is stabilised by K 1 . By definition, |Re(ν 0 )| ε with probability , so (i) has probability (1 − )/2. Let q denote the probability of stabilisation by K 1 in case (ii). Summing over the N − 1 non-outlier eigenvalues of K 0 that can be stable or unstable, the probability p = P(K * stable) is for N 1, using the binomial theorem and (1 + 1/x) x ∼ e x for x 1. A similar expression determines p = P(J * stable), with q replaced q . Eq. (6) shows that J 1 acts on K 0 as the negative definite matrix −κI, so q > q [5]. It follows that confirming the trend observed for small N in Fig. 2: the full model is much more likely to be stable than the averaged model. Hence stochastic switching to a rare phenotype such as the persister phenotype can enhance the stability of a community. Switching to an abundant phenotype destabilises the community, however, since, as argued earlier, the introduction of such a phenotype essentially increases the effective number of species in the system. The detailed analysis of Eq. (5) above has emphasised that this effect relies on both the spectral distribution (the structure of which allows overcoming an asymptotic separation of switching rates to stabilise or destabilise the system) and the detailed structure of the system (which can suppress this mechanism). Ref. [6] recently introduced a family of models with explicit resource dynamics for which any feasible equilibrium is stable. Switching to an abundant phenotype also destabilises a phenotype-resolving version of the model of Ref. [6] provided that the difference in switching rates is large enough [5], confirming that this effect is generic. These conclusions are therefore likely to be relevant for the stability of microbial communities such as the microbiome, for which competitive interactions are known to play an important, stabilising role [22].
Linear stability analysis cannot elucidate the effect of large perturbations on the coexistence state, however. Biologically, perturbations of this kind might correspond to antiobiotic treatments, which bacterial communities can survive by forming persister cells [12]. In the final part of this Letter, we therefore explore such perturbations numerically. Here, we do not  (5) and the corresponding averaged model as a function of N. Possible behaviours are convergence back to coexistence of all species, to coexistence of some, but not all species, or non-convergence (e.g. for cases in which a limit cycle arises). In the remaining cases (#), numerical solution failed. Probabilities were estimated from up to 10 4 random systems. (b) Distribution of the number n of remaining species for systems and initial conditions converging to a new coexistence state not involving all species, for N = 10 and for both the full and averaged models. (c) Probability that the full and averaged models converge to different coexistence states, for systems and initial conditions for which both the full and averaged models converge to a coexistence state. The contributions from systems for which one model converges to the full coexistence state of all species are highlighted. In all panels, error bars indicate 95% confidence intervals. model the dynamics of antibiotic treatments in detail [5]. Instead, we suppose that, at the end of such a treatment, bacterial and persister abundances are both small, and we ask: does the community evolve back into its coexistence state, or do some species disappear from the community? Do the answers differ between the averaged model and the full persister model? To answer these questions, we consider random systems for which the coexistence state is stable in both the averaged and full models, and compute numerical evolutions of both systems from consistent random small initial conditions using the stiff solver ode15s of M (The MathWorks, Inc.). The distribution of possible outcomes (convergence back to the coexistence state of all N species, convergence to a new coexistence state of n < N species, or absence of a steady state) is shown in Fig. 4(a). As N increases, the probability that some species are killed off increases [ Fig. 4(a)]. The distribution of the number n of species that survive the treatment if not all species do [ Fig. 4(b)] shows that n N/2 is somewhat more likely than n N/2, i.e. at least half the community tends to survive the antibiotic treatment if not the whole community does. (It is easy to see that the respective trivial equilibria A = 0 and B = C = 0 of the averaged and full models are unstable, and so 1 n < N.) The averaged and full persister models give comparable outcome distributions: this does not contradict our earlier result that a rare phenotype stabilises the community since here, we only consider states stable in both the averaged and full models so that results can be compared meaningfully. In fact, as N increases, individual realisations of full and averaged models are increasingly likely to give different outcomes [ Fig. 4(c)]. In most cases, different results arise because the system converges back to the full coexistence state of all species in one model only; in a small number of cases, different species die out in the full and averaged models [ Fig. 4(c)]. These observations thus extend our earlier results for linear perturbations to large perturbations.
In this Letter, we have revealed the strong effects of subpopulation structure on the stability of competing microbial communities and the surprising stabilising effect of a rare phenotype. While the competitive interactions considered here are important in systems such as the human microbiome [22], the type of interactions and the modularity of the community are known to affect its stability [15,23]. These studies, in the spirit of May's seminal work [14], are however based on the analysis of random Jacobians, assumed to correspond to feasible equilibria of otherwise unspecified dynamical systems. By contrast, here, we could not avoid specifying dynamical systems, since we had to establish a correspondence between phenotype-resolving and averaged models (and indeed our analysis has shown that the detailed structure of the models can matter for the questions considered here). In particular, this prepends the question of feasibility to the question of stability. Recent work [24] has shown that these questions can be disentangled in simple systems, but that analysis does not apply for the models considered here. While feasibility can be understood in a statistical sense [25], this does not easily yield an understanding of the feasibility of an individual system. For these reasons, we only studied competitive interactions, for which we could avail ourselves of the mathematical trick described in the first paragraph. Nonetheless, the fact that our explanation of the observed behaviour only relies on generic properties of the spectral distribution suggests that our results carry over to more general interactions.
The authors gratefully acknowledge support from the Engi-