The ratio of single to co-colonization is key to complexity in interacting systems with multiple strains.

The high number and diversity of microbial strains circulating in host populations have motivated extensive research on the mechanisms that maintain biodiversity. However, much of this work focuses on strain-specific and cross-immunity interactions. Another less explored mode of pairwise interaction is via altered susceptibilities to co-colonization in hosts already colonized by one strain. Diversity in such interaction coefficients enables strains to create dynamically their niches for growth and persistence, and "engineer" their common environment. How such a network of interactions with others mediates collective coexistence remains puzzling analytically and computationally difficult to simulate. Furthermore, the gradients modulating stability-complexity regimes in such multi-player endemic systems remain poorly understood. In a recent study (Madec & Gjini, Bulletin of Mathematical Biology, 82), we obtained an analytic representation for N-type coexistence in an SIS epidemiological model with co-colonization. We mapped multi-strain dynamics to a replicator equation using timescale separation. Here, we examine what drives coexistence regimes in such co-colonization system. We find the ratio of single to co-colonization, µ, critically determines the type of equilibrium and number of coexisting strains, and encodes a trade-off between overall transmission intensity R 0 and mean interaction coefficient in strain space, k. Preserving a given coexistence regime, under fixed trait variation, requires balancing between higher mean competition in favorable environments, and higher cooperation in harsher environments, and is consistent with the stress gradient hypothesis. Multi-strain coexistence tends to steady-state attractors for small µ, whereas as µ increases, dynamics tend to more complex attractors. Following strain frequencies, evolutionary dynamics in the system also display contrasting patterns with µ, interpolating between multi-stable and fluctuating selection for cooperation and mean invasion fitness, in the two extremes. This co-colonization framework could be applied more generally, to study invariant principles in collective coexistence, and to quantify how critical shifts in community dynamics get potentiated by mean-field and environmental gradients.

INTRODUCTION

Rich ecosystems comprise many species interacting together in a myriad of ways and on multiple temporal and spatial scales. Understanding the scope and consequences of such interactions has been the focus of countless theoretical ecology studies, starting with the seminal work by Lotka (1926) and Volterra (1926) on mathematical models of the population dynamics of interacting species. This model has been later extended and sophisticated by many other theoretical studies (May 1972; Pascual et al., 2006), and is currently extensively used to characterize interaction networks in empirical microbiome communities (Bucci et al., 2016; Stein et al., 2013). Theoretically, a crucial question has been to study stability and coexistence patterns in such Lotka–Volterra multi‐species communities, analyzing both structured ecological networks and random networks (Serván et al., 2018; Song & Saavedra, 2018). Modeling efforts seek to understand organizing principles for species composition, including the balance between competition and cooperation (Mougi & Kondoh, 2012).

Overall, analysis of such models with arbitrarily high dimensionality has been and continues to remain difficult. In particular, beyond the complexity‐stability debate which represents a major force in ecology (Landi et al., 2018; May 1972; McCann, 2000), many studies are increasingly addressing the problem of deriving collective dynamics from pairwise outcomes between species, which constitutes another major challenge, both at an analytical (Levine et al., 2017; Momeni et al., 2017) and empirical (Friedman et al., 2017) level.

The challenge of high‐dimensionality in ecological microbial networks parallels a similar challenge in the epidemiology of polymorphic pathogen systems, where understanding the mechanisms and forces that maintain diversity among interacting strains, is also an area of active research (Cobey & Lipsitch, 2012; Gupta & Anderson, 1999; Lipsitch et al., 2009; Wearing & Rohani, 2006). While it is well recognized that population patterns of infection are to a large extent determined by susceptibility to infection, most multi‐strain SIR models, inspired from influenza, dengue, and malaria parasites, have focused on cross‐immunity between strains as driver of population structure (Gog & Grenfell, 2002; Gomes et al., 2002; Gupta et al., 1998; Lin et al., 1999). Yet, other factors, besides persistent host immunity, may make strains compete or cooperate with each other, and it remains unclear which environmental variables also contribute to their epidemiologic fitness.

Here, we bridge between multi‐species ecology and multi‐strain epidemiology, revisiting coexistence and diversity in a new context. We explore another mode of strain interactions, namely altered susceptibilities to coinfection, whereby N strains compete in SIS endemic scenarios of no persistent immunity and no virulence. Coinfection models with up to two strains have described vulnerability to coinfection with a single parameter (Alizon et al., 2013; Davies et al., 2019; Gaivão et al., 2017; van Baalen & Sabelis, 1995), two coefficients (Lipsitch, 1997) or four coefficients (Gjini et al., 2016) depending on model structure and aims, but very few analytical investigations have been done for a larger number of interacting strains (Adler & Brunet, 1991), recognizing the difficulties of including within‐and between‐strain details for such coefficients (Mosquera & Adler, 1998). Moreover, analytic solutions for strain frequency dynamics in coinfection models remain rare, due to nonlinearities even for N = 2.

In a recent co‐colonization (coinfection) SIS model framework, with N‐strains, we have simplified the complex ecology embedded in N 2 epidemiological variables (Madec & Gjini, 2020). Using timescale separation, we obtained a model reduction from the matrix of pairwise coinfection vulnerabilities between strains. This coincides with a special replicator equation (Cressman & Tao, 2014; Hofbauer & Sigmund, 2003) by which we can predict explicitly multi‐strain frequency evolution. This N‐dimensional model reduction makes the entire epidemiology more accessible to analysis, and relates emergent collective dynamics to the ensemble of pairwise competitive outcomes, not only qualitatively but moreover in an explicit quantitative manner.

In the present article, we harness the simplicity of this co‐colonization model framework (Madec & Gjini, 2020) to investigate coexistence, stability, and evolution of such multi‐strain systems with variable co‐colonization susceptibility coefficients among strains. We start by studying the behavior of the system for different global variables such as total transmission intensity R 0 and mean interaction coefficient in the pool of available strains k. We then study coexistence through random co‐colonization interactions, where the matrix coefficients are drawn from fixed distributions, and can range from competitive to cooperative links. We ask what is the number of strains that can coexist when starting from a pool of N strains, and in which diversity–stability configuration. We uncover rich transient and asymptotic behavior of such systems, where steady states, limit cycles, multi‐stability, and chaotic attractors are possible.

