Subpopulations and Stability in Microbial Communities

TL;DR

This paper investigates how phenotypic diversity within microbial species affects ecological stability, revealing that phenotypic variation can both destabilize and stabilize communities, challenging classical stability theories.

Contribution

It introduces a detailed analysis of phenotypic switching effects on stability, highlighting the stabilizing role of rare phenotypes like bacterial persisters.

Findings

Phenotypic variation can destabilize microbial communities.

Rare phenotypes such as persisters can stabilize communities.

Phenotypic diversity influences stability under large perturbations.

Abstract

In microbial communities, each species often has multiple, distinct phenotypes, but studies of ecological stability have largely ignored this subpopulation structure. Here, we show that such implicit averaging over phenotypes leads to incorrect linear stability results. We then analyze the effect of phenotypic switching in detail in an asymptotic limit and partly overturn classical stability paradigms: abundant phenotypic variation is linearly destabilizing but, surprisingly, a rare phenotype such as bacterial persisters has a stabilizing effect. Finally, we extend these results by showing how phenotypic variation modifies the stability of the system to large perturbations such as antibiotic treatments.

Equations

$$\displaystyle\boldsymbol{\dot{A}}=\boldsymbol{A}\bigl{(}\boldsymbol{\alpha}-\mathsf{D}\cdot\boldsymbol{A}\bigr{)},$$

$$\displaystyle\boldsymbol{\dot{B}}=\boldsymbol{B}\bigl{(}\boldsymbol{\beta}\!-\!\mathsf{E}\!\cdot\!\boldsymbol{B}\!-\!\mathsf{F}\!\cdot\!\boldsymbol{C}\bigr{)},\;\;\boldsymbol{\dot{C}}=\boldsymbol{C}\bigl{(}\boldsymbol{\gamma}\!-\!\mathsf{G}\!\cdot\!\boldsymbol{B}\!-\!\mathsf{H}\!\cdot\!\boldsymbol{C}\bigr{)},$$

$${\mathbf{A}}_{\ast }$$

$${\mathbf{A}}_{\ast }\mathsf{D}{\mathbf{A}}_{\ast }$$

$$\dot{\mathbf{B}}$$

$$\dot{\mathbf{C}}$$

$$\dot{\mathbf{B}}$$

$$\dot{\mathbf{C}}$$

$$\displaystyle\mathsf{J_{0}}=\left(\!\begin{array}[]{c|r}-\boldsymbol{B_{0}}\mathsf{E}&\boldsymbol{\lambda}\mathsf{I}\\ \hline\cr\mathsf{O}&-\boldsymbol{\lambda}\mathsf{I}\end{array}\!\right)\!,\;\mathsf{J_{1}}=\left(\!\begin{array}[]{r|c}-\boldsymbol{\kappa}\mathsf{I}&-\boldsymbol{B_{0}}\mathsf{F}\\ \hline\cr\boldsymbol{\kappa}\mathsf{I}&(\boldsymbol{\gamma}-\mathsf{G}\cdot\boldsymbol{B_{0}})\mathsf{I}\end{array}\!\right)\!,$$

$${\mathsf{K}}_0$$

$$p$$

$$={\left(\frac{1}{2}+\frac{c(2q-1)\epsilon }{2\sqrt{N}}\right)}^{N-1}\sim \frac{\exp \left[c(2q-1)\epsilon \sqrt{N}\right]}{2^{N-1}},$$

$$\frac{p^{\prime }}{p}\sim \exp \left[2c(q^{\prime }-q)\epsilon \sqrt{N}\right]\to \infty \text{as }N\to \infty ,$$

Full text

Over forty years ago, May suggested that equilibria of large ecological communities are overwhelmingly likely to be linearly unstable [1]. His approach did not specify the details of the dynamical system that describes the full population dynamics, but rather assumed that the linearized dynamics near the fixed point were represented by a random Jacobian matrix. Invoking results from random matrix theory, he concluded that unstable eigenvalues are more likely to arise as the number of interacting species increases. Actual large ecological communities certainly seem stable, and a major research theme in theoretical ecology has been to identify those features of the population dynamics that stabilize them [2, 3, 4, 5, 6, 7].

