Erida Gjini
email: egjini@igc.gulbenkian.pt

Sten Madec
email: sten.madec@univ-tours.fr

The key to complexity in interacting systems with multiple strains

Keywords: SIS co-colonization model, quasi-neutrality, altered susceptibilites, single-to-coinfection ratio, cooperation-competition, stress gradient hypothesis, social interactions

Ecological community structure, persistence and stability are shaped by multiple forces, acting on multiple scales. These include patterns of resource use and limitation, spatial heterogeneities, drift and migration. Pathogen strains co-circulating in a host population are a special type of an ecological community. They compete for colonization of susceptible hosts, and sometimes interact via altered susceptibilities to co-colonization. Diversity in such pairwise interaction traits enables the multiple 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 systems remain poorly understood. In a recent study, we presented an analytic framework for N -type coexistence in an SIS epidemiological system with co-colonization interactions. The multi-strain complexity was reduced from O(N 2 ) dimensions of population structure to only N equations for strain frequency evolution on a long timescale. Here, we examine the key drivers of coexistence regimes in such a system. We find the ratio of single to co-colonization µ critically determines the type of equilibrium for multi-strain dynamics. This key quantity in the model 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, can only be achieved from a balance between higher competition in favourable environments, and higher cooperation in harsher environments, consistent with the stress gradient hypothesis in ecology. Multi-strain coexistence regimes are more stable when µ is small, whereas as µ increases, dynamics tends to increase in complexity. There is an intermediate ratio that maximizes the existence and stability of a unique coexistence equilibrium between strains. This framework provides a foundation for linking invariant principles in collective coexistence across biological systems, and for understanding critical shifts in community dynamics, driven by simple and random pairwise interactions but 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 [START_REF] Lotka | Elements of physical biology[END_REF] and Volterra [START_REF] Volterra | Fluctuations in the abundance of a species considered mathematically[END_REF] on mathematical models of the population dynamics of interacting species. This model has been later extended and sophisticated by many other theoretical studies [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Pascual | Ecological networks: linking structure to dynamics in food webs[END_REF], and is currently extensively used to characterize interaction networks in empirical microbiome communities [START_REF] Bucci | Mdsine: Microbial dynamical systems inference engine for microbiome time-series analyses[END_REF][START_REF] Stein | Ecological modeling from time-series inference: Insight into dynamics and stability of intestinal microbiota[END_REF]. 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 [START_REF] Serván | Coexistence of many species in random ecosystems[END_REF]. Modeling efforts seek to understand organizing principles for species composition, such as the balance between competition and cooperation [START_REF] Mougi | Diversity of interaction types and ecological community stability[END_REF], which many are increasingly arguing that should be studied together in integrated and context-dependent frameworks [START_REF] Bascompte | Mutualism and biodiversity[END_REF][START_REF] Coyte | Understanding competition and cooperation within the mammalian gut microbiome[END_REF]. Overall, analysis of such models with arbitrarily high dimensionality has been and continues to remain difficult.

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 [START_REF] Cobey | Niche and neutral effects of acquired immunity permit coexistence of pneumococcal serotypes[END_REF][START_REF] Gupta | Population structure of pathogens: the role of immune selection[END_REF][START_REF] Lipsitch | No coexistence for free: neutral null models for multistrain pathogens[END_REF][START_REF] Wearing | Ecological and immunological determinants of dengue epidemics[END_REF]. It is well recognised that population patterns of infection are to a large extent determined by susceptibility to infection. Most multi-strain SIR compartmental models, inspired for example from influenza, dengue, and malaria parasites, have focused on cross-immunity thresholds and thus competition between antigenically similar strains as the main driver of strain interactions and pathogen population structure [START_REF] Gupta | Chaos, persistence, and evolution of strain structure in antigenically diverse infectious agents[END_REF]. To deal with the inherent mathematical complexity, such studies have adopted model reduction, traveling wave solutions [START_REF] Gog | Dynamics and selection of many-strain pathogens[END_REF], and aggregation of variables [START_REF] Gomes | On the determinants of population structure in antigenically diverse pathogens[END_REF][START_REF] Lin | Dynamics of influenza a drift: the linear three-strain model[END_REF], which have led to key insights and analytical results. Studies have shown that even in a simple transmission model with a single infectious agent, where immunity by infection and vaccination reduce susceptibility to subsequent infections, altered susceptibility induces a 'reinfection threshold', above which reinfection dominates transmission, leading to higher infection levels and vaccination failure [START_REF] Gomes | The reinfection threshold promotes variability in tuberculosis epidemiology and vaccine efficacy[END_REF]. It is unsurprising that if multiple strains are involved, very complex behavior can emerge in infectious disease systems when such nonlinearities multiply and interact. Yet, it remains unclear which other factors, besides persistent host immunity, make strains compete or cooperate with each other, and whether there are environmental factors that can shift their fitness balance epidemiologically.

Here, we set out to study another mode of strain interactions, that has received, albeit less, modeling attention, namely interactions between strains via altered susceptibilities, not after, but during infection, in endemic scenarios of no immunity (SIS). Co-colonization vulnerability in SIS models with up to 2 strains, has been shown to be important in several studies [START_REF] Alizon | Multiple infections and the evolution of virulence[END_REF][START_REF] Davies | Within-host dynamics shape antibiotic resistance in commensal bacteria[END_REF][START_REF] Gaivão | Transmission fitness in co-colonization and the persistence of bacterial pathogens[END_REF][START_REF] Gjini | How direct competition shapes coexistence and vaccine effects in multi-strain pathogen systems[END_REF][START_REF] Lipsitch | Vaccination against colonizing bacteria with multiple serotypes[END_REF], but very few analytical investigations have been done for a larger number of interacting strains [START_REF] Adler | The dynamics of simultaneous infections with altered susceptibilities[END_REF]. In a recent mathematical study, extending the cocolonization interaction space to N -strains, we have provided a new analytical front for understanding a complex SIS system of multiple strains interacting upon co-colonization [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF]. Using timescale separation, we obtained a model reduction, where we showed how from the matrix of multi-type interactions arising in SIS cocolonization, we can obtain a special replicator equation [START_REF] Cressman | The replicator equation and other game dynamics[END_REF][START_REF] Hofbauer | Evolutionary game dynamics[END_REF] for strain frequencies. This reduced model formulation makes the entire dynamics more accessible to analysis, and links directly epidemiological SIS dynamics to Lotka-Volterra systems [START_REF] Bomze | Lotka-volterra equation and replicator dynamics: new issues in classification[END_REF] and evolutionary dynamics [START_REF] Nowak | Evolutionary dynamics of biological games[END_REF].

In the present article, we harness the simplicity of this co-colonization model framework to investigate coexistence and stability of such multi-strain systems, without immunity, but with variable co-colonization susceptibility coefficients among strains. We start by studying the behavior of the system for different global environmental 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 diversitystability configuration. We uncover rich transient and asymptotic behaviour of such systems, where steady states, limit cycles, multi-stability and chaotic attractors are possible.

We argue that the ratio of single-to co-colonization is a critical factor in collective dynamics, by modulating the asymmetry in pairwise mutual invasibility outcomes between types, and consequently, the dynamic complexity of the system as a whole. We argue that the analytically explicit form of this ratio in our model enables direct connection with the stress gradient hypothesis (SGH) in ecology [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Callaway | Competition and facilitation: a synthetic approach to interactions in plant communities[END_REF]. 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 such multi-strain systems and point to key global modulators of collective polymorphic coexistence in nature.

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 [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF]. 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 Materials and Methods), we describe the proportion of hosts in several compartments: susceptibles, S, hosts colonized by one type I i , and co-colonized hosts I ij , with two types of each combination, independent of the order of their acquisition. Fitness differences between types are encoded in how strains interact with each other upon co-colonization (K ij ), -whether there is facilitation between resident and co-colonizer (K ij > 1) or inhibition (K ij < 1), -assuming equivalence in transmission β and clearance rate γ. The structure of the K ij is central to multi-strain dynamics. The co-colonization interaction matrix between N strains in the system can arbitrarily include cooperative and competitive entries, but subject to small deviation from neutrality K ij = k + εα ij (ε small). 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 strains in the host population on temporal trends in individual strain frequencies.

In an earlier mathematical investigation [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF], we have derived in detail two timescales in this system: a fast one, given by the neutral model, and a slow one governed by the variation in co-colonization coefficients K ij . During fast dynamics the system stabilizes the total prevalence of single and co-colonization at the endemic persistence equilibrium, whereas over the slow timescale, under conserved global quantities, strain selective dynamics unfold. In the N equations, describing frequency evolution, the changing mean fitness of the collective appears explicitly, highlighting a key environmental feedback on all strains (see Box 1).

Box 1. Key model features for N -strain cocolonization and dynamics

• Deviation from neutrality in interaction trait space

The basis of the framework [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF] is to conceptualize each interaction coefficient between closelyrelated strains in co-colonization, as a mean value k plus some deviation from neutrality, thus re-writing it as:

Kij = k + εαij (1)
where 0 < ε < 1 is small. Thus, A = (αij) is the normalized interaction matrix, relative to the reference k. The parameter k is an indicator of mean interaction in co-colonization between any two strains, which when above 1 describes facilitation, and when below 1 describes competition. Another key epidemiological parameter in the model is the basic reproductive number, R0, describing the transmissibility of infection, defined as the average number of secondary infections that arise from single case in a naive population [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations[END_REF]).

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

On the fast time-scale (ε = 0), the system tends towards neutrality and reaches the manifold where total prevalence of susceptibles S, total prevalence of single colonization I and total prevalence of co-colonization D are conserved, and depend only on mean epidemiological parameters:

S = 1 R0 , I = N i=1 I * i = R0 -1 R0 + R0k(R0 -1) , D = 1≤i,j≤N I * i,j = 1 (R0 -1)k I, (2)
but, where each individual strain frequency is neutrally stable and free to vary zi ∈ [0, 1]. On a second time scale, which is the slow time scale τ = εt, consistent with weak selection, exact asymmetries between strains play out, and while the conservation law above still holds, individual strain frequencies zi are not free anymore. They obey deterministic slow dynamics, given the N -dimensional system:

   d dτ zi = Θzi •   j =i λ j i zj -Q(z)   , i = 1, • • • , N (3)
where Q(z) = 1≤k =j≤N λ k j zjz k , is a quadratic term symmetric on all strains, encapsulating the effect of the system as a whole on each individual strain. We have i zi = 1 and the rate Θ is given by: Θ = β (1-S)ID 2(1-S) 2 -ID .

• From dyads to collective dynamics

In the above equations, which drastically reduce the system from N (N -1)/2 + N to only N dimensions, the slow dynamics of strain frequencies are a direct function of pairwise invasion fitnesses λ j i between strains, which in our model, for each dyad (i, j), have the simple form:

λ j i = µ(αji -αij) + αji -αjj. ( 4
)
The quantities λ j i denote the initial rate of growth of strain i in an exclusion equilibrium where only strain j is resident, and depend on the ratio between single and co-colonization µ which is given by

µ = I D = 1 (R0 -1)k . ( 5
)
• Invariant principles in non-equilibrium multi-strain dynamics On the slow timescale, at all times, and for all strains, the following relationships hold:

zi = Ii I = Di D , (6)
stating that zi, as a measure of dominance of strain i in the total prevalence of carriage, occupies an equal relative frequency in single (Ii) and co-colonization (Di = j (Iij + Iji)/2)). After solving for strain frequencies zi, the Eqs. 3 can inform the epidemiological variables as follows:

Ii(t) = Izi(τ ); Iij(t) := Dzi(τ )zj(τ ). ( 7
)
The influence of mean-field parameters on tempo and type of coexistence

According to this model, global mean-field parameters, such as those affecting transmission rate, β, basic reproduction number R 0 , and mean interaction coefficient k affect explicitly the multi-strain dynamics over long time scales (εt). The equation for strain frequencies z i contains information for how they set the speed (Θ) and mode of strain stabilization (µ). More specifically, the relative dominance of single colonization in the system µ, because it appears as a factor in pairwise invasion fitness λ j i , modulates the importance of cross-strain asymmetry in co-colonization trait comparison between two strains: the higher µ is, the higher the dominance of single colonization, the more important the asymmetry between i and j becomes (a ji -a ij ). While this first component of invasion fitness is sensitive to feedbacks from global transmission and mean parameters between strains, the other component of the λ j i is independent of such global quantities and depends only on characteristics of the resident, namely on the relative strength of inter-vs. intra-strain interaction of the resident (a ji -a jj ).

Notice, that for fixed normalized variation among strains, hence rescaled matrix A = (a ij ), increasing µ in the system, amplifies asymmetries and should lead to more competitive exclusion between strains. We illustrate this for the case of N = 2 in Fig. 1a, where for randomly generated rescaled interaction matrices A (Normal distribution with mean 0 and variance ε), increasing µ stretches the region occupied by λ j i towards the competitive exclusion zone (λ 2 1 > 0, and λ 1 2 < 0). This phenomenon escalates also to the case of higher N , where effectively the multi-strain network is derived from pairwise invasion fitnesses between any two constituent strains (Fig. 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 subject to more constraints.

Although the mean number of strains that can coexist (n) apparently increases with µ (See Supplements for N = 10), the possibilities for stable coexistence tend to be more and more restricted to only special numbers of strains (odd values), as expected generically in special cases of the replicator equation. To illustrate further the role of µ, we show the dynamics of strain frequencies as a function of µ in Fig. 2, for the case of N = 3 and N = 4. Under fixed asymmetries A, changing the ratio between single and co-colonization, shifts the multi-strain dynamics from stable coexistence towards limit cycles and ultimately towards May-Leonard type oscillations, with effectively only one strain persisting for long periods of time, to be subsequently replaced by another one, and so on. These extreme regimes of behavior are expected in the replicator equation for special cases, but here the novelty is that we have identified a global environmental quantity, the ratio between single and co-colonization, as a tuning parameter for moving our epidemiological system between such extremes (see Box 2).

Thanks to the explicit formula for this critical ratio, we can see that µ = I/D = 1/k(R 0 -1) can change in two ways: i) either by changing basic reproduction number R 0 , or by changing k, the mean interaction coefficient between strains in co-colonization. If R 0 increases µ decreases, and similarly, if k increases (strains tending towards facilitation) µ also decreases. From the expression for λ j i , and the illustrations above, µ determines completely the type of equilibrium, for a given set of standardized variation in strain interactions. However if µ changes, either via shifts in R 0 or in k, the speed of the actual selective dynamics between strains may also change. In the case of R 0 effects on µ, this depends on how R 0 changes, via transmission rate or duration of infection changes. Thus we have 3 cases, whereby changes in µ can lead to changes in Θ (See Supplementary Figure 2). When changes in µ are mediated by changes in k, for high values of µ, the speed of selection goes down. When changes in µ are mediated by R 0 , for µ small (i.e. high R 0 ) there may be a small region where increasing µ (either via decreasing β or increasing m) can actually speed up the dynamics. This then would imply, that for relatively high R 0 , reduction in R 0 (via faster population turnover, shorter duration of carriage, or lower transmission rate) not only acts to increase the complexity of the dynamics, but also speed up of the multi-strain selection, thus contributing to quicker settlement of the competitive hierarchies between entities at an equilibrium.