We find that the ratio of single‐ to co‐colonization is a critical factor in collective dynamics, by modulating the asymmetry in pairwise invasion fitness between types, and consequently, the dynamic complexity of the system as a whole. This ratio is key to observe the emergent context dependence of strain interactions (Bascompte, 2019; Coyte & Rakoff‐Nahoum, 2019) in our model. We argue that the analytically explicit form of this ratio in our formalism enables direct connection with the stress gradient hypothesis (SGH) in ecology (Bertness & Callaway, 1994; Callaway & Walker, 1997). This hypothesis postulates that as stress increases, the importance of positive facilitative effects increases in a community, whereas in benign environmental conditions, competitive effects are higher; a finding that emerges also from our results. In support of complex higher‐order dynamics emergent from simple pairwise interactions, with critical links between mean and variance, we uncover the exact formulation for why the sum as a collective is much more than its parts. Our results invite a deeper understanding of the biology of endemic multi‐strain systems and point to key global modulators of collective polymorphic coexistence in nature.

METHODS

N‐strain SIS model with co‐colonization

We study the epidemiological model (Madec & Gjini, 2020) given by the following ordinary differential equations:where gives the force of infection of strain , and the variables , and refer to the proportion of susceptible hosts, singly colonized by strain and co‐colonized by strains and . Notice that , thus the dimension of the system is effectively . A key assumption in the model is that hosts co‐colonized with different strains and transmit either with equal probability. In the above notation, , is the infected host turnover rate, encapsulating both clearance rate of colonization episodes and recruitment rate of susceptible hosts (balanced by natural mortality rate ).

This model has also been described in detail previously (Gjini & Madec, 2017) for the case of . The model allows for each strain to interact differently with other strains upon co‐colonization, altering the susceptibility of an already‐colonized host to incoming strains. The magnitude and type of such interactions are described by the matrix , where values above 1 indicate facilitation between strains, and values of below 1 indicate inhibition or competition between strains. We do not make any specific assumptions on the mechanisms underlying such interactions, but a key feature of this formulation is the explicit incorporation of intra‐strain and inter‐strain interactions. The derivation of the slow‐fast dynamics decomposition (Madec & Gjini, 2020) lies on the assumption that the variance of such interaction coefficients is low, thus deviation from neutrality of , with respect to a reference , is small. In our particular simulations here, we assume a normal distribution for such interaction coefficients. We can write every as: under the assumption that is small.

Although many different parametrizations are possible, we define as the statistical mean of the co‐colonization interaction coefficients between all possible strains in the available pool. We define the deviation from neutrality as the standard deviation of (K)1 ≤  ,  ≤  . This leads to the normalized interaction matrixto have the same distribution as , but with mean 0 and variance 1, at the start of dynamics in the unpruned community. We adopt analysis and simulations to understand the behavior of the epidemiological model for different assumptions on the co‐colonization interaction matrix, and for different system size . We use the reduced system (of equations) for dynamics on the slow timescale (Madec & Gjini, 2020) to obtain the steady‐states of the system and analyze their properties and stability (see Box 1).

RESULTS

In the present model, the entire system is structured as a collection of hosts that can be in different colonization states: susceptible, singly colonized, or co‐colonized. We consider a multi‐type infection, transmitted via direct contact, following susceptible‐infected‐susceptible (SIS) epidemiological dynamics with co‐colonization (Madec & Gjini, 2020). Adopting a general formulation for diversity, we assume there are N types, without specifying the mechanisms for their definition. Thus, with an ordinary differential equations model (see Section 2), we describe the proportion of hosts in several compartments: susceptibles, S, hosts colonized by one type I, and co‐colonized hosts I, with two types of each combination, independent of the order of their acquisition. The model structure follows that of Gjini et al. (2016); van Baalen and Sabelis (1995) allowing also for same strain coinfection (I). Hosts in the mixed coinfection compartment (I) transmit either strain with equal probability. Fitness differences in this system are encoded in how strains interact with each other upon co‐colonization (K), whether there is facilitation between resident and co‐colonizer (K > 1) or inhibition (K < 1), and its exact magnitude, – assuming equivalence in transmission β and clearance rate γ. The structure of the K, subject to small deviation from perfect symmetry, is central to multi‐strain dynamics. We take no account of other biological details such as mutation, seasonality in transmission, or heterogeneity in the host population, which may influence the transmission dynamics of particular strains. These exclusions serve our purpose to assess the impact of selection imposed by co‐colonization interactions between given strains in the host population on temporal trends in individual strain frequencies.

In an earlier mathematical investigation (Madec & Gjini, 2020), we have derived in detail two timescales in this multi‐strain system: a fast one, given by the neutral model, and a slow one governed by the variation in co‐colonization coefficients K. During fast dynamics, the system stabilizes conserved quantities such as the endemic prevalence of single and co‐colonization whereas over the slow timescale, strain selective dynamics unfold (see Box 1).

Deviation from symmetry in interaction trait space: the basis of the framework

We conceptualize each altered susceptibility to co‐colonization (i.e., interaction coefficient) between closely related strains, as a mean value k plus some deviation from symmetry, thus re‐writing it as:where is small. Thus, is the normalized interaction matrix, relative to the reference . The parameter in this case encodes mean interaction in co‐colonization between any two strains, which, if describes mutual facilitation, and if describes mutual inhibition. Another key epidemiological parameter is the basic reproduction number, , in our model given by β/m (Section 2), describing the intensity of transmission, and defined as the average number of secondary infections that arise from a single case in a naive population (Diekmann et al., 1990).

Fast and slow dynamics: derivation of the replicator equation in a new context

On the fast time‐scale (), strains behave as neutral (). The system reaches the central manifold where total prevalence of susceptibles , total prevalence of single colonization , and total prevalence of co‐colonization are conserved, and depend only on mean epidemiological parameters:but, where each individual strain frequency is neutrally stable and free to vary . On the slow time scale , consistent with weak selection, exact asymmetries between strains play out, and under the conservation law above, individual frequencies obey deterministic slow dynamics, given by an ‐dimensional replicator equation:where , is a quadratic term symmetric on all strains, encapsulating the effect of the system as a whole on each individual strain. The quantity is conserved and the rate is given by: . This replicator equation links directly epidemiological SIS dynamics to Lotka–Volterra systems (Bomze, 1995) and multi‐strain co‐colonization processes to ecology and evolution (Nowak & Sigmund, 2004).

From pairs to collective dynamics

The above equations drastically reduce the system from to dimensions. Furthermore dynamics are a direct function of pairwise invasion fitnesses between strains, which, for each pair , have the form:

The quantities denote the initial rate of growth of strain in an exclusion equilibrium where only strain is resident, and depend on the ratio between single and co‐colonization which is given by

Invariant principles in nonequilibrium multi‐strain dynamics

