Towards a mathematical understanding of invasion resistance in multispecies communities — TEI-tagged view

Erida GjinipersName, Sten MadecpersName (2023date). “Towards a mathematical understanding of invasion resistance in multispecies communities.” Royal Society Open Sciencetitle·j.

Object of study — invasion resistanceterm·dom of multispeciesterm·key microbial communities modelled by a replicator equationterm·mth

Research question — how colonization resistanceterm·key can be described mathematically through mean invasion fitnessterm·key in a replicator system of N species

Contribution — a replicator equationterm·mth framework in which invasion resistanceterm·dom of a system appears explicitly and can be studied

Full text — tagged

Towards a mathematical understanding of invasion resistanceterm·dom in multispeciesterm·key communities Abstract Multispeciesterm·key community composition and dynamics are key to health and disease across biological systems, a prominent example being microbial ecosystems. Explaining the forces that govern diversity and resilience in the microbial consortia making up our body’s defences remains a challenge. In this, theoretical models are crucial, to bridge the gap between species dynamics and underlying mechanisms and to develop analytic insight. Here we propose a replicator equationterm·mth framework to model multispeciesterm·key dynamics where an explicit notion of invasion resistanceterm·dom of a system emerges and can be studied explicitly. For illustration, we derive the conceptual link between such replicator equationterm·mth and N microbial species’ growth and interaction traits, stemming from micro-scale environmental modification. Within this replicator framework, mean invasion fitnessterm·key arises, evolves dynamically, and may undergo critical predictable shifts with global environmental changes. This mathematical approach clarifies the key role of this resident system trait for invader success, and highlights interaction principles among N species that optimize their collective resistance to invasion. We propose this model based on the replicator equationterm·mth as a powerful new avenue to study, test and validate mechanisms of invasion resistanceterm·dom and colonization in multispeciesterm·key microbial ecosystems and beyond.

Keywords: colonization resistanceterm·key, microbial ecology, multispeciesterm·key community, pairwise invasion fitnessterm·key matrix, replicator equationterm·mth, system invasibility Introduction Deciphering the rules and mechanisms governing dynamics of living systems is a daunting task. Theoretical models contribute to meet this challenge, as they provide formalisms for precise quantitative descriptions of phenomena, clear hypotheses formulation and cause-and-effect predictions within a set of pre-stated assumptions. One of the very widely used models in ecology today is the Lotka–Volterra modelling framework for N species, which abstractly describes how the abundances of a set of species change with time as a function of their growth rates and mass-action interaction parameters, typically describing competition [ 1 – 3 ]. In this model, the reproductive rate (fitness) of each species depends linearly on the abundance of other species. This model has a long history of study, refinement and applications in many areas of biology and ecology [ 3 – 7 ].

Another classical model originating from evolutionary game theoryterm·key is the replicator equationterm·mth formalism [ 8 ], tracking the relative evolution of N competing strategies, and how they spread in a population, and is mostly used in population genetics, economics, theoretical physics and evolutionary dynamics. Although it is known that the replicator equationterm·mth is closely related to Lotka–Volterra model (topologically equivalent, to be mathematically strict) [ 9 , 10 ], most models of multispeciesterm·key dynamics utilize the Lotka–Volterra frameworkterm·mth, while applications of the replicator equationterm·mth in ecology remain rarer.

Here, we propose a replicator equationterm·mth framework to model multispeciesterm·key dynamics where species frequencies can be tracked explicitly from the N variables of the system. The key aspect of evolutionary game theoryterm·key, explicit in the replicator equationterm·mth, is frequency-dependent selectionterm·dom, whereby the fitness of an individual depends on the frequency of other ‘strategies’ in the population. We take advantage of this explicit nature to highlight that furthermore, within the replicator equationterm·mth framework, an explicit notion of invasion resistanceterm·dom of the system emerges naturally as a term in the equation and can be studied. In contrast, it is not natural to see mean invasion fitnessterm·key of a system explicit in a typical Lotka–Volterra frameworkterm·mth.

Invasion resistanceterm·dom is important in all ecosystems [ 3 , 4 , 11 , 12 ], and the replicator framework is flexible to be applied to many modelling contexts, but in the following, without loss of generality, we will first elaborate an illustration focused on a polymorphic microbial ecosystem, and then synthesize some broader qualitative insights for multi-species invasion resistanceterm·dom under replicator dynamics and various special cases.

From empirical colonization resistanceterm·key to mathematical invasion resistanceterm·dom . A pertinent example of a microbial ecosystem where invasion, interaction and colonization processes among several species are important is the commensal microbiotaterm·key, which plays a crucial role in health and disease. Recent studies are increasingly uncovering one of its clearest contributions, across organisms: protection against invading pathogens, known as colonization resistanceterm·key [ 13 , 14 ]. In humans, this protection is important when enteric bacterial pathogens challenge the gastrointestinal tract, but applies also in the nasopharynx, skin, and other environments which are populated by rich consortia of commensal microbes. The composition of gut microbiotaterm·key affects several collective metabolic functions and feedbacks with the host’s immune system, leading to or sometimes impairing immune homeostasis. Disruption of healthy microbiotaterm·key (dysbiosis), for example through antibioticsterm·dom, or cytotoxic chemotherapy, can result in loss of colonization resistanceterm·key and increased susceptibility to pathogens [ 15 , 16 ], as well as long-term shifts in microbial ecosystem composition [ 17 ] and other diseases [ 18 ]. Conversely, reconstitution of normal microbiotaterm·key (e.g. through faecal transplants) has been demonstrated to restore host protection, and help cure patients from recurring pathogenic infections [ 19 ]. Such dynamic processes of loss and gain of colonization resistanceterm·key are particularly important in the intestinal tract, making stability, diversity and resilience of our gut microbial ecosystem an intense topic of investigation [ 14 , 20 – 22 ].

Despite the progress in experimentally revealing the beneficial roles of the greater than 100 species (trillions of microbes) populating the human intestine [ 23 ], understanding the microbiotaterm·key from a quantitative perspective is a hugely difficult challenge with difficulties in collecting useful data, and then identifying scalable mathematical models to tackle what can quickly become incredibly complex highly dimensional problems. Within this highly complex and fast-growing research field, a theoretical understanding of how colonization resistanceterm·key emerges as a collective trait from a community network, is maintained as a dynamic process through host lifespan, and reacts to various perturbations remains elusive. With this study, we hope to provide some new tools and concepts to meet this challenge.

In particular, to contribute to this question, we synthesize a new conceptual and mathematical perspective to describe colonization resistanceterm·key at the level of a single host, based on the mathematical quantity of mean invasion fitnessterm·key emerging from a replicator equationterm·mth system for N species.

Modelling invasion resistanceterm·dom based on replicator dynamics The classical N -dimensional replicator equationterm·mth [ 8 ] can be written as follows: where z = ( z i ) i =1, …, N are the frequencies of N competing strategies in a population with ∑ i z i = 1 , f i ( z ) is the fitness of strategy i and f ¯ ( z ) = ∑ i = 1 N z i f i ( z ) the mean fitness over all strategies in the population. Of course, f i may be time- and system-dependent. In the most general form of this equation, there are no constraints on how linear or nonlinear f i ( z ) is as a function of z . The equations describe how the frequencies z i change over time as a function of other strategies in the population; essentially, one strategy will go up in success at any point in time if its fitness is higher than the mean fitness of all strategies present in the system.