Figure 1: Collective coexistence from pairwise invasion fitnesses and the effect of single-to cocolonization ratio µ. a. Effects of µ on the repartition of the pairwise interactions. In our model, there are four possible outcomes between two strains (depending on the values of λ j i ): extinction of 1, extinction of 2, coexistence or bistability. The coefficients α ij of the matrix A = (α ij ) are generated following a normal distribution N (0, 1). When µ → 0, the partitioning between the four outcomes is the same. When µ increases, the strain pairs, for the same A matrix, follow more likely competitive exclusion dynamics. b. Effect of µ on the multi-strain pairwise λ networks. Red line: coexistence, blue line: bistability, and a grey arrow is a competitive exclusion (the arrows shows the winner between the two). Here we fixed the normalized interaction coefficients between strains (matrix A) and N = 10 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, with the grey edges becoming more common. This feature makes stable coexistence harder for µ large.

Box 2. Single to co-colonization ratio as a tuning parameter for complexity One limit for µ may be achieved by either increasing cooperation k → ∞, or increasing transmission intensity R 0 = β m → ∞. The other limit can be achieved by increasing competition 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 Θ lim in these limits, see also Supplementary Text S2, and Figures S2-S3. A specific example with dynamics close to these limits is provided in Figure S4.

• Limit 1: Co-colonization, µ → 0 (R0 ↑ or k ↑)
Quality of the dynamics Recall that λ j i > 0 implies that the strain i (the invader) may invade the system with only strain j (the resident) present, starting from a very small frequency. In the limit µ → 0, this pairwise invasion fitness reduces to:

λ j i = αj,i -αj,j.
Importantly, in this limit, the fitness λ j i depends only on the coefficients of the resident, where αj,i measures how much the resident j contributes to mixed cocolonization with strains i and j: Dij (transmission and fitness of other strains), and αj,j how much the resident j contributes to self co-colonization: Djj (transmission and fitness of itself). Hence, the relative fitness of the invader i depends only on the resident j. Invasion is possible if and only if αj,i > αj,j.

This phenomenon in a system of N players can increase the availability of niches. For example, if all intra-specific coefficients are higher than inter-specific coefficients: αj,j > αi,j for all j = 1, • • • N and i = j then all strains may be stable residents when alone, and we would have at least N 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 αjj = 0 for all j, then Λ = A T , and the dynamics could be richer, in line with generic replicator systems [START_REF] Yoshino | Rank abundance relations in evolutionary dynamics of random replicators[END_REF]. In general, for random A, with no special structure, qualitatively, the system likely has many stable steady states wherein only few species coexist.

Speed of the dynamics

i) Θ lim = 0. If k → ∞, we have Θ ∼ m 2k → 0.
For fixed population turnover (m: clearance and death rate), the selective dynamics happen very slowly, i.e. the dynamics are effectively neutral. Similarly if R0 → ∞ via decreases in m (m → 0), the multi-strain system stays very close to neutrality.

ii

) Θ lim > 0. If R0 → ∞, via increases in trans- mission rate, β → ∞, then Θ → m
2k > 0 which is far from zero. In this case, the model reduction captures well the real dynamics and the qualitative study of µ → 0 is appropriate.

• Limit 2: Single colonization, µ → ∞ (R0 ↓ or k ↓)
Quality of the dynamics Recall the matrix Λ is the matrix of all pairwise invasion fitnesses between N strains in the system. In this limit, we have generically Λ → +∞ and Θ → 0. In order to keep a bounded matrix Λ we rewrite the system as

żi = µΘzi µ -1 Λ z i -z t µ -1 Λ z , (8)
where i = 1, • • • , N . Then the speed is given by µΘ and the qualitative behaviour by µ -1 Λ. Notice that when µ → +∞, the matrix

µ -1 Λ → A T -A.
Importantly, this matrix is skew symmetric, for which there are known mathematical results of the replicator equation in zero-sum games (see [START_REF] Allesina | A competitive network theory of species diversity[END_REF][START_REF] Chawanya | Large-dimensional replicator equations with antisymmetric random interactions[END_REF][START_REF] Fisher | Optimal strategies for random tournament games[END_REF]. The quadratic term in Eq. 8 is zero:

Q(z) = z T Λz = 0.
So the N -strain frequency dynamics reduce to:

żi = Θ lim zi (A T -A)z)i . (9)
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 n the number of strains coexisting at equilibrium, for this skew-symmetric case of the replicator equation, it has been shown that only an odd number of strains may coexist (see also Fig. S1,S5). More specifically, the probability to observe n = k strains, out of a total pool of N is

P (n = k) = 0 if k is even N k 2 1-N if k is odd.
Speed of the dynamics i) Θ lim = 0. If R0 → 1 decreases transmission intensity, either via lower β or higher m, then Θ → 0, and selective dynamics are too slow, making the system effectively neutral.

ii)

Θ lim > 0 If k → 0, µΘ → m 2 (R0 -1) > 0.
Thus the selective dynamics are well defined and occur on a feasible timescale.

Coexistence, stability and diversity

The community dynamics resulting 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. While there is an ongoing debate in ecology about the relationship between stability and diversity in ecosystems [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Mccann | The diversity-stability debate[END_REF][START_REF] Odum | Fundamentals of ecology[END_REF][START_REF] Tilman | Biodiversity and stability in grasslands[END_REF], here we set out to examine this Here we fixed the normalized interaction coefficients between strains (matrix A) and we varied µ. Because µ affects the λ j i , for each µ we have a different effective competition network between strains modulated by top-down factors. The plots show the strain frequency (z i ) dynamics following the equations on the slow manifold. For N = 3 and N = 4, for small values of µ, the final attractor is a stable steady state, while for a large value of µ, the stable attractor is a union of heteroclines of 3 strains (May Leonard type). For intermediate values of µ (µ = 0.5) we obtain a stable spiral of coexistence for N = 3 and a limit cycle of coexistence for N = 4, both consistent with an increase in rate of strain turnover. relationship in our system, assuming random interaction structure between strains, namely the matrix A and fixing the total pool size of strains N . Thus, we study the model's rich behaviour for a larger pool of strains in the system, namely N = 10 (Figure 3). We sample the normalized interaction matrix A randomly from a normal distribution. For each interaction matrix, we compute numerically the feasible equilibria of the system (z i ≥ 0, ∀i), and for each equilibrium found, we evaluate its local stability and its associated Shannon entropy

H = - N i z * i log(z * i )
. By simulating different randomly generated structures of the matrix, 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 [START_REF] Mccann | The diversity-stability debate[END_REF], although there are also other notions of general stability related to permanence [START_REF] Law | Permanence and the assembly of ecological communities[END_REF] that could be informative.

The first aspect of coexistence we focus on is the probability of reaching a stable n-strain coexistence, starting with a pool of N = 10 strains. Using our model, under random allocation of co-colonization interaction coefficients, we find that the probability of stable n-strain coexistence decreases with n, thus confirming that the more strains there are, with random mutual interactions in co-colonization, the harder it is for all of them to coexist (colored dots in Fig. 3a). However, notice that in contrast to the stable steady states, there are many more unstable steady states in the system (98%), spread at any diversity and stability level, with no apparent order or hierarchy (black dots in Fig. 3a). This suggests that unstable coexistence is much more likely and that the number of coexisting strains in that case is not a limiting factor. For a fixed number of strains coexisting n, out of a total pool of N , the higher the diversity (Shannon entropy) of the equilibrium, the higher the stability, as measured by the negative part of the dominant eigenvalue of the Jacobian matrix evaluated at that equilibrium.

Then, to quantify in more detail this relationship, accounting for the feasibility of the steady state itself, we go on to examine how the probability of stable coexistence for a given n depends on the diversity of the equilibrium. Using Shannon entropy as a measure of diversity of the coexistence equilibrium, we find that among stable equilibria, the proportion of systems at a given entropy level decreases with entropy, and this pattern is consistent for all n. This implies that under randomly sampled interactions, stable coexistence equilibria with equally-distributed strain frequencies, and thus more diversity, are harder to obtain.

Thus, on one hand we find more diversity is associated with more stability, in support of other ecological studies [START_REF] Mccann | The diversity-stability debate[END_REF], but also that higher-diversity equilibria are less feasible and harder to reach under random sampling of interaction coefficients from a pool of N strains. When combining these two forces: diverse Figure 3: Diversity-stability relationships for n-strain coexistence starting with a pool of N = 10 strains. a. Summary of 100.000 simulations with a random normalized interaction matrix A, with each entry α ij drawn from the Normal distribution N (0, 1). µ = 0.05. 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 n. c. Plot of the rank-order frequency distribution of strains at the "optimal diversity" for n = 4 strains coexisting out of a pool of 10. d. Plot of the rank-order frequency distribution of strains at the "optimal diversity" for n = 8 strains coexisting out of a pool of 10.

coexistence is rarer to obtain, but generally more stable, we find that this trade-off between probability of stable coexistence and stability, give rise to an optimal intermediate diversity for n-strain coexistence (Fig. 3b)), where stable coexistence is on one hand sufficiently feasible, and the resilience of that equilibrium to perturbations is sufficiently high. Thus, mean stability of n-strain coexistence is effectively maximized at particular, typically intermediate, values of Shannon entropy (Fig. 3c-d). This implies that in practice, independently of which n strains ultimately coexist in a system out of a pool of N , for a given n, we should tend to observe such system around the same rank-frequency distribution of strains. Such "optimal" rank-order frequency distribution will display higher variance for low number of coexisting strains (for example n = 4), but be more robust for higher n, (e.g. n = 8), and generally independent of community composition.

