Predicting N-Strain Coexistence from Co-colonization Interactions: Epidemiology Meets Ecology and the Replicator Equation

Multi-type infection processes are ubiquitous in ecology, epidemiology and social systems, but remain hard to analyze and to understand on a fundamental level. Here, we study a multi-strain susceptible-infected-susceptible model with coinfection. A host already colonized by one strain can become more or less vulnerable to co-colonization by a second strain, as a result of facilitating or competitive interactions between the two. Fitness differences between N strains are mediated through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^2$$\end{document}N2 altered susceptibilities to secondary infection that depend on colonizer-cocolonizer identities (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{ij}$$\end{document}Kij). By assuming strain similarity in such pairwise traits, we derive a model reduction for the endemic system using separation of timescales. This ‘quasi-neutrality’ in trait space sets a fast timescale where all strains interact neutrally, and a slow timescale where selective dynamics unfold. We find that these slow dynamics are governed by the replicator equation for N strains. Our framework allows to build the community dynamics bottom-up from only pairwise invasion fitnesses between members. We highlight that mean fitness of the multi-strain network, changes with their individual dynamics, acts equally upon each type, and is a key indicator of system resistance to invasion. By uncovering the link between N-strain epidemiological coexistence and the replicator equation, we show that the ecology of co-colonization relates to Fisher’s fundamental theorem and to Lotka-Volterra systems. Besides efficient computation and complexity reduction for any system size, these results open new perspectives into high-dimensional community ecology, detection of species interactions, and evolution of biodiversity. Electronic supplementary material The online version of this article (10.1007/s11538-020-00816-w) contains supplementary material, which is available to authorized users.

1 Illustrative model parameters for a biological system Table S1: Model parameters, some illustrative values, and their biological interpretation (e.g. in pneumococcus bacteria transmission). Notice that the model is general enough to deal with any combination of parameters, provided they lead to a global endemic steady state and close values of co-colonization traits between strains. R 0 Basic reproduction number, R 0 = β/m, equal among strains R 0 < 1 Disease extinction R 0 > 1 Persistence 2-3 (Dietz, 1993;Gjini et al., 2016) Gjini, E., Valente, C., Sá-Leão, R., and Gomes, M. G. M. (2016). How direct competition shapes coexistence and vaccine effects in multi-strain pathogen systems. Journal of Theoretical Biology, 388:50-60.

Derivation details for the slow-fast dynamics
First, we rewrite the system (3)-(5). With the perturbation from neutrality K ij = k + εα ij and the aggregated variables the strain-specific force of infection becomes F i = βJ i and the system (1) reads For ε = 0 we have the neutral model with the limits when t → +∞ : If ε > 0, the first two limits remain true. It follows that in order to apply the Tikhonov theorem on slow fast dynamics (See Tikhonov (1952) Systems of differential equations containing small parameters in the derivatives. This reference and a few of more recent references on this subject are also cited in the main text), we may take S = S * and T = T * (these limits being reached uniformly fast in ε). In contrast, the equation of I(t) depends on ε. Hence we use the ansatz: Finally, remark that the equations on I ij have no influence on the equations for I, I i and J i . With these considerations and notations, we obtain wherein we have denoted: The final slow-fast form is obtained by using the new variables H i and z i defined in (8) in the main paper. Hence, define: Finally, using these new unknowns, one obtains the expression of the slow-fast system (S3): where δ = β 2(T * ) 2 −I * D * is a positive constant and −ξ is the non-zero eigenvalue of matrix A 0 (see Eq. (7) in main text).

Characterization of the model via invasion fitnesses
We focus on the trivial solutions of (10), which are the extreme competitive exclusion equilibria In particular, we note λ j i the invasion fitness of the strain i with respect to the trivial solution E j . More precisely, λ j i is the per-capita growth rate of the strain i evaluated at E j . Since the i-th line of the system reads d dτ z i = z i h i (z 1 , · · · , z n ), we have the explicit formula: With the notations introduced above, this can be rewritten as: Finally, the fact that Plugging this into (S9), we obtain a new version of system (S7), enterely parametrized by the fitnesses λ j i : (S10) or more shortly by noting λ i i = 0 and the pairwise invasion fitness matrix (S11)

