- Email: [email protected]
Analysis of metapopulation epidemic process on arbitrary networks*
4th IFAC Conference on Analysis and Control of Chaotic Systems 4th IFAC Conference on Analysis Control of Chaotic Systems 4th IFAC Conference on and August 2015. Tokyo, Japan and 4th IFAC26-28, Conference on Analysis Analysis and Control Control of of Chaotic Chaotic Systems Systems August 26-28, 2015. Tokyo, Japan Available online at www.sciencedirect.com August 26-28, 2015. Tokyo, Japan August 26-28, 2015. Tokyo, Japan
ScienceDirect IFAC-PapersOnLine 48-18 (2015) 141–145
Analysis of metapopulation epidemic Analysis of metapopulation epidemic Analysis of metapopulation epidemic process on arbitrary networks process on arbitrary networks process on arbitrary networks
Taro Takaguchi ∗∗ Renaud Lambiotte ∗∗ ∗∗ Taro Takaguchi ∗∗ Renaud Lambiotte ∗∗ Taro Taro Takaguchi Takaguchi Renaud Renaud Lambiotte Lambiotte ∗∗ ∗ National Institute of Informatics, Tokyo 101-8430, Japan (e-mail: ∗ ∗ Institute of Informatics, Tokyo 101-8430, Japan (e-mail: ∗ National National Tokyo [email protected]) National Institute Institute of of tInformatics, Informatics, Tokyo 101-8430, 101-8430, Japan Japan (e-mail: (e-mail: t [email protected]) ∗∗ t [email protected]) [email protected]) Department of Mathematics and naXys, University of Namur, 5000 t ∗∗ ∗∗ Department of Mathematics and naXys, University of Namur, 5000 ∗∗ Department Mathematics naXys, Namur, of Belgium (e-mail:and [email protected]) Department of Mathematics and naXys, University University of of Namur, Namur, 5000 5000 Namur, Belgium (e-mail: [email protected]) Namur, Belgium (e-mail: [email protected]) Namur, Belgium (e-mail: [email protected]) Abstract: Epidemic process on networks is considered. The system is discribed as a metapopAbstract: Epidemic process on networks issubpopulation considered. The system is discribed as a metapopAbstract: Epidemic process on networks system is discribed as metapopulation network in which a node (e.g., city or and Abstract: Epidemic process on represents networks is is considered. considered. The The system is school) discribed asisaa connected metapopulation network in which a node represents subpopulation (e.g., city or school) and is connected ulation network in which a node represents subpopulation (e.g., city or school) and is connected to other nodes via undirected links. Particles represent the subject of infection (e.g., individuals) ulation network in which a node represents subpopulation (e.g., city or school) and is connected to nodes via undirected links. represent the subject of individuals) to other nodes viaeach undirected links. Particles Particles represent the from subject of infection infection (e.g., individuals) andother interact with other within nodes while migrating nodes to nodes(e.g., in the manner of to other nodes via undirected links. Particles represent the subject of infection (e.g., individuals) and interact with each other within nodes while migrating from nodes to nodes in the mannersize of and interact with each other within nodes while migrating from nodes to nodes in the random diffusion. The nonlinear dependence of contact ratefrom within a node on itsinpopulation and interact with each other within nodes while migrating nodes to nodes the manner manner of of random diffusion. The nonlinear dependence of contact rate within a node on its population size random diffusion. The nonlinear dependence of contact rate within a node on its population size is introduced, according to the recent findingofbased onrate emprical The impacts random diffusion. The nonlinear dependence contact withinphone-call a node on data. its population size is according to the finding based on emprical The impacts is introduced, according to recent finding on phone-call data. The impacts of introduced, the nonlinear dependence arerecent investigated for three aspects ofphone-call epidemic data. process: is introduced, according to the the recent finding based based on emprical emprical phone-call data. Theepidemic impacts of the nonlinear dependence are investigated for three aspects of epidemic process: epidemic of the nonlinear nonlinear dependence are investigated investigated fortransient three aspects aspects of epidemic epidemic process: process: epidemic epidemic threshold, infection size at stationary state, and dynamics. of the dependence are for three of threshold, infection size at stationary state, and transient dynamics. threshold, threshold, infection infection size size at at stationary stationary state, state, and and transient transient dynamics. dynamics. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Complex systems; Networks; Differential equations; Eigenmode analysis; Monte Keywords: Complex Keywords: Complex systems; systems; Networks; Networks; Differential Differential equations; equations; Eigenmode Eigenmode analysis; analysis; Monte Monte Carlo simulation Keywords: Complex systems; Networks; Differential equations; Eigenmode analysis; Monte Carlo simulation Carlo simulation Carlo simulation 1. INTRODUCTION recently observed in various emprical studies (Smith, et al. 1. recently observed various emprical (Smith, et al. 1. INTRODUCTION INTRODUCTION recently observed in emprical studies (Smith, et (2009); Schl¨ apfer,in al. (2014)). To studies get better insight of 1. INTRODUCTION recently observed inetvarious various emprical studies (Smith, et al. al. (2009); Schl¨ a pfer, et al. (2014)). To get better insight of (2009); aapfer, et (2014)). To better of process will insight introduce (2009); Schl¨ Schl¨ pfer,under et al. al.realistic (2014)).settings, To get get we better insight of Understanding and controlling spreading process of infec- epidemic epidemic process under realistic settings, we will introduce Understanding and controlling spreading process of infecepidemic process under realistic settings, we will introduce Understanding and controlling spreading process of infecthese two factors into standard model and theoretically epidemic process under realistic settings, we will introduce tious diseases are of much interstspreading in the research fields such these two factors into standard model and theoretically Understanding and controlling process of infectious diseases of interst in research fields such two into tious diseases are are public of much muchhealth, interstand in the the research fields (Ansuch these investigate its properties. these two factors factors into standard standard model model and and theoretically theoretically as epidemiology, network science tious diseases are of much interst in the research fields such investigate its properties. as epidemiology, public health, and network science (Aninvestigate its properties. as epidemiology, public health, health, and network network science (An- investigate its properties. derson & May (1991)). A standard framework to model as epidemiology, public and science (Anderson & (1991)). standard to 2. MODEL derson & May May (1991)). A A standard isframework framework to model model epidemics in structrured population so-called metapopderson & May (1991)). A standard framework to model 2. MODEL MODEL epidemics in structrured population is so-called metapop2. epidemics in structrured population is so-called metapop2. MODEL ulation model, in which nodes represent subpopulations epidemics in structrured population is so-called metapopulation model, in which nodes represent subpopulations suppose that the network consists of undirected and ulation model, in which which nodes represent subpopulations and links connectes nodenodes pairs. represent The subject of epidemic We ulation model, in subpopulations We suppose links. that the consists of undirected and links connectes pairs. subject of epidemic We network consists of and We network also suppose that the networkand is and links connectes node node pairs. The The subject of with epidemic We suppose suppose that that the the network consists of undirected undirected and process is represented by particles which interact each unweighted and links connectes node pairs. The subject of epidemic unweighted links. We also suppose that the network is process is represented by particles which interact with each unweighted links. We also suppose that the network is connected (i.e., one can move between any node pairs via process is represented by particles which interact with each unweighted links. We also suppose that the network is other within a node and between in awith random process is represented by move particles whichnodes interact each connected (i.e., one can move between any node pairs via other within a node and move between nodes in a random connected (i.e., one can move between any node pairs via its links) and simple (i.e., no self-loops or multiple links). other within a node and move between nodes in a random connected (i.e., one can move between any node pairs via manner. For example, in airport networks, nodes are cities other within a node and move between nodes in a random its links) and simple (i.e., no self-loops or multiple links). manner. For example, airport nodes are and (i.e., multiple links). Letlinks) us define N nodes inno theself-loops networkor np particles manner. For and example, in airport networks, networks, nodes are cities cities its links) and simple simple (i.e., no self-loops orand multiple links). or countries links in represent connections between air- its manner. For example, in airport networks, nodes are cities Let us define N nodes in the network and n particles p or countries and links represent connections between airLet us define N nodes in the network and n particles p on it. Therefore, the average number of particles on a or countries and links represent connections between airLet us define N nodes in the network and np particles ports via flights. In this example, a particle between represents or countries and links represent connections air-a on it. Therefore, the average number of particles on ports via flights. In this example, a particle represents a on it. Therefore, the average number of particles on aaa node is given by ρ ≡ n /N . We define the rules governing ports via flights. In this example, a particle represents a on it. Therefore, the average number of particles on p passenger traveling cities by air. Therepresents theoreticala node is given by ρ ≡ n /N . We define the rules governing ports via flights. In between this example, a particle p passenger traveling between cities by air. The theoretical node is given by ρ ≡ n /N . We define the rules governing the states of by particles according to susceptible-infectedpassenger traveling betweenmodels cities by by air. The theoretical is given ρ ≡ npp /N . We define the rules governing studies of traveling metapopulation seek forThe thetheoretical impact of node passenger between cities air. the states of particles according to susceptible-infectedstudies metapopulation models seek the impact the particles according to (SIS) dynamics which has been investigated in studies oftopology metapopulation models seek for for on the large-scale impact of of susceptible the states states of of particles according to susceptible-infectedsusceptible-infectednetworkof and model parameters studies of metapopulation models seek for the impact of susceptible (SIS) which has investigated in network topology susceptible (SIS) dynamics dynamics which has been beenSalda˜ investigated in previous studies (Colizza, et al. (2007); na (2008); network topology and and model model parameters parameters on on large-scale large-scale susceptible (SIS) dynamics which has been investigated in outbreaks. network topology and model parameters on large-scale previous studies (Colizza, et al. (2007); Salda˜ n a (2008); outbreaks. previous&studies studies (Colizza, et Juher, al. (2007); (2007); Salda˜ naa Salda˜ (2008); Colizza Vespignani (2008); et al. (2009); n a outbreaks. previous (Colizza, et al. Salda˜ n (2008); outbreaks. & Vespignani (2008); Juher, et al. (2009); Salda˜ n aa In this paper, we take into accout two new aspects: Colizza Colizza & Vespignani (2008); Juher, et al. (2009); Salda˜ n (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. Colizza & Vespignani (2008); Juher, et al. (2009); Salda˜ n a In this take into accout new aspects: Masuda (2010); Iggidr, et al. (2012); Lund, et al. In this paper, paper, wefree take into accout two two and newnonlinear aspects: (2010); network topologywe from approximation (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. In this paper, we take into accout two new aspects: (2013); Juher & Ma˜ n osa (2014)). In SIS dynamics, the (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. network topology free from approximation and nonlinear & n (2014)). SIS dynamics, network topology free free contact from approximation approximation and recent nonlinear population-dependent rate. First, while re- (2013); (2013); Juher & Ma˜ Ma˜ nosa osa to (2014)). In SIS compartments: dynamics, the the network topology from and nonlinear state of Juher a particle belongs either In of two (2013); Juher & Ma˜ n osa (2014)). In SIS dynamics, the population-dependent contact rate. First, while recent restate of a particle belongs to either of two compartments: population-dependent contact rate. First, while recent research on metapopulation models has unveiled the impact state of a particle belongs to either of two compartments: population-dependent contact rate. First, while recent re- state susceptible (S) and infected (I). The states of particles of a particle belongs to either of two compartments: search on models has impact (S) and infected (I). states of particles search on metapopulation metapopulation models has unveiled unveiled the impact of heterogeneous distribution of degree (i.e., thethe number of susceptible susceptible (S) infected The states of search on metapopulation models has unveiled the impact evolve through the infection andThe recovery as susceptible (S) and and infected (I). (I). The states processes of particles particles of heterogeneous distribution of degree (i.e., the number through the infection and recovery processes as of heterogeneous distribution of degree (i.e., the number ofa evolve links connected todistribution a node) (Colizza, et (i.e., al. (2007); Salda˜ nof evolve through the infection and recovery processes as of heterogeneous of degree the number of follows: evolve through the infection and recovery processes as links connected to a node) (Colizza, et al. (2007); Salda˜ n a { follows: links connected to a node) (Colizza, et al. (2007); Salda˜ n a (2008); Colizza & (2008);etJuher, et al.Salda˜ (2009); follows: links connected to Vespignani a node) (Colizza, al. (2007); na follows: { S + I → I + I with rate β , i { (2008); Colizza & Vespignani (2008); Juher, et al. (2009); {S (2008); Colizza &Juher Vespignani (2008); Juher, the et al. al.majority (2009); + I → II + I with rate β , (1) Salda˜ naColizza (2010);& & Ma˜ nosa (2014)), (2008); Vespignani (2008); Juher, et (2009); + with rate rate µ. β S→ + IIS→ →I+ + II with βiii ,, (1) IS Salda˜ n & n (2014)), the (1) Salda˜ naaa (2010); (2010); Juher & Ma˜ Ma˜ nosa osatype (2014)), the majority majority of them assumedJuher the mean-filed approximation of (1) Salda˜ n (2010); Juher & Ma˜ n osa (2014)), the majority I → S with rate µ. II → S with rate µ. of them assumed the mean-filed type approximation of → S with rate µ. of them assumed the mean-filed type approximation of varies for different We assume that the infection rate β network equate alltype nodes with the same of them topology assumed which the mean-filed approximation of We assume that the infection rate βi varies for different network topology which equate all with same varies for that different We assume assume that the infection rate β βiii be network topologywhile which equate all nodes nodesstudies with the the same nodes i (1 that ≤ i the ≤ infection N ). It should noted the varies for different We rate degree. Second, most of previous presume network topology which equate all nodes with the same nodes i (1 ≤ i ≤ N ). It should beat noted noted thatnode. the degree. Second, while most of previous studies presume nodes i (1 ≤ i ≤ N ). It should that degree. Second, while most of of particles previous within studiesapresume presume infection occurs particlesbe the same i (1 ≤ ionly ≤ between N ). It should be noted that the the that theSecond, contact while rate between node is nodes degree. most previous studies infection occurs only between particles at the same node. that the contact rate between particles within a node is infection occurs only between particles at the same node. that the contact contact rate proportional between particles particles within node at is infection In addition to theonly reactions, particles move between occurs between particles at the same nodes node. irrelevant or linearly to the population that the rate between within aa node is addition to the reactions, particles between nodes irrelevant or proportional to the at In addition to reactions, particles move between nodes irrelevant or linearly linearly proportional to population the population population at In in manner of simple diffusion withmove diffusion coefficient In the addition to the the reactions, particles move between nodes the node, nonlinear dependence on the size are irrelevant or linearly proportional to the population at in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are D and D for S and I particles, respectively. In general, in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are DS and DI for S and I particles, respectively. In general, S II for S and D and II particles, respectively. general, This work was partly supported by Bilateral Joint Research S there is aD value of β (denoted by βc In and called D Dthreshold particles, respectively. In general, S and I for S and there is a threshold value of β (denoted by β and called c This work was partly supported by Bilateral Joint Research there is a threshold value of β (denoted by β and endemic If β of>β β of This work work wasJSPS, partly supported by Bilateral Bilateral Joint Research Research there is a threshold). threshold value (denoted by βcc fraction and called called Projects between Japan, and F.R.S.-FNRS, Belgium. T.T. was c a positive This was partly supported by Joint endemic threshold). If β > β a positive fraction of Projects between JSPS, Japan, and F.R.S.-FNRS, Belgium. T.T. was If β β fraction of Projects JSPS, Japan, and T.T. supported by JST, ERATO, Kawarabayashi LargeBelgium. Graph Project. Iendemic particlesthreshold). remains at state (called endemic threshold). If the β > >stationary βccc aa positive positive fraction the of Projects between between JSPS, Japan, and F.R.S.-FNRS, F.R.S.-FNRS, Belgium. T.T. was was supported by JST, ERATO, Kawarabayashi Large Graph Project. I particles remains at the stationary state (called the supported by by JST, JST, ERATO, ERATO, Kawarabayashi Kawarabayashi Large Large Graph Graph Project. Project. particles remains remains at at the the stationary stationary state state (called (called the the supported II particles Copyright 2015 IFAC 141 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright 2015 IFAC 141 Copyright © 2015 IFAC 141 Peer review© of International Federation of Automatic Copyright ©under 2015 responsibility IFAC 141Control. 10.1016/j.ifacol.2015.11.026
IFAC CHAOS 2015 142 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
endemic equilibrium). Otherwise, if β < βc , the fraction of I particles is zero in the long-term limit (called the extinction equilibrium). Dependence of βc on network structure is one of the main concerns in this study. We employ the continuous-time formulation of metapopulation SIS dynamics (Salda˜ na (2008)) which assumes the simultaneity of interaction between particles and diffusion of particles. Taking advantage of this formulation, we can describe the average behavior of the system by the deterministic reaction-diffusion equation (Masuda (2010)) as ∑ β ∂t ρS,i = µρI,i − i ρS,i ρI,i − DS Lij ρS,j , ρi j (2) ∑ βi ∂ ρ = −µρ + ρ ρ − D L ρ , t I,i I,i S,i I,i I ij I,j ρi j
where ρS,i and ρI,i represent the number of S and I particles at node i, respectively. The total number of particles at node i is given by ρi ≡ ρS,i + ρI,i . Matrix L denotes the random-walk Laplacian whose elements obey { 1 (i = j), Lij = (3) −Aij /kj (i = j),
where Aij is the (i, j)-element of adjacency matrix A and takes a binary value, i.e., unity if nodes i and j are adjacent and zero otherwise. Degree ki of node i is the total number of links associated to the node. Since we suppose the network is undirected, A is a symmetric matrix but L is asymmetric.
We assume that the infection rate βi at node i depends on ρi as follows: (4) βi = βργi (γ ≥ 0). This nonlinear dependence is motivated by recent empirical findings (Smith, et al. (2009); Schl¨ apfer, et al. (2014)) described in Sec. 1. Equation (4) was previously mentioned in (Anderson & May (1982); McCallum, et al. (2001)) and was recently considered in combination with meanfield networks (Juher & Ma˜ nosa (2014)). Equation (4) interpolates the two standard assumption: the so-called frequency interaction recovers if γ = 0 and the so-called mass interaction (Colizza, et al. (2007); Masuda (2010)) recovers if γ = 1. 3. RESULTS In this section, we report following three aspects of the metapopulation SIS model defined in Sec. 2. First, in Sec. 3.1, we investigate the endemic threshold βc . We cannot derive the exact βc value in general case, however, we derive its upper bound which appears to be fairly tight according to numerical simulations on some network models. Second, in Sec. 3.2, we investigate the infection size at the stationary state. In particular, we assume that diffusion of particles as the perturbation to the steady state without diffusion, and derive the infection size as a perturbed solution of the reaction-diffusion equation (Eq. (2)). Third, in Sec. 3.3, we investigate the impact of the population-dependent infection rate (Eq. (4)) on the transient process in which the configuration of particles converges from the initial to the stationary distribution. We show that the the population-dependent infection rate 142
induces nonmonotonical behavior in the convergene of population. The metapopulation SIS model (Eq. (2)) contains four parameters related to its time scale: β, µ, DS , and DI . Without loss of generality, we set µ = 1 and DS = DI for all the numerical simulaitons reported in this section. 3.1 Endemic threshold The results reported in this section are essentially in parallel with previous studies (Salda˜ na (2008); Juher, et al. (2009); Salda˜ na (2010); Masuda (2010); Lund, et al. (2013); Juher & Ma˜ nosa (2014)). Our results differ from these studies because we consider arbitrary networks free from the heterogeneous mean-field (HMF) approximation of networks, while these studies are based on this approximation, except for (Masuda (2010)). With the HMF approximation, the nodes with the same ki are supposed to behave in the same way on average (Salda˜ na (2008, 2010)). To be more precise, a pair of nodes with degree k and k are stochastically connected to each other with probability kk P (k)P (k )/k2 , where k is the average degree of the network, in an uncorrelated network. In addition, the population of particles at a node with k is determined by (ρS,k , ρI,k ) regardless of the nodes’ indices i. We are interested in endemic threshold βc . Since the extinction equilibrium ρ∗I,k = 0 is the stable if β < βc , βc is given as the β value with which the extinction equilibrium becomes unstable. By analysing the eigenvalue of the Jacobian at the equilibrium, one obtain the upper bound of βc (Salda˜ na (2008)) as )−γ ( kmax ρ β HMF ≡ (µ + DI ) (5) k In words, βc < β HMF holds true under the HMF approximation.
In the following, we prove that the upper bound β HMF also holds true for arbitrary networks beyond the HMF approximation. As we described for the derivation under the HMF approximation, we consider the stability of the extinction equilibrium (ρ∗S,i , ρ∗I,i ) = (ki /kρ, 0) (1 ≤ i ≤ N ) of the model (Eq. (2)). The allocation ρ∗S,i = ki /kρ streightforwardly follows as the stationary distribution of diffusion of single-type particles. The linear stability of the extinction equilibrium is determined by the Jacobian matrix at the equilibrium: ( (1) (3) ) J J J= (6) O J (2) All submatrices J (1) , J (2) , J (3) , and O are with size N ×N . Matrix O is zero matrix. The element of J (1) and J (2) is defined by (1)
Jij = −DS Lij , (7) )γ ) ( ( ki (2) ρ − µ δij − DI Lij . Jij = β (8) k When the largest eigenvalue of J is positive, the extinction equilibrium is linearly unstable. Since the eigenvalues of L are in 0 ≤ λ(L) ≤ 2 (Chung (1997)), the eigenvalues of J (1)
IFAC CHAOS 2015 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
are in −2DS ≤ λ(J (1) ) ≤ 0 and not positive regardless of parameter values. Therefore, to derive the condition for the endemic equilibrium, we investigate the largest eigenvalue of J (2) in the following. There are two special cases and we can exactly derive βc . As the first case, when γ = 0, the element of J (2) is equal to (2) Jij = (β − µ) δij − DI Lij . (9)
Therefore, its largest eigenvalue λmax (J (2) ) is given by (β − µ) and it leads to βc = µ. As the second case, when the network is regular and ki = k (1 ≤ i ≤ N ), the element of J (2) is equal to (2) Jij = (βργ − µ) δij − DI Lij . λmax (J (2) ) is given by (βργ − µ) and
(10)
Therefore, it leads to βc = µρ−γ . It should be noted that DI (> 0) is irrelevant to βc for these two cases. We now move to general cases with γ > 0 on a nonregular network. For the convenience for calculation, we symmetrize J (2) via the following similarity transformation (Masuda (2010)): ) ( 1 1 1 J (2) Jˆ(2) ≡diag √ , √ , . . . , √ k1 k2 kN (√ √ √ ) (11) × diag k1 , k2 , . . . , kN , where diag(a1 , a2 , . . . , aN ) represents a diagonal matrix in which the (i, i)-element is given by ai . With this transformation, the element of J (2) is converted into )γ ) ( ( Aij ki (2) ˆ ρ − µ − DI δij + DI √ . (12) Jij = β k ki kj
If λmax (Jˆ(2) ) is positive, the extinction equilibrium is unstable and the endemic equilibrium arises. We can represent λmax (Jˆ(2) ) by using the Rayleigh quotient as x Jˆ(2) x . (13) λmax (Jˆ(2) ) = max x x |x|=1
We show that βc < β HMF holds true for arbitrary networks as follows. Let us define x = (0, 0, . . . , 1, . . . , 0) in which only the i-th element is equal to unity and the rest equal to zero. With this x, )γ ( ki ρ − µ − DI . (14) x Jˆ(2) x = β k Therefore, )γ ( ki ρ − µ − DI ≤ λmax (Jˆ(2) ), (15) max β 1≤i≤N k and we obtain )−γ ( kmax β > (µ + DI ) ρ = β HMF . (16) k Equation (16) states that, if β > β HMF , λmax (Jˆ(2) ) is positive and the endemic threshold arises. In other words, βc < β HMF holds true. We show an improved upper bound of βc . We observed that the upper bound β HMF is achieved via approximation of the principal eigenvector of Jˆ(2) with x = (0, 0, . . . , 1, . . . , 0) in which only the element corresponding to a single node with kmax is equal to unity. Therefore, 143
143
we elaborate this approximation so that the elements of x corresponding to the nodes with the largest degree have nonzero values. Let us define x such that (ki < kc ), 0√ (17) xi = ki (ki ≥ kc ), Kc takes an integer value in [1, kmax ] where free parameter k c ∑ and Kc ≡ j:kj ≥kc kj is a normalization constant. With this choice of x, )γ ( ) ( ∑ 2Ec ki ki ˆ(2) ρ , β − µ − DI 1 − x J x= k Kc Kc i:ki ≥kc
(18) )γ ( ) 2Ec kc ρ − µ − DI 1 − , (19) ≥β k Kc where Ec represents the total number of links between nodes with ki ≥ kc . Therefore, x Jˆ(2) x λmax (Jˆ(2) ) = max x x |x|=1 )γ [ ( ] 2Ec kc ρ + DI − (µ + DI ) . ≥ max β kc k Kc (20) (
A sufficient condition for the endemic equilibrium is that there exists kc which satisfies )γ ( ) ( 2Ec kc ρ − 1− DI − µ > 0. (21) β k Kc In terms of β, this is equivalent to {[ )−γ } ( ) ]( 2Ec kc ρ = β. (22) DI β > min µ + 1 − kc Kc k This sufficient condition improves the previous one (given by Eq. (16)) because β ≤ β HMF and the equality holds true if kc = kmax . 3.2 Infection size at the steady state In the previous section, we discussed endemic threshold βc . Another important property of metapopulation SIS dynamics is the infection size at the stationary state when β > βc . In particular, we want to derive the infection size when diffusion coefficient DS = DI = ε is very small, in order to understand the impact of diffusion process on epidemic spreading. To this aim, we approximate the solution of Eq. (2) with ε > 0 by perturbing its solution with ε = 0 with small ε. The basic idea of the analysis reported in this section was originally proposed in (Sanders, et al. (2013)). We focus on the infection size at the stationary state, while the previous study (Sanders, et al. (2013)) considered its time evolution from the initial state. We assume that the initial number of particles at node i is proportional to ki , that is, ρi = ki ρ/k. When ε = 0, the time evolution of ρI,i is governed by ∂t ρI,i = −µρI,i + βργ−1 (ρi − ρI,i )ρI,i , i
(23)
where we substitute ρS,i = ρi − ρI,i . The steady state (i.e., the solution of ∂t ρI,i = 0) without diffusion is
ρ∗I,i
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
−γ 0 ( ) (β < µρi ), µ = ρi 1 − ρ−γ (β > µρ−γ i ). β i
(24)
We introduce small diffusion coefficient ε and assme that ε is so small that ρi is constant over time. We can formally write ρ∗I,i in the expansion form with ε as ρ∗I,i =
∞ ∑
(m)
(25)
εm ρI,i .
m=0
Substituting this expansion form to ∂t ρI,i = −µρ∗I,i + βργ−1 (ρi − ρ∗I,i )ρ∗I,i = 0 leads to i ∑ ∑ (m) (m ) (m) 0 = (βργi − µ) εm ρI,i − βργ−1 εm+m ρI,i ρI,i i −
∑
m,m
m
εm+1
m
∑
(m)
(26)
Lij ρI,j .
j
First, let us take the terms of O(ε0 ). These terms satisfy )2 ( (0) (0) 0 = (βργi − µ) ρI,i − β ρI,i (27) ) ( (0) which has the solutions ρI,i = 0 and ρi 1 − βµ ρ−γ . They i
recover the non-perturbed solution given by Eq. (24). Second, the terms of O(ε) satisfy ∑ (1) (0) (1) (0) 0 = (βργi − µ) ρI,i − βργ−1 · 2ρ ρ − Lij ρI,j . (28) i I,i I,i j
µρ−γ i
If β < (i.e., node i is at the extinction equilibrium (0) with ε = 0), ρI,i = 0 and ∑ 1 (1) (0) Lij ρI,j . (29) ρI,i = γ βρi − µ j
If β > µρ−γ (i.e., node i ( i is at )the endemic equilibrium (0) and with ε = 0), ρI,i = ρi 1 − βµ ρ−γ i ∑ 1 (1) (0) Lij ρI,j . (30) ρI,i = − γ βρi − µ j We can extend the same derivation to the terms of O(εm ) with m ≥ 2.
In the following, we demonstrate the results of the perturbation analysis when the network is the wheel graph with N nodes. The wheel graph with N nodes consists of a cycle with N − 1 nodes and a single hub connected to the rest of the nodes. Nodes on the cycle have degree 3 and the node at the center has degree N − 1. The average degree of the N -wheel graph is equal to k = 4(N − 1)/N . We refer to the nodes on the cycle as fringe nodes and to the node at the center as hub node. We take the wheel graph as an example because L comprises only two type of connection, i.e., (hub-fringe) and (fringe-fringe) and we can explicitly (m) write down ρI,i . At the initial state, the allocation of particles is set to ρfringe = 3N ρ/4(N − 1) and ρhub = N ρ/4. When β satisfy −γ µρ−γ hub < β < µρfringe ,
(31)
the hub is at endemic equilibrium with ε = 0 whereas the fringe nodes are at the extinction equilibrium. Under this parameter setting, the O(ε) term of perturbed solution for a fringe node is given by 144
0.8 fraction of infected particles
IFAC CHAOS 2015 144 August 26-28, 2015. Tokyo, Japan
ρI,hub ρI,fringe
0.6 0.4 0.2 0
0
0.1
0.2
0.3
0.4
diffusion coefficient ε
Fig. 1. Fractions of I particles at the hub node (red squares) and a fringe node (black circles) as a function of diffusion coefficient ε. The fraction of I particles at fringe nodes are averaged over all N − 1 fringe nodes. The solid and dashed lines represent the theoretical solutions up to the order of O(ε5 ), for the hub and fringe nodes, respectively. We set (N, ρ, β, γ) = (100, 100, 0.06, 0.5). )γ ( ( )−1 Nρ 3N ρ (1) ρI,fringe = − β −µ 4(N − 1) 4(N − 1) ( ( )−γ ) µ Nρ > 0. (32) × 1− β 4
In addition, the O(ε) term of perturbed solution for the hub node is given by ( )−γ+1 1 Nρ (1) ρI,hub = − < 0. (33) β 4 Therefore, in the order of O(ε), the number of I particles at fringe nodes increases due to diffusion whereas that at the hub node decreases. In Fig. 1, we confirm the accuracy of the perturbation anal∑5 k (k) ysis by comparing theoretical solutions k=0 ε ρI,fringe ∑5 (k) and k=0 εk ρI,hub with ρI,fringe and ρI,hub obtained from numerical simulations. As expected from the analysis described above, ρI,fringe increases with ε for small ε > 0 and ρI,hub decreases. The theoretical solutions agree well with the numerical results in the range of 0 < ε < 0.3. 3.3 Transient dynamics As we observed in Sec. 3.1, only the principal eigenvector of L matter as long as we focus on the properties at the steady state such as endemic threshold. However, other eigenvalues should have impact on short-term properties, that is, transient dynamics from initial to steady states. Because the metapoplation SIS model conserves the total number of particles and particles move between nodes in a random manner, ρi converges to the stationary distribution {ki / kρ}i=1,2,...,N with time scale λ2 · max[DS , DI ], where λ2 is the smallest nonzero eigenvalue of L. It should be noted that λ2 is real for arbitrary undirected networks. If diffusion processes are sufficiently faster than infections, i.e., λ2 · max[DS , DI ] βi , the number of particles converges to ki / kρ before interaction occurs and there is no characteristic transient dynamics. Otherwise, if λ2 ·
IFAC CHAOS 2015 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
(a)
(b)
8000 6000
6000
(N-1)ρI, fringe ρhub (N-1)ρfringe
4000
4000
2000 0
2000
0
20
40
60
80
100
0
time
(c)
8000
6000
6000
4000
4000
2000
2000
0
20
40
0
20
40
60
80
100
0
0
60
80
100
60
80
100
time
(d)
8000
0
REFERENCES
8000
ρI,hub
145
20
time
40 time
Fig. 2. Transient dynamics on the wheel graph with (a) D = 1 and (b,c,d) D = 0.1. Solid lines correspond to ρI,hub (red) and (N − 1)ρI,fringe (black), respectively. Dashed lines correspond to ρhub (magenta) and (N − 1)ρfringe (blue), respectively. For (a,b,c), we set other parameters (N, ρ, γ, β) = (100, 100, 0.5, 0.06). Only for (d), we set (γ, β) = (0, 1.5) while other parameters are the same as those in (a,b,c). The initial number of particles on the hub are (a,b,d) ρhub (t = 0) = N ρ/2 and (c) ρhub (t = 0) = N ρ/5. max[DS , DI ] βi , infection occurs before particles reach to the stationary distribution. In Fig. 2, time series of the transient dynamics are shown for four different parameter settings. We use the wheel graph with N = 100 again as an example. Here we assume D = DS = DI for the sake of simplicity. When D is large (see Fig. 2(a)), ρi rapidly converges to the stationary value ρ∗i and ρI,i (t) monotonically reaches to the stationary state with small fluctuation. In contrast, when D are comparable with βi (see Fig. 2(b)), infection occurs before ρi converges to the stationary distribution. In the case shown in Fig. 2(b), since ρhub (t = 0) = N ρ/2 > ρ∗hub = N ρ/4, ρI,hub initially increases rapidly and then turn to decrease as ρhub decreases due to the diffusion. The overshoot behavior of ρI,hub shown in Fig. 2(b) depends on initial condition ρhub (t = 0). As shown in Fig. 2(c), if ρhub (t = 0) < ρ∗hub , ρI,hub monotonically reaches to the steady state even if D is the same as the value in Fig. 2(b). It also should be noted that these varieties of transient dynamics observed in Figs. 2(a),(b),(c) are realized essentially due the population-dependent infection rate βi = βργi with γ = 0. If γ = 0 and the contact rate is independent from population size, there is no characteristic behavior of ρI,hub even when ρhub (t = 0) > ρ∗hub and D is sufficiently small (see Fig. 2(d)). 4. CONCLUSIONS We considered the metapopulation SIS model with the nonlinear population-dependent contact rate on arbitrary networks. We showed the impacts of network topology and the contact rate on endemic threshold, infection size at the stationary state, and transient dynamics. 145
R. M. Anderson and R. M. May, Directly transmitted infections diseases: control by vaccination. Science 215:1053–1060, 1982. R. M. Anderson, and R. M. May, Infectious diseases of humans: Dynamics and control. Oxford University Press. 1991. F. R. K. Chung, Spectral Graph Theory. American Mathematical Society, 1997. V. Colizza, R. Pastor-Satorras, and A. Vespignani Reaction-diffusion processes and metapopulation models in heterogeneous networks. Nat. Phys. 3:276–282, 2007. V. Colizza and A. Vespignani, Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations. J. Theor. Biol. 251:450–67, 2008. S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, k-Core organization of complex networks. Phys. Rev. Lett. 96:040601, 2006. A. Iggidr, G. Sallet, and B. Tsanou, Global stability analysis of a metapopulation SIS epidemic model. Math. Popul. Stud. 19:115–129, 2012. D. Juher, J. Ripoll, and J. Salda˜ na, Analysis and Monte Carlo simulations of a model for the spread of infectious diseases in heterogeneous metapopulations. Phys. Rev. E 80:041920, 2009. D. Juher and V. Ma˜ nosa, Spectral properties of the connectivity matrix and the SIS-epidemic threshold for mid-size metapopulations. Math. Model. Nat. 9:108– 120, 2014. H. Lund, L. Lizana, and I. Simonsen, Effects of citysize heterogeneity on epidemic spreading in a metapopulation: a reaction-diffusion approach. J. Stat. Phys. 151:367–382, 2013. N. Masuda, Effects of diffusion rates on epidemic spreads in metapopulation networks. New J. Phys. 12:093009, 2010. H. McCallum, N. Barlow, and J. Hone, How should pathogen transmission be modelled? Trends Ecol. Evol. 16:295–300, 2001. J. Salda˜ na, Continuous-time formulation of reactiondiffusion processes on heterogeneous metapopulations. Phys. Rev. E 78:012902, 2008. J. Salda˜ na, Modelling the spread of infectious diseases in complex metapopulations. Math. Model. Nat. 5:22–37, 2010. L. P. Sanders, B. S¨oderberg, D. Brockmann, and T. Ambj¨ornsson, Perturbative solution to susceptibleinfected-susceptible epidemics on networks. Phys. Rev. E 88:032713, 2013. M. Schl¨apfer, L. M. A. Bettencourt, S. Grauwin, M. Raschke, R. Claxton, Z. Smoreda, G. B. West, and C. Ratti, The scaling of human interactions with city size. J. R. Soc. Interface 11:20130789, 2014. M. J. Smith, S. Telfer, E. R. Kallio, S. Burthe, A. R. Cook, X. Lambin, and M. Begon, Host-pathogen time series data in wildlife support a transmission function between density and frequency dependence. Proc. Natl. Acad. Sci. U.S.A. 106:7905–7909, 2009.
ScienceDirect IFAC-PapersOnLine 48-18 (2015) 141–145
Analysis of metapopulation epidemic Analysis of metapopulation epidemic Analysis of metapopulation epidemic process on arbitrary networks process on arbitrary networks process on arbitrary networks
Taro Takaguchi ∗∗ Renaud Lambiotte ∗∗ ∗∗ Taro Takaguchi ∗∗ Renaud Lambiotte ∗∗ Taro Taro Takaguchi Takaguchi Renaud Renaud Lambiotte Lambiotte ∗∗ ∗ National Institute of Informatics, Tokyo 101-8430, Japan (e-mail: ∗ ∗ Institute of Informatics, Tokyo 101-8430, Japan (e-mail: ∗ National National Tokyo [email protected]) National Institute Institute of of tInformatics, Informatics, Tokyo 101-8430, 101-8430, Japan Japan (e-mail: (e-mail: t [email protected]) ∗∗ t [email protected]) [email protected]) Department of Mathematics and naXys, University of Namur, 5000 t ∗∗ ∗∗ Department of Mathematics and naXys, University of Namur, 5000 ∗∗ Department Mathematics naXys, Namur, of Belgium (e-mail:and [email protected]) Department of Mathematics and naXys, University University of of Namur, Namur, 5000 5000 Namur, Belgium (e-mail: [email protected]) Namur, Belgium (e-mail: [email protected]) Namur, Belgium (e-mail: [email protected]) Abstract: Epidemic process on networks is considered. The system is discribed as a metapopAbstract: Epidemic process on networks issubpopulation considered. The system is discribed as a metapopAbstract: Epidemic process on networks system is discribed as metapopulation network in which a node (e.g., city or and Abstract: Epidemic process on represents networks is is considered. considered. The The system is school) discribed asisaa connected metapopulation network in which a node represents subpopulation (e.g., city or school) and is connected ulation network in which a node represents subpopulation (e.g., city or school) and is connected to other nodes via undirected links. Particles represent the subject of infection (e.g., individuals) ulation network in which a node represents subpopulation (e.g., city or school) and is connected to nodes via undirected links. represent the subject of individuals) to other nodes viaeach undirected links. Particles Particles represent the from subject of infection infection (e.g., individuals) andother interact with other within nodes while migrating nodes to nodes(e.g., in the manner of to other nodes via undirected links. Particles represent the subject of infection (e.g., individuals) and interact with each other within nodes while migrating from nodes to nodes in the mannersize of and interact with each other within nodes while migrating from nodes to nodes in the random diffusion. The nonlinear dependence of contact ratefrom within a node on itsinpopulation and interact with each other within nodes while migrating nodes to nodes the manner manner of of random diffusion. The nonlinear dependence of contact rate within a node on its population size random diffusion. The nonlinear dependence of contact rate within a node on its population size is introduced, according to the recent findingofbased onrate emprical The impacts random diffusion. The nonlinear dependence contact withinphone-call a node on data. its population size is according to the finding based on emprical The impacts is introduced, according to recent finding on phone-call data. The impacts of introduced, the nonlinear dependence arerecent investigated for three aspects ofphone-call epidemic data. process: is introduced, according to the the recent finding based based on emprical emprical phone-call data. Theepidemic impacts of the nonlinear dependence are investigated for three aspects of epidemic process: epidemic of the nonlinear nonlinear dependence are investigated investigated fortransient three aspects aspects of epidemic epidemic process: process: epidemic epidemic threshold, infection size at stationary state, and dynamics. of the dependence are for three of threshold, infection size at stationary state, and transient dynamics. threshold, threshold, infection infection size size at at stationary stationary state, state, and and transient transient dynamics. dynamics. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Complex systems; Networks; Differential equations; Eigenmode analysis; Monte Keywords: Complex Keywords: Complex systems; systems; Networks; Networks; Differential Differential equations; equations; Eigenmode Eigenmode analysis; analysis; Monte Monte Carlo simulation Keywords: Complex systems; Networks; Differential equations; Eigenmode analysis; Monte Carlo simulation Carlo simulation Carlo simulation 1. INTRODUCTION recently observed in various emprical studies (Smith, et al. 1. recently observed various emprical (Smith, et al. 1. INTRODUCTION INTRODUCTION recently observed in emprical studies (Smith, et (2009); Schl¨ apfer,in al. (2014)). To studies get better insight of 1. INTRODUCTION recently observed inetvarious various emprical studies (Smith, et al. al. (2009); Schl¨ a pfer, et al. (2014)). To get better insight of (2009); aapfer, et (2014)). To better of process will insight introduce (2009); Schl¨ Schl¨ pfer,under et al. al.realistic (2014)).settings, To get get we better insight of Understanding and controlling spreading process of infec- epidemic epidemic process under realistic settings, we will introduce Understanding and controlling spreading process of infecepidemic process under realistic settings, we will introduce Understanding and controlling spreading process of infecthese two factors into standard model and theoretically epidemic process under realistic settings, we will introduce tious diseases are of much interstspreading in the research fields such these two factors into standard model and theoretically Understanding and controlling process of infectious diseases of interst in research fields such two into tious diseases are are public of much muchhealth, interstand in the the research fields (Ansuch these investigate its properties. these two factors factors into standard standard model model and and theoretically theoretically as epidemiology, network science tious diseases are of much interst in the research fields such investigate its properties. as epidemiology, public health, and network science (Aninvestigate its properties. as epidemiology, public health, health, and network network science (An- investigate its properties. derson & May (1991)). A standard framework to model as epidemiology, public and science (Anderson & (1991)). standard to 2. MODEL derson & May May (1991)). A A standard isframework framework to model model epidemics in structrured population so-called metapopderson & May (1991)). A standard framework to model 2. MODEL MODEL epidemics in structrured population is so-called metapop2. epidemics in structrured population is so-called metapop2. MODEL ulation model, in which nodes represent subpopulations epidemics in structrured population is so-called metapopulation model, in which nodes represent subpopulations suppose that the network consists of undirected and ulation model, in which which nodes represent subpopulations and links connectes nodenodes pairs. represent The subject of epidemic We ulation model, in subpopulations We suppose links. that the consists of undirected and links connectes pairs. subject of epidemic We network consists of and We network also suppose that the networkand is and links connectes node node pairs. The The subject of with epidemic We suppose suppose that that the the network consists of undirected undirected and process is represented by particles which interact each unweighted and links connectes node pairs. The subject of epidemic unweighted links. We also suppose that the network is process is represented by particles which interact with each unweighted links. We also suppose that the network is connected (i.e., one can move between any node pairs via process is represented by particles which interact with each unweighted links. We also suppose that the network is other within a node and between in awith random process is represented by move particles whichnodes interact each connected (i.e., one can move between any node pairs via other within a node and move between nodes in a random connected (i.e., one can move between any node pairs via its links) and simple (i.e., no self-loops or multiple links). other within a node and move between nodes in a random connected (i.e., one can move between any node pairs via manner. For example, in airport networks, nodes are cities other within a node and move between nodes in a random its links) and simple (i.e., no self-loops or multiple links). manner. For example, airport nodes are and (i.e., multiple links). Letlinks) us define N nodes inno theself-loops networkor np particles manner. For and example, in airport networks, networks, nodes are cities cities its links) and simple simple (i.e., no self-loops orand multiple links). or countries links in represent connections between air- its manner. For example, in airport networks, nodes are cities Let us define N nodes in the network and n particles p or countries and links represent connections between airLet us define N nodes in the network and n particles p on it. Therefore, the average number of particles on a or countries and links represent connections between airLet us define N nodes in the network and np particles ports via flights. In this example, a particle between represents or countries and links represent connections air-a on it. Therefore, the average number of particles on ports via flights. In this example, a particle represents a on it. Therefore, the average number of particles on aaa node is given by ρ ≡ n /N . We define the rules governing ports via flights. In this example, a particle represents a on it. Therefore, the average number of particles on p passenger traveling cities by air. Therepresents theoreticala node is given by ρ ≡ n /N . We define the rules governing ports via flights. In between this example, a particle p passenger traveling between cities by air. The theoretical node is given by ρ ≡ n /N . We define the rules governing the states of by particles according to susceptible-infectedpassenger traveling betweenmodels cities by by air. The theoretical is given ρ ≡ npp /N . We define the rules governing studies of traveling metapopulation seek forThe thetheoretical impact of node passenger between cities air. the states of particles according to susceptible-infectedstudies metapopulation models seek the impact the particles according to (SIS) dynamics which has been investigated in studies oftopology metapopulation models seek for for on the large-scale impact of of susceptible the states states of of particles according to susceptible-infectedsusceptible-infectednetworkof and model parameters studies of metapopulation models seek for the impact of susceptible (SIS) which has investigated in network topology susceptible (SIS) dynamics dynamics which has been beenSalda˜ investigated in previous studies (Colizza, et al. (2007); na (2008); network topology and and model model parameters parameters on on large-scale large-scale susceptible (SIS) dynamics which has been investigated in outbreaks. network topology and model parameters on large-scale previous studies (Colizza, et al. (2007); Salda˜ n a (2008); outbreaks. previous&studies studies (Colizza, et Juher, al. (2007); (2007); Salda˜ naa Salda˜ (2008); Colizza Vespignani (2008); et al. (2009); n a outbreaks. previous (Colizza, et al. Salda˜ n (2008); outbreaks. & Vespignani (2008); Juher, et al. (2009); Salda˜ n aa In this paper, we take into accout two new aspects: Colizza Colizza & Vespignani (2008); Juher, et al. (2009); Salda˜ n (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. Colizza & Vespignani (2008); Juher, et al. (2009); Salda˜ n a In this take into accout new aspects: Masuda (2010); Iggidr, et al. (2012); Lund, et al. In this paper, paper, wefree take into accout two two and newnonlinear aspects: (2010); network topologywe from approximation (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. In this paper, we take into accout two new aspects: (2013); Juher & Ma˜ n osa (2014)). In SIS dynamics, the (2010); Masuda (2010); Iggidr, et al. (2012); Lund, et al. network topology free from approximation and nonlinear & n (2014)). SIS dynamics, network topology free free contact from approximation approximation and recent nonlinear population-dependent rate. First, while re- (2013); (2013); Juher & Ma˜ Ma˜ nosa osa to (2014)). In SIS compartments: dynamics, the the network topology from and nonlinear state of Juher a particle belongs either In of two (2013); Juher & Ma˜ n osa (2014)). In SIS dynamics, the population-dependent contact rate. First, while recent restate of a particle belongs to either of two compartments: population-dependent contact rate. First, while recent research on metapopulation models has unveiled the impact state of a particle belongs to either of two compartments: population-dependent contact rate. First, while recent re- state susceptible (S) and infected (I). The states of particles of a particle belongs to either of two compartments: search on models has impact (S) and infected (I). states of particles search on metapopulation metapopulation models has unveiled unveiled the impact of heterogeneous distribution of degree (i.e., thethe number of susceptible susceptible (S) infected The states of search on metapopulation models has unveiled the impact evolve through the infection andThe recovery as susceptible (S) and and infected (I). (I). The states processes of particles particles of heterogeneous distribution of degree (i.e., the number through the infection and recovery processes as of heterogeneous distribution of degree (i.e., the number ofa evolve links connected todistribution a node) (Colizza, et (i.e., al. (2007); Salda˜ nof evolve through the infection and recovery processes as of heterogeneous of degree the number of follows: evolve through the infection and recovery processes as links connected to a node) (Colizza, et al. (2007); Salda˜ n a { follows: links connected to a node) (Colizza, et al. (2007); Salda˜ n a (2008); Colizza & (2008);etJuher, et al.Salda˜ (2009); follows: links connected to Vespignani a node) (Colizza, al. (2007); na follows: { S + I → I + I with rate β , i { (2008); Colizza & Vespignani (2008); Juher, et al. (2009); {S (2008); Colizza &Juher Vespignani (2008); Juher, the et al. al.majority (2009); + I → II + I with rate β , (1) Salda˜ naColizza (2010);& & Ma˜ nosa (2014)), (2008); Vespignani (2008); Juher, et (2009); + with rate rate µ. β S→ + IIS→ →I+ + II with βiii ,, (1) IS Salda˜ n & n (2014)), the (1) Salda˜ naaa (2010); (2010); Juher & Ma˜ Ma˜ nosa osatype (2014)), the majority majority of them assumedJuher the mean-filed approximation of (1) Salda˜ n (2010); Juher & Ma˜ n osa (2014)), the majority I → S with rate µ. II → S with rate µ. of them assumed the mean-filed type approximation of → S with rate µ. of them assumed the mean-filed type approximation of varies for different We assume that the infection rate β network equate alltype nodes with the same of them topology assumed which the mean-filed approximation of We assume that the infection rate βi varies for different network topology which equate all with same varies for that different We assume assume that the infection rate β βiii be network topologywhile which equate all nodes nodesstudies with the the same nodes i (1 that ≤ i the ≤ infection N ). It should noted the varies for different We rate degree. Second, most of previous presume network topology which equate all nodes with the same nodes i (1 ≤ i ≤ N ). It should beat noted noted thatnode. the degree. Second, while most of previous studies presume nodes i (1 ≤ i ≤ N ). It should that degree. Second, while most of of particles previous within studiesapresume presume infection occurs particlesbe the same i (1 ≤ ionly ≤ between N ). It should be noted that the the that theSecond, contact while rate between node is nodes degree. most previous studies infection occurs only between particles at the same node. that the contact rate between particles within a node is infection occurs only between particles at the same node. that the contact contact rate proportional between particles particles within node at is infection In addition to theonly reactions, particles move between occurs between particles at the same nodes node. irrelevant or linearly to the population that the rate between within aa node is addition to the reactions, particles between nodes irrelevant or proportional to the at In addition to reactions, particles move between nodes irrelevant or linearly linearly proportional to population the population population at In in manner of simple diffusion withmove diffusion coefficient In the addition to the the reactions, particles move between nodes the node, nonlinear dependence on the size are irrelevant or linearly proportional to the population at in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are D and D for S and I particles, respectively. In general, in the manner of simple diffusion with diffusion coefficient the node, nonlinear dependence on the population size are DS and DI for S and I particles, respectively. In general, S II for S and D and II particles, respectively. general, This work was partly supported by Bilateral Joint Research S there is aD value of β (denoted by βc In and called D Dthreshold particles, respectively. In general, S and I for S and there is a threshold value of β (denoted by β and called c This work was partly supported by Bilateral Joint Research there is a threshold value of β (denoted by β and endemic If β of>β β of This work work wasJSPS, partly supported by Bilateral Bilateral Joint Research Research there is a threshold). threshold value (denoted by βcc fraction and called called Projects between Japan, and F.R.S.-FNRS, Belgium. T.T. was c a positive This was partly supported by Joint endemic threshold). If β > β a positive fraction of Projects between JSPS, Japan, and F.R.S.-FNRS, Belgium. T.T. was If β β fraction of Projects JSPS, Japan, and T.T. supported by JST, ERATO, Kawarabayashi LargeBelgium. Graph Project. Iendemic particlesthreshold). remains at state (called endemic threshold). If the β > >stationary βccc aa positive positive fraction the of Projects between between JSPS, Japan, and F.R.S.-FNRS, F.R.S.-FNRS, Belgium. T.T. was was supported by JST, ERATO, Kawarabayashi Large Graph Project. I particles remains at the stationary state (called the supported by by JST, JST, ERATO, ERATO, Kawarabayashi Kawarabayashi Large Large Graph Graph Project. Project. particles remains remains at at the the stationary stationary state state (called (called the the supported II particles Copyright 2015 IFAC 141 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © 2015, IFAC (International Federation of Automatic Control) Copyright 2015 IFAC 141 Copyright © 2015 IFAC 141 Peer review© of International Federation of Automatic Copyright ©under 2015 responsibility IFAC 141Control. 10.1016/j.ifacol.2015.11.026
IFAC CHAOS 2015 142 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
endemic equilibrium). Otherwise, if β < βc , the fraction of I particles is zero in the long-term limit (called the extinction equilibrium). Dependence of βc on network structure is one of the main concerns in this study. We employ the continuous-time formulation of metapopulation SIS dynamics (Salda˜ na (2008)) which assumes the simultaneity of interaction between particles and diffusion of particles. Taking advantage of this formulation, we can describe the average behavior of the system by the deterministic reaction-diffusion equation (Masuda (2010)) as ∑ β ∂t ρS,i = µρI,i − i ρS,i ρI,i − DS Lij ρS,j , ρi j (2) ∑ βi ∂ ρ = −µρ + ρ ρ − D L ρ , t I,i I,i S,i I,i I ij I,j ρi j
where ρS,i and ρI,i represent the number of S and I particles at node i, respectively. The total number of particles at node i is given by ρi ≡ ρS,i + ρI,i . Matrix L denotes the random-walk Laplacian whose elements obey { 1 (i = j), Lij = (3) −Aij /kj (i = j),
where Aij is the (i, j)-element of adjacency matrix A and takes a binary value, i.e., unity if nodes i and j are adjacent and zero otherwise. Degree ki of node i is the total number of links associated to the node. Since we suppose the network is undirected, A is a symmetric matrix but L is asymmetric.
We assume that the infection rate βi at node i depends on ρi as follows: (4) βi = βργi (γ ≥ 0). This nonlinear dependence is motivated by recent empirical findings (Smith, et al. (2009); Schl¨ apfer, et al. (2014)) described in Sec. 1. Equation (4) was previously mentioned in (Anderson & May (1982); McCallum, et al. (2001)) and was recently considered in combination with meanfield networks (Juher & Ma˜ nosa (2014)). Equation (4) interpolates the two standard assumption: the so-called frequency interaction recovers if γ = 0 and the so-called mass interaction (Colizza, et al. (2007); Masuda (2010)) recovers if γ = 1. 3. RESULTS In this section, we report following three aspects of the metapopulation SIS model defined in Sec. 2. First, in Sec. 3.1, we investigate the endemic threshold βc . We cannot derive the exact βc value in general case, however, we derive its upper bound which appears to be fairly tight according to numerical simulations on some network models. Second, in Sec. 3.2, we investigate the infection size at the stationary state. In particular, we assume that diffusion of particles as the perturbation to the steady state without diffusion, and derive the infection size as a perturbed solution of the reaction-diffusion equation (Eq. (2)). Third, in Sec. 3.3, we investigate the impact of the population-dependent infection rate (Eq. (4)) on the transient process in which the configuration of particles converges from the initial to the stationary distribution. We show that the the population-dependent infection rate 142
induces nonmonotonical behavior in the convergene of population. The metapopulation SIS model (Eq. (2)) contains four parameters related to its time scale: β, µ, DS , and DI . Without loss of generality, we set µ = 1 and DS = DI for all the numerical simulaitons reported in this section. 3.1 Endemic threshold The results reported in this section are essentially in parallel with previous studies (Salda˜ na (2008); Juher, et al. (2009); Salda˜ na (2010); Masuda (2010); Lund, et al. (2013); Juher & Ma˜ nosa (2014)). Our results differ from these studies because we consider arbitrary networks free from the heterogeneous mean-field (HMF) approximation of networks, while these studies are based on this approximation, except for (Masuda (2010)). With the HMF approximation, the nodes with the same ki are supposed to behave in the same way on average (Salda˜ na (2008, 2010)). To be more precise, a pair of nodes with degree k and k are stochastically connected to each other with probability kk P (k)P (k )/k2 , where k is the average degree of the network, in an uncorrelated network. In addition, the population of particles at a node with k is determined by (ρS,k , ρI,k ) regardless of the nodes’ indices i. We are interested in endemic threshold βc . Since the extinction equilibrium ρ∗I,k = 0 is the stable if β < βc , βc is given as the β value with which the extinction equilibrium becomes unstable. By analysing the eigenvalue of the Jacobian at the equilibrium, one obtain the upper bound of βc (Salda˜ na (2008)) as )−γ ( kmax ρ β HMF ≡ (µ + DI ) (5) k In words, βc < β HMF holds true under the HMF approximation.
In the following, we prove that the upper bound β HMF also holds true for arbitrary networks beyond the HMF approximation. As we described for the derivation under the HMF approximation, we consider the stability of the extinction equilibrium (ρ∗S,i , ρ∗I,i ) = (ki /kρ, 0) (1 ≤ i ≤ N ) of the model (Eq. (2)). The allocation ρ∗S,i = ki /kρ streightforwardly follows as the stationary distribution of diffusion of single-type particles. The linear stability of the extinction equilibrium is determined by the Jacobian matrix at the equilibrium: ( (1) (3) ) J J J= (6) O J (2) All submatrices J (1) , J (2) , J (3) , and O are with size N ×N . Matrix O is zero matrix. The element of J (1) and J (2) is defined by (1)
Jij = −DS Lij , (7) )γ ) ( ( ki (2) ρ − µ δij − DI Lij . Jij = β (8) k When the largest eigenvalue of J is positive, the extinction equilibrium is linearly unstable. Since the eigenvalues of L are in 0 ≤ λ(L) ≤ 2 (Chung (1997)), the eigenvalues of J (1)
IFAC CHAOS 2015 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
are in −2DS ≤ λ(J (1) ) ≤ 0 and not positive regardless of parameter values. Therefore, to derive the condition for the endemic equilibrium, we investigate the largest eigenvalue of J (2) in the following. There are two special cases and we can exactly derive βc . As the first case, when γ = 0, the element of J (2) is equal to (2) Jij = (β − µ) δij − DI Lij . (9)
Therefore, its largest eigenvalue λmax (J (2) ) is given by (β − µ) and it leads to βc = µ. As the second case, when the network is regular and ki = k (1 ≤ i ≤ N ), the element of J (2) is equal to (2) Jij = (βργ − µ) δij − DI Lij . λmax (J (2) ) is given by (βργ − µ) and
(10)
Therefore, it leads to βc = µρ−γ . It should be noted that DI (> 0) is irrelevant to βc for these two cases. We now move to general cases with γ > 0 on a nonregular network. For the convenience for calculation, we symmetrize J (2) via the following similarity transformation (Masuda (2010)): ) ( 1 1 1 J (2) Jˆ(2) ≡diag √ , √ , . . . , √ k1 k2 kN (√ √ √ ) (11) × diag k1 , k2 , . . . , kN , where diag(a1 , a2 , . . . , aN ) represents a diagonal matrix in which the (i, i)-element is given by ai . With this transformation, the element of J (2) is converted into )γ ) ( ( Aij ki (2) ˆ ρ − µ − DI δij + DI √ . (12) Jij = β k ki kj
If λmax (Jˆ(2) ) is positive, the extinction equilibrium is unstable and the endemic equilibrium arises. We can represent λmax (Jˆ(2) ) by using the Rayleigh quotient as x Jˆ(2) x . (13) λmax (Jˆ(2) ) = max x x |x|=1
We show that βc < β HMF holds true for arbitrary networks as follows. Let us define x = (0, 0, . . . , 1, . . . , 0) in which only the i-th element is equal to unity and the rest equal to zero. With this x, )γ ( ki ρ − µ − DI . (14) x Jˆ(2) x = β k Therefore, )γ ( ki ρ − µ − DI ≤ λmax (Jˆ(2) ), (15) max β 1≤i≤N k and we obtain )−γ ( kmax β > (µ + DI ) ρ = β HMF . (16) k Equation (16) states that, if β > β HMF , λmax (Jˆ(2) ) is positive and the endemic threshold arises. In other words, βc < β HMF holds true. We show an improved upper bound of βc . We observed that the upper bound β HMF is achieved via approximation of the principal eigenvector of Jˆ(2) with x = (0, 0, . . . , 1, . . . , 0) in which only the element corresponding to a single node with kmax is equal to unity. Therefore, 143
143
we elaborate this approximation so that the elements of x corresponding to the nodes with the largest degree have nonzero values. Let us define x such that (ki < kc ), 0√ (17) xi = ki (ki ≥ kc ), Kc takes an integer value in [1, kmax ] where free parameter k c ∑ and Kc ≡ j:kj ≥kc kj is a normalization constant. With this choice of x, )γ ( ) ( ∑ 2Ec ki ki ˆ(2) ρ , β − µ − DI 1 − x J x= k Kc Kc i:ki ≥kc
(18) )γ ( ) 2Ec kc ρ − µ − DI 1 − , (19) ≥β k Kc where Ec represents the total number of links between nodes with ki ≥ kc . Therefore, x Jˆ(2) x λmax (Jˆ(2) ) = max x x |x|=1 )γ [ ( ] 2Ec kc ρ + DI − (µ + DI ) . ≥ max β kc k Kc (20) (
A sufficient condition for the endemic equilibrium is that there exists kc which satisfies )γ ( ) ( 2Ec kc ρ − 1− DI − µ > 0. (21) β k Kc In terms of β, this is equivalent to {[ )−γ } ( ) ]( 2Ec kc ρ = β. (22) DI β > min µ + 1 − kc Kc k This sufficient condition improves the previous one (given by Eq. (16)) because β ≤ β HMF and the equality holds true if kc = kmax . 3.2 Infection size at the steady state In the previous section, we discussed endemic threshold βc . Another important property of metapopulation SIS dynamics is the infection size at the stationary state when β > βc . In particular, we want to derive the infection size when diffusion coefficient DS = DI = ε is very small, in order to understand the impact of diffusion process on epidemic spreading. To this aim, we approximate the solution of Eq. (2) with ε > 0 by perturbing its solution with ε = 0 with small ε. The basic idea of the analysis reported in this section was originally proposed in (Sanders, et al. (2013)). We focus on the infection size at the stationary state, while the previous study (Sanders, et al. (2013)) considered its time evolution from the initial state. We assume that the initial number of particles at node i is proportional to ki , that is, ρi = ki ρ/k. When ε = 0, the time evolution of ρI,i is governed by ∂t ρI,i = −µρI,i + βργ−1 (ρi − ρI,i )ρI,i , i
(23)
where we substitute ρS,i = ρi − ρI,i . The steady state (i.e., the solution of ∂t ρI,i = 0) without diffusion is
ρ∗I,i
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
−γ 0 ( ) (β < µρi ), µ = ρi 1 − ρ−γ (β > µρ−γ i ). β i
(24)
We introduce small diffusion coefficient ε and assme that ε is so small that ρi is constant over time. We can formally write ρ∗I,i in the expansion form with ε as ρ∗I,i =
∞ ∑
(m)
(25)
εm ρI,i .
m=0
Substituting this expansion form to ∂t ρI,i = −µρ∗I,i + βργ−1 (ρi − ρ∗I,i )ρ∗I,i = 0 leads to i ∑ ∑ (m) (m ) (m) 0 = (βργi − µ) εm ρI,i − βργ−1 εm+m ρI,i ρI,i i −
∑
m,m
m
εm+1
m
∑
(m)
(26)
Lij ρI,j .
j
First, let us take the terms of O(ε0 ). These terms satisfy )2 ( (0) (0) 0 = (βργi − µ) ρI,i − β ρI,i (27) ) ( (0) which has the solutions ρI,i = 0 and ρi 1 − βµ ρ−γ . They i
recover the non-perturbed solution given by Eq. (24). Second, the terms of O(ε) satisfy ∑ (1) (0) (1) (0) 0 = (βργi − µ) ρI,i − βργ−1 · 2ρ ρ − Lij ρI,j . (28) i I,i I,i j
µρ−γ i
If β < (i.e., node i is at the extinction equilibrium (0) with ε = 0), ρI,i = 0 and ∑ 1 (1) (0) Lij ρI,j . (29) ρI,i = γ βρi − µ j
If β > µρ−γ (i.e., node i ( i is at )the endemic equilibrium (0) and with ε = 0), ρI,i = ρi 1 − βµ ρ−γ i ∑ 1 (1) (0) Lij ρI,j . (30) ρI,i = − γ βρi − µ j We can extend the same derivation to the terms of O(εm ) with m ≥ 2.
In the following, we demonstrate the results of the perturbation analysis when the network is the wheel graph with N nodes. The wheel graph with N nodes consists of a cycle with N − 1 nodes and a single hub connected to the rest of the nodes. Nodes on the cycle have degree 3 and the node at the center has degree N − 1. The average degree of the N -wheel graph is equal to k = 4(N − 1)/N . We refer to the nodes on the cycle as fringe nodes and to the node at the center as hub node. We take the wheel graph as an example because L comprises only two type of connection, i.e., (hub-fringe) and (fringe-fringe) and we can explicitly (m) write down ρI,i . At the initial state, the allocation of particles is set to ρfringe = 3N ρ/4(N − 1) and ρhub = N ρ/4. When β satisfy −γ µρ−γ hub < β < µρfringe ,
(31)
the hub is at endemic equilibrium with ε = 0 whereas the fringe nodes are at the extinction equilibrium. Under this parameter setting, the O(ε) term of perturbed solution for a fringe node is given by 144
0.8 fraction of infected particles
IFAC CHAOS 2015 144 August 26-28, 2015. Tokyo, Japan
ρI,hub ρI,fringe
0.6 0.4 0.2 0
0
0.1
0.2
0.3
0.4
diffusion coefficient ε
Fig. 1. Fractions of I particles at the hub node (red squares) and a fringe node (black circles) as a function of diffusion coefficient ε. The fraction of I particles at fringe nodes are averaged over all N − 1 fringe nodes. The solid and dashed lines represent the theoretical solutions up to the order of O(ε5 ), for the hub and fringe nodes, respectively. We set (N, ρ, β, γ) = (100, 100, 0.06, 0.5). )γ ( ( )−1 Nρ 3N ρ (1) ρI,fringe = − β −µ 4(N − 1) 4(N − 1) ( ( )−γ ) µ Nρ > 0. (32) × 1− β 4
In addition, the O(ε) term of perturbed solution for the hub node is given by ( )−γ+1 1 Nρ (1) ρI,hub = − < 0. (33) β 4 Therefore, in the order of O(ε), the number of I particles at fringe nodes increases due to diffusion whereas that at the hub node decreases. In Fig. 1, we confirm the accuracy of the perturbation anal∑5 k (k) ysis by comparing theoretical solutions k=0 ε ρI,fringe ∑5 (k) and k=0 εk ρI,hub with ρI,fringe and ρI,hub obtained from numerical simulations. As expected from the analysis described above, ρI,fringe increases with ε for small ε > 0 and ρI,hub decreases. The theoretical solutions agree well with the numerical results in the range of 0 < ε < 0.3. 3.3 Transient dynamics As we observed in Sec. 3.1, only the principal eigenvector of L matter as long as we focus on the properties at the steady state such as endemic threshold. However, other eigenvalues should have impact on short-term properties, that is, transient dynamics from initial to steady states. Because the metapoplation SIS model conserves the total number of particles and particles move between nodes in a random manner, ρi converges to the stationary distribution {ki / kρ}i=1,2,...,N with time scale λ2 · max[DS , DI ], where λ2 is the smallest nonzero eigenvalue of L. It should be noted that λ2 is real for arbitrary undirected networks. If diffusion processes are sufficiently faster than infections, i.e., λ2 · max[DS , DI ] βi , the number of particles converges to ki / kρ before interaction occurs and there is no characteristic transient dynamics. Otherwise, if λ2 ·
IFAC CHAOS 2015 August 26-28, 2015. Tokyo, Japan
Taro Takaguchi et al. / IFAC-PapersOnLine 48-18 (2015) 141–145
(a)
(b)
8000 6000
6000
(N-1)ρI, fringe ρhub (N-1)ρfringe
4000
4000
2000 0
2000
0
20
40
60
80
100
0
time
(c)
8000
6000
6000
4000
4000
2000
2000
0
20
40
0
20
40
60
80
100
0
0
60
80
100
60
80
100
time
(d)
8000
0
REFERENCES
8000
ρI,hub
145
20
time
40 time
Fig. 2. Transient dynamics on the wheel graph with (a) D = 1 and (b,c,d) D = 0.1. Solid lines correspond to ρI,hub (red) and (N − 1)ρI,fringe (black), respectively. Dashed lines correspond to ρhub (magenta) and (N − 1)ρfringe (blue), respectively. For (a,b,c), we set other parameters (N, ρ, γ, β) = (100, 100, 0.5, 0.06). Only for (d), we set (γ, β) = (0, 1.5) while other parameters are the same as those in (a,b,c). The initial number of particles on the hub are (a,b,d) ρhub (t = 0) = N ρ/2 and (c) ρhub (t = 0) = N ρ/5. max[DS , DI ] βi , infection occurs before particles reach to the stationary distribution. In Fig. 2, time series of the transient dynamics are shown for four different parameter settings. We use the wheel graph with N = 100 again as an example. Here we assume D = DS = DI for the sake of simplicity. When D is large (see Fig. 2(a)), ρi rapidly converges to the stationary value ρ∗i and ρI,i (t) monotonically reaches to the stationary state with small fluctuation. In contrast, when D are comparable with βi (see Fig. 2(b)), infection occurs before ρi converges to the stationary distribution. In the case shown in Fig. 2(b), since ρhub (t = 0) = N ρ/2 > ρ∗hub = N ρ/4, ρI,hub initially increases rapidly and then turn to decrease as ρhub decreases due to the diffusion. The overshoot behavior of ρI,hub shown in Fig. 2(b) depends on initial condition ρhub (t = 0). As shown in Fig. 2(c), if ρhub (t = 0) < ρ∗hub , ρI,hub monotonically reaches to the steady state even if D is the same as the value in Fig. 2(b). It also should be noted that these varieties of transient dynamics observed in Figs. 2(a),(b),(c) are realized essentially due the population-dependent infection rate βi = βργi with γ = 0. If γ = 0 and the contact rate is independent from population size, there is no characteristic behavior of ρI,hub even when ρhub (t = 0) > ρ∗hub and D is sufficiently small (see Fig. 2(d)). 4. CONCLUSIONS We considered the metapopulation SIS model with the nonlinear population-dependent contact rate on arbitrary networks. We showed the impacts of network topology and the contact rate on endemic threshold, infection size at the stationary state, and transient dynamics. 145
R. M. Anderson and R. M. May, Directly transmitted infections diseases: control by vaccination. Science 215:1053–1060, 1982. R. M. Anderson, and R. M. May, Infectious diseases of humans: Dynamics and control. Oxford University Press. 1991. F. R. K. Chung, Spectral Graph Theory. American Mathematical Society, 1997. V. Colizza, R. Pastor-Satorras, and A. Vespignani Reaction-diffusion processes and metapopulation models in heterogeneous networks. Nat. Phys. 3:276–282, 2007. V. Colizza and A. Vespignani, Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations. J. Theor. Biol. 251:450–67, 2008. S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, k-Core organization of complex networks. Phys. Rev. Lett. 96:040601, 2006. A. Iggidr, G. Sallet, and B. Tsanou, Global stability analysis of a metapopulation SIS epidemic model. Math. Popul. Stud. 19:115–129, 2012. D. Juher, J. Ripoll, and J. Salda˜ na, Analysis and Monte Carlo simulations of a model for the spread of infectious diseases in heterogeneous metapopulations. Phys. Rev. E 80:041920, 2009. D. Juher and V. Ma˜ nosa, Spectral properties of the connectivity matrix and the SIS-epidemic threshold for mid-size metapopulations. Math. Model. Nat. 9:108– 120, 2014. H. Lund, L. Lizana, and I. Simonsen, Effects of citysize heterogeneity on epidemic spreading in a metapopulation: a reaction-diffusion approach. J. Stat. Phys. 151:367–382, 2013. N. Masuda, Effects of diffusion rates on epidemic spreads in metapopulation networks. New J. Phys. 12:093009, 2010. H. McCallum, N. Barlow, and J. Hone, How should pathogen transmission be modelled? Trends Ecol. Evol. 16:295–300, 2001. J. Salda˜ na, Continuous-time formulation of reactiondiffusion processes on heterogeneous metapopulations. Phys. Rev. E 78:012902, 2008. J. Salda˜ na, Modelling the spread of infectious diseases in complex metapopulations. Math. Model. Nat. 5:22–37, 2010. L. P. Sanders, B. S¨oderberg, D. Brockmann, and T. Ambj¨ornsson, Perturbative solution to susceptibleinfected-susceptible epidemics on networks. Phys. Rev. E 88:032713, 2013. M. Schl¨apfer, L. M. A. Bettencourt, S. Grauwin, M. Raschke, R. Claxton, Z. Smoreda, G. B. West, and C. Ratti, The scaling of human interactions with city size. J. R. Soc. Interface 11:20130789, 2014. M. J. Smith, S. Telfer, E. R. Kallio, S. Burthe, A. R. Cook, X. Lambin, and M. Begon, Host-pathogen time series data in wildlife support a transmission function between density and frequency dependence. Proc. Natl. Acad. Sci. U.S.A. 106:7905–7909, 2009.