Here, we propose the following special version of the replicator equationterm·mth, with a linear fitness function, to model explicit N -species frequency dynamics: with ∑ i z i = 1 . In our formulation, λ i j denote pairwise invasion fitnessesterm·key between species i and j and Θ is a speed constant (typically assumed 1 in the classical replicator equationterm·mth). Hence in this replicator equationterm·mth, f i ( z ) = ∑ j ≠ i λ i j z j and f ¯ ( z ) is the quadratic term Q = ∑ 1 ≤ k , j ≤ N ⁡ λ k j z j z k . The mutual invasion fitnessterm·key, appearing in this fitness function, λ i j , is a key quantity at the heart of adaptive dynamics [ 24 ], and depicts the initial demographic performance, namely the invasion growth rate of a species i in an equilibrium set by a single resident j . Naturally, it is defined only for two distinct species, with λ k k = 0 .

The first feature to remark is that pairwise invasion fitnessterm·key is typically derived from the underlying ecology of the model or details of the competition process between species. The second feature to remark is that every abstract replicator system with a linear payoff matrix A , through basic algebraic operations, can be brought into the special form ( 2.2 ) where the new payoff matrix A ′ has zeros on the diagonal, a transformation that preserves exactly the dynamics, and in addition benefits from the extra attribute that the resulting payoff matrix could be read as a matrix of pairwise invasion fitnessesterm·key between the competing entities.

An emergent exact definition of invasion resistanceterm·dom In this replicator equationterm·mth (equation ( 2.2 )), from the very meaning of the variables and the parameters, the dynamic quadratic term Q = ∑ 1 ≤ k , j ≤ N ⁡ λ k j z j z k (representing f ¯ ( z ) in equation ( 2.1 )), equals precisely mean invasion fitnessterm·key of the system. Indeed, noting that λ k k = 0 and λ ¯ k = ∑ j = 1 N λ k j z j , we have which reduces to In other words, Q sums the global effect of the system on each member, and reflects mathematically whether the multispeciesterm·key consortium as a whole adds ( Q < 0) or detracts from the growth of each individual member ( Q > 0) [ 25 ]. For any entity i in the system, the rate of growth of its frequency ( z i ) relative to other entities is ultimately given by λ i ¯ − λ ¯ , namely by how it performs in pairwise invasion ‘games’ with any other member of the system, and also by how the entire system as a whole ( λ ¯ ) acts to promote or reduce net individual growth. By the same arguments, Q ( z ) is also a quantity that explicitly embodies the invasion resistanceterm·dom of the system towards outsiders.

Invasion by outsiders depends on systemic invasion resistanceterm·dom Mathematically, it is straightforward to see that the initial growth rate of any invader in a system at state z ( τ ) is given by The first summation term in r invader depends on invader traits, λ invader j , namely how this species invades any existing species currently in the system. The second term instead, − Q , is independent of the invader, and crucially describes how the system invasion fitnessterm·key by itself detracts from or adds to invader’s initial growth, if Q > 0 or Q < 0, respectively. One can thus immediately use Q as a global indicator of systemic resistance to invasion by outsider strains, that applies at any point in time, independent of equilibrium or non-equilibrium dynamics.

Matching the generalized Lotka–Volterra and replicator equationterm·mth framework It is a well-known result that every generalized Lotka–Volterra system of dimension n − 1 ( y i , i = 1, …, n − 1 with y = ∑ i = 1 n − 1 y i ) can be written as a replicator equationterm·mth of dimension n ( z i , i = 1, … n ). The equivalence can be shown by the transformation z i = y i /( y + 1) ( i = 1, …, n − 1) and z n = 1/( y + 1) [ 8 ]. Although the two models are topologically equivalent, they are not identical [ 8 , 9 ], and it has been argued that sometimes one formalism can be preferred over the other. One or the other formalism may be more easily connected to data, may be more easily amenable to analysis and hard proofs, or sometimes it may be a matter of convenience, practicality and personal conceptual choice of the modeller for the system of study. In our perspective, the above-mentioned mathematical link is the most important one for bridging between our work here using the replicator equationterm·mth and the studies based on generalized Lotka–Volterra modelling of multispeciesterm·key dynamics.

While discussions are old and ongoing in the literature about the merits and limitations of the Lotka–Volterra approach [ 26 – 29 ], one of which being the lack of explicit characterization of the underlying limiting resource in the system, further comparison between the Lotka–Volterra frameworkterm·mth and the replicator equationterm·mth falls beyond the scope of this paper, but see [ 8 – 10 ].

For us, it is natural to look at the replicator equationterm·mth ( 2.2 ), how it arises by itself in an ecological context, noting that the key quantities modelled are not the same as those in the Lotka–Volterra equations, and the meaning of the parameters involved is not the same. Most crucially, invasion resistanceterm·dom of the system emerges naturally within the replicator equationterm·mth, while it is not common to see mean invasion fitnessterm·key of a system explicit in a typical Lotka–Volterra frameworkterm·mth.

Invasion fitnessesterm·key from first principles and specific models The pairwise invasion fitnessterm·key coefficients are at the core of the replicator dynamics (equation ( 2.2 )) and reflect initial demographic performance in dyadic i → j invasion scenarios between two species, starting from low frequency of i in an equilibrium set by j [ 24 ]. There are at least two key challenges regarding invasion fitnessesterm·key: (i) the first one is how to derive them bottom-up mechanistically from the details of the underlying ecology and biology of species dynamics and interactions; (ii) the second one is how to make strong qualitative and quantitative predictions about their effect on global dynamics once their structure and magnitudes are known.

In this paper, we address both of these challenges. First we show one possible route for the derivation of λ i j , bottom-up, exploiting an analogy with an SIS transmission dynamics framework proposed for multiple strains [ 25 , 30 , 31 ]. Application of this model to new contexts, for example, to in-host microbiotaterm·key dynamics, requires re-interpretation of original variables and parameters, and a conceptual shift of scales. In §§3.1–3.3, we present the key ingredients for such re-interpretation, and develop the analogy for translation of the model to the within-host level.

Secondly, in §4, we illustrate how this model can be used and which predictions can be made based purely on knowledge of the Λ matrix structure between N species.

The analogy: from SIS to within-host microbial dynamics A host is a system of ‘free’ and ‘occupied’ micro-niches for colonization . We define the host as a well-mixed system, and the dimensionality of the system as the number of similar microbial types (entities) interacting in colonization of such a system, hereafter denoted as species. In the re-framing following [ 25 ], we model the within-host environment as a number of potential micro-‘niches’ depicting generic growth units (resources) to be used, and we track the proportion of such niches that remain free or susceptible ( S ), those that get singly colonized by any species ( I i ), and those that get co-colonized, either twice by the same species ( I ii ) or by two different species ( I ij ) (see §3.2, and figure 1 a ). By taking suitably small such growth units within-host, we can truncate the multiplicity of ‘infection’ (MOI) to at most two co-colonizing species per unit. It is assumed that such growth units (similar to hosts in the epidemiological setting) homogeneously mix in the system, and when they are free, they are equally accessible to propagules of all species, whereas when they are already colonized, they become more or less available to co-colonization/co-utilization by others, depending on the identities of the interacting partners.

