Abstract
Microbial communities routinely have several possible species compositions or community states observed for the same environmental parameters. Changes in these parameters can trigger abrupt and persistent transitions (regime shifts) between such community states. Yet little is known about the main determinants and mechanisms of multistability in microbial communities. Here, we introduce and study a consumerresource model in which microbes compete for two types of essential nutrients each represented by multiple different metabolites. We adapt gametheoretical methods of the stable matching problem to identify all possible species compositions of such microbial communities. We then classify them by their resilience against three types of perturbations: fluctuations in nutrient supply, invasions by new species, and small changes of abundances of existing ones. We observe multistability and explore an intricate network of regime shifts between stable states in our model. Our results suggest that multistability requires microbial species to have different stoichiometries of essential nutrients. We also find that a balanced nutrient supply promotes multistability and species diversity, yet make individual community states less stable.
eLife digest
In nature, different species of bacteria and fungi often live together in stable microbial communities. Exactly which species are present in the group and in which proportion may vary between communities. Changes in the environment, and in particular in the availability of nutrients, can trigger abrupt, extensive, and longlasting changes in the composition of a community: these events are known as regime shifts. For instance, when bodies of water receive large quantities of phosphorus and nitrogen, certain algae can start to multiply uncontrollably and take over other species. A given community can have different stable species compositions, but it was unclear exactly how variations in nutrients can influence regime shifts.
To examine this problem, Dubinkina, Fridman, Pandey and Maslov harnessed mathematical techniques used in game theory and economics and modeled all the possible stable compositions of a community. They could then predict which environmental conditions – in this case, the amount of specific nutrients – were necessary for each stable composition to exist. These models also showed which conditions could trigger a regime shift. Finally, how resilient the communities were to different types of perturbations – for instance, an invasion by new species or changes in nutrient supply – was examined.
The results show that if competing species require different quantities of the same nutrients, then the community can have several possible stable compositions and it is more likely to go through regime shifts. In addition, a small number of keystone species were identified which can drive regime shifts by preventing other microbes from invading the community. Ultimately, these results suggest ways to control microbial communities in our environment, for example by manipulating nutrient supplies or introducing certain species at the right time. More work is needed however to verify the predictions of the model in real communities of microbes.
Introduction
Recent metagenomics studies revealed that microbial communities living in similar environments are often composed of rather different sets of species (Zhou et al., 2007; Lahti et al., 2014; Lozupone et al., 2012; Zhou et al., 2013; Pagaling et al., 2017; Gonze et al., 2017). It remains unclear to what extent such alternative species compositions are deterministic as opposed to being an unpredictable outcome of communities’ stochastic assembly. Furthermore, changes in environmental parameters may trigger abrupt and persistent transitions between alternative species compositions (Shade et al., 2012; Rocha et al., 2018; Scheffer and Carpenter, 2003). Such transitions, known as ecosystem regime shifts, significantly alter the function of a microbial community and are difficult to reverse. Understanding the mechanisms and principal determinants of alternative species compositions and regime shifts is practically important. Thus, they have been extensively studied over the past several decades (Sutherland, 1974; Holling, 1973; May, 1977; Tilman et al., 1997; Schröder et al., 2005; Fukami and Nakajima, 2011; Bush et al., 2017).
Growth of microbial species is affected by many factors, with availability of nutrients being among the most important ones. Thus, the supply of nutrients and competition for them plays a crucial role in determining the species composition of a microbial community. The majority of modeling approaches explicitly taking nutrients into account are based on the classic MacArthur consumerresource model and its variants (MacArthur and Levins, 1964; MacArthur, 1970; Huisman and Weissing, 1999; Tikhonov and Monasson, 2017; Posfai et al., 2017; Goldford et al., 2018; Goyal et al., 2018; Butler and O'Dwyer, 2018). This model assumes that every species coutilizes several substitutable nutrients of a single type (e.g. carbon sources). However, nutrients required for growth of a species exist in the form of several essential (nonsubstitutable) types including sources of C, N, P, Fe, etc. Real ecosystems driven by competition for multiple essential nutrients have been extensively experimentally studied (see recent papers; Fanin et al., 2015; Browning et al., 2017; Camenzind et al., 2018 and references therein). The theoretical foundation for all existing consumerresource models capturing this type of growth has been laid in Tilman (1982), where a model with two essential resources has been introduced and studied. Future studies extended Tilman’s approach to three and more essential resources, where it has been shown to sometimes result in oscillations and chaos (Huisman and Weissing, 1999; Huisman and Weissing, 2001; Shoresh et al., 2008). However, all the previously studied models accounted for just a single metabolite per each essential nutrient.
Here, we introduce and study a new consumerresource model of a microbial community supplied with multiple metabolites of two essential types (e.g. C and N or N and P). This ecosystem is populated by microbes selected from a fixed pool of species. We show that our model has a very large number of possible steady states classified by their distinct species compositions. Using gametheoretical methods adapted from the wellknown stable marriage (or stable matching) problem (Gale and Shapley, 1962; Gusfield and Irving, 1989), we predict all these states based only on the ranked lists of competitive abilities of individual species for each of the nutrients. We further classify these states by their dynamic stability, and whether they could be invaded by other species in our pool. We then focus our attention on a set of steady states that are both dynamically stable and resilient with respect to species invasion.
For each state, we identify its feasibility range of all possible environmental parameters (nutrient supply rates) for which all of state’s species are able to survive. We further demonstrate that for a given set of nutrient supply rates, more than one state could be simultaneously feasible, thereby allowing for multistability. While the overall number of stable states in our model is exponentially large, only very few of them can be realized for a given set of environmental conditions defined by nutrient supply rates. The principal component analysis of predicted microbial abundances in our model shows a separation between the alternative stable states reminiscent of reallife microbial ecosystems. We further explore an intricate network of regime shifts between the alternative stable states in our model triggered by changes in nutrient supply. Our results suggest that multistability requires microbial species to have different stoichiometries of two essential resources. We also find that wellbalanced nutrient supply rates matching the average species’ stoichiometry promote multistability and species diversity yet make individual community states less structurally and dynamically stable. These and other insights from our consumerresource model may help to understand the existing data and provide guidance for future experimental studies of alternative stable states and regime shifts in microbial communities.
Results
Microbial community growing on two types of essential nutrients represented by multiple metabolites
Our consumerresource model describes a microbial ecosystem colonized by microbes selected from a pool of $S$ species. Growth of each of these species could be limited by two types of essential resources, to which we refer to as ‘carbon’ and ‘nitrogen’. In principle, these could be any pair of resources essential for life: C, N, P, Fe, etc. A generalization of this model to more than two types of essential resources (e.g. C, N and P) is straightforward. Carbon and nitrogen resources exist in the environment in the form of $K$ distinct metabolites containing carbon , and $M$ other metabolites containing nitrogen. For simplicity, we ignore the possibility of the same metabolite providing both types. We further assume that each of the $S$ species in the pool is a specialist, capable of utilizing only a single pair of nutrients, that is one metabolite containing carbon and one metabolite containing nitrogen.
We assume that for given environmental concentrations of all nutrients, a growth rate of a species $\alpha $ is limited by a single essential resource via Liebig’s law of the minimum (de Baar, 1994):
Here, $c$ and $n$ are the environmental concentrations of the unique carbon and nitrogen resources consumed by this species. The coefficients ${\lambda}_{\alpha}^{(c)}$ and ${\lambda}_{\alpha}^{(n)}$ are defined as speciesspecific growth rates per unit of concentration of each of two resources. They quantify the competitive abilities of the species $\alpha $ for its carbon and nitrogen resources, respectively. Indeed, according to the competitive exclusion principle, if two species are limited by the same resource, the one with the larger value of $\lambda $ wins the competition. Note that according to Liebig’s law, if the carbon source is in short supply so that ${\lambda}_{\alpha}^{(c)}c<{\lambda}_{\alpha}^{(n)}n$, it sets the value for this species growth rate. We refer to this situation as $c$source limiting the growth of the species $\alpha $. Conversely, when ${\lambda}_{\alpha}^{(c)}c>{\lambda}_{\alpha}^{(n)}n$, the $n$source is limiting the growth of this species. Thus, each species always has exactly one growthlimiting resource and one nonlimiting resource.
In our model, microbes grow in a wellmixed chemostatlike environment subject to a constant dilution rate $\delta $ (see Figure 1A for an illustration). The dynamics of the population density, ${B}_{\alpha}$, of a microbial species $\alpha $ is then governed by:
where ${c}_{i}$ and ${n}_{j}$ are the specific pair of nutrients defining the growth rate ${g}_{\alpha}$ of this species according to the Liebig’s law (Equation 1). These nutrients are externally supplied at fixed rates ${\varphi}_{i}^{(c)}$ and ${\varphi}_{j}^{(n)}$ and their concentrations follow the equations:
Here, ${Y}_{\alpha}^{(c)}$ and ${Y}_{\alpha}^{(n)}$ are the growth yields of the species $\alpha $ on its $c$ and $n$resources respectively. Yields quantify the concentration of microbial cells generated per unit of concentration of each of these two consumed resources. The yield ratio ${Y}_{\alpha}^{(n)}/{Y}_{\alpha}^{(c)}$ determines the unique C:N stoichiometry of each species.
A steady state of the microbial ecosystem can be found by setting the right hand sides of Equations 23 to zero and solving them for environmental concentrations of all nutrients ${c}_{i}$, and ${n}_{j}$, and abundances ${B}_{\alpha}$ of all species. We choose to label all possible steady states by the list of species present in the state and by the growthlimiting nutrient ($c$or $n$) for each of these species. Thus, two identical sets of species, where at least one species is growth limited by a different nutrient are treated as two distinct states of our model. Conversely, our definition of a steady state does not take into account species’ abundances. Examples of such states in a system with two carbon, two nitrogen nutrients and four species (one species for every pair of carbon and nitrogen nutrients) with specific values of species’ competitive abilities ${\lambda}_{\alpha}^{(c)}$ and ${\lambda}_{\alpha}^{(n)}$ and yields ${Y}_{\alpha}^{(c)}$ and ${Y}_{\alpha}^{(n)}$ (see Supplementary files 1,2 for their exact values) are shown in Figure 1B. For the sake of brevity we refer to this model as 2C × 2N × 4S.
Because each of the $S$ species in the pool could be absent from a given state, or, if present, could be limited by either its $c$ or its $n$resource, the theoretical maximum of the number of distinct states is 3^{S} (equal to 81 in our 2C × 2N × 4S example). However, the actual number of possible steady states is considerably smaller (equal to 34 in this case). Indeed, possible steady states in our model are constrained by a variant of the competitive exclusion principle (Gause, 1932) (see Materials and methods for details). One of the universal consequences of this principle is that the number of species present in a steady state of any consumerresource model cannot exceed K + M − the total number of nutrients.
We greatly simplified the task of finding all steady states in our model by the discovery of the exact correspondence between our system and a variant of the celebrated stable matching (or stable marriage) problem in game theory and economics (Gale and Shapley, 1962; Gusfield and Irving, 1989). The matching in our model connect pairs of C and N resources via microbial species using both of them. Unlike in the traditional stable marriage model, a given resource can be involved in more than one matching but cannot be limiting for more than one microbe. Thus, the competitive exclusion principle provides a number of constraints on the set of possible matchings and their stability, which are described in detail in Materials and methods and Appendix 3.
Three criteria for stability of microbial communities
Each of the steady states identified in the previous chapter can be realized only for a certain range of nutrient supply rates. These ranges can be calculated using the steady state solutions of Equations 2, 3, governing the dynamics of microbial populations and nutrient concentrations, respectively (see Materials and methods). Among all formal mathematical solutions of these equations we select those, where populations of all species and all nutrient concentrations are nonnegative. This imposes constraints on nutrient supply rates, thereby determining their feasible range for a given steady state (shown in green in Figure 1C). The volume of such feasible range has been previously used to quantify the socalled structural stability of a steady state (Rohr et al., 2014; Grilli et al., 2017; Butler and O'Dwyer, 2018). States with larger feasible volumes generally tend to be more resilient with respect to fluctuations in nutrient supply.
Stability of a steady state could be also disturbed by a successful invasion of a new species (see Figure 1D). We can test the resilience of a given state in our model with respect to such invasions. A state is called uninvadable if none of the other species from our pool can survive in the environment shaped by the existing species. Figure 1B shows all seven states that are uninvadable in our variant of the 2C × 2N × 4S model. Whether or not a given state is uninvadable is determined by the specific choice of parameters ${\lambda}_{\alpha}^{(c)}$, ${\lambda}_{\alpha}^{(n)}$. For example, for parameters listed in the Supplementary file 1 the state in which ${B}_{12}$ is limited by carbon ${c}_{1}$, and ${B}_{22}$  by carbon ${c}_{2}$ could be invaded by the species ${B}_{11}$. Indeed, ${\lambda}^{(c)}$ of ${B}_{11}$ is larger than that of ${B}_{12}$, and the nitrogen concentration ${n}_{1}$ is not limited by any species. Hence, this state is not shown in Figure 1B. However, the same state may turn out to be uninvadable for a different combination of parameters. The onetoone correspondence between our model and a variant of the stable matching problem (Gale and Shapley, 1962) allows us to identify all uninvadable steady states for a given choice of ${\lambda}_{\alpha}^{(c)}$, ${\lambda}_{\alpha}^{(n)}$ describing species competitiveness for resources (see Materials and methods and Appendix 3 ).
Note that, in the regime of our model, where the supply of all nutrients is high, that is ${\varphi}^{(c,n)}\gg {\delta}^{2}/{\lambda}_{\alpha}^{(c,n)}$, invadability of individual states does not depend on supply rates. Indeed, in this regime the outcome of an attempted invasion is fully determined by the competition between species, which in turn depends only on the rankorder of competitive abilities $\lambda $ of the invading species relative to the species currently present in the ecosystem (see Materials and methods for details).
In addition to structural and invasion types of stability described above, there is also a notion of dynamic stability of a steady state actively discussed in the ecosystems literature (see e.g. May, 1972; Allesina and Tang, 2012; Butler and O'Dwyer, 2018). Dynamic stability can be tested by exposing a steady state to small perturbations in populations of all species present in this state (see Figure 1E). The state is declared dynamically stable if after any such disturbance the system ultimately returns to its initial configuration (see Materials and methods for details of the testing procedure used in our study).
We classify all the steady states in our model according to these three types of stability. The example of this classification for our realization of 2C × 2N × 4S model is summarized in Figure 1F. Note, that in general, one type of stability does not imply another. Out of 34 possible steady states realized for different ranges of nutrient supply rates there are only seven uninvadable ones. In the 2C × 2N × 4S model only one of the states (labelled seven in Figure 1B) turned out to be dynamically unstable, while for the remaining 33 states small perturbations of microbial abundances present in the state do not trigger a change of the state. Unlike two other types of stability, the structural stability has a continuous range. It could be quantified by the fraction of all possible combinations of nutrient supply rates for which a given state is feasible (referred to as state’s normalized feasible range). We estimated the normalized feasible ranges of all states in the 2C × 2N × 4S model using a Monte Carlo procedure described in Materials and methods. The results are reflected in the second column of Figure 1F, where a structurally stable state is defined as that whose normalized feasible range exceeds 0.1 (an arbitrary threshold). In general we find that normalized feasible ranges of uninvadable states in our model have a broad lognormal distribution (see Figure 1—figure supplement 1 for details).
It is natural to focus our attention on steady states that are simultaneously uninvadable and dynamically stable. Indeed, such states correspond to natural endpoints of the microbial community assembly process. They would persist for as long as the nutrient supply rates do not change outside of their structural stability range. Therefore, they represent the states of microbial ecosystems that are likely to be experimentally observed. From now on, we concentrate our study almost exclusively on those states and refer to them simply as stable states.
Regime shifts between alternative stable states
The feasible ranges of nutrient supply of different stable states may or may not overlap with each other (see Figure 2A–B for a schematic illustration of two different scenarios). Whenever feasible ranges of two or more states overlap (see Figure 2B)  multistability ensues. Note that the states in the overlapping region of their feasibility ranges constitute true alternative stable states defined and studied in the ecosystems literature (Sutherland, 1974; Holling, 1973; May, 1977; Fukami and Nakajima, 2011; Bush et al., 2017). The existence of alternative stable states goes handinhand with regime shifts manifesting themselves as large discontinuous and hysteretic changes of species abundances (Scheffer and Carpenter, 2003).
Every pair of states with overlapping feasibility ranges in our model corresponds to a possible regime shift between these states as illustrated in Figure 2D (note abrupt changes in the population ${B}_{11}$ at the boundary of the overlapping region). In general, a discontinuous regime shift happens in our model when one of the species (${B}_{12}$ in this example) changes its growthlimiting nutrient thereby making the state invadable. It is then promptly invaded by the species present in the new state (${B}_{11}$ and ${B}_{22}$ in our example) which may lead to immediate changes in populations of multiple species. Conversely, when feasible ranges of a pair of states do not overlap with each other but share a boundary (Figure 2A), the transition between these states is smooth and nonhysteretic (Figure 2C). It manifests itself in continuous changes in abundances of all microbial species at the boundary between states. Such continuous transitions happen in our model when the growth rate of one of the species (${B}_{21}$ in this example) falls below the dilution rate $\delta $. This species then slowly disappears from the ecosystem thereby changing its state. When this boundary is crossed in the opposite direction, the same species (${B}_{21}$) gradually appears in the ecosystem.
As expected for regime shifts, dynamically unstable states always accompany multistable regions in our model (Scheffer and Carpenter, 2003) (see below for a detailed discussion of the interplay between multistability and dynamically unstable states). We observed that dynamically unstable state #7 in our 2C × 2N × 4S is feasible in the overlapping region between states #1 and #2 in Figure 2B. The population ${B}_{11}$ in this state is shown as dashed line in Figure 2D.
We identified all possible regime shifts in the 2C × 2N × 4S model by systematically looking for overlaps between the feasible ranges of nutrient supply of all six uninvadable dynamically stable states. These regime shifts can be represented as a network in which nodes correspond to community’s stable states and edges connect states with partially overlapping feasible ranges (see thick red edges in Figure 2E). One can see that regime shifts are possible only for of nine pairs of uninvadable states. We performed additional simulations (see Materials and mthods) looking for shared boundaries (continuous transitions) between uninvadable states and identified additional four pairs of states bordering each other (thin blue edges in Figure 2E). The pairs of states #5  #6 and #3  #4 do not directly transition to each other either continuously or discontinuously. This indicates that their feasible ranges are too far apart from each other, so that they do not have any overlaps or common boundaries.
Combining the information in Figure 1B and Figure 2E one can find that all states connected by a discontinuous regime shift in our 2C × 2N × 4S model have two distinct sets of keystone species: ${B}_{11}$${B}_{22}$ in one state and ${B}_{12}$${B}_{21}$ in another. This is because all regime shifts are driven by the same bistable switch in which these pairs of species compete and mutually exclude each other. The dynamically unstable state #7 is formed by the union of all four keystone species and, when perturbed, collapses into a state with either one or another keystone set. Conversely, states connected by a continuous transition share the same pair of keystone species. One of the ‘satellite’ species, that is species distinct from the keystone, gradually goes extinct when the boundary between these states is crossed. When the nutrient supply is changed in the opposite direction this species gradually invades the system.
Figure 2F shows a much larger network of 8633 regime shifts between 893 uninvadable dynamically stable states in the 6C × 6N × 36S realization of our model. In this model the microbial community is supplied with six carbon and six nitrogen nutrients and colonized from a pool of 36 microbial species (one for each pair of C and N nutrients) (see Supplementary files 3, 4, 5, 6 for the values of $\lambda $’s and yields). For simplicity, we did not show the remaining 165 uninvadable stable states that have no possible regimes shifts to any other states. The size of a node is proportional to its degree (i.e. the total number of other states it overlaps with) ranging between 1 and 164 with average around 20 (the degree distribution is shown in Figure 2—figure supplement 1).
The network modularity analysis (see Materials and methods for details) revealed seven network modules indicating that pairs of states that could possibly undergo a regime shift are clustered together in the multidimensional space of nutrient supply rates. This modular structure suggests the existence of distinct sets of keystone species driving regime shifts within each module. However, the complexity of the 6C × 6N × 36S model does not allow a straightforward identification of these drivers (paired sets of keystone species and dynamically unstable states).
Patterns of multistability
In a general case, the number of stable states that are simultaneously feasible for a given set of nutrient supply rates can be more than two. Furthermore, as the number of nutrients increases, the multistability with more than two stable states becomes progressively more common. In Figure 3A, we quantify the frequency of multistability with $V$ stable states occur in our 6C × 6N × 36S model across all possible nutrient supply rates (see Materials and methods for details of how this was estimated). $V1$ approximately follows a Poisson distribution (dashed line in Figure 3A) with $\lambda =0.063$. Note that for some supply rates up to five stable states can be simultaneously feasible. However, the probability to encounter such cases is exponentially small.
We further explored the factors that determine whether multistability is possible in resourcelimited microbial communities. Like in a simple special case of regime shift between two microbial species studied in Tilman (1982), multistability in our model is only possible if individual microbial species have different C:N stoichiometry. This stoichiometry is given by the ratio of species’ nitrogen and carbon yields. Our numerical simulations and mathematical arguments show that when all species have exactly the same stoichiometry ${Y}_{\alpha}^{(n)}/{Y}_{\alpha}^{(c)}$, there is no multistability or dynamical instability in our model (see Appendix 5 ) . That is to say, in this case for every set of nutrient supply rates the community has a unique uninvadable state, and all these states are dynamically stable.
A complementary question is whether multistable states are more common around particular ratios of carbon and nitrogen supply rates. Figure 3B shows this to be the case: the likelihood of multistability has a sharp peak around the wellbalanced C:N nutrient supply rates. In this region multiple stable states are present for roughly 15% of nutrient supply rate combinations. Note that the average C:N stoichiometry of species in our model is assumed to be 1:1. In case of an arbitrary C:N stoichiometry, by redefining the units of nutrient concentrations and supply rates one can transform any ecosystem to have 1:1 nutrient ratio. Hence, in general, we predict that the highest chance to observe multistability will be when the ratio of nutrient supply rates is close to the average C:N stoichiometry of species in the community.
To illustrate how multistable states manifest themselves in a commonly performed Principal Component Analysis (PCA) of species’ relative abundances, we picked the environment with $V=5$ simultaneously feasible stable states in our 6C × 6N × 36S model. In natural environments, nutrient supply usually fluctuates both in time and space. To simulate this we sampled a ±10% range of nutrient supply rates around this chosen environment (see Materials and methods) and calculated species’ relative abundances in each of the uninvadable states feasible for a given nutrient supply. To better understand the relationship between dynamically stable and unstable states we included the latter in our analysis. Figure 3C shows the first vs the second principal components of relative microbial abundances sampled in this fluctuating environment. (two more examples calculated for different multistable neighborhoods are shown in Figure 3—figure supplement 1A–B). One can see five distinct clusters, each corresponding to a single dynamically stable uninvadable state. Interestingly, in the PCA plot these states are separated by $V1=4$ dynamically unstable ones. Furthermore, all states are aligned along a quasi1D manifold with an alternating order of stable and unstable states. It is tempting to conjecture that some variant of our model may explain similar arrangements of clusters of microbial abundances, commonly seen in PCA plots of real ecosystems. If this is the case, the gaps between neighboring clusters would correspond to dynamically unstable states of the ecosystem, which may be experimentally observable as long transients in community composition.
Patterns of diversity and structural stability of states
Above we demonstrated that multistable states are much more common for balanced nutrient supply rates, that is to say, when the average ratio of carbon and nitrogen supply rates matches the average C:N stoichiometry of species in the community (see Figure 3B). Interestingly, a balanced supply of nutrients also promotes species diversity. In Figure 4A, we plot the average number of species in a stable state, referred to as species richness, as a function of the average balance between carbon and nitrogen supplies for 6C × 6N × 36S model. The species richness is the largest (around 10.5) for balanced nutrient supply rates, while dropping down to the absolute minimal value of six in two extreme cases of very large imbalance of supply rates, where the nutrient supplied in excess becomes irrelevant in competition. In this case, only six species that are teh top competitors for carbon metabolites (if nitrogen supply is plentiful) or, respectively nitrogen metabolites (if carbon is large) survive, while the rest of less competitive species are never present in uninvadable states. The number of distinct community states also has a sharp peak at balanced nutrient supply (see 3orders of magnitude difference in Figure 4—figure supplement 1).
For balanced nutrient supply rates the relationship between species’ competitiveness and its prevalence in the community is much less pronounced than for imbalanced ones. It is shown in Figure 4B, where we plot the prevalence of the species as a function of its average competitiveness. Here, the average competitiveness rank of a species is defined as the mean of its ranks of competitive abilities ($\lambda $ parameters of the model) for its carbon and nitrogen resources. The rank 1 being assigned to the most competitive species for a given resource (the species with the largest value of $\lambda $), while the rank 6  to the least competitive species for this resource. Species prevalence is given by the fraction of all environments where it can survive. Note that all 36 species in our pool are present in some of the environments.
In general, more competitive species tend to survive in a larger subset of environments (see the dashed curve in Figure 4B). For example, in our pool there is one species which happens to be the most competitive for both its carbon and nitrogen sources. This species is present in all of the states in every environment. However, we also find that some of the least competitive species (those at the right end of the xaxis in Figure 4B) survive in a broad range of environments. For example, one species with average competitiveness rank of 5.5 corresponding to the last and next to last rank for its two resources still has relatively high prevalence of around 20%. This illustrates complex ways in which relative competitiveness of all species in the pool shapes their prevalence in a broad range of environments.
We also explore the relationship between species richness of a state (i.e. its total number of surviving species) and its other properties. Figure 4C shows an exponential increase of the number of uninvadable states as a function of species richness. In our 6C × 6N × 36S model all uninvadable states with less than 10 species are dynamically stable (solid line in Figure 4C), while those with 10 or more species can be both stable or unstable (dashed line in Figure 4C). Overall, the fraction of stable states to dynamically unstable ones decreases with species richness. In other words, the probability for a state to be dynamically unstable increases with the number of species. In this aspect, our model behaves similar to the gLV model in Robert May’s study (May, 1972).
In Figure 4D, we show a negative correlation between the species richness of a stable state and its feasible range of nutrient supplies. Thus in our model the number of species in an ecosystem has detrimental effect on the structural stability of the community quantifying its robustness to fluctuating nutrient supply (Rohr et al., 2014). The empirically observed exponential decay of state’s feasible range with its number of species is well described by a twofold decrease per each species added (see Serván et al., 2018 and Grilli et al., 2017 for related results in the gLV model). Note that the observed decrease in feasible range with species richness goes handinhand with an increase in the overall number of states. Thus, in wellbalanced environments a large number of states are carving all possible combinations of nutrient supply into many small and overlapping ranges.
Overall, the results of our model with a large number of nutrients suggest the following picture. In nutrientbalanced environments, we expect to observe a high diversity of species in the existing communities. These species can form a very large number of possible combinations (uninvadable states). Each of these states could be realized only for a narrow range of nutrient supply rates indicating their low structural stability. Moreover, in such environments we predict common appearance of multistability between some of these states.
Discussion
The inspiration for our model was the common appearance of alternative stable states in ecosystems in general, and microbial communities in particular (Sutherland, 1974; Tilman et al., 1997; Schröder et al., 2005; Fukami and Nakajima, 2011; Bush et al., 2017; Pagaling et al., 2017; Gonze et al., 2017). For example, eutrophication of shallow lakes caused by algal competition for N and P is one of the best studied examples of alternative stable states and regime shifts (Scheffer and Jeppesen, 2007). To the best of our knowledge our model is the first consumerresource model capable of multistability between several states, each characterized by a high diversity of species. We extend Tilman’s scenario (Tilman, 1982) in which the growth of two species is limited by a pair of essential resources to the case of multiple nutrients of each type. This allows us to assemble complex communities with large number of coexisting species and provides additional insights into patterns of multistability in such communities.
Multistability requires diverse species stoichiometry
We find that multistability in our model requires a mix of species with different nutrient stoichiometries. In this aspect it is similar to both the Tilman model (Tilman, 1982), and the MacArthur model (MacArthur and Levins, 1964; MacArthur, 1970; Chesson, 1990). Common variants of the MacArthur model assume identical biomass yields of different species growing on a given nutrient (Tikhonov and Monasson, 2017; Posfai et al., 2017; Goldford et al., 2018; Goyal et al., 2018; Butler and O'Dwyer, 2018). In this case, the absence of multistability is guaranteed by a convex Lyapunov function (MacArthur, 1970) guiding any dynamical trajectory of the system to its unique minimum. However, a MacArthur model with different nutrient yields of different species is capable of multistability. For some growth functions $g(C,N)$ multistability is possible even in a community of species with identical nutrient yields/stoichiometries. For example, the growth function $g(C,N)=\lambda \cdot C\cdot N$ has been numerically studied in the context of autocatalytic polymer growth and shown to be capable of multistability (Tkachenko and Maslov, 2018). This model had 1:1 stoichiometry: a ligation always eliminates one left end and one right end of two polymer chains and generates one autocatalytic polymer segment. Another type of growth function with two essential resources has been shown to have bistable solutions even for identical species stoichiometries (see Figure 5C in Marsland et al., 2019b). The Minimum Environmental Perturbation Principle introduced in this study may provide additional insights on the necessary conditions for multistability in consumer resource models.
Given that multistability in our model is impossible in communities of species with identical stoichiometries, it is reasonable to expect that the larger is the variation of C:N ratio of individual microbes, the higher is the likelihood to observe multistability. We investigated this question in the 2C × 2N × 4S model and summarized the results in Figure 5. It shows that the likelihood of finding nutrient supply rates with multistability systematically decreases with standard deviation of the logarithm of species stoichiometry. We also found that about half of the combinations of species stoichiometries yielded no multistable states at all. Multistability in our model is caused by a complex interplay between species’ competitiveness abilities $\lambda $ and their C:N stoichiometries $\frac{{Y}^{(n)}}{{Y}^{(c)}}$. Appendix 4 explains why multistability is impossible for half of yield combinations for which the enumerator in Equation S11 exceeds its denominator.
Nutrient stoichiometry of phytoplankton species in marine ecosystems has been known to be relatively universal with C:N:P ≃ 106:16:1 known as Redfield ratio (Redfield, 1958). Thus speciestospecies variability of C:N ratio for phytoplankton is rather small with logarithmic standard deviation estimated to be around 0.05 based on data from Finkel et al. (2016). In this limit, the multistability in our model is rather unlikely (observed in ∼3% of nutrient supply combinations, see green arrow in Figure 5). The likelihood of multistability is also low (∼1%) for the mammalian gut microbiome, where variability of the logarithm of C:N ratio in different 'keystone’ gut species studied in Reese et al. (2018) is around 0.03. The chances of multistability increase in terrestrial ecosystems such as soil, where significant deviations from the Redfield ratio have been reported (Cleveland and Liptzin, 2007). For example, using the data for the microbial species from grassland leaf litter community reported in Mouginot et al. (2014), with log(C:N) variability of 0.12 we predict the likelihood of multistability to be around 10%.
Multistability requires balanced nutrient supply matching the average species stoichiometry
Another important factor favoring multistability in our model is the balanced supply of two essential nutrients (see Figure 3B). It occurs when the average ratio of supply rates of two essential nutrients matches the average C:N stoichiometry of community’s species (see Figure 3B). When nutrient supplies are balanced, microbial community multistability is relatively common. Furthermore, for balanced nutrients the community can be in one of many different states, characterized by different combinations of limiting nutrients. These states tend to have high species diversity (Figure 4A) – a trend consistent with lake ecosystems in Interlandi and Kilham (2001) – and relatively small range of feasible supply rates (Figure 4D). Hence, regime shifts can be easily triggered by changes in nutrient supply. The balanced region is characterized by a complex relationship between species competitiveness and survival, so that even relatively poor competitors could occasionally have high prevalence (species in the upper right corner of Figure 4B).
In the opposite limit, the supply of nutrients of one type (say nitrogen) greatly exceeds that of another type (say carbon). For such imbalanced supply, the community has a unique uninvadable state, where every carbon nutrient supports the growth of the single most competitive species. Nitrogen nutrients are not limiting the growth of any species and thus have no impact on species survival and community diversity. As a consequence, the average diversity of microbial communities in such nutrientimbalanced environments is low (about one half of that for balanced supply conditions). This is in agreement with many experimental studies showing that addition of high quantities of one essential nutrient (e.g. as nitrogen fertilizer) tends to decrease species diversity. This has been reported in numerous experimental studies cited in the chapter 'Resource richness and species diversity’ of Tilman (1982) as well as in recent experiments in microbial communities (Mello et al., 2016).
Multistability and the total number of states are affected by tradeoffs
Species in our model are characterized by their competitiveness abilities ${\lambda}^{(c)}$, ${\lambda}^{(n)}$ and nutrient yields ${Y}^{(c)}$ and ${Y}^{(n)}$. As we showed above, the rank order of the former fully defines the total number of stable states and their invadability. On the other hand, multistability highly depends on combination of species’ nutrient yields. While in the current version of our model we did not assume any specific correlations between these parameters, imposing such correlations due to various biologically motivated tradeoffs may affect multistability and the total number of states of the ecosystem.
One possibility is a negative correlation between the competitive abilities of a given species for different nutrients. Such tradeoff may exist due to a limited amount of internal resources (such as the overall number of transporters) this species can allocate for consumption of all nutrients. This type of tradeoff was shown to result in an increased species diversity in wellbalanced environments, but does not lead to multistability (Posfai et al., 2017; Tikhonov, 2016). Similar negative (positive) correlations to increase (decrease) the number of stable states in a very different consumerresource model based on the stable marriage problem (Goyal et al., 2018). We expect these results to also apply to our model with tradeoff between ${\lambda}_{\alpha}^{(c)}$ and ${\lambda}_{\alpha}^{(n)}$. Negative correlations between species’ competitive abilities for carbon and nitrogen are expected to increase the total number of stable states in our model, while positive correlations  to decrease it.
Another possibility is a negative correlation between species’ competitive ability and its yield for the same nutrient. It is known as a ‘growthyield tradeoff’, which states that microbial species with faster growth on a given nutrient tend to use it less efficiently (have a smaller yield) (Pfeiffer et al., 2001; Beardmore et al., 2011; Novak et al., 2006). Growthyield tradeoff is expected to increase the likelihood of multistability in our model. It could be demonstrated already in the model of Tilman (1982) with two species competing for two essential resources. If the species with the higher growth rate on, say, carbon source has a smaller yield on this resource than the other species  bistability always ensues. Note that, while growthyield tradeoff is known to be common among microorganisms, the macroscopic (e.g. plant) ecosystems, which are the main focus of Tilman (1982), have the opposite correlation in which species’ yield $Y$ is proportional to its competitive ability $\lambda $. This type of tradeoff leads to a relative scarcity of multistability in macroscopic ecosystems. Conversely, multistability is expected to be more common in microbial ecosystems due to the growthyield tradeoff.
Interplay between diversity and stability in ecosystems with multiple essential nutrients
Ever since Robert May’s provocative question ‘Will a large complex system be stable?’ (May, 1972) the focus of many theoretical ecology studies has been on investigating the interplay between dynamic stability and species diversity in real and model ecosystems (Ives and Carpenter, 2007). May’s prediction that ecosystems with large number of species tend to be dynamically unstable needs to be reconciled with the fact that we are surrounded by complex and diverse ecosystems that are apparently stable. Thus, it is important to understand the factors affecting stability of ecosystems in general and microbial ecosystems in particular .
Here, we explored the interplay between diversity and stability in a particular type of microbial ecosystems with multiple essential nutrients. We discussed three criteria for stability of microbial communities shaped by the competition for nutrients: (i) how stable is the species composition of a community to fluctuations in nutrient supply rates; (ii) the extent of community’s resilience to species invasions; and (iii) its dynamical stability to small stochastic changes in abundances of existing species. Naturallyoccurring microbial communities may or may not be stable according to either one of these three criteria (Ives and Carpenter, 2007). The degree of importance of each single criterion is determined by multiple factors such as how constant are nutrient supply rates in time and space and how frequently new microbial species migrate to the ecosystem.
Our model provides the following insights into how these three criteria are connected to each other. First, as evident from Figure 1F, the three types of stability are largely independent from each other. Second, communities growing on a well balanced mix of nutrients tend to have high species diversity (see peak in Figure 4A). The similar effect was demonstrated in other consumer resource models (Posfai et al., 2017; Tikhonov and Monasson, 2017; Taillefumier et al., 2017; Marsland et al., 2019a). However, each of the community states in this regime tends to have a low structural stability with respect to nutrient fluctuations. In environments with highly variable nutrient supplies the community will frequently shift between these states. That is to say, some of the species will repeatedly go locally extinct and the vacated niches will be repopulated by others. Furthermore, many of the steady states in this regime are dynamically unstable giving rise to multistability and regime shifts. In this sense our model follows the general trend reported in May (1972). Conversely, microbial communities growing on an imbalanced mix of essential nutrients have relatively low diversity (Figure 4A) but are characterized by a high degree of structural and dynamic stability (see Figure 4D and Figure 4C respectively). We expect these trends to apply to a broad variety of consumerresource models.
The existence of dynamically unstable states always goes hand in hand with multistability (Scheffer and Carpenter, 2003) (see the dashed line in Figure 2D for an illustration of this effect in our model). Interestingly, in our model we always find $V1$ dynamically unstable states coexisting with $V$ dynamically stable ones for the same environmental parameters (see Figure 3C and Figure 3—figure supplement 1 for some examples). All states (both dynamically stable and unstable) shown in Figure 3C are positioned along some onedimensional curve in PCA coordinates. This arrangement hints at the possibility of a nonconvex onedimensional Lyapunov function whose $V$ minima (corresponding to stable states) are always separated by $V1$ maxima (unstable stable states) as dictated by the Morse theory (Milnor, 1963). This should be contrasted with convex multidimensional Lyapunov functions used in MacArthur (1970), Case and Casten (1979) and Chesson (1990).
Extensions of the model
Our model can be extended to accommodate several additional properties of reallife microbial ecosystems. First, one could include generalist species capable of using more than one nutrient of each type. The growth rate of such species is given by:
Here, the sum over $i$ (respectively $j$) is carried out over all carbon (nitrogen, respectively) sources that this species is capable of converting to its biomass. One may also consider the possibility of diauxic shifts between substitutable nutrient sources. In this case, each generalist species is following a predetermined preference list of nutrients and uses its carbon and nitrogen resources oneatatime, as modelled in Goyal et al. (2018). Since at any state each of the species is using a ‘specialist strategy’, that is to say, it is growing on a single carbon and a single nitrogen source, we expect that many of the results of the current study would be extendable to this model. Interestingly, the stable marriage problem can be used to predict the stable states of microbial communities with diauxic shifts between substitutable resources (Goyal et al., 2018) and those in communities growing on a mix of two essential nutrients as in this study. It must be pointed out that these models use rather different variants of the stable marriage model.
It is straightforward to generalize our model to Monod’s growth equation and to take into account nonzero death rate (or maintenance cost) of individual species (see Appendix 1).
One can extend our model to include crossfeeding between the species. In this case some of the nutrients are generated as metabolic byproducts by the species in the community. These byproducts should be counted among nutrient sources and thus would allow the number of species to exceed the number of externally supplied resources.
Above we assumed a fixed size of the species pool. This constraint could be modified in favor of an expanding pool composed of a constantly growing number of species. These new species correspond to either migrants from outside of the community or mutants of the species within the community. This variant of the model would allow one to explore the interplay between ecosystem’s maturity (quantified by the number of species in the pool) and its properties such as multistability and propensity to regime shifts.
Control of microbial ecosystems exhibiting multistability and regime shifts
In many practical situations we would like to be able to control microbial communities in a predictable and robust manner. That is to say, we would like to be able to reliably steer the community into one of its stable states and to maintain it there for as long as necessary. Alternative stable states and regimes shifts greatly complicate the task of manipulation and control of microbial ecosystems. Indeed, multistability means that the environmental parameters alone do not fully define the state of the community. In order to get it to a desired state, one needs to carefully select the trajectory along which one changes the environmental parameters (nutrient supply rates). Changing these parameters could lead to disappearance (local extinction) of some microbial species and open the ecosystem for colonization by others thereby changing its state. However, not all the states could be directly converted to each other in one step due to them being restricted to vastly different environments. Thus, densely interconnected networks of regime shifts shown in Figure 2E–F can be viewed as maps guiding the selection of the optimal environmental trajectory leading to the desired stable species composition.
These maps also suggest that microbial ecosystems described by our model might have a relatively small number of key drivers of regime shifts roughly corresponding to network modules (see Figure 2F). Regime shifts in each of the modules are driven by the competition between two mutually exclusive sets of keystone species. In addition to these keystone species, states also include ‘satellite’ species that do not generally affect the bistable switch. The exploration of different manipulation strategies of microbial ecosystems and the role of keystone and peripheral species in regime shifts is the subject of our future research (Maslov et al., unpublished).
Materials and methods
Identification of all states and classification of them as invadable or uninvadable
Request a detailed protocolThe competitive exclusion principle states that, in general, two species competing for the same growthlimiting nutrient cannot coexist with each other. Accounting for nonlimiting nutrients present in our model, the competitive exclusion principle can be reformulated as the following two rules:
Rule 1: In a given steady state each nutrient (either carbon or nitrogen) limits the growth of no more than one species.
Rule 2: Any number of species can use a given nutrient in a non growthlimiting fashion. However, each of such species needs to be able to survive at the steady state concentration of this nutrient set by the growthlimited species. That means that for every nutrient each of the non growthlimited species $\beta $ needs to be more competitive than the growlimited species $\alpha $ for the same resource: ${\lambda}_{\alpha}^{(c)}<{\lambda}_{\beta}^{(c)}$ (or ${\lambda}_{\alpha}^{(n)}<{\lambda}_{\beta}^{(n)}$ in case of a nitrogen nutrient).
Note that in any state of our model every species has a unique nutrient limiting its growth. By the virtue of the Rule 1, if a nutrient is limiting the growth of any species at all, such species is also unique. Hence, in a given state the relationship between surviving species and their growthlimiting nutrients (marked as shaded ovals in Figure 1A) is an example of a matching on a graph of resource utilization. Rule two imposes additional limitations on this matching. As we show in the Appendix 2 , uninvadable states correspond to stable matchings in a variant of the celebrated stable marriage problem (Gale and Shapley, 1962; Gusfield and Irving, 1989).
Just like in the MacArthur model (MacArthur and Levins, 1964) or any other consumerresource model for that matter, the number of species present in a steady state of the community cannot exceed the total number of nutrients they consume. Any community constructed using Rules 1 and 2 represents a steady state of the ecosystem feasible for a certain range of nutrient supply rates. This state can be either invadable or uninvadable, and either dynamically stable or not.
For simplicity in our simulations we use equal numbers of C and N resources ($L$ carbons and $L$ nitrogens), with one unique species capable of utilization of every pair of resources (${L}^{2}$ species in total). One must reiterate that our theory is not restricted to the specific values of $K$, $M$, and $S$. We first selected the values of ${\lambda}_{(i,j)}^{(c)}$ and ${\lambda}_{(i,j)}^{(n)}$ from a uniform random distribution between 10 and 100. Note that in the regime of high nutrient supply (${\varphi}^{(c,n)}\gg \frac{{\delta}^{2}}{{\lambda}^{(c,n)}}$) all steady states of the community can be identified and tested for invadability using only the relative rank order of species’ competitiveness for nutrients. For this we used the following exhaustive search algorithm:
Step 1  Select the subset of species whose growth is limited by C (Climited species). For every carbon nutrient there are $L$ ways to choose a species using this nutrient. There is also an additional possibility that this nutrient is not limiting the growth of any species. Thus, the total number of possibilities is $L+1$ for each of $L$ carbon nutrients. There are ${(L+1)}^{L}$ ways to choose the set of Climited species and our algorithm will exhaustively investigate each of these potential steady states onebyone.
Step 2  Given the set of Climited species selected in Step 1, we now select the set of Nlimited species. We first eliminate from our search any species that doesn’t have enough carbon to grow. That is to say, we go over all carbon nutrients onebyone and eliminate all species whose ${\lambda}^{(c)}$ is smaller than that of the Climited species (if any) for this carbon nutrient. Among the remaining species we go over the nitrogen nutrients onebyone and look for all possible ways to add a species limited by a given nitrogen source ${n}_{j}$ that satisfy the Rule 2. More specifically, we identify all species that use ${n}_{j}$ and can grow on their carbon sources (those are the only species that remained after the elimination procedure described above). We then compare ${\lambda}^{(n)}$s of these species to ${\lambda}^{(n)}$s of all Climited species using ${n}_{j}$. To satisfy the Rule 2 for each ${n}_{j}$ we can add at most one Nlimited species and its ${\lambda}^{(n)}$ has to be smaller than ${\lambda}^{(n)}$s of all Climited species using ${n}_{j}$. Let ${M}_{j}$ be the number of such species (${M}_{j}=0$ if there are no such species for a given ${n}_{j}$). The total number of possible steady states of our model for a given combination of Climited species selected in Step 1 is given by ${\prod}_{j=1}^{L}({M}_{j}+1)$. Here the $+1$ factor in ${M}_{j}+1$ takes into account an additional possibility to not add any Nlimited species for ${n}_{j}$.
The unique way to construct an uninvadable state by following this algorithm is to go over all nitrogen sources onebyone and for each of them add the Nlimited species with the largest ${\lambda}^{(n)}$ among all species using this resource, whose growth is allowed by carbon constraints. If for every ${n}_{j}$ this species is allowed by the Rule 2, that is to say, if its ${\lambda}^{(n)}$ is smaller than ${\lambda}^{(n)}$ of all Climited species using ${n}_{j}$, we successfully constructed a unique uninvadable state for a given set of Climited species. Indeed, all possible invading species that are allowed to grow by their carbon nutrients will be blocked by their nitrogen nutrients. If, however, for any of ${n}_{j}$, the species with the largest ${\lambda}^{(n)}$ is not allowed by the Rule 2, that is to say, if its ${\lambda}^{(n)}$ is larger than ${\lambda}^{(n)}$ of at least one of the Climited species, this species would make a successful invader of any state we construct. In this case, there is no uninvadable state for the set of Climited species selected during the Step 1.
We used the above procedure to identify all possible steady states and to classify them as invadable and uninvadable for different numbers of resources used in our 2C × 2N × 4S and 6C × 6N × 36S examples. Note that, while this method is computationally feasible for a relatively small number of nutrients (we were able to successfully use it for up to 9 nutrients of each type), for larger systems one should rely on computationally more efficient algorithms based on the stable marriage problem (Gale and Shapley, 1962; Gusfield and Irving, 1989) as described in the Appendix 3.
MonteCarlo sampling of nutrient supply rates to identify feasible ranges of states
Request a detailed protocolGiven the parameters defining all species (i.e., the set of their $\lambda $s and $Y$s) and the chemostat dilution constant $\delta $, each state $p$ is feasible within a finite region in the nutrient supply space (a $K+M$ dimensional space $\overrightarrow{\varphi}=\{{\varphi}_{i}^{(c)},{\varphi}_{j}^{(n)}\}$), where all microbial populations and nutrient concentrations are nonnegative and the limiting nutrients of every surviving species do not change. It is easy to show that in a steady state our system satisfies mass conservation laws for each of the nutrients:
To simplify the process of calculating the feasible volumes of all states we worked in the limit of high nutrient supply where ${\varphi}_{i}^{(c)}\gg \frac{{\delta}^{2}}{{\lambda}_{\alpha}^{(c)}}$ and ${\varphi}_{j}^{(n)}\gg \frac{{\delta}^{2}}{{\lambda}_{\alpha}^{(n)}}$ for all species $\alpha $. In this case the concentration $\delta /{\lambda}_{\alpha}^{(c,n)}$ of any nutrient limiting growth of some species ($\alpha $ in this case) is negligible compared to its ‘abiotic concentration’${\varphi}_{i}^{(c,n)}/\delta $, that is to say, its concentration before any microbial species were added to the chemostat. In this case one can ignore the terms ${c}_{i}$ and ${n}_{j}$ in Equation 4 for all nutrient limiting growth of some species and leave only the ones that are not limiting the growth of any species. It is convenient to introduce the K + M − dimensional vector ${\overrightarrow{X}}_{p}$ of microbial abundances and nonlimiting nutrient concentrations in a given state $p$. For example, for the uninvadable state #5 in the 2C × 2N × 4S model we have: ${\overrightarrow{X}}_{5}=\{{B}_{(1,1)},{B}_{(1,2)},{B}_{(2,2)},{n}_{2}\}$.
The mass conservation laws (Equation 4) can be used to obtain the feasible volumes of all states and can be represented in a compact matrix form for each state $p$:
where $\varphi $ is the vector of $K+M$ nutrient supply rates and $\widehat{{R}_{p}}$ is a matrix composed of inverse yields ${Y}^{1}$ of surviving species and '1’ for each of the nonlimiting nutrients in a given state $p$. For example, for the state #5 in our 2C × 2N × 4S model the Equation 5 expands to:
Using Equation 5, it is easy to check if a given state is feasible at a particular nutrient supply rate $\overrightarrow{\varphi}$ by multiplying ${\widehat{{R}_{p}}}^{1}$ (the inverse of the matrix $\widehat{{R}_{p}}$) with $\overrightarrow{\varphi}$. If all of the elements of the resulting vector ${\overrightarrow{X}}_{p}$ are positive, then the state $p$ is feasible at $\overrightarrow{\varphi}$. If the matrix $\widehat{{R}_{p}}$ is not invertible that is $det(\widehat{{R}_{p}})=0$, the state is feasible only on a lowdimensional subset of nutrient supply rates. This is not possible for a general choice of yields $Y$ and is not considered in our study.
The parameters $\lambda $ and $Y$ for 2C × 2N × 4S and 6C × 6N × 36S realizations of the model were drawn from the uniform random distributions and are listed in Supplementary files 1,2 and 3,4,5,6. For $\lambda $’s the distribution ranges between 10 and 100. For $Y$’s it is between 0.1 and 1.
In our numerical simulations, for each model realization, we sampled 10^{6} random nutrient supply rate combinations $\overrightarrow{\varphi}$. Supply rates of each individual nutrient were independently selected from a uniform distribution on the [10, 1000] interval. We refer to this procedure as Monte Carlo sampling. The lower bound ensures that the system is always in the limit of high nutrient supply since max $(\frac{{\delta}^{2}}{{\lambda}_{\alpha}})=0.1\ll 10.$.
Then we checked the feasibility of each of the 33 possible steady states (both invadable and uninvadable) in the 2C × 2N × 4S model and each of the 1211 uninvadable steady states in the 6C × 6N × 36S model. That is to say, for every set of nutrient supply rates $\overrightarrow{\varphi}$ and for every state $p$ we checked whether all elements of ${\overrightarrow{X}}_{p}$ are positive. The feasible range of nutrient supply rates of each state was estimated as the fraction of nutrient supply rate combinations (out of 1 million vectors $\overrightarrow{\varphi}$ sampled by our Monte Carlo algorithm) where this state was found to be feasible.
The network of regime shifts from overlaps of feasible ranges
Request a detailed protocolTwo stable states are said to be capable of a regime shift if their feasibility ranges overlap with each other, that is if there exists at least one nutrient supply rate combination at which both these states are feasible. We used the data obtained by the MonteCarlo sampling to look for such cases and to construct networks shown in Figure 2E, Figure 2F.
We performed additional simulations to look for boundaries between uninvadable states in our realization of the 2C × 2N × 4S model. In order to do that, for each state we generated a large ensemble of random vectors of bacterial abundances of surviving species and concentrations of all notlimiting nutrients. The population of one of the species (also randomly selected) was set to be a small negative number (−0.01). This represents continuous gradual extinction of this species upon crossing of the boundary, The populations of state’s other surviving species and the concentrations of its nonlimiting nutrients were drawn from the uniform distribution between (0,1]. Using Equation 5, we calculated the nutrient supply rates $\varphi $ corresponding to this case. For these nutrient supply rates lying just across the feasibility boundary of the originally selected state, we checked the feasibility of other five uninvadable dynamically stable states. If any of these states ended up being feasible, we assumed that this state shares a boundary with the originally selected one. Using these procedure we found four bordering pairs of states shown as red edges in Figure 2E.
We used Gephi 0.9.2 software package to visualize the network in Figure 2F and to perform its modularity analysis. Seven densely interconnected clusters shown with different colors in Figure 2F were identified using Gephi’s builtin moduledetection algorithm (Blondel et al., 2008) with the resolution parameter set to 1.5.
Dynamic stability of states
Request a detailed protocolWe checked the dynamic stability of every 33 possible steady state (both invadable and uninvadable) (for the 2C × 2N × 4S model) and each of the 1211 uninvadable states (for the 6C × 6N × 36S model) using the following two algorithms:
Small perturbation analysis For 2C × 2N × 4S example and each of the 33 states we selected many supply rate combination where this state is feasible. For each of these supply rates we generated the populations in our ecosystem to be equal to their steady state values. We then subjected them to small perturbations of of all nutrient concentrations and of all microbial populations of species present in the state. If after some transient period all populations and concentrations returned to their steady state values  the state was declared to be dynamically stable. If they drifted to these in some other steady state  the original state was declared to be dynamically unstable. Based on our numerical simulations the dynamical stability of the state was independent of the nutrient supply rates at which this numerical experiment has been performed. We choose to perturb only the populations of the species present in the state because an invadable state, by definition, would always be dynamically unstable against an addition of a very small population of at least one invading species from the species pool. This instability should not render it dynamically unstable. The numerical integration of the system dynamics following a perturbation was done in C programming language using the CVODE solver library of the SUNDIALS package (Hindmarsh et al., 2005).
Inference of state’s dynamic stability from the pattern of its overlaps with other states The number of uninvadable states (1211) in our 6C × 6N × 36S model was too large to be directly tested for dynamic stability as we did for the 2C × 2N × 4S model. Their dynamic stability was instead inferred from our MonteCarlo simulations listing all feasible uninvadable states for every sampled nutrient supply rate combination. We first identified 1022 uninvadable states, which were the unique feasible state for at least one nutrient supply point. All such states should be dynamically stable, since for every nutrient supply rate there should be at least one uninvadable dynamically stable state representing the end point of system’s dynamics. The remaining 173 uninvadable states, which were feasible for at least one of 1 million sampled nutrient supply rates were labelled as potentially dynamically unstable. Note that in our MonteCarlo analysis we only sampled a finite (albeit large) number of nutrient supply rate combinations. Thus it is entirely possible that we missed some crucial supply rate combinations for which one of these states was the only uninvadable state. Any such point would have rendered this state as dynamically stable. Such false assignments might lead to a violation of the basic empirical rule in our model stating that $V$ uninvadable stable states are always accompanied by $V1$ uninvadable dynamically unstable states ($V/(V1)$ rule) for some sampled nutrient supply rates. In our Monte Carlo simulations of the 6C × 6N × 36S model the $V/(V1)$ rule was violated for only 370 nutrient supply rates combinations out of 1,000,000 sampled points. We believe that these violations were caused by an incorrect identification of dynamically unstable states mentioned above. To iteratively refine the lists of stable and unstable states, we went over all potentially unstable states onebyone and checked whether reclassifying the state involved in the largest number of violations as stable would reduce the overall number of violations. If it did, we reclassified this state as stable and recalculated the number of violations for all remaining points. By the end of this iterative procedure we were able to completely eliminate violations by reassigning 36 potentially unstable states as dynamically stable. This left us with 1022 + 36 = 1058 dynamically stable and 173 − 36 = 137 dynamically unstable uninvadable states in the 6C × 6N × 36S model. The remaining 1211 − 1058 − 137 = 16 uninvadable states were not feasible for any of 1,000,000 sampled nutrient supply rates. Hence their dynamic stability remains unidentified. Both 36 reassigned states and 16 undetected states are expected to have very small ranges of feasible nutrient supply rates.
Multistability as a function of variation in stoichiometric ratios of different species
Request a detailed protocolTo investigate how multistability in our model depends on variation in stoichiometric ratios of different species, we simulated 4000 variants of the 2C × 2N × 4S model. In these variants, we kept the same choice of species competitiveness (quantified by their $\lambda $’s) but reassigned their yields $Y$. To cover a broad range of standard deviations of N:C stoichiometry of different species (their ${Y}_{\alpha}^{(c)}/{Y}_{\alpha}^{(n)}$) we randomly sampled yield combinations from gradually expanding intervals. First we simulated 1000 model variants, where yields of four species were independently drawn from uniform distribution $U(0.45,0.55)$. These simulations were followed by 1000 model variants, where yields of four species were drawn from $U(0.3,0.7)$, 1000 model variants with yields from $U(0.1,0.9)$ and, finally, 1000 model variants with yields from $U(0.01,1.0)$. In each variant of the model with a particular set of yields of four species, we calculated the fraction of multistable points among 10^{5} nutrient supply rate combinations as described in the section 5.2 MonteCarlo sampling of nutrient supply rates to identify feasible ranges of states of Materials and methods. The results are shown in Figure 5.
GitHub repository of the code used in our project
Request a detailed protocolThe PCA analysis, plots and statistical tests were implemented using R version 3.4.4. Other simulations were carried out in C (using compiler gcc version 5.4.0) and Python 3.5.2. Matlab analysis was done using MATLAB and Statistics Toolbox Release 2018a, The MathWorks, Inc, Natick, Massachusetts, United States. The code for both our simulations and statistical analysis can be downloaded from: https://github.com/ssm57/CandN (Dubinkina and Maslov, 2019; copy archived at https://github.com/elifesciencespublications/CandN).
Appendix 1
General form of growth laws
It is straightforward to generalize our model to allow a more general functional form for growth laws than Liebig’s law, $\mathrm{min}({\lambda}_{\alpha}^{(c)}{c}_{i},{\lambda}_{\alpha}^{(n)}{n}_{j})$. Microbial growth on two essential substrates is thought to normally follow Monod’s equation for the ratelimiting nutrient: ${g}_{\alpha}^{(m)}\mathrm{min}({c}_{i}/({K}_{\alpha}^{(c)}+{c}_{i}),{n}_{j}/({K}_{\alpha}^{(n)}+{n}_{j}))$ (See KovárováKovar and Egli, 1998 for a discussion of limitations of Monod’s law). For low concentrations of the ratelimiting nutrient, say carbon source, the Monod’s law simplifies to the linear growth law used throughout this study: ${g}_{\alpha}={\lambda}_{\alpha}^{(c)}{c}_{i}$. Microbes’ competitive abilities, also known as their specific affinities towards each substrate, are related to the parameters of Monod’s law via
In another variant of growth laws, two essential nutrients at low concentrations jointly affect the growth rate of the microbe: ${g}_{\alpha}^{(m)}{c}_{i}\cdot {n}_{j}/[({K}_{\alpha}^{(c)}+{c}_{i})\cdot ({K}_{\alpha}^{(n)}+{n}_{j})]$ (see Bader, 1978 for a discussion of these and other forms of doublesubstrate growth law). For simplicity of mathematical calculation we limited this study to Liebig’s law. However, many of the essential results we obtained (e. g. possible multistability in a system where species have different yields) hold for any growth laws listed above. In fact, the low concentration version of the previous growth law, where ${g}_{\alpha}^{(m)}{c}_{i}\cdot {n}_{j}$ has been studied by one of us in the context of autocatalytic growth of heteropolymers (Tkachenko and Maslov, 2018). Instead of exponentially replicating microbial species Tkachenko and Maslov (2018) considers pairs of mutually catalytic (and thus exponentially growing) complementary ‘2mers’ (a specific sequence of two consecutive monomers anywhere within a polymer chain). This minor difference complicates the math, while leaving the basic properties unchanged. Just like in our system, where up to $2L$ species (out of ${L}^{2}$ candidates) may simultaneously survive in the steady state of an ecosystem grown on of $L$ carbon and $L$ nitrogen sources, the polymer systems have no more than $2Z$ 2mer ‘species’ (out of ${Z}^{2}$ candidates) surviving in the steady state with polymers having $Z$ possible monomers on their right ends and $Z$ possible monomers on their left ends. Many (but not all) results of this paper are largely consistent with the present study. Note that for polymers the yields of all ‘species’ are equal to 1, that is to say, one new 2mer is formed upon ligation of one left end of a polymer with one right of another polymer chain. Yet, the model in Tkachenko and Maslov (2018) is capable of (at least) bistability. At present, it is not clear if this is due to autocatalytic cycles having length two or this property would survive in a simpler version of the model in which instead of the Equation (1) of Tkachenko and Maslov (2018) one has
and the overall fluxes of left and right ends are independent from each other (instead of both being equal to ${c}_{i}={\varphi}_{i}/\delta $ as in Tkachenko and Maslov, 2018).
Another variant of the model is where each species $\alpha $ has its own unique ‘death’ or ‘maintenance’ rate ${\delta}_{\alpha}$, playing the role of the same dilution rate $\delta $. The steady states of this model (but not the dynamics leading to these states) can be calculated by dividing both sides of Equations 2 by ${\delta}_{\alpha}$. This is equivalent by redefining the competitiveness parameters to $\stackrel{~}{{\lambda}_{\alpha}}={\lambda}_{\alpha}/{\delta}_{\alpha}$ and setting the chemostat dilution rate to $\stackrel{~}{\delta}=1$. All of our results in the highflux regime $\varphi \gg {\delta}^{2}/\lambda $ would remain unchanged.
From (Equation 5Equation 6) one can see that when all species have the same C:N stoichiometry, the maximal number of microbialspecies in a state is equal to the number of nutrients minus 1. Indeed, one can show that a state $p$ with ${S}_{surv}=K+M$ has $det(\widehat{{R}_{p}})=0$, which means that the feasible volume of any such state is zero. These states are only possible on a lowerdimensional manifold in the $(K+M)$dimensional space of supply rates (these results have been already discussed by Tilman in his special case; Tilman, 1982).
Multistability is also possible in a variant of the MacArthur model (MacArthur and Levins, 1964; MacArthur, 1970; Chesson, 1990) in which different species have different yields on individual carbon sources (Maslov, unpublished). A convex Lyapunov function (MacArthur, 1970) precluding multistability does not exist in this case.
Appendix 2
Constraints on steady states from microbial and nutrient dynamics
A steady state of equations describing the microbial dynamics (Equation 2) is realized when either ${B}_{\alpha}=0$ (the species was absent from the system from the start or subsequently went extinct) or when its growth rate ${g}_{\alpha}$ is exactly equal to the chemostat dilution rate $\delta $. This imposes constraints on steady state nutrient concentrations with the number of constraints equal to the number of microbial species present with nonzero concentrations. Since, in general, the number of constraints cannot be larger than the number of constrained variables, no more than $K+M$ of species could be simultaneously present in a steady state of the ecosystem. For Liebig’s growth law used in this study, each resource can have no more than one species for which this resource limits its growth, that is to say, which sets the value of the minimum in $\mathrm{min}({\lambda}_{\alpha}^{(c)}{c}_{i},{\lambda}_{\alpha}^{(n)}{n}_{j})$ The steady state concentrations of these resources are given by ${c}_{i}^{(*)}=\delta /{\lambda}_{\alpha}^{(c)}$ (if the growth is limited by the carbon source) and ${n}_{j}^{(*)}=\delta /{\lambda}_{\alpha}^{(n)}$ (if the growth is limited by the nitrogen source). Here $\alpha $ is the species whose growth is ratelimited by the resource in question. In a general case, no more than one species can be limited by the same resource (carbon in our example), since the species with the largest ${\lambda}^{(c)}$ would outcompete other species with smaller values of ${\lambda}^{(c)}$ by making the steady state concentration ${c}_{i}^{(*)}$ so low that other species can no longer grow on it. Note however, that multiple species $\beta $ could consume the same resource as the ratelimiting species $\alpha $, as long as their growth is not limited by the resource. Each of these species must then be limited by their other nutrient (a nitrogen source in our example). However, their survival requires that carbon concentration set by species $\alpha $ is sufficient for their growth. Thereby, any species growing on a resource in a nonlimited fashion must have ${\lambda}_{\beta}^{(c)}>{\lambda}_{\alpha}^{(c)}$.
Mathematically, it cane be proven by observing that, since species $\beta $ is limited by its nitrogen resource, one must have ${\lambda}_{\beta}^{(c)}{c}_{i}^{(*)}>{\lambda}_{\beta}^{(n)}{n}_{j}^{(*)}$. At the same time in a steady state, the concentrations of all ratelimiting resources are determined by the dilution rate $\delta $ via ${\lambda}_{\beta}^{(n)}{n}_{j}^{(*)}=\delta $, and ${\lambda}_{\alpha}^{(c)}{c}_{i}^{(*)}=\delta $. Combining the above three expressions one gets: ${\lambda}_{\beta}^{(c)}{c}_{i}^{(*)}>{\lambda}_{\beta}^{(n)}{n}_{j}^{(*)}=\delta ={\lambda}_{\alpha}^{(c)}{c}_{i}^{(*)}$, or simply ${\lambda}_{\beta}^{(c)}>{\lambda}_{\alpha}^{(c)}$. The constraints on competitive abilities $\lambda $ for species present in a steady state in our model are then:
Exclusion Rule 1: Each nutrient (either carbon or nitrogen source) can limit the growth of no more than one species $\alpha $. From this it follows that the number of species coexisting in any given steady state cannot be larger than $K+M$, the total number of nutrients.
Exclusion Rule 2: Each nutrient (e.g. specific carbon source) can be used by any number of species in a nonratelimiting fashion (that is to say, where it does not constrain species growth in Liebig’s law). However, any such species $\beta $ has to have ${\lambda}_{\beta}^{(c)}>{\lambda}_{\alpha}^{(c)}$, where ${\lambda}_{\alpha}^{(c)}$ is the competitive ability of the species whose growth is limited by this nutrient. In case of a nitrogen nutrient, the constraint becomes ${\lambda}_{\beta}^{(n)}>{\lambda}_{\alpha}^{(n)}$.
Note that the steady state solutions of equations Equation 2 do not depend on populations ${B}_{\alpha}$ of surviving species. Their steady state populations ${B}_{\alpha}^{(*)}$ are instead determined by Equation 3. Taking into account that, in a steady state, the growth rate of each surviving species is exactly equal to the dilution rate $\delta $ of the chemostat, after simplifications one gets:
As described above, the steady state concentration of resources are given by $\delta /{\lambda}_{\alpha}^{(\text{c or n})}$, where $\alpha $ are the species ratelimited by each resource. In the absence of such species, the concentration of a resource is given by anything left after it being consumed by surviving species in a nonratelimiting manner. One can show that in this case, the resource (e.g. carbon) concentration has to be larger than $\delta /{\lambda}_{\beta}^{(c)}$, where ${\lambda}_{\beta}^{(c)}$ is the smallest affinity among microbes utilizing this resource.
One convenient approximation greatly simplifying working with Equation S2 is the ‘highflux limit’ in which ${\varphi}_{i}^{(c)}\gg {\delta}^{2}/{\lambda}_{\alpha}^{(c)}$ and ${\varphi}_{j}^{(n)}\gg {\delta}^{2}/{\lambda}_{\alpha}^{(n)}$. In this approximation one can approximately set to zero the steady state concentrations of all resources that have a species ratelimited by them. The steady state concentrations of the remaining resources can take any value as long as it is positive. Hence, in this limit the Equation S2 can be viewed as a simple matrix test of whether a given set of surviving species limited by a given set of resources is possible for a given set of nutrient fluxes. Indeed, my multiplying the vector of fluxes with the inverse of the matrix $\widehat{R}$ composed of inverse yields of surviving species and one for nutrients not limiting the growth of any species one formally gets the only possible set of steady state species abundances, ${B}_{\alpha}^{(*)}$, and a subset of nonlimiting resource concentrations ${c}_{i}^{(*)}$ and ${n}_{j}^{(*)}$. If all of them are strictly positive  the steady state is possible. If just one of them enters the negative territory  the steady state cannot be realized for these fluxes of nutrients.
The above rule can be modified to apply even below the highflux limit with the following modifications: 1) Instead of ${\varphi}^{(c)}$ (or ${\varphi}^{(n)}$), one uses their ‘effective values’ ${\stackrel{~}{\varphi}}^{(c)}$ (or ${\stackrel{~}{\varphi}}^{(n)}$) introduced in Goyal et al. (2018), determined as
where $\alpha (i)$ is the (unique) species limited by the nutrient $i$. If the nutrient is not limiting for any os the species in the steady state, $\alpha (i)$ is the species using the nutrient in a nonlimited fashion, which has the smallest value of $\lambda $. This last rule comes from the observation that in order for a nonlimiting resource not to become limiting for a species $\beta $ currently using it in a nonlimiting fashion, its concentration cannot fall below $\delta /{\lambda}_{\beta}^{(x)}$. Thus, when checking the feasibility of a given state, the concentration of a nonlimiting resource can be written as $\delta /{\lambda}_{\beta}^{(x)}+$ a positive number, or (more conveniently) the influx of this resource can be offset as described in Equation S3
Appendix 3
Stable matching approach for identification and classification of steady states
First we describe the exact onetoone mapping between all uninvadable steady states (UIS) in our model and the complete set of ‘stable marriages’ in a variant of a wellknown stable marriage or stable allocation problem developed by Gale and Shapley in the 1960s (Gale and Shapley, 1962) and awarded the Nobel prize in economics in 2012. This mapping provides us with constructive algorithms to identify and count all uninvadable steady states in our ecosystem.
We start by considering a special case of our problem with $L$ carbon and $L$ nitrogen sources and a pool of ${L}^{2}$ species, such that for every pair of sources ${c}_{i}$ (carbon) and ${n}_{j}$ (nitrogen) there is exactly one microbe ${B}_{ij}$ capable of using them. For the sake of simplicity we have switched the notation from ${B}_{\alpha}$ to ${B}_{ij}$, where $\alpha =(ij)$ is the unique microbe in our pool capable of growing on ${c}_{i}$ and ${n}_{j}$. Having considered this simpler situation we will return to the most general case of unequal numbers of carbon ($K$) and nitrogen ($M$) resources and any number of microbes from a pool of $S$ species competing for a given pair of resources.
This is where we need to revise in certain ways the network representation of a steady state used in the main text (see Figure 1A). In the marriage game related theory, the notion of (stable or unstable) matching explicitly refers to a bipartite graph with two distinct sets of vertices and edges arranged in such a way that each one may join only a pair of elements belonging to different sets. In our case, it is natural to consider two sets of resource nodes (vertices), one including all carbon nodes and the other one containing all nitrogen nodes. An edge, or link, will appear between a carbon ${c}_{i}$ and nitrogen ${n}_{j}$ node if the microbe ${B}_{ij}$ using these two nutrients is present in the state represented by this particular bipartite network.
Furthermore, the specifics of our version of 'marriage game’, or rather 'residents vs hospitals’, problem requires us to consider directed bipartite graphs as the steady state representations in our model. For any species ${B}_{ij}$ present in a given state, we choose the direction of the edge joining node ${c}_{i}$ with node ${n}_{j}$ to be pointing from ${c}_{i}$ to ${n}_{j}$ if the microbe is limited by its carbon nutrient (and the other way around, from ${n}_{j}$ to ${c}_{i}$, on the case of ${B}_{ij}$ being nitrogenlimited). Appendix 3—figure 1A shows the directed bipartite graph representation of the state #5 of a particular example of 2C × 2N × 4S system considered in the 'Results’ section of the main text.
In what follows we will refer to a resource as occupied if in a given steady state there is a microbe for which this resource is ratelimiting. In our bipartite network representation occupied resources have an outgoing edge (their outdegree is equal to 1), while unoccupied resources have outdegree equal to 0.
Review of results about stable matchings in the hospitals/residents problem
The hospitals/residents problem (Gale and Shapley, 1962) is known in various settings. The one directly relevant to our problem is the following. There are $L$ applicants for residency positions in $H\le L$ hospitals. A hospital number $i$ has ${V}_{i}$ vacancies for residents to fill, ${V}_{i}$ ranging from zero to $L$, $\sum {V}_{i}=L$. Each hospital has a list of preferences in which residency applicants are strictly ordered by their ranks, from 1 (the most desirable) to $L$ , (the least desirable). These lists are generally different for different hospitals. Each applicant has a ranked list of preferred hospitals ranging from 1 (the most desirable) to $H$ (the least desirable). Those lists can also vary between applicants. A matching is an assignment of applicants to hospitals such that all applicants got residency and all hospital vacancies are filled. A matching is unstable if there is at least one applicant $a$ and hospital $h$ to which $a$ is not assigned such that:
Condition 1. Applicant $a$ prefers hospital $h$ to his/her assigned hospital;
Condition 2. Hospital $h$ prefers applicant $a$ to at least one of its assigned applicants.
If such a pair $(a,h)$ exists, it is called ‘a blocking pair’ or ‘a pair that blocks the matching’. A stable matching by definition has no blocking pairs. Gale and Shapley proved that for any set of applicant/hospital rankings and hospital vacancies there is at least one stable matching (Gale and Shapley, 1962). Generally the number of stable matching is larger than one. For example, for stable marriages and random rankings the average number of stable matchings is given by $L/e\mathrm{log}L$(Gusfield and Irving, 1989). To the best of our knowledge, the dependence of this number on the distribution of hospital vacancies has not been investigated. The fact that the actual number of uninvadable states is rather close to its lower bound indicates that, at least for $L\le 9$, the number of stable matchings averaged over all possible indegree allocations is rather close to 1.
Gale and Shapely not only proved the existence of at least one stable matching, but also proposed a constructive algorithm on how to find it. Listed below are the main steps in this algorithm optimized for for applicants. Each applicant first submits his/her application to the hospital ranking 1 in his/her preference lists. Each hospital considers all applications it received so far and accepts all of the applicants if their number is less or equal than hospital’s announced number of vacancies, ${L}_{i}$. If the number of applicants exceeds ${L}_{i}$, the hospital gives a conditional admission to the bestranking ${L}_{i}$ applicants according to hospital’s own preference list. Each applicant not admitted to their top hospital goes a step down on his/her preference list and applies to the secondbest hospital. The latter admits this applicant if (1) this hospital has not yet filled all of its vacancies or (2) all vacancies are filled, but among the conditionally admitted applicants there is at least one who ranks lower (according to hospital’s list) than the new applicant. Such lowerranked applicants are declined admission and replaced with better ones. They subsequently lower their expectations and apply to the next hospital on their list. After a number of iterations all applicants are admitted and all vacancies are filled so that this process stops. As Gale and Shapley (1962) proved, the resulting matching is stable. Furthermore, the theorem states that in this matching every applicant gets admitted to the best hospital among all stable matchings, while every hospital gets the worst set of residents among all stable matchings. Later research described in Gusfield and Irving (1989) describe more complex constructive algorithms allowing one to efficiently find all of the stable matchings starting with the applicantoptimal one.
Well developed mathematical apparatus of stable matching problem provides an invaluable help in the task of identifying all uninvadable states in microbial ecosystems. Indeed, without its assistance this task would require exponentially long time. To connect the problem of finding all uninvadable states to that of finding all stable matchings between hospitals and residents, we start with the following three observations:
In any uninvadable steady state, either all carbon sources or all nitrogen sources (or both) are occupied. Indeed, if in a steady state a carbon source ${c}_{i}$ and a nitrogen source ${n}_{j}$ are notlimiting to any microbes, then microbe ${B}_{ij}$ can always grow and thereby invade this state. Thus uninvadable states can be counted separately: one first counts the states where all nitrogen sources are occupied, and then counts those in which all carbon sources are occupied. Double counting happens when both carbon and all nitrogen sources are occupied. We will keep the possibility of double counting in mind and return to this problem later.
For a pool of species, where for every pair of resources there is exactly one microbe using each (carbon, nitrogen) pair, one can think of each of $L$ carbon (alternatively, nitrogen) sources as if it had a list of ‘preferences’ ranking all nitrogen (correspondingly carbon) sources. Indeed, the ranking of competitive abilities ${\lambda}_{ik}^{(c)}$ of different microbes using the same carbon source ${c}_{i}$ but different nitrogen sources ${n}_{k}$ can be viewed as the ranking of nitrogen sources $k$ by the carbon source $i$. Conversely, the ranking of ${\lambda}_{mj}^{(n)}$ with the same ${n}_{j}$ but variable ${c}_{m}$ can be thought of as ranking of carbon sources ${c}_{m}$ by the nitrogen source ${n}_{j}$.
Consider a steady state in which all nitrogen sources are occupied. In our network representation it corresponds to every nitrogen source sending an outgoing link to some carbon source. Let ${L}_{i}$ be the number of microbes using the carbon source $i$ in a nonlimiting fashion (the indegree of these outgoing links ending on ${c}_{i}$, see Appendix 3—figure 1C). Then, obviously, $L=\sum {L}_{i}$ (note that some of the terms in this sum might be equal to zero).
One can prove that if the state is uninvadable, then the matching given by all edges going from nitrogen sources to carbon sources must be stable in the GaleShapley sense. To prove this, let’s think of nitrogen sources as ‘applicants’ and carbon sources as ‘hospitals’ with their numbers of ‘vacancies’ given by ${L}_{i}$. Indeed, any unstable matching has at least one blocking pair $({n}_{j},{c}_{i})$ such that:
Condition 1. The nitrogen source (‘applicant’) ${n}_{j}$‘prefers’ the carbon source (‘hospital’) ${c}_{i}$ to its currently assigned carbon source (the one used by the current microbe ${B}_{kj}$ limited ${n}_{j}$). This means that ${\lambda}_{ij}^{(n)}>{\lambda}_{kj}^{(n)}$. Thus the microbe ${B}_{ij}$ can grow on its nitrogen source (provided that it can also grow on its carbon source).
Condition 2. The carbon source (‘hospital’) ${c}_{i}$‘prefers’ the nitrogen source (‘applicant’) ${n}_{j}$ to at least one of ${L}_{i}$ of its currently assigned carbon sources (the set of microbes using ${c}_{i}$ in a nonratelimiting fashion). Thereby ${\lambda}_{ij}^{(c)}$ must be larger than the smallest ${\lambda}^{(c)}$ among these microbes. According to the Exclusion Rule 2, this smallest ${\lambda}^{(c)}$ is still larger than ${\lambda}^{(c)}$ of the microbe limited by ${c}_{i}$ (if it exists). Thus the microbe ${B}_{ij}$ can also grow on its carbon source.
This proves that the microbe ${B}_{ij}$ corresponding to any blocking pair can grow on both its carbon and its nitrogen sources, and thereby can successfully invade the steady state. This finishes the proof that any uninvadable state has to be a stable matching in the GaleShapley sense.
However, this does not prove that any stable matching corresponds to exactly one uninvadable state. To prove this we first notice that, up to this point, our candidate uninvadable state contained only the nitrogenlimited species (See Appendix 3—figure 1B) . We will now supplement it with carbonlimited species in such a way that (1) added species do not violate the exclusion rule 2; (2) added species render the state completely uninvadable. Let is introduce a new notation (applicable to our case in which all nitrogen sources are occupied). Let ${\lambda}_{\text{min}}^{(c)}(i)$ denote the smallest ${\lambda}^{(c)}$ among all species using ${c}_{i}$ in a nonratelimiting fashion. The GaleShapley theorem only guarantees the protection of our state from invasion by a species $(i,j)$ with ${\lambda}_{ij}^{(c)}$ larger than ${\lambda}_{\text{min}}^{(c)}(i)$ (see the Condition 2 above). To ensure that our state is uninvadable by the rest of the species, one needs to add some carbonlimited species to this state. In order to do this in a systematic way, for each ${c}_{i}$ we compile the list of all species using this carbon source with ${\lambda}^{(c)}<{\lambda}_{\text{min}}^{(c)}(i)$. Each of these species is a potential invader. Some species could be crossed off from the list of potential invaders because they cannot grow on their nitrogen source. These species have ${\lambda}^{(n)}$ below that of the (unique) species limited by their nitrogen source. Among the species that remained on the list of invaders after this procedure, we select that with the largest ${\lambda}^{(c)}$ and add it to our steady state as a $C\to N$ directed edge, that is to say, as a carbonlimited species. This will prevent all other potential invaders on our list, since they have smaller ${\lambda}_{(c)}$ and thus, following the addition of our top carbonlimited species, they would no longer be able to grow based on their carbon source. We will go over all ${c}_{i}$ and add such carbonlimited species if they are needed. The only scenario when such species is not needed if our list of potential invaders would turn up to be empty. In this case we will leave this carbon source unoccupied. See Appendix 3—figure 1C for the illustration of an uninvadabe state constructed by the above procedure in $6C\times 6N\times 36S$ model. Since for each carbon source the above algorithm selects the carbonlimited species (or selects to add no such species) in a unique fashion, there is a single uninvadable state for every stable matching in the GaleShapley sense. We are now in a position to predict and enumerate all uninvadable states in our model.
Lower bound on the number of uninvadable states
To count the number of partitions $({L}_{1},{L}_{2},\mathrm{\dots},{L}_{L})$ such that $\sum {L}_{i}=L$, one can use a well known combinatorial method. According to this method, one introduces $L1$ identical ‘separators’ (marked with ) which are placed between $L$ identical objects (marked ·) separating them into $L$ (possibly empty) partitions. For example, for $L=4$ a partition $0,1,0,3$ would be denoted as $\cdot \mathrm{\cdots}$. The combinatorial number of all possible arrangements of separators and objects is obviously $\left(\genfrac{}{}{0pt}{}{2L1}{L}\right)$. For every such partition the GaleShapley theorem guarantees at least one stable matching (that is, at least one uninvadable steady state). The lower bound on the number of uninvadable steady states has to be doubled to account for reversal of roles of carbons and nitrogens. There is a small possibility that we double counted one partition $(1,1,\mathrm{\dots},1)$. Indeed, the unique uninvadable stable state corresponding to this partition could in principle be counted both when we start from nitrogen sources and when we start from carbon sources. While such uninvadable steady states can and do occur, if only under rather special circumstances described below, they cannot be obtained by the classical Gale and Shapley resident(nitrogen)proposing algorithm. Indeed, at the last step of the nitrogenproposing algorithm all carbon sources except just one carbon sources have all their vacancies already filled. Since the algorithm makes all nitrogen sources to go down their preference lists, at this point none of the carbon source whose outgoing arrow has already been accepted would prefer the last nitrogen with a vacancy to the one they have sent an arrow to. Thus, in terms of our adaptation of the GaleShapley algorithm to identify all uninvadable states, the last carbon to receive the incoming link from the last nitrogen does not need to be occupied (send a link from C to N) to prevent an invasion, as there is no one willing to invade it because of nitrogen preferences. Thus, in the uninvadable state identified by the Gale and Shapley nitrogenproposing algorithm, this carbon source would have an incoming arrow but no outgoing arrow. Such states would not be double counted and would not affect our lower bound on the total number of uninvadable states in the model. A partition in which each carbon and each nitrogen source have exactly one vacancy could have more than two (possibly identical) stable states found by the Gale and Shapley nitrogenproposing and carbonproposing algorithms could lead to double counting. Hence, in general one may need to merge a handful of identical states independently found starting from carbon and nitrogen sides.
Then we have ${N}_{UIS}\ge 2\left(\genfrac{}{}{0pt}{}{2L1}{L}\right)=\left(\genfrac{}{}{0pt}{}{2L}{L}\right)$. The Sterling approximation for this expression is ${2}^{2L}/\sqrt{\pi L}$. Thus the overall lower bound for the number of uninvadable stable states is given by
More generally, the number of carbon sources, $K$, is not equal to the number of nitrogen sources, $M$. The resource type with a larger number will always have at least one resource left without input. Thus here one never needs to correct for double counting. Using the same reasoning as for $K=M=L$, the lower bound on the number of resources in this case is given by $\left(\genfrac{}{}{0pt}{}{K+M1}{K1}\right)+\left(\genfrac{}{}{0pt}{}{K+M1}{M1}\right)=\left(\genfrac{}{}{0pt}{}{K+M}{K}\right)$. Here, the first term counts the uninvadable steady states in which all nitrogen sources are occupied and the partition divides $M$ edges sent by nitrogen sources among $K$ carbon sources, which requires $K1$‘dividers’. The second term counts the number of uninvadable steady states in which all carbon sources are occupied. Denoting the fraction of carbon resources among all resources as $p=K/(K+M)$ and using the Stirling approximation one gets
In the case of multiple microbial species using the same pairs of resources, our version of the GaleShapley residentoriented algorithm must be further updated. Let $M$ be the number of nitrogen sources, and $K$ — the number of carbon sources in the ecosystem, $S$ the number of species in our pool, each requiring a pair of resources to grow. As now there may be more than one microbe that uses a given pair of resources ${c}_{i}$ and ${n}_{j}$, we introduce the notation ${B}_{ij}^{(r)}$ for the $r$th microbe using the same pair of sources ${c}_{i}$ and ${n}_{j}$. On average, each nitrogen (carbon) source has $S/K$ ($S/M$) microbes, which are capable of using it. As in the traditional GaleShapley algorithm, each nitrogen (carbon) source ranks all microbes capable of using it by their ${\lambda}^{(n)}$ (${\lambda}^{(c)}$).
The way to identify all uninvadable stable states in this case is determined by a variant of the stable marriage problem (or rather the hospital/resident problem) in which every man (and every woman) may have more than one way to propose marriage to the same woman (man). In our model, this corresponds to more than one microbe (a type of marriage) capable of growing on the same pair of carbon (corresponding to, say, men) and nitrogen (corresponding to women) sources. You may think of it as if each participant has several different ways to propose to the person of the opposite sex (send flowers, take to a restaurant, etc). Each of these proposals is ranked by both parties independent of other ways. As far as we know, this variant has not been considered in the literature yet. However, all of the results of the usual stable marriage (or hospitalresident) problem remain unchanged.
One can easily see that our lower bound (Equation S5) on the number of uninvadable states (equal to the number of stable marriages in all partitions) remains unchanged. Indeed, it is given by the number of partitions and hence depends only on $K$ and $M$ and not on $S$. However, for $S\gg K\cdot M$ one expects to have many more stable marriages for each partition. Thus the lower bound we have established is likely to severely underestimate the actual number of UIS in the ecosystem.
Appendix 4
Conditions of multistability in the 2C × 2N × 4S ecosystem
Below we consider the general case of a 2C × 2N × 4S ecosystem. Let’s assume that the selected set of $\lambda $ parameters allow potentially unstable state in which all four species are present. These species form a single loop in the network representation (like state S7 in our 2C × 2N × 4S example). Without loss of generality we may assume that $\lambda$–parameters satisfy ${\lambda}_{11}^{(n)}>{\lambda}_{21}^{(n)};{\lambda}_{21}^{(c)}>{\lambda}_{22}^{(c)};{\lambda}_{22}^{(n)}>{\lambda}_{12}^{(n)};{\lambda}_{12}^{(c)}>{\lambda}_{11}^{(c)}$. The 4species loop is then formed by the links ${C}_{1}\to {N}_{1}$, ${N}_{1}\to {C}_{2}$, ${C}_{2}\to {N}_{2}$, and ${N}_{2}\to {C}_{2}$ In all other cases we may rename C and N resources until the direction of the loop is as stated above. The system also has two uninvadable steady states (A) in which two microbes (${N}_{1}\to {C}_{1}$ and ${N}_{2}\to {C}_{1}$) are limited by their nitrogen sources and (B) in which two other microbes are limited by their carbon sources (${C}_{1}\to {N}_{2}$ and ${C}_{2}\to {N}_{1}$). These two states have their regions of feasibility in the influx space. In order for these regions to overlap with each other, thereby resulting in bistability within the overlapping region, the yield parameters of species and nutrient supply rates have to satisfy the following conditions.
For the state (A), the conservation laws read as
The last two relations in (Equation S6) define microbe concentrations as ${B}_{11}={Y}_{11}^{(n)}\phantom{\rule{thinmathspace}{0ex}}(\phantom{\rule{thinmathspace}{0ex}}\frac{{\varphi}_{1}^{(n)}}{\delta}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\frac{\delta}{{\lambda}_{11}^{(n)}})$, ${B}_{22}={Y}_{22}^{(n)}\phantom{\rule{thinmathspace}{0ex}}(\phantom{\rule{thinmathspace}{0ex}}\frac{{\varphi}_{2}^{(n)}}{\delta}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\frac{\delta}{{\lambda}_{22}^{(n)}})$. Substituting these expressions into the first two relations in (Equation S6) and invoking the requirements $c}_{1}>\frac{\delta}{{\lambda}_{11}^{(c)}},\phantom{\rule{thickmathspace}{0ex}}{c}_{2}>\frac{\delta}{{\lambda}_{22}^{(c)}$ (guaranteeing that neither of two carbons limits microbes’ growth), in the high flux limit we obtain
The inequality conditions above are only natural given that the microbes use up their nitrogen fluxes very thoroughly in the state (A) while (at least in the highflux limit) they not getting even close to consuming all of their carbon supply rates .
The state (B) in the high flux limit will require different conditions, although obtained in a perfectly similar way:
Combining Equation S7 and Equation S8 one gets
Hence, for the carbon fluxes ratio, one would have
which implies that for multistability to be possible at least for some ratio of fluxes the yields have to satisfy the following inequality:
Note that the Equation S11 is both necessary and sufficient for bistability between (A) and (B) for some set of supply rates. Indeed, if Equation S11 is satisfied, ${\varphi}_{2}^{c}$ can be chosen arbitrarily (the only thing one would have to mind here is the highflux limit requirements), then ${\varphi}_{1}^{c}$ should be chosen in accordance with Equation S10, and any nitrogen fluxes satisfying Equation S9. All the procedures are legitimate whenever Equation S11 holds. Once chosen in the way described above, the point $({\varphi}_{1}^{(c)},{\varphi}_{2}^{(c)},{\varphi}_{1}^{(n)},{\varphi}_{2}^{(n)})$ of the influx space will make both (A) and (B) steady states feasible.
In a more general case of an arbitrary number of nutrients of each type, the conditions allowing for bistability or even multistability can be expressed by simple inequalities connecting yields and fluxes in combinations dictated by network topology of potentially bistable states in a very similar way to the simple case presented above. Thus the solution to the puzzle of why roughly half of all possible yield combinations has no multistability whatsoever becomes intuitively clear. Indeed, these yields and fluxes must come in ‘dimensionless’ combinations so that any inequality can be written as a function of only C:N stoichiometry ${S}_{ij}^{(C:N)}={\displaystyle \frac{{Y}_{ij}^{(n)}}{{Y}_{ij}^{(c)}}}$ for all microbial species present in any of the set of potentially multistable states. Note that should nitrogen and carbon yields exchange places for each of these microbes, the key inequality similar to Equation S11 would be reversed, thus prohibiting multistability where it was permitted and vice versa.
In the yield space, the proposed swap of carbon and nitrogen yields of all species ${Y}^{(c)}\to {Y}^{(n)},{Y}^{(n)}\to {Y}^{(c)}$ is a volume preserving transformation. This means that, for each set of potentially multistable states, the fraction of the yield space favouring multistability is exactly the same as the fraction prohibiting multistability. In a possible (albeit unlikely) scenario when the suggested permutation affects not only the potential multistability in question, but also some other multistable conditions, ‘turning off’ one multistability might in some cases ‘turn on’ others. This is why the ultimate empirical probability of multistability for a given combination of species’ yields might somewhat deviate from 1/2.
Appendix 5
Proof of the absence of multistability when all species have identical nutrient stoichiometries
For randomly selected nutrient supply rates, the chance for the system to end up in a given steady state is proportional to the volume of its region of feasibility. In a $(K+M)$ dimensional space, where each point corresponds to the set of K carbon and M nitrogen nutrient supply rates, regions of feasibility for some uninvadable states may intersect. In this case, the area they have in common defines the initial conditions at which multistability can occur. Naturally, for the effect to be observable, that shared area has to be $(K+M)$ dimensional, else its volume and, correspondingly, the chance to encounter such a situation in a real ecosystem would be negligibly small.
In the general case, when the microbes differ in terms of their stoichiometry, the steady state conditions dictate, for carbon supply rates:
where $\frac{{\varphi}_{i}^{(c)}}{\delta}$ is the abiotic value of the concentration of the nutrient ${c}_{i}$, $i$ ranging from 1 to $K$, ${B}_{ik}$, as usual, is the concentration of the microbe growing on ${c}_{i}$ and ${n}_{k}$ (equal to zero if the microbe is absent in the state), ${Y}_{ik}^{(c)}$ — the yield of the microbe with respect to its carbon source.
Similarly, for nitrogen supply rates one would have:
here ${\varphi}_{j}^{(n)}$ being the abiotic value of the concentration of the nutrient ${n}_{j}$, ${Y}_{lj}^{(n)}$ – the yield of the microbe using ${c}_{l}$ and ${n}_{j}$ with respect to its nitrogen source, $j$ ranging from 1 to $M$.
When multistability occurs, there exist at least two different sets of values of ${B}_{ij}$, ${c}_{l}$ and ${n}_{m}$ such that either of them can fit both Equation S12 and Equation S13 for the same values of nutrient supply rates. After marking one of them with the symbol ∼ for carbon sources one gets:
and for nitrogen sources:
In the general case, where the microbes differ in their yields, Equation S14 and Equation S15 taken together for all relevant $i$, $j$, describe a $(K+M)$–dimensional area. Indeed, the nutrient concentration on any occupied resource is fixed to $\delta /{\lambda}_{\alpha}^{(c,n)}$, where $\alpha $ is the species limited by this resource. Conversely, for any unoccupied resource nutrient concentration is free to vary as long as it is larger than $\delta /{\lambda}_{\beta}^{(c,n)}$, where $\beta $ is the species using this resource in a given state which has the smallest value of ${\lambda}^{(c,n)}$. In other words each unoccupied resource contributes one degree of freedom, while occupied resources do not. There is also a positive microbe concentration value for any occupied nutrient source. Thus, the total number of the unknowns is being equal to $K+M$ on the left side and $K+M$ on the right side of the system of the equations. The total number of the unknowns is thus $2K+2M$. The number of restrictions provided by Equation S14 and Equation S15 being $K+M$, the area defined by the system of the equations turns out to be $(K+M)$ dimensional.
If the stoichiometry (yields) of all species are equal to each other, without loss of generality one can assume them all to be equal to unity. Indeed, in this case the units of resource concentrations or bacterial abundances can be universally changed to match each other so that one unit of resource is used to make exactly one unit of bacterial abundance.
For identical stoichiometry of species there is an additional constraint on our $2K+2M$ variables. Indeed, after adding up all equations Equation S14 and then, separately, all equations Equation S15, and comparing the results of these additions we get
This reflects the fact that for identical stoichiometries both sides of Equation S16 are equal to the total biomass concentration ${\sum}_{ij}{B}_{ij}$ of all surviving microbes. Hence, if the two sets of concentrations satisfying both Equation S14 and Equation S15 exist, under the degenerate stoichiometry condition they must also satisfy an additional constraint
as the difference between the total concentration of carbon nutrients and that of nitrogen nutrients must be the same in both states by the virtue of Equation S16.
The additional constraint Equation S17 reduces the dimensionality of the multistability area to $K+M1$, so in the $(K+M)$ dimensional space of nutrient concentrations it has zero volume.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files. The code for both our simulations and statistical analysis can be downloaded from: https://github.com/ssm57/CandN (copy archived at https://github.com/elifesciencespublications/CandN).
References

Analysis of doublesubstrate limited growthBiotechnology and Bioengineering 20:183–202.https://doi.org/10.1002/bit.260200203

Fast unfolding of communities in large networksJournal of Statistical Mechanics: Theory and Experiment 2008:P10008.https://doi.org/10.1088/17425468/2008/10/P10008

Stability criteria for complex microbial communitiesNature Communications 9:2970.https://doi.org/10.1038/s4146701805308z

Nutrient limitation of soil microbial processes in tropical forestsEcological Monographs 88:4–21.https://doi.org/10.1002/ecm.1279

Global stability and multiple domains of attraction in ecological systemsThe American Naturalist 113:705–714.https://doi.org/10.1086/283427

MacArthur's consumerresource modelTheoretical Population Biology 37:26–38.https://doi.org/10.1016/00405809(90)90025Q

Von Liebig's law of the minimum and plankton ecology (1899–1991)Progress in Oceanography 33:347–386.https://doi.org/10.1016/00796611(94)900221

College admissions and the stability of marriageThe American Mathematical Monthly 69:9–15.https://doi.org/10.1080/00029890.1962.11989827

Experimental studies on the struggle for existence: I. mixed population of two species of yeastJournal of Experimental Biology 9:389–402.

Multistability and the origin of microbial community typesThe ISME Journal 11:2159–2166.https://doi.org/10.1038/ismej.2017.60

Feasibility and coexistence of large ecological communitiesNature Communications 8:14389.https://doi.org/10.1038/ncomms14389

Sundials: suite of nonlinear and differential/algebraic equation solversACM Transactions on Mathematical Software 31:363–396.

Resilience and stability of ecological systemsAnnual Review of Ecology and Systematics 4:1–23.https://doi.org/10.1146/annurev.es.04.110173.000245

Growth kinetics of suspended microbial cells: from singlesubstratecontrolled growth to mixedsubstrate kineticsMicrobiology and Molecular Biology Reviews : MMBR 62:646–666.

Tipping elements in the human intestinal ecosystemNature Communications 5:4344.https://doi.org/10.1038/ncomms5344

Species packing and competitive equilibrium for many speciesTheoretical Population Biology 1:1–11.https://doi.org/10.1016/00405809(70)900390

Elemental stoichiometry of Fungi and Bacteria strains from grassland leaf litterSoil Biology and Biochemistry 76:278–285.https://doi.org/10.1016/j.soilbio.2014.05.011

Experimental tests for an evolutionary tradeoff between growth rate and yield in E. coliThe American Naturalist 168:242–251.https://doi.org/10.1086/506527

Metabolic TradeOffs promote diversity in a model ecosystemPhysical Review Letters 118:028103.https://doi.org/10.1103/PhysRevLett.118.028103

The biological control of chemical factors in the environmentAmerican Scientist 46:205–221.

Microbial nitrogen limitation in the mammalian large intestineNature Microbiology 3:1441–1450.https://doi.org/10.1038/s4156401802677

Catastrophic regime shifts in ecosystems: linking theory to observationTrends in Ecology & Evolution 18:648–656.https://doi.org/10.1016/j.tree.2003.09.002

Coexistence of many species in random ecosystemsNature Ecology & Evolution 2:1237–1242.https://doi.org/10.1038/s4155901806036

Fundamentals of microbial community resistance and resilienceFrontiers in Microbiology 3:417.https://doi.org/10.3389/fmicb.2012.00417

Multiple stable points in natural communitiesThe American Naturalist 108:859–873.https://doi.org/10.1086/282961

Collective phase in resource competition in a highly diverse ecosystemPhysical Review Letters 118:048103.https://doi.org/10.1103/PhysRevLett.118.048103

Resource competition and community structureMonographs in Population Biology 17:1–296.

Onset of natural selection in populations of autocatalytic heteropolymersThe Journal of Chemical Physics 149:134901.https://doi.org/10.1063/1.5048488
Decision letter

Wenying ShouReviewing Editor; Fred Hutchinson Cancer Research Center, United States

Ian T BaldwinSenior Editor; Max Planck Institute for Chemical Ecology, Germany

Wenying ShouReviewer; Fred Hutchinson Cancer Research Center, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Thank you for submitting your article "Multistability and regime shifts in microbial communities explained by competition for essential nutrients" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Wenying Shou as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Ian Baldwin as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary
"Multistability and regime shifts in microbial communities explained by competition for essential nutrients" by Dubinkina et al. explores the presence of diverse, multiple stable states for diverse ecosystems. The model builds on the groups earlier work using results from game theory on the stable marriage problem to ask under what conditions it is possible to have multistability in diverse communities. The authors employ an idealized model of competition for essential nutrients to study patterns of multistability and diversity in microbial ecosystems. By computationally enumerating steady states and classifying them by stability, the authors establish that multistability is a typical feature of resource competition between specialized consumers. Specifically, they show that multistability emerges robustly when nutrient supplies are balanced and competitors have different C:N stoichiometries. The topic of how diversity is maintained in microbial ecosystems is likely to be of broad interest to the readership of eLife. While the current model is quite simplified compared to any real multispecies community, the authors do a good job of highlighting general lessons from the model that make testable predictions for natural and experimental systems. Overall, all three reviewers found the work interesting.
Essential revisions
Essential revisions fall under three categories:
1) Present your model more clearly;
2) Add mechanisms and intuitions (rather than simply describing your data);
3) Discuss generalizability of the results, including sensitivity to parameter choices.
To help you revise, we attach the reviews.
Reviewer #1:
My main suggestions are as following:
1) The authors define states by the presence/absence of species, and if present, which nutrient it would be limited in. Out of the theoretical possible states, the authors found all steady states by noting a correspondence of their problem and the stable marriage problem. This correspondence needs to be explained in an accessible form (perhaps in a figure).
2) Put more emphasis on going through the data and explain what they mean. For example, in Figure 1B: why do we get these states? For the twospecies 11, 22 cases, one could imagine a different state where 11 is limited for carbon, and 22 is limited for nitrogen. Why is this case not seen? Does this have something to do with your parameter choices? For Figure 3B, why does the likelihood of multistability peak at C:N supply rates similar to the average C:N stoichiometry? If you change the average C:N stoichiometry ratio by 10fold, does the peak change by 10fold? Seems that the authors were simply describing results rather than also giving intuition on why the results arose.
3) Subsection “Three criteria for stability of microbial communities” paragraph four: not clear why invasion does not depend on nutrient supply rates, given that nutrient supply rates could change growth rates.
4) Abstract: Multistability requires different stoichiometries of essential nutrients. How different does "different" have to be? It would be good to supplement with previous experimental data on these ratios and discuss.
Reviewer #2:
My main concern is really about how much the conclusions from this particular model used by the authors generalize to a more general setting. I think the assumptions of the model are quite stringent (and need to be in order to really use the technical apparatus of game theory/stable marriage problem they use):
Each species can use only one C and one N resource; a single metabolite cannot provide both C and N (even though we know for example bacteria can and do use amino acids as C sources).
As the authors point out, this gives rise to nongeneric behavior such as the fact:
"Unlike other consumer resource models, in our model the dynamic stability of a state with respect to species invasions does not depend on nutrient supply rates."
At the minimum, the authors have to argue what depends on these assumptions and what does not. Comparing the results to more standard Consumer Resource Models (CRMs) would go a long way at clarifying this. Despite this criticism, I want to emphasize that I think the paper explores and interesting question and gives rise to interesting hypothesis that are interesting and potentially testable.
I now give more detailed comments:
Technical Comments:
The section on MonteCarlo sampling is extremely unclear. What exactly is the distribution that is sampled? Is this done only in the highflux regime? What is assumed about the R matrix being invertible? Is this somehow an assumption that the Y are "generic"? It is very unclear to me why one cannot be feasible on a lowdimensional manifold. This is usually a generic situation in most CRMs so why is it excluded here. I would appreciate a much more technical discussion, perhaps with pseudo code so that I can better understand if and when the arguments generalize.
The equations in subsection “Conditions of multistability in the 2C x 2N x 4S ecosystem” have typos (extra factors of Y).
I did not understand in the section on lowerbound, why "double counting" can be discarded (either when # of C = # of N or otherwise).
It is very unclear how sensitive are the results to different "draws" of the λ and Y. There is a particular realization given in SI Tables but it would be nice to have a sense of the fluctuations rather than just for one network.
Another general way of modeling essential resources is to use a function of the from
where the $c}_{I\alpha$ are constants with the same units as $R}_{\alpha$ that encode ratios of resources required for building biomass (see Taillefumier et al., 2017). Recently, it was shown in arxiv:1901.09673 that such functions can also give rise to multistability using a Minimum Environmental Perturbation Principle based on measuring distances from the input flux vectors with appropriate distance metric. This suggests it should be possible to generically construct input fluxes into the system with multiple steadystates by finding input vectors that are equally distant from the intersection of all ZGNIs for a different, more traditional class of CRMs where species can eat many resources. This again raises questions about what generalizes from this simple but interesting model. It seems that this would be more in line with the generalization in equation in subsection “Extensions of the model” but it may qualitatively change the results.
Conceptual Comments:
I really enjoyed the idea of how by continuously changing the environments one cannot connect very different stable states. This reminds me of old arguments in evolutionary theory on neutralness in adaptation by Fontana and collaborators (PNAS 1996) where very different evolutionary minima could be connected even though they changed differently. This is just a thought/comment not a suggestion for investigation.
The idea that equalizing resources promotes diversity is now something that has been found in various CRMs (Posfai et al., 2017; Tikhonov et al; Marsland et al., 2019). The idea is that one gets neutral like dynamics. Is there a way to reconcile these ideas with the multistability? It is curious that both these dynamics like equalizing resources and I would like to understand what, if any, the authors think the connection is.
I am really confused as stated above how these results depend on the available resource species pool. In particular, what if I changed the statistics of things were generated, yields, etc.? I think it would be nice to at least understand these things at least semiquantitatively.
Reviewer #3:
1) Existing literature on competition for essential nutrients.
The authors state that "while reallife ecosystems driven by competition for multiple essential nutrients have been studied experimentally Fanin et al., 2016; Browning et al., 2017; Camenzind et al., 2018, the resourceexplicit models capturing this type of growth are not so well developed beyond the foundational work by Tilman, 1982." Although this is generally true, it would be helpful to reference the work that does exist, to orient the reader as to whether the current results are surprising. Huisman's paper (Nature 402 1999) is relevant, as is the plankton literature more generally.
2) Assumptions and parameters of the model.
The assumption of complete specialization to a single carbon and nitrogen resource is very limiting. More care should be taken to justify it in the context of real ecosystems. Would the results hold qualitatively if this assumption were relaxed? Why or why not? Furthermore, each model seems to feature only one realization of the growth coefficient λs. The authors don't need to redo everything for new parameters, but they should confirm that their results are not specific to these parameters. Finally, in the current model a species' λs and Ys (growth yields) are uncorrelated. However, it would be more biologically plausible if these were subject to tradeoffs for each species (e.g. high growth rates correspond to low yields). It would be informative if the authors would check the case where λY=constant to confirm that their results do not depend on λ and Y being uncorrelated.
3) Explanation of Figure 2EF
The authors should try to provide an intuitive, mechanistic understanding of these computational results. Why does state 2 have the largest feasibility volume, and why are regime shifts possible between some states and not others? If I understood the simple case of 2e, I would have a better appreciation of 2f. A discussion of the ecological implications would also be welcome here.
4) 1D manifold in Figure 3C
In the same spirit as point #3, the authors should try to explain this result and its implications. This is certainly an interesting observation, but as of now it's asserted without any mechanistic understanding.
5) Competitive rank and diversity
The "competitive rank" metric is unintuitive, at least as presented. The authors should explain more clearly why ranks run from 16 and why they're a species property independent of the environment. The surprising result that a rank 5.5 species is uninvadable 20% of the time should be explained. What is a specific scenario in which this occurs, and why?
6) Volume in abundance space
The authors focus on the feasibility volume of a state in the space of nutrients. However, it is also natural to consider the volume of the basins of attraction in the space of species abundances for a particular state in a given environment, when multiple stable states are possible. Does the volume of the basin of attraction for a state relate to its feasibility volume?
https://doi.org/10.7554/eLife.49720.sa1Author response
Reviewer #1:
My main suggestions are as following:
1) The authors define states by the presence/absence of species, and if present, which nutrient it would be limited in. Out of the theoretical possible states, the authors found all steady states by noting a correspondence of their problem and the stable marriage problem. This correspondence needs to be explained in an accessible form (perhaps in a figure).
We thank the reviewer for this suggestion. To address this we added a paragraph with additional details clarifying the general ideas behind the correspondence between our model and the stable marriage problem (see subsection “Microbial community growing on two types of essential nutrients represented by
85 multiple metabolites” paragraph six and “Three criteria for stability of microbial communities” paragraph two). These are further explained in the Materials and methods section "Identification of all states and classification of them as invadable or uninvadable”. We also believe that knowing the exact correspondence to the stable marriage problem is not essential for understanding of our main results and may even be distracting to our readers.
We have thought long and hard if this correspondence could be illustrated in an additional figure in the main text of the article and (after long deliberations) decided to leave it for supplementary materials (see Appendix 3figure 1). Unlike a straightforward mapping to the stable marriage problem outlined in our previous paper (Goyal et al., 2018), the relation to the stable marriage problem described in this study, while mathematically exact, is a bit more formal. Indeed, as explained in the section "Stable matching approach for identification and classification of steady states" of the supplementary notes, "marriages" in the C:N model are unidirectional and represented by directed arrows between C and N sources. It would be confusing to readers to work with two different graphical representations of the same problem (compare Figure 1A and Appendix 3figure 1 in the current version of the manuscript).
2) Put more emphasis on going through the data and explain what they mean. For example, in Figure 1B: why do we get these states? For the twospecies 11, 22 case, one could imagine a different state where 11 is limited for carbon, and 22 is limited for nitrogen. Why is this case not seen? Does this have something to do with your parameter choices?
We thank the reviewer for bringing this up and helping us to clarify our results. Indeed, the list of uninvadable states illustrated in Figure 1B is determined by our specific choice of model parameters λ quantifying competitiveness of individual microbial species (see Supplementary file 1). For this specific case, the state where the species 11 is limited by carbon, and the species 22 – by nitrogen is always invadable by the species 21 and hence not listed among uninvadable states in Figure 1B. We added the explanation of another case to the main text to provide readers with a specific example and additional understanding of how the relative ranking of lambdas defines the invadability of a state (paragraph two subsection “Three criteria for stability of microbial communities”).
For Figure 3B, why does the likelihood of multistability peak at C:N supply rates similar to the average C:N stoichiometry? If you change the average C:N stoichiometry ratio by 10fold, does the peak change by 10fold? Seems that the authors were simply describing results rather than also giving intuition on why the results arose.
The short answer is yes: for a 10:1 stoichiometry the multistability peak in Figure 3B will shift to a 10fold higher $\varphi}^{(c)}/{\varphi}^{(n)$ ratio. In general, the C:N stoichiometry averaged across all species can be renormalized to the 1:1 ratio used in our study by redefining the nutrient concentrations and supply rates. For example if the average stoichiometry is 10:1 one can measure all carbon concentrations c_{i} in 10fold larger units than nitrogen concentrations n_{j} thereby renormalizing the stoichiometry to 1:1. We have added a brief explanation of this in paragraph three of subsection “Patterns of multistability”
3) Subsection “Three criteria for stability of microbial communities” paragraph four: not clear why invasion does not depend on nutrient supply rates, given that nutrient supply rates could change growth rates.
Indeed, nutrient supply rates change nutrient concentrations, which in turn, change the growth rates of species. However, they do so proportionally for all species using the same nutrient. Throughout the manuscript and in all our numerical simulations we are working in the high nutrient supply regime. In this regime, all individual microbes would be able to grow as isolates. Hence, when multiple species are introduced, their survival is solely due to the competition between them. In this case, the invadability of a state is fully determined by the rankorder of lambdas of the invading species and all species present in the state. In the opposite regime, where at least some of the supply rates are low the set of uninvadable states may expand beyond what we predict. However, our uninvadable states would remain uninvadable even in the low nutrient supply regime.
To clarify this argument we have rephrased the text. It now reads:
"Note that, in the regime of our model, where the supply of all nutrients is high, i.e. $\varphi}^{(c,n)}\gg {\delta}^{2}/{\lambda}_{\alpha}^{(c,n)$, invadability of individual states does not depend on supply rates. Indeed, in this regime the outcome of an attempted invasion is fully determined by the competition between species, which in turn depends only on the rankorder of λ's of the invading species relative to the species currently present in the ecosystem“. We also clarified it in the respective subsection in Materials and methods.
4) Abstract: Multistability requires different stoichiometries of essential nutrients. How different does "different" have to be? It would be good to supplement with previous experimental data on these ratios and discuss.
Multistability is possible for arbitrarily small variation of stoichiometry, but in this case, it would be confined to a very narrow region of nutrient supply rates. We have previously illustrated it in what was Figure S3 of the original submission. Inspired by this comment and comments of other reviewers we decided to move this figure to the main text of the manuscript (see Figure 5 in the revised submission). It shows that the smaller is the variance of stoichiometries (the xaxis) – the less likely one is to encounter multistability for randomly selected combination of nutrient supply rates (the yaxis of Figure 5). We have added a discussion of this along with comparison with experimental data in subsection “Multistability requires diverse species stoichiometry”.
We also showed where the reallife soil microbial ecosystem stoichiometry variation lies on this plot.
We also added the experimentally determined stoichiometry variance for real ecosystems (mammalian gut, phytoplankton, and soil). This suggests that multistability might not be very rare (up to 10% likelihood for soil). We also wrote supplementary note (Appendix 5) that proves that multistability is impossible if all species have identical nutrient stoichiometries.
Reviewer #2:
1) My main concern is really about how much the conclusions from this particular model used by the authors generalize to a more general setting. I think the assumptions of the model are quite stringent (and need to be in order to really use the technical apparatus of game theory/stable marriage problem they use):
Each species can use only one C and one N resource; a single metabolite cannot provide both C and N (even though we know for example bacteria can and do use amino acids as C sources).
As the authors point out, this gives rise to nongeneric behavior such as the fact:
"Unlike other consumer resource models, in our model the dynamic stability of a state with respect to species invasions does not depend on nutrient supply rates."
At the minimum, the authors have to argue what depends on these assumptions and what does not. Comparing the results to more standard Consumer Resource Models (CRMs) would go a long way at clarifying this. Despite this criticism, I want to emphasize that I think the paper explores and interesting question and gives rise to interesting hypothesis that are interesting and potentially testable.
We agree with the reviewer that the assumptions of our model are both stringent and crucial for the exact mapping to the Stable Marriage Problem. However, the conceptual framework in which we reduce multistability between a pair of states to a problem of overlaps of their feasibility regions is general and applicable to a broad variety of Consumer Resource Models. The extent to which our findings hold in other CRM models is the subject of ongoing research in our lab. We added a paragraph comparing one of the key results in our model to that in other CRMs (subsection “Multistability requires diverse species stoichiometry”). We also emphasized the universality (or lack thereof) of our results in other CRMs in subsection “Interplay between diversity and stability in ecosystems with multiple essential nutrients”. Finally, we cited arxiv:1901.09673, which has a number of additional insights into multistability in CRMs.
In the end of our manuscript we propose two generalizations of our model that can be investigated numerically: the one in which $g(C,N)=min\left(\sum _{i}{\lambda}_{i}^{(c)}{C}_{i},\sum _{i}{\lambda}_{j}^{(c)}{N}_{j}\right)$ (a hybrid of our model with the MacArthur model) and another one in which microbes diauxically shift between resources (a hybrid of our model with our previously published model (Goyal et al., 2018) also mapped onto the Stable Marriage Problem). We expect our mathematical techniques to be at least partially extendable to both of these models.
I now give more detailed comments:
Technical Comments:
2) The section on MonteCarlo sampling is extremely unclear. What exactly is the distribution that is sampled? Is this done only in the highflux regime?
We MonteCarlo (MC) sampled our flux space randomly (and uniformly in linear space) between 10 and 1000. For each supply combination point, all nutrients were sampled independently as described in Materials and methods. The lower bound was selected in such a way that all our MC sampling was done for the high nutrient supply regime (indeed our λ’s are >10 and δ=1 giving the lower bound of the high flux regime <0.1 well below the lowest nutrient flux of 10. It is now described in more detail in the Materials and methods section.
3) What is assumed about the R matrix being invertible? Is this somehow an assumption that the Y are "generic"?
Yes, to instantiate our model we randomly generated Y, which all but guaranteed that the matrices $R}_{p$ in Equations 5,6 are invertible. Thus given nutrient supply rates we can always calculate bacterial abundances and concentrations of all nutrients, which are not limiting to any microbes. The problem is made nontrivial by the feasibility constraint that all of these have to be positive.
4) It is very unclear to me why one cannot be feasible on a lowdimensional manifold. This is usually a generic situation in most CRMs so why is it excluded here. I would appreciate a much more technical discussion, perhaps with pseudo code so that I can better understand if and when the arguments generalize.
One common strategy in CRMs is to work with nutrient supply rates artificially constrained to a lowdimensional manifold (e.g. by adding a constraint that $\sum _{i}{\varphi}_{i}={\varphi}_{total}$). Such constraints would not change our results except by reducing all of the dimensions used in our study by 1. For example, it would turn K^{+}M dimensional space of $\varphi$ into K^{+}M1 dimensional space of constrained $\varphi$. At the same time the dimensionality of the boundaries between states would go down from K^{+}M1 to K^{+}M2. Hence, a shared boundary between states still would not imply multistability.
Another argument on why we ignore states feasible only on a lowerdimensional manifold is as follows.
We evaluate the feasibility of a state separately for each set of nutrient supply rates $\varphi$. If a state is feasible only on a lowdimensional subset of $\varphi$, any small perturbation of these will take the system out of this state. Such a state would be nearly impossible to detect in our Monte Carlo simulations. It is also unlikely to be experimentally observed in real microbial ecosystems unless there is a special reason why nutrient supply rates would satisfy this constraint. This treatment of feasibility volumes is standard (see e.g. Grilli et al., 2016 for GLV model or Butler and O'Dwyer, 2018 for the MacArthur model).
The equations in subsection “Conditions of multistability in the 2C x 2N x 4S ecosystem” have typos (extra factors of Y).
These have been corrected.
5) I did not understand in the section on lowerbound, why "double counting" can be discarded (either when # of C = # of N or otherwise).
We tried to better explain our proof that double counting is impossible for stable states identified by the GaleShapley men (nitrogen) proposing algorithm. These are the only states we use in our calculations of the lower bound on the number of uninvadable states. The modified explanation is in Appendix 3.
6) It is very unclear how sensitive are the results to different "draws" of the λ and Y. There is a particular realization given in SI Tables but it would be nice to have a sense of the fluctuations rather than just for one network.
We agree with the reviewer that the question of how parameters λ and Y affect the set of states in our model is an interesting topic for future investigation. During our preparation of this manuscript, we experimented with different sets of lambdas and Y's. We noticed that while the exact set of dynamically stable and uninvadable states depends on our choice of parameters, their number and statistical properties remain the same. We also noticed that the choice of Y’s has a much stronger influence on multistability, which is the main focus of this study. This influence was illustrated in what was the Figure 3S (now moved as Figure 5 in the main text) with different draws of Y and the same set of λ's in the 2C x 2N x 4S model. Similar parameter sensitivity analysis of the 6C x 6N x 36S model was prohibitively computationally expensive.
7) Another general way of modeling essential resources is to use a function of the from
$g=(\sum _{\alpha}{c}_{i\alpha}/{R}_{\alpha}{)}^{1}$ where the $c}_{I\alpha$ are constants with the same units as $R}_{\alpha$ that encode ratios of resources required for building biomass (see Taillefumier et al., 2017). Recently, it was shown in arxiv:1901.09673 that such functions can also give rise to multistability using a Minimum Environmental Perturbation Principle based on measuring distances from the input flux vectors with appropriate distance metric. This suggests it should be possible to generically construct input fluxes into the system with multiple steadystates by finding input vectors that are equally distant from the intersection of all ZGNIs for a different, more traditional class of CRMs where species can eat many resources. This again raises questions about what generalizes from this simple but interesting model. It seems that this would be more in line with the generalization in equation in subsection “Extensions of the model” but it may qualitatively change the results.
We have now cited arxiv:1901.09673 and the Minimum Environmental Perturbation Principle (MEMP), which may provide new insights on general question of multistability in CRMs. We do not fully understand the connection between our methods and MEMP. For example, our model does not have multistability if the yields of all microbes on a given nutrient are equal to each other (see proof in Appendix 5). That is to say, in this case the feasibility regions of different uninvadable states do not overlap with each other, except on a lowdimensional manifold of nutrient supply rates. The model from arxiv:1901.09673 with $g=(\sum _{\alpha}{c}_{i\alpha}/{R}_{\alpha}{)}^{1}$ can demonstrably be multistable (see their Figure 5C). However, it is unclear whether this multistability is limited to a lowdimensional manifold composed of "input vectors that are equally distant from the intersection of all ZGNIs". If it does, our results agree with arxiv:1901.09673. Overall, we believe that these questions could be fruitfully explored in future extensions of our model and are out of scope for this study.
Conceptual Comments:
8) I really enjoyed the idea of how by continuously changing the environments one cannot connect very different stable states. This reminds me of old arguments in evolutionary theory on neutralness in adaptation by Fontana and collaborators (PNAS 1996) where very different evolutionary minima could be connected even though they changed differently. This is just a thought/comment not a suggestion for investigation.
We thank the reviewer for this positive feedback. We thought about the analogy between addition of species in our model and addition of mutations allowed by natural selection in models similar to Fontana, 1996. We explained this observation in more detail in our manuscript (see subsection “Control of microbial ecosystems exhibiting multistability and regime shifts”). We also have partially implemented this analogy in a forthcoming manuscript where we explore possible changes in the community state along different trajectories in nutrient supply rates and species additions.
9) The idea that equalizing resources promotes diversity is now something that has been found in various CRMs (Posfai et al., 2017; Tikhonov et al; Marsland et al., 2019). The idea is that one gets neutral like dynamics. Is there a way to reconcile these ideas with the multistability? It is curious that both these dynamics like equalizing resources and I would like to understand what, if any, the authors think the connection is.
We have rather good conceptual understanding of this issue. In our model, multistability and species diversity are two sides of the same coin. In the region where nutrient supplies are equalized (wellbalanced region in our notation), there is an excess of states with high diversity of species (see our Figure 4A). Some of these states are dynamically stable – in this case, high diversity is maintained. Dynamically unstable states on the other hand break up into bistable pairs of lowerdiversity states.
It is easy to understand this relation using the example of Tilman model with a single C and a single N resource. In this case, the wellbalanced region is located near the diagonal of the CN plane. For such nutrient supplies one has either the coexistence between two species (a highdiversity state) or bistability between two lowdiversity ones in which either one or the other species survive. The choice between these two scenarios is dictated by the combination of yields and competitiveness of two species, which determine dynamical stability of the twospecies steady state. Similar arguments apply to our model with many C and many N resources. However, since the number of states that could be dynamically unstable is large, there is no mutual exclusivity between multistability and high species diversity in the wellbalanced region. We have cited this connection between our model and other CRMs in subsection “Interplay between diversity and stability in ecosystems with multiple essential
nutrients”.
10) I am really confused as stated above how these results depend on the available resource species pool. In particular, what if I changed the statistics of things were generated, yields, etc.? I think it would be nice to at least understand these things at least semiquantitatively.
We plan to explore some of these topics in future papers. Mapping to the Stable Marriage Problem proves that the statistics (meaning PDFs) of lambdas does not matter for our model. The relative ranks of lambdas fully determine the set of steady states in our model and their invadability. However, as discussed in the paper and shown in Figure 5 (formerly Figure S3) the standard deviation of yields has a strong impact on the extent of multistability in our model. Moreover, we suspect that tradeoffs or correlations between λ^{(c)} and λ^{(n)} or λ’s and yields might affect the number of stable states in the model. Indeed, we have demonstrated this in our previous work (Goyal et al., 2018) in a different model also based on the stable marriage problem. We now speculate about potential consequences of different tradeoffs in a dedicated section of the Discussion in subsection “Multistability and the total number of states are affected by tradeoffs”.
Reviewer #3:
1) Existing literature on competition for essential nutrients.
The authors state that "while reallife ecosystems driven by competition for multiple essential nutrients have been studied experimentally Fanin et al., 2016; Browning et al., 2017; Camenzind et al., 2018, the resourceexplicit models capturing this type of growth are not so well developed beyond the foundational work by Tilman, 1982." Although this is generally true, it would be helpful to reference the work that does exist, to orient the reader as to whether the current results are surprising. Huisman's paper (Nature 402 1999) is relevant, as is the plankton literature more generally.
We expanded our literature list with several experimental papers quantifying the variability of nutrient stoichiometry in several real ecosystems (phytoplankton in oceans, microbes in soil, and mammalian gut). These references helped us to put a number on likelihood of multistability in each of these systems, which is now shown in Figure 5 (see subsection “Multistability requires diverse species stoichiometry” for the discussion of this figure). We also added an important reference to the Huisman's paper, which we previously overlooked (see paragraph two of the Introduction).
2) Assumptions and parameters of the model.
The assumption of complete specialization to a single carbon and nitrogen resource is very limiting. More care should be taken to justify it in the context of real ecosystems. Would the results hold qualitatively if this assumption were relaxed? Why or why not?
The mathematical mapping to the stable marriage problem depends on species' complete specialization to use a single carbon and a single nitrogen resource. However, the conceptual framework in which the multistability between uninvadable states is connected to overlaps of their feasibility regions of nutrient supply rates is completely general. We expect the following findings of our model to hold in a broad variety of consumer resource models and real ecosystems: 1) positive correlation between the likelihood of finding multistable nutrient supply combinations and the variability of nutrient stoichiometry (Figure 5 and discussion in subsection “Multistability requires diverse species stoichiometry”); 2) positive correlation between multistability (Figure 3C), species diversity (Figure 4A) and balanced nutrient supply (see the discussion in response to point 9 by reviewer #2 and subsection “Interplay between diversity and stability in ecosystems with multiple essential nutrients”). The extent to which these and other findings can be extended to other CRM models is the subject of ongoing research in our lab as outlined in Discussion.
Furthermore, each model seems to feature only one realization of the growth coefficient λs. The authors don't need to redo everything for new parameters, but they should confirm that their results are not specific to these parameters.
Mapping to the Stable Marriage Problem proves that the statistics of lambdas does not at all matter for our model. The relative ranks of lambdas fully determine the set of steady states in our model and their invadability in the high nutrient supply regime (see subsection “Three criteria for stability of microbial communities” for explanation). We suspect that tradeoffs or correlations between different lambdas might affect the number of stable states in the model (as discussed in subsection “Multistability and the total number of states are affected by tradeoffs”). Indeed, we have demonstrated this in our previous work (Goyal et al., 2018) in a different model also based on stable marriage problem.
During the preparation of our manuscript we experimented with different sets of λ's and Y's. We noticed that while the exact set of dynamically stable and uninvadable states depends on our choice of parameters, their statistical properties remain the same. We also noticed that the choice of Y’s has a much stronger influence on multistability, which is the main focus of this study. This influence was illustrated in what was the Figure 3S (now moved as Figure 5 in the main text) with different draws of Y and the same set of λ's in the 2C x 2N x 4S model (see subsection “Multistability requires diverse species stoichiometry” for the discussion of this figure). Similar parameter sensitivity analysis of the 6C x 6N x 36S model was prohibitively computationally expensive.
Finally, in the current model a species' λs and Ys (growth yields) are uncorrelated. However, it would be more biologically plausible if these were subject to tradeoffs for each species (e.g. high growth rates correspond to low yields). It would be informative if the authors would check the case where λY=constant to confirm that their results do not depend on λ and Y being uncorrelated.
We thank the reviewer for this comment as it lead us to an important discussion. While we do not study any correlated sets of parameters in the current variant of our model, we added a subsection in the Discussion section that summarizes our ideas on this topic (see subsection “Multistability could be affected by tradeoffs”).
3) Explanation of Figure 2EF
The authors should try to provide an intuitive, mechanistic understanding of these computational results. Why does state 2 have the largest feasibility volume, and why are regime shifts possible between some states and not others? If I understood the simple case of 2e, I would have a better appreciation of 2f. A discussion of the ecological implications would also be welcome here.
We thank the reviewer for this comment, as it helped us to improve the manuscript. The states #1 and #2 have the lowest species diversity (2 species as opposed to 3 species in states #36), thus they can be realized for a larger number of nutrient supply combinations (this diversity/feasibility dependence is illustrated in Figure 4C and discussed in the manuscript). The fact that the volume of the state #2 is larger than that of the state #1 is purely due to chance and is determined by our specific parameter choice.
We also added a new type of edges (marked blue) to Figure 2E, which mark continuous, smooth transitions between states without discontinuous regime shifts. These appear due to shared boundaries (instead of overlaps) between the pairs of states. Using the network shown in Figure 2E we developed some mechanistic intuition behind which pairs of states are connected by a bistable regime shifts, which have a shared boundary, and which have no direct transition at all. As one can see, the driver of all bistable switches shown in Figure 2E are the mutual exclusivity of pairs of keystone species (1,1), (2,2) and (1, 2), (2,1) correspondingly. Once one of these pairs is established, it cannot be invaded by another pair. The state #7 with all 4 species present is dynamically unstable and for some nutrient supply rates it spontaneously collapses to one of two states with different keystone species giving rise to bistability. Conversely, a shared boundary is possible for some of the states sharing the pair of keystone species.
Using the example of our 2C x 2N x 4S model, we can now conjecture the following intuitive view of regime shifts in more complex models including our 6C x6N x 36S model. Every ecosystem has a relatively small number of "drivers of bistability’’ generated by certain dynamically unstable states (like state #7). Each of these states can collapse into mutually exclusive pairs of stable states each with its own set of "keystone" species. In addition to keystone species the states could also include ‘’satellite’’ species that do not generally affect the bistable switch between states.
Our detailed explanations of these effects were added to the manuscript (see Results and Discussion sections).
4) 1D manifold in Figure 3C
In the same spirit as point #3, the authors should try to explain this result and its implications. This is certainly an interesting observation, but as of now it's asserted without any mechanistic understanding.
We currently do not have a mechanistic interpretation of these results. We suspect that it is a consequence of some special type of Lyapunov function in our problem (see the Discussion) but we have not found this function and thus cannot answer this question. We added the following text describing possible ecological implications of our results in real microbial ecosystems:
"It is tempting to conjecture that some variant of our model may explain similar arrangements of clusters of microbial abundances commonly seen in PCA plots of real ecosystems. If this is the case, the gaps between neighboring clusters would correspond to dynamically unstable states of the ecosystem, which may be experimentally observable as long transients in species composition."
5) Competitive rank and diversity
The "competitive rank" metric is unintuitive, at least as presented. The authors should explain more clearly why ranks run from 16 and why they're a species property independent of the environment.
We clarified our definition of the species competitiveness rank in the manuscript (see subsection “Patterns of diversity and structural stability of states”). The competitive rank of a species for carbon (nitrogen) is given by the rank order of its λ^{(c)} (λ^{(n)}) among all species using the same carbon (nitrogen) source. Since in our example there are L=6 species using each carbon (nitrogen) resource, the rank varies between 1 for the most competitive species in our pool and 6 for the least competitive species. The average competitiveness rank is the mean of species' carbon and nitrogen ranks. Since the growth rate of a species limited by a resource c is given by λ^{(c)}*c, with c – being its environmental concentration, the species with a higher λ^{(c)} outcompetes the species with a lower λ^{(c)} for the resource limiting to both of them.
The surprising result that a rank 5.5 species is uninvadable 20% of the time should be explained. What is a specific scenario in which this occurs, and why?
The factors deciding species prevalence are complex. In fact, only one of the 3 species with the mean competitiveness rank of 5.5 has high prevalence. It is decided by the interplay between competitiveness of all species in our pool. This particular species using c_{1} and n_{1} by chance was present as a ``satellite'' species in a state (or several states) with high prevalence. The goal of Figure 4B was to illustrate that the outcome of the competition in our model highly dependent on the composition of the entire species pool. In contrast, in a model with only one type of resource (say C), the species with the rank 1 would always win and all species with lower ranks would have zero prevalence.
6) Volume in abundance space
The authors focus on the feasibility volume of a state in the space of nutrients. However, it is also natural to consider the volume of the basins of attraction in the space of species abundances for a particular state in a given environment, when multiple stable states are possible. Does the volume of the basin of attraction for a state relate to its feasibility volume?
There is a close relationship between state's feasibility volume in the space of nutrients and its basin of attraction in the space of species abundances. They are identical up to the determinant of the matrix R_{p} in Equation 5 in Materials and methods. It could be understood using the specific example in Equation 6 describing the state #5 in our 2C x 2N x 4S model. To calculate the feasibility volume of this state we search for all vectors of ${\varphi}_{i}\text{}({\varphi}_{i}{\varphi}_{max})$ for which it is feasible, that is to say, when all values of B and n in the left hand side of this equation are positive. To calculate the basin of attraction of species abundances, one goes in the opposite direction and scans all positive values of B_{11}, B_{12}, B_{22}, and n_{2} between, say 0 and some upper bound, and constrain them by the condition $\varphi}_{i}<{\varphi}_{max$ for all nutrients. The elementary linear algebra dictates that the volumes of these two searches are connected via det(R_{p}). In our manuscript we choose to use the feasibility volume since 1) it is more relevant to possible experiments, where it is easy to control nutrient supply rates but hard – bacterial abundances; 2) it has been previously studied in the literature (Grilli et al., 2016, Butler and O'Dwyer, 2018) and even given a special name "structural stability".
https://doi.org/10.7554/eLife.49720.sa2Article and author information
Author details
Funding
No external funding was received for this work.
Acknowledgements
Part of this work has been carried out at the University of Padova, Italy, in August 2018, during a scientific visit by one of us (SM). The authors wish to thank Prof. James O’Dwyer, Prof. Seppe Kuehn and Akshit Goyal for a critical reading of an earlier version of the manuscript. The authors also thank reviewers and editors for their comments, that helped to improve the manuscript.
Senior Editor
 Ian T Baldwin, Max Planck Institute for Chemical Ecology, Germany
Reviewing Editor
 Wenying Shou, Fred Hutchinson Cancer Research Center, United States
Reviewer
 Wenying Shou, Fred Hutchinson Cancer Research Center, United States
Publication history
 Received: June 27, 2019
 Accepted: October 1, 2019
 Version of Record published: November 22, 2019 (version 1)
Copyright
© 2019, Dubinkina et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,330
 Page views

 194
 Downloads

 8
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.