key: cord-0726777-k6fk13ry
authors: Cota, Wesley; Soriano-Panos, David; Arenas, Alex; Ferreira, Silvio C.; G'omez-Gardenes, Jes'us
title: Infectious disease dynamics in metapopulations with heterogeneous transmission and recurrent mobility
date: 2021-02-21
journal: New Journal of Physics
DOI: 10.1088/1367-2630/ac0c99
sha: 429d8cc92bbf7e2c11119971779dc19e5a9d6415
doc_id: 726777
cord_uid: k6fk13ry

Human mobility, contact patterns, and their interplay are key aspects of our social behavior that shape the spread of infectious diseases across different regions. In the light of new evidence and data sets about these two elements, epidemic models should be refined to incorporate both the heterogeneity of human contacts and the complexity of mobility patterns. Here, we propose a theoretical framework that allows accommodating these two aspects in the form of a set of Markovian equations. We validate these equations with extensive mechanistic simulations and derive analytically the epidemic threshold. The expression of this critical value allows us to evaluate its dependence on the specific demographic distribution, the structure of mobility flows, and the heterogeneity of contact patterns, thus shedding light on the microscopic mechanisms responsible for the epidemic detriment driven by recurrent mobility patterns reported in the literature.

The proliferation and accessibility of large data sets describing the essential aspects of human behavior is being crucial to reveal the influence that our social habits have on the development of epidemics as well as providing useful insights to design nonpharmaceutical containment strategies. Human mobility is one of the aspects of our social behavior determining the form and speed of the transmission of infectious diseases.

In this sense, the recent availability of data about the mobility patterns of individuals at different levels [1, 2, 3] , from global to urban, demands to revisit epidemic models, in particular those studying the geographical spread of pathogens leveraging the mobility of hosts [4] .

Data-driven models are developed to improve the spatio-temporal accuracy of predictions of real epidemic outbreaks by using a large amount of real data as inputs [5, 6, 7, 8, 9, 10, 11] . However, agent-based and mechanistic models based on large-scale stochastic Monte Carlo simulations have as a counterpart the impossibility of performing analytical treatments that shed light on the role played by the different aspects of our sociability in the transmission of communicable diseases. To fill the gap between accurate epidemic forecasting systems and mathematical models, theoretical frameworks should be refined in order to be able to incorporate as much social data as possible.

The most usual way to incorporate mobility patterns into epidemic models is the use of metapopulations. In this case, individuals are considered to live in a set of subpopulations (or patches) whereas flows of individuals happen among these patches. Within this framework, the spread of diseases is characterized by local reactions inside each patch [12, 13, 14, 15, 16] that mimic the interactions between individuals giving rise to the transmission of the pathogen. This reaction process within each patch interplays with the global diffusion of agents that captures the mobility patterns at work.

The first metapopulation frameworks were built by considering assumptions that simplify their mathematical analysis while limiting their direct application in real situations. However, with the advent of the XXI century and the massive use of online platforms, real data capturing individual flows between different geographical areas were incorporated into metapopulation frameworks [17, 18, 19, 20] in an attempt of increasing their accuracy while preserving the ability to perform analytical predictions. Still, the first models in this line assumed simple mobility patterns such as random diffusion [21, 22] or continuous models of commuting flows [23, 24, 25] , that allowed analytical studies about the influence of mobility on the epidemic threshold [17] . The next step in the search for more reliable and accurate metapopulation models was to get rid of the simplifying assumptions about human diffusion and find ways to take into account aspects such as the recurrent nature [26, 27, 28, 29, 30] and high order memory of human displacements [31] , the coexistence of different mobility modes [32] , or the correlation between the time-scales associated to human mobility and that of infection dynamics [33] . These models, apart from yielding important insights about the role that human behavior has on the unfolding of epidemic states, have turned to be useful tools to reproduce the real prevalence distribution of endemic diseases [34] and the advance of real epidemic outbreaks [35, 36] , thus showing a versatile and hybrid facet as mathematical yet informative models.