A conceptual analogy for within-host microbiotaterm·key dynamics (here illustrated for N = 2). One mechanistic path to arrive at the replicator equationterm·mth (equation ( 2.2 )) is from the population structure model, based on epidemiological SIS dynamics with co-infection [ 25 ]. ( a ) Conceptualization of the within-host environment, where micro-compartments can be free (uncolonized S ), singly colonized ( I i ) or doubly colonized ( I ii , I ij , I jj ) by microbes. Full equations are provided in §3.2. ( b ) Model structure for compartmental state transitions, including colonization (black arrows), co-colonization (colour-shaded arrows) and clearance (dashed arrows). Natural birth/death is denoted by green arrows. ( c ) Singly colonized compartments by one species ( i ) may have a reduced or increased susceptibility to becoming co-colonized by a second species ( j ), noted by coefficients K ij below or above 1. For similar species, we may write: K i j = k + ϵ α i j , with the parameter 0 < ϵ ≪ 1 being small enough. ( d ) For N = 2, there are four outcomes ( table 1 ), depending on species interactions: exclusion of 1 or 2, coexistence, or bistabilityterm·dom. This is mathematically encoded in the mutual signs of the invasion fitnessesterm·key λ i j , whether (+ , −), (− , +), (+ , +) or (− , −). ( e ) Context-dependence of mutual invasion fitnessesterm·key, in this mechanistic model, where λ i j = α j i − α j j − 1 k ( R 0 − 1 ) ( α i j − α j i ) . If k or R 0 change, with α ij fixed, the pairwise invasion fitnessesterm·key change and generate different ecology between 2 species. In this example, starting from coexistence (shaded region with both λ 1 2 and λ 2 1 positive), lowering overall colonization ( R 0 ↓ ), e.g. via antibioticsterm·dom, or increasing it drastically, may lead to competitive exclusionterm·dom of an opposite species in each direction. Going from stable coexistence of 2 species to competitive exclusionterm·dom with a single winner reduces invasion resistanceterm·dom of the system from a positive value to 0. In an N -species network such shifts can be more complex but yet predictable [ 30 ].

Within each single growth unit, microbial propagules can grow, then upon release, can be transmitted to other units. Microbial propagation dynamics, thus, under homogeneous mixing of propagules in the in-host milieu, can be described, to a first-order approximation, through simple mass-action kinetics ( figure 1 b ). The growth of each microbial species depends on their transmission rate β and frequency of niches emitting propagules of that species. The transmission parameter, β , in this analogy, can be taken as a net parameter, encapsulating the chain of events including local growth happening in the micro-niche of ‘origin’, followed by migration, and diffusion of infectious propagules in the system, multiplied by the effective probability of settling in other ‘destination’ niches. In its basic version, the model assumes β is equal for all i , yielding reproductive fairness in the system (but see [ 31 ] for a formal proof and generalization to a higher number of traits varying among relatively similar species).

Interactions among species happen upon co-colonization, where the colonizer ( i ) and co-colonizer ( j ) microbes may compete or facilitate each other, each ordered pair with its own coefficient K ij ( figure 1 c ). We do not specify the underlying explicit ecological mechanisms, which may be synergistic or antagonistic, and include cross-feeding, antagonistic warfare via bacteriocins, local nutrient alteration or depletion, acidification of the micro-environment, substrate availability and others [ 32 – 37 ]. What counts essentially at an abstract level is that the micro-scale presence of one species alters the value of that same resource-growth unit complex for individuals of that same/other species, making it somewhat better or worse. Such altered susceptibilities to co-colonization K ij ( > 1 : facilitationterm·dom , < 1 : competition ) comprise a large potential spectrum ( figure 1 c ), and for N = 2, form a 2-by-2 matrix which fully determines one of four possible outcomes between 2 species [ 38 ]: exclusion of 1, exclusion of 2, stable coexistence or bistabilityterm·dom ( figure 1 d ). Bistabilityterm·dom means that one species or the other may be the sole winner in the system depending on initial frequency. However, the co-colonization interaction matrix in general is higher-dimensional N × N , when considering N possible ecological partners in the community. Singly and co-colonized units are assumed to produce an equal number of propagules, and mixed co-colonization units I ij produce and transmit i / j propagules with probability 1/2 (a symmetry assumption that is relaxed in subsequent model extensions [ 31 , 39 ]).

Free growth units are replenished at rate r (in-flow), assumed equal to the natural decay rate (out-flow), keeping thus constant the carrying capacity of the system. Within each available micro-niche, microbes can grow for a mean duration of 1/ γ time units, hence yielding a turnover rate of m = r + γ . The basic reproduction number R 0 denotes how many new colonizations each colonized micro-niche will produce over its lifetime as occupied, if occurring in a totally colonization-free environment: R 0 = β / m , and if R 0 > 1 there is a globally stable endemic state with 1 − 1/ R 0 colonization prevalence, a classical feature of SIS models. If R 0 < 1, the system is colonization-free.

Overarching dynamics of colonization and co-colonization In the original model for N species, described in detail in [ 25 ], we track the variables S , I i and I ij , that here refer to the proportion of uncolonized niches, those singly colonized by species i and those co-colonized by species i and j . The dynamics of colonization and co-colonization are given by the following equations: where F i = β ( I i + ∑ j = 1 N ( 1 / 2 ) ( I i j + I j i ) ) gives the force of propagation of species i in the system, summing the contribution of all micro-niches emitting propagules of that species. Denoting by I = ∑ i I i the total prevalence of single colonization, and D = ∑ i , j I i j , the total prevalence of co-colonization, we have: S + I + D = 1. In this model in micro-niches co-colonized with different species i and j , the opportunity of each species to grow and be propagated is equal, thus represented with the 1/2 probability. In the system, the composite parameter m = γ + r represents the growth unit turnover rate, encapsulating both the clearance rate γ of each colonization episode, after which growth units become immediately available for new colonizations, and the natural replenishment rate of new universal units r . It is assumed that natural replenishment is balanced by natural out-flow rate r = d . Species interactions via local environmental modification are represented by the matrix K , where values K ij above 1 indicate pairwise facilitationterm·dom, while values of K ij below 1 reflect inhibition between propagules of i and j , and where the order of colonizer and co-colonizer landing in a growth unit matters.

From population structure to replicator dynamics. The mathematical step that connects the population structure of the colonization with the frequency of individual species is time-scale separation [ 25 , 31 ], based on some similarity arguments for between-species traits.