While Figure 3 refers to a scenario of low µ where co-colonization is dominant, in the Supplementary Figure S5 we explore also the case of a high µ tending to the opposite limit of single colonization dominance (µ = 10 is more realistic for a system like pneumococcus serotypes). Strikingly the fundamental feature that diversity increases stability, is preserved. Likewise, the trade-off between probability of reach and stability of coexistence equilibria is also preserved. This means the ecological forces leading to an optimal effective stability in an intermediate entropy range are robust, independently of µ. In this limit, we also verify that for the same random structures of interaction coefficients A, typically more strains coexist, an odd number of strains is more likely to coexist [START_REF] Chawanya | Large-dimensional replicator equations with antisymmetric random interactions[END_REF], and that their coexistence dynamics are more complex (Box 2).

The role of µ in multi-stability

We revisit the relationship between complexity of the dynamics and ratio of single to co-colonization (µ = 1/k(R 0 -1)) in our model, now with simulations over a larger range of µ. Using 100.000 random system simulations with N = 10 for different values of µ, we confirm the pattern observed for N = 2 and illustrated for N = 4 in Figure 2, that increasing µ (via decreasing either R 0 or k) 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 (µ → 0) as shown in Fig. 4a. 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 N and normalized interaction matrix A between strains.

Figure 4: System complexity for large N and multi-stability. a. For several values of µ, 100.000 systems of N = 10 strains have been generated (using a normal distribution of A). For each value of µ, we show the proportion of systems with a given number of stable steady states (in {0, • • • , 6}). As µ increases, the proportion of systems without any stable steady state increases, reflecting larger potential for complex coexistence dynamics, where the attractors may be limits cycles, May Leonard trajectories or even strange attractors. For µ < 0.5 (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 µ = 0, and for each of the 100.000 systems, we compare the strain overlapping index SOI (see Materials and Methods) 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 0. This shows that in multi-stability, the attractors contain mostly disjoint sets of strains.

In particular, when zooming-in on the multi-stable equilibria 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 Materials and Methods), we find that multi-stability is largely characterized by non-overlapping subsets of strains (Fig. 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 N 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 µ = 10.

When quantifying system outcomes for N strains, in three broad categories: monostability, multistability, and complexity, with an even higher resolution for µ, we find that there are effectively 3 regions (Figure 5): for small µ, stable coexistence and multi-stability is more likely, for intermediate µ monostability is more likely, and for high µ, (i.e. very small R 0 or k) complexity, May-Leonard type coexistence and unstable coexistence are more probable outcomes. This result highlights the importance of the environmental gradient in shifting the qualitative dynamics: the more cooperation there is between all strains on average (k higher), or the higher the overall transmission intensity (higher R 0 ), i.e. the lower µ is, the weaker the selection imposed by co-colonization asymmetries, and hence stable coexistence (in multiple combinations) is more probable. In contrast, when µ increases, either via more competition between strains on average (k lower), or lower transmission intensity (R 0 closer to 1), the stronger the importance of asymmetries, and thus the stronger the selection imposed by co-colonization trait differences between strains. In that extreme, because the strains are so inter-dependent, coupled more strongly via the available dynamic host resources in the single infection compartment (implying dynamic niche creation and destruction along at most N axes), the probability of unstable coexistence between them is largest. In one extreme, we expect low variability in a setting over time, in the other extreme we expect high variability and complex strain turnover over time.

Constructing multi-strain communities with random co-colonization interaction strengths (in normalized matrix A form), we find that the global epidemiological quantity µ = 1 (R0-1)k , namely the single to co-colonization ratio modulates regimes of system behavior and provides insights on the diversity-stability debate [START_REF] Mccann | The diversity-stability debate[END_REF]. For a fixed µ and a given set of interaction coefficients between strains however, as their frequencies Figure 5: Dynamic regimes as a function of singe to co-colonization ratio µ. We plot the number of steady states as a function of µ ∈ [0.01, 100] for many randomly-generated interaction matrices A. Since µ depends on R 0 and k, we kept their combination implicit, and randomly sampled the rescaled interactions matrix A = (α ij ). For any system (A), 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, in principle every system with a unique global attractor is monostable, but there may be other more complex unique attractors like stable limit cycles. 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 N = 10, and for each value of µ we generated randomly 1000 matrices A leading to 1000 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 (grey line). For a small µ, very few systems have zero stable steady states, while nearly 80% of them are multistable. For a large value of µ the multistability is nearly impossible while more than 70% of interaction systems having no steady state and thus displaying a complex dynamics with no stationary attractors. change epidemiologically, the epidemiological dynamics and the evolutionary dynamics interact. Evolution happens on two levels in this system. First, frequency dynamics (z i ) can be linked with selective dynamics on mean interaction in co-colonization trait space (see [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF] where k ef f ective = k + εq(t). Secondly, another trait changing among community members is mean invasibility of the system, Q, 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 the product between 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 particular, in the µ → 0 limit, where multi-stability is common, systems can evolve towards more competition or more cooperation in co-colonization depending on initial conditions. Conversely in the other limit µ → ∞, where complex oscillatory coexistence is more probable, systems preserve polymorphism in interaction trait space, with alternating periods between competitors and cooperators. 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 increase in mean resistance to invasion of the multi-strain system (Q evolutionary dynamics) does not necessarily always correlate with increased cooperation or competition (higher or lower q), and in the oscillatory regime the two can also be synchronous and asynchronous (Figure S6).

Box 3. Stability-diversity-complexity and evolutionary dynamics of interaction traits

In multi-strain communities with random co-colonization interaction strengths (A), the global epidemiological quantity µ = 1 (R 0 -1)k , describing the single to co-colonization ratio modulates regimes of system behavior, shaping its stability, complexity and diversity properties.

Together with epidemiological dynamics, evolution happens in real time on two traits in this system: i) mean interaction in co-colonization trait space (q) where k ef f ective = k + εq(t), ii) mean invasibility of the system Q in mutual invasion fitness space (Eq.3).

Evolutionary dynamics in the microbial community reflect the changing mean fitness landscape, and depend on frequency dynamics [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF] as follows:

q(z) = 1≤k,j≤N α jk z k zj and Q(z) = 1≤k =j≤N λ k j zjz k .
Notice that q increasing means the system tends to increase cooperation on average over time, whereas q decreasing means the multi-strain community tends to increase competition 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. • Mean trait evolution on 2 levels: q and Q over time

(R0 = 3)
• Either competition or cooperation can be selected • Evolution of "specialist" communities • High diversity between host populations

• Competition or low transmission (µ high)

-Many species coexist -Unstable coexistence common (no multi-stability)