Recent advances in our understanding of large natural microbial communities such as the human microbiome have emphasized the important link between stability and function: adult individuals typically carry the same microbiome composition for long periods of time and disturbances thereof are often associated with disease [8, 9, 10]. Moreover, while genetically identical organisms may exhibit different phenotypes [11, 12, 13, 14] and despite the known ecological importance of phenotypic variation [15], studies of stability have largely ignored the existence of such subpopulations within species. Most models are therefore implicit averages over subpopulations.

We show here that this averaging yields incorrect stability results. With stochastic switching between phenotypes as an example of subpopulation structure, we show that while multiple abundant phenotypes are destabilizing, a rare phenotype can be stabilizing. This surprising result partly overturns May’s paradigm, and stresses the importance of phenotypic variation in ecological stability.

Our starting point is the Lotka–Volterra model [16], one of the most studied in population dynamics: $N$ species with abundances $\boldsymbol{A}$ compete as

[TABLE]

where $\boldsymbol{\alpha}$ is a vector of birth rates and $\mathsf{D}$ a matrix of non-negative competition strengths 111Throughout, we imply elementwise multiplication of vectors and (rows or columns of) matrices by writing symbols next to each other. Dots denote tensor contractions.; this is the competitive (as opposed to predator-prey) flavor of the model. If $\det{\mathsf{D}}\not=0$, Eq. (1) has a unique equilibrium $\boldsymbol{A_{\ast}}=\mathsf{D}^{-1}\cdot\boldsymbol{\alpha}$ of coexistence of all $N$ species. This equilibrium is feasible (i.e. $\boldsymbol{A_{\ast}}>\boldsymbol{0}$) if and only if $\boldsymbol{\alpha}$ lies in the positive span of the columns of $\mathsf{D}$.

If the $N$ species each have two subpopulations with respective abundances $\boldsymbol{B},\boldsymbol{C}$, then

[TABLE]

where $\boldsymbol{\beta},\boldsymbol{\gamma}$ are birth rates, and $\mathsf{E},\mathsf{F},\mathsf{G},\mathsf{H}$ are competition strengths [*[For$N=2$, systemsofthisformhavebeenconsideredinanalysesofevolutionarystability, seee.g.][]cressman03]. The dynamics of the sum $\boldsymbol{B}+\boldsymbol{C}$ derived from the subpopulation-resolving “full” system (2) are not of the averaged form (1). However, to a coexistence equilibrium $(\boldsymbol{B_{\ast}},\boldsymbol{C_{\ast}})$ of Eq. (2) we can associate an equilibrium $\boldsymbol{A_{\ast}}$ of Eq. (1), determined by the requirement that, at equilibrium, population sizes, births, and competition be equal [Fig. 1(a)], i.e. that

[TABLE]

These consistency conditions are a property of the model, based on the interpretation of its terms. Given Eqs. (2), they uniquely define an equilibrium $\boldsymbol{A_{\ast}}$ and an averaged model of the form (1), and this $\boldsymbol{A_{\ast}}$ is an equilibrium of this averaged model, and feasible if $(\boldsymbol{B_{\ast}},\boldsymbol{C_{\ast}})$ is.

We select random averaged and subpopulation-resolving systems by sampling model parameters from a uniform distribution 222See Supplemental Material at [url to be inserted], which includes Refs. [6, 20, 21, 22, 5, 23, *armstrong76, 25, 26], for (i) a detailed description of the selection of random systems, (ii) details of the asymptotic calculations leading to Eqs. (10) and (11), (iii) details of the spectral analysis leading to Eq. (13), including a note on the stability of sums of stable matrices, (iv) a discussion of models with explicit resource dynamics, based on those of Ref. [6], and (v) a brief discussion of the dynamics of antibiotic treatments., and analyze the stability of their coexistence equilibria by computing the eigenvalues of their Jacobians. As the number of species $N$ increases, stable equilibria of the averaged model (1) are increasingly likely to be unstable in the full model (2) [Fig. 1(b)]. This is because the full model effectively has $2N$ species, and stable equilibria become increasingly rare as the number of species increases [1, 3]. It is therefore all the more striking that, as $N$ increases, stable equilibria of the full model are also increasingly likely to be unstable in the averaged model [Fig. 1(b)]. The full model (2) can be extended to species with $M\geqslant 2$ subpopulations each, but increasing $M$ at fixed $N=2$ does not significantly affect the probability that a stable equilibrium of the full model destabilises in the averaged model [Fig. 1(c)].