In the most basic model [ 25 ], with variable co-colonization interaction coefficients K ij , the features of the full slow–fast approximation are summarized in table 1 . Species similarity in traits K ij means with ϵ suitably small (e.g. the standard deviation of K ij ) and k a benchmark, e.g. the mean of K i j : = ∑ 1 ≤ i , j ≤ N K i j / N 2 . With this formulation, the full dynamics of the system can be decomposed into: fast dynamics, which are neutral and serve to equilibrate global population structure of colonization and co-colonization variables (see S , I , D in table 1 ), and slow dynamics, which are effectively the N -species frequency dynamics over the timescale τ = ϵ t , captured by a replicator equationterm·mth under some particular values of λ i j .

Key features of the replicator dynamics emerging from a co-colonization model for similar species [ 25 ]. But see [ 31 , 39 ] for the model with 5-trait dimensions where species can vary.

What is important in this separation of time-scales is that global variables from the fast dynamics stay constant during the slow dynamics, and can play a role in determining the magnitudes of the invasion fitnessesterm·key between species. During the slow dynamics λ i j do not change, only species frequencies change obeying equation ( 2.2 ). This feature is robust also in extensions of the model [ 31 , 39 ].

How species traits inform invasion fitnessesterm·key λ i j Special case: species vary only in K ij . In the model with variation only in co-colonization coefficients among species [ 25 ], species pairwise invasion fitnessterm·key is given by and is a direct function of the phenotypic variation in co-colonization interactions (the α ij ’s), and of the ratio of single to co-colonization in the system: What is key in this λ i j is that the invasion performance of i into a j - only equilibrium depends on how much the resident entity j relatively favours the invader i over itself (trait difference α ji − α jj in λ i j ) and also on the net benefit of i from transitions to mixed co-colonization with j (− μ ( α ij − α ji ) term) [ 25 , 30 ]. This second term depends on relative availability of singly colonized micro-niches to co-colonize in the system, hence on the quantity μ .

Thus, there are two effects in this trait-based λ i j : (i) comparative advantage from non-self to self-interaction of the resident, and this is not dependent on co-colonization prevalence in the system and (ii) relative benefit from asymmetric investment in mixed co-colonization, and this decreases with prevalence of co-colonization in the system.

Notice that in this model, whether the species are competing or facilitating each other in their common growth on average ( k → 0, or k large) has an effect on λ i j via the parameter μ , which is an explicit decreasing function of R 0 and k . The dependence on μ happens only for non-symmetric co-colonization interactions, i.e. only if α ji ≠ α ij and is stronger if the α ′ s are more anti-symmetric. For fixed α ’s, the linear dependence of λ i j on μ , implicating R 0 and k , exemplifies an elegant contextual feedback on pairwise invasion fitnessterm·key ( μ high or low), which is key to understand analytically special cases of the replicator dynamics in these two extremes, and as a consequence, the limit cases of system invasion resistanceterm·dom [ 30 ].

General case: species are more complex and vary along five fitness dimensions. More generally, in a scenario when more traits are allowed to vary between species [ 31 , 39 ], the trait-based pairwise invasion fitnessterm·key expression is more complex, but the same replicator dynamics can be derived. We outline in appendix A that the replicator dynamics (equation ( 2.2 )) are general, and can be applied also to cases where species have different proliferation rates, clearance rates, and biases in local growth and priority effects; the pairwise invasion fitnessterm·key in that case becomes a weighted average of the asymmetries along all phenotype axes (see eqn (2.22) in [ 31 ]).

Whatever their exact expression, together the λ i j between all pairs drive collective dynamics in the system, and the definition of system invasion resistanceterm·dom via Q = λ ¯ remains robust.

Using the replicator equationterm·mth to study system invasion fitnessterm·key Invasion resistanceterm·dom: system versus outsider species It is clear from equation ( 2.2 ), that to survive in the system, an individual/species must always ‘play games’ with co-occurring partners, and the payoff depends not only on an individual’s own (vector) strategy ( λ i j , j = 1 , . . . , N ), but also on strategies and abundances of its opponents. Because Q sums products of frequencies, it is non-zero only when there are more than 1 species coexisting, and it is exactly zero when λ j k = − λ k j for all pairs ( j , k ). In the simplest system with only two species, N = 2 ( table 1 ), it is easy to see that at the coexistence steady state, Q ( z 1 , z 2 ) is positive and maximized for z 1 = z 2 = 1/2. Thus, a 2-species system is mostly protected against a third invader, when its two constituent members are equiabundant ( figure 2 ).

System invasion resistanceterm·dom Q in the N = 2 system with 2 species coexisting. Here we visualize a three-dimensional plot of Q ( x , y ) = xy /( x + y ) as a function of mutual invasion fitnessesterm·key x = λ 1 2 and y = λ 2 1 when the system is at the coexistence equilibrium with species frequencies ( z 1 ∗ , z 2 ∗ ) = ( λ 1 2 / ( λ 2 1 + λ 1 2 ) , λ 2 1 / ( λ 2 1 + λ 1 2 ) ) . The maximum of Q is reached at λ 1 2 = λ 2 1 , i.e. for equiabundance of two resident species, and it increases with the absolute magnitude of the λ ′s.

Depending on how many species make up the resident species pool, N , and how they interact, one models a different system size and a different network of mutual invasion fitnessesterm·key (directed edges with precise magnitudes). This will generate system-specific frequency dynamics eventually leading up to a pruned community with n ≤ N species, and simultaneously the evolution of system resistance to invasion by outsiders ( table 1 ).

We tested several invasions numerically with model simulations, and how the success of an invader’s initial growth depends on system features. When comparing different multispeciesterm·key systems in their resistance to invasion, for fixed parameters, we find that the higher the value of Q , the harder the invasion by outsider species ( figure 3 a – g ). By contrast, we find no significant role of system diversity (Shannon, Gini–Simpson indices; figure 3 i – k ), but a positive effect of the number of species for invader success ( figure 3 h ). The final outcome of any given invasion however, as captured in equation ( 2.4 ), besides system Q , depends also on other details of resident dynamics, and crucially on the traits of the invader, relative to the resident community, at the invasion time. This does not necessarily imply restrictions on cooperation or competition [ 40 ], but simply that the mean invasion fitnessterm·key of the invader must exceed the mean invasion fitnessterm·key of the resident system, whatever the path to satisfying such criterion may be.

System invasion resistanceterm·dom, as defined by our model, and outsider invasion. In ( a – d ), we fixed the multispeciesterm·key K ij interactions in susceptibilities to co-colonization in two systems with N = 6, β = 1, k = 0.3, R 0 = 2 but with different underlying species dynamics. We simulated 100 realizations of invasion by an outsider species, sampling its interaction traits from the same distribution as the resident N ( k , ϵ 2 ) , and introducing it at the mid-point of the interval. We counted invasion successful if the invader grew from its low initial frequency (here 10 −4 ). System 2 (higher Q ) was superior in preventing invasion to system 1: experiencing only about 30% successful invasions compared to more than 60% successful invasions in system 1. In ( e ), we repeated the procedure for 100 randomly generated systems (blue dots) and examined system Q at invasion and invader success, finding a significant negative relationship between invasion resistanceterm·dom and invasibility ( p < 10 −5 ). The final outcome of invasion whether the invader ultimately coexists, or how original species are affected was not studied. In ( f – k ), we examined invasion success for N = 20. Simulating 100 random systems as in ( e ), we regressed outsider invasion probability on different system characteristics: Q ( p < 10 −5 ) in ( f – g ) and the number of species in the system ( p < 0.01) in ( h ) were significant, while diversity indices ( i – k ) were not.