The former refinements have focused on the way real human mobility patterns are incorporated into metapopulation frameworks, but continue using simple mixing rules for the interaction of individuals within each patch. These simplifying hypotheses include well mixing assumptions and explore scenarios where the number of contacts inside each patch is homogeneous and usually determined by some demographic aspects such as the density of the patch or its age distribution. However, human contact patterns are known to be highly heterogeneous and this attribute plays a central role in the transmission of some communicable diseases [37] . In fact, the analysis of the propagation of recent coronavirus such as SARS-CoV-1 [38, 39, 40] , MERS-CoV [41, 42] and SARS-CoV-2 [43, 44, 45, 46, 47, 48] , reveals that a small proportion of cases were responsible for a large fraction of the infections. This empirical evidence supports the existence of superspreading events [49] , an attribute of transmission chains that cannot be captured by models in which the contacts of individuals, and hence their infectiousness, are assumed to be homogeneous.

There have been some attempts in the literature to account for the impact of individual diversity in metapopulation modeling [50, 35] . However, they usually rely on the stratification of the population into different age-groups [51] , which are assumed to be homogeneous, and the introduction of mixing matrices governing the interactions among them. Therefore, a general formalism able to accommodate heterogeneous subpopulations with any arbitrary degree distribution is still missing in the literature. In this paper, we aim at filling this gap and including the heterogeneity of social contact patterns in the body of a metapopulation model, in particular that presented in reference [30] and used in subsequent works [32, 34, 35] .

The most important result found in these works was the detrimental effect of human daily recurrent mobility for the emergence of epidemic outbreaks. Nonetheless, the mean-field assumption included within each subpopulation in these formalisms precludes getting any microscopic explanation about the mechanism triggering this phenomenon. The model presented here is therefore a step forward towards a metapopulation formalism that includes concomitantly the demographic distribution of real populations, the recurrent nature of human displacements, and the heterogeneity of social contacts and sheds light on the unexpected phenomena arising from their interplay. In fact, the most important finding in this new framework is that the detrimental effect of human daily recurrent mobility is recovered despite the fact that the number of interactions does not depend on the number of agents that meet inside each patch. Thus, individual interactions appear here as an intensive parameter, rather than an extensive one as in reference [30] , shedding light on the microscopic roots of the epidemic detriment phenomenon.

Let us start the construction of the metapopulation framework by describing the interaction rules that govern the mixing of individuals across and within patches. We consider a metapopulation network with Ω patches, each one of population n i Figure 1 . Example of a metapopulation with two patches, both having the same average connectivity k = 5. The first is a heterogeneous patch with resident individuals of connectivity 1 or 20, and the second is a homogeneous patch in which all residents have the same connectivity 5.

(i = 1, . . . , Ω), thus accumulating a total of N = i n i individuals. Each individual is associated with a single residence (one of the patches) and can travel to another patch according to some mobility rules. The flow of individuals from a patch i to another j is described by a directed and weighted network of patches, in which the weight W ij is the number of individuals from i that commute to j daily. The matrix W ij is also called origin-destination (OD) matrix and allows us to define the probability that, when an individual living in i decides to move, she or he goes to patch j as

where Ω l=1 W il = s i is the total number of trips observed from patch i. According to the framework presented in reference [30] , mobility and interactions are iterated in consecutive rounds of a process that involves three stages: Mobility, Interaction, and Return (MIR). Namely, first the agents with residence in a patch i decide to move with probability p (or they stay in i with probability 1 − p). If they move, their destination j is chosen with probability R ij , given by equation (1) . Once all the agents in each patch have been assigned to their new locations (either their residence or a new destination chosen according to the matrix R) the interaction on the assigned patch takes place with the rest of agents in the same subpopulation. Finally, once the interaction stage has finished, agents are placed in the original population, i.e., they come back to their corresponding residence. Now we propose a modification to consider heterogeneous contacts inside each patch. In reference [30] , all individuals inside a patch interact with all others with the same probability thus following a homogeneous mixing hypothesis. Here we propose a model in which each individual in a patch has a different social degree or connectivity k as shown in figure 1. In this way, each patch i has n [k] i individuals with connectivity k, so that the population of patch i can be written as:

where P i (k) is the probability that a randomly chosen individual living inside i has a connectivity k:

In the following, we assume that individuals with social connectivity k will preserve this value when traveling to another patch, i.e., we assume that sociability is an intrinsic individual attribute that does not depend on their location. This later hypothesis captures the biological and behavioural aspect of hosts that can turn them into superspreaders, i.e., individuals that are highly efficient in transmitting the disease due to a high viral shedding [52] or because they have a high contact rate due to a pronounced social behavior. However, other causes that are inherently related to the location, such as the existence of high-risk scenarios related to work or leisure, are not captured by the former assumption.