-Complex attractors (limit or heteroclinic cycles)

-High variability over time in one population

• Mean trait evolution on 2 levels: q and Q over time

(R0 = 2)
• Competition and cooperation fluctuate! • Evolution of "generalist" communities • Lower diversity between host populations

The R 0 vs. k trade-off in µ and the stress gradient hypothesis

The structure of a multi-strain system is the result of complex interactions among different members, where the relative importance of each interaction can also be a function of the global environmental context, and of the biological interactions among the rest of the community. This context-dependence in effects of ecosystem engineers for the rest of the community has long been highlighted in ecological systems [START_REF] Hastings | Ecosystem engineering in space and time[END_REF].

In our model, we have seen that strains can mutually engage in more or less competitive interactions and more or less cooperative interactions, and differential abilities in such traits allow each strain to potentially act as an engineer. Thanks to this similarity assumption for dynamics decomposition, we can disentangle very clearly the effect of mean interaction k and the variation around that mean in what concerns intra-specific and inter-specific interactions, A = (α ij ), and we can also see how mean 'environmental' variables contribute to relative fitness between system members λ ij . Until now, we have seen that in our system, the ratio of single to co-colonization µ = 1/(R 0 -1)l) is a key modulator of the qualitative complexity of the dynamics (Figure 5, Box 2). Thanks to the simplicity of this formula, it is clear that both transmission intensity R 0 and mean interaction coefficient in cocolonization k decrease this ratio, and the only way for µ to be held constant, is if R 0 and k are traded-off against one another. In other words, in order to keep µ constant when transmission opportunities for all strains are reduced, the mean interaction coefficient in the system must increase, i.e. if all strains cooperate more in co-colonization. And viceversa, if the environment gets more favourable, i.e. R 0 increases, then µ can be kept constant only if k decreases, thus if all strains increase competition. Such exact quantitative finding in our model makes an explicit link with the stress-gradient hypothesis (SGH) postulated in ecology [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Callaway | Competition and facilitation: a synthetic approach to interactions in plant communities[END_REF]. The SGH predicts that positive interactions should be more prevalent in stressful environments, while more favourable environments should favour competition. This hypothesis has been supported and tested mainly in plant systems [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF][START_REF] Eränen | Increasing intraspecific facilitation in exposed environments: consistent results from mountain birch populations in two subarctic stress gradients[END_REF][START_REF] Pugnaire | Changes in plant interactions along a gradient of environmental stress[END_REF], and more recently in microbial communities [START_REF] Fetzer | The extent of functional redundancy changes as species? roles shift in different environments[END_REF][START_REF] Hoek | Resource availability modulates the cooperative and competitive nature of a microbial cross-feeding mutualism[END_REF][START_REF] Lawrence | Evolution of resource use along a gradient of stress leads to increased facilitation[END_REF][START_REF] Mccluney | Shifting species interactions in terrestrial dryland ecosystems under altered water availability and climate change[END_REF][START_REF] Piccardi | Toxicity drives facilitation between 4 bacterial species[END_REF]).