Symmetric invasion matrix
In this case, starting with λ j i = λ i j , for all pairs of strains the following holds: This case resembles the classical condition for coexistence in ecology where the difference in withinstrain interactions being bigger than that in between-strain interactions, favours coexistence between the strains, independently of whether the mean trait is cooperation or competition: the important principle is about deviation from the mean. If two strains cooperate more between than within-strain, then they will coexist. Similarly, if two strains compete more within strain than between strains, they will coexist. In general, in this case, we have that Q always increases:Q ≥ 0. Indeed, denoteλ i = j λ j i z j . We haveQ We have Q = z i λ i (i.e. mean invasibility trait in the system) and then • Special case of all λs negative The system will be characterized by multistability and competitive exclusion. Q will be negative and will increase towards zero. (similar too the case of defense-driven invasion).
• Will only the pairs with λ j i = λ i j > 0 coexist? If λ j i = λ i j < 0 then the coexistence steady state E i,j exists but it is unstable. For a system with many strains, the outcomes can be complex.

Invader-driven invasion matrix
For this case, we find that there are many possible coexistence solutions subject to the constraint: This leads to z i = 1 − Q λ i at steady state. Since z i < 1 this yields that for all the strains i which coexist, Q and λ i have the same sign. So the invasion fitnesses of the coexisting strains have the same sign. If there is a steady state with a negative Q, and at least one member strain j with positive invasion fitness, then the steady state is unstable for the growth rate of j being in that case λ j − Q > 0. This implies that for the steady state to be stable, either all the strains have negative fitness or the coexisting strains have positive fitness. In the first case, we have a multistability of competitive exclusion steady state and Q goes to 0. In the second case, then from z i > 0, we have λ i > Q at equilibrium, and this means only the strains with a large enough intrinsic invasion fitness λ i will coexist.
In detail, at steady state with n coexisting strains, one obtains Q = n−1 j 1 λ j , and thus for coexisting equilibrium strain frequencies:

Resident-driven invasion matrix
In this case the equation for z i can be written as: If a strain has negative invasion fitness, e.g. λ 1 < 0 (good defender of his equilibrium when alone), then the corresponding trivial steady state E 1 = (1, 0, · · · , 0) is asymptotically stable (because the linear growth rates of all other strains evaluated at such equilibrium, will always be negative. Hence, this situation leads often to a scenario of competitive exclusion with a single winner (and thus Q(t) → 0). Several strains may have this property which leads to multistability of exclusion. Depending on initial conditions, there may exist some strain, for which j =i λ j z j will be maximal compared to other strains in the system (the space of 'favourable' permissive territory is largest), and together with its initial frequency, this will drive the disproportionate growth of that strain in the system, leading to competitive exclusion, thus making Q = 0 in the long term. If the values of resident-driven fitnesses span the positive and negative range, and initial frequencies vary, different competitive exclusion scenarios may emerge, thus reinforcing the possibility of multi-stability.
However, there are exceptions. If a subset Ω of n strains have positive invasion fitnesses, they may coexist at a steady state E I and from Q = j =i λ j z j and z j = 1 we infer that at such steady state: . But if this subset Ω is strict (n < N ), then this steady state E Ω is unstable because the other strains i 0 / ∈ Ω, at such equilibrium, will experience the positive growth rate For example, if all the strains have a positive invasion fitness (everybody is a bad 'defender', thus everybody is invadable) then a stable coexistence of all N types in the system is possible.
To conclude, for random λ values which can be positive or negative, this type of matrix typically leads to a multistability of competitive exclusion steady states, although for very special cases, coexistence may arise.

Anti-symmetric matrix
Here we have for each couple of strains in the system λ j i = −λ i j . This leads to : meaning that the cumulative (summed) deviation from the mean in within vs. between-strain interaction coefficients is equal for any two strains in the system, making the dynamics more difficult to settle, and thus complexity is more likely to arise in the multi-player community. For Q, we may write: In that case, it becomes immediately clear that Q = 0 at any time, independently on the values of type frequencies z i .

Supplementary figures and movies of model dynamics
Supplementary Movie S1. See separate file movie1.mp4. Here we illustrate the dynamic trajectory of the system (Figure 2 of the paper). A. Network of pairwise invasion fitnesses between strains, where edges can be of 4 types: coexistence (red), bistability (blue), exclusion of i, exclusion of j (directed gray arrows toward the winner). We plot the trajectory by considering the z = (z 1 , · · · , z 6 ) as the barycentric coordinate of a point M into the regular polygon determined by the nodes of the graph. B. Trajectory of the strain frequencies over time. Figure S1: Illustration of more complex dynamics. Here we simulate the model in its reduced form for another set of parameters, leading to limit cycle coexistence between 5 strains in an N = 6 system. See associated movie S2. The precise co-colonization interaction matrix is the following :  Figure S1). A. Network of pairwise invasion fitnesses between strains, where edges can be of 4 types: coexistence (red), bistability (blue), exclusion of i, exclusion of j (directed gray arrows toward the winner). The system trajectory z i is shown as a movement in a 6-dimensional space. B. Trajectory of the strain frequencies over time. The constant Θ is the rate that sets the tempo of multi-type "motion" on the slow manifold towards an equilibrium. Θ depends specifically on the absolute transmission rate of the pathogen β, which calibrates the slow dynamics, but also on mean traits R 0 and k, via the conserved aggregated quantities T * , I * , D * , where R 0 and k appear in nonlinear coupling with each other. There is an intermediate value of transmission intensity that maximizes the speed of strain frequency dynamics, and this increases with more competition between strains (k → 0). b) Θ as a function of both R 0 and k. The speed of slow dynamics is maximal for intermediate values of R 0 and relatively low values of k. It is when competition dominates in co-colonization processes, that prevalence of singly-colonized hosts increases most in the system, speeding up the competitive dynamics between strains in the longer time scale. Here β = 2. Figure S3: Illustration of quasi-neutrality principles on the slow manifold for N = 20. a) Strain frequencies tend to equalize in single and co-colonization, for all strains and for all time when on the slow manifold. b) Prevalence of co-colonization I ij tends to scale with the product of strain prevalences in single colonization (I i , I j ), for all strain pairs and for all time, when on the slow manifold. This example is simulated using a random matrix K, with N = 20.
Each trajectory corresponds to a given strain in the system (a), or a given strain pair (b). Figure S4: Accuracy of the slow-fast approximation over time We confirm there exists t 0 > 0 such that for any T > t 0 and any ε > 0 small enough, there is M > 0 such that for any τ ∈ (t 0 , T ), the error between the two is small: Error between the neutral model and the full system, starting from random initial conditions. We summarize simulations for different ε, where for each ε, 20 co-colonization matrices of different size (N ∈ [2, 15]) were randomly generated, together with initial conditions. On a fast time-scale (o(1/ε)) the neutral model is a very good approximation to the real system. b) Error of the neutral and the slow dynamics approximation, starting from the slow manifold (9). The error (mean over all model variables) between the original system trajectories and the neutral vs. slow dynamics. We generated random co-colonization matrices (N ∈ [2, 15]) for 3 values of ε (10 −1 ,10 −2 ,10 −3 ) while keeping mean k = 1. The K ij were drawn from the normal distribution N (k, ε 2 ). On the slow manifold, the slow-dynamics approximate well the original system, unlike the neutral model. Figure S5: Computation speed of the approximation. Here, we verify that the slowfast approximation enables exploration of dynamics in a much more computationally efficient manner over a very large number of strains N , which is difficult with the original system. Computation time increases as a function of number of strains N in the full system (black dashed line) but not in the slow-fast dynamic aproximation (blue solid line) for T = 5000 time units. The line shows the mean over 20 random matrices for each N , with ε = 0.1, k = 1 and R 0 = 2. For details see Figure S6. Figure S6: Computation speed of the slow-dynamics approximation in detail. Here we plot the computation time as a function of number of strains N , in terms of number of seconds on a PC to simulate the slow-fast dynamic aproximation (boxplots) for T = 5000 time units starting from the slow manifold. The line shows the mean over 20 random matrices for each N , where ε = 0.1 and mean co-colonization coefficient k = 1 and R 0 = 2. Note that we use the default representation of the function boxplot of matplotlib.pyplot. In particular, for each value of N , the box is delimited by the first and third quartile, Q1 and Q3. The orange line inside the box correspond to the median Q2. The whiskers are going from Q1-1.5(Q3-Q1) to Q3+1.5(Q3-Q1). The values outside this interval are represented with dots.