Under the former hypothesis about the invariance of the connectivity k under mobility and assuming that those individuals with connectivity k move with probability p k , we can calculate the effective population of a patch i,ñ i , after the movement stage has been performed, as the sum of the effective number of agents with connectivity k:

In the latter equation,ñ

[k] i is calculated considering the number of individuals with connectivity k that travel from any patch j to i:

where

Another quantity that can be evaluated is the effective connectivity distribution of a patch,P i (k), defined as the probability of finding an individual of connectivity k in patch i after the mobility stage. This probability is given by:

From the effective connectivity distribution of a patch i we can measure the effective moments as:

The coupling of interaction and mobility patterns of agents produces, for a given set of mobility probabilities {p k }, a variation of the main structural attributes of the patches, as shown by the expressions of the effective population, equations (4)- (5) , and the effective connectivity distribution, equation (7) . These variations occur once the mobility step is performed and become crucial when the spreading process (the interaction step of the MIR model) enters into play.

Here the interaction stage is incorporated as a Susceptible-Infected-Susceptible (SIS) spreading dynamics. To this aim, we denote the number of infected individuals residing in i that have connectivity k as I [k] i , implying that the total number of infected residents in i is

i . Thus, the probability that an agent with residence in patch i and connectivity k is infected is given by:

The probabilities {ρ

i } (with i = 1, . . . , Ω and k = 1, . . . , k max ) constitute our dynamical variables. From these variables we can compute the fraction of infected individuals with residence in patch i:

or the fraction of infected individuals in the whole metapopulation:

To derive the corresponding Markovian evolution equations of the probabilities {ρ [k] i } corresponding to the SIS dynamics we make use of the so-called heterogeneous mean-field theory (HMF) in the annealed regime [53] . Thus, after the movement stage, each susceptible agent with connectivity k that is placed in patch j connects randomly with k individuals in the same patch and, for each infected contact, the susceptible agent will become infected and infectious with probabilityλ. In addition, those infected agents at time t will recover and become susceptible again with probabilityμ. Following these simple rules, the equations for the time evolution of the probabilities {ρ

where Π i (t) is the probability that a healthy individual with connectivity k and residence in patch i becomes infected at time t:

where π i (t) is the probability that an individual of connectivity k placed in patch i becomes infected at time t and reads:

In the former expression,P i (k |k) is the probability that an agent with connectivity k placed in patch i is connected with another agent with k placed in the same patch. In addition,ρ [k] i is the effective fraction of infected individuals with connectivity k placed in patch i:ρ

where the denominator is given by (5) and the numerator is the number of infected individuals that are in patch i. In the following we will consider that the contact networks created at each interaction step are completely uncorrelated. This way, the probabilityP i (k |k) can be written in terms of the effective connectivity distribution of patch i as:

which is the probability of selecting an edge from an individual with connectivity k placed in patch i, independent of k.

The derived Markovian equations are general for a set of Ω patches, their population n i , degree distribution P i (k), and OD matrix elements W ij , (i, j = 1, . . . , Ω). We now study the impact of heterogeneous distributions of individual contacts by using synthetic metapopulations to validate these equations by comparing the results obtained by the iteration of equations (12)- (14) with the results of mechanistic Monte Carlo (MC) simulations in which we keep track of the dynamics of each agent.

Although the formalism presented can accommodate any arbitrary mobility network and set of connectivity distributions, we restrict our analysis, as in reference [30] , to synthetic star-like metapopulation networks. Our choice is rooted in their versatility for, despite being simplistic structures, star-like metapopulations exhibit a wide variety of regimes caused by the non-uniform distribution of the population across patches and the asymmetry in th mobility patterns connecting them. This kind of synthetic metapopulation, shown in figure 2 , is composed by a central patch (the hub) connected to κ patches (the leaves). The hub h has a population of n h individuals, while each leaf l has a fraction α ∈ [0, 1] of the hub population, n l = αn h . The mobility towards leaves of individuals with residence in the hub is uniform, given by: Figure 2 . Example of a star-like metapopulation network with κ + 1 patches. In this example, the leaves and the hub have the same number of individuals, n l = αn h , with α = 1, while the hub is a heterogeneous patch with resident individuals of connectivity 1 with probability η, or k max = 20 with complementary probability 1 − η, and each leaf is a homogeneous patch with residents of same connectivity k l = β k h , with β = 1 and k h = 5. The flow from hub to a leaf happens with probability R hl = κ −1 , from leaves to hub with R lh = δ, and between adjacent leaves with R l,l+1 = 1 − δ, in counterclockwise direction.

while the mobility of those residents in the leaves is controlled by a parameter δ. This way, a resident in a leave l that decides to move will go to the hub with probability δ,

or move to the next (counterclockwise direction) leave with probability

Note that the choice of the direction of movements among leaves is not relevant as long as it is uniform across all the leaves, for they are statistically equivalent. Up to this point, the design of the metapopulation is identical to that presented in reference [30] , being characterized by two parameters α and δ. However, the synthetic metapopulations used here get rid of the assumption of homogeneous (all-to-all) contact patterns in the patches. To this aim, and keeping the symmetry of the original star-like metapopulations, we consider that the residents of the central patch (the hub) have a contact distribution P h (k) that is different from that of the residents in the leaves, P l (k). A particular case of this setting used along the manuscript is to consider that the connectivity distribution of the individuals belonging to the hub is bimodal:

i.e., agents in the hub have connectivity 1 with probability η and connectivity k max with probability (1 − η). This way, the n-th moment of the hub's connectivity distribution is:

In their turn, those individuals belonging to leaves have the same number of contacts ( k l ):

Note that the values of η and k max are correlated if we impose the additional constraint that the hub has an average connectivity k h fixed. In this case, given a value k max , the value of η that allows it is given by:

In this simple configuration, the heterogeneous nature of the contacts is twofold. From a microscopic point of view, the bimodal distribution existing inside the central node induces local heterogeneities in the contacts made by residents there, which are controlled by parameters η and k max . In its turn, another global connectivity heterogeneity emerges driven by the asymmetry existing between the connectivity of residents of the hub and the leaves. In particular, we will assume throughout the manuscript that k l = β k h , with β ∈ [0, 1]. According to this formulation, the star-like metapopulation shown in figure 2 has k h = 5, k max = 20, and α = β = 1.

To check the validity of the Markovian equations, we define a MC algorithm for the stochastic simulation of the SIS model on top of a metapopulation with heterogeneous contact patterns. As in the case of Markovian equations, equations (12)- (14) , the proposed process is also a discrete-time dynamics. At each time step t, each individual is tested to move with probability p k (being k the number of contacts assigned to this individual). If accepted, it moves to a patch j with probability R ij . Then, each susceptible individual with connectivity k chooses randomly k individuals in the patch they currently occupy and are infected with probabilityλ if the contacted individual is infectious. Once all the potential infections events have been simulated, healing happens with probabilityμ for each infected individual at time t − 1. In this sense, we perform a synchronous update of the state of the entire metapopulation.

First, a fraction ρ ini of the population is randomly infected as the initial condition and the simulation procedure in a give time step t can be summarized as follows:

(i) For each patch i, each individual with connectivity k resident in i is tested to move with probability p k . If she or he moves, a patch j is chosen proportionally to R ij .

(ii) Each susceptible individual with connectivity k selects k contacts at random in patch i. For each attempt, it can be infected with probability:

or remains susceptible with complementary probability. These attempts stop when the individual becomes infected and reproduce the annealed regime proposed in section 2, since all edges are available for each individual in the same time step.

(iii) Each individual with infected state at time step t − 1 heals in time step t with probabilityμ.

(iv) Finally, all individuals return to their residences and time step t + 1 starts in (i).

To avoid the absorbing state, we infect a small fraction ρ pump = 2 × 10 −4 of individuals at random when this state is reached [54, 55] . This keeps the dynamics always active and the equilibrium state is defined after comparing averages over sequential time windows of size T = 100, and accepting if the absolute difference is smaller than ρ cvg = 10 −6 .

The comparisons between MC and Markovian equations are performed in star-like metapopulations with κ = 10 and α = 1, i.e., in which all patches (hubs and leaves) contain the same number of individuals (n l = n h = 10 4 individuals per patch), to focus on the effect of contact heterogeneity. Furthermore, for the same reason, we focus on the case that mobility is independent of the connectivity of individuals, p k = p ∀ k.

First we neglect local heterogeneities and consider that contact heterogeneity only happens between patches. In mathematical terms, this assumption implies that the population of the hub has an homogeneous contact distribution (η = 0) although its mean connectivity k h = k max is different from that of the leaves k l = β k h , with β = 1. In particular, in figure 3 we plot the mean epidemic prevalence ρ * in the equilibrium state as a function of the infection probabilityλ scaled by the epidemic threshold in the case of null mobilityλ 0 ≡λ c (p = 0). To derive the latter quantity, we realize that the absence of flows among the patches precludes the interaction among the residents in different areas, so the epidemic threshold corresponds to the well-known expression provided by HMF equations [53] for the most vulnerable patch. Therefore,

We consider that k h = 100 while leaves have k l = 10 (β = 0.1) and explore two different mobility patterns. In particular, in (a) we set δ = 0.1 so that most of the residents of leaves move circularly, i.e., passing from one leave to another and avoiding the hub. In this case, the so-called epidemic detriment by mobility shows up so that the epidemic state is delayed as the mobility p increases, with the exception of very large values of p. However, note that, at variance with reference [30] , here both the hubs and the leaves are equally populated; we will explore the roots of this detriment below. Second, in panel (b), we set δ = 0.9 so that the situation is the opposite and the residents of leaves tend to visit the hub. In this case, the epidemic detriment is also evident although this behavior is restricted to values p < 0.5, while for p > 0.5 the increase of mobility produces a progressive decrease of the epidemic threshold. In both cases, the agreement with MC simulations is almost perfect. Next we analyze a star-like metapopulation that generalizes the contact heterogeneity of the first one. In this case the hub is very heterogeneous, containing a power-law distribution, P h (k) ∼ k −γ h withγ h = 2.3, while leaves have also a powerlaw distribution P l (k) ∼ k −γ l withγ l = 3.5, both with k ∈ [3, 100], the hub being the most heterogeneous one. The cases explored in figure 4 are again (a) δ = 0.1 and (b) δ = 0.9, showing similar qualitative behaviors with the mobility, namely the emergence of epidemic detriment, to those found in figure 3 . Quantitatively, it is worth stressing that the existence of strong local heterogeneities within both hub and leaves in absence of mobility will lead to an activation described by the HMF theory, in which the epidemic prevalence approaches zero close to the epidemic threshold as ρ ∼ (λ −λ c )β whereβ > 1 if the degree exponent is smaller than 4 [56] , and valid for large population sizes (thermodynamic limit). The convexity of the prevalence curve approaching the transition in the finite-size population of the investigated patches is reminiscent of this behavior. Again, the agreement with MC is good, except around the epidemic threshold due to difficulties in avoiding the absorbing state. 

Figures 3 and 4 reveal that the epidemic detriment emerges even when dealing with uniformly distributed populations, contrarily with reference [30] , in which increasing mobility in homogeneous populations favors epidemic spreading by reducing the epidemic threshold,λ c , here defined as as the minimum infectivity per contact,λ, such that an epidemic state can be stable. Therefore, the emergence of epidemic detriment here should be rooted in the interplay among contact heterogeneities and human mobility. In this section, we aim at deriving an analytical expression of the epidemic threshold,λ c for general configurations, to shed light on the mechanisms giving rise to the behavior shown above. Let us assume that the dynamics has reached its steady state, so that ρ 

. Under this assumption, equation (12) reads:

Furthermore, forλ values close to the epidemic threshold, the fraction of infected individuals is negligible, which means that ρ * k) . This fact allows us to linearize the equations characterizing the steady state of the dynamics by neglecting all the terms O(¯ 2 ). In particular, the probability that an individual with connectivity k and placed in i contracts the disease, π * i [k] , can be approximated by