On the slow timescale, at all times, and for all strains, the following relationships hold:stating that , as a measure of dominance of strain in the total prevalence of carriage, occupies an equal relative frequency in single () and co‐colonization (). After solving for strain frequencies , the Equation (4) can inform the epidemiological variables of the original system (1) as follows:

How mean‐field parameters determine type and tempo of selection

What type of coexistence?

According to this model, global mean‐field parameters, such as those affecting transmission rate, , basic reproduction number , and mean interaction coefficient , explicitly impact the multi‐strain dynamics over long time scales (). The equation for strain frequencies contains information for how they set the speed () and mode of strain stabilization (). More specifically, the relative dominance of single colonization in the system appears as a factor in pairwise invasion fitness () of when invading a only equilibrium. The parameter modulates the importance of cross‐strain asymmetry in co‐colonization trait comparison between two strains. The higher is, the higher the prevalence of single colonization, thus the more important the asymmetry between how and invest in the mixed co‐colonization compartment . Because hosts transmit both and with equal probability, the relative superiority of will be reduced from such “altruistic” investment; thus, it appears with a negative sign in invasion fitness. While this component of invasion fitness of is sensitive to feedbacks from global transmission and mean parameters between strains, the other fitness component of in depends only on characteristics of the resident, namely on the relative strength of inter‐ versus intra‐ strain interaction of the resident (), which indicates how much the resident promotes its own coinfection, compared to its vulnerability to coinfection by .

Notice, that for fixed normalized variation among strains, hence rescaled matrix , increasing in the system, amplifies asymmetries among all pairs of strains. Mathematically this is reflected in the correlation between and tending to ‐1, which leads to more pairwise competitive exclusion. We illustrate this for the case of in Figure 1a, where for randomly generated rescaled interaction matrices (Normal distribution with mean 0 and variance 1), increasing stretches the region occupied by towards the competitive exclusion zone (, and and viceversa). This phenomenon applies also to the case of higher , where effectively the multi‐strain network is derived from pairwise invasion fitnesses between any two constituent strains (Figure 1b). As the ratio of single to co‐colonization, , increases, the proportion of edges between constituent strains, that lead to competitive exclusion increases, making collective coexistence more probable as a result.

FIGURE 1

Collective coexistence from pairwise invasion fitnesses and the effect of single‐to co‐colonization ratio . (a) Effects of on the repartition of the pairwise invasion fitnesses. In our model, there are four possible outcomes between any two strains (depending on the values of and ): extinction of , extinction of , coexistence or bistability of the exclusion equilibria, the latter being known also as priority effects (Ke & Letten, 2018). The coefficients of the matrix are randomly generated from a normal distribution . When , the partitioning (probabilities) between the four outcomes is the same. When increases, the strain pairs, for the same matrix, follow more likely competitive exclusion dynamics. Among 10,000 random simulations of A, as we increase , we find 51% of pairwise competitive exclusion for , 64% exclusion for and 92% exclusion for . (b) Effect of on the multi‐strain invasion (pairwise ) networks. Red line: coexistence, blue line: bistability, and a gray arrow indicates competitive exclusion (the arrow points to the winner between the two). Here we fixed the normalized interaction coefficients between strains (matrix ) and strains, and we varied . As increases even though the actual normalized interaction coefficients remain fixed, the effective outcomes between each pair of strains, being ‐dependent, tend to exclusion for , with the gray edges becoming more common

To illustrate concretely how the ratio impacts ecological coexistence, we visualize multi‐strain dynamics as a function of in Figure 2, for the case of and . Under fixed interaction asymmetries , changing the ratio between single and co‐colonization, shifts multi‐strain dynamics from stable coexistence towards limit cycles and ultimately towards heteroclinic cycle type oscillations, with effectively only one strain persisting for long periods of time, to be subsequently replaced by another one, and so on. We find that for and random , we have a case of a general replicator equation whose steady‐state analysis when may be reduced to Generalized Lotka–Volterra (GLV) dynamics with constant growth rates and random interactions. Whereas for the case of , dynamics converge to hyper‐tournament dynamics studied by Allesina and Levine (2011) for and Grilli, Barabás, et al. (2017) for the case in which . These extreme regimes of behavior are expected in the replicator equation for special cases, but the novelty in our framework is that we have identified a global system quantity, the ratio between single and co‐colonization, as a tuning parameter for moving our multi‐strain dynamics between such extremes (see Box 2).

FIGURE 2

Increasing the ratio of single to co‐colonization () increases the complexity of multi‐strain dynamics. We illustrate two examples with fixed interaction matrix and shifting for two system sizes: (a) and (b) . Here, we fixed the normalized interaction coefficients between strains (matrix ) and we varied from 0 to 0.5 and 1.5. Because affects the (Equation 5), for each , we obtain a different effective competition network between strains modulated by top‐down factors. The plots show the strain frequency dynamics following the replicator equation on the slow manifold. For and , for small value of closer to zero, the final attractor is a stable steady state, while for a larger value of , the stable attractor is a union of heteroclines of three strains (May Leonard type or rock‐paper‐scissors dynamics). For intermediate values of () in this case we obtain a stable spiral of coexistence for and a limit cycle of coexistence for , both consistent with an increase in rate of strain turnover

One limit for may be achieved by either increasing global facilitation , or increasing transmission intensity . The other limit can be achieved by increasing competition in co‐colonization or decreasing transmission intensity. Although in reality, these limits may never be literally attained biologically, here we describe the mathematical trends of the system, which will most likely reside in the intermediate range. For a more detailed analysis of system behavior and speed of the dynamics in these limits, see also Text S2, and Figures S2 and S3. A specific example with dynamics close to these limits is provided in Figure S4.

Limit 1: Co‐colonization, ( or )

Quality of the dynamics

Recall that implies that the strain (the invader) may invade the system with only strain (the resident) present, starting from a very small frequency. In the limit , this pairwise invasion fitness reduces to:

Importantly, in this limit, the fitness depends only on the coefficients of the resident, where measures how much the resident contributes to mixed co‐colonization with strains and : (transmission and fitness of other strains), and how much the resident contributes to self co‐colonization: (transmission and fitness of itself). Hence, the relative fitness of the invader depends only on the resident . Invasion is possible if and only if .

This phenomenon in a system of players can increase the availability of independent niches. For example, if for all and then all strains may be stable residents when alone, leading to at least stable monomorphic equilibria. Cycles or more complex dynamics are rare.

