Human mobility networks and persistence of rapidly mutating pathogens

Rapidly mutating pathogens may be able to persist in the population and reach an endemic equilibrium by escaping hosts’ acquired immunity. For such diseases, multiple biological, environmental and population-level mechanisms determine the dynamics of the outbreak, including pathogen's epidemiological traits (e.g. transmissibility, infectious period and duration of immunity), seasonality, interaction with other circulating strains and hosts’ mixing and spatial fragmentation. Here, we study a susceptible-infected-recovered-susceptible model on a metapopulation where individuals are distributed in sub-populations connected via a network of mobility flows. Through extensive numerical simulations, we explore the phase space of pathogen's persistence and map the dynamical regimes of the pathogen following emergence. Our results show that spatial fragmentation and mobility play a key role in the persistence of the disease whose maximum is reached at intermediate mobility values. We describe the occurrence of different phenomena including local extinction and emergence of epidemic waves, and assess the conditions for large-scale spreading. Findings are highlighted in reference to previous studies and to real scenarios. Our work uncovers the crucial role of hosts’ mobility on the ecological dynamics of rapidly mutating pathogens, opening the path for further studies on disease ecology in the presence of a complex and heterogeneous environment.

SM, 0000-0001-6202-3302; CP, 0000-0002-4051-1716 Rapidly mutating pathogens may be able to persist in the population and reach an endemic equilibrium by escaping hosts' acquired immunity. For such diseases, multiple biological, environmental and population-level mechanisms determine the dynamics of the outbreak, including pathogen's epidemiological traits (e.g. transmissibility, infectious period and duration of immunity), seasonality, interaction with other circulating strains and hosts' mixing and spatial fragmentation. Here, we study a susceptible-infected-recovered-susceptible model on a metapopulation where individuals are distributed in sub-populations connected via a network of mobility flows. Through extensive numerical simulations, we explore the phase space of pathogen's persistence and map the dynamical regimes of the pathogen following emergence. Our results show that spatial fragmentation and mobility play a key role in the persistence of the disease whose maximum is reached at intermediate mobility values. We describe the occurrence of different phenomena including local extinction and emergence of epidemic waves, and assess the conditions for large-scale spreading. Findings are highlighted in reference to previous studies and to real scenarios. Our work uncovers the crucial role of hosts' mobility on the ecological dynamics of rapidly mutating pathogens, opening the path for further studies on disease ecology in the presence of a complex and heterogeneous environment.