where we have used O(ρ) = O(¯ ) as shown by equation (15) . In particular, plugging (15)-(16) into the last expression leads to:

where

is the effective number of edges in patch i. Note that i Q i = k j kP j (k)n j is the total number of edges in the system, a conserved quantity. After introducing (29) and some algebra, equation (27) transforms into:

Finally, if we introduce these values into equation (26) and retain only linear terms in¯ , we arrive to the following expression

that defines an eigenvalue problem. According to its definition, the epidemic threshold is thus given by:λ

The elements of matrixM given by (32) represent four types of interactions in the metapopulation. Namely, the elementM jk ik represents the probability that a resident of patch i with connectivity k is in contact with another individual of patch j and connectivity k . The first term accounts for interactions of residents of the patch, that do not move. In second term, an individual of i stays and interacts with a traveler from patch j in patch i, that arrived with probability p k R ji . A similar event happens in the third term, in which an individual of i travels to patch j and interact there with a resident of j with probability p k R ij . Finally, in the forth term, both individuals of patches i and j travel to a patch l, arriving there with probability p k p k R il R jl . In computational terms, each row or column identifies individuals from one degree class living inside a patch. Therefore, the dimension of the matrix corresponds with the sum of the different degree classes observed within each patch.

Equation (34) computes the exact expression of the epidemic threshold in presence of heterogeneous contact patterns. However, its computation involves solving the spectrum of a matrix whose dimension is determined by the number of connectivity classes and patches in the metapopulation. In particular, in presence of highly heterogeneous populations with fine spatial resolution, this problem can be computationally very hard due to a large number of elements of the critical matrix. For this reason, in what follows, we assume that mobility is independent of the connectivity so that p k = p which will considerably reduce the complexity of the problem as proved below.