The final outcomes can vary because the structure of the matrix is free in principle. In the special case for all , then , and the dynamics could be richer, in line with generic replicator systems (Yoshino et al., 2008). In that special case, under random i.i.d. symmetric about 0, the probability of having a feasible equilibrium is exactly , where is the number of coexisting strains (Morrison, 2013; Serván et al., 2018). In general, for random , with no special structure, qualitatively, the system likely has many stable steady states wherein only few species coexist.

Speed of the dynamics

. If , we have . For fixed population turnover (m: clearance and death rate), the selective dynamics happen very slowly, that is, the dynamics are effectively neutral. Similarly if via lower (), the multi‐strain system stays close to neutrality.

. If , via increases in transmission rate, , then which is far from zero. In this case, the model reduction captures well the real dynamics and the qualitative study of is appropriate.

Limit 2: Single colonization, ( or )

Recall the matrix is the matrix of all pairwise invasion fitnesses between strains in the system. In this limit, we have generically ||Λ|| → +∞ and . In order to keep a bounded matrix , we rewrite the system as

Then, the speed is given by and the qualitative behavior by . Notice that when , the matrix

Importantly, this matrix is skew symmetric, for which there are known mathematical results of the replicator equation in zero‐sum games (see Allesina & Levine, 2011; Chawanya & Tokita, 2002; Fisher & Reeves, 1995). The quadratic term in Equation (9) is zero: . The ‐strain dynamics reduce to:

There exists exactly one nonnegative linearly stable equilibrium (in particular multistability is impossible). In practice, similar to the classical Lotka–Volterra model, there is a one‐parameter family of limit cycles parametrized by initial conditions. This type of coexistence is not structurally stable and will be lost for a large but finite value of leading in general to heteroclinic limit cycles among strains.

Denoting the number of strains coexisting, for this skew‐symmetric case, it has been shown that only an odd number of strains may coexist (see also Figures S1 and S5). The probability to observe strains, out of a total pool of is

In this case, at equilibrium , hence is an eigenvector associated with a zero eigenvalue. For a general skew‐symmetric matrix, a zero eigenvalue is expected only when its size is odd (see also Grilli, Barabás, et al., 2017).

If decreases transmission intensity, either via lower or higher , then , and selective dynamics are too slow, making the system effectively neutral.

. If , . Thus dynamics are well defined and occur on a feasible timescale.

How many strains coexist?

Using simulations with random strain interactions and varying , we find that the mean number of strains that can coexist () increases with (see full distributions in Figure S1 for ). In the limit and with the additonal assumption that for each , the probability of a feasible (positive) steady‐state of strains in our replicator equation is exactly the same as that of a positive steady state in GLV dynamics with equal growth rates and random interactions: . As increases, the number of coexisting strains tends to , with where denotes the total pool size. In particular, when is low, lower numbers of coexisting strains are progressively more probable, but for large values of the ratio , favoring single to co‐colonization, stable coexistence tends to become restricted to only an odd number of strains, a feature expected in special cases of the replicator equation (Chawanya & Tokita, 2002). Such case of perfect pairwise exclusion between species has emerged as the best‐studied stabilizing competitive network involving intransitive competition among species, known as rock‐paper‐scissors tournament games (Allesina & Levine, 2011; Kerr et al., 2002). In our system, this extreme, obtained in the limit , emerges as a special case of a far more general gradient mediated by different values of , where edges between any pairs of strains can vary among coexistence, exclusion and bistability, and with quantitative subtleties like in hypertournament games.

How fast is the dynamics?

Thanks to the explicit formula for this critical ratio, one can immediately see that can change in two ways, either by changing basic reproduction number , or by changing , the mean interaction coefficient between strains in co‐colonization. If overall colonization increases (), or if facilitation between strain increases (), decreases, and viceversa: when total prevalence is lower, or strains compete more in co‐colonization, increases. So far, we have seen that the ratio determines completely the type of multi‐strain equilibrium, for a given standardized variation in strain interactions. But how fast is this equilibrium reached? This is determined by , also explicit, in our replicator equation describing multi‐strain dynamics over the slow time . As changes, not only the selection “movie” changes, but also the effective speed at which such “movie” is played. Essentially, decreases with , hence selection is slower with increases in or . But there are in general three cases, illustrated in Figure S2. When changes in are obtained via , if is small (i.e., large) it is possible that increasing increases . This would be obtained either from lower or higher in . In such case, reducing colonization prevalence from relatively high to lower levels (via lower strain transmission/growth, or faster host birth‐death, shorter colonization episodes) not only acts to increase the complexity of the dynamics, but also to speed up multi‐strain selection in the biological co‐colonization “game”.

Coexistence, stability, and diversity

The community dynamics resulting purely from co‐colonization interactions between strains can be very complex, and ranges from simple coexistence equilibria to limit cycles and even possible wildly oscillatory dynamics. Noting that dynamics among strains are bounded in our system, failure to identify existence of a steady state necessarily implies the existence of another long‐time attractor, of a more complex nature. We find the dimension of the attractor is 0 for small and becomes 1, and may even exceed 1, for larger values of this parameter. This is what we refer to as complex coexistence.

As the important debate in ecology about the relationship between stability and diversity in ecosystems is still ongoing (Landi et al., 2018; May, 1972; McCann, 2000; Odum & Barrett, 1971; Tilman & Downing, 1994), here we set out to examine this relationship in our system, assuming random interaction structure between strains, described in the matrix , and fixing the total pool size of strains . Thus, we study the model's rich behavior for a relatively larger pool of strains in the system, namely (Figure 3). We sample the normalized interaction matrix randomly from a normal distribution . For each interaction matrix, we compute numerically the feasible equilibria of the system (verifying ), and for each equilibrium found, we evaluate the number of strains coexisting, , the local stability of the steady state (see Appendix), and its associated Shannon entropy also known as species evenness in ecology. By simulating different randomly generated normalized interaction matrices , we can explore regimes of exclusion, multi‐stability, multi‐strain coexistence, limit cycles and even chaos. In this analysis we use the definition of equilibrium stability (McCann, 2000), although there are also other notions of general stability related to permanence (Law & Morton, 1996) that could be informative.

FIGURE 3

Diversity‐stability relationships for ‐strain coexistence starting with a pool of strains. (a) Summary of 100,000 simulations with a random normalized interaction matrix , with each entry drawn from the Normal distribution . . The equilibria of each system are mapped as dots in this plot, after computing their local stability and Shannon entropy. (b) Effective stability is calculated as a product between the mean dominant eigenvalue and proportion of steady states at each entropy level, stratified by number of strains coexisting . (c) Plot of the rank‐ order frequency distribution of strains at the “optimal intermediate evenness” for strains coexisting out of a pool of 10. (d) Plot of the rank‐ order frequency distribution of strains at the “optimal intermediate evenness” for strains coexisting out of a pool of 10

