Multi-strain disease dynamics on metapopulation networks

A much-updated version of this work is now available published


Introduction
Many of the most impactful infectious diseases that affect humans (influenza, malaria, human papillomavirus, etc.), livestock (porcine reproductive and respiratory syndrome, foot-and-mouth disease, etc.), and wildlife (anthrax, plague, etc.) have clusters in their population-genetic variability that we classify as strains.This variation in pathogen genotype is often associated with differences in phenotype, which can have profound effects on the efficacy of host immune defenses.While the human immune system is usually capable of preventing re-infection i.e. infection with something to which it has been previously exposed sufficient, divergent evolution among pathogen strains can reduce the ability of the host to recognize, and thus mount an immunological response to, subsequent exposures.In some cases, this change is not sufficient to completely avoid recognition by the host's immune system, yielding an immune response that is neither as strong as it would be in the case of re-exposure to the same strain, nor as weak as in the case of exposure to a novel pathogen.This partial cross-reactive immunity can likewise lead to reduced transmissibility, affecting disease dynamics across the population.
While the study of multi-strain diseases goes back decades (1; 2), the resulting modelling framework has not yet been generalized to a collection of sub-populations connected through host movement, i.e. a metapopulation (but see (3)).Initially introduced through the concepts of island biogeography (4), the network approach of metapopulations can be applied to a variety of systems, including human movement between cities, livestock transport between farms, and wildlife living in fragmented natural habitats.In each case, there exist relatively high-density areas which are connected to one another through a network of individuals' movement.A metapopulation framework allows the application of network analyses to characterize patterns of connection within the larger system, and can provide unique insights across scales.
Historically, metapopulation studies have been been divided into two main camps: those that model within-population dynamics and "cell occupancy" models.The latter of which, where only the presence or absence of a given species within a population is recorded in a given timepoint (5), has received much more theoretical attention.Importantly, cell occupancy models rest on an assumption of temporal separation in which local dynamics occur on a timescale that can be treated as instantaneous relative to that of the between-population dynamics (6).When considering pathogens in systems with relatively high migration rates, however, this assumption rarely holds, and the presence-absence approach can significantly limit model accuracy (7; 8; 9).The presence of metapopulation structure has been repeatedly associated with increased stability (10; 11; 12).This is due in part to the ability of migration between asynchronous populations to rescue temporarily low density populations from extinction (13).This is particularly relevant when populations are undergoing cyclical or chaotic dynamics, where repeated instances of low density are generally considered to be at greater risk of extinction than a population maintaining steady state dynamics (14; 15).
Here, we build on the strain theory of host-pathogen systems proposed by (16), considering a scenario where a collection of populations undergoing local dynamics are furthermore interconnected through the movement of individuals between populations.We simulate disease dynamics on this system, characterizing the effects of parametrization and network structure on these dynamics.This work is divided into three sections: first, we explore a case of interconnected populations which have been parametrized to display identical dynamics in the absence of host migration.Second, we consider cases where parameters differ between populations.Finally, we explore the role of network structure on disease dynamics in larger networks of connected populations.