These toy models thus show that implicit averaging of subpopulations leads to incorrect stability results, and hence underline their importance. Mathematically, this result is not fundamentally surprising: the determinant of the sum of two matrices is not the sum of their determinants, and so the linear relations between the Jacobians of the two systems resulting from Eqs. (3) cannot be expected to lead to simple relations between their stability.

We now specialize by taking phenotypic variation in microbial communities as an instance of subpopulation structure. It is useful to focus on one particular biological example: bacterial persisters [27, 25]. 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 or exposure to antibiotics [25, 28, 29]. Infections can thus be difficult to treat even in the absence of genetic antibiotic resistance; for this reason, this phenotype has great biomedical relevance [25, 28].

Adding switching between normal cells and persisters to Eqs. (2) leads to a phenotype-resolving model,

[TABLE]

where $\boldsymbol{\kappa},\boldsymbol{\lambda}$ are rates of (stochastic) switching. This form has previously been used to study phenotypic switching of a single species without competition [30, 26]. Since the rates are balanced, given an equilibrium of Eqs. (4), the consistency conditions (3) still define a corresponding “averaged” model (1) without phenotypic variation.

To analyze the effect of this switching, we again compare the stability of the full and averaged models. Steady states of Eqs. (4) cannot in general be found in closed form. To sample random systems, we must therefore sample parameters indirectly [19]. With increasing number of species, random stable coexistence states of either this full model or the corresponding averaged model again become increasingly likely to be unstable in the other model (Fig. 2).

Eqs. (4) do not however take into account the weak growth and competition of the persisters or the large separation of switching rates [30]. Adding a small parameter $\varepsilon\ll 1$, we therefore modify Eqs. (4) into

[TABLE]

Hence $\boldsymbol{B}$ are normal cells and $\boldsymbol{C}$ are persisters. For wild-type E. coli, $\varepsilon\approx 10^{-5}$ [30], but here, we take $\varepsilon=0.01$ for numerical convenience. We justify this below by confirming the numerical results by an asymptotic analysis of Eqs. (5). A more intricate asymptotic separation of parameters arises in the hipQ mutant of E. coli [30], but we do not pursue this further.

The separation of growth and competition terms and switching rates in Eqs. (5) allows for steady states to be found by expansion in $\varepsilon$. Writing \boldsymbol{B_{\ast}}=\boldsymbol{B_{0}}+\varepsilon\boldsymbol{B_{1}}+O\bigl{(}\varepsilon^{2}\bigr{)} and \boldsymbol{C_{\ast}}=\boldsymbol{C_{0}}+\varepsilon\boldsymbol{C_{1}}+O\bigl{(}\varepsilon^{2}\bigr{)}, we find $\boldsymbol{B_{0}}=\mathsf{E}^{-1}\cdot\boldsymbol{\beta}$ and $\boldsymbol{C_{1}}=\boldsymbol{\kappa B_{0}}/\boldsymbol{\lambda}$, but $\boldsymbol{B_{1}}=\boldsymbol{C_{0}}=\boldsymbol{0}$ [19]. This is the expected asymptotic separation of the two population sizes: few cells are persisters, at least under laboratory conditions [30]. This asymptotic solution enables direct sampling of all model parameters [19].

To analyze the stability of equilibria of Eqs. (5), we expand its Jacobian, \mathsf{J_{\ast}}=\mathsf{J_{0}}+\varepsilon\mathsf{J_{1}}+O\bigl{(}\varepsilon^{2}\bigr{)}, finding

[TABLE]

with $\mathsf{I}$ the identity and $\mathsf{O}$ the zero matrix [19]. The averaged model has Jacobian \mathsf{K_{\ast}}=\mathsf{K_{0}}+\varepsilon\mathsf{K_{1}}+O\bigl{(}\varepsilon^{2}\bigr{)}, with

[TABLE]