Introduction
Many pathogens are able to persist in a host population by repeatedly reinfecting individuals who lose the immunity acquired during infection [1]. This mechanism is particularly important for the case of rapidly mutating pathogens when mutation is associated with antigenic drift-e.g. for influenza virus [2] and respiratory syncytial virus [3]. The probability of pathogen persistence is then determined by the interplay between the transmissibility and infection time scales (i.e. duration of infection and immunity period), along with ecological factors like seasonality, interaction with other circulating pathogens, mixing patterns among individuals, their spatial distribution and mobility.
Waning of immunity and partial immunization have been accounted for by several modelling studies [4][5][6][7][8][9] addressing detailed immunological mechanisms [10], interaction with other strains [7,11], human contact structure [6], environmental factors [12], seasonality and temporal variation in transmission [13][14][15]. The fate of an emerging pathogen in terms of persistence/extinction was also studied, both as a theoretical problem [8,9] and through applications to the emergence of a novel influenza strain [2,16]. However, the vast majority of studies on waning of immunity and persistence have disregarded spatial structure so far.
Recently, it has been shown that the geographical distribution of the host population and the rate of travelling among distant locations determine the time scale of human mixing, hence affecting the probability of having large-scale epidemic spread [17][18][19][20] and the incidence profile [21][22][23][24]. For the case of childhood diseases (e.g. measles), previous studies showed that these dynamical effects may impact disease persistence and alter the critical community size [25][26][27][28][29][30][31], i.e. the population threshold above which a disease is able to persist [25]. Results reported so far, however, suggest that the persistence probability is determined by the interplay between the different spatial ingredients-e.g. distance and frequency of travels, environmental influences on transmission-in a non-trivial way. Indeed, different contrasting hypotheses have been suggested, with persistence maximized for intermediate levels of spatial coupling [28,30], for maximum levels of coupling [31] or either of the two depending of the level of spatial heterogeneity in transmission [29].
The spatial structure of human population could have important and non-trivial effects on the spreading of diseases also when waning of immunity is at work, affecting the possible dynamical regimes following an outbreak. In this paper, we address this problem by means of a stochastic discrete spatially explicit metapopulation model. We reconstruct the dynamics of an emerging infection through extensive numerical simulations. We characterize the role played by the spatial fragmentation, the coupling, and the structure of the mobility network through which individuals travel from one sub-population to another. We focus on the regime of parameters that best describes influenza-like infections. Our metapopulation framework realistically reproduces the statistical properties of the human population, i.e. strong heterogeneities and correlations in the population distribution and mobility. We perform a detailed analysis of the mechanisms shaping the probability of pathogen persistence and the spatial dynamics following emergence in the short and long term. Our approach offers the theoretical and computational framework to further explore the effects of other features not taken into account here (e.g. multiple strains, seasonality and others) in order to gain a full understanding of pathogens' persistence and be able to perform realistic data-driven simulations.