Analyzing coexistence patterns under random allocation of co‐colonization interaction coefficients (Figure 3), there are multiple phenomena to note. First (and trivially) the higher the number of strains coexisting (), out of a pool of strains, the higher the maximal Shannon entropy (evenness) in coexistence, delimited mathematically by , a function that increases with (dashed lines in Figure 3a). Second, within a given , the higher the evenness in coexistence, the higher the stability of such steady‐state, as measured by the negative part of the dominant eigenvalue of the Jacobian matrix evaluated at that equilibrium, a result similar to Song and Saavedra (2018) using a Lotka–Volterra framework. In contrast, across different , as maximal entropy increases, the maximal stability of coexistence steady states tends to decrease (color peaks decreasing diagonally down from left to right in Figure 3a). This implies that under randomly sampled interactions, hence in randomly assembled communities, stable coexistence equilibria with more strains, and more evenly distributed strain frequencies, become harder to obtain.

Next, when formally computing the probabilities of reaching a stable ‐strain coexistence in this model, we find that the relative probability of stable ‐strain steady‐state decreases with , thus confirming that the more strains there are, randomly interacting, the harder it is for all of them to coexist (relative comparison among colored dots frequencies in Figure 3a). In contrast to the overall stable steady states (2%), there are many more unstable steady states in the system (98%), spread at any entropy (evenness) and stability level, with no apparent order or hierarchy (black dots in Figure 3a). This indicates that stable coexistence at a steady state is rare under random co‐colonization interactions, and that the attractors are often of dimension 1 or more. This also suggests that the number of coexisting strains in that case is not a limiting factor.

Quantification and stratification of these probabilities jointly by , stability and evenness, combines the feasibility of a given state with its expected stability and strain composition. Then, to visualize in more detail this distribution, among steady states with a given , for each entropy level , we multiply the existence probability of a stable steady state (feasibility) with the absolute value of mean dominant eigenvalue (mean stability) of such steady state and plot the resulting effective stability for each (Figure 3b). This integrates two opposing forces: on one hand higher‐evenness equilibria are less feasible and harder to reach under random sampling of interaction coefficients from a pool of strains, that is, decreases with ; on the other hand, more evenness among strains is associated with more stability in coexistence, increases with , similar to (McCann, 2000). As can be seen in the emergent peaked curves in Figure 3b, such trade‐off between probability of stable coexistence and stability gives rise to an optimal intermediate evenness for ‐strain coexistence, for any , where stable coexistence is sufficiently feasible in the first place, and secondly where the stability of that equilibrium to perturbations is sufficiently high. Thus, mean stability of ‐strain coexistence is effectively maximized at particular, typically intermediate, values of Shannon entropy (Figure 3c,d). In practice, independently of particular strain composition, out of a total pool of , for a given number coexisting, , the rank‐frequency distribution of strains should be preserved. Such “optimal” rank‐order frequency distribution will display higher variance for low number of coexisting strains (for example ), but be more robust for higher (e.g., ).

In Figure S5, we explore the same phenomena for the case of a higher , tending to dominant single colonization ( more likely applies to a system like pneumococcus serotypes, Gjini et al., 2016). We find the ecological forces leading to an optimal effective stability at an intermediate entropy range remain robust also in this case. This limit also verifies that for the same randomization of interaction coefficients , typically more strains coexist than for lower , and in particular, tending to an odd number of them (Chawanya & Tokita, 2002), while their joint dynamics become more complex (Box 2).

The role of the ratio in multi‐stability

With the richness of potential dynamics in mind for a fixed (Figure 3), we next revisit the relationship between complexity and ratio of single to co‐colonization, with model simulations where is varied. Using 100,000 random system simulations with for different values of , we confirm the pattern observed for and illustrated for in Figure 2, that increasing (via decreasing either or ) increases strain turnover and the complexity of multi‐strain coexistence. Computing the number and stability properties of equilibria for each system, we can see that multi‐stability of steady‐states is very common, especially in the limit of single colonization dominance () as shown in Figure 4a. By multi‐stability we broadly refer to multiple “alternative stable steady‐states” in the system, which may contain different strains, a number of which may also overlap. Pure multi‐stability, with 100% overlapping constituent strains, but different frequencies, between two steady states, is impossible. In contrast, as increases, multi‐stable equilibria become less common but the proportion of systems without any stable steady state increases, pointing to an increase in dynamic complexity for the same and normalized interaction matrix between strains.

FIGURE 4

System complexity as a function of and multi‐stability. (a) For several values of , 100,000 systems of strains have been generated (using a normal distribution of ). For each value of , we show the proportion of systems with a given number of stable steady states (in ). As increases, the proportion of systems without any stable steady state increases, reflecting larger potential for complex coexistence dynamics, where the attractors may be limit cycles, heteroclinic cycle or even strange attractors. For (high transmission or competition between strains), stable multi‐strain coexistence is very probable, and most of the systems have several stable steady states, thus multi‐stability is common. (b) In the particular limit case of , and for each of the 100,000 systems, we compare the strain overlapping index SOI (see Appendix) within all steady states and within only the stable steady states. We see that without the constraint of stability, two steady states of our system share on average 20% of their coexisting strains, while within the set of stable steady states, this overlapping index is nearly always . This shows that in multi‐stability, the attractors contain mostly disjoint sets of strains

In particular, when zooming‐in on the alternative stable steady‐states of each system, we can quantify how different are the subsets of strains coexisting in each of those equilibria. Defining an index of strain overlap, based on the Jaccard index between two sets (see Appendix), we find that multi‐stability is largely characterized by nonoverlapping subsets of strains (Figure 4b), in contrast to the higher average strain overlap (about 20%) between any two generic steady states. This means that depending on initial conditions and founder effects, the same global pool of strains with random co‐colonization interactions, can lead to entirely different community composition of strains in different epidemiological settings, an effect that is rather common (60% probability) for values of around 0.1, but becomes improbable for .