Model framework for one population
We work from a system of ordinary differential equations which delineate a population into classes based on current and past exposure to different strains of a pathogen.Pathogens with strain structure can differ in both the number of strains and the level of cross-reactive immunity afforded by past exposure to similar strains.To model the number of strains, we signify a strain i = {x 1 , x 2 , . . ., x n } as a set of n loci, each of which can take on a finite number of alleles.For instance, a pathogen with two loci (a and b) and two alleles at each loci has a total of four potential strains: {a For cross-reactive immunity, we use a parameter γ which indicates the degree of reduced susceptibility a host has to strains that are similar to (i.e.strains that share at least one allele with one another) past exposures.Importantly, in this framework, the number of strains is fixed and finite.While strains may go extinct over time, there is no process for the generation of new strains or to re-introduce strains that had previously gone extinct (but see ( 16)).
The model consists of sets of three nested equations (one set for each strain i): y, z, and w.See (17) for a more comprehensive discussion of the model framework, including a graphical representation.y i represents the proportion of the population currently infectious with strain i. z i represents the proportion of the population that has been exposed to strain i.These individuals harbor complete immunity to future infections with strain i and include those currently infected, i.e. y i , those that have recovered but were previously infectious, and those that were exposed, but protected from becoming infectious due to partial cross-protective immunity.Finally, w i represents the proportion of the population which has been exposed to any strain j which has at least one allele in common with strain i (including strain i itself), i.e. j ∩ i = Ø.These individuals have at least partial immunity to strain i. N.b. these equations are nested such that any individual in y i is also in z i and any individual in z i is also in w i , and y i ≤ z i ≤ w i ∀ strain i.In traditional Susceptible-Infected (SI), Susceptible-Infected-Recovered (SIR), etc. single-strain mathematical frameworks: the y class is analogous to the I class, while w and z are composed of combinations of I and R classes.The susceptible population is not modelled explicitly in this framework.
Explicitly, these three equations (for a given strain i) are: As above, we denote strains as subscripts and, in the equation for w i , we sum over all strains j which share at least one allele with the focal strain i. β, σ, and μ are the infection, recovery, and death rates, respectively.γ (as mentioned above) is an indicator of the level of cross-reactive immunity gained by prior exposure to alleles in the target strain.Note that while we depict only one value per demographic parameter (i.e.all strains are functionally equivalent) for clarity of notation, these values could also be written to vary by strain (e.g. Immunity in this framework is non-waning: exposure to a strain yields consistent protection from future infection over the lifespan of the individual.Moreover, this protection is trichotomous: an individual can either have no protection from a given strain (it has not seen any of the alleles before), complete protection (it has seen this exact combination of alleles before), or a set point in-between according to the parameter γ (it has seen at least one, but not all alleles before).Put another way, we do not distinguish between loci, assuming that sharing an allele at one locus is functionally identical to sharing an allele at any other locus, or indeed all other loci except one.

Extensions to more than one population
Following (18), we model movement between populations using a dispersal matrix where A is the weighted adjacency matrix containing elements A kl indicating the proportion of population k (row) moving to population l (column) per unit time and E is a diagonal matrix representing emigration, where each entry E kk = n k=1 A kl where n is the number of populations.
Thus, the whole system can be depicted by a set of three equations per strain i per population k: Where Δ T signifies the transpose of Δ, and each equation from Section 2.1 is now additionally indexed according to population and has an additional term to account for migration between populations.While in principle the elements of Δ can take any value [0, 1], signifying a (continuous) movement of between 0 and 100% of individuals, for simplicity we use a constant value δ for the strength of each movement, i.e. for each non-zero off-diagonal element of Δ. Sensitivity to this value is explored in the supporting information (Fig S2).
This framework can be applied to a metapopulation of arbitrary size and complexity, with the number of equations being linearly related to the number of populations.The dynamics of each population are governed by a set of three equations per pathogen strain, and these equations are interlinked within populations by partial, cross-reactive immunity, and between populations through a movement network.The total number of differential equations for any given system will be three times the number of strains multiplied by the number of populations in the metapopulation.

Simulation Procedure
All simulations were carried out in Julia version 1.3.0(19), with graphics produced using the ggplot2 package (20) in R version 3.6.1 (21).For simplicity of presentation, we fix the values of all variables other than β (the infection rate) and Δ (the network of movement information) to be identical for all populations in the metapopulation.To demonstrate the variety of dynamics obtainable in this modeling framework, we vary where E kk = -Δ kk , as noted above, signifies the total outgoing movement from the population of interest.We add a ˜over R 0 to denote that this is an approximation of the true reproductive number, the precise form of which would additionally take into account the inflow of infectious individuals from other populations.We additionally vary Δ according to the number and interconnectedness of the populations.For the figures of the main text, we utilize a strain structure of two loci, each with two alleles.Sensitivity to these parameter choices is explored in the supporting information (Fig S3).

Populations with identical parametrizations
To assess the effect of migration on population dynamics, we first consider the simplest case of a set of populations sharing the same disease parametrization: β such that R0 = 2, σ = 8, μ = 0.1, and γ = 0.66.We use a movement network described by a chain of populations, i.e.A→B→C→D or , where δ = 0.05, and ask how the dynamics of populations further down the chain (i.e.B, C, D) differ from those of the origin population (i.e.A), recalling that, without migration, all populations would have identical dynamics.

Populations with varying parametrizations
We next consider the case where parameters differ between connected populations, we restrict our consideration to a system of two populations, identical in all respects other than the parameter β, which is set to either induce a steady state of coexisting strains (β such that R0 = 5 in population A) or cyclically coexisting strains (β such that R0 = 2 in population B).We then display three potential patterns of connection: no migration (left column), A→B (middle column), and B→A (right column).Specifically, we set respectively, again with δ = 0.05, σ = 8, μ = 0.1, and γ = 0.66.
To address the case of multiple origin populations feeding into a single destination population, we consider a system of three populations: A→C←B, or , where populations A and C have β corresponding to an R0 = 5, but population B has β corresponding to an R0 = 2; all other parameters as above.

Larger network structure
Finally, we characterize the role of global network structure through considering the impact of degree distribution on a few summary statistics of system-wide disease burden: the mean proportion of infectious individuals (area under the currently infectious (i.e.y) curve), the mean level of strainspecific immunity (average z value), and the mean time between epidemic peaks (i.e. between local maxima in y) over the course of the final 33% of the simulation.We omit the initial period of the simulation to reduce the impact of transient dynamics.
We perform 100 simulations for each of five generic network ensembles each with 25 populations and a connectedness of approximately 0.15.Specifically, we examine Erdős-Rényi (links randomly assigned between populations), stochastic block (a metapopulation consisting of two groups of populations which have high migration within the group, but low migration to populations in the other group), tree-like (where there are many chains of populations and no potential for cycles), Barabasi-Albert (a scale-free network in which there tends to be a few populations with very many links, and many populations with few links), and Watts-Strogatz (a small-world network structure which is produced by partially re-wiring a spatially connected grid of populations) network structures.To generate these networks, we utilize functions from the tidygraph R package (22), except in the case of the tree and Watts-Strogatz configuration for which we use custom algorithms.Note that we use a parameter of attachment of 4 for the Barabasi-Albert random graphs.This allows for comparable connectences to the other random graphs as well as distinguishing these randomizations from trees (as would result from a default parameter of attachment of 1).In all cases, each migration strength is set to a constant δ = 0.01, only the pattern of connections varies.Each population is assigned a random β value corresponding to a R0 between [1,6].These results are qualitatively similar if instead every population is assigned the same value of β.

Results
In the following sections, we provide figures to demonstrate the effect of metapopulation structure on disease dynamics.In these figures, we plot a time series for each of three subsets of the population: those currently infected with a particular strain of the pathogen, those having (complete) specific immunity against the focal strain, and those who have at least partial cross-reactive immunity to the focal strain, due to past exposure to a similar strain (see Methods).Populations differ in their β (and thus R 0 ) value.This can be considered, for example, as differences in population density, which affects the probability of disease transmission.We only depict one representative strain in each plot for visual clarity and parametrize the model such that all strains are functionally equivalent (i.e. they all have the same transmission and recovery rates within any given population).

Cyclical dynamics are dampened along chains in the metapopulation network
We found that even when all populations share the same parametrizations and initial conditions, that populations further along network chains have reduced proportions of currently infectious individuals and dampened oscillatory dynamics compared to those they would exhibit in isolation (Fig 1).This is due to the movement of (partially) immune individuals between the populations, increasing the proportion of individuals with specific and cross-reactively immunity in populations further along the chain.While infectious individuals move at an equal rate, the proportion of the population that is currently infectious at any given time is much smaller than the proportion with immunity.

Dynamics propagate through metapopulation networks
We found that in the case of a simple chain of populations, the dynamics of destination populations can be overridden by the dynamics of origin populations (Fig 2).Interestingly, this is true both of cyclical dynamics overruling stable dynamics and vice versa, though the required amount of migration differs according to the origin and destination dynamics (see supporting information The issue of dynamics propagation gets more complicated when there are multiple, varying origin populations for a given destination population.We found that there is a hierarchy of dynamics in their propagation through the network: when there are origin populations with both cyclical and steady state dynamics, the destination population inherits the cyclical dynamics (Fig 3), albeit dampened from what they would have been without migration from a steady state population.
This asymetrical inheritance is robust to imbalance in the relative contributions of the origins.Put another way, if just one of many origin populations (or a small proportion of the total movement) has cyclical dynamics, the destination population will also have cyclical dynamics.

Degree distribution affects pathogen prevalence and immunity
These simple patterns in the effects of origin population dynamics on those in the destination population have clear implications when pieced together into larger network structures.For instance, the propagation of immune individuals through the metapopulation suggests that populations further "up the chain" will tend to have higher on-average disease incidence and also greater variability.
The inheritance of dynamical regimes combined with a hierarchy of dynamics in that inheritence suggests that chaos and cycles should be more common, especially in populations further "down the chain."That is, except in cases where the ultimate origin populations are all disposed toward steady states, in which case the stabilizing effect could overrule downstream local parametrizations, leading to an overall stable system.
In Fig 4, we report the effect of various network structures on three summary statistics of pathogen prevalence (and levels of immunity) using five common network ensembles.Depending on the system being explored, empirical network structures might have elements in common with one or more of these ensembles, for instance, many social networks are considered to be "small-world" in structure like Watts-Strogatz random graphs, while ecological networks are often commented on for their formation of "modules" or clusters of more densely interacting species as in stochastic 13 Posted on 28 Oct 2022 -The copyright holder is the author/funder.All rights reserved.No reuse without permission.-https://doi.org/10.22541/au.156026839.96630781/v2-This a preprint and has not been peer reviewed.Data may be preliminary.
block random graphs.Networks were parametrized to have approximately equal connectance and size in order to reduce uninformative variation (see Section 2.3.3).This is because metapopulation size and connectance have known effects on pathogen persistence, independent of further network structure (23; 24; 25).
We found that the network configurations with higher variation in indegree (i.e. the number of other populations each population receives migration from) distributions (supporting information Fig S4 ), such as those found in the tree and Barabasi-Albert networks, tend to have higher levels of infection over time, despite similar levels of immunity as the other three network types (Fig 4).
We saw similar patterns to those in infectious individuals when looking at time between epidemic peaks across network types.While fewer of the populations ended up with cyclical dynamics in the Barabasi-Abert graphs, the mean period of the cycles tended to be slightly higher and have higher variance, but this was not robust to alternative parametrizations (supporting information Fig S5 ).
Figure 4: The effect of network structure on pathogen prevalence and levels of immunity through time.In the top row, we depict a representative network from each of the five ensembles.The second row shows the distributions of each of three response variables for prevalence of the pathogen and specific immunity over the course of the simulation, with the units for the horizontal axes given by the panel headings.We depict one point for each randomized network structure and box-plots indicating the median and inter-quartile range of each network-type's distribution.Network generating algorithms were tuned to produce networks of the same size and approximate connectance and model parameters were either the same for all populations and across simulations (σ = 8, μ = 0.1, δ = 0.01, and γ = 0.66) or randomized for each population in each simulation (initial densities of infectious and immune individuals [0, 1] and β value corresponding to a R0 within [1,6] for each population).

Discussion
Both metapopulation (26) and strain (16) structure have long been known to be important to disease dynamics and are increasingly being recognized as ubiquitous.Yet the combination of these two areas of theory has been underexplored.We show here that this lacuna can have real consequences for our understanding of disease dynamics in empirical systems.
In probing the relationship between origin and destination dynamics in simple metapopulations, we have demonstrated several patterns that expand our understanding of disease dynamics in these systems.By directly incorporating a movement network into our model framework, we have constructed a very general approach that lends itself to arbitrarily large and complex systems.This is noteworthy, as more and more natural systems are being thought of in terms of networks of interacting components (e.g.separate species in ecological communities (27) or host individuals exchanging parasites ( 28)).By adjusting the scale of our metapopulation, we can ask and answer different questions about the forces influencing disease dynamics.
We found that the dynamics of prevalence and immunity among migrationally connected populations are not independent, and that even very small rates of population movement can have profound effects on a population's disease dynamics: from reducing pathogen prevalence to changing the dynamical regime of destination populations entirely.Our findings regarding the reduction in cycle amplitude (Section 3.1) echo results in dispersal networks in ecology, where population dynamics were dampened following the introduction of migration (29).
Contrary to prior focus in the literature on the role of migrating infectious individuals (30; 31; 26), we found that the migration of immune individuals can be equally (or even more) important.This is noteworthy, as the few previous studies relating multi-strain diseases and metapopulation structure only allow pathogen transmission between populations, not the movement of individuals explicitly (3; 32)-an approach that is more mathematically tractable, but omits the potentially influential transmission of immune individuals.
Finally, we show that larger network structure also has a part to play in disease dynamics, resulting in significant differences in pathogen prevalence across network types.Our results are in agreement with previous results suggesting increased epidemic size in scale-free network structures (such as those found in Barabasi-Albert random graphs) when the spreading rate is sufficiently slow (33; 34) due to the high-degree nodes serving as "super-spreaders" (35; 25).Along these lines, there has been some previous research indicating that node degree (the number of other populations a given population is connected to) is directly related to pathogen prevalence in that focal population ((36), but see ( 37)), however a complete investigation into network structure at the node-level is beyond the scope of this work.A comprehensive investigation of the role of more complex network structures in disease dynamics, however, remains a topic for further investigation.
In this work, we have utilized a relatively simple model for disease dynamics in an effort to maximize interpretability.Such simplification inevitably comes with a cost, and several of our assumptions can be critiqued as unrealistic.Perhaps foremost is the assumption of continuous movement.While continuous movement might be appropriate for very large populations with frequent, relatively small migrations between them, when any of these three components is not present, we would expect deviation from these predictions.Future work could explore the importance of discrete movement regimes on these patterns.Likewise, in this work we omit strain mutation and recombination ( 17) (yet the latter is included in the original framework of ( 16)).The generation of novel strains is likely important to the global persistence of diseases in humans (38) and animals (39).Finally, in representing movement by adding a proportion of the origin population to the destination, we introduce an assumption that the two populations in question are of approximately the same size.
In a metapopulation where populations vary widely in size, the proportion leaving one population would not correspond to the proportion entering another.
This work should not be seen as an attempt at comprehensive categorization of the role of metapopulation structure on the dynamics of multi-strain diseases, but rather as an initial step in exploring the complex interplay between the population structure of hosts and strain structure of pathogens.Our results suggest there may be simple rules underlying this relationship, at least for a wide range of parameter values, but it remains to be seen where networks based on empirical data fall in these parameter regimes, as well as how such systems might deviate from theoretical expectations.
We would like to thank José Lourenço for helpful discussions regarding the mathematical framework of this work.This work was supported by the USDA National Institute of Food and Agriculture, Animal Health project #MINV-62-057.

S1 Supporting Information
Figure S1: Considering all the dynamics of all four strains from population A in Fig 1 .Note that lines are colored according to strain rather than population.Strains can be divided into two discordant sets of non-overlapping alleles: {1, 1} and {2, 2}, and {1, 2} and {2, 1}.Each strain of a discordant set behaves identically due to identical parametrization and no interaction between strains that do not share at least one allele, but discordant sets interact with one another due to partial cross-reactive immunity.Thus, when one set is abundant, the other is rare and vice versa.
We highlight the maximum value of each discordant set's cycle with a vertical line in order to facilitate comparisons between strains and sets.
Figure S2: The effect of variable migration rate on transference of dynamical regime from origin to destination in a simple metapopulation of two populations linked by unidirectional movement.
In the top panel are cases where the origin population has cyclical dynamics and the destination has steady state dynamics when in isolation, while in the bottom panel the opposite is true.As the movement rate (δ· horizontal axis) increases, there is a phase-transition at which point the destination population's dynamics (indicated by the vertical axis) switch to match those of the origin.We fit a binomial spline to highlight this transition point.We see that even with very small rates of migration, a stable population can be converted to a cyclical one (top panel).Yet, it is more difficult to convert a cyclical population to one with steady state dynamics (bottom panel).Three parametrizations are recorded here (color; see legend), with additional parameter μ = 0.15 being the same for both populations.Finally, note that the two values of R0 listed in the legend correspond to the two populations, with the larger value corresponding to the steady state population and the smaller value the cyclical population (when in isolation) .While the observed differences in total infected are robust, note that here the mean time between epidemic peaks is approximately equal across randomizations.

Fig S2 ) .
Fig S2).This migration can also allow for strain coexistence even in populations where the local parameters would suggest extinction of one or more strains.

Figure 1 :
Figure1: Connecting multiple populations with the same parameters results in reduced pathogen prevalence and dampened cycles in populations further down the chain.Here, populations are connected such that A → B → C → D. Each column indicates a population, while each row is one of the three population classes laid out above and in the Methods, i.e. those currently infectious with the given strain, those with (complete) specific immunity to the given strain, and those with partial cross-protective immunity to the given strain.The mean level of immunity (both specific (middle row) and cross-reactive (bottom row)) increases in each sequential population, while the mean level of currently infectious individuals (top row) decreases.All populations have parameters σ = 8, μ = 0.1, δ = 0.05, γ = 0.66 and a β chosen to make R0 = 2 for all populations.We use a two-loci, two-allele strain structure, but show only one strain for clarity (but see supporting information FigS1).

Figure 2 :
Figure 2: Destination populations tend to inherit origin population dynamics when linking populations with different model parametrizations.As in Fig 1, rows correspond to population classes, but here, columns indicate network structure.While in isolation (left column), population A has cyclical dynamics and population B has steady-state dynamics, when the two populations are linked by migration, the destination population inherits the dynamics of the origin population (center and right columns).This is true regardless of the direction of the movement (depending on the level of migration; see supporting information Fig S2).Populations A and B have parameters σ = 8, μ = 0.1, δ = 0.05, and γ = 0.66 in common and β chosen to reflect R0 = 2 and 5, respectively.We use a two-loci, two-allele strain structure, but show only one strain for clarity.

Figure 3 :
Figure 3: When multiple origin populations differ in their dynamics, the destination population inherits cycles over steady states.As in Fig 1, rows indicate population classes, and columns the component populations.Here, we have populations A and B feeding into population C at the same rate of δ = 0.05.Populations A and C are parametrized to produce steady state dynamics in the absence of migration, with σ = 8, μ = 0.1, γ = 0.66, and β corresponding to a R0 = 2. Population B shows cyclical dynamics with β corresponding to a R0 = 5 and all other parameters the same.Note that, even though the parameters of population C would lead to a steady state in the absence of migration, we see cyclical dynamics being inherited from population B. We use a two-loci, two-allele strain structure, but show only one strain for clarity.

Figure S3 :
Figure S3: Effect of parametrization on dynamic regime for a population in isolation.Here, columns indicate values of γ and rows indicate values of σ.Depending on the combination of β, σ, μ, and γ, a population can exhibit a range of dynamics including steady states for all strains (pink), cyclical or chaotic dynamics for all strains (orange), or partial extinction of some strains (blue).Green cells indicate a numerical failure in integration.All simulations here utilize a two-loci, two-allele strain structure.See (16) for a similar figure for alternative parameterizations.

Figure S4 :
Figure S4: Summary statistics for the degree distributions of each randomized network used for Fig 4 in the main text.Networks were constructed to have the same size and approximate connectance, but with the network structure (which populations are connected to which other popuations)otherwise generated according to one of five algorithms: Erdős-Rényi, Barabasi-Albert, and Watts-Strogatz, stochastic block, and tree (see Section 2.3.3 of the main text).Some algorithms allowed perfect matching of connectance (Erdős-Rényi, Barabasi-Albert, and Watts-Strogatz), while others necessitated some minor variation (stochastic block and tree).

Figure S5 :
Figure S5: Similar to the lower row of Fig 4,but with γ and σ equal to 0.55 and 32, respectively.All other parameters are equal to or set randomly as in Fig 4. All other parameters are equal to or set randomly as in Fig 4.While the observed differences in total infected are robust, note that here the mean time between epidemic peaks is approximately equal across randomizations.