Material and methods
We use a stochastic individual-based metapopulation model [32][33][34], where mobility and infection of individuals are explicitly accounted for as discrete-time processes. Individuals mix homogeneously within the local communities (also called sub-populations, patches or nodes of the metapopulation system), whereas at the global level the coupling between these sub-populations is introduced by individuals that travel along mobility connections (figure 1). The model is thus represented in terms of a network of links, i.e. mobility flows, connecting different nodes representing local populations.
Massive data and extensive research studies have uncovered the statistical properties of the distribution of human populations and their associated mobility, pointing out their network structure, characterized by large variabilities and non-trivial correlations [35][36][37][38]. In particular, the population of urban areas and the number of connections between locations (considering e.g. air-transportation or commuting) span several orders of magnitude, and more-populated centres are also often the hubs of the mobility network. We account for these properties by assuming a power law distribution for the patches' connectivity and a linear relationship between connectivity and population. Specifically, we considered a network of V patches, each of them (node i of the network) characterized by a population of size N i , and a degree k i that represents the number of sub-populations through which node i is linked via mobility connections. The nodes' degree distribution is assumed to be of the form P(k) ∼ k −γ with . Schematic representation of the stochastic metapopulation model. Panel (a) represents the disease dynamics that follows a SIRS scheme, where individuals change states with given transition probabilities (Susceptible (S), Infected (I), Recovered (R) and back to S when the acquired immunity is lost). Panel (b) illustrates the structure of the metapopulation network, which consists of sub-populations that are connected by individuals travelling from one to another through a mobility network.
no degree-degree correlation present; a random network with these characteristics was generated by means of the uncorrelated configuration model [39] with γ = 3.0. Mobility fluxes are assigned to each link by assuming that, in each sub-population, individuals may travel at each time step to neighbouring sub-populations with a probability p. Departing individuals from node i choose at random one of the available k i links, so that the probability of travelling along a connection is given by p/k i . This process yields a population distribution that at equilibrium is proportional to the number of connections, namely, N i = k i N / k ; where N and k are the average population per node and average degree of the metapopulation network, respectively [20]. The infection dynamics is implemented through the Susceptible-Infected-Recovered-Susceptible (SIRS) compartmental scheme [4]. The evolution of the dynamics results from the following transition rules iterated at each time step, corresponding to 1 day: (i) Susceptible individuals within each sub-population i may be infected with probability where β represents the transmission probability, and N i and I i denote, respectively, the total population and the number of infected individuals in node i; (ii) Infected individuals enter with probability μ the recovered compartment, that represents a state in which individuals are immune to the disease; here immunity is assumed to be complete; (iii) Eventually, recovered individuals lose their immunity and get back to the susceptible state.
This occurs with probability λ, that is the inverse of the average immunity period, L.
Infection transition rules combined with travelling can be summarized, for each sub-population, by the following discrete master equations: .
where υ i represents the set of neighbouring patches of i. The described process disregards vital dynamics. This choice simplifies the dynamical characterization of the model and it is motivated by our interest in diseases with durations of both infection and immunity periods much shorter than the average life duration (e.g. influenza). Starting with a fully susceptible population, the disease is seeded in a randomly chosen node by infecting 0.5% of its population, in such a way that stochastic runs differ for both the initial seed and the stochastic realization of the dynamics. The disease may then propagate inside the local population and spread to neighbouring populations due to host movements. Progression of the disease and host movements are stochastic events that are numerically modelled through multinomial processes. For each scenario considered, 10 3 independent epidemics are simulated starting from different seeded sub-populations. We consider influenza-like infections spreading in space as a prototypical example for our study. We fix the duration of infection μ −1 = 2.5 days and the basic reproductive number R 0 = β/μ = 2.0. The population is assumed to be composed of 10 8 individuals divided in V patches. We explored scenarios characterized by: (i) different spatial fragmentation of the host population, considering V = 10 2 , 10 3 , 10 4 ; (ii) a varying coupling between patches expressed by host mobility, considering p ∈ [10 −6 , 1]; (iii) different topological structures of mobility connections, comparing the heterogeneous network described by a power-law degree distribution to a homogeneous Erdős-Rényi graph [40] described by a Poisson degree distribution with the same average value; (iv) a varying immunity period, considering L ∈ [1, 800] days.
For each set of parameters, we numerically monitor the global prevalence (i.e. the total fraction of infectious individuals) and the number of infected patches D(t) in time. We also estimate the probability that the pathogen persists in the population at the global level, P gl , by computing the fraction of stochastic runs for which the epidemic reaches the endemic state, once a given set of parameter values is considered. The same quantity can be defined at the local level to measure the probability of local persistence P loc . This can be computed considering a closed population, or a population integrated into a metapopulation framework, thus open to incoming and outgoing fluxes of individuals. At the local level, we also define an outbreak as a single epidemic (infecting at least 0.1% of the population) that continuously persists in the population. Repeated outbreaks instead are considered as distinct waves following local extinction. These behaviours are characterized for the different scenarios explored, in order to provide a coherent picture of the mechanisms for virus persistence and endemic equilibrium in the spatially structured population.

Results and discussion
In the following, we present first the results regarding the persistence of the pathogen in the population and the mechanisms behind it, followed by a detailed characterization of the spatial dynamics across the scenarios explored.

