Infectious disease dynamics in metapopulations with heterogeneous transmission and recurrent mobility
Wesley Cota, David Soriano-Paños, Alex Arenas, Silvio C. Ferreira, Jesús Gómez-Gardeñes
IInfectious disease dynamics in metapopulations withheterogeneous transmission and recurrent mobility
Wesley Cota , , David Soriano-Pa˜nos , , A. Arenas , Silvio C.Ferreira , , Jes´us G´omez-Garde˜nes , , Departamento de Fisica, Universidade Federal de Vi¸cosa, 36570-900 Vi¸cosa, MinasGerais, Brazil GOTHAM Lab – Institute of Biocomputation and Physics of Complex Systems(BIFI), University of Zaragoza, E-50018 Zaragoza, Spain Department of Condensed Matter Physics, University of Zaragoza, E-50009Zaragoza, Spain Department of Computer Science and Mathematics, Universitat Rovira i Virgili,43007 Tarragona, Spain National Institute of Science and Technology for Complex Systems, 22290-180 Riode Janeiro, Brazil Center for Computational Social Science, Kobe University, Kobe 657-8501, Japan.February 2021
Abstract.
Human mobility, contact patterns, and their interplay are key aspects ofour 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 modelsshould be refined to incorporate both the heterogeneity of human contacts and thecomplexity of mobility patterns. Here, we propose a theoretical framework that allowsaccommodating these two aspects in the form of a set of Markovian equations. Wevalidate these equations with extensive mechanistic simulations and derive analyticallythe epidemic threshold. The expression of this critical value allows us to evaluate itsdependence on the specific demographic distribution, the structure of mobility flows,and the heterogeneity of contact patterns.
Keywords : Epidemic models, human mobility, metapopulations, complex networks
1. Introduction
The proliferation and accessibility of large data sets describing the essential aspectsof human behavior is being crucial to reveal the influence that our social habits haveon the development of epidemics as well as providing useful insights to design non-pharmaceutical containment strategies. Human mobility is one of the aspects of oursocial 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 individualsat different levels [1, 2, 3], from global to urban, demands to revisit epidemic models, in a r X i v : . [ phy s i c s . s o c - ph ] F e b isease dynamics with heterogeneous transmission and recurrent mobility patterns isease dynamics with heterogeneous transmission and recurrent mobility patterns
2. Metapopulation model
Let us start the construction of the metapopulation framework by describing theinteraction rules that govern the mixing of individuals across and within patches.We consider a metapopulation network with Ω patches, each one of population n i ( i = 1 , . . . , Ω), thus accumulating a total of N = (cid:80) i n i individuals. Each individual isassociated with a single residence (one of the patches) and can travel to another patchaccording 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 alsocalled origin-destination matrix and allows us to define the probability that, when anindividual living in i decides to move, she or he goes to patch j as R ij = W ijN (cid:88) l =1 W il , (1)where (cid:80) Nl =1 W il = s i is the total number of trips observed from patch i .According to the framework presented in [30], mobility and interactions are iteratedin consecutive rounds of a process that involves three stages: Mobility, Interaction, andReturn (MIR). Namely, first the agents with residence in a patch i decide to move withprobability p (or they stay in i with probability 1 − p ). If they move, their destination j is chosen with probability R ij , given by Eq. (1). Once all the agents in each patchhave been assigned to their new locations (either their residence or a new destinationchosen according to the matrix R ) the interaction on the assigned patch takes placewith the rest of agents in the same subpopulation. Finally, once the interaction stage isease dynamics with heterogeneous transmission and recurrent mobility patterns Figure 1.
Example of a metapopulation with two patches, both having the sameaverage connectivity. The first is a heterogeneous patch with resident individuals ofconnectivity 1 or 20, and the second is a homogeneous patch with all residents withthe same connectivity 5. In both patches, the average connectivity is (cid:104) k (cid:105) = 5. has finished, agents are placed in the original population, i.e. , they come back to theircorresponding residence.Now we propose a modification to consider heterogeneous contacts inside eachpatch. In Ref. [30], all individuals inside a patch interact with all others with thesame probability thus following a homogeneous mixing hypothesis. Here we propose amodel in which each individual in a patch has a different social degree or connectivity k as shown in figure 1. This way, each patch i has n [ k ] i individuals with connectivity k ,so that the population of patch i can be written as: n i = (cid:88) k n [ k ] i = (cid:88) k n i P i ( k ) , (2)where P i ( k ) is the probability that a randomly chosen individual living inside i has aconnectivity k : P i ( k ) = n [ k ] i n i . (3)In the following, we assume that individuals with social connectivity k will preservethis value when traveling to another patch, i.e. , we assume that sociability is an intrinsicindividual attribute that does not depend on their location. This later hypothesiscaptures the biological and behavioural aspect of hosts that can turn them into super-spreaders, i.e. , individuals that are highly efficient in transmitting the disease due to ahigh viral shedding [50] or because they have a high contact rate due to a pronouncedsocial behavior. However, other causes that are inherently related to the location, suchas the existence of high-risk scenarios related to work or leisure, are not captured bythe former assumption.Under the former hypothesis about the invariance of the connectivity k undermobility and assuming that those individuals with connectivity k move with probability p k , we can calculate the effective population of a patch i , ˜ n i , after the movement stagehas been performed, as the sum of the effective number of agents with connectivity k :˜ n i = (cid:88) k ˜ n [ k ] i . (4) isease dynamics with heterogeneous transmission and recurrent mobility patterns n [ k ] i is calculated considering the number of individuals withconnectivity k that travel from any patch j to i :˜ n [ k ] i = (cid:88) j n [ k ] j → i , (5)where n [ k ] j → i = [(1 − p k ) δ ij + p k R ji ] n j P j ( k ) . (6)Another quantity that can be evaluated is the effective connectivity distribution ofa patch, ˜ P i ( k ), defined as the probability of finding an individual of connectivity k inpatch i after the mobility stage. This probability is given by:˜ P i ( k ) = ˜ n [ k ] i ˜ n i . (7)From the effective connectivity distribution of a patch i we can measure the effectivemoments as: (cid:103) (cid:104) k n (cid:105) i = (cid:88) k k n ˜ P i ( k ) . (8) The coupling of interaction and mobility patterns of agents produces, for a given set ofmobility probabilities { p k } , a variation of the main structural attributes of the patches,as shown by the expressions of the effective population, Eqs. (4)-(5), and the effectiveconnectivity distribution, Eq. (7). These variations occur once the mobility step isperformed and become crucial when the spreading process (the interaction step of theMIR 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 individualsresiding in i that have connectivity k as I [ k ] i , implying that the total number of infectedresidents in i is I i = (cid:80) k I [ k ] i . Thus, the probability that an agent with residence inpatch i and connectivity k is infected is given by: ρ [ k ] i = I [ k ] i n [ k ] i . (9)The probabilities { ρ [ k ] i } (with i = 1,...,Ω and k = 1 , ..., k max ) constitute our dynamicalvariables. From these variables we can compute the fraction of infected individuals withresidence in patch i : ρ i = (cid:88) k P i ( k ) ρ [ k ] i , (10)or the fraction of infected individuals in the whole metapopulation: ρ = (cid:88) i n i N ρ i . (11) isease dynamics with heterogeneous transmission and recurrent mobility patterns { ρ [ k ] i } corresponding to the SIS dynamics we make use of the so-called heterogeneousmean-field theory (HMF) in the annealed regime [51]. Thus, after the movement stage,each susceptible agent with connectivity k that is placed in patch j connects randomlywith k individuals in the same patch and, for each infected contact, the susceptibleagent will become infected and infectious with probability ¯ λ . In addition, those infectedagents at time t will recover and become susceptible again with probability ¯ µ . Followingthese simple rules, the equations for the time evolution of the probabilities { ρ [ k ] i } read: ρ [ k ] i ( t + 1) = (1 − ¯ µ ) ρ [ k ] i ( t ) + (cid:104) − ρ [ k ] i (cid:105) Π [ k ] i ( t ) , (12)where Π [ k ] i ( t ) is the probability that a healthy individual with connectivity k andresidence in patch i becomes infected at time t :Π [ k ] i ( t ) = (1 − p k ) π [ k ] i ( t ) + p k N (cid:88) j =1 R ij π [ k ] j ( t ) , (13)where π [ k ] i ( t ) is the probability that an individual of connectivity k placed in patch i becomes infected at time t and reads: π [ k ] i ( t ) = 1 − (cid:89) k (cid:48) (cid:104) − ¯ λ ˜ P i ( k (cid:48) | k ) ˜ ρ [ k (cid:48) ] i ( t ) (cid:105) k . (14)In the former expression, ˜ P i ( k (cid:48) | k ) is the probability that an agent with connectivity k placed in patch i is connected with another agent with k (cid:48) placed in the same patch. Inaddition, ˜ ρ [ k ] i is the effective fraction of infected individuals with connectivity k placedin patch i : ˜ ρ [ k ] i = ˜ I [ k ] i ˜ n [ k ] i = (cid:88) j I [ k ] j → i ˜ n [ k ] i = (cid:88) j n [ k ] j → i ρ [ k ] j ( t )˜ n [ k ] i , (15)where the denominator is given by Equation (5) and the numerator is the number ofinfected individuals that are in patch i .In the following we will consider that the contact networks created at eachinteraction step are completely uncorrelated. This way, the probability ˜ P i ( k (cid:48) | k ) canbe written in terms of the effective connectivity distribution of patch i as:˜ P i ( k (cid:48) | k ) = k (cid:48) ˜ P i ( k (cid:48) ) (cid:102) (cid:104) k (cid:105) i = k (cid:48) ˜ n [ k (cid:48) ] i (cid:88) k (cid:48)(cid:48) k (cid:48)(cid:48) ˜ n [ k (cid:48)(cid:48) ] i , (16)which is the probability of selecting an edge from an individual with connectivity k (cid:48) placed in patch i . isease dynamics with heterogeneous transmission and recurrent mobility patterns Figure 2.
Example of a star-like metapopulation network with κ + 1 patches. In thisexample, 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 connectivity1 with probability η , or k max = 20 with complementary probability 1 − η , and eachleaf is a homogeneous patch with residents of same connectivity (cid:104) k (cid:105) l = β (cid:104) k (cid:105) h , with β = 1 and (cid:104) k (cid:105) h = 5. The flow from hub to a leaf happens with probability R hl = κ − ,from leaves to hub with R lh = δ , and between adjacent leaves with R l,l +1 = 1 − δ , incounterclockwise direction.
3. Metapopulations with heterogeneous subpopulations
Once the Markovian equations have been derived, we now study the impact ofheterogeneous distributions of individual contacts by using synthetic metapopulations.This will allow us to validate the Markovian equations derived by comparing the resultsobtained by the iteration of Eqs. (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 andset of connectivity distributions, we restrict our analysis, as in [30], to synthetic star-like metapopulation networks. Our choice is rooted in their versatility for, despite beingsimplistic structures, star-like metapopulations exhibit a wide variety of regimes causedby the non-uniform distribution of the population across patches and the asymmetry inth mobility patterns connecting them. This kind of synthetic metapopulation, shown infigure 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] ofthe hub population, n l = αn h . The mobility towards leaves of individuals with residence isease dynamics with heterogeneous transmission and recurrent mobility patterns R hl = 1 κ , (17)while the mobility of those residents in the leaves is controlled by a parameter δ . Thisway, a resident in a leave l that decides to move will go to the hub with probability δ , R lh = δ, (18)or move to the next (counterclockwise direction) leave with probability R l,l +1 = 1 − δ. (19)Note that the choice of the direction of movements among leaves is not relevantas 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 inRef. [30], being characterized by two parameters α and δ . However, the syntheticmetapopulations used here gets rid of the assumption of homogeneous (all-to-all) contactpatterns in the patches. To this aim, and keeping the symmetry of the original star-likemetapopulations, we consider that the residents of the central patch (the hub) havea 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 thatthe connectivity distribution of the individuals belonging to the hub is bimodal: P h ( k ) = ηδ k + (1 − η ) δ kk max , (20) i.e. , agents in the hub have connectivity 1 with probability η and connectivity k max withprobability (1 − η ). This way, the n -th moment of the hub’s connectivity distributionis: (cid:104) k n (cid:105) h = (cid:88) k k n P h ( k ) = η + (1 − η ) k n max . (21)In their turn, those individuals belonging to leaves have the same number of contacts( (cid:104) k (cid:105) l ): P l ( k ) = δ k (cid:104) k (cid:105) l . (22)Note that the values of η and k max are correlated if we impose the additionalconstraint that the hub has an average connectivity (cid:104) k (cid:105) h fixed. In this case, given avalue k max , the value of η that allows it is given by: η = k max − (cid:104) k (cid:105) h k max − . (23)In this simple configuration, the heterogeneous nature of the contacts is two-fold. From a microscopic point of view, the bimodal distribution existing inside thecentral node induces local heterogeneities in the contacts made by residents there,which are controlled by parameters η and k max . In its turn, another global connectivityheterogeneity emerges driven by the asymmetry existing between the connectivity ofresidents of the hub and the leaves. In particular, we will assume throughout themanuscript that (cid:104) k (cid:105) l = β (cid:104) k (cid:105) h , with β ∈ [0 , (cid:104) k (cid:105) h = 5, k max = 20, and α = β = 1. isease dynamics with heterogeneous transmission and recurrent mobility patterns To check the validity of the Markovian equations, we define a MC algorithm for thestochastic simulation of the SIS model on top of a metapopulation with heterogeneouscontact patterns. As in the case of Markovian equations, Eqs. (12)-(14), the processproposed is also a discrete-time dynamics. At each time step t , each individual is testedto 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 individualwith connectivity k chooses randomly k individuals in the patch they currently occupyand are infected with probability ¯ λ if the contacted individual is infectious. Once allthe potential infections events have been simulated, healing happens with probability¯ µ for each infected individual at time t −
1. In this sense, we perform a synchronousupdate of the state of the entire metapopulation.The simulation procedure in a give time step t can be summarized as follows:(i) First, a fraction ρ ini of the population is randomly infected as the initial condition.(ii) For each patch i , each individual with connectivity k resident in i is tested to movewith probability p k . If she or he moves, a patch j is chosen proportionally to R ij .(iii) Each susceptible individual with connectivity k selects at random at maximum k contacts in patch i . For each attempt, it can be infected with probability:¯ λ (cid:88) k k ˜ I [ k ] i (cid:88) k k ˜ n [ k ] i , (24)or remains susceptible with complementary probability. These attempts stop whenthe individual becomes infected.(iv) Each individual with state infected at time step t − t withprobability ¯ µ .(v) 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 × − ofindividuals at random when this state is reached [52, 53]. This keeps the dynamicsalways active and the equilibrium state is defined after comparing averages oversequential time windows of size T = 100, and accepting if the absolute difference issmaller than t cvg = 10 − . The comparisons between MC and Markovian equations are performed in star-likemetapopulations with κ = 10 and α = 1, i.e. , in which all patches (hubs and leaves)contain the same number of individuals ( n l = n h = 10 individuals per patch), to focuson the effect of contact heterogeneity. Furthermore, for the same reason, we focus onthe case that mobility is independent of the connectivity of individuals, p k = p ∀ k . isease dynamics with heterogeneous transmission and recurrent mobility patterns . . . . . . . ¯ λ c / ¯ λ . . . . . ρ ∗ . . . . . . p . . . . . . . ¯ λ c / ¯ λ . . . . . ρ ∗ . . . . . . p (a) (b) Figure 3.
Equilibrium regimes of the Markovian equations (lines) and MC simulations(symbols) for a star-like metapopulation with n h = n l = 10 and κ = 10. The hubcontains individuals with connectivity (cid:104) k (cid:105) h = 100 ( η = 0 and k max = 100), and theleaves (cid:104) k (cid:105) l = 10 ( β = 0 . δ = 0 . .
9. A fraction ρ pump = 2 × − and 10 stochastic samples were used for the MCsimulations. First we neglect local heterogeneities and consider that contact heterogeneity onlyhappens between patches. In mathematical terms, this assumption implies that thepopulation of the hub has an homogeneous contact distribution ( η = 0) although itsmean connectivity (cid:104) k (cid:105) h = k max is different from that of the leaves (cid:104) k (cid:105) l = β (cid:104) k (cid:105) h , with β (cid:54) = 1. In particular, in figure 3 we plot the mean epidemic prevalence ρ ∗ in theequilibrium state as a function of the infection probability ¯ λ scaled by the epidemicthreshold in the case of null mobility ¯ λ c ( p = 0) = ¯ λ . To derive the latter quantity,we realize that the absence of flows among the patches precludes the interaction amongthe residents in different areas, so the epidemic threshold corresponds to the well-knownexpression provided by HMF equations [51] for the most vulnerable patch. Therefore,¯ λ = ¯ µ min (cid:26) (cid:104) k (cid:105) h (cid:104) k (cid:105) h , (cid:104) k (cid:105) l (cid:104) k (cid:105) l (cid:27) . (25)We consider that (cid:104) k (cid:105) h = 100 while leaves have (cid:104) k (cid:105) l = 10 ( β = 0 .
1) and explore twodifferent mobility patterns. In particular, in (a) we set δ = 0 . i.e. , passing from one leave to another and avoidingthe hub. In this case, the so-called epidemic-detriment by mobility shows up so that theepidemic state is delayed as the mobility p increases, with the exception of very largevalues of p . However, note that, at variance with [30], here the populations of hubsand leaves are identical; we will explore the roots of this detriment below. Second, inpanel (b), we set δ = 0 . p < .
5, while for p > . isease dynamics with heterogeneous transmission and recurrent mobility patterns . . . . . . . ¯ λ c / ¯ λ . . . . . . . ρ ∗ . . . . . . p . . . . . . . ¯ λ c / ¯ λ . . . . . . . ρ ∗ . . . . . . p (a) (b) Figure 4.
Equilibrium regimes of the Markovian equations (lines) and MC simulations(points) for n h = n l = 10 and κ = 10. The patches contain individuals with power-law connectivity distributions P i ( k ) ∼ k − ¯ γ i , k ∈ [3 , γ h = 2 . γ l = 3 . δ = 0 . .
9. Afraction ρ pump = 2 × − and 10 stochastic samples were used for the MC simulations. Next we analyze a star-like metapopulation that generalizes the contactheterogeneity of the first one. In this case the hub is very heterogeneous, containinga power-law distribution, P h ( k ) ∼ k − ¯ γ h with ¯ γ h = 2 .
3, while leaves have also a power-law distribution P l ( k ) ∼ k − ¯ γ l with ¯ γ l = 3 .
5, both with k ∈ [3 , δ = 0 . δ = 0 .
9, showing similar qualitative behaviors with the mobility, namely theemergence of epidemic detriment, to those found in figure 3. Quantitatively, it is worthstressing that the existence of strong local heterogeneities within both hub and leaves inabsence of mobility will lead to an activation described by the HMF theory, in which theepidemic prevalence approaches zero close to the epidemic threshold as ρ ∼ (¯ λ − ¯ λ c ) ¯ β where ¯ β >
4. Epidemic threshold
Figures 3 and 4 reveal that the epidemic detriment emerges even when dealing withhomogeneously distributed populations, at variance with [30], in which increasingmobility in homogeneous populations favors epidemic spreading by anticipating theepidemic threshold, ¯ λ c , here defined as as the minimum infectivity per contact, ¯ λ ,such that an epidemic state can be stable. Therefore, the emergence of epidemicdetriment here should be rooted in the interplay among contact heterogeneities andhuman mobility. In this section, we aim at deriving an analytical expression of theepidemic threshold, ¯ λ c , to shed light on the mechanisms giving rise to the behavior isease dynamics with heterogeneous transmission and recurrent mobility patterns ρ [ k ] i ( t + 1) = ρ [ k ] i ( t ) = ρ ∗ i [ k ] . Under this assumption, Eq. (12) reads:¯ µρ ∗ i [ k ] = (cid:104) − ρ ∗ i [ k ] (cid:105) Π ∗ i [ k ] (26)with Π ∗ i [ k ] = (cid:34) (1 − p k ) π ∗ i [ k ] + p k N (cid:88) j =1 R ij π ∗ j [ k ] (cid:35) . (27)Furthermore, for ¯ λ values close to the epidemic threshold, the fraction of infectedindividuals is negligible, which means that ρ ∗ i [ k ] = ¯ (cid:15) ik (cid:28) ∀ ( i, k ). This fact allows usto linearize the equations characterizing the steady state of the dynamics by neglectingall the terms O (¯ (cid:15) ). In particular, the probability that an individual with connectivity k and placed in i contracts the disease, π ∗ i [ k ] , can be approximated by π ∗ i [ k ] = 1 − (cid:89) k (cid:48) (cid:104) − ¯ λ ˜ P i ( k (cid:48) | k ) ˜ ρ ∗ [ k (cid:48) ] i (cid:105) k (cid:39) ¯ λk (cid:88) k (cid:48) ˜ P i ( k (cid:48) | k ) ˜ ρ ∗ [ k (cid:48) ] i , (28)where we have used that O ( ˜ ρ ) = O (¯ (cid:15) ) as shown by Eq. (15). In particular, pluggingEqs. (15-16) into the last expression leads to: π ∗ i [ k ] = ¯ λkQ i (cid:88) k (cid:48) k (cid:48) (cid:88) j [(1 − p k (cid:48) ) δ ij + p k (cid:48) R ji ] n j P j ( k (cid:48) )¯ (cid:15) jk (cid:48) , (29)where Q i ≡ (cid:88) k k (cid:88) j [(1 − p k ) δ ij + p k R ji ] n j P j ( k ) (30)is the effective number of edges in patch i . Note that (cid:80) i Q i = (cid:80) k (cid:80) j kP j ( k ) n j is thetotal number of edges in the system, a conserved quantity. After introducing Eq. (29)and some algebra, Eq. (27) transforms into:Π ∗ i [ k ] = ¯ λ (cid:88) j (cid:88) k (cid:48) ¯ M jk (cid:48) ik ¯ (cid:15) jk (cid:48) , (31)where ¯ M jk (cid:48) ik = kk (cid:48) P j ( k (cid:48) ) (cid:20) (1 − p k )(1 − p k (cid:48) ) δ ij Q i + (1 − p k ) p k (cid:48) R ji Q i + p k (1 − p k (cid:48) ) R ij Q j + p k p k (cid:48) (cid:88) l R il R jl Q l (cid:35) n j . (32)Finally, if we introduce these values into Eq. (26) and retain only linear terms in (cid:15) ,we arrive to the following expression¯ µ ¯ (cid:15) ik = ¯ λ (cid:88) j (cid:88) k (cid:48) ¯ M jk (cid:48) ik ¯ (cid:15) jk (cid:48) , (33) isease dynamics with heterogeneous transmission and recurrent mobility patterns λ c = ¯ µ Λ max ( ¯M ) . (34)The elements of matrix ¯M given by Equation (32) represent four types ofinteractions in the metapopulation. Basically, the element ¯ M jk (cid:48) ik represents theprobability that a resident of patch i with connectivity k is in contact with anotherindividual of patch j and connectivity k (cid:48) . The first term accounts for interactions ofresidents of the patch, that do not move. In second term, an individual of i stays andinteracts with a traveler from patch j in patch i , that arrived with probability p k (cid:48) 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 (cid:48) R il R jl . Equation (33) computes the exact expression of the epidemic threshold in presence ofheterogeneous contact patterns. However, its computation involves solving the spectraof a matrix whose dimension is determined by the number of connectivity classes andpatches in the metapopulation. In particular, in presence of highly heterogeneouspopulations with fine spatial resolution, this problem can be computationally very harddue 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 willconsiderably reduce the complexity of the problem as proved below.Before going ahead, it is convenient to make the transformation ¯ (cid:15) ik (cid:55)→ k(cid:15) ik inEquation (33), resulting in ¯ µ(cid:15) ik = ¯ λ (cid:88) j (cid:88) k (cid:48) M jk (cid:48) ik (cid:15) jk (cid:48) , (35)where the elements of the new matrix M read as M jk (cid:48) ik = k (cid:48) P j ( k (cid:48) ) (cid:20) (1 − p k )(1 − p k (cid:48) ) δ ij Q i + (1 − p k ) p k (cid:48) R ji Q i + p k (1 − p k (cid:48) ) R ij Q j + p k p k (cid:48) (cid:88) l R il R jl Q l (cid:35) n j . (36)If p k = p , Equation (35) becomes independent of k , which allows a dimensionalityreduction of the matrix. In particular, Equation (35) reads:¯ µ(cid:15) i = ¯ λ (cid:88) j M ij (cid:15) j , (37)and the elements of the reduced matrix M are given by: M ij = (cid:104) k (cid:105) j (cid:34) (1 − p ) δ ij Q i + p (1 − p ) (cid:18) R ji Q i + R ij Q j (cid:19) + p (cid:88) l R il R jl Q l (cid:35) n j , (38) isease dynamics with heterogeneous transmission and recurrent mobility patterns . . . . . . p ¯ λ c / ¯ λ − − − − − ρ ∗ . . . . . . p ¯ λ c / ¯ λ (50 , . , . , .
9) (100 , . , . , . (a) (b) Figure 5.
Dependence of the epidemic threshold on the mobility parameter p . Allpatches have the same population ( α = 1), where the hub has agents with a bimodalconnectivity distribution with k = 1 or k max , fixing (cid:104) k (cid:105) h = 5, and the leaves haveagents with connectivity (cid:104) k (cid:105) l = 5 ( β = 1). (a) Comparison of the theoretical epidemicthreshold obtained using Equation (40) (solid line), scaled by its value for p = 0,and the steady values of the prevalence ρ obtained from Eqs. (12)-(14). (b) Relativeepidemic threshold for different configurations ( k max , δ ), shown in the legends, withsolid and dashed lines for k max = 50 and 100, respectively, and different values of δ . where the effective number of edges Q i is now expressed as Q i = (cid:88) j (cid:104) k (cid:105) j [(1 − p ) δ ij + pR ji ] n j . (39)Once matrix M is constructed the epidemic threshold is computed as¯ λ c = ¯ µ Λ max ( M ) . (40)To test the accuracy of the former expression for the epidemic threshold, we compareits value computed according to Equation (40) with the heat map of the steady state ofthe dynamics obtained from the iteration of Eqs. (12)-(14). Figure 5(a) reveals that thetheoretical prediction of the epidemic threshold by Equation (40) is very accurate andcaptures the dependence of the epidemic threshold on the mobility p . This thresholdincreases while promoting mobility until it reaches a maximum at p = p ∗ since theinfection is gradually reduced in the hub as p increases, and the activation is thentriggered 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 forEq. (40) retrieving, as expected, the expression for the epidemic threshold providedby HMF equations on contact networks. Moreover, to quantify the effects of promotingmobility among disconnected patches, we perform a perturbative approach to the latterthreshold which holds for small p values in Appendix C. Interestingly, at variance withthe perturbative analysis carried out for (non-structured) well mixed metapopulationsin [30], here the linear correction of the epidemic threshold strongly depends on thetopological properties of the metapopulation. isease dynamics with heterogeneous transmission and recurrent mobility patterns . . . . . . γ . . . . . . α ¯ λ c ( p ∗ ) / ¯ λ . . . . . β . . . . . . α . . . . . ¯ λ c ( p ∗ ) / ¯ λ (a) (b) Figure 6.
Heat maps of the relative magnitude of the peak of the epidemic threshold¯ λ c ( p ∗ ) / ¯ λ as a function of α , β , and γ , with (cid:104) k (cid:105) h = 5. In (a), all patches have the sameaverage connectivity, with β = 1 while the local heterogeneity of the hub is modulatedby γ . Dashed lines correspond to the values of γ for k max = 100, 50, and 20, fromleft to right. The plot in (b) considers a fixed value of γ ≈ . k max = 100 when β = 1, tuning the connectivity of the leaves with β . The populationasymmetry is modulated by α for all cases. In what follows, to shed light on the nature of the epidemic detriment, we aim atquantifying the impact of the different components of the formalism, namely theunderlying metapopulation structure and the contact heterogeneities existing among itspopulation, on the relative magnitude ¯ λ c ( p ∗ ) / ¯ λ . To simplify this analysis, we will focuson the case of mobility independent of k , p k = p , and consider the configuration definedin Sec. 3.1, in which the hub has individuals with connectivity either 1 or k max , withfixed average connectivity (cid:104) k (cid:105) h , and the ones of the leaves have the same connectivity (cid:104) k (cid:105) l = β (cid:104) k (cid:105) h . For the sake of clarity, let us also express (cid:104) k (cid:105) l = γ (cid:104) k (cid:105) h . Note that inthis configuration the values of η and k max are correlated by Eq. (23), while γ is alsocorrelated with β and k max via γ = β (cid:104) k (cid:105) h (cid:104) k (cid:105) h ( k max + 1) − k max . (41)First, we fix α = β = 1, so that n l = n h and (cid:104) k (cid:105) h = (cid:104) k (cid:105) l , to study the effects ofvarying either the local heterogeneity existing in the hub by tuning k max or the flowsfrom leaves to the hub with δ in figure 5(b). Fixing k max = 50 and changing δ , itbecomes clear that the increase of δ leads to a decrease of p ∗ as a consequence of thehigher mixing among individuals from the central node and the leaves, but does notchange the relative magnitude ¯ λ c ( p ∗ ) / ¯ λ .The former beneficial effect is rooted in the homogenization of the connectivitydistribution driven by the mixing among individuals from the hub and the leaves.Interestingly, the position of the peak p ∗ remains unaltered when keeping δ constant,suggesting that the qualitative dependence of the epidemic threshold is governed by isease dynamics with heterogeneous transmission and recurrent mobility patterns p . Moreover, for small values of p , the behavior does notdepend on the local heterogeneities of the patches, as shown by a pertubative analysisin Appendix C. Quantitatively, it becomes clear that increasing the degree heterogeneityin the central node boosts the beneficial effect of the mobility, since the homogenizationeffect gains more relevance due to the higher vulnerability of the central node.Finally, we extend our analysis to cover populations distributed heterogeneouslyacross the metapopulation. In particular, we are interested in determining howthe population asymmetry α and the local connectivity heterogeneity η shape therelative magnitude of the peak of the epidemic threshold. To this aim, we represent¯ λ c ( α, β, γ ; p ∗ ) / ¯ λ ( β, γ ) in figure 6, for n l = αn h , (cid:104) k (cid:105) l = β (cid:104) k (cid:105) h , and (cid:104) k (cid:105) l = γ (cid:104) k (cid:105) h , inwhich γ is given by Eq. (41) for the constraints imposed in Sec. 3.1. We can observe that,as in figure 5(b), increasing the local heterogeneity of the hub (lowering γ ) increases thebeneficial effect of the population mixing, as shown in figure 6(a). Interestingly, if wefix γ and study the dependence of ¯ λ c ( α, β, γ ; p ∗ ) / ¯ λ ( β, γ ) with α and β , as shown infigure 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 reducingthe population of the periphery nodes, i.e. decreasing α , agents in the leaves are notable to substantially modify the connectivity distribution of residents in the hub, thushindering the detriment effect in all investigated cases.
5. Conclusions
Driven by the advance of data mining techniques in mobility and social patterns [1, 55,56, 57], epidemic models are continuously refined to bridge the gap existing betweentheir theoretical predictions and the outcomes of real epidemic scenarios. In particular,within the very diverse realm of epidemic models, the proliferation of data setscapturing human movements across fine spatial scales have prompted the evolution ofmetapopulation frameworks, which constitute the usual approach to study the interplaybetween human mobility and disease spreading. In this sense, the first theoreticalframeworks assuming the population to move as random walkers across syntheticmetapopulations [17] have given rise to models incorporating the recurrent nature ofhuman mobility [30, 58, 59, 28], the socio-economic facets of human movements [32, 60]or high-order mobility patterns [31].While most of the advances previously described have been focused on capturingthe mobility flows more accurately, less attention has been paid to improve the contactpatterns within each subpopulation. With few exceptions, such as the model recentlyproposed in [61] incorporating the time varying nature of social contacts, humaninteractions are usually modeled using well-mixing hypothesis that do not capture theheterogeneous nature of human interactions and the role that this social heterogeneityhas on the so-called super-spreading events. In this work, we tackle this challenge andadapt the metapopulation model presented in [30] to account for the heterogeneity isease dynamics with heterogeneous transmission and recurrent mobility patterns isease dynamics with heterogeneous transmission and recurrent mobility patterns Acknowledgments
W.C. acknowledges financial support from the Coordena¸c˜ao de Aperfei¸coamento dePessoal de N´ıvel Superior, Brazil – CAPES, Finance Code 001. A.A. acknowledgesfinancial support from Spanish MINECO (grant PGC2018-094754-BC21), Generalitatde Catalunya (grant No. 2017SGR-896), and Universitat Rovira i Virgili (grant No.2019PFR-URVB2-41). A.A. also acknowledges support from Generalitat de CatalunyaICREA Academia, and the James S. McDonnell Foundation (grant 220020325). W.C.and S.C.F. acknowledge financial support from CAPES (grant No. 88887.507046/2020-00), Conselho Nacional de Desenvolvimento Cient´ıfico e Tecnol´ogico – CNPq (grants No.430768/2018-4 and 311183/2019-0) and Funda¸c˜ao de Amparo `a Pesquisa do Estado deMinas Gerais – FAPEMIG (grant No. APQ-02393-18). D.S.P. and J.G.G. acknowledgefinancial support from MINECO (projects FIS2015-71582-C2 and FIS2017-87519-P),from the Departamento de Industria e Innovaci´on del Gobierno de Arag´on y FondoSocial Europeo (FENOL group E-19), and from Fundaci´on Ibercaja and Universidad deZaragoza (grant 224220).
Appendix A. Exact evaluation of the epidemic threshold for a star-likemetapopulation
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 − : 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 theseexpressions in Equation (38), we have isease dynamics with heterogeneous transmission and recurrent mobility patterns M hh = (cid:104) k (cid:105) h (cid:20) (1 − p ) Q h + p κ Q l (cid:21) n h , (1.1 a ) M lh = (cid:104) k (cid:105) h (cid:20) (1 − p ) p (cid:18) κ Q l + δ Q h (cid:19) + p (1 − δ ) κ Q l (cid:21) n h , (1.1 b ) M hl = (cid:104) k (cid:105) l (cid:20) (1 − p ) p (cid:18) δ Q h + 1 κ Q l (cid:19) + p (1 − δ ) κ Q l (cid:21) n l , (1.1 c ) M ll = (cid:104) k (cid:105) l (cid:20) (1 − p ) Q l + p (1 − δ ) Q l + p δ Q h (cid:21) n l , (1.1 d ) M l,l +1 = (cid:104) k (cid:105) l (cid:20) (1 − p ) p (1 − δ ) Q l + p δ Q h (cid:21) n l , (1.1 e ) M l,l − = (cid:104) k (cid:105) l (cid:20) (1 − p ) p (1 − δ ) Q l + p δ Q h (cid:21) n l , (1.1 f ) M ln = (cid:104) k (cid:105) l (cid:18) p δ Q h (cid:19) n l . (1.1 g )Again, by evaluating Equation (37), we have, for the hub,¯ µ(cid:15) h = ¯ λ (cid:88) j M hj (cid:15) j = ¯ λM hh (cid:15) h + κ ¯ λM hl (cid:15) l , (1.2)while for a leaf we have¯ µ(cid:15) l = ¯ λ (cid:88) j M lj (cid:15) j = ¯ λM lh (cid:15) h + ¯ λM ll (cid:15) l + ¯ λM l,l +1 (cid:15) l + ¯ λM l,l − (cid:15) l + ¯ λ ( κ − M ln (cid:15) l , (1.3)in which the factor 3 in the last term is since there are κ − R ln = 0). The statistical equivalence of the leaves allows usto recast the computation of the epidemic threshold in a eigenvalue problem of a 2 × M = (cid:32) M hh κM hl M lh M ll + M l,l +1 + M l,l − + ( κ − M ln , (cid:33) . (1.4)The leading eigenvalue will be given by Λ max = Tr M + (cid:112) (Tr M ) − M Appendix B. Epidemic threshold in the static case
To check the consistency of these equations, let us consider the static case in whichall individuals stay in their patches and do not move: p k = 0 ∀ k . So, Equation (38)becomes M ij | p k =0 = (cid:104) k (cid:105) j δ ij Q i | p k =0 n j , isease dynamics with heterogeneous transmission and recurrent mobility patterns Q i | p k =0 = n i (cid:104) k (cid:105) i , that after being used in Equation (37) results in ¯ µ(cid:15) i = ¯ λ (cid:104) k (cid:105) i (cid:104) k (cid:105) i (cid:15) i . This case consists of isolated subpopulations in an annealed regime in which the epidemicthreshold will be given by the first subpopulation in the active state, if its populationis not so small compared to other patches. Indeed, the usual epidemic threshold knownin the HMF theory is obtained, ¯ λ c = ¯ µ min i (cid:26) (cid:104) k (cid:105) i (cid:104) k (cid:105) i (cid:27) . (2.1)Therefore, the epidemic threshold of the metapopulation corresponds to the individualepidemic threshold of the most vulnerable patch in the static case. Appendix C. Perturbative analysis of the epidemic threshold
We proceed by making a perturbative analysis of the eigenvalues of the matrix M up tofirst order on p to complement the discussions of the main text. First, it is convenientto rewrite Equation (38) to split the terms with different order in p : M ij = (cid:104) k (cid:105) j (cid:26) δ ij Q i + p (cid:20) R ji Q i + R ij Q j − δ ij Q i (cid:21) + p (cid:34) δ ij Q i − R ji Q i − R ij Q j + (cid:88) l R il R jl Q l (cid:35)(cid:41) n j . (3.1)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 (cid:104) k (cid:105) i . The first derivative of Q i isd Q i d p (cid:12)(cid:12)(cid:12)(cid:12) p =0 = (cid:88) j (cid:104) k (cid:105) j ( R ji − δ ij ) n j . Let us define r i ≡ (cid:88) j ( − R ji + δ ij ) n j (cid:104) k (cid:105) j , so that dd p (cid:18) Q i (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) p =0 = r i ( n i (cid:104) k (cid:105) i ) . Next, keeping only terms up to order 1, we have1 Q i = 1 n i (cid:104) k (cid:105) i + p r i ( n i (cid:104) k (cid:105) i ) + O ( p ) . Substituting the last expression in Equation (3.1) we get, after some algebra, M ij = (cid:103) M (0) ij + p (cid:103) M (1) ij + O ( p ) , (3.2) isease dynamics with heterogeneous transmission and recurrent mobility patterns . . . . . . p . . . . . . . . Λ / Λ κ = 10 δ = 0 . δ = 0 . δ = 0 . δ = 1 . − (1 + κδ ) p . . . . . p . . . . . Λ / Λ δ = 0 . κ = 5 κ = 10 κ = 20 κ = 501 − (1 + κδ ) p (a) (b) Figure C1.
Normalized leading eigenvalue of matrix M as a function of the mobilityfor different values of the number of leaves κ and the agents flow from leaves to thehub governed by δ , with the same number of individuals ( α = 1). Solid lines show theexact values whereas dotted lines corresponds to the linear correction estimated by theperturbative approach via Equation (3.6). We fix the number of leaves κ and modify δ (color code) in (a) and present the complementary analysis in (b). where (cid:103) M (0) ij = δ ij (cid:104) k (cid:105) i (cid:104) k (cid:105) i , (3.3 a ) (cid:103) M (1) ij = (cid:20) R ij n j (cid:104) k (cid:105) j + R ji n i (cid:104) k (cid:105) i + δ ij n i (cid:104) k (cid:105) i (cid:18) r i n i (cid:104) k (cid:105) i − (cid:19)(cid:21) n j (cid:104) k (cid:105) j . (3.3 b )From the static case, we know that there are Ω unperturbed eigenvalues Λ (0) i = (cid:104) k (cid:105) i / (cid:104) k (cid:105) i , for p = 0, with normalized eigenvectors (cid:126)(cid:15) i = { (cid:15) j } and (cid:15) j = δ ij ; see Equation(2.1). Assuming that the eigenvalues are not degenerate, the new eigenvalues will begiven by [63] Λ i ≈ Λ (0) i + p Λ (1) i , (3.4)where Λ (0) i = (cid:104) k (cid:105) i (cid:104) k (cid:105) i , (3.5 a )Λ (1) i = (cid:126)(cid:15) i (cid:103) M (1) (cid:126)(cid:15) i . (3.5 b )Substituting Equation (3.3 b ) in Equation (3.5 b ), after some algebra we get the firstcorrection to the eigenvalue,Λ (1) i Λ (0) i = R ii − − (cid:88) j (cid:54) = i R ji n j (cid:104) k (cid:105) j n i (cid:104) k (cid:105) i . (3.6)Interestingly, unlike the original MIR model, the first order correction dependson the underlying topology. To check the accuracy of this correction, we represent infigure C1 the leading eigenvalues of the matrix M along with the linear correctionprovided by the perturbative analysis, finding a remarkable agreement in the lowmobility regime p (cid:28) isease dynamics with heterogeneous transmission and recurrent mobility patterns References [1] Guimera R, Mossa S, Turtschi A and Amaral L N 2005
Proc. Natl. Acad. Sci.
Nature
Phys. Rep.
Epidemics Nature et al.
BMC Med. et al. Science
J. Infect. Dis.
S375–S379[9] Zhang Q, Sun K, Chinazzi M, Pastore y Piontti A, Dean N E, Rojas D P, Merler S, Mistry D,Poletti P, Rossi L, Bray M, Halloran M E, Longini I M and Vespignani A 2017
Proc. Natl. Acad.Sci.
E4334–E4343[10] Kraemer M U G, Yang C H, Gutierrez B, Wu C H, Klein B, Pigott D M, du Plessis L, Faria N R,Li R, Hanage W P, Brownstein J S, Layan M, Vespignani A, Tian H, Dye C, Pybus O G andScarpino S V 2020
Science
Proc. Natl. Acad.Sci.
Math. Biosci.
Math. Biosci.
J. Theoret. Biol.
Trends Ecol. Evol. Ecol. Lett. Nat. Phys. Phys. Rev. Lett. J. Theoret. Biol.
Proc. Natl. Acad.Sci.
Phys. Rev. E (4) 042820[22] Silva D H and Ferreira S C 2018 Chaos Nature
PLoS ONE e60069[25] Masucci A P, Serras J, Johansson A and Batty M 2013 Phys. Rev. E Nat. Phys. Phys. Rev. X Eur. Phys. J. B J. Theoret. Biol.
Nat. Phys. J. R. Soc. Interface Phys. Rev. X (3) 031039[33] Soriano-Pa˜nos D, Ghoshal G, Arenas A and G´omez-Garde˜nes J 2020 J. Stat. Mech. Theory Exp.
Phys. Rev. Research Phys. Rev. X isease dynamics with heterogeneous transmission and recurrent mobility patterns [36] Costa G S, Cota W and Ferreira S C 2020 Phys. Rev. Research Phys. Rev. Lett. (14) 3200–3203[38] Shen Z, Fang N, Weigong Z, Xiong H, Changying L, Chin D P, Zonghan Z and Schuchat A 2004 Emerg. Infect. Dis. Nature
Int. J. Infect. Dis. Cell Host Microbe Lancet
Emerg. Infect. Dis. New Sci.
Int. J. Infect. Dis. PLOS Biology Science eabe2424[48] Althouse B M, Wenger E A, Miller J C, Scarpino S V, Allard A, H´ebert-Dufresne L and Hu H2020 Stochasticity and heterogeneity in the transmission dynamics of SARS-CoV-2 (
Preprint arXiv:2005.13689 )[49] Meyerowitz E A, Richterman A, Gandhi R T and Sax P E 2020
Ann. Intern. Med.
M20–5008[50] Woolhouse M E J, Dye C, Etard J F, Smith T, Charlwood J D, Garnett G P, Hagan P, Hii J L K,Ndhlovu P D, Quinnell R J, Watts C H, Chandiwana S K and Anderson R M 1997
Proc. Natl.Acad. Sci. Rev. Modern Phys. Phys. Rev. E (4) 042308[53] Cota W and Ferreira S C 2017 Comput. Phys. Commun.
Phys. Rev. E (6) 066117[55] Chowell G, Hyman J M, Eubank S and Castillo-Chavez C 2003 Phys. Rev. E Proc. Natl. Acad. Sci.
Netw. Spat. Econ. Phys. Rev. E Phys. Rev. E (2) 022306[60] Bosetti P, Poletti P, Stella M, Lepri B, Merler S and De Domenico M 2020
Proc. Natl. Acad. Sci.
J. Roy. Soc. Interface PeerJ Comput. Sci. e103[63] Marcus R 2001 J. Phys. Chem. A105