When quantifying system outcomes for strains, in three broad categories: those leading to a unique stable steady state, “monostability”; those having multiple alternative stable steady states, “multistability,” and those leading to complex attractors but no stable steady state, as we apply an even higher resolution for , we find that there are effectively three regions as a function of single to co‐colonization prevalence (Figure 5). For small , stable coexistence and multi‐stability are more likely, for intermediate , a single stable steady state is more likely, and for high , it is very rare for the system to have any stable steady state, a regime we refer to as no stable steady state, implying complexity. In this regime, rock‐paper‐scissors type coexistence and unstable coexistence are the more, if not the only, probable outcomes. This result highlights the importance of an environmental mean‐field gradient in shifting the qualitative dynamics: moving from higher endemic prevalence or more facilitation in co‐colonization towards lower prevalence or more competition, the inter‐dependence between strains increases, the role of asymmetric investment in co‐colonization increases, and thus, the propensity for complex high‐dimensional coexistence becomes higher. In one extreme, we expect low variability in a system over time, while in the other, we expect high variability and complex strain turnover over time.

FIGURE 5

Dynamic regimes as a function of single to co‐colonization ratio . We plot the number of steady states as a function of for many randomly generated interaction matrices . Since depends on and , we kept their combination implicit, and randomly sampled the rescaled interactions matrix . For any system (), if there are two or more stable steady states, the system is considered multistable and the dynamics depend on the initial conditions. If there is only one stable steady state the dynamics can be considered simpler. Indeed every system with a unique global attractor is monostable, but the converse is not true for a monostable system may also have non stationary attractors. If there is no steady state, the dynamics are more complex, meaning the attractors are necessarily neither stationary like limit cycles, heteroclinic cycles or even strange attractors. For the simulations, we used a total strain pool size of , and for each value of we generated randomly matrices leading to different interaction systems for each value of . For a given value of , we then counted the proportion of systems with no stable steady state (red line), exactly one stable steady state (black line) or several stable steady states (gray line). For a small , very few systems have zero stable steady states, while nearly of them are multistable. For a large value of the multistability is nearly impossible while more than of interaction systems having no steady state and thus displaying a complex dynamics with no stationary attractors

Contrasting patterns of evolution along the gradient

Until now, we have seen how for a fixed and a given set of interaction coefficients between strains, their relative abundances change at the population level. Next, we consider how the epidemiological dynamics and the evolutionary dynamics interact. Evolution happens on several levels in this system, made explicit by dynamics and various possible trait definitions for each system member . Here we focus on two important levels. First, frequency dynamics () can be linked with selective dynamics on mean interaction in co‐colonization trait space (see Madec & Gjini, 2020) where k effective = k + εq(t). Secondly, another trait changing in the community is mean invasibility of the system, , in mutual invasion fitness space (see Box 1). These two quantities reflect the changing mean fitness landscape, and depend on frequency dynamics as quadratic terms involving summation over products of strain pair frequencies over time. When varying , as the multi‐strain selective dynamics unfold in different ways, so do the evolutionary dynamics of mean traits in the system (Box 3).

In multi‐strain communities with random co‐colonization interaction strengths (), the global epidemiological quantity , describing the single to co‐colonization ratio in the host population, modulates regimes of system behavior, shaping its stability, complexity and diversity properties.

Alongside transmission dynamics among hosts, we may track evolution in real time on at least two traits in this system:

mean interaction in co‐colonization trait space (q) where k effective = k + εq(t),

mean invasibility of the system (Q) in mutual invasion fitness space among strains (Equation 4).

Evolutionary dynamics in the microbial community reflect the changing mean fitness landscape, and depend on strain frequency dynamics (Madec & Gjini, 2020) as follows:

Notice that q increasing means the system tends to increase facilitation in co‐colonization on average over time, whereas q decreasing means the multi‐strain community tends to increase inhibition in co‐colonization over time. Similarly Q increasing implies the system becomes less invadable over time, and Q decreasing means the community of strains becomes easier to invade from outsider strains. As frequency dynamics can be complex, so can mean trait dynamics. An indirect measure of the diversity of the system is given by the magnitude of Q. A dynamic measure of the “selfishness” in the system is given by the difference:which indicates the relative strength in the system of self‐positive feedbacks in co‐colonization. If the system selects for strains which are better at favoring self than others, and viceversa. If , the system selects for strains which, on average, are better at favoring others than self (“altruists”). A special case are competitive exclusion scenarios (Madec & Gjini, 2020), where a decrease in () will be associated with a positive dynamics of , illustrating how extreme selfishness leads to selection of a single strain in the system.

Mean trait evolution on 2 levels: and over time

Resistance to invasion tends to be stable

Mean trait evolution on 2 levels: and over time

Resistance to invasion tends to oscillate

In particular, in the limit, where multi‐stability is common, systems can slowly evolve towards higher mean competition (), or higher facilitation () in co‐colonization depending on initial conditions, and on the strains that get selected in the pruned community to stably coexist. Conversely, in the limit, where complex oscillatory coexistence is more probable, systems are more likely to preserve many strains in oscillatory dynamics, and thus fluctuating polymorphism in interaction trait space, with alternating periods between higher average competition and higher average facilitation ( oscillatory). This result highlights two routes to coexistence or maintenance of biodiversity: little variability within one population but diversity between populations in the first case, and more variability within one population but low diversity between populations in the second. Notice also that an increase in mean resistance to invasion of the multi‐strain system ( dynamics) does not necessarily always correlate with increased facilitation or competition (higher or lower ), and in the oscillatory regime the two can also be synchronous and asynchronous (Figure S6). As can be taken as a relative measure of selfishness versus altruism in the system (Box 3), in the asynchronous oscillatory regimes for high , we can expect fluctuating prevalence of “altruist” and “selfish” strains, all maintained over time.

Linking to the stress gradient hypothesis

Our main result is that the ratio of single to co‐colonization is a key modulator of the qualitative complexity of the dynamics (Figure 5, Box 2). While both transmission intensity and mean interaction coefficient in cocolonization decrease this ratio, a natural corollary is that the only way for to be held constant, is if and are traded‐off against each other. In other words, when colonization opportunities for all strains are reduced, in order to keep constant, the mean interaction coefficient must increase, thus enhancing facilitation among strains in co‐colonization. And viceversa, if the environment gets more favorable, that is, increases, then can only be be kept constant if decreases accordingly, thus if all strains become less vulnerable to co‐colonization. Such quantitative feature of our model makes an explicit link with the stress‐gradient hypothesis (SGH) postulated in ecology (Bertness & Callaway, 1994; Callaway & Walker, 1997; Chamberlain et al., 2014). The SGH predicts that positive interactions should be more prevalent in stressful environments, while more favorable environments should favor competition. This hypothesis has been supported and tested mainly in plant systems (Callaway et al., 2002; Eränen & Kozlov, 2008; He et al., 2013; Pugnaire & Luque, 2001), and more recently in microbial communities (Fetzer et al., 2015; Hoek et al., 2016; Lawrence & Barraclough, 2016; McCluney et al., 2012; Piccardi et al., 2019).