In our system, which is an epidemiological system at basis, but with microbial interactions embedded in co-colonization dynamics, we can verify via the single-to co-colonization ratio, an interpretation of the stressgradient hypothesis. In order for a set of species with a fixed set of normalized interactions (A=(α ij )) between them, to maintain a given configuration of coexistence (µ) under changing environmental conditions (e.g. R 0 ), the only way to do it is by universal adaptation in mean-field interaction, increasing competition in favourable environments and increasing cooperation in harsher environments, respectively (see Figure 6). However, it is worth noting that even if µ is kept constant, i.e. even if the species respond to harsh environmental changes by increasing cooperation, and thus preserve their qualitative equilibrium configuration, the speed of their collective dynamics may still be affected, typically reduced with R 0 , although the absolute magnitude will depend on whether and how much they are mediated through altered transmission rate vs. altered duration of carriage (see Supplementary Figures S3-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, and one way to keep that balance in check would be by quick adaptation of the mean.

Figure 6: Transmission intensity R 0 and mean interaction in cocolonization determine µ in a tradeoff manner. While the type of equilibrium where the system will tend to, is strongly driven by µ (the ratio of single to co-colonization), hence directly by R 0 and k, it is clear from the expression µ = 1 (R0-1)k that the only way for µ to be held constant, is if R 0 and k are traded-off against one another. This principle reflects the expectation from the stress-gradient hypothesis (SGH): when transmission opportunities for all strains are reduced (harsher environment), the mean interaction coefficient in the system must increase, i.e. if all strains should cooperate more in co-colonization, to keep the same µ and hence qualitative coexistence regime. And viceversa, if the environment gets more favourable, i.e. R 0 increases, then µ can be kept constant only if k decreases, thus if all strains increase competition. Note that keeping µ fixed only preserves the qualitative nature of the dynamics between strains, but does not guarantee that its speed Θ will remain constant. Changes in global R 0 and k and the underlying mechanisms giving rise to such changes, will be important drivers of changes in selection speed.

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 cocolonization, which broadly include both competition and facilitation. By clearly delineating the role of mean interaction coefficient, basic reproduction number R 0 , and biases in pairwise coefficients relative to the mean, we obtain a model reduction for strain frequency evolution for any N . The questions we addressed within such a system are similar to a long-standing ecological quest on multispecies dynamics [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Pascual | Ecological networks: linking structure to dynamics in food webs[END_REF][START_REF] Serván | Coexistence of many species in random ecosystems[END_REF]: 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 [START_REF] Dawid | The blp bacteriocins of Streptococcus pneumoniae mediate intraspecies competition both in vitro and in vivo[END_REF][START_REF] Leggett | War and peace: social interactions in infections[END_REF][START_REF] Lysenko | Within-host competition drives selection for the capsule virulence determinant of streptococcus pneumoniae[END_REF][START_REF] Riley | The ecological role of bacteriocins in bacterial competition[END_REF][START_REF] Shen | Pneumococcal quorum sensing drives an asymmetric owner-intruder competitive strategy during carriage via the competence regulon[END_REF] quantitatively we can predict very important features of system behavior with a simple mathematical framework derived under quasi-neutrality and strain similarity assumptions [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF]. This model captures selective dynamics between strains over long time, and coincides with an instance of the replicator equation [START_REF] Hofbauer | Evolutionary game dynamics[END_REF], parametrized in our case, in terms of pairwise invasion fitnesses between strains, a central quantity in adaptive dynamics [START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF]. Studying this equation and its epidemiological consequences in detail, here we find that there is a critical ratio that tunes asymmetry and dynamic regimes in such multi-strain interacting system: namely the ratio of single to co-colonization. This ratio µ, given by the inverse of the product of basic reproduction number and mean permissiveness to co-colonization µ = 1/((R 0 -1)k), increases or decreases the importance of asymmetry in interaction between strains as they encounter each-other epidemiologically in resident-co-colonizer combinations, and lends support to the context-dependence [START_REF] Coyte | Understanding competition and cooperation within the mammalian gut microbiome[END_REF] of relative fitnesses in a coupled microbial community.

In our epidemiological system, individual strain frequency evolution over long times depends on how that strain interacts with other strains, but also on the total net interaction between extant strains. Such global, symmetric and temporally-varying, environmental feedback on all strains emerges as a key feature of this cocolonization model, and can be used to summarize the evolution of diversity in the system. It 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 and chaos. Each strain is competing for susceptible hosts, and for (at most) N singly-colonized hosts that can be co-colonized. These N types of available hosts are the resources upon which the strains can grow epidemiologically, but these resources are not independent from each other. With each additional strain, the landscape of interactions can change dramatically for any system composition, because each additional strain offers a new type of resource upon which the competition of the residents can be manifested, while also modulating directly the relative success of those strains by taking over space from their colonization. This dynamism embedded in the inter-connected and inter-dependent system is the source of dynamic complexity, recognized generally in biological games [START_REF] Nowak | Evolutionary dynamics of biological games[END_REF], and here uncovered and quantified in an explicit epidemiological context.

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 non-overlapping strain subsets can coexist, depending on initial conditions. This result is consistent with earlier multi-strain SIRS models with cross-immunity, that in order to minimize competition and maximize chances of persistence, polymorphic pathogens self-organize in non-overlapping subsets [START_REF] Gupta | Population structure of pathogens: the role of immune selection[END_REF]. Here, within a much simpler, SIS framework, we find a similar principle applies, despite the interactions between strains being mediated via altered susceptibilities to co-colonization, without persistent immune memory.

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, hydrology and economics are calling attention on critical transitions, which occur when natural systems drastically shift from one state to another [START_REF] Scheffer | Anticipating critical transitions[END_REF]. 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 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.

In a recent theoretical study focused on malaria, without explicitly modeling multiple strains and host immune memory, superinfection (concurrent multiple infection) has been shown to consistently create tipping points that give rise to hysteresis in responses to control or seasonal variation in vector abundance [START_REF] Alonso | Critical transitions in malaria transmission models are consistently generated by superinfection[END_REF]. Our present work, developed in the framework of many interacting strains, lends support to a similar perspective, but more generally relevant to high-dimensional conservative SIS systems and coexistence regimes, rather than prevalence tipping points. Cooperation and competition among strains, studied here under the same integrated context, 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 (Fig. 5), but coexistence of more strains becomes possible (Fig. S1). Although at first sight this may seem to suggest, somewhat contrary to previous expectation from lower-dimensional models [START_REF] Chen | Fundamental properties of cooperative contagion proc[END_REF][START_REF] Hébert-Dufresne | Complex dynamics of synergistic coinfections on realistically clustered networks[END_REF], that average cooperation in co-colonization is stabilizing and average competition is destabilizing, our result should be read in the broader view of the underlying resource dynamics shaped dynamically by N strains, and the opportunities for small trait differences between them to get manifested in a homogeneously mixing host population. When contact structure heterogeneities are accounted for, it would be interesting to test, how much of this behavior will still be preserved.

Taken together, our results show that even though differences in pairwise interactions can be small and random, as expected for closely related strains, a large community of such strains can collectively coexist, and self-organize, albeit often in complex and unstable fashion, simply thanks to the critical interplay between global transmission potential and mean permissiveness to co-colonization. We found that the requirement of being simultaneously stable and feasible tends to push coexistence regimes toward intermediate diversity, independently of precise composition. We argue that the stress gradient hypothesis (SGH) [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Callaway | Competition and facilitation: a synthetic approach to interactions in plant communities[END_REF] can be used as a framework to help interpret the critical role of the single-to colonization ratio in shaping epidemiological multi-strain systems.

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 cocolonization 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-strain pathogen systems. For example, empirical data in polymorphic S. pneumoniae bacteria, have been consistent with estimates of about 90% mutual inhibition between co-colonizing serotypes [START_REF] Gjini | How direct competition shapes coexistence and vaccine effects in multi-strain pathogen systems[END_REF][START_REF] Lipsitch | Estimating rates of carriage acquisition and clearance and competitive ability for pneumococcal serotypes in Kenya with a Markov transition model[END_REF], and R 0 values around 2. This implies that for this system, µ ≈ 10, and multi-stability is highly unlikely but a single stable equilibrium point is almost as likely as complex unstable coexistence, which may reconcile the secular trends observed in some settings consistently over many years [START_REF] Ekdahl | Bacteraemic pneumococcal infections in southern sweden 1981-96: trends in incidence, mortality, age-distribution, serogroups and penicillin-resistance[END_REF][START_REF] Feikin | Historical changes in pneumococcal serogroup distribution: implications for the era of pneumococcal conjugate vaccines[END_REF][START_REF] Fenoll | Evolution of streptococcus pneumoniaeserotypes and antibiotic resistance in spain: Update (1990 to 1996)[END_REF]. Such secular trends can interfere with vaccine introductions and need to be accounted for when estimating impact [START_REF] Moore | Rethinking replacement and resistance[END_REF][START_REF] Moore | Emergence of nonvaccine serotypes following introduction of pneumococcal conjugate vaccine: cause and effect?[END_REF]. Our model makes explicit predictions about the timescale and qualitative aspects of such secular trends (Box 2, Fig. S4), under the plausible assumption that they are driven by co-colonization interactions.

The consistency of the optimal diversity (entropy) configuration for a given number of strains, is a model prediction that has in fact been observed empirically for pneumococcus serotypes across geographical settings [START_REF] Hanage | Evidence that pneumococcal serotype replacement in massachusetts following conjugate vaccination is now complete[END_REF]. While the serotype composition has changed in different populations before and after vaccination, the rank-frequency distribution of serotypes, related to entropy in our model, has been conserved. Our model, offers a new perspective on such empirical observations, in terms of a null expectation for how such pattern emerges and is maintained, without needing to invoke strain-specific transmission rate dominance or duration of carriage dominance [START_REF] Weinberger | Pneumococcal capsular polysaccharide structure predicts serotype prevalence[END_REF], but simply based on co-colonization interactions between serotypes, which when random balance out in a particular optimal entropy configuration. Since this optimal rank-order abundance distribution depends on the context µ, model expectations for such dependence can be tested with data from different R 0 settings. Finally, the link with the stress gradient hypothesis, suggests that to preserve a certain "optimal" single-to cocolonization ratio (optimal complexity/coexistence balance), independently of community size N , facilitation in co-colonization between microbial strains should be more common in epidemiological settings with low prevalence, and this is an interesting prediction to test 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 revisit and interpret multi-strain dynamics in other systems of relevance for public health, for example influenza [START_REF] Yang | Dynamic interactions of influenza viruses in hong kong during 1998-2018[END_REF]), dengue (Mier-y Teran-Romero et al., 2013), malaria [START_REF] Alonso | Critical transitions in malaria transmission models are consistently generated by superinfection[END_REF][START_REF] Gupta | Exploring the evolution of diversity in pathogen populations[END_REF] or human papilloma viruses [START_REF] Murall | Revising ecological assumptions about Human papillomavirus interactions and type replacement[END_REF].

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 vs. equalizing forces in biodiversity.

Methods

N-strain SIS model with co-colonization

The model is given by the following equations:

     Ṡ = m(1 -S) -S N j=1 F j , İi = F i S -mI i -I i j K ij F j , 1 ≤ i ≤ N İ ij = I i K ij F j -mI ij , 1 ≤ i, j ≤ N (10)
where

F i = β I i + N j=1 1 2 (I ij + I ji )
gives the force of infection of strain i, and S, I i and I ij refer to the proportion of susceptible hosts, singly-colonized by strain i and co-colonized by strains i and j. Notice that S = 1 -I i + I ij , thus the dimension of the system is effectively N + N (N -1)/2. A key assumption in the model is that hosts co-colonized with different strains i and j transmit either with equal probability. In the above notation m = γ + r, is a population turnover rate, encapsulating both clearance rate γ of colonization episodes and recruitment rate of susceptible hosts r (balanced by natural mortality rate r = d).

This model has been described in detail elsewhere [START_REF] Gjini | A slow-fast dynamic decomposition links neutral and non-neutral coexistence in interacting multi-strain pathogens[END_REF] for the case of N = 2. The model allows for each strain to interact differently with other strains upon co-colonization, altering the suceptibility of an already-colonized host to incoming strains. The magnitude and type of such interactions is described by the matrix K, where values K ij above 1 indicate facilitation between strains, and values of K ij 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 lies on the assumption that the variance of such interaction coefficients is low, thus deviation from neutrality of K ij with respect to a reference k is small. In our particular simulations here, we assume a normal distribution for such interaction coefficients. We can write every K ij as: K ij = k + εα ij under the assumption that ε is small. Although many different parametrizations are possible, we define k = 1≤i,j≤N Kij N 2 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 ij ) 1≤i,j≤N . This leads to the normalized interaction matrix A = (α ij ) = K ij -k ε to have the same distribution as K, but with mean 0 and variance 1. We use 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 N . We use the reduced system (of N equations) for dynamics on the slow timescale εt [START_REF] Madec | Co-colonization interactions drive explicit frequency-dependent dynamics among n types[END_REF] to obtain the steady-states of the system and analyze their properties and stability.

Strain overlap index (SOI) for two coexistence equilibria of the same system Consider a set S = {E 1 , E 2 , • • • , E p } with all steady states of a given system (interaction matrix), where E j ⊂ {1, • • • , N } and E j = ∅. In our applications, E j is the index of feasible strains (with z i > 0) in the steady state number j of a given system of N = 10 strains. We note P(S) the set of each subsets of S. As instance, if S = {E 1 , E 2 } then P(S) = {∅, {E 1 }, {E 2 }, {E 1 , E 2 }} . Recall that |P(S)| = 2 |S| = 2 p .

We denote the strain overlapping index by SOI and define it as a mean over the Jaccard indices calculated for each pair of coexistence equilbria: We assume m = 0.5 (clearance + natural mortality rate) so that R 0 = 2β and a maximum time for simulation T max = 500. For each row the qualitative dynamics remain the same as shown on the graph. However, the speed of this dynamics depends on β. a. On the first row µ = 0.1, we have a competitive exclusion scenario.

Only strain 5 persists. When β = 1 (and then k = 10 to compensate via the trade-off for µ constant) all the strains are still present in the system at T max . For β = 5 (then k = 100/9), dynamics are much quicker and selection of finally only one strain occurs very rapidly. The values of Θ are respectively: 0.043, and 0.39. b. In the second row µ = 10. For this value of µ, the same matrix A yields more complex dynamics in general, in line with the result that qualitatively the dynamics are determined by µ. However, when β = 2, the dynamics seems to approach a cycle with dominance of strain 5 (grey) and extinction of strains 1, (blue) 4 (green) and 6 (purple). For β = 10 (so k = 1/90), the dynamics over the same timescale show a very complex behavior. First strain 6 appears to come back, and 5 goes to extinction, then strain 1 comes back again and ultimately strain 5 dominates in coexistence with 2 and 3. The values of Θ respectively are: 0.043, and 0.39. Thus increasing transmission rate has increased the speed about ten fold. (which is more realistic for pneumococcus). a. Relation between Shannon entropy (x-axis) and stability (yaxis).for difference number of n coexisting strains within a pool of N = 10. Due to this large value of µ = 10, it is very rare to observe an even number of coexisting strains n (see figure S1). b. However, the shape of the effective stability continues to be similar to the case µ = 0.05 for both even and odd n. For n ≥ 3, there is an "optimal diversity" which leads to the optimal balance between feasibility and stability of the steady state. c. and d. To be consistent, we illustrate the rank order at the "optimal diversity" for the odd number of coexisting species n = 3 and n = 7. where the dynamics of q and Q can be correlated or anti-correlated. Evolution of system resilience to invasion (Q high) can be obtained with evolution of more cooperation or more competition in the system, depending on initial conditions. c-d. Two cases of R 0 -2, where since µ is higher, the dynamics are more complex, and so is evolution of mean traits in the system. Remarkably q and Q can be synchronous or asynchronous. This regime corresponds more likely to maintenance of polymorphism in interaction trait space within the same population.

Figure 2 :

 2 Figure 2: Increasing the ratio of single to co-colonization (µ) increases the complexity of multistrain dynamics. We illustrate two examples of fixed interaction matrix A and shifting µ. a. for N = 3 and b. for N = 4.Here we fixed the normalized interaction coefficients between strains (matrix A) and we varied µ. Because µ affects the λ j i , for each µ we have a different effective competition network between strains modulated by top-down factors. The plots show the strain frequency (z i ) dynamics following the equations on the slow manifold. For N = 3 and N = 4, for small values of µ, the final attractor is a stable steady state, while for a large value of µ, the stable attractor is a union of heteroclines of 3 strains (May Leonard type). For intermediate values of µ (µ = 0.5) we obtain a stable spiral of coexistence for N = 3 and a limit cycle of coexistence for N = 4, both consistent with an increase in rate of strain turnover.

R0

  ↑ or k ↑ R0 ↓ or k ↓ • Cooperation or high transmission (µ low) -Few species coexist -Multistability is common -Simple coexistence (fixed point) -Low variability over time in one population

Figure S3 :

 S3 Figure S2: Speed of dynamics depends on µ and other parameters. The ratio of single to co-colonization µ can vary in three different ways. When k is varied: β = 2, m = 1, R 0 = 2, µ decreases with k (blue line).When R 0 is varied: via β: k = 1, m = 1, µ decreases with R 0 (red line). When R 0 is varied via m: k = 1, β = 1, µ again decreases with R 0 (black line). The effects of the same change in µ on the parameter Θ are different for these three different cases.

Figure S4 :

 S4 FigureS4: Keeping µ constant preserves the same qualitative dynamics in principle but there are timescale effects. For two values of µ (rows a and b) we compute the dynamics and we show the multi-strain frequency trajectories for two different choice of transmission rate β (columns 2 and 3). µ and β being fixed, we have k = 1/µ(R 0 -1). For the normalized matrix of interactions A we assume the same matrix in each figure

Figure S5 :

 S5 Figure S5: Diversity-Stability relationship for µ = 10. This figure is similar to figure 3 but for µ = 10(which is more realistic for pneumococcus). a. Relation between Shannon entropy (x-axis) and stability (yaxis).for difference number of n coexisting strains within a pool of N = 10. Due to this large value of µ = 10, it is very rare to observe an even number of coexisting strains n (see figureS1). b. However, the shape of the effective stability continues to be similar to the case µ = 0.05 for both even and odd n. For n ≥ 3, there is an "optimal diversity" which leads to the optimal balance between feasibility and stability of the steady state. c. and d. To be consistent, we illustrate the rank order at the "optimal diversity" for the odd number of coexisting species n = 3 and n = 7.

Figure S6 :

 S6 Figure S6: Evolutionary dynamics of mean traits in the N -strain system. a.-b Two cases of R 0 = 3,where the dynamics of q and Q can be correlated or anti-correlated. Evolution of system resilience to invasion (Q high) can be obtained with evolution of more cooperation or more competition in the system, depending on initial conditions. c-d. Two cases of R 0 -2, where since µ is higher, the dynamics are more complex, and so is evolution of mean traits in the system. Remarkably q and Q can be synchronous or asynchronous. This regime corresponds more likely to maintenance of polymorphism in interaction trait space within the same population.

The quadratic part is the scalarQ(z) = z T Λz. From Λ = -Λ we obtain Q(z) = Q(z) T = z T Λ T z = -Q(z), thus Q(z) = 0.

Acknowledgement

EG acknowledges support by the University of Tours and EURAXESS for an invited researcher stay at IDP, in Tours during 2019, which was instrumental to complete this collaborative work.

Supporting information

S1 Supplementary figures

Figure S1: Probability for stable n-strain coexistence out of a pool of N = 10 strains.We plot the distribution of the number of coexisting strains n for 6 values of µ. For each value of µ, we generated 100 000 matrices A (normalized interactions between strains) and computed their stable steady states. The data are the same as those in the figures 3 and 5. We find that the distribution of n is very dependent on µ. In general, the probability is hard to compute explicitly but the limit µ → ∞ is explicitly known using the fact that the pairwise invasion fitness matrix becomes skew symmetric in that limit (see Box 2). One remarkable fact is that for small µ the probability of finding a stable steady state is greater for n = 1 and decreases very quickly with number of strains n, while for large values of n, the distribution looks like a binomial distribution centered at N 2 for odd values of n, and zero for even values of n. Hence, the probability to have a large number of strains coexisting at a given steady state is small for small µ and increases for large values of µ.

S2 Mathematical details for the limits of µ → 0 and µ → ∞

We use the following notations, for the basic reproduction number and for the global steady state obtained via the neutral system:

We define the important ratio µ of single and dual colonization by

.

(S1)

The reduced system reads żi

with the constraint

The qualitative behavior depends on the fitness N × N matrix Λ = (λ j i ) i,j whose i, j term is the fitness of strain i with respect to the strain j and is explicitly given by λ j i = µ(α j,i -α i,j ) + α j,i -α j,j . Here, A = (α i,j ) = ((K ij -k)/ε) is the normalized interaction matrix, accounting for deviation from neutrality. Recall that λ j i denotes the growth rate of the strain i when introduced to the steady state of the strain j alone :

The time scale of the selective dynamics between multiple strains, in this system, depends basically on the constant

Our goal here is to give a description of the system in the limit cases µ → 0 and µ → +∞. There are 3 parts in the discussion.

(i) Which parameters change? µ depends on k, β and m. This point is straightforward from the formula (S1).

(ii) The qualitative behaviour. This depends only on Λ (up to a positive multiplicative constant) which depends only on µ.

Hence the point (i) has no importance here.

(iii) How does the speed of the dynamics (Θ Λ ) vary? This speed is especially important when it is close to 0. Remember that this system has been obtained by focusing of the first order perturbation in the quasi-neutral approximation (the zero order perturbation, ε = 0, yields exactly the neutral model). If the speed of the dynamics is below some threshold O(ε), then we cannot remove the next order perturbation and our model does not describe well the dynamics. In other words, if the speed is too slow then the system will stay very close to neutrality for a very long time. We will say that the system is effectively neutral. Hence, the main question in this section is the following:

May the speed of the dynamics go to zero when µ changes? If yes, what is the relevance of the parameters in point (i)?

The qualitative dynamics. Recall that λ j i > 0 means that the strain i (the invader) may invade the system from a very small concentration within a system where only the strain j (the resident) is present.

Passing to the limit µ → 0, we obtain the matrix Λ lim with the (i, j) coefficient being simply

In general, very little may be said, because the structure of the matrix Λ lim can be anything. As an example, in the special case α jj = 0 for all j, then Λ lim = A T . Since A may have any structure, so can Λ lim , and very few results are known for the general replicator system. However, it is useful to see A as a random matrix and to check the probability of a given type of dynamics (see [START_REF] Yoshino | Rank abundance relations in evolutionary dynamics of random replicators[END_REF]). We have the following interpretation.

• The fitness λ j i depends only on the coefficients α j,i and α j,j which are characteristics of the resident. α j,i measures how much the resident j contributes to the pool of co-colonization by the strains i and j : D ij , whereas α j,j measures how much the resident j increases the pool of self co-colonization: D jj .

• Hence, the relative fitness of the invader i depends only of the resident j. The invader i may invade the resident j if and only if the resident reinforces more the co-colonization D ij than D jj .

• This phenomenon provides a lot of niches for multiple strains to coexist. As an example, if A has a strong diagonal : α j,j > α i,j for any j = 1, • • • N and i = j then all the strains may be stable residents when alone, and we will have at least N stable monomorphic steady states.

• In conclusion, let A be chosen randomly. The system is likely to have a lot of stable steady states wherein only few strains coexist. It is rare to see cycles or more complex dynamics.

The speed of the dynamics. Since µ → 0 then µ is bounded and far from 0, and so is Λ . The speed of the dynamics is then measured by

• If k → +∞, from the second expression in (S3) we have: Θ ∼ m 2k → 0 = Θ lim . This is a problematic situation. If the second infection occurs very quickly then Θ is very small and the model does not capture the real dynamics.

• If m → 0 then again the same expression give Θ ∼ m 2k → Θ lim = 0. If the natural mortality or the clearance of infection happen very slowly (which means that the rate of population turnover is very slow) then this model do not capture the real dynamics.

• If β → +∞ then the same equation above shows that Θ → m2k = Θ lim > 0 which, in this case, is far from zero.

Thus, if the pathogen transmission rate is very large, then the system captures well the real dynamics and the qualitative study of µ → 0 is appropriate.

In summary, the only way for Θ lim to be positive as µ → 0 is when β → +∞.

S2.2 Limit 2 : single colonization, µ → +∞ : k → 0 or R 0 → 1.

In that cases, we have generically Λ → +∞ and Θ → 0. In order to keep a bounded matrix Λ we rewrite the system as

The speed is then given by µθ and the qualitative behavior by µ -1 Λ.

The qualitative dynamics. When µ → +∞, the matrix µ -1 Λ → Λ lim = A T -A This matrix of pairwise invasion fitnesses becomes skew symmetric.

As opposed to the previous limit, there exist some results for the skew symmetric case. In the particular case where the coefficients are ±1 the system is known as the Tournament game which has been studied in detail in in [START_REF] Fisher | Optimal strategies for random tournament games[END_REF] and has been used recently in [START_REF] Allesina | A competitive network theory of species diversity[END_REF] in coexistence theory. Some of these results remain true for a general skew matrix [START_REF] Chawanya | Large-dimensional replicator equations with antisymmetric random interactions[END_REF].

• The quadratic term in the frequency evolution equation is always 1 zero. So the replicator equation system for N strains reduces to

• Generically, there exists exactly one non-negative linearly stable equilibrium (in particular multistability is impossible). Without loss of generality, we can assume that this equilibrium reads

• From the linear point of view, ζ is a center: it is linearly stable but not asymptotically linearly stable. Indeed, the N -n last eigenvalues are negative and the n first eigenvalues are purely imaginary.

• The function H(z) = N i=1 ζ i ln(z i ) is bounded and increasing. We have Ḣ → 0 and Ḣ(z) = 0 if and only z k = 0 for each k > n.

• It follows that z k → 0 for each k > n, and that the sub system of the n first strains is strongly persistent.

• In practice, exactly like the classical Lotka-Volterra model of prey and predators, we have a one-parameter family of cycles, parametrized by initial conditions, more precisely by the value of H(z(0)).

• Denote Λ n lim the sub-matrix of the first n rows and n columns of Λ ∞ and ζn = (ζ 1 , • • • , ζn). We have Λ n lim ζn = 0 In particular it is necessary that |Λ n lim | = 0. Since λ n lim is skew symmetric, we have

On the contrary, if n is even Λ n lim has no reason to be zero (except in very specific cases which have a zero probability of occurrence in the sense below) and there will be no steady state of n coexisting strains.

• A more complete description may be given. Let the coefficient of Λ lim be chosen randomly. Denote n the number of strains coexisting, then the probability to have a steady state of n = k strains over a total pool of N strains is

We illustrate this case in the figure S1.

We finish this section by highlighting the fact that the above-described dynamics is not structurally stable. This structure will be lost for a large but finite value of µ. As µ decreases from +∞, the center will become either a stable focus or an unstable focus. In the first case, we will observe a global attractor, while in the second case, which is the most probable as shown in Figure 5, we will observe more complex dynamics like heteroclinic cycle or even chaos.

The speed of the dynamics As shown below, the speed, in this second limit, is given by µΘ

• If R 0 → 1 then 2 µΘ → 0 = Θ lim and the model is effectively neutral.

If the transmission intensity of the infection is too low, then this model do not capture well the dynamics.

• If k → 0 then µΘ → Θ lim = m 2 (R 0 -1) > 0 which is far from zero.

In summary : the only way for Θ lim to be positive as µ → +∞ is to have k → 0, thus in the limit of extreme competition between strains in co-colonization. In that case, the compartment of co-colonization D becomes very small, and the dynamics are well described by the approximation µ → +∞, which depends only on the asymmetries α j,i -α i,j .

2 Except if β → +∞ and R 0 = β m → 1, like by example if β = m + a. This corresponds to a situation where both the mortalityclearance rate and the infection rate become simultaneously very large. In that case, we rewrite µΘ µ 2 (β-m) 2(µ+1) 2 -µ → θ lim = a 2 which is positive.