Before going ahead, it is convenient to make the transformation¯ ik → k ik in equation (33) . Note that this represents a similarity transformation which does not alter the spectrum of the matrix. After doing such transformation, equation (33) 

where the elements of the new matrix M read as

If p k = p, equation (35) becomes independent of k, which allows a dimensionality reduction of the matrix. In particular, equation (35) reads:

and the elements of the reduced matrix M are given by:

where the effective number of edges Q i is now expressed as

Once matrix M is constructed the epidemic threshold is computed as

To test the accuracy of the former expression for the epidemic threshold, we compare its value computed according to equation (40) with the heat map of the steady state of the dynamics obtained from the iteration of equations (12)- (14) . Figure 5 the theoretical prediction of the epidemic threshold by equation (40) is very accurate and captures the dependence of the epidemic threshold on the mobility p. This threshold increases while promoting mobility until it reaches a maximum at p = p * since the infection is gradually reduced in the hub as p increases, and the activation is then triggered in the leaves since hub's residents spend longer times there. For the sake of completeness, in Appendix B, we analyze the case p = 0 for equation (40) retrieving, as expected, the expression for the epidemic threshold provided by HMF equations on contact networks. Moreover, to quantify the effects of promoting mobility among disconnected patches, we perform a perturbative approach to the latter threshold which holds for small p values in Appendix C. Interestingly, at variance with the perturbative analysis carried out for (non-structured) well mixed metapopulations in reference [30] , here the linear correction of the epidemic threshold strongly depends on the topological properties of the metapopulation.

In what follows, to shed light on the nature of the epidemic detriment, we aim at quantifying the impact of the different components of the formalism, namely the underlying metapopulation structure and the contact heterogeneities existing among its population, on the relative magnitudeλ c (p * )/λ 0 . To simplify this analysis, we will focus on the case of mobility independent of k, p k = p, and consider the configuration defined in section 3.1, in which the hub has individuals with connectivity either 1 or k max , with fixed average connectivity k h , and the ones of the leaves have the same connectivity k l = β k h . For the sake of clarity, let us also express k 2 l = γ k 2 h . Note that in this configuration the values of η and k max are correlated by equation (23), while γ is also correlated with β and k max via

First, we fix α = β = 1, so that n l = n h and k h = k l , to study the effects of varying either the local heterogeneity existing in the hub by tuning k max or the flows from leaves to the hub with δ in figure 5 (b). Fixing k max = 50 and changing δ, it becomes clear that the increase of δ leads to a decrease of p * as a consequence of the higher mixing among individuals from the central node and the leaves, but does not change the relative magnitudeλ c (p * )/λ 0 .

The former beneficial effect is rooted in the homogenization of the connectivity distribution driven by the mixing among individuals from the hub and the leaves. Interestingly, the position of the peak p * remains unaltered when keeping δ constant. Moreover, for small values of p, the behavior does not depend on the local heterogeneities of the patches, as shown by a perturbative analysis in Appendix C. Quantitatively, it becomes clear that increasing the degree heterogeneity in the central node boosts the beneficial effect of the mobility, since the homogenization effect gains more relevance due to the higher vulnerability of the central node. Mathematically, the invariance of p * , when introducing local contact heterogeneities without varying the mobility patterns, implies that the spatial distribution of cases close to the epidemic threshold -controlled by the components of the eigenvector of matrix Mis ruled by the structure of the underlying mobility network. We also observe that the value of the epidemic threshold at the peak p * is independent of the mobility network but is instead determined by the local heterogeneities, the difference in mixing of the subpopulations.