Invading one system under different conditions: effect of global environmental perturbations When the invasion fitnessesterm·key in the replicator equationterm·mth have a mechanistic derivation (e.g. equation ( 3.2 )) where a dependence with global quantities is immediately apparent (e.g. μ = 1/( R 0 − 1) k ), it is possible to study in some detail the system’s response to a change in such mean-field parameter. In the special case of a microbial ecosystem where invasion fitnessesterm·key are defined by ( 3.2 ), perturbations affecting global context, such as basic reproduction number R 0 and mean interaction in co-colonization k , directly change the mutual invasion fitnessesterm·key λ i j [ 25 , 30 ]. In highly diverse communities (many species N > 2, and/or more traits varying), effects of context-altering perturbations can be very complex, as we have shown in some cases [ 30 , 31 , 39 ].

However, among N species varying randomly only in co-colonization interaction phenotypes (susceptibilities K ij ), a key modulator of coexistence regimes is the ratio of single- to co-colonization in the system, μ [ 30 ], decreasing with mean growth potential R 0 and mean facilitationterm·dom k . Perturbing these mean-field parameters can have predictable effects on invasion resistanceterm·dom of the system.

How antibioticsterm·dom alter colonization resistanceterm·key of a microbial ecosystem: R 0 effect. One such perturbation example could be broad-spectrum antibioticsterm·dom, leading to alteration of R 0 in the system, likely reducing the colonization of all microbial species in the system, and shifting their internal competitive balance. Even when acting in the same way on system members, increasing their clearance rate or reducing proliferation, more antibioticsterm·dom in our mechanistic model can drive the system from coexistence to competitive exclusionterm·dom (see the case of N = 2 in figure 1 e ), hence toward Q = 0, promoting outsider invasion. Through our explicit λ i j [ 25 ], one can clearly see that this R 0 impact drastically varies, depending on the interaction network between species; for example it is stronger when k is lower, i.e. when competition is higher in the system (see electronic supplementary material, S3–S4). In particular, even when starting at the same baseline invasion fitnessesterm·key between two species (i.e. also the same coexistence scenario in terms of relative abundances), the stronger the mean competition in micro-scale co-colonization (the lower k ), the more drastic will be the selection in the system induced by the antibioticterm·dom. Such explicit individualization of perturbation effects, based on underlying interactions, may be easily related to data and to context-dependence of outcomes in microbial consortia.

How species facilitationterm·dom alters colonization resistanceterm·key of a microbial ecosystem: k effect. In [ 30 ], we have shown in detail how in this model, μ = 1/(( R 0 − 1) k ) acts as an axis of tunable contextual feedback on the multispeciesterm·key dynamics: for standard normally distributed rescaled interaction strengths ( α ij ), when μ → ∞, and all else is kept fixed, the mutual invasibility network tends to contain mostly competitive exclusionterm·dom edges, which makes more species coexist, but in unstable and oscillatory fashion, contrasting the regime of μ → 0 displaying mostly stable steady states. The extreme of μ → 0 can be obtained for increased facilitationterm·dom between resident species (higher values of k ), which typically would lead to higher invasion resistanceterm·dom of the system. This is similar to a result about cooperative microbiotaterm·key obtained via simulations in other studies [ 40 ]. In our model, the phenomenon is entirely predictable analytically, and a similar effect could be obtained with amplifying R 0 , all else fixed. Conversely, decreasing k , hence decreasing facilitationterm·dom, in the multispeciesterm·key system the ratio of single- to co-colonization μ would increase, thereby promoting lower and oscillatory invasion resistanceterm·dom via fluctuating multispeciesterm·key dynamics [ 30 ].

As the pivotal quantity in λ i j is the product between R 0 − 1 and k , an equivalent effect could be obtained by decreasing R 0 in such a multispeciesterm·key scenario. Lower overall colonization rate and higher μ could lead possibly to a fluctuating Q and could make a given system more vulnerable to opportunistic invasion, particularly when Q hits low or negative values, even if transiently. This seems to mimic what has been observed empirically, for example in individuals who following antibioticterm·dom therapy experienced domination of their gut microbiotaterm·key by pathogenic E. faecalis or C. difficile bacteria, and succumbed to infection [ 15 , 41 ].

The role of pairwise invasion structures among N species Above we considered a case of mechanistic knowledge of λ coefficients and their context-dependence, and how this would make them change with global parameters, leading to a predictable shift in system invasion resistanceterm·dom. What can we say about the species dynamics and invasion resistanceterm·dom of a system when the Λ matrix is quantitatively known but without explicit contextual drivers? How can we compare two systems that simply vary in their Λ matrix structure? We have previously sketched several canonical mutual invasion structures for the matrix of λ i j in a system with N members [ 25 ], namely, invader-driven (equal columns); resident-driven (equal rows); symmetric ( Λ = Λ T ); anti-symmetric ( Λ = − Λ T ); and random.

In the special mechanistic case ( 3.2 ), these structures could be mapped to special assumptions about the co-colonization interaction matrix with K ij . For example, symmetric K , such that K ij = K ji = K ii + K jj − k leads to invader-driven Λ ; whereas symmetric K , such that K ij = k for i ≠ j , with the diagonal entries free to vary, leads to resident-driven Λ (see table 2 in [ 25 ] for more details).

However, there are in principle infinitely many possible microscopic K ij or other trait combinations and model details, leading eventually to the same macroscopic λ i j architecture. That is why next we revisit these Λ structures, more generally, from a top-down perspective, to explore their direct implications for invasion resistanceterm·dom in a completely closed multispeciesterm·key ecosystem.

For symmetric mutual invasion ( λ i j = λ j i ), it can be shown that Q always increases (d Q /d t > 0), a case related to Fisher’s fundamental theorem [ 42 ] (for details see electronic supplementary material, section 4 in [ 25 ]). In such instance, cycles among species are impossible, and their dynamics converge typically to stable coexistence fixed points, and competitive exclusionterm·dom is unlikely.

In the opposite extreme, for anti-symmetric mutual invasion ( λ i j = − λ j i ), the system colonization resistanceterm·key is exactly zero over all time, Q = 0, corresponding to zero-sum games in evolutionary game-theoretic terms [ 43 ]. Dynamics in that case tend to a (structurally unstable) family of cycles around a centre, like in Lotka–Volterra prey–predator models [ 44 ], where an odd number of species coexist. More special cases, such as the species-centric invader-driven and resident-driven invasion, and explicit K – Λ matrix links are analysed in [ 25 ].

Investigating community dynamics under these pairwise invasion architectures, we find that keeping all else equal, and starting from random initial frequencies, Q increases fastest over time in a system regulated by invader-driven mutual invasion ( figure 4 ), i.e. where species are defined by their active invasiveness trait. This case typically leads to coexistence among several species and rarely to exclusion. In the resident-driven mutual invasion case, where species are defined by how they allow invasion by others, hence by their invasibility, dynamics tend most likely to competitive exclusionterm·dom, i.e. toward Q = 0, but following an initial selection period during which systemic Q is negative.

