Multi-strain SIS dynamics with coinfection across heterogeneous patches

Multi-strain SIS dynamics with coinfection across heterogeneous patches | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Multi-strain SIS dynamics with coinfection across heterogeneous patches Miguel B. Marôco , Louis Pailloncy , Sten Madec , View ORCID Profile Erida Gjini doi: https://doi.org/10.1101/2025.10.26.684704 Miguel B. Marôco a Center for Computational and Stochastic Mathematics, Instituto Superior Tecnico, University of Lisbon , Lisbon, Portugal Find this author on Google Scholar Find this author on PubMed Search for this author on this site Louis Pailloncy b ENSTA , Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sten Madec c Institut Denis Poisson, University of Tours , France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Erida Gjini a Center for Computational and Stochastic Mathematics, Instituto Superior Tecnico, University of Lisbon , Lisbon, Portugal Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Erida Gjini For correspondence: erida.gjini{at}tecnico.ulisboa.pt Abstract Full Text Info/History Metrics Preview PDF Abstract We study a Susceptible-Infected-Susceptible (SIS) model with coinfection and multiple interacting strains where hosts move between a set of inter-connected patches. Under strain similarity and slow migration, we obtain a discrete model, following the corresponding continuous space model derived in ( Le and Madec, 2023 ). In this model, the fast variables are total prevalence of susceptibles, single-infected and co-infected hosts in each patch ( S, I, D ). The slow variables are strain frequencies ( z ) in each patch. These local strain frequencies in each patch are governed by a replicator-like equation, where an additional contribution arises from migration, scaled explicitly by emergent patch heterogeneity. In our model, the strains can vary epidemiologically within- and between-patches along several traits, including transmission rates, clearance rates, priority effects and pairwise susceptibilities to coinfection; a complexity that totally defies straightforward prediction of their ecological outcome. However, harnessing the analytical advantage of this model reduction, we can investigate several key scenarios for the outcome of the coupled P-patch N-strain system. We focus on the global and local factors that promote or hinder coexistence, and on regimes that favour spatial segregation of strains, or ultimately competitive exclusion. The analytical tractability of this meta-population model, illustrated with several examples for N = 2 strains and 2-3 patches, makes it a useful framework for application to general multi-species co-colonization systems with migration across heterogeneous environments. 1. Introduction Population dynamics in heterogeneous environments have many applications in biology, ecology, epidemiology and medicine ( Hanski, 1999 ; Okubo et al., 2001 ; Keeling et al., 2004 ). Typical modeling frameworks for exploring population dynamics across different habitats involve logistic growth models plus diffusion ( Arditi et al., 2018 ), Lotka-Volterra systems with migration between patches ( Levin, 1974 ), chemostat models ( Tan et al., 2019 ), consumer-resource models with symmetric/asymmetric diffusion ( Zhang et al., 2017 ), epidemiological compartmental models with different connected environments such as community-hospital, hospital-hospital, city X with city Y, rural with urban area, disease spread on networks, etc. ( Smith et al., 2004 ). In within-host infectious disease models, the question of population dynamics across different habitats also emerges when studying virus dynamics in different within-host niches ( Wong and Yukl, 2016 ), or bacterial spatiotemporal dynamics across different organs ( Grant et al., 2008 ). On the theoretical side, the questions arising in these models are those concerning the global persistence of the system ( Butler et al., 1986 ), the influence of diffusion on the global stability and magnitude of positive equilibrium ( Arditi et al., 2015 ), the comparison with the total abundance/dynamics in non-coupled system ( Poggiale et al., 2005 ), the influence of environmental heterogeneity for persistence and stability, and when different strains or species are modeled, how the competitive outcome between these players is altered by diffusion and environmental heterogeneity. On the more practical side, the application of such models to real data ( Trindade et al., 2022 ; Grant et al., 2008 ; Korevaar et al., 2020 ; Price et al., 2017 ), involves the challenges of suitable parametrization, parameter inference, uncertainty quantification and limits of predictability or control in such systems. In a broader sense, models of population dynamics in heterogeneous environments have a clear and direct bearing on global efforts towards conservation of biodiversity, health and stability of ecosystems and modern coexistence theory. In epidemiology, the study of spatiotemporal dynamics of infectious disease has a long history of study, first with reaction-diffusion type models which consider space as a continuous variable (see e.g. ( Noble, 1974 ) for early models of plague and ( Källén et al., 1985 ; Murray et al., 1986 ) for early models on the spread of rabies), and more recently with the application of multi-patch models (modeling space in discrete units), for example studying SIS dynamics ( McCormack and Allen, 2007 ; Arrigoni and Pugliese, 2002 ), application of aggregation methods under fast migration assumptions ( Kouokam et al., 2008 ; Marvá et al., 2012 ), vector-host epidemic models with multiple strains in a patchy environment ( Qiu et al., 2013 ) and multi-strain diseases on networks of interconnected patches ( Michalska-Smith et al., 2022 ). In general, spatial and non-spatial models differ, and often even the mere addition of a spatial dimension and random diffusion as a homogeneous variable to the non-spatial dynamics can have critical implications for the persistence of a growing diffusive species ( Skellam, 1951 ). Indeed, when space (or multiple patches) is heterogeneous, and diffusion is also space- or density-dependent, and when other parameters of growth and transmission can vary with space, the infectious disease dynamics become inevitably much more complex. In this work we will assume space is discrete. The motivation comes from the realization that many populations of humans, livestock, and wildlife consist of densely occupied subpopulations, or “patches”, connected by the movement of individuals, thus giving rise to a metapopulation. Today, most of the world’s population lives in cities, farm animals are concentrated in specific areas, and wildlife often cluster where their habitats are fragmented. In the past, people moved so little that new diseases usually stayed near where they began—like the Black Death in medieval Europe or smallpox centuries ago. Now, with our world more connected than ever, infections such as influenza, Ebola, and COVID-19 can spread faster and farther. As people travel more and science tracks pathogens in greater detail, disease models need to capture the complex links between how different areas, hosts, and strains interact. These connections ultimately shape how diseases spread and evolve, both locally and around the world. Here, we contribute to this challenge by considering a multi-patch multi-strain SIS system with coinfection, strain interactions, and symmetric host migration between patches. The basic model structure follows the Suscetible-Infected-Susceptible coinfection framework for N interacting strains proposed by Madec and Gjini (2020) , generalized by ( Le et al., 2023 ) and extended to continuous space in ( Le and Madec, 2023 ). In the non-spatial model ( Madec and Gjini, 2020 ; Le et al., 2023 , 2022 ), it was shown that under strain similarity, strain frequencies follow a slow dynamics given by the replicator equation: where denote pairwise invasion fitnesses between any two strains, and Θ gives the speed of the dynamics. In the subsequent epidemiological model extension ( Le and Madec, 2023 ), allowing for spatial diffusion of hosts, a similar model reduction was obtained, and in the case of low-diffusion a replicator-like equation was again derived for strain frequencies over space. In the present work, we will consider the same epidemiological multi-strain model but over discrete space, where we assume hosts can move between a set of fully-connected heterogeneous patches. Our aims are threefold: i) to present the equation governing strain frequencies over discrete patches and its properties, ii) to illustrate and analyze special cases of these dynamics for N = 2 in the case of P = 2 or 3 patches, iii) to suggest ways in which this model representation can be theoretically studied and extended. In a companion purely theoretical paper ( Madec and Gjini, 2025 ), we show the same model arrived here to via discretization of the spatial PDE ( Le and Madec, 2023 ), can be obtained by application from scratch of the slow-fast method on the epidemiological multi-strain multi-patch model. Indeed, this replicator equation, contains a local selection term, plus an explicit contribution term from host migration coupled with the environmental difference between patches, making it amenable to analysis and more efficient computation. 2. The model - metacommunity of hosts across P patches and N co-circulating strains 2.1. The SIS model with coinfection and host migration across patches We describe the model over a discrete set of P separate patches. Susceptible individuals are infected upon contact with infected individuals at a rate β . Single infected can become co-infected upon further contact with infected individuals, but with an altered susceptibility coefficient, relative to a susceptible host. The altered susceptibility parameters between strains are carried on, like in the earlier formulations ( Madec and Gjini, 2020 ), and introduced in the form of a matrix of coefficients k i j , which when above 1 indicate facilitation between strains i and j and when below 1, indicate competition. Single infected individuals and coinfected individuals recover from infection back to the susceptible class at equal global rate γ , but there can be differences in clearing infection by each strain γ i , or composite infections γ i j (for coinfection). There is no virulence so population size remains constant, with susceptibles having recruitment rate r , equal to the population per capita death rate. Strain dissimilarity is modelled through a parameter ϵ , which applies to all strain-specific parameters of the model. This formulation together with the biological meaning of all parameters is summarized in Table 1 . View this table: View inline View popup Download powerpoint Table 1: Strain similarity along 5 biological parameters. Notice that the model is general enough to encapsulate variation between patches both in the neutral model parameters β, γ, k, r, m (environmental mean-field heterogeneity) and in the non-neutral parameters , namely the relative variation between strains. However most of the time, we will assume explicitly only environmental heterogeneity between patches and study the net variation between strains across patches in terms of 4 effective λ coefficients, as shown in Eq. 12 . Note that all variables and parameters are to be understood as P-dimensional column vectors, with each entry being associated to one of the patches. With this, the equations of the model, describing host proportions in different compartments, translated to discrete P-dimensional space, are: with 𝒟 = ( d pk ) p,k ∈ ℝ P × P a Metzler matrix relating the connectivity between patches. For the total population in each patches to be constant (and renormalized to 1) we assumed that 𝒟(1 …, 1) T = (0, …, 0) T which reads for all p ∈ ⟦1, P ⟧: The force of infection for each strain i in a patch p is given by where: with ℙ ( j,i )→ i + ℙ ( j,i )→ j = 1, sums all transmission contributions of all hosts that can potentially transmit strain I within a patch k . The total single infection prevalence in patch p is given by . The total co-infection prevalence in patch p is given by . The total infection prevalence in patch p is given by T p = I p + D p . Like for the spatially continuous version of this model (see ( Le and Madec, 2023 )) we focus on a particular regime of the strain-dependent parameters and on the migration speed δ. Let 0 < ε ≪ 1 be a small parameter. We require two conditions for the model reduction to work. H1 Quasi-Neutral assumption In each patch p , the strain-dependent parameters are assumed to be ε -close. For instance, we write . See Table 1 for the notation on all the parameters. Note that this is an assumption of similarity between strains. There is no assumption of the variability of the parameters between patches. H2 Slow migration assumption We assume that the migration rate of hosts is ε -small : δ = εd for a fixed d > 0. Note that it is essential that the same ε appears in H1 and H2. 2.2. Discretization of the slow-fast approximation: explicit N strain-frequency dynamics across P patches In the spatially continuous system, under the assumption H1 and H2, through a slow-fast decomposition, it has been shown in ( Le and Madec, 2023 ) that, as ϵ → 0 each aggregated variable in each point in space approaches the equilibrium set by the neutral model ( ϵ = 0): with as the basic reproduction number. Note that these aggregated variables are indeed space-dependent. The model relies on R 0 ( x ) > 1 in each point in space. For the slow-variables z i , Le and Madec (2023) arrived at the following PDE equation, denoting the rates of change in strain frequencies over space and time z i ( x, t ): where the z i represent the frequencies of strain i in each point in space, with and with also space-dependent. Hence the vector affecting the advection term in the frequency equation, keeps a trace of the heterogeneity over space in mean-field transmission and other epidemiological parameters. We will now discretize the spatiotemporal equations in the case of N strains with slow diffusion δ = ϵd with ϵ ≪ 1, considering movement among P patches. During discretization, the continuous space variable x is replaced with the integer p ∈ [1, P ]. This way, z i ( t, x ) is replaced with , which gives the value of the frequency z for strain i in patch p (the t generally being implicit). Patch-specific global epidemiological variables Here we present the discrete version of the aggregated variables, ∀ p ∈ ⟦1, P ⟧, when t → +∞: From these global variables, we also have several other key patch-specific quantities: ∀ p (see Table 1 ). Patch-specific selection forces Now we can define the patch-specific pairwise invasion fitnesses between strains: ∀ i , ∀ j , ∀ p , and (pairwise invasion fitness matrix experienced in patch p ). Notice that these invasion fitnesses are weighted averages of strain variation along each of the 5 key dimensions, like in ( Le et al., 2023 ; Le and Madec, 2023 ). Considering the connectivity matrix 𝒟 between all patches, we can write the correspondence between the continuous replicator model and the discrete replicator equation: ∀ p ∈ ⟦1, P ⟧ and ∀ i ∈ ⟦1, N ⟧ following the steps in Table 2 . Therefore, the discretizations of (6) shows that the dynamics of the frequency of strain i in patch p is given by: which can be simplified as: View this table: View inline View popup Download powerpoint Table 2: Elements of the discretization of the PDE system, from continuous space to discrete space. In the above replicator-like equation plus the effect of movement, we have that the patch heterogeneity explicitly modulates the net contribution of host movement to strain dynamics, via the parameters ν pk : ∀ p , ∀ k measuring the relevant epidemiological ‘gradient’ between any two patches p and k in the system. The speed of local dynamics Θ p is also patch-specific (see Appendix A for more details on these parameters). Note that we have the later inequality coming from the fact that for any k . Let us denote the connectivity matrix whose coefficients measure the effective strength of links between patches. The facts that 𝒟 = ( d kp ) kp is an irreducible Metzler matrix and that 1 + ν kp > 0 imply that ℳ is an irreducible Metzler matrix too. With these notations, the multi-strain dynamics system reads compactly, for all strains i and patches p : This is the system that we study in the following sections. Population-level dynamics in each patch An immediate consequence of this fast-slow approximation framework for the SIS model with coinfection, similar to ( Madec and Gjini, 2020 ; Le and Madec, 2023 ), is that the final epidemiological compartmental variables in each patch p can be obtained after strain frequencies as: for single and co-infection respectively, which in our case will be patch-specific, first because the global prevalence compartments ( S p , I p , D p ) can be different between the patches, secondly because the strain frequencies can be different between the patches ( Eq. 6 ). Here, we have presented the discretization of ( Le and Madec, 2023 ), because it is easier to follow, but the same final equation can be formally obtained by deriving explicitly the slow-fast approximation and its regularization, directly from the SIS multi-strain multi-patch system. For a complete bottom-up derivation of this slow-fast approximation, see companion paper ( Madec and Gjini, 2025 ). Box 1. Results on the general N = 2 strains , P -patch model When N = 2 we have . Hence, denoting and z = ( z p ) p , the system (8) reads shortly for p ∈ ⟦1, P ⟧: We have for any p ≠ k . Then the systeme is cooperative. By the theorem 1.1 p.56 in Smith (1995) , the irreductibility of ℳ implies that the system is strongly monotonous (see theoreme 1.1 p.56). As a consequences, in the particular case of N = 2, we obtain some general properties that are inherited from the strongly monotone structure. This fact leads to several important properties. We refer to Smith and Waltman (1995) , where Chapter 6 presents results closely related to our current setting, and Appendix C provides a concise rationale for the properties of monotone systems. For further details on monotone systems, see Smith (1995) . Without being exhaustive, let us now state 3 key properties of system (10). There are no attracting limit cycles. Nearly all the trajectories converge to a steady state. To state the last properties, we note the simple fact that the trivial solutions 0 = (0, …, 0) T ∈ ℝ P and 1 = (1, …, 1) T ∈ ℝ P are both stationary solutions of the system (10), corresponding respectively to for all p , and for all p . If both 0 and 1 are linearly unstable, then there exist two steady states z * and z ** satisfying 0 < z * ≤ z ** < 1 such that the set attracts any trajectory z ( t ) with initial state z (0) ∈ [0, 1] P { 0, 1 }. In particular, if z * = z ** then this coexistence solution is unique and globally stable. If both 0 is unstable and 1 is stable then either their exists a non trivial coexistence solution or 1 is globally stable. Similar results hold under different stability assumptions on 0 and 1 . See Smith and Waltman (1995) , Chapter 6, for examples of such results in a similar framework. It is important to highlight the fact that the strong monotonicity of this system essentially allows one to reduce the study of its dynamics to the mere existence and stability of steady states. However, even though this reduces the complexity of the analysis, due to spatial heterogeneity, an exhaustive study of system (10) remains out of reach. This is why we restrict our study to the cases of P = 2 and P = 3 patches. The case of some identical patches By identical patches, we mean the same in all aspects. Two patches are called identical if they verify the same initial conditions, have the same parameters and strain species. It is an interesting case because it allows to perform a system reduction. If patches 1, 2, … m are identical, (for m ∈ ⟦2, P ⟧) from 8 it remains only to solve: ∀ p ∈ { m , …, P } and with And the m − 1 first patches follow the same dynamics as patch m : ∀ i, z i ,1 = z i,m −1 = z i,m . This shows that the N-strain P-patch network with m equal patches is equivalent to an N-strain (P-m)-patch reweighted network, allowing to reduce problem dimension. 2.3. The ecological scenario of 2 strains: N = 2 across P = 2 and P = 3 patches The case of 2 co-circulating strains is the easier one to tackle analytically ( Box 1 ) and often quite relevant for applications, so here, we will focus on this case in a metapopulation composed of 2 or 3 patches, or environments, where hosts can move. We will generally assume in this work that movement is uniform across patches, d pk = d for all p, k , however the model is general and encapsulates cases beyond this symmetric movement assumption. Regarding patch heterogeneity arising i) from mean-field global parameters ( ν pk ≠ 0, Θ p patch-specific), and from ii) patch-specific strain invasion fitnesses, we can say that this interplay of different heterogeneity layers will be the main focus of our investigation ( Figure 1 ). Download figure Open in new tab Figure 1: Schematic illustrating the context of the epidemiological SIS coinfection model for N = 2 strains and P = 2 patches. We show an epidemiological model with two patches A and B and 2 interacting strains, where hosts move in space and follow broadly SIS dynamics with coinfection. Assuming strain similarity,following ( Le and Madec, 2023 ) we can decompose the dynamics into a fast and slow component, and consider the special case of the reduced model in the limit relatively slow diffusion. In the discrete patch case, we arrive at a replicator-like equation for strain frequency in each patch. Besides local selection, this equation contains explicitly the contribution of patch heterogeneity and global migration, to the competitive dynamics between strains. 2.3.1. The model for P = 2 patches and N = 2 strains The replicator-like equation gets written in terms of the frequencies of one of the strains (e.g. strain , and for the case of 2 patches,it becomes: where we add the subscripts A and B, to the strain invasion fitnesses in patch A and B, respectively. 2.3.2. The model for P = 3 patches and N = 2 strains 3. Results for N = 2 It is clear that the role of ν pk , and the role of d pk are crucial for the difference with the ‘non-spatial’ model, or single-patch model. The equations for focal strain frequency in each patch have two clear components, a first component which is just the classical replicator dynamics (purely local selection dynamics in each patch), and a second component, which is the net benefit from global movement experienced by the focal strain in the current patch and in the rest of the patches (net metapopulation effect). While the speed parameter Θ is inconsequential in the single-patch model for determining final selection outcome, in the coupled P-patch system, the speed of the dynamics in each patch matters for final selection outcome: Θ p determine often both the magnitude and stability of the equilibria. Another immediate thing to notice is that under equal neutral mean-field parameters between patches, there is no global spatial heterogeneity in the system, ν pk = 0 for all p, k , and Θ p = Θ, so the contribution of host movement to strain dynamics in that case is purely driven by diffusion. Regarding the complexity of dynamics, we can say that while in the single-patch system, the magnitude of the unique coexistence equilibrium is very easy to compute ( Table 3 ) and its stability can be unambiguously decided on the basis of studies of the marginal exclusion equilibria ( 0 = (0, …, 0) and 1 = (1, …, 1)), in the coupled P-patch system, the coexistence equilibria have no straightforward analytic solutions, there can be multiple such interior coexistence equilibria, and hence the stability of interior coexistence points cannot be straightforwardly decided from studying simply the behaviour of the marginal exclusion equilibria 0 and 1 . Even though analysis of this system is harder, we can still advance in some special cases and derive explicit criteria. View this table: View inline View popup Download powerpoint Table 3: The final ecological outcomes for N = 2, in each patch, when uncoupled, depend purely on the mutual invasion fitnesses between 2 strains. 3.1. A sufficient condition of competitive exclusion for N = 2 strains and any number of patches P As explained in Box 1 , for N = 2 the system is strongly monotone, which is a powerful property. In this part, we focus on the case where the same strain wins in competitive exclusion in each patch. Let’s assume that without migration, strain 1 is dominant in each patch, that is: We will show that in that case 1 is globally stable. Indeed, under assumption (14), it is easy to show that 0 is unstable and that 1 is stable. By point ( iv ) of Box 1 , in order to prove that 1 is globally stable, it suffices to show that there is no non-trivial (interior) equilibrium. The system reads where we denote We proceed by contradiction. Assume that there exists a coexistence solution z ∈ (0, 1) P . Then z is a positive solution of wherein we have set 𝒜 ( z ) = d ℳ + diag( R p ( z p )). Like the matrix ℳ, the matrix 𝒜 ( z ) is an irreducible Metzler matrix. Firstly, (16) shows that z ∈ (0, 1) P is a positive eigenvector of 𝒜 ( z ) for the eigenvalue 0. By the PerronFrobenius theorem for Metzler matrices (see Bullo (2024) , Theorem 10.2), this implies that 0 is the principal eigenvalue of 𝒜 ( z ). Secondly, we also have ℳ 1 = 0 , which yields From assumption (14), we see that R p ( z p ) < 0 for each p , and thus 𝒜 ( z ) < 0. By Theorem 14.1 (iv) of Bullo (2024) , this implies that the principal eigenvalue of 𝒜 ( z ) is negative, which contradicts (16). This proves that when the same strain wins in competitive exclusion in each patch, it is the global winner also in the multi-patch system. 3.2. Results and simulations on the strains, patch model Although the equations governing strain selection dynamics across two patches look very simple, they can generate very complex behavior, - where up to nine possible equilibria are possible - which is not present in the single-patch system. Our scope is not to provide an exhaustive exploration of all cases, but rather to illustrate a few key scenarios where the coupling between patches can unequivocally favour exclusion or coexistence, and maintain or overturn the final outcome predicted by each single patch alone. We will focus only on analysis of selection dynamics between strains, given fixed parameters ( Eq. 6 ), as subsequent translation to epidemiology is straightforward via Eq.9. View this table: View inline View popup Download powerpoint Table 4: Possible outcomes for the ecological dynamics of 2 strains circulating across two heterogeneous patches connected via host movement. Outcome relative to a focal strain (here strain 1) frequency in patch A and patch B. We denote by (0) and (1) the outcome of single-patch dynamics favouring exclusion of strain 1 or persistence of strain 1. We denote by (0,0), (1,1) the case where the coupled patch dynamics favours competitive exclusion or competitive dominance of strain 1 in both patches. The scenarios of strain coexistence in single patch or in both patches are denoted by ( z ) and ( z A , z B ) respectively. Coupling (1) − (1) patches Without loss of generality, this is the simplest case to analyze, (equivalent to (0)-(0) coupling). In this case the same strain is a competitive single winner in either patch when uncoupled (monostability of same exclusion in each patch). By the general result in 3.1; (1, 1) is globally stable. Coupling ( z ) − ( z ′ ) coexistence patches This is also a simple case to analyze. In this case the two strains coexist in either patch when uncoupled (stability of of coexistence in each patch). This means both strains have positive individual invasion fitnesses in each patch: in patch A,and in patch B. It is straightforward to verify that that in this case, the two marginal equilibria (1,1) and (0,0) are both unstable, hence the only stable equilibrium of the system, independently of other patch-specific parameters, must inevitably be some interior coexistence point . Imposing the positive determinant and negative trace stability criteria on the Jacobian evaluated at each marginal equilibrium, leads to a contradiction, hence proving the impossibility of such scenario under this parameter regime. Coupling (0) − (1) patches: patch heterogeneity, invasion fitnesses and migration rate drive final outcome This is the most interesting case. Each patch alone favours the exclusion of a different strain. For illustration, we will study here a simplified scenario, assuming the invasion fitnesses are fixed free parameters and do not depend on mean-field parameters of the system. We assume an arbitrary fixed Λ matrix guaranteeing these patch-specific regimes. We examine numerically the stability of the two marginal equilibria (0,0) and (1,1) in the coupled system. Fixing values for λ i, j in each patch, and migration parameter d , we can broadly assume the following relations for epidemiological parameters between patches: R 0, B = cR 0, A and µ B = ωµ A , allowing to vary relative heterogeneity (see Appendix A for detailed parameter consequences). Under a certain value of c we can plot resulting qualitative regions of stability in the 2-patch system on the ( R 0 , µ )-plane, based on evaluating the stability properties of the marginal equilibria (0, 0) and (1, 1) ( Figure 2 ). Each region represents one of the possible outcomes for the system at each point ( R 0 , µ ) which are the parameters of patch A. Thus, in figure 2 we observe how the coupled 2-patch system can arrive at four possible broad selection outcomes depending on the basic reproduction number of patch A, R 0 , and the ratio of single to co-colonization µ = I / D (here assumed equal between patches). Download figure Open in new tab Figure 2: Outcome regions for the coupled 2-patch 2-strain system where alone patch A favours the exclusion of strain 2 and patch B the exclusion of strain 1, but patch B has higher total infection prevalence . Patch A parameters: , and strain invasion fitnesses . Patch B parameters: µ B = µ and , and invasion fitness coefficients . Migration rate between patches: d = 1, and demographic parameter for host turnover equal in each patch: m = 1.5. Notice that every time we change parameters of patch A, this changes also the parameters of patch B , despite the relative ratio between the two is held fixed. Inevitably, the precise shape of these regions (based on stability of marginal equilibria) does heavily depend on the choice of fixed parameters but they still provide good insight about the outcome of the system and can be used as a starting point for further numerical or analytical exploration. As a a function or patch variation in R 0 and µ we can obtain different regimes of system behavior. And within each of these regimes, there is a potential for underlying complexity, such as multiple coexistence points and multi-stability of several equilibria. See also Appendix D. For high values of the ratio µ coexistence is impossible, and the 2-patch system always approaches an exclusion, be it mono-stable or bistable. The magnitude of R 0 , which scales the total infection prevalence, on the other hand, decides the winning strain. In figure 2 , high values of R 0 benefit strain 1, whereas low values benefit strain 2. What is interesting to notice here is that unlike the simple case of exactly anti-symmetric invasion fitnesses within each patch (full analysis in Appendix B), for arbitrary four λ coefficients across patches, bistability of (0,0) and (1,1) also becomes possible in the coupled system. In Appendix D, we also provide more details and numerical explorations of the case in Figure 2 , including the effect of varying host movement rate d . Connecting two originally bistable patches: the (0/1) − (0/1) case This type of system naturally leads to more complex behavior, as the strains final outcome is frequency-dependent on initial conditions in each patch. In this case, it can be proven that the (0, 0) and (1, 1) global exclusion equilibria are both stable in the 2-patch system, implying the existence of unstable interior coexistence points, but without excluding the possibility of eventual stable coexistence points too. We illustrate such cases by plotting the nullclines of z A and z B for different assumptions of parameter variation between patches ( Figures 3 - 4 ). These figures show that the resulting behavior of the system depends heavily on the global patch heterogeneity, be it in basic reproduction number R 0 or single to coinfection ratio µ . In the R 0 variation case ( Figure 3 ), we see the most complex dynamics happens when R 0 in patch B is 1.15 times higher than the R 0 in patch A (middle scenario); there are 9 equilibria for this system, including 7 interior coexistence points. The marginal equilibria are still stable, however this does not guarantee global exclusion of either strain, because there are two interior coexistence points that are stable and could be reached depending on initial conditions. The other two values for c confirm the exclusive bistability of the two marginal equilibria. Download figure Open in new tab Figure 3: Connecting two bistable patches, under different R 0 gradient. The parameter c , different in each plot, varies the R 0 gradient between patches, . Left: c = 0.8 (2 stable equilibria) Middle: c = 1.15 (4 stable equilibria) Right: c = 2 (2 stable equilibria). In these plots the pairwise invasion fitnesses between the two strains are equal in each patch (equal ‘strength’ of bistability). The fixed parameter values are: . Stable/unstable steady states are denoted by the filled/empty circles, respectively. In the coinfection gradient case ( Figure 4 ), we also see a sensitivity to the relative magnitude of coinfection in each patch, where depending on parameter values, and initial conditions, bistability in each patch may drive ultimately interior stable coexistence everywhere in the system, instead of the classically expected exclusion. Download figure Open in new tab Figure 4: Connecting two bistable patches, under different coinfection gradient. The parameter w , different in each plot, varies the µ gradient between patches (ratio of single to co-infection) µ B = wµ A , keeping R 0 constant, here . Left: w = 0.5 (4 stable equilibria) Middle: w = 1 (4 stable equilibria) Right: w = 2 (2 stable equilibria). In these plots the pairwise invasion fitnesses between the two strains are equal in each patch (equal ‘quality’ of bistability). The fixed parameter values are: . In the scenarios above, we considered identical λ i, j values in both patches, meaning the bistability in each original patch is of equal magnitude - this led to a certain symmetry in the phase-plots, visible in the intersection structure between the two nullclines. When we relax this symmetry assumption, and instead allow for the magnitudes of the λ ’s to vary between strains and between patches, while still remaining negative, we see a broken symmetry in the corresponding phase-plots and having 2 additional stable interior coexistence points becomes more improbable. In such a scenario, we have also checked the effect of host movement rate on the final outcome of the system ( Figure 5 ). We see that in this asymmetric λ parameter regime and equal global epidemiological parameters, increased movement rate between patches reduces the chance strain coexistence, by reducing the number of alternative stable steady states in the system and pushing the system to the classical marginal equilibria: either one strain or the other everywhere. Download figure Open in new tab Figure 5: Effect of host movement rate d on strain dynamics between 2 locally bistable patches. Left: d = 0.8 (3 stable equilibria) Middle: d = 1 (3 stable equilibria) Right: d = 5 (2 stable equilibria). In these plots the pairwise invasion fitnesses between the two strains are unequal in patch A, and in patch B. The fixed parameter values are: . 3.3. Results and simulations on the strains, P = 3 patch model This case is more complex to study, as there is one extra environment where the strains can compete for host resources, and this will affect their net fitness and dynamics in the system. In the 3-patch case, there appears to be the two usual trivial equilibria: (1,1,1) and (0,0,0), and perhaps other coexistence ones, which are very complex to solve for since it comes to intersecting three cubics in 3D. Substitutions for and lead to finding the zeros of a degree 27 polynomial in , so there are at most 25 coexistence equilibria. Despite this complexity and high-dimensionality, one can notice straightforward outcomes of the dynamics in special cases. For example, in the case where pairwise strain fitnesses are proportionally-related by the same positive constant in each patch, and for some α > 0, the coexistence point is also an equilibrium of the system. Another result is that mixed trivial outcomes are not possible, for instance (1,0,1) or (0,0, z 3 ) because of the movement terms. Besides, some outcomes should be easily predicted. For instance, when connecting patches that have a stable exclusion steady-state on their own, coupled dynamics will lead to a global exclusion everywhere of that strain: (1), (1), (1) → (1, 1, 1) and (0), (0), (0) → (0, 0, 0), similar to the 2-patch case. Equal patches case For the sake of simplicity and illustration, we consider here the case of equal patches, in the sense of the neutral model. Some analysis can be performed in a straightforward manner (see Appendix C) and we can obtain a criterion for the stability of the marginal equilibria. Even if these results are partial, the stability of the exclusion state 0 and 1 says a lot about the dynamics thanks to the strong monotonicity of the system. Introducing aggregated parameters: , we have: encoding a balance between local stability and diffusion rate. Similarly, for the stability of (1, 1, 1), we introduce: . We then reach a condition very similar to the previous one: We simulate 3 key scenarios to illustrate system outcomes depending on the above stability criteria: i) (0,0,0) stable and (1,1,1) unstable ii) (0,0,0) stable and (1,1,1) stable hence bistability, and iii) (0,0,0) unstable and (1,1,1) unstable, leading to an interior coexistence point. In Figure 6A we show an illustration of Scenario 1, satisfying the stability criteria for one equilibrium but violating those for the other, and the dynamics converge to an exclusion state of strain 1 in all patches. In Figure 6B we show an illustration of Scenario 2, satisfying the stability of both exclusion equilibria, so depending on initial conditions, the dynamics converge to an exclusion of one strain or of the other, in all patches. In Figure 6C , we show an illustration of Scenario 3, not satisfying the stability criteria for either marginal equilibrium, leading to an interior coexistence state, albeit at a different frequency, in all patches. Download figure Open in new tab Figure 6: Examples of 2 strain dynamics across 3 patches. A . (0, 0, 0) stable and (1, 1, 1) unstable. Dynamics with ( z A (0), z B (0), z C (0)) = (0.9, 0.97, 0.95) B . Bistability of the marginal equilibria. Solid lines: Dynamics with ( z A (0), z B (0), z C (0)) = (0.61, 0.92, 0.43). Dashed lines: Dynamics with( z A (0), z B (0), z C (0)) = (0.61, 0.22, 0.43) C . (0,0,0) and (1,1,1) unstable. Dynamics with ( z A (0), z B (0), z C (0)) = (0.01, 0.12, 0.03). The system tends to a coexistence between the 2 strains in all three patches. Basic parameter values, equal among patches: β = 5, γ = 0.5, k = 0.1, r = 0.6, d = 1, Θ = 0.84, R 0 = 4.54. Pairwise invasion fitness between two strains are patch-specific: (A) Overall the simple numerical and analytical illustrations above for low-dimensional cases of the N -strain P −patch SIS model with co-colonization, serve to highlight the potential for rich nonlinear dynamics of such spatially-extended ecosystems, accessible via replicator-like dynamics ( Eq. 8 ) - where environmental heterogeneity together with local selective forces and host movement combine to generate final competitive outcomes between co-circulating strains. Ultimately space may favour or hinder coexistence, the answer is: it depends, and the study of this diversity entangled with environmental gradients should get more complex and fascinating for high number of strains. 4. Discussion Metacommunity ecology is increasingly recognizing and calling attention on the importance of local and regional dynamics for community outcomes ( Holyoak et al., 2005 ; Leibold et al., 2004 ; Leibold and Chase, 2018 ). Despite research efforts, the understanding of how different spatial scales interplay to shape metacommunities remains limited. Furthermore the inclusion of environmental heterogeneity into a unified theory of metapopulation dynamics, and the generality of modeled interactions beyond competition ( Leibold et al., 2004 ) are two key areas calling for a long time for mathematical analytical developments, beyond simulation ( Thompson et al., 2020 ). In this work, we have addressed this challenge, building an N − strain P − patch metacommunity model, based on SIS dynamics with co-colonization within patches and host migration between patches, incorporating patch heterogeneity and several dimensions of strain variability and interactions. Our model is fully-deterministic, and exploits the methods of time-scale separation under strain similarity ( Madec and Gjini, 2020 ; Le et al., 2023 ; Le and Madec, 2023 ). Although inspired by epidemiology, it is easily applicable to other endemic multi-species co-colonization scenarios, developing over a network of interconnected habitats. We studied the model in the case of of low diffusion, where global population structure endemic variables stabilize on a fast timescale (susceptibles S , singly-colonized I and co-colonized D host prevalences), and a replicator-like equation emerges for local strain frequencies in the slow timescale. The ODE model presented was based on discretization of the PDE framework, originally derived by Le and Madec (2023) , but it can also be obtained from first principles and re-application of the slow-fast approach on the multi-patch multi-strain SIS system ( Madec and Gjini, 2025 ). The power of this approach yields a formulation where patch heterogeneity, in the form of a relevant epidemiological gradient, elegantly appears in the master equation ( Eq. 8 ) as a net effective modulation of the movement term in the strain frequency ODE. While the replicator-like equation regulates local strain frequencies within each patch, it also ‘feels’ global processes operating at larger spatial scales, namely across all the patches in the network, where different prevalence of colonization or co-colonization can magnify or reduce effective strain fluxes across the metacommunity of hosts. Although we presented the general model formulation for any number of strains and patches, we illustrated the analytical tractability of this model through focusing on the N = 2 strain special case. In such case, the system is strongly monotonous, which considerably simplifies the study of the dynamics ( Box 1 , Section 3.1 ). In a network of 2 or 3 patches, the equations look apparently simple, but these cases are already able to generate very interesting and non-trivial dynamics, depending on isolated patch selection outcomes between 2 strains, connectivity, and heterogeneity parameters. While metacommunity frameworks typically focus on the situation of local competitive exclusion ( Leibold and Chase, 2018 ), here we have covered all possible local outcomes, from same-strain competitive exclusion, to opposite-strain competitive exclusion, both-patch coexistence, and bistability. The details of local selection become critical for global outcome especially if the connection is between two opposite-outcome patches or two bistable patches. Indeed, the final outcome depends in nonlinear ways on the balance between magnitude of local outcome stability, speed of local dynamics and its comparison with movement rate and its net modulation from other patch heterogeneity. Limitations and future directions Our model is deterministic, this means it cannot capture how demographic stochasticity can impact metacommunity dynamics ( Lerch et al., 2023 ; Capitán et al., 2017 ). Next model extensions could account explicitly for different population size between patches and investigate in detail the role of host stochasticity for multi-strain selection dynamics both locally and regionally. We expect in the limit of large population sizes, the model behavior will be similar to the trends and phenomena captured here, but we expect also major differences in small population scenarios, especially when founder-effects are operating at each patch (bistable-bistable connection), and in general a shrinking of coexistence regions in parameter regimes predicted by the deterministic model ( Neuhauser and Pacala, 1999 ). We also studied local pairwise invasion fitnesses between strains, within a simplified paradigm, namely as totally independent of the global mean-field parameters. However, if as in the original derivation by Madec and Gjini (2020) ; Le et al. (2023) , the invasion fitnesses are mechanistic model-based, and derived from the epidemiological traits between strains, the feedbacks between patch heterogeneity and strain selection could be different. In particular, we expect selection outcomes to be more constrained by such mechanistic couplings, where the λ ’s are not entirely free to vary. We have explored some of these scenarios, assuming the invasion fitnesses depend on mean-field parameters of the system, as outlined in Equation 7 , in Appendix E. What is fixed is the relative variation between strains. Indeed, if the pairwise invasion fitnesses between strains are trait-based, tracking selection of specific traits in a multi-patch model would become possible in the present framework. Because strain frequency evolution is explicitly governed by a replicator-like equation, we can link such selection dynamics directly to the evolution of certain traits locally in each patch and ultimately in the metapopulation; one of which could also be system resistance to invasion. Although the main goal of this paper is to provide a theoretical underpinning for multi-strain dynamics in a metacommunity of hosts in different patches, our models make also testable predictions in the form of bifurcation diagrams and critical shifts with certain parameters (e.g. Fig.2 . These predictions could be tested in laboratory microcosms or with epidemiological field data from co-colonization or co-infection systems with many strains. Indeed, a key challenge in the field of community ecology remains to understand to what extent observational data is informative and how it should be used to infer the underlying community assembly processes ( Ovaskainen et al., 2019 ). Examples of biological systems for applications could be dynamics of antibiotic-resistant and antibiotic-sensitive strains between two interconnected habitats, dynamics of two virus variants across two body compartments, pneumococcus serotype dynamics across geographic locations, co-colonizing species dynamics across two locations differing in temperature, humidity or other physical parameters. In such contexts, the questions of habitat fragmentation and biodiversity, climate change and biodiversity, antibiotic pressure gradient and selection in metapopulations could be explored further. The role of movement, or dispersal for biodiversity could be studied for specific cases, and the model framework provided here could help reconcile conflicting patterns, namely that sometimes movement could help promote biodiversity at the regional scale and other times it could hinder it. We believe mechanisms of global vs. local coexistence, captured here, can be put in practice for effective management of biodiversity and ecosystem health in a wide variety field of applications in the future, from agriculture to biomedicine. Acknowledgement We thank Thi Minh Thao Le for helpful discussions at the beginning of this work. The study was supported by the Portuguese Foundation for Science and Technology (FCT grant number 2022.03060.PTDC), by CEMAT (fellowship grant under UIDP/04621/2020 to M.M.) and by the ERASMUS EU programme (fellowship to L.P). Appendix A. Patch heterogeneity in mean-field parameters and ν pk Explicit parameter values: Θ, ν for each patch. In the 2-patch case, assuming without loss of generality that demographic turnover is equal in both patches ( m ) and the mean-field environmental difference between patches is given by basic reproduction number variation and single - to coinfection ratio variation: in the system above the speed parameters can be expressed as: and the global migration contribution parameters nu AB and nu BA are given by: which clarifies their dependence on global mean-field parameters in each patch and their ratios, but defies a straight-forward immediate understanding or prediction of effects. That is why we resort to simulations. Appendix B. Examples for 2 patch system for N = 2 Case of anti-symmetric invasion fitnesses in each patch. Let’s consider a simple case of: with t, x > 0 so that we keep the necessary fitness structure for opposite exclusion. We can control the stability difference of opposite strain exclusion in each patch by varying the coefficient t . We assume patches have equal ratios of single to co-infection µ ( ω = 1). The case c = 1: Patches have equal R 0 . Assuming equal patches in both R 0 and µ parameters: c = 1, ω = 1, we have the following. This becomes a very simple system of equations and can be completely described analytically. Besides the trivial global exclusion equilibria, there exists one feasible coexistence equilibrium which can be understood as a perturbation from the average point and if t = 1 then . The conditions for existence and stability of this point are identical and can be described as: So for fixed patch-specific parameters, there is a critical value of migration between patches d above which we go from a system of 2 exclusion equilibria (0,0) and (1,1), to a system with 3 equilibria including coexistence in both patches. This critical value is given by: . And for the marginal equilibria, where the same strain gets fixed in both patches, we have: and The conditions can be interpreted as comparing the relative strengths of stability of each exclusion equilibrium in each singular patch when uncoupled and identifying a critical lower bound for migration rate between patches. (0, 0) is the equilibrium where Patch A trend dominates ultimately in both patches in the coupled system, and is achieved for relatively higher stability of the z = 0 equilibrium in A vs. the z = 1 equilibrium in B, namely t > 1. This, together with a sufficiently high diffusion rate guarantees strain 1 loses in both patches, and strain 2 being the sole strain that persists in the coupled system. Similar arguments apply to the (1, 1) equilibrium: it is achieved for t < 1 and again its stability requires , These conditions, taken together indicate that strong diffusion between two originally opposite exclusion patches (sink - source) favours exclusion, eventually one or the other strain will be the sole winner in both patches. Conversely, for intermediate rate of diffusion between such patches, coexistence of both strains in both patches is possible and stable: In such coexistence, there can be biases in relative abundances, driven mainly by the asymmetric stabilities of each original exclusion in the respective patches under no coupling ( t ≠ 1). The above conditions imply that bistability between the marginal equilibria (0, 0) and (1, 1) is impossible under these assumptions. Furthermore, once there is coexistence, the marginal equilibria become unstable so there is also no bistability between a coexistence and an exclusion equilibrium. This means that the entire qualitative behaviour of the system can be completely described by the verification of the previous three conditions relating t, x, d and Θ. The case c > 1: Patch B has higher R 0 . When c > 1, hence when patch B has a bigger R 0 than patch A, we are still able to find closed form solutions for the coexistence equilibrium, but it is no longer straightforward to obtain the corresponding conditions for stability. The complications come from ν AB ≠ 0 and ν BA ≠ 0 anymore. Similarly, as before, this equilibrium is given as a perturbation from by and it should always exist. Because the direct conditions for its stability are not immediately tractable analytically, we can resort to reasoning about the stability of the marginal equilibria. Here, we have: and And so the outcome of the system might be predicted from these two conditions. These criteria demonstrate that the ultimate fate of the E1-E2 system will be decided by an interplay of all parameters together, including the migration rate d , the difference between patches in global mean-field epidemiological parameters manifested as ν AB and ν BA but also from the exact magnitudes of stability of exclusion in each patch, namely the parameters x and t . Notice also that even when the two patches are equal in some parameter, here µ , and as a consequence in Θ, this parameter, may still have an effect on global outcome; in this case Θ represents the speed of local dynamics and it interacts with d the speed of global mixing of the two strains in the system. Appendix C. The 2-strain 3-patch system with homogeneous patches The analysis performed in this special case requires the following additional hypothesis: All 3 patches have the same neutral model. In other terms, all patches are identical for fast timescales : β A = β B = β C , γ A = γ B = γ C … Only pairwise invasion fitnesses and slow variables z are now patch patch-specific. This framework enables in particular to take out the advection terms ν p = 0. In this context, the equation for the 3-patch system is: We now propose to study the stability of the two trivial equilibria of this system. Our two Jacobians are written on these points: We get two symmetric matrixcs. And this is where the additional hypothesis made in the beginning of the section makes sense (without it, the Jacobians would not be symmetric). Now, we want to see, under which conditions, each of these two Jacobians has strictly negative eigenvalues, which is equivalent to them being symmetric definite negative. Noticing this enables us to use an algebraic criterion, which is going to simplify a lot the study of stability : Sylvester’s criteria. Sylvester’s criteria Sylvester’s criteria: “Let 𝔸 = ( a i j ) 1≤ i, j ≤ n be a real symmetric matrix. Then, 𝔸 is positive definite if and only if the determinants of the n leading principal submatrices 𝔸 p = ( a i j ) 1≤ i, j ≤ p for p = 1, …, n are all strictly positive, that is, if all the leading principal minors of 𝔸 are strictly positive.” Proof . We will work in a vector space of dimension n 𝔼. Since 𝔸 is a real symmetric matrix, it is diagonalizable with transition matrix ℚ. Let us introduce its real eigenvalues ( λ ) i,i ∈1,… n , and their eigenvectors ( e i ) i ∈1,… n . This means we have 𝔸 = ℚ T ≱ℚ, where ≱ = Diag ( λ 1 , λ 2 , … λ n ). For each e i , let us say 𝔼 i = Vect ( e i ). Since 𝔸 is diagonalisable one can write: 𝔼 = 𝔼 1 ⊕ 𝔼 2 ⊕ 𝔼 n 𝔸 positive definite enables us to write: ∀ i ∈ ⟦ 1, n ⟧, λ i > 0. Now we are going to consider the reductions of endomorphism A on the vector spaces 𝔼 1 , 𝔼 1 ⊕ 𝔼 2 ,…𝔼 1 ⊕ 𝔼 2 ⊕… 𝔼 n . Each of them will guarantee the positivity of the determinant of a leading submatrix. For instance let us look at 𝔸 p = ( a i j ) 1≤ i, j ≤ p , for p ⟦1, n ⟧. In order to prove that its determinant is > 0, we have to look at the reduction of endomorphism 𝔸 on 𝔼 1 ⊕ 𝔼 2 ⊕ ….𝔼 p . Then we have that : det (𝔸 p ) = det (ℚ p ) 2 Π i ∈[1, n ] λ i > 0, where ℚ p = ( q i j ) 1≤ i, j ≤ p and the result is proven. By considering the same reductions of endomorphism 𝔸, we get that λ 1 > 0, λ 1 λ 2 > 0, … i ∈ ⟦1, n ⟧ λ i > 0. This leads to: ∀ i ∈ ⟦1, n ⟧, λ i > 0, which leads to 𝔸 definitive positive. Stability of (0,0,0) We are going to use this result in the following way : an equilibrium is stable if and only if all the leading principal minors of the opposite Jacobian matrix at this point are strictly positive. Here is our Jacobian: Its opposite is written: The criteria guarantees (0,0,0) to be stable if and only if: Now let us introduce: We have: (0,0,0) is stable for the network iff Now the following table predicting the outcome of a 2-strain scenario allows to make a relevant remark View this table: View inline View popup Download powerpoint Corollary . In particular if in all three isolated patches (0) is stable then is in , then (Φ A , Φ B , Φ C ) is in ℝ 3 , then all three conditions are verified, which implies that (0,0,0) is stable. In short our triple criteria enables to show the intuitive: ‘(0) stable in every isolated patch’ ⇒ ‘(0,0,0) is stable for the network’ Remarks It is worth mentionning that the three patches are interchangeable, meaning that their labels A, B and C can be taken arbitrarily. Besides, it can be interesting to see the impact of network connectivity on the “stability volume”. – Strongly connected networks : If d → + ∞ then the volume defined tends to the semi-space x+y+z¡0. (0,0,0) will be stable iff , that is iff (0) is “stable on average” – Weakly connected networks : If d → 0 then the volume defined tends to . At the limit the patches aren’t connected and so (0,0,0) is stable iff (0) is stable for every patch. Stability around (1,1,1) Now the Jacobian studied is:: In the same spirit as last section, we write its opposite in: Where we have introduced : We then reach a condition very similar than the previous one. (1,1,1) is stable for the network iff Corollary . In particular if in all three isolated patches (0) is stable then is in , then (Φ A , Φ B , Φ C ) is in , then all three conditions are verified, which implies that (1,1,1) is stable for the network. In short our triple criteria enables to show the intuitive: ‘(1) stable in every isolated patch’ ⇒ ‘(1,1,1) is stable for the network’. AppendixD. Supplementary figures and further illustrations Here, we test the accuracy of our predictions in Figure 2 of the paper, by picking points in each region and solving the system numerically ( Figure D.1 ). Download figure Open in new tab Figure D.1: Solutions of the system associated to different point of the ( R 0 , µ )-plane for parameters as in figure. The regions of bistability can arrive at either exclusion, depending on initial conditions. The Cases of (1, 1) and (0, 0) behave as expected, although the resulting dynamics could be more complicated. The coexistence region always maintains the diversity of strains over both patches, although the relative frequencies of each strain will change for the points inside the region. Finally, we can illustrate a sample solution for one of these points in the regions where more than one coexistence point exists ( Figure D.3 ), for R 0 = 3.2 and µ = 1.4. We can study the effect of host movement rate on this particular system by keeping every parameter as before but changing d ( Figure D.4 ). In this example, we can see that increasing the magnitude of movement d hinders diversity, in this particular system. The regions for coexistence shrink as d increases. These regions naturally do not allow us to understand whether the relative frequencies in these coexistence regions become more homogenized as d increases but it is clear that higher movement does not, for this set of invasion fitnesses, favour these two strains to coexist. Download figure Open in new tab Figure D.2: Regions of competitive outcomes for a special system with a test on uniqueness of a coexistence point. In the textured areas, there exists more than one coexistence point, meaning the behavior of the system can be more complicated. It is over these regions that marginal equilibria stability is not sufficient to capture all dynamics. Download figure Open in new tab Figure D.3: Verification of the complex behavior of the system when there exists more than coexistence point. (1, 1) is stable however there is another stable coexistence point. Depending on initial conditions (left: z A (0) = z B (0) = 0.2, right: z A (0) = 0.93, z B (0) = 0.89), the system will either fixate strain 1 or approach the stable interior coexistence. Strain 2, however, can never be fixated. Download figure Open in new tab Figure D.4: Effect of movement on the system with (0)-(1) patch connection, and parameters as in Figure 2 of paper. Left: d = 0.5, middle panel: d = 1, right: d = 1.5. As host movement increases the region of unique stable coexistence decreases. AppendixE. Mechanistic invasion fitnesses and predicting 2-strain outcomes in the ( c, ω )-plane In the main bulk of the paper, we have considered the invasion fitness to be 4 free parameters from their specific patch’s parameters. This is, however, not the most realistic definition, when recalling the derivation from the SIS model. Here we will analyze the effect that connecting invasion fitness to the environment has on the resulting dynamics of the system. The mechanistic definition of invasion fitness ( Madec and Gjini, 2020 ), assuming only trait variation in α i j , is and recall that in the non-spatial case the relative sign of the invasion fitness’ is enough to determine the competitive behavior of the system ( Figure E.5 ). With the usual framework, by fixing µ A = µ and varying ω we will be able to shift the system into heterogeneous regimes of competitive preference, that is, the patches will have different competitive outcome preference as we vary µ . This allows us to connect the invasion fitness’ to “environmental” properties of patches, in particular, to the ratio of single- to co-infected prevalcne µ . Of course, we still have to fix the set of α i j over both patches, but ultimately this can be done in a systematic manner directed at having different types of dynamics arise from sweeping the values of ω , that is, changing µ B in such that patch B changes its separate uncoupled outcome. Recall still that α i j in the full epidemiological model 4 is associated to the neutral altered factor of susceptibility k between strains and the strain dissimilarity ϵ As before, we will assume We decided to fix R 0 , as it does not appear explicitly in and μ ; and let c and ω vary. The motivation for this approach comes from the last section: in fixing µ we can keep the parameters of patch A fixed and, by varying ω , we are able to possibly change patch B into a different parameter regime ( Figure E.6 ). Furthermore, the type of selection change happens in a way which is connected to the “environmental” parameters and not taken arbitrarily, as before. That is, for every point in the ( c, ω )-plane there are a set of corresponding invasion fitness’; whereas in the paper (e.g. Figure 2 ) the four invasion fitness where fixed and constant over the whole ( R 0 , µ )-plane. We split the analysis into various qualitative cases of interest that arise from certain choices of α 11 , α 12 , α 21 and α 22 . This constraint will restrict the quality of patches we can connect. It is for example, impossible to connect a coexistence patch with a bistable one. Under this framework, not much analysis can be carried out analytically. Instead, we have used Mathematica to evaluate the usual Routh-Hurtwiz criteria for the stability of the marginal equilibria for points in the ( c, ω )-plane. We need to fix R 0 , µ, d ; and we will look into specific choices of α i j so as to enable different selection scenarios. We have also taken α i j ’s to be equal between the two patches. Throughout these various cases, we have kept the same values for the fixed parameters. These are resumed in Tables below: View this table: View inline View popup Download powerpoint Table E.7: Fixed parameters for each patch. (0) in patch A and patch B transitions from (0) to (1) . Here we have looked at a case where the invasion fitness functions allow the patches to shift (0), (0) to (0), (1) as a function of µ . Patch A is always of type (0) but as we vary ω we force patch B to become of type (1) ( Figure E.7 ). For a set of fixed We can take a closer look at how the transition from (0, 0) to (1, 1) happens in the ( ω, c )-plane and see there is actually a small boundary region of coexistence ( Figure E.8 .) Patch A is multistable and patch B transitions from multistability to (0) . This invasion fitness structure is similar as before, but patches move from a multistable competitive outcome (0 / 1), (0 / 1) to a (0 / 1), (0) competitive outcome ( Figure E.9 ). Here There is coexistence in patch A and patch B transitions from coexistence to (1). Taking the symmetric invasion fitness from the previous case allows to shift the patches from a coexistence coupling ( z ), ( z ′ ) to a ( z ), (1) coupling ( Figure E.10 ). That is (0) in patch A and patch B transitions from 3 outcomes. And finally, perhaps more interestingly, we can define invasion fitness function that allow for three different scenarios depending on the choice of ω . Now we have taken There are clearly 3 regions for µ , associated to three different types of coupling in Figure E.11 . µ ∈]0, 1.5[we have (0), (0) coupling; 2 µ ∈]1.5, 2[we have (0), (0 / 1) coupling; µ ∈]2, +∞[we have (0), (1) coupling; Depending on which of these regions of µ we are placed in the ( ω, c )-plane, the effect of c will be significant or not in the global coupled outcome. To conclude, the mechanistic connection of invasion fitness to environmental parameters as in the strict model derivation ( Madec and Gjini, 2020 ) reduces the possible outcomes of the system. Primarily Not all types of patch couplings can be achieved, that is, we can only connect certain types of patches whose invasion fitness functions shift from one competitive regime to another and functions do not allow for all possible regimes to be connected via µ .; For the possible couplings, variation in c and ω is enough to obtain different competitive outcomes The maximum number of possible outcomes found for a given system was 3. The 5-dimensional trait variation governing Λ as derived in ( Le et al., 2023 ) is interesting to be studied in the future. Download figure Open in new tab Figure E.5: Relation between the competitive outcomes of a single uncoupled patch and the invasion fitness functions of µ , for fixed α i j ’s. Download figure Open in new tab Figure E.6: Respective values for mechanistic invasion fitness under full model derivation, in µ and ω framework. Left: ω = 2.5, right: ω = 0.6 Download figure Open in new tab Figure E.7: Possible outcomes arising from coinfection susceptibility variation A . Left: Invasion fitness curves. Right:Regions of stability in the ( c, ω )-plane. The system can go from a (0, 0) to a (1, 1) outcome, meaning that the fixation of strains will depend on the size of ω and c . Download figure Open in new tab Figure E.8: The transition from (0, 0) to (1, 1) is done via a thin region of coexistence. Download figure Open in new tab Figure E.9: Possible outcomes for the chosen A . Left: Invasion fitness curves. Right: Regions of stability in the ( c, ω )-plane. Both multistability and 0, 0 are possible, depending on ω and c . Download figure Open in new tab Figure E.10: Possible outcomes from fixing A . Left: invasion fitness curves. Right. Regions of stability in the ( c, ω )-plane. We now can arrive at either or (1, 1 ). Download figure Open in new tab Figure E.11: Three possible outcomes from coupling. Left: Invasion fitness curves. Right: Regions of stability in the ( c, ω )-plane. Depending on ω and c we can arrive at either (0 , 0), multistability or (1, 1 ). Increasing c and ω favors the fixation of the first strain. Between the two mono-exclusions there is a region of multistability. Funder Information Declared Fundação para a Ciência e Tecnologia, https://ror.org/00snfqn58 , 2022.03060.PTDC , UIDP/04621/2020 European Commission, https://ror.org/00k4n6c32 , Erasmus Mobility Programme References ↵ Arditi , R. , Lobry , C. , Sari , T. , 2015 . Is dispersal always beneficial to carrying capacity? new insights from the multi-patch logistic equation . Theoretical population biology 106 , 45 – 59 . OpenUrl PubMed ↵ Arditi , R. , Lobry , C. , Sari , T. , 2018 . Asymmetric dispersal in the multi-patch logistic equation . Theoretical population biology 120 , 11 – 15 . OpenUrl PubMed ↵ Arrigoni , F. , Pugliese , A. , 2002 . Limits of a multi-patch sis epidemic model . Journal of mathematical biology 45 , 419 – 440 . OpenUrl PubMed ↵ Bullo , F. , 2024 . Lectures on Network Systems, version 1.7 Edition . CreateSpace , http://motion.me.ucsb.edu/book-lns . ↵ Butler , G. , Freedman , H. I. , Waltman , P. , 1986 . Uniformly persistent systems . Proceedings of the American Mathematical Society 96 ( 3 ), 425 – 430 . OpenUrl ↵ Capitán , J. A. , Cuenda , S. , Alonso , D. , 2017 . Stochastic competitive exclusion leads to a cascade of species extinctions . Journal of Theoretical Biology 419 , 137 – 151 . OpenUrl PubMed ↵ Grant , A. J. , Restif , O. , McKinley , T. J. , Sheppard , M. , Maskell , D. J. , Mastroeni , P. , 2008 . Modelling within-host spatiotemporal dynamics of invasive bacterial disease . PLoS biology 6 ( 4 ), e74 . OpenUrl CrossRef PubMed ↵ Hanski , I. , 1999 . Metapopulation ecology . Oxford University Press . ↵ Holyoak , M. , Leibold , M. A. , Holt , R. D. , 2005 . Metacommunities: spatial dynamics and ecological communities . University of Chicago Press . ↵ Källén , A. , Arcuri , P. , Murray , J. , 1985 . A simple model for the spatial spread and control of rabies . Journal of theoretical biology 116 ( 3 ), 377 – 393 . OpenUrl CrossRef PubMed Web of Science ↵ Keeling , M. J. , Bjørnstad , O. N. , Grenfell , B. T. , 2004 . Metapopulation dynamics of infectious diseases . In: Ecology, genetics and evolution of metapopulations . Elsevier , pp. 415 – 445 . ↵ Korevaar , H. , Metcalf , C. J. , Grenfell , B. T. , 2020 . Structure, space and size: competing drivers of variation in urban and rural measles transmission . Journal of the Royal Society Interface 17 ( 168 ), 20200010 . OpenUrl PubMed ↵ Kouokam , E. , Auger , P. , Hbid , H. , Tchuente , M. , 2008 . Effect of the number of patches in a multi-patch sirs model with fast migration on the basic reproduction rate . Acta Biotheoretica 56 , 75 – 86 . OpenUrl PubMed ↵ Le , T. M. T. , Gjini , E. , Madec , S. , 2023 . Quasi-neutral dynamics in a coinfection system with n strains and asymmetries along multiple traits . Journal of Mathematical Biology 87 ( 3 ), 48 . OpenUrl PubMed ↵ Le , T. M. T. , Madec , S. , 2023 . Spatiotemporal evolution of coinfection dynamics: a reaction–diffusion model . Journal of Dynamics and Differential Equations , 1 – 46 . ↵ Le , T. M. T. , Madec , S. , Gjini , E. , 2022 . Disentangling how multiple traits drive 2 strain frequencies in sis dynamics with coinfection . Journal of Theoretical Biology 538 , 111041 . OpenUrl PubMed ↵ Leibold , M. A. , Chase , J. M. , 2018 . Metacommunity ecology , volume 59 . Princeton University Press . ↵ Leibold , M. A. , Holyoak , M. , Mouquet , N. , Amarasekare , P. , Chase , J. M. , Hoopes , M. F. , Holt , R. D. , Shurin , J. B. , Law , R. , Tilman , D. , et al. , 2004 . The metacommunity concept: a framework for multi-scale community ecology . Ecology letters 7 ( 7 ), 601 – 613 . OpenUrl CrossRef Web of Science ↵ Lerch , B. A. , Rudrapatna , A. , Rabi , N. , Wickman , J. , Koffel , T. , Klausmeier , C. A. , 2023 . Connecting local and regional scales with stochastic metacommunity models: Competition, ecological drift, and dispersal . Ecological Monographs 93 ( 4 ), e1591 . OpenUrl CrossRef ↵ Levin , S. A. , 1974 . Dispersion and population interactions . The American Naturalist 108 ( 960 ), 207 – 228 . OpenUrl CrossRef Web of Science ↵ Madec , S. , Gjini , E. , 2020 . Predicting n-strain coexistence from co-colonization interactions: epidemiology meets ecology and the replicator equation . Bulletin of Mathematical Biology 82 ( 142 ). ↵ Madec , S. , Gjini , E. , 2025 . Derivation of a spatial replicator system with environmental heterogeneity from a co-colonization sis model with n strains and p patches. submitted . ↵ Marvá , M. , de La Parra , R. B. , Poggiale , J.-C. , 2012 . Approximate aggregation of a two time scales periodic multi-strain sis epidemic model: A patchy environment with fast migrations . Ecological Complexity 10 , 34 – 41 . OpenUrl ↵ McCormack , R. K. , Allen , L. J. , 2007 . Multi-patch deterministic and stochastic models for wildlife diseases . Journal of biological dynamics 1 ( 1 ), 63 – 85 . OpenUrl PubMed ↵ Michalska-Smith , M. , VanderWaal , K. , Craft , M. E. , 2022 . Asymmetric host movement reshapes local disease dynamics in metapopulations . Scientific reports 12 ( 1 ), 9365 . OpenUrl PubMed ↵ Murray , J. D. , Stanley , E. A. , Brown , D. L. , 1986 . On the spatial spread of rabies among foxes . Proceedings of the Royal society of London. Series B. Biological sciences 229 ( 1255 ), 111 – 150 . OpenUrl CrossRef ↵ Neuhauser , C. , Pacala , S. W. , 1999 . An explicitly spatial version of the lotka-volterra model with interspecific competition . The Annals of Applied Probability 9 ( 4 ), 1226 – 1259 . OpenUrl ↵ Noble , J. V. , 1974 . Geographic and temporal development of plagues . Nature 250 ( 5469 ), 726 – 729 . OpenUrl CrossRef PubMed Web of Science ↵ Okubo , A. , Levin , S. A. , et al. , 2001 . Diffusion and ecological problems: modern perspectives . Vol. 14 . Springer . ↵ Ovaskainen , O. , Rybicki , J. , Abrego , N. , 2019 . What can observational data reveal about metacommunity processes? Ecography 42 ( 11 ), 1877 – 1886 . OpenUrl ↵ Poggiale , J.-C. , Auger , P. , Nerini , D. , Manté , C. , Gilbert , F. , 2005 . Global production increased by spatial heterogeneity in a population dynamics model . Acta Biotheoretica 53 ( 4 ), 359 – 370 . OpenUrl PubMed ↵ Price , D. J. , Breuzé , A. , Dybowski , R. , Mastroeni , P. , Restif , O. , 2017 . An efficient moments-based inference method for within-host bacterial infection dynamics . PLOS Computational Biology 13 ( 11 ), e1005841 . OpenUrl ↵ Qiu , Z. , Kong , Q. , Li , X. , Martcheva , M. , 2013 . The vector–host epidemic model with multiple strains in a patchy environment . Journal of Mathematical Analysis and Applications 405 ( 1 ), 12 – 36 . OpenUrl CrossRef ↵ Skellam , J. G. , 1951 . Random dispersal in theoretical populations . Biometrika 38 ( 1/2 ), 196 – 218 . OpenUrl CrossRef PubMed Web of Science ↵ Smith , D. L. , Dushoff , J. , Perencevich , E. N. , Harris , A. D. , Levin , S. A. , 2004 . Persistent colonization and the spread of antibiotic resistance in nosocomial pathogens: resistance is a regional problem . Proceedings of the National Academy of Sciences 101 ( 10 ), 3709 – 3714 . OpenUrl Abstract / FREE Full Text ↵ Smith , H. , 1995 . Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems: An Introduction to the Theory of Competitive and Cooperative Systems . Mathematical surveys and monographs. American Mathematical Society . URL https://books.google.pt/books?id=vOfNAwAAQBAJ ↵ Smith , H. , Waltman , P. , 1995 . The Theory of the Chemostat: Dynamics of Microbial Competition . Cambridge Studies in Mathematical Biology . Cambridge University Press . URL https://books.google.pt/books?id=wFLdVo89vq8C ↵ Tan , C. , Wang , Y. , Wu , H. , 2019 . Population abundance of a two-patch chemostat system with asymmetric diffusion . Journal of Theoretical Biology 474 , 1 – 13 . OpenUrl PubMed ↵ Thompson , P. L. , Guzman , L. M. , De Meester , L. , Horváth , Z. , Ptacnik , R. , Vanschoenwinkel , B. , Viana , D. S. , Chase , J. M. , 2020 . A process-based metacommunity framework linking local and regional scale community ecology . Ecology letters 23 ( 9 ), 1314 – 1329 . OpenUrl CrossRef PubMed ↵ Trindade , S. , De Niz , M. , Costa-Sequeira , M. , Bizarra-Rebelo , T. , Bento , F. , Dejung , M. , Narciso , M. V. , López-Escobar , L. , Ferreira , J. , Butter , F. , et al. , 2022 . Slow growing behavior in african trypanosomes during adipose tissue colonization . Nature communications 13 ( 1 ), 7548 . OpenUrl PubMed ↵ Wong , J. K. , Yukl , S. A. , 2016 . Tissue reservoirs of hiv . Current Opinion in HIV and AIDS 11 ( 4 ), 362 . OpenUrl ↵ Zhang , B. , Kula , A. , Mack , K. M. , Zhai , L. , Ryce , A. L. , Ni , W.-M. , DeAngelis , D. L. , Van Dyken , J. D. , 2017 . Carrying capacity in a heterogeneous environment with habitat connectivity . Ecology letters 20 ( 9 ), 1118 – 1128 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted October 27, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Multi-strain SIS dynamics with coinfection across heterogeneous patches Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Multi-strain SIS dynamics with coinfection across heterogeneous patches Miguel B. Marôco , Louis Pailloncy , Sten Madec , Erida Gjini bioRxiv 2025.10.26.684704; doi: https://doi.org/10.1101/2025.10.26.684704 Share This Article: Copy Citation Tools Multi-strain SIS dynamics with coinfection across heterogeneous patches Miguel B. Marôco , Louis Pailloncy , Sten Madec , Erida Gjini bioRxiv 2025.10.26.684704; doi: https://doi.org/10.1101/2025.10.26.684704 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Ecology Subject Areas All Articles Animal Behavior and Cognition (7618) Biochemistry (17633) Bioengineering (13856) Bioinformatics (41841) Biophysics (21399) Cancer Biology (18529) Cell Biology (25422) Clinical Trials (138) Developmental Biology (13352) Ecology (19860) Epidemiology (2067) Evolutionary Biology (24282) Genetics (15582) Genomics (22462) Immunology (17700) Microbiology (40295) Molecular Biology (17140) Neuroscience (88419) Paleontology (666) Pathology (2823) Pharmacology and Toxicology (4813) Physiology (7632) Plant Biology (15107) Scientific Communication and Education (2042) Synthetic Biology (4284) Systems Biology (9808) Zoology (2267)