In our epidemiological system, with microbial interactions embedded in co‐colonization dynamics, the overarching single‐to‐co‐colonization ratio can be considered as a conceptual and mathematical formalization of the stress‐gradient hypothesis. In order for a set of species with a fixed set of normalized interactions (A = ()) between them, to maintain a given configuration of coexistence () under changing environmental conditions (e.g., ), the only way is by universal adaptation in mean‐field interaction, increasing antagonism towards co‐colonizers in favorable environments and increasing facilitation in harsher environments, respectively (see Figure 6). However, even if is kept constant, the speed of collective dynamics may still be affected, typically reduced with , although the absolute magnitude will depend on whether and how they are mediated through transmission rate versus duration of carriage (see Figures S3 and S4). In summary, while the fitness landscape is shaped by the phenotypic distribution of all the community members, as the external environment changes, the relative balance of their interactions may shift, due to ultimate context‐dependence in mutual invasion fitness. One solution to keep that balance in check would be by quick reversal adaptation of the mean.

FIGURE 6

The versus trade‐off in and the stress gradient hypothesis. While the type of equilibrium the system will tend to, is strongly driven by (the ratio of single to co‐colonization), hence directly by and , it is clear from the expression that the only way for this ratio to be held constant, is if and are traded‐off against one another. This principle reflects a core expectation of the stress‐gradient hypothesis (SGH): when transmission opportunities for all strains are reduced (harsher environment), the mean interaction coefficient in the system must increase, that is, if all strains cooperate more in co‐colonization, to keep the same and coexistence regime. And viceversa, if the environment gets more favorable, that is, increases, then can be kept constant only if decreases, thus if all strains increase competition. Note that keeping fixed only preserves the qualitative multi‐strain dynamics, but does not guarantee its speed will remain unchanged. Different global and and the underlying route to such changes may effectively impact selection speed

Taken together, our results show that even though differences in pairwise co‐colonization interactions can be small and random, as expected for closely related strains (or species), they are enough to drive selection in a large community of such members, which can collectively coexist, adapt, and self‐organize, albeit often in complex and unstable fashion. The critical influence of average growth potential and mean “willingness to share” co‐colonization, analyzed here, highlights the importance of global context, and calls for detailed system‐specific investigations in the future.

DISCUSSION

A central question in microbial ecology is whether community members compete or cooperate with one another, and how such interactions mediate community stability, resilience and function. In our model, we study an epidemiological multi‐strain system, where members interact with each other via altered susceptibilities to co‐colonization, which broadly include both competition and facilitation. By clearly delineating the role of mean interaction coefficient, basic reproduction number , and biases in pairwise coefficients relative to the mean, we obtain a model reduction for strain frequency evolution for any . The questions we addressed within such a system are similar to a long‐standing ecological quest on multispecies dynamics (Bunin, 2017; May, 1972; Pascual et al., 2006; Serván et al., 2018; Song & Saavedra, 2018): what governs coexistence regimes, diversity and stability in such systems?

Although we do not specify the molecular mechanisms that can mediate positive or negative interactions between strains in co‐colonization (Dawid et al., 2007; Leggett et al., 2014; Lysenko et al., 2010; Riley & Gordon, 1999; Shen et al., 2019), this renders the system generic and broadly applicable. We can quantitatively predict very important features of system behavior with a simple mathematical framework derived under quasi‐neutrality and strain similarity assumptions (Madec & Gjini, 2020). This model reduction captures selective dynamics between strains over long time, and coincides with an instance of the replicator equation from evolutionary game theory (Hofbauer & Sigmund, 2003; Nowak & Sigmund, 2004). Studying this equation and its biological consequences in detail, here we find that in our system, there is a critical ratio that tunes complexity and dynamic regimes in such multi‐strain contagion context: namely the ratio of single to co‐colonization. This ratio, , is given by the inverse of the product of the basic reproduction number and mean vulnerability to co‐colonization . We show that it amplifies the importance of asymmetry in interaction between strains, lending theoretical support to the principle of context‐dependence (Coyte & Rakoff‐Nahoum, 2019) of relative fitnesses in a coupled microbial community.

As strain relative abundances fluctuate, the resource and fitness landscape change dramatically and feed back on the system. This dynamic interdependence is on one hand, the source of complexity, but also a key driver of evolving mean fitness between extant strains. A global, symmetric and temporally varying, environmental feedback on all strains emerges naturally in co‐colonization dynamics ( in Equation 4), describing mean invasibility and following the evolution of diversity. Mean invasion fitness will tend to increase when more strains are stably maintained over time, and it will tend to zero in the extreme cases of competitive exclusion. Depending on parameter values, there is a gradient for number of resources in the system. Strains compete more strongly for susceptible hosts (one limiting resource) in the extreme of co‐colonization (), whereas they compete for (at most) singly colonized hosts, open to co‐colonization ( resources) in the opposite extreme of single colonization dominance ().

We find that for small values of the ratio between single and co‐colonization , which means high transmission intensity or high cooperation among strains on average, the system dynamics is characterized by multi‐stability, where typically nonoverlapping strain subsets can coexist, depending on initial conditions, a result that resonates with earlier multi‐strain SIRS models with cross‐immunity (Gupta & Anderson, 1999). In our model, we find a similar principle applies, despite the interactions between strains being mediated via altered susceptibilities to co‐colonization, without persistent immune memory. While in the present setup this gradient emerges naturally from the intrinsic structure of infection‐mediated interactions among all strains, similar gradients have been found also in ecological communities studied with generalized Lotka–Volterra models, when the correlation between and interactions was explicitly varied (Bunin, 2017).

Our study brings together several themes of interest across many multi‐type systems shaped by interactions and higher‐order feedbacks between their members. Many fields including ecology, geophysics, and economics are calling attention on critical transitions, which occur when natural systems drastically shift from one state to another (Scheffer et al., 2012). Critical transitions in the epidemiology of infectious diseases are of relevance to the emergence of new pathogens and escape from control, such as vaccines. The critical transitions analyzed in this paper relate global and mean‐field environmental variables to the manifestation of competitive hierarchies between multiple strains interacting in co‐colonization. We have made explicit how a gradient emerges from the epidemiological ratio of single to co‐colonization, and how it tunes effectively the diversity, stability and complexity of the coexistence between strains. Such gradient can mediate critical transitions in collective dynamics, when the normalized interaction coefficients between members are held fixed. These transitions may underlie and potentially enhance (or counteract) efforts to control and eliminate multi‐type infectious pathogens, as via vaccines or drugs, or in the face of climate change. Other studies have shown that concurrent multiple infection in malaria creates tipping points that give rise to hysteresis in responses to control or seasonal variation in vector abundance (Alonso et al., 2019). Our work supports a similar perspective, but more generally relevant to interacting systems with multiple strains, and coexistence regimes, rather than prevalence tipping points. Mean facilitation and competition among strains, affecting , appear as two sides along a continuum for the system, which particular strain compositions or environmental drivers (seasonality, general host immunity, population turnover) may tip towards one or the other extreme.