Mutual invasion structure of the species in the ecosystem, dynamics and system mean invasion fitnessterm·key. ( a ) Illustration of a mutual invasion network between 10 species, whose structure and pairwise invasion fitnessterm·key values lead to multispeciesterm·key dynamics on the right (equation ( 2.2 )), with ultimately five species coexisting. ( b ) Colonization resistanceterm·key for special invasion structures. On the left, we plot Q dynamics (mean: lines, ±s.d.: shading) over 30 stochastic realizations for six cases of Λ between N = 30 species (shown on the right). For each case, random uniform λ i j ∈ [ − 1 , 1 ] , from the same distribution and range, are drawn and fed each into equation ( 2.2 ) (assuming Θ = 1 ), which are then solved with random initial relative abundances z i (0), as in [ 25 ]. Symmetric Λ yields dynamics described by Fisher’s fundamental theorem, where mean invasion fitnessterm·key always increases and stable coexistence is likely. Invader-driven Λ implies large coexistence potential among multiple species. Resident-driven Λ implies large potential for competitive exclusionterm·dom, and favours a single species becoming the only persistent one of the community. Antisymmetric Λ makes Q ≡ 0, corresponding to zero-sum games among species and complex structurally unstable coexistence dynamics. Almost-antisymmetric Λ typically produces complex but more regular dynamics, for example limit cycles, leading to periodic Q . Random mutual invasion allows rich multispeciesterm·key dynamics to unfold. Cases a–c can be obtained for three particular instances of symmetric co-colonization interactions matrix K [ 25 ].

A vast complexity of possible co-colonization interaction networks and emergent invasion architectures between N species may yield coexistence scenarios of more complex nature than fixed points, including multistable attractors, limit cycles and chaos [ 30 ]. Yet, in all these cases, the net invasion resistanceterm·dom that arises and changes as an adaptive trait of an evolving ecosystem can be mapped directly to the relative frequencies of constituent species. Thus, depending on how frequency dynamics unfold, subject to intrinsic or extrinsic drivers, Q dynamics can be higher or lower, monotonic or oscillatory, and feed back on the system, in a fully analytically explicit manner ( table 1 ).

Invading one system at different times A corollary of this mathematical understanding of invasion is that following the dynamics of systemic mean invasion fitnessterm·key Q ( τ ) = λ ¯ , the timing of invasion (by the same invader) will likely have a different success and effect on a given system, depending on the actual magnitude of Q . We performed simulations of invasion in two qualitatively different systems to illustrate this phenomenon: (i) in a system where multispeciesterm·key coexistence tends to a stable fixed point ( figure 5 a , b ) and (ii) in a system where underlying species coexistence is oscillatory ( figure 5 g , h ). In figure 5 c , d , we show one illustration of dynamics of a particular successful invader in multispeciesterm·key system 1 (invader coloured in brown), and the corresponding dynamics of systemic invasion resistanceterm·dom Q after invasion. In figure 5 i , j , we show one illustration of dynamics of another particular successful invader in multispeciesterm·key system 2 (invader again coloured in brown), and the corresponding dynamics of systemic invasion resistanceterm·dom Q after this given invasion.

Two qualitatively different resident systems and sensitivity of invasion outcome to invasion timing. Here we fix two resident systems with N = 5 ( a–f and g–l ) and study the effect of invasion time over 200 random invaders. ( a , g ) The baseline dynamics of species frequencies in the absence of invasion (system with stable steady state and oscillatory dynamics, respectively). ( b,h ) The dynamics of mean invasion fitnessterm·key of each system in the absence of invasion. ( c,i ) Example of species frequencies for a particular invader (in brown) and a particular invasion time (dashed line). ( d,j ) Corresponding mean invasion fitnessterm·key dynamics for each system with invasion shown in ( c , i ). ( e,k ) Summary of many invasion outcomes as a function of invasion timing in the same system. We sample 200 invaders for each time point, simulate invasion starting from low initial frequency of 1%, and compute the proportion that results in initial growth (local invasion) and probability of ultimate persistence in the system (global success) using a dynamic horizon of T = 1000. In a system with monotonic Q ( τ ) = λ ¯ , mean invader success tends to a constant, whereas in the system with oscillatory λ ¯ invader success tends to oscillate. ( f,l ) One aspect of system outcome with and without invasion is the final number of coexisting species (here calculated as those above a given relative abundance threshold z i > 1 % during the last 10% of the simulation time interval, and averaged over many different simulations). On the left, we observe that invasion may increase slightly number of coexisting species if invasion happens early in the course of community dynamics, when the system is far from equilibrium. If invasion happens closer to equilibrium, the qualitative number of final coexisting species tends to be the same, but with one of the species effectively replaced by the invader. On the right, we observe that invasion on average reduces diversity slightly below its originally expected level.

Then we repeat such invasion experiments over different invasion time points, and multiple random invaders. First, we confirm that temporal success of outsider invasion on average follows clearly the dynamics of Q of the underlying system: it tends to a constant when Q is monotonic ( figure 5 e ), and oscillates periodically when Q is periodic ( figure 5 k ). In particular, the less Q changes (i.e. when the underlying system tends to stable coexistence fixed points), the weaker the timing effect, except for a short transient whenever the resident community may be far from its equilibrium ( figure 5 e ), and it may seem that the probability of successful invasion decreases with invasion time, until it saturates. Conversely, the more oscillatory invasion resistanceterm·dom Q is (i.e. the system is characterized by limit cycles), the stronger the timing effect on average, and invasion success can vary more along time ( figure 5 k ), albeit being overall relatively higher.

These patterns seem to persist when computed both locally (invader growth after invasion) and globally (invader persistence after invasion), with the relative difference between the two correlated probabilities being about 20%. This suggests that on average the final fate of invaders may be quite predictable based on their initial invasion dynamics: if they invade, most likely they will persist (70–80% of invaders that do experience some growth in the system, will persist).

Our simulations also suggest that mean invasion fitnessterm·key tends to increase after invasion ( figure 5 d , j ), a plausible result when considering that a prerequisite for successful invasion is λ ¯ invader > Q = λ ¯ system , so if an invader eventually manages to stay in the system it does so via augmenting systemic Q . Under stable coexistence, the successful invader tends to replace one species in a given niche—thereby not affecting final number of coexisting species n ( figure 5 f ), but increasing the stability of this coexistence. In the oscillating case, such increase in stability comes from the successful invader more often decreasing the number of coexisting species ( figure 5 g ).

Invading qualitatively different systems Random invader success can also be studied in given community contexts, fixing invasion time, but now varying systematically the underlying mutual invasion matrix Λ of the resident species ( figure 6 ). In other words, how does the probability or random invader success depend on the ‘quality of the invasion game’ among the residents? For this, we performed 50 simulations of random outsider invasion, for 200 realizations of each canonical Λ structure (invader-driven, resident-driven, symmetric, anti-symmetric and random), and computed two indicators of successful invasion: (i) the probability of initial invader growth (local invasion) and (ii) the probability of final invader persistence in the system (global invasion).