Interplay between waning of immunity and mobility
The probability of persistence is found to be very high and equal to 1 for short immunity periods; once immunity becomes longer, a reduction of the probability of persistence is observed (figure 2a). This behaviour is common to all host mobility rates explored. When the waning of immunity is very fast, there is a high rate of renewal of susceptible individuals that can be newly infected by the pathogen, thus ensuring its survival in the population. A slower rate of renewal makes persistence less likely and the specific behaviour of the decrease in persistence probability depends on the coupling between patches. Reduction may occur through a drop to zero persistence (e.g. p = 10 −6 , 10 −5 , etc. and p = 10 −2 , 10 −1 , etc.), or a drop to an intermediate positive value of persistence probability that is maintained for a large range of immunity period, before reaching P gl = 0 (e.g. p = 10 −4 , 2 × 10 −4 , etc.).
Both panels show results obtained with V = 10 4 , N = 10 4 and power law degree distribution.
In the low mobility regime, persistence at the global level is mainly determined by the probability that the pathogen would persist locally in the seeded population. If p is low enough so that no infectious individual travels out of the infected patch during the first wave of the epidemic, pathogen survival following the first wave depends on the length of the immunity period. A long enough immunity period prevents the replenishment of healthy individuals in the short term, thus leading to the extinction of the epidemic in the seeded patch, with no chances for global spread. If the availability of new susceptible hosts is instead provided fast enough by a short immunity period, the epidemic can survive the first wave in the seeded population, reaching a local endemic equilibrium. This, in its turn, can slowly seed the epidemic in neighbouring patches even when small mobility rates are considered, thanks to the local equilibrium condition. This same process is then occurring in each newly seeded patch, so that propagation unfolds at the global level and a high probability of persistence is observed.
The condition distinguishing between an epidemic trapped locally in a patch and one that spreads at the global level has been quantified in previous studies in terms of the global invasion threshold [17,20,[41][42][43][44][45]. Its analytical expression, provided in [20], provides an upper bound for the low mobility regime Once informed with the demographic and topological features of our spatially structured population, this condition is found to provide a good approximation to the limit of the low mobility regime (red arrow in figure 2b). The spatial structure of the population becomes less important once the system is found in the high mobility regime. In this case, the probability of moving from one patch to another is so high as to largely increasing the opportunity of mixing of the populations belonging to different patches. The system behaves effectively as a single population of the size of the full metapopulation (i.e. in our case 10 8 individuals) affected by a SIRS dynamics. The conditions for global persistence are therefore those described before for the seeded population in the low mobility regime. Differently from before, here local persistence is equivalent to global persistence, due to the large degree of mixing across space. The values of the immunity period at which a given persistence probability is observed (e.g. L 300 for P gl = 5%) differ from the case of low mobility regime (L 100), because of the difference in the population size, i.e. a single patch in the low mobility regime and effectively the full metapopulation size in the high mobility regime.
The intermediate mobility regime displays a higher persistence probability for a given duration of the immunity period, compared with low and high mobility regimes (figure 2b). The contour line at a fixed P gl (e.g. P gl = 5% or 95%) exhibits a maximum for intermediate values of the travelling probability p that is reached for high values of L. In addition, a smoother variation of the persistence probability is observed by varying L for a given p, compared with the faster transitions recovered in the low and high mobility regimes. The intermediate degree of mixing provided by intermediate mobility seems to enhance stochastic effects responsible for local extinction, thus having an important impact on the persistence at the global level.
Maximum persistence for intermediate level of coupling was also reported in the context of childhood diseases, where SIR models with continuous external introduction of the infection were considered [28,30]. No consensus was, however, found in this context, with other studies reporting different findings. For instance, Hagenaars et al. found an increase in the persistence with increasing levels of the coupling among sub-populations [31]. Rozhnova et al. confirmed this finding when homogeneous transmission across patches is considered, and additionally showed that heterogeneity in transmission among patches selects dynamical regimes with different spatial effects [29]. Different results obtained from different models show how the detailed mechanism of susceptible replenishment, together with infectious transmission and recovery, and spatial mixing, have an impact on the resulting persistence. Prompted by this, we assess in the following the role played by the spatial fragmentation of the population and by the mobility network on our model where the replenishment of individuals is obtained through immunity waning.

Role of spatial fragmentation and mobility network
Spatial fragmentation is explored by subdividing the total population of N = 10 8 hosts into V = 10 2 , 10 3 , 10 4 patches. Global persistence in the three scenarios shows the same qualitative behaviour described in the previous section (figure 3a). The separation between different mobility regimes, however, depends on the number of patches considered. More precisely, it is a demographic effect due to a varying average population size per patch. Referring to the transition from low mobility to intermediate mobility regime, we have indeed that the estimate expressed in equation (3.1) depends on N . Increasing the number of patches from V = 10 2 to V = 10 3 , 10 4 , the average size of the population in a given patch decreases by one and two orders of magnitude, respectively, thus inducing a corresponding increase in the mobility value of the transition from low to intermediate regimes. As a result, the maximum of the global persistence is obtained for correspondingly increasing values of the travelling probability.
It is interesting to note that the behaviour observed in the high mobility regime is the same across the three scenarios. This confirms that the resulting dynamics when p is high enough corresponds effectively to the local persistence dynamics of a single population of 10 8 individuals, as discussed in the previous subsection.
If the topological structure connecting the patches via mobility links is altered and homogenized, we find that higher values of mobility need to be reached for the transition from low to intermediate mobility to occur (figure 3b). As expected, heterogeneous mobility patterns (as the ones adopted in figure 2) favour hosts' mobility resulting in a larger global persistence for lower travelling probabilities. This is consistent with what is predicted by the global invasion of equation (  same. Mixing is high enough so that the structured population behaves effectively as a single population regardless of the topology through which such mixing occurs.

Spatial dynamics
In this section, we go more in depth into the understanding of the spatial dynamics following the emergence of a pathogen in the fully susceptible and spatially structured population. Starting from the global persistence diagram reported in figure 2, we characterize the epidemic diffusion in space considering a single value of reference for p in each of the three mobility regimes (p = 3 × 10 −5 in the low, p = 6 × 10 −4 in the intermediate, and p = 10 −2 in the high mobility regime), associated with different values of the immunity period L (figure 4).

Low mobility coupling
We first analyse the low mobility regime, by characterizing the spatial epidemic dynamics obtained for p = 3 × 10 −5 and L = 30 leading to P gl = 1.0 (figure 4g), and p = 3 × 10 −5 and L = 80 leading to P gl = 0.50 (figure 4d). As discussed previously, global persistence is determined by local persistence of the disease at the patch level. The very high persistence observed for p = 3 × 10 −5 and L = 30 is ensured by an endemic equilibrium that is reached in all patches after a given time, as signalled by 100% of patches experiencing an outbreak (i.e. D(t)/V = 1 after a given t). Even with a low probability of moving from one patch to another, the initially seeded population is able to transfer infected hosts to neighbouring populations in the long time limit thanks to the endemic equilibrium reached. In their turn, newly infected patches can  further propagate to the infection reaching all patches of the system, thus achieving a high probability of persistence of the pathogen at the global level, fuelled by each endemic dynamics. The same travelling probability leads to lower P gl once the immunity period increases, from L = 30 to L = 80 (figure 4d). Here, a similar mechanism of local persistence and spatial diffusion ensured by the endemic equilibrium is at play, but it does not occur in all patches. At the end of our simulations, we find indeed that only a fraction of patches is experiencing an outbreak (D/V 80%). The different dynamics occurring across patches depends on the patch population size. For a fixed value of the immunity period L, there exists a minimum population N min needed in order to achieve a certain probability of local persistence (figure 5) [8,25,46,47]. The population size of the patch of our system varies approximately in the range [5 × 10 3 , 2.7 × 10 5 ], due to the heterogeneous topology connecting the patches. If L = 30, we find that all population values are above the minimum population needed to achieve local persistence and mobility can act as spatial spreader to reach global persistence. If L = 80, we find that a portion of the patches do not have enough population to achieve local persistence. In particular, to obtain a local persistence probability higher than 95% the minimum population size is 5 × 10 4 . The epidemic is able to spread at the global scale; however, smaller peripheral patches do not meet the condition for local endemic equilibrium. This results in a slower spatial spread, compared with the case of L = 30 (as displayed by the time axis in figure 4d,g), and an epidemic that is not able to endemically reach all sub-populations (D/V < 1).
The dynamics for the L = 80 case reflects a spatial hierarchy similar to the one studied in the context of childhood diseases-notably, for measles epidemics in the UK [23,24,48]. At equilibrium, highly  Minimum population size N min needed to achieve a local persistence probability P loc = 5% and P loc = 95% as function of the duration of immunity L for the case of an isolated population and for populations of the metapopulation system. The shaded area denotes the possible population values in our system. connected and populated patches are steadily infected, while the peripheral ones experience repeated outbreaks on the rare occasions in which infectious travellers arrive from the hubs and find enough susceptible individuals to induce a local outbreak. To characterize this hierarchy, we analyse how local persistence changes across classes of patches and quantify the number of outbreaks each class experiences throughout the duration of the simulated epidemic. Classes of patches are defined by their connectivity, i.e. we focus on nodes having the same number k of neighbouring sub-populations. Owing to the correlation between degree of a patch and its population size, these connectivity classes correspond effectively to population classes. Increasing degree leads therefore to higher local persistence, because of the relation to the population size (figure 6a for the low mobility regime). P loc reaches 1 in high-degree nodes for small enough values of the immunity period (i.e. up to L = 100 days among the cases shown in the figure 6a). In addition, we find that in these conditions hubs consistently experience only one outbreak (figure 6b), corresponding to an epidemic that persists in the population after the first seeding event. Low-degree nodes instead experience repeated outbreaks seeded from neighbouring nodes interspersed with local extinctions. For higher values of L (e.g. L = 150), the maximum local persistence achieved even in highly connected nodes is smaller than 1 and global persistence is approximately P gl = 0.03. Our results show that when disease persistence occurs, it is driven by repeated outbreaks taking place in different locations. Patches with intermediate degrees have the highest turnover of consecutive epidemics (about three outbreaks measured in the explored timeframe; figure 6b). Hubs present fewer outbreaks, as these last longer because of their larger population size and due to the continuous reseeding from neighbouring patches that lead to multiple consecutive waves. Peripheral nodes instead have less outbreaks because, being less reachable, they are not frequently reseeded after local extinctions by other nodes.

Intermediate mobility coupling
For intermediate values of the travelling probability p, the spatial dynamics displays multiple waves followed by an endemic behaviour, if the probability of global persistence is larger than 5% ( figure 4b,c,h). For smaller persistence (e.g. P gl = 0.01, figure 4a), the epidemic goes extinct after an initial wave. Interestingly, the amplitude of the initial wave is determined only by the mobility, regardless of the final persistence or extinction outcome. We explain this by noticing that the first wave refers to the process of spatial invasion from the initially seeded patch, that can be quantified by the global invasion threshold [20]: measuring the number of neighbouring patch that a seeded patch can infect. This is the analogue of the basic reproductive number at the sub-population level. Though this estimate was obtained for a SIR dynamics, we consider it here as we focus on the initial spatial invasion following the seeding, before the mechanism of immunity waning becomes important. The invasion condition R * > 1 indeed  figure 4a,b,c,h), we find indeed the same initial invasion, regardless of the duration of the immunity period. Following the initial wave, dampened oscillations in the number of infected patches are observed at the spatial level, for large enough persistence probability. At the local level, oscillations in the incidence are caused by the waning of immunity and have a frequency that decreases with L and the patch population size, as predicted by the SIRS theory [4]. These local oscillations then interact at the spatial level because of the mobility coupling, translating in multiple waves in the spatial propagation. Those may be the result of travelling waves departing from well-connected nodes [23,24,48,49]. Figure 6c indeed shows that if the pathogen is able to persist in the population, after the initial transient, the epidemic has more chances to persist in highly connected nodes. The more peripheral patches experience a sequence of sporadic outbreaks, similarly to the low mobility case (figure 6d). Differently from before, here the local persistence is enhanced by the spatial coupling provided by the intermediate mobility. The continuous reseeding by infectious individuals travelling between sub-populations, in a situation of low level of synchronization between local epidemics, fuel the circulation of the infection, hence allowing for a non-zero persistence even at high values of L, corresponding to the maximum discussed in §3.1. This rescue effect was previously introduced in the context of childhood diseases [26][27][28]. Here, we find that it is driven by the increase in local persistence in highly connected nodes as illustrated in figures 5 and 6. In figure 5, the comparison of this mobility regime with the low mobility one shows that the rescue effect lowers the population threshold needed to have persistence P gl ≥ 0.05 and P gl ≥ 0.95. Moreover, comparing the intermediate and low mobility regimes (blue curves in figure 6a,c) for the case L = 150, we see that while local persistence is not achieved when patches are fairly isolated, this is not the case for a large class of nodes as soon as mobility increases.