Finally, we extend our analysis to cover populations distributed heterogeneously across the metapopulation. In particular, we are interested in determining how the population asymmetry α and the local connectivity heterogeneity η shape the relative magnitude of the peak of the epidemic threshold. To this aim, we represent λ c (α, β, γ; p * )/λ 0 (β, γ) in figure 6 , for n l = αn h , k l = β k h , and k 2 l = γ k 2 h , in which γ is given by equation (41) for the constraints imposed in section 3.1. We can observe that, as in figure 5(b) , increasing the local heterogeneity of the hub (lowering γ) increases the beneficial effect of the population mixing, as shown in figure 6(a) . Interestingly, if we fix γ and study the dependence ofλ c (α, β, γ; p * )/λ 0 (β, γ) with α and β, as shown in figure 6 (b), we observe that the detriment effect becomes stronger for larger values of β since k max increases so to keep γ constant. In the opposite direction, when reducing the population of the periphery nodes, i.e., decreasing α, agents in the leaves are not able to substantially modify the connectivity distribution of residents in the hub, thus hindering the detriment effect in all investigated cases.

Driven by the advance of data mining techniques in mobility and social patterns [1, 57, 58] , epidemic models are continuously refined to bridge the gap existing between their theoretical predictions and the outcomes of real epidemic scenarios. In particular, within the very diverse realm of epidemic models, the proliferation of data sets capturing human movements across fine spatial scales have prompted the evolution of metapopulation frameworks, which constitute the usual approach to study the interplay between human mobility and disease spreading. In this sense, the first theoretical frameworks assuming the population to move as random walkers across synthetic metapopulations [17] have given rise to models incorporating the recurrent nature of human mobility [30, 59, 60, 28] , the socio-economic facets of human movements [32, 61] or high-order mobility patterns [31] .

While most of the advances previously described have been focused on capturing the mobility flows more accurately, less attention has been paid to improve the contact patterns within each subpopulation. With few exceptions, such as the model recently proposed in [62] incorporating the time varying nature of social contacts, human interactions are usually modeled using well-mixing hypothesis that do not capture the heterogeneous nature of human interactions and the role that this social heterogeneity has on the so-called super-spreading events.

In this work, we tackle this challenge and adapt the metapopulation model presented in reference [30] to account for the heterogeneity in the number of contacts made by individuals. We describe a complete set of Markovian equations for a discrete-time Susceptible-Infected-Susceptible dynamics on subpopulations with recurrent mobility patterns. These equations characterize the spatio-temporal evolution of the number of infected individuals across the system and show a good agreement with extensive agent-based simulations results. Computationally, iterating the equations of our formalism is orders of magnitude faster than performing the simulations because the latter should account for each microscopic stochastic process occurring in the population at each time step. Apart from the computational advantages, our formalism allows for deriving analytical results on the interplay between epidemics, mobility, and the structure of contacts within the metapopulation. Specifically, the linearization of these equations yields an accurate expression for the epidemic threshold, which is a crucial indicator for the design of interventions aimed at mitigating emerging outbreaks.

Our most important finding here is the emergence of the epidemic detriment when enhancing mobility, despite the fact that the individuals preserve their number of interactions independently of the visited locations. This result cannot be explained following the macroscopic arguments proposed in reference [30] and shed light on the microscopic nature of the epidemic detriment phenomenon. In particular, it becomes clear that this phenomenon is inherent to the variation of the contact structure of the population driven by redistribution of its individuals. Specifically, close to the epidemic threshold, the outbreak is mainly sustained by super-spreaders and the ties existing among them, which are weakened due to the homogenization of the underlying connectivity distributions caused by human mobility. Interestingly, the epidemic detriment observed in critical regimes is reversed in the super-critical regimes, where mobility increases epidemic prevalence, for it increases the average number of potentially infectious contacts made by scarcely connected individuals.