Mutual invasion structure of the resident microbial consortium and random invader success. We generated 200 random Λ matrix structures of size N = 5, with entry coefficients in the range [−1, 1], for each of six canonical architectures: symmetric, invader-driven, resident-driven, anti-symmetric, almost anti-symmetric and random. For each random invader (random uniform λ invader j , j = 1 , . . . , 5 ), we simulated invasion in 200 systems of each type ( Λ matrices), starting at the same initial frequency of the invader of 1% and equal relative frequencies of the resident five species totalling 99%. Over 200 systems, we computed the probability of outcome for each invader, summarized in the coloured boxplots, which show entire variation over 50 invaders for each type of Λ matrix. Systems were run for a time horizon of T = 1000. The invader was introduced at two timepoints: ( a,b ) early in the development of the system close to initial equiabundance ( τ = 10), and ( c,d ) later in the development of the system, closer to its intrinsic equilibrium ( τ = 500), which contains more signature of Λ -driven selection. Local invasion counts all those instances where the invader grew from its initial frequency at some point from introduction up to T = 1000. Global invasion counts all those instances where the random invader is still present in the system at final time T = 1000. The ranking observed between Λ structures can be very useful when engineering multispeciesterm·key communities for the purpose of maximizing colonization resistanceterm·key.

As expected from our previous analysis [ 25 ] and the dynamics of Q ( τ ) in figure 4 , we find that the mutual invasion structure in the resident community strongly impacts invader success. Importantly, the highest probability of invasion, both local and global, is obtained for the anti-symmetric invasion architecture of the resident species, yielding Q = 0 as the most favourable invasion context for outsiders. By contrast, the lowest probability of invader growth and invader persistence is found for symmetric, invader-driven and random mutual invasion of the resident species, indicating that these communities contain the least favourable (i.e. highest Q ) colonization resistanceterm·key for outsider invasion. This ranking is robust to random variation in the resident (200 different realizations of the same structure) and invader invasion coefficients (50 different invaders), pointing to a key signature of system stability and invasibility encoded precisely in the Λ matrix type. Although if introduced early ( τ = 10), the probability of outsider species initial growth is lowest in a community characterized by invader-driven (row-only variation) Λ architecture ( figure 6 a ), ultimate outsider persistence in the system ( figure 6 b ), and in general both outsider growth and persistence for later invasions are minimized if the underlying system displays symmetric mutual invasion structure (blue boxplot in figure 6 c , d ).

Crucially, the ranking of these mutual invasion structures ( Λ matrix types) suggests that the most resistant, externally stable communities are those whose constituent species display symmetric or hierarchical invader-driven pairwise invasion. Such Λ matrices, in our trait-based model, may be both obtained via special cases of symmetric co-colonization susceptibilities K ij (see table 2 in [ 25 ]). This suggests that stronger collective coexistence may indeed arise when species pairs derive variable but symmetric benefits from shared resources.

Ultimately, the type of replicator dynamics critically depends on the species invasion matrix, whether the λ i j are derived from particular biological traits and processes [ 25 , 31 ], or represented via qualitative structural constraints. While capturing all the realistic complexities and feedbacks between species traits, context, and systemic invasion resistanceterm·dom may be difficult at best, the conceptual perspective and results presented here should provide a new key to study system invasibility: precisely the architecture and magnitudes of mutual invasion coefficients among constituent species. This opens an avenue for exciting future developments in the study of invasion resistanceterm·dom using the replicator equationterm·mth.

Discussion Colonization resistanceterm·key in microbiology [ 14 ] is closely linked with the theory of invasibility in ecology [ 11 , 45 , 46 ], long studied but rarely integrated between the disciplines. Despite our initial motivation, over the course of this study, it became increasingly clear that what we are advocating for goes in fact beyond microbiology, and constitutes an approach to study invasion resistanceterm·dom in multispeciesterm·key systems via the replicator equationterm·mth formalism. The microbiotaterm·key analogy developed is useful to connect the dots from biological details to invasion fitnessterm·key and replicator dynamics, and serves as an illustration of a very relevant context of potentially immediate applications, where this approach to study system invasibility can be adopted.

In the late 1950s, Elton argued that complex communities should be more resistant to invasion by new species [ 47 ]. Later May [ 3 , 4 ] showed that complex ecological communities tend to be less stable, i.e. in these systems it is harder to return to equilibria from small perturbations to existing species. Since then, the complexity–stability debate began, inspiring many investigations [ 11 ]. Crucial to study invasibility and stability of ecosystems have been simulations, using random interaction networks or varying their properties, e.g. the percentage of cooperative versus competitive links [ 48 ], or examining the relative interaction in invader versus within-resident community [ 40 ]. More recently, theoretical approaches are also invoking trait-based invasion fitnessterm·key concepts [ 12 ].

Two notions of stability of a multispeciesterm·key system, internal versus external, have been studied, with the latter receiving much less attention. It is external stability that relates most with invasion resistanceterm·dom. While the stabilizing versus destabilizing role of mutualism and cooperative interactions among species is still theoretically debated [ 48 , 49 ], as are diversity effects [ 50 ], recent work is attempting to reconcile Elton’s and May’s views [ 51 ] and seeking deeper syntheses. However until now, to our knowledge, no unifying formal analytical quantity has emerged for invasion resistanceterm·dom per se , that integrates species interactions along the competition–cooperation spectrum, their traits and resource availability.

In this paper, by calling attention on the replicator equationterm·mth formalism, we propose a fresh conceptual framework to address invasibility, and in particular, to contribute to the formalization of colonization resistanceterm·key concepts in microbiology [ 52 ]. Despite being the central equation of evolutionary game theoryterm·key [ 8 ], the replicator equationterm·mth is surprisingly little used in ecology or microbiology, in favour of the predominant use of generalized Lotka–Volterra frameworksterm·mth for total abundances [ 5 , 7 , 49 ]. In our setup, we favour the use of the replicator equationterm·mth because it precisely models the frequencies of multiple species, which are typically more readily available empirically, e.g. only relative abundances between (pseudo-) species of interest are most often empirically accessible via metagenomic sequencing data. In our model, chosen a ‘species’ resolution level, these relative abundances can be fully predicted via replicator dynamics.

Central to our replicator equationterm·mth (equation ( 2.2 )) is the measure of mutual invasion fitnessterm·key between any two species, λ i j , and the various topologies of the matrix or contextual dependencies it generates. This higher-level trait inevitably arises as a function of lower-level traits and species asymmetries, including cooperation/competition (in our particular case co-colonization susceptibilities K ij ), but also mean-field resource availability indicators (in our case R 0 , k )—dependencies which can be sometimes known or derived a priori from an underlying biological model. Within the replicator equationterm·mth, these pairwise invasion fitnessesterm·key undergo selection and evolution in their own right following species frequency dynamics. Such evolution can be mapped to a collective mean system trait, Q = λ ¯ , applicable to both equilibrium and non-equilibrium scenarios. This replicator framework unifies ecological, evolutionary dynamics and emergent system properties in multispeciesterm·key communities. The same spirit of reconciling simplicity with complexity in functional biology and of improving predictability of models has also been advocated by others [ 53 , 54 ].