High mobility coupling
To characterize the unfolding of the epidemic in the high mobility regime, we consider two points of figure 4e having the same travelling probability p = 10 −2 , and different lengths of the immunity period, namely L = 150 and L = 250, corresponding to P gl = 1.0 and P gl = 0.60, respectively. The epidemic dynamics in these cases shows an oscillatory behaviour followed by an endemic equilibrium ( figure 4f,i), similar to what was previously discussed in the intermediate mobility regime. Here, frequencies of oscillations are, however, higher. This may be induced by an increase in synchronization among patches, enabled by the stronger level of spatial coupling [31]. Larger travelling fluxes lead, indeed, to almost all patches being infected during the first wave (e.g. D/V 1 for p = 10 −2 , figure 4i, compared with D/V 0.6 for p = 6 × 10 −4 , figure 4h, for the same value of L). Over time, these fluxes guarantee a large degree of seeding across patches along with increased mixing. For L 360 and in a situation of global persistence, local persistence is found to be much higher across degree classes compared with intermediate mobility-see e.g. curves in figure 6e versus figure 6c for the same value of L and consistently the values of N min in figure 5 corresponding to the p = 10 −2 case. Spatial heterogeneities are therefore progressively reduced, with a larger number of patches behaving similarly and epidemic differences across patches almost disappearing.
If persistence is enhanced locally for all patches, the resulting more synchronized system is, however, more prone to extinction of the epidemic for large enough immunity periods. This is the same type of phenomenon observed for measles in [26][27][28], i.e. a decrease in spatial heterogeneity reduces the effectiveness of the rescue effect by enhancing the level of synchronization. To test this hypothesis, we compared the local and importation dynamics occurring in the most connected patch in the intermediate and high mobility regimes. We considered a single stochastic run to highlight the critical role of stochastic effects within this mechanism (figure 7). Immunity is maintained equal in both scenarios, for comparison. When p = 6 × 10 −4 (intermediate mobility, figure 7a), importation of cases from neighbouring patches, desynchronized with respect to the local dynamics, allows the epidemic to survive the first wave and then quickly reach an endemic equilibrium, i.e. the rescue effect referred to before. If mobility is higher (p = 0.01, figure 7b), importations are synchronized with the local dynamics, i.e. a peak of importations is occurring during the peak of the local epidemic wave, and both profiles fade out at the same time. In this way, there is no importation from outside the patch that could help sustain or relaunch the outbreak in the patch, and the epidemic goes extinct.