Since $\mathsf{J_{0}}$ is block-upper-triangular, its eigenvalues are those of $-\boldsymbol{\lambda}\mathsf{I}$, which are stable, and those of $-\boldsymbol{B_{0}}\mathsf{E}=\mathsf{K_{0}}$. Hence any unstable eigenvalues of $\mathsf{J_{\ast}}$ and $\mathsf{K_{\ast}}$ are equal to lowest order in the expansion. Equivalently, at $\varepsilon=0$, the full phenotype-resolving model (5) is stable if and only if the corresponding averaged model is stable. This result is not borne out however by numerics at finite $\varepsilon$: as $N$ increases, the probability that a random stable equilibrium of the full model is destabilized in the 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 destabilized in the full model is vastly reduced (Fig. 2, inset). This is all the more surprising as we argued earlier that the opposite behavior was to be expected since larger systems are more likely to be unstable. We have also sampled exact equilibria of Eqs. (5), similarly to our analysis of Eqs. (4) above, yielding results in qualitative agreement with those based on the asymptotic equilibria (Fig. 2, inset). This justifies basing the detailed analysis of the destabilization mechanism on the asymptotic results.

To explain these surprising results, we analyze the spectra of the Jacobians in more detail. With the exception of a single outlier eigenvalue that is large and negative, the eigenvalues of $\mathsf{K_{\ast}}$ lie approximately within a circle [Fig. 3(a)], as expected from the circular law of random matrices [31]. The spectral distribution of $\mathsf{J_{\ast}}$ is the sum of this distribution and the (uniform) distribution of eigenvalues of $-\boldsymbol{\lambda}\mathsf{I}$ [Fig. 3(b)]. The outlier eigenvalue can be analyzed in great generality [32], but heuristics suffice here: denoting by $-\mu_{0}<0$ the mean of the distribution of entries of $\mathsf{K_{0}}$ and neglecting correlations between entries, each row of $\mathsf{K_{0}}$ has approximate sum $-N\mu_{0}$, and so $\mathsf{K_{0}}$ has an approximate eigenvector $\mathbf{1}=(1,1,\dots,1)$ with eigenvalue $-N\mu_{0}<0$, as argued in the Supplemental Material of Ref. [3]. The other eigenvalues of $\mathsf{K_{0}}$ [Fig. 3(c)] are uniformly distributed on a disk of radius $r\sim\sqrt{N}$ for $N\gg 1$ by the circular law [31]. (Hence, by the Perron–Frobenius theorem [*[][Chap.1.3, pp.44--56, Chap.2.5, pp.100--112, andChap.8.2, pp.495--503.]linalg], the outlier eigenvalue is indeed real.) An eigenvalue $\nu_{0}$ of $\mathsf{K_{0}}$ has $|\mathrm{Re}(\nu_{0})|\lesssim\varepsilon$ with probability \varpi\sim(\varepsilon r)\big{/}r^{2}\sim\varepsilon\big{/}\sqrt{N}, i.e. \varpi=c\varepsilon\big{/}\sqrt{N}, for some $c=O(1)$. The average distance $\delta$ between eigenvalues is determined by $N\delta^{2}\sim r^{2}$, so $\delta=O(1)$ and the eigenvalues of $\mathsf{K_{0}}$ are pairwise different at order $O(1)$ [Fig. 3(c)]. Hence, if $\nu_{0}$ is an eigenvalue of $\mathsf{K_{0}}$, then $\mathsf{K_{\ast}}$ has an eigenvalue $\nu_{\ast}=\nu_{0}+O(\varepsilon)$ [19]. Thus $\nu_{\ast}$ is stable if either (i) $\mathrm{Re}(\nu_{0})<0$ and $|\mathrm{Re}(\nu_{0})|\gtrsim\varepsilon$ or (ii) $|\mathrm{Re}(\nu_{0})|\lesssim\varepsilon$ and the small real part of $\nu_{\ast}$ is stabilized by $\mathsf{K_{1}}$. By definition, $|\mathrm{Re}(\nu_{0})|\lesssim\varepsilon$ with probability $\varpi$, so (i) occurs with probability $(1-\varpi)/2$. Let $q$ denote the probability of stabilization by $\mathsf{K_{1}}$ in case (ii). Summing over the $N-1$ non-outlier eigenvalues of $\mathsf{K_{0}}$, the probability $p=\mathbb{P}(\mathsf{K_{\ast}}\,\text{stable})$ is