By illustrating a mechanistic path for such replicator equationterm·mth, through an analogy with an epidemiological model [ 25 ], we simplify the complexity of multispeciesterm·key dynamics, while including an abstract growth unit, integrating cooperation and competition in the same continuum via the tunable parameter k , and showing how the magnitude of the invasion fitnessterm·key [ 24 ] can be bottom-up derived and trait-based. Other particular mechanistic or non-mechanistic realizations of such invasion fitnessterm·key matrix between species could be studied in the future.

The whole-community behaviour exhibited in our model is a generic consequence of dynamic resource availability, similar to [ 55 ], where, in our case, abstract resources are amplified or diminished in value, by the simple co-presence of other ecological partners. In our closed system, essentially, species are simultaneously the ‘active’ entities, striving to survive, and the very ‘substrates’ upon which their life depends. There is no requirement for cooperative interactions in order for Q to evolve as a community-level function; only the possibility of species micro-scale co-occurrence, k > 0, is needed, coupled with the assumption that organisms transform their local environment in slightly different ways. Then it becomes another application of the basic notion that fitness is context-dependent [ 8 , 56 – 58 ].

Many theoretical phenomena, like those explored above and beyond, directly accessible via this formalism, could be empirically matched when comparing gradients in microbiotaterm·key dynamics and composition, caused by antibioticsterm·dom, diet or other stressors in multispeciesterm·key communities, and relating such gradients to mathematical instances of classical replicator dynamics [ 30 , 44 ]. Empirically, these gradients could include variation across human or mammalian hosts of different age or conditions, or environments under the effect of nutrients and abiotic factors such as moisture, light, temperature (e.g. microbial consortia in soil, cheese rinds, natural biofilms).

Naturally, a conceptual and empirical challenge arises when wishing to access invasion or colonization resistanceterm·key of a system. An idea can come from equation ( 2.4 ). Turning this equation around, one can see an indirect way of estimating Q = λ ¯ in a given system can bypass the precise knowledge of the N ( N − 1)/2 pairwise invasion fitnessterm·key coefficients of the constituent species, and just rely on a single observation of initial growth rate of any invader into the composite system (defined by z j ), coupled with empirical observation of N pairwise invasion experiments of such invader versus each of the species ( j = 1, ..., N ) in the system. Further averaging of such experimental procedure over a few invaders can augment statistical power for estimation of Q . The framework should flexibly accommodate protocols of invasion experiments that consider all factors governing Q , namely species number and composition, their frequencies, and their time-dependence.

With theoretical extensions ongoing, to increase model scope and generality, for example adding space [ 59 ], or more trait dimensions for species differences [ 31 , 39 ], including growth and clearance biases besides interactions, possibly relevant to more biological scenarios [ 60 ], we show the backbone of the replicator equationterm·mth formalism stays the same. This constitutes a big advantage as results on replicator dynamics [ 8 ] can be readily harnessed. The precise notion of invasion resistanceterm·dom, stemming from this model as mean invasion fitnessterm·key of a community, can enable deeper study of the role of species interactions in the ecology of antibioticterm·dom resistance [ 61 ], of general community-level cohesion independently of cooperation or competition [ 55 ], and provide more precise quantifiable links with data, from both in vitro and in vivo systems [ 21 , 52 , 62 ].

In the light of increasing calls for models to better describe microbial consortia and predict system-level properties [ 53 , 54 , 63 ], here we outline a new path for explicit mathematical investigation of colonization resistanceterm·key and mean invasion fitnessterm·key in a microbial ecosystem. This replicator equationterm·mth framework may help synthesize a new understanding of the ecology, development and evolution of host health and resilience to disease, mediated by microbiotaterm·key. Its extended applications could also provide broader insights on invasibility [ 11 , 45 , 46 , 64 ] and invasion resistanceterm·dom in other multispeciesterm·key ecosystems.

Acknowledgements We thank the several anonymous reviewers who provided challenging and useful feedback on earlier versions of the manuscript.

Appendix A A.1. N colonizing species varying along five traits also follow replicator dynamics The more general colonization model for variation in other trait dimensions between similar species can be found in [ 31 ] where again the replicator equationterm·mth with invasion fitnessesterm·key is derived mechanistically for the species frequency dynamics, using slow–fast dynamics arguments. Now, the mechanistic mutual invasion fitnessterm·key λ i j depends on all other traits, besides the matrix K = ( K ij ), describing micro-scale susceptibilities to co-colonization [ 25 ]. The five trait axes studied in [ 31 , 39 ] are: (i) species-specific transmission/proliferation rates β i , (ii) species-specific clearance rates in single colonization γ i , (iii) pair-specific clearance rates in co-colonization ( γ ij ), (iv) species onward propagation biases from mixed co-colonization, depending on their order of arrival ( p i j i ≠ p i j j ≠ 1 / 2 ), also known as ‘priority effects’, and (v) species pairwise susceptibilities to co-colonization K ij .

In this more complex scenario, the definition of ‘global context’ changes, as it may comprise also mean-field parameters in other traits, besides R 0 , k (relevant for the K ij -only model [ 25 ]). The N = 2-species system, in the case of 5-dimensional trait variation, is studied in [ 39 ], where the effect of the single to co-colonization ratio μ is shown to be very important, but this time influencing λ i j in a nonlinear fashion. The dynamic invasion resistanceterm·dom of the N -species community emerges in this model again as a mean invasion fitnessterm·key over all pairs in the system Q = λ ¯ , dynamically shaped by their changing frequencies. This model generalization shows that studying system dynamics and evolution in terms of the mutual invasion coefficients between constituent species is robust and widely applicable.

Ethics This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility Data and codes related to this study can be found on the Dryad Digital Repository: https://doi.org/10.5061/dryad.0k6djhb60 [ 65 ].

Supplementary files are provided in electronic supplementary material [ 66 ].

Declaration of AI use We have not used AI-assisted technologies in creating this article.

Authors' contributions E.G.: conceptualization, formal analysis, investigation, methodology, project administration, software, visualization, writing—original draft, writing—review and editing; S.M.: data curation, formal analysis, funding acquisition, investigation, methodology, software, visualization, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration We declare we have no competing interests.

Funding This study benefited from PESSOA grant no. 5666/44637YE awarded by Portuguese Foundation for Science and Technology (FCT) (to E.G. and S.M.) and cooperative agreement 2020-2001-230-Y20V7 by Le Studium, Loire Valley Institute of Advanced Studies (E.G.). E.G. acknowledges support also by CEECIND/03051/2018 and funding via FCT project 2022.03060.PTDC.

References Associated Data Data Citations Data Availability Statement Data and codes related to this study can be found on the Dryad Digital Repository: https://doi.org/10.5061/dryad.0k6djhb60 [ 65 ].

Supplementary files are provided in electronic supplementary material [ 66 ].