Conclusion
As has been previously discussed in the literature, the epidemiological mechanism of waning of immunity alters the fate of an emerging infection and its chance to persist in the population. Here, we have characterized in depth this effect once a spatially structured population with varying coupling is considered. By using a stochastic individual-based metapopulation model that explicitly accounts for travelling fluxes among sub-populations, we have found that there are three distinct mobility regimes. Our results show that, when mobility is low, the persistence conditions for the disease in the metapopulation system are the same as if the patches were isolated. When mobility is very high the epidemic behaves as if the whole metapopulation were well mixed. Eventually, at intermediate levels of mobility, the coupling between sub-populations creates rescue effects due to non-synchronized epidemics ongoing in the patches that maximize disease persistence. These effects are able to reseed the epidemics locally once it became extinct.
Different persistence outcomes are associated with a variety of short-term dynamical behaviours. These are due to the interplay between local replenishment of susceptible individuals and the level of epidemic synchronization across patches. For a very rapid waning of immunity, persistence occurs regardless of the level of mobility coupling, with local endemic equilibria allowing for the epidemic to propagate from one patch to another even at very low travelling probabilities. For longer durations of the immunity period, the probability of having local persistence is shaped by the high heterogeneity in the population distribution (large patches may be above the critical size needed for local persistence, while small ones may not) and the rescue effect. These two mechanisms are important at intermediate levels of travelling coupling and are responsible for the peak of persistence observed in this regime.
Our findings are similar in many aspects to the ones obtained for the case of measles and, more in general, for an SIR model with vital dynamics and external introductions [23,24,27,28,30,31,48,[50][51][52]. This exemplifies that different mechanisms of susceptible replenishment can lead to common dynamical features. Distinct and apparently contradicting behaviours have been shown in previous studies [28][29][30][31], highlighting the important role that mobility can have in triggering these mechanisms. Here, we have provided a detailed description of the epidemic dynamics for the case of a rapidly spreading acute infection with waning of immunity, where the main statistical features of human space distribution and mobility are accounted for in a realistic way. We showed that the main result, i.e. that persistence is maximized at intermediate mobility couplings, is valid for different levels of spatial fragmentation and also for a homogeneous distribution of population and connectivity.
In order to better focus the analysis on the interplay between mobility and time scale of waning of immunity, we considered throughout the manuscript a single value of R 0 . We expect different values of R 0 yielding the same qualitative picture with an increase in the parameter determining a shift of the persistence peak toward smaller values of p as predicted by equation (3.1). We also disregarded ingredients known to impact infectious disease dynamics in many real cases. In particular, we assumed no seasonality in transmission and no multi-strain interference. These two ingredients are known to alter the spread of diseases like influenza and respiratory syncytial virus, thus they are usually included in realistic data-driven modelling of these infections [2,7,53]. The inclusion of seasonality in the studies of childhood diseases and its interplay with demography showed the emergence of annual and biannual cycles in the dynamics of the disease [23,24,51]. In the systematic exploration conducted here for the SIRS case, accounting for these ingredients would be complicated by the large increase in the complexity of the dynamics, the number of parameters and the variety of possible outcomes in terms of persistence, transient dynamics and equilibrium state [21,22,29,[54][55][56], thus hindering a clear theoretical understanding. The chosen approximations may, however, limit the applicability of our work directly to real-case situations, such as the case of the emerging influenza strains [2,[57][58][59].
Other aspects could alter the results, such as human response and intervention to the epidemice.g. vaccination, case isolation and travel restrictions. In particular, phenomena like travel bans or spontaneous reductions of travel often occur during the initial stage of an outbreak as a consequence of the public concern about the emerging infection. They were shown, however, to have limited impact on the early stage invasion [19], given that reductions of orders of magnitude needed to significantly affect spatial spread are unfeasible in reality [60,61]. For the same reason we argue that their effect would be minimal also on pathogen persistence, unless the system is in the intermediate mobility regime, where probability of persistence is more sensitive to the level of spatial coupling (figure 2). Vaccination, on the other hand, by boosting the level of immunity in the population increases the chance of pathogen extinction and may profoundly alter the spreading dynamics, affecting the frequency of epidemic waves [23]. For the case of rapidly mutating pathogens the quantification of vaccine efficacy and its impact on disease eradication is complicated by the fact that immunity caused by vaccination also decays in time when novel virus variants emerge, and vaccine match with the circulating pathogen is often compromised due to the difficulty in anticipating virus mutation in the vaccine design and production [58].
Starting from the complete characterization of the epidemic dynamics proposed here, we plan to address these points in the future in order to better understand the epidemic dynamics resulting from a complex interplay of different factors.