The formalism here presented constitutes a step forward to account for the interplay between contact and flow structures and thus present several limitations. First of all, we assume that the number of interactions of each individual is constant and depends on the features of her residence patch, regardless of the place to which they move. Although this assumption can be interpreted as the preservation of the sociability of individuals, it prevents us from accounting for super-spreading events [46] associated to events or particular gatherings in which social connectivity is punctually amplified. In addition, as remarked in the former paragraph, the results here obtained rely on assuming uncorrelated connectivity distributions within each patch. In this context, the effect of degree-degree correlations inside the patches deserves to be investigated; for example, one could expect the epidemic detriment to lose relevance in assortative populations, where ties connecting super-spreaders are strengthened and less likely to be influenced by the mobility. Finally, although we have explored the physics of the interplay between contact heterogeneity and recurrent mobility with simple synthetic metapopulation networks, the model represents a general framework that can accommodate any arbitrary set of degree distributions within a population and any mobility network structure. In this sense, when data is available, the model can be investigated using a data-driven approach in the sense that one can easily include real data of demographics, mobility, and contact patterns to describe more realistic situations. threshold will be given by the first subpopulation in the active state, if its population is not so small compared to other patches. Indeed, the usual epidemic threshold known in the HMF theory is obtained,λ c =μ min

(2.1) Therefore, in the static case the epidemic threshold of the metapopulation corresponds to the individual epidemic threshold of the most vulnerable patch.

We proceed by making a perturbative analysis of the eigenvalues of the matrix M up to first order on p to complement the discussions of the main text. First, it is convenient to rewrite equation (38) to split the terms with different order in p:

Since Q i is also a function of p, we must perform a Taylor expansion around p = 0, knowing that Q i | p=0 = n i k i . The first derivative of Q i is dQ i dp p=0 = j k j (R ji − δ ij ) n j .

Let us define r i ≡ j (−R ji + δ ij ) n j k j , so that d dp

Next, keeping only terms up to order 1, we have

Substituting the last expression in (3.1) we get, after some algebra,

where From the static case, we know that there are Ω unperturbed eigenvalues Λ (0) i = k 2 i / k i , for p = 0, with normalized eigenvectors i = { j } and j = δ ij ; see equation (2.1) . Assuming that the eigenvalues are not degenerate, the new eigenvalues will be given by [64] Λ i ≈ Λ 

Interestingly, unlike the original MIR model, the first order correction depends on the underlying topology. To check the accuracy of this correction, we represent in figure C1 the leading eigenvalues of the matrix M along with the linear correction provided by the perturbative analysis, finding a remarkable agreement in the low mobility regime p 1.

Generalitat de Catalunya (grant No. 2017SGR-896 and 2020PANDE00098), and Universitat Rovira i Virgili (grant No. 2019PFR-URVB2-41). A.A. also acknowledges support from Generalitat de Catalunya ICREA Academia, and the James S. McDonnell Foundation (grant 220020325). W.C. and S.C.F. acknowledge financial support from CAPES (grant No

acknowledge financial support from MINECO (projects FIS2015-71582-C2 and FIS2017-87519-P), from the Departamento de Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo

Stochasticity and heterogeneity in the transmission dynamics of SARS-CoV-2

In this case, we have to evaluate seven different terms:• M hh : contact of two individuals residing in the hub;• M lh : contact of one resident from a leaf with another from the hub;• M hl : contact of one resident from the hub with another from a leaf;• M ll : contact of two individuals residing in the same leaf;• M l,l+1 : contact of one resident from a leaf with another from its adjacent leaf;• M l,l−1 : contact of one resident from the adjacent leaf with one from the other leaf;• M ln : contact of two residents from different and not adjacent leaves;The mobility matrix elements R ij are expressed in eqs. (17) to (19) . Applying these expressions in (38) , we have Disease dynamics with heterogeneous transmission and recurrent mobility patterns 20Again, by evaluating equation (37), we have, for the hub, in which the factor 3 in the last term is since there are κ − 3 other leafs not directly connected to a single leaf (R ln = 0). The statistical equivalence of the leaves allows us to recast the computation of the epidemic threshold in a eigenvalue problem of a 2 × 2 matrixThe leading eigenvalue will be given by Λ max = Tr M + (Tr M ) 2 − 4 det M 2 , that was solved using SymPy [63] to get the results shown in the main text.

To check the consistency of these equations, let us consider the static case in which all individuals stay in their patches and do not move: p k = 0 ∀k. So, equation (38) becomeswhere Q i | p k =0 = n i k i , that after being used in (37) results inμ i =λ k 2 i k i i . This case consists of isolated subpopulations in an annealed regime in which the epidemic