We find that when tends to favor co‐colonization, for example in the limit of more facilitation between strains on average, the system tends to multi‐stability and stable coexistence of a few strains in simple dynamics. In contrast, when tends to favor single colonization, for example in the limit of more competition between strains on average, the system tends to more complexity and unstable equilibria (Figure 5), but coexistence of more strains becomes possible (Figure S1). Although at first sight this may seem to suggest, somewhat contrary to previous expectation from lower‐dimensional models (Chen et al., 2017; Hébert‐Dufresne & Althouse, 2015), that average cooperation in co‐colonization is stabilizing and average competition is destabilizing, our result should be related to the fact that in our model is not a measure of cooperation in the classical sense, whereby cooperation is exclusively defined as a between‐strain phenomenon. In our model includes self and nonself‐interaction, and is critical to the global resource dynamics shaped dynamically by strains.

Our study has also several limitations. Although we explored the qualitative aspects of the dynamics, including the complexity and stability of steady states, their dimensionality and associated entropy, we did not study which strains ultimately coexist. This question relates to optimal strategies among players in a game theoretic context. How the interaction traits of each member determine its persistence or exclusion from the system, and what is the role of , requires further investigation and algorithmic optimization (Bárány et al., 2007). We also did not develop all the links with the Lotka–Volterra modeling literature in microbial ecology, although special cases of our model are similar to particular cases of GLV dynamics. The metaphor of co‐colonization, as adopted here, could be applied to more general cases of ecological dynamics between microbial species or ecotypes. We did not study alternative distributions of rescaled interaction strengths , only focusing on a symmetric distribution around 0. It would be interesting in the future to test our findings against nonrandom topologies, and empirical interaction networks, as studied for example by Grilli, Adorisio, et al. (2017). In the limit the patterns analyzed here still hold independently of the distribution governing , but in the limit , the distribution of matters, and different distributions may lead to slightly different numerical predictions.

Finally, this model makes several predictions which can be tested empirically. First, the invariant principles in the slow time scale (Box 1) suggest that dominance patterns in single and co‐colonization of particular strains should be the same. The other finding that stable (multi‐stable) coexistence between types through co‐colonization is more likely when mean interactions tend towards cooperation, and that a single stable coexistence between types becomes more probable at intermediate values of , and ultimately only unstable coexistence is possible for large values of could be tested in endemic multi‐type microbial ecosystems. For example, empirical data in polymorphic Streptococcus pneumoniae bacteria, have been consistent with estimates of about 90% mutual inhibition between co‐colonizing serotypes (Gjini et al., 2016; Lipsitch et al., 2012), and values around 2. This implies that for this system, , and multi‐stability is highly unlikely but a single stable equilibrium point is almost as likely as complex unstable coexistence. This may reconcile the secular trends observed in some settings consistently over many years (Ekdahl et al., 1998; Feikin & Klugman, 2002; Fenoll et al., 1998). Such secular trends can interfere with vaccine introductions and need to be accounted for when estimating impact (Moore, 2009; Moore & Whitney, 2008). Our model makes explicit predictions about the timescale and qualitative aspects of such secular trends (Box 2, Figure S4), under the plausible assumption that they are driven by co‐colonization interactions.

The requirement of being simultaneously stable and feasible tends to push coexistence regimes toward intermediate entropy, independently of precise strain composition. The consistency of the optimal evenness (Shannon entropy) configuration for a given number of strains has been observed empirically for pneumococcus serotypes across geographical settings, before and after vaccination (Hanage et al., 2010). Our model, offers a new perspective on such observations, thanks to co‐colonization interactions between serotypes. As this optimal rank‐order abundance distribution depends on the context , model expectations for such dependence could be tested with multi‐site data from different endemicity levels.

We postulate that the stress gradient hypothesis (SGH) (Bertness & Callaway, 1994; Callaway & Walker, 1997; Chamberlain et al., 2014) could help interpret the critical role of multiple infection in shaping epidemiological multi‐strain systems. Our link suggests that to preserve a certain “optimal” single‐to co‐colonization ratio (optimal complexity/coexistence balance), independently of community size , facilitation in co‐colonization between microbial strains should be more common in settings with low prevalence, a prediction to be tested in the future. It becomes intriguing to verify, beyond pneumococcus, to what extent this model and its insights (Box 3) can be used as an analytic backbone to interpret multi‐strain dynamics in other systems of relevance for public health, for example influenza (Yang et al., 2019), dengue (Mier‐y Teran‐Romero et al., 2013), malaria (Alonso et al., 2019; Gupta & Maiden, 2001) or human papilloma viruses (Murall et al., 2014).

Disentangling multi‐strain interactions and their role in community function at the epidemiological level remains challenging, but can be made more accessible analytically using frameworks such as the one proposed here. With the simplicity and deep insights afforded by this model, we can address better the role of mean fitness of the microbial system as a whole, trait variance, and the role of environmental gradients for stabilizing versus equalizing forces in biodiversity.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Erida Gjini: Conceptualization (lead); Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (lead); Software (supporting); Visualization (supporting); Writing‐original draft (lead); Writing‐review & editing (lead). Sten Madec: Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (supporting); Software (lead); Visualization (lead); Writing‐original draft (supporting); Writing‐review & editing (supporting).

Supplementary Material

Click here for additional data file.

High prevalence or high facilitation (μ low)	Low prevalence or high competition (μ high)
Co‐colonization dominates Few strains coexist Multistability is common Simple coexistence (fixed point) Low variability over time in one population Mean trait evolution on 2 levels: q and Q over time Resistance to invasion Q tends to be stable Either competitors or facilitators can be selected Evolution of “specialist” communities High diversity between host populations	Single colonization dominates Many strains coexist No multi‐stable steady states Complex attractors (limit or heteroclinic cycles) High variability over time in one population Mean trait evolution on 2 levels: q and Q over time Resistance to invasion Q tends to oscillate Competition and facilitation fluctuate Evolution of “generalist” communities Lower diversity between host populations