[TABLE]

for $N\gg 1$, using the binomial theorem and $(1+1/x)^{x}\sim\mathrm{e}^{x}$ for $x\gg 1$. A similar expression determines $p^{\prime}=\mathbb{P}(\mathsf{J_{\ast}}\,\text{stable})$, with $q$ replaced $q^{\prime}$. Eq. (10) shows that $\mathsf{J_{1}}$ acts on $\mathsf{K_{0}}$ as the negative definite matrix $-\boldsymbol{\kappa}\mathsf{I}$, so $q^{\prime}>q$ [19]. It follows that

[TABLE]

confirming the trend in Fig. 2: the full model is much more likely to be stable than the averaged model.

Hence, switching to a rare phenotype such as persisters can enhance the stability of a community. The detailed analysis of Eqs. (5) above has emphasized that this effect relies on both the spectral distribution (which allows terms beyond leading order to change the stability of the system) and the detailed structure of the system (which can suppress or enhance this mechanism). By contrast, switching to an abundant phenotype destabilizes the community: the introduction of such a phenotype essentially increases the effective number of species, which is destabilizing [1, 3]. Reference [6] recently introduced a family of models with explicit resource dynamics for which any feasible equilibrium is stable. Switching to an abundant phenotype also destabilizes a phenotype-resolving version of the model of Ref. [6] provided that the difference in switching rates is large enough [19], 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, stabilizing role [10].

Linear stability analysis cannot elucidate the effect of large perturbations on coexistence. These might arise from antibiotic treatments, which bacterial communities can survive by forming persisters [25]. In the final part of this Letter, we therefore explore such perturbations numerically. Rather than modelling the dynamics of antibiotic treatment in detail [19], we suppose that it reduces the abundances of both normal cells and persisters, and we ask: does the community converge back to coexistence, or do some species disappear from the community? Do the answers from the averaged and full models differ? To answer these questions, we reconsider the exact stable equilibria of the full model (5) that are also stable in the averaged model, and evolve both systems from consistent random small initial conditions using the stiff solver ode15s of Matlab (The MathWorks, Inc.).

Fig. 4(a) shows the distribution of possible outcomes: (1) convergence back to coexistence of all $N$ species, (2) convergence to a new coexistence state of $n<N$ species, and (3) convergence to a limit cycle. (The trivial equilibria $\boldsymbol{A}=\boldsymbol{0}$ and $\boldsymbol{B}=\boldsymbol{C}=\boldsymbol{0}$ of the averaged and full models are clearly unstable, so $1\leqslant n<N$.) The probability of outcome (2) increases with $N$ [Fig. 4(a)] and the distribution of $n$ in Fig. 4(b) shows that $n\gtrapprox N/2$ is somewhat more likely than $n\lessapprox N/2$. Thus if the whole community does not survive, then at least half does. The averaged and full persister models give comparable outcome distributions. This does not contradict our earlier result that a rare phenotype stabilizes 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 realizations of full and averaged models are increasingly likely to give different outcomes [Fig. 4(c)]. In most cases, outcomes differ because the system converges back to coexistence 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 results for linear perturbations to large perturbations.

Here, we have revealed the strong effects of subpopulation structure on the stability of competing microbial communities and the surprising stabilizing effect of stochastic switching to a rare phenotype. Very recently, Ref. [33] similarly emphasized the stabilizing effect of phenotypic variation using a different model, without phenotypic switching. While the competitive interactions considered here are important in systems such as the human microbiome [10], future work will need to explore the phenotypic structure in more detail. The interaction structure of ecological communities without phenotypic variation is known to affect their stability [3, 4], but these and related studies, in the spirit of May’s seminal work [1], are based on the analysis of random Jacobians. By contrast, here, we could not avoid specifying explicit dynamical systems, since we had to establish a correspondence between full and averaged models (and indeed our analysis has shown that the details of the model structure can matter here). In particular, this prepends the question of feasibility to the question of stability. This question of feasibility can be treated in a statistical sense [34, 35] but cannot be eschewed in general. Nonetheless, our analysis only relying on generic properties of the spectral distribution suggests that our conclusions apply not only to competitive, but also to more general interactions.