Effects of diffusion rates on epidemic spreads in metapopulation networks

# Effects of diffusion rates on epidemic spreads in metapopulation networks

Naoki Masuda

Department of Mathematical Informatics,
The University of Tokyo,
7-3-1 Hongo, Bunkyo, Tokyo 113-8656, Japan
PRESTO, Japan Science and Technology Agency,
4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan

masuda@mist.i.u-tokyo.ac.jp
###### Abstract

It is often useful to represent the infectious dynamics of mobile agents by metapopulation models. In such a model, metapopulations form a static network, and individuals migrate from one metapopulation to another. It is known that heterogeneous degree distributions of metapopulation networks decrease the epidemic threshold above which epidemic spreads can occur. We investigate the combined effect of heterogeneous degree distributions and diffusion on epidemics in metapopulation networks. We show that for arbitrary heterogeneous networks, diffusion suppresses epidemics in the sense of an increase in the epidemic threshold. On the other hand, some diffusion rates are needed to elicit epidemic spreads on a global scale. As a result of these opposing effects of diffusion, epidemic spreading near the epidemic threshold is the most pronounced at an intermediate diffusion rate. The result that diffusion can suppress epidemics contrasts with that for diffusive SIS dynamics and its variants when individuals are fixed at nodes on static networks.

## 1 Introduction

Infectious diseases are transmitted through social contacts between individuals. The relationships between the structure of social networks and the number of infected individuals have been extensively studied in mathematical epidemiology and network sciences. Such studies have established that heterogeneity in contact rates enhances epidemics [1, 2, 3, 4].

The contact networks underlying, for example, sexually transmitted diseases of humans and computer viruses on the Internet are considered to be static on the time scale of epidemics. However, the social networks underlying humans’ prevailing infectious diseases such as influenza are presumably dynamic even during a day. The dynamics of networks are induced by the migration of individuals among residences, workplaces, places for social activities, and so on. Other animals also migrate between habitats. Metapopulation models are useful in describing epidemics and ecological invasions in such a situation [5, 1, 2, 6]. A node in such a model represents a metapopulation or a habitat, and not an individual. A link represents a physical pathway connecting a pair of metapopulations. Individuals travel from one node to another. Simultaneously, interactions between individuals, such as infection, can occur in each metapopulation.

The basic reproduction number and the condition for the occurrence of global epidemics have been theoretically determined for the susceptible-infected-suscepetible (SIS) dynamics and its variants on metapopulation networks with general connectivity profiles [7, 2, 8, 9, 10]. Real metapopulation networks relevant to epidemics, such as networks of urban metapopulations and networks of airports, often have heterogeneous degree distributions; some metapopulations are highly connected as compared to many others [3, 4, 11]. Colizza et al. put important efforts for understanding infectious dynamics on complex metapopulation networks [12, 13, 14, 15, 16]. In particular, they showed that heterogeneous degree distributions enhance epidemics in uncorrelated networks of metapopulations (also see [17]). Similar results have been established for models of epidemics on static networks of individuals (see [3, 4] for reviews). However, the relationship between the epidemic threshold of the infection rate (or population density) above which large epidemics can occur and the degree distribution differs in these two classes of models.

For metapopulation models, effects of stochasticity in infection, recovery, and migration dynamics were also analyzed in a recent paper [18]. Other authors also investigated the threshold for infection and epidemic dynamics in uncorrelated heterogeneous networks of metapopulations [19, 20].

In the present work, we investigate the effect of diffusion on epidemics in metapopulation networks on which individuals diffuse. Diffusion increases the possibility of epidemics for the SIS model [21, 22] and its variants [23, 24] in conventional static networks of individuals. This is presumably because diffusion helps infected individuals to contact new susceptibles to be infected. The effect of traffic-driven diffusion on the epidemic threshold was also studied recently [25]. We show that for epidemics in metapopulation networks, diffusion prohibits rather than enhances epidemic spreads in the sense of the epidemic threshold. Although this result has been implicitly recognized in previous literature [20], we show it by developing an analytical method to calculate the epidemic threshold for arbitrary networks and diffusion rates. We then support this result by numerical simulations and further show that an intermediate diffusion rate magnifies epidemics near the epidemic threshold. Qualitatively identical behavior is observed in numerical simulations of the susceptible-infected-recovered (SIR) dynamics on metapopulation networks.

## 2 Model

We assume connected and undirected networks of metapopulations; extending the following results to the case of weighted networks is straightforward. Each node represents a metapopulation or habitat and accommodates any number of individuals. For a network having nodes, the by adjacency matrix is denoted by ; when node and node are adjacent, and otherwise. is symmetric and we set ().

We consider the SIS dynamics of diffusive individuals on this network using the continuous-time formulation developed in [19, 20]. The use of the original discrete-time formulation of the model [13] does not essentially change the following results. The network contains individuals that independently diffuse on the network. Each individual stays at a node and takes either a susceptible or an infectious state. The SIS dynamics with infection rate and recovery rate occur in each metapopulation. We assume mass interaction and therefore do not normalize the infection rate by the number of individuals in the metapopulation [13]. By doing so, we focus on situations in which metapopulations that host relatively many individuals would be a source of epidemics once an infected individual is introduced to such a metapopulation. The infection event at node occurs at a rate of , where and are the numbers of susceptible and infected individuals at node , respectively. This assumption implies all-to-all interaction within each metapopulation.

At the same time, individuals perform a simple random walk. In infinitesimal time , a susceptible (infected) individual at node with degree moves to one of its neighboring nodes with equal probability (), where () is the diffusion rate for the susceptible (infected) individual and is the degree of node .

## 3 Epidemic threshold for arbitrary metapopulation networks

We derive the epidemic threshold above which the endemic state (i.e., survival of infected individuals) can occur. Although general solutions have been formulated [7, 2, 8, 9, 10] and solutions are known for uncorrelated random networks [13, 19, 20, 18], we are concerned with the effect of the diffusion rate on the epidemic threshold in general heterogeneous networks.

The master equations are given by

 dρS,idt =−βρS,iρI,i+μρI,i−DSρS,i+DSN∑j=1AjikjρS,j, (1) dρI,idt =βρS,iρI,i−μρI,i−DIρI,i+DIN∑j=1AjikjρI,j. (2)

To derive the epidemic threshold, we consider the steady state in which () is infinitesimally small. In this situation, Eq. (1) represents the master equation for the continuous-time pure diffusion of susceptible individuals because and . Given , the steady state for the number of susceptible individuals is given by

 ρ∗S,i=ki⟨k⟩ρ, (3)

where is the average degree and , the total population density. Substituting Eq. (3) in Eq. (2) yields

 dρI,idt≈βρ⟨k⟩kiρI,i−μρI,i−DIρI,i+DIN∑j=1AjikjρI,j. (4)

Eq. (4) has () as the disease-free solution. The destabilization of the disease-free solution implies an endemic state, that is, () [7, 10]. We focus on the threshold value of above which the endemic state is induced, which is denoted by . We note that the threshold obtained below can also be expressed in terms of , as is done in previous literature [13, 19, 20], because and appear only as the multiple in Eq. (4).

In the strong diffusion limit , Eq. (4) implies

 ρI,i=ki⟨k⟩ρI, (5)

where is the total density of the infected individual. By substituting Eq. (5) in Eq. (4) and taking the summation over , we obtain

 dρIdt=(βρ⟨k2⟩⟨k⟩2−μ)ρI, (6)

which reproduces the threshold known for uncorrelated networks [13]. When , heterogeneous degree distributions decrease the epidemic threshold in arbitrary metapopulation networks.

When , Eq. (4) represents decoupled dynamics, and is equal to the epidemic threshold at an isolated node that hosts the largest number of individuals among the nodes. We obtain

 βc=⟨k⟩μρkmax, (7)

where is the degree of node [20]. When , the nodes that satisfy can eventually have infected individuals, whereas the other nodes can not. In heterogeneous networks, we obtain . Therefore, the value of for is smaller than that of for . Although this fact apparently indicates that the endemic state occurs more easily for than for , this statement is somewhat misleading because implies that the infection does not travel across metapopulations. When , the infection is confined in the metapopulations that contain index patients.

Assume that and are arbitrary. We develop a method to calculate the epidemic threshold for an arbitrary structure of networks and diffusion rates. The Jacobian of the dynamics represented by Eq. (4) around the disease-free solution is given by

 Jij=(βρki⟨k⟩−μ−DI)δij+DIAjikj,(1≤i,j≤N) (8)

where is the Kronecker delta. When the real parts of all the eigenvalues of are negative, the disease-free solution is stable. Otherwise, the endemic state emerges. Similar treatments were carried out for this model on uncorrelated random networks [20] and the SIS dynamics on static networks of individuals [26, 27, 28, 29].

is isospectral to

 J′=(J′ij)≡diag(1√k1,…,1√kN)Jdiag(√k1,…,√kN), (9)

where denotes the diagonal matrix with the th diagonal element equal to . Because

 J′ij=(βρki⟨k⟩−μ−DI)δij+DIAji√kikj (10)

is symmetric, all the eigenvalues of and hence those of are real. Let denote the maximum eigenvalue of .

Fix and express using the Rayleigh quotient:

 λmax=max|x|=1x⊤J′x,(x∈RN) (11)

where denotes the transpose. Suppose that the maximum of Eq. (11) is realized by , where . When is replaced by , we obtain , where is the largest eigenvalue before is replaced by . Therefore, monotonically increases with ; this guarantees that the minimum value of satisfying is equal to . The condition is equivalent to

 detJ=0. (12)

With regard to the value of , Eq. (12) is equivalent to

 0 =⟨k⟩ρdet[diag(1k1,…,1kN)J] =det(βδij−μ+DIki⟨k⟩ρδij+DI⟨k⟩ρAijkikj)=0. (13)

Equation (13) indicates that the minimum eigenvalue of the symmetric matrix , where

 Mij=(μ+DIkiδij−DIAijkikj)⟨k⟩ρ, (14)

is equal to . For arbitrary networks, we can reliably compute the minimum eigenvalue of in time by combining the Cholesky decomposition, which is applicable to symmetric matrices, and the inverse power method.

For uncorrelated networks, substituting in Eq. (14) yields

 Mij=(μ+DIkiδij−DI⟨k⟩N)⟨k⟩ρ. (15)

The inverse of Eq. (15) is given by

 M−1ij=ρ⟨k⟩(μ+DI)(kiδij+DIkikjμ⟨k⟩N). (16)

Because all the eigenvalues of are positive, all the eigenvalues of are also positive. Therefore, we can rapidly calculate by applying the standard power method to .

## 4 Epidemic threshold increases with the diffusion rate

In this section, we show that monotonically increases with in arbitrary heterogeneous networks. To this end, we decompose as

 M=(B+DIL)⟨k⟩ρ, (17)

where

 B=μ×diag(1k1,…,1kN) (18)

and

 L=diag(1k1,…,1kN)−Aijkikj. (19)

is a Laplacian matrix such that its minimum eigenvalue is 0. All the other eigenvalues of are positive for connected networks. The interlacing eigenvalue theorem (also called Weyl’s inequality) implies that the minimum eigenvalue of the summation of two symmetric matrices is not smaller than the sum of the minimum eigenvalues of the two individual matrices [30, 31]. Because the minimum eigenvalue of is equal to for and that of is equal to 0, the minimum eigenvalue of , i.e., for a given , is not smaller than for . More generally, if , we can apply the abovementioned arguments to the sum of and , i.e., . Therefore, monotonically increases with .

For regular graphs, that is, networks in which all the nodes have the same degree, both for and . Because monotonically increases with , holds true for any . The diffusion rate affects the epidemic threshold for heterogeneous metapopulation networks but not for homogeneous networks.

## 5 Numerical results

The results presented in the previous section imply that diffusion suppresses epidemics on heterogeneous metapopulation networks in the sense of an increased epidemic threshold. However, if the diffusion rate is too small, infected individuals may be localized within the metapopulations having index patients. We investigate the combined effect of the heterogeneity and the diffusion rate by carrying out direct numerical simulations of infectious dynamics.

We set the recovery rate and population density without loss of generality and vary the infection rate . For simplicity, we set . We start with one infected individual selected from the population with equal probability . The other individuals are initially susceptible. We measure the fraction of infected individuals in the steady state, denoted by . We calculate as the number of infected individuals after transient divided by and averaged over multiple realizations. The relaxation time is governed by for small and otherwise. Therefore, we set the transient length to 500 for and , and for .

We first examine the SIS dynamics on scale-free networks generated by the Barabási-Albert (BA) model [32] with and mean degree 6 (i.e., 3 links are added per new node). To generate a network, nodes are sequentially added to a triangle according to the preferential attachment. The degree distribution of the network is the power law [32]. For a generated scale-free network, in Fig. 1(a), we compare obtained theoretically and numerically for a range of . The theoretical value of is equal to the minimum eigenvalue of matrix given by Eq. (14). We obtain a numerical estimate of as the value of above which an infected individual survives in at least one of the 5000 realizations. The numerical results are in a good agreement with the theoretical results.

Next, we examine the SIS dynamics for suprathreshold . For fixed and , we generate 400 realizations of the scale-free network and run the SIS dynamics 5 times. The mean fraction of infected individuals is calculated on the basis of 2000 runs. The dependence of , normalized by , on and are shown in Fig. 1(b). In agreement with Fig. 1(a), the epidemic threshold increases with the diffusion rate.

However, when , is small for any . This is because infected individuals do not migrate to other metapopulations, and they are localized within the single metapopulation to which the index patient belongs even above the epidemic threshold. When , the infection reaches multiple metapopulations above the epidemic threshold, but only to a limited extent because of a small diffusion rate. When , above the epidemic threshold, the infection seems to spread to more metapopulations and individuals than when . A visible fraction of infected individuals () survives between and , for which larger diffusion rates do not allow the endemic state. When , is large for large values of , presumably because diffusion helps the infected individuals to reach metapopulations with small degrees. We stress that the epidemic threshold for is nonetheless larger than that for and , which is consistent with Fig. 1(a).

To quantify the localization of infected individuals, in particular when an appropriately small (i.e., ) yields relatively many infected individuals (i.e., ), we plot the inverse participation ratio defined by (see [33] for an exemplary use of for epidemic dynamics on networks). When is of the order of unity, the infected individuals are localized in a small number of metapopulations. When , the infected individuals have spread to metapopulations. In Fig. 1(c), are plotted for the values of used in Fig. 1(b). For , irrespective of , as expected. For , is relatively large for . For , in this range of becomes smaller to approach the values for (note the logarithmic scale in Fig. 1(c)).

In infectious diseases with which metapopulation models are concerned, such as influenza, the SIR model seems to have a wider applicability than the SIS model [12, 14, 34]. In the SIR model, the individuals that recover from the infected state do not return to the susceptible state but transit to the recovered state. When the basic reproduction number is assumed to be the same for each metapopulation, analytical results for the SIR model on heterogeneous metapopulation networks have been established [16, 15, 18]. However, the analysis seems difficult when the basic reproduction number depends on the metapopulation. In our model, highly populated metapopulations, which presumably make contact with many other metapopulations, have large basic reproduction numbers.

We carry out numerical simulations starting from an arbitrarily chosen infected individual until there is no infected individual. Then, we measure the final fraction of the recovered individuals averaged over 2000 realizations. The results shown in Fig. 2(a) are qualitatively the same as those for the SIS model shown in Fig. 1(b). To confirm that an appositely small facilitates spreads of infection to many metapopulations, we measure . For the SIR dynamics, we need to modify the definition of because runs without any secondary infections lead to . Runs with a minor fraction of infected individuals also yield a spuriously large . Therefore, in the calculation of , we exclude the runs in which less than 1% of the individuals are eventually infected. The dependence of on and shown in Fig. 2(b) are qualitatively the same as those for the SIS model shown in Fig. 1(c).

As an example of real heterogeneous networks, we simulate SIS and SIR dynamics on the network of US airports [11]. This network has nodes, 2126 links, and a long-tailed degree distribution. We ignore the link weight. A airport is considered to be a metapopulation that hosts traveling individuals. Figure 3(a) indicates that the values of obtained from direct numerical simulations are in a good agreement with the theoretical values. The suprathreshold behavior of the SIS and SIR dynamics shown in Figs. 3(b) and 3(c), respectively, is qualitatively the same as that on the BA model (Figs. 1(b) and 2(a)).

The combined effect of the heterogeneity and the diffusion rate is absent for homogeneous networks. To demonstrate this, we carry out numerical simulations on the regular random graph with and degree 6; all the nodes have degree 6 and are wired randomly under the condition that the generated network is connected. The theory obtained in Sec. 4 indicates that the epidemic threshold is independent of for the regular random graph. As shown in Fig. 4(a), the numerically obtained is largely independent of and close to the theoretical value, albeit there is some disagreement for small . For a range of , the normalized and for the SIS and SIR models are shown in Figs. 4(b) and 4(c), respectively. Because little depends on , a large diffusion rate results in a large fraction of infected individuals for any .

## 6 Discussion and Conclusions

We have shown that diffusion increases the epidemic threshold of the SIS dynamics in arbitrary heterogeneous networks. Nevertheless, a certain amount of diffusion is necessary to enable the transmission of epidemics across metapopulations. Numerical simulations of the SIR dynamics also yielded similar results. Similar conclusions are expected for other population dynamics used in ecology [6]. Indeed, dispersal decreases the population growth rate in ecological invasion models in which each metapopulation carries heterogeneous birth and death rates [35].

Our numerical results indicate that excessive diffusion suppresses epidemics in metapopulation networks near the epidemic threshold. This effect of diffusion contrasts with that on the diffusive SIS model and variants in static networks of individuals [21, 22, 23, 24]. We note that, in the metapopulation model, an increased diffusion rate can also suppress the spreading speed of epidemics. In heterogeneous static networks of individuals, the spreading speed in an early stage of transmission is controlled by the largest eigenvalue of the matrix that represents transmission of infection [33]. The counterpart in the metapopulation model is the matrix given by Eq. (8). The largest eigenvalue of the matrix given by Eq. (8) decreases with the diffusion rate; this can be shown by adapting the proof of the monotonic dependence of the epidemic threshold on the diffusion rate (Sec. 4).

The nonmonotonic dependence of the fraction of infected individuals on the diffusion rate near the epidemic threshold is derived from the fact that epidemics confined in a single metapopulation are irrelevant on the global scale. A similar distinction between global and local epidemics is also important when analyzing epidemics in static networks of individuals with community structure; epidemics confined in a single community are not usually regarded as global epidemics [36].

## Acknowledgments

We thank Takehisa Hasegawa for providing valuable comments on the manuscript. N.M. acknowledges the support through Grants-in-Aid for Scientific Research (Nos. 20760258 and 20540382, and Innovative Areas “Systems Molecular Ethology”) from MEXT, Japan.

## References

• [1] Anderson R M and May R M 1991 Infectious Diseases of Humans (Oxford: Oxford University Press)
• [2] Diekmann O and Heesterbeek J A P 2000 Mathematical Epidemiology of Infectious Diseases (Chichester: John Wiley & Sons, Ltd.)
• [3] Newman M E J. 2003 SIAM Rev. 45 167
• [4] Barrat A, Barthélemy M and Vespignani A 2008 Dynamical Processes on Complex Networks (Cambridge: Cambridge University Press)
• [5] Rvachev L A and Longini Jr. I M 1985 Math. Biosci. 75 3
• [6] Hanski I and Ovaskainen O 2000 Nature 404 755
• [7] Diekmann O, Heesterbeek J A P and Metz J A J 1990 J. Math. Biol. 28 365
• [8] Arino J and van den Driessche P 2003 LNCIS 294 135
• [9] Arino J, Davis J R, Hartley D, Jordan R, Miller J M and van den Driessche P 2005 Math. Med. Biol. 22 129
• [10] Allen L J S, Bolker B M, Lou Y and Nevai A L 2007 SIAM J. Appl. Math. 67 1283
• [11] Batagelj V and Mrvar A 2006 Pajek datasets.
• [12] Colizza V, Barrat A, Barthélemy M and Vespignani A 2006 Proc. Natl. Acad. Sci. USA 103 2015
• [13] Colizza V, Pastor-Satorras R and Vespignani A 2007 Nat. Phys. 3 276
• [14] Colizza V, Barrat A, Barthelemy M, Valleron A-J and Vespignani A 2007 PLoS Med. 4 e13
• [15] Colizza V and Vespignani A 2008 J. Theor. Biol. 251 450
• [16] Colizza V and Vespignani A 2007 Phys. Rev. Lett. 99 148701
• [17] Fulford G R and Roberts M G 2002 Theor. Pop. Biol. 61 15
• [18] Barthélemy M, Godrèche C and Luck J-M 2010 arXiv:1004.1553v1
• [19] Saldaña J 2008 Phys. Rev. E 78 012902
• [20] Juher D, Ripoll J and Saldaña J 2009 Phys. Rev. E 80 041920
• [21] Jensen I and Dickman R 1993 J. Phys. A 26 L151
• [22] Konno N 1995 J. Theor. Prob. 8 833
• [23] Dickman R and Tomé T 1991 Phys. Rev. A 44 4833
• [24] Fiore C E and de Oliveira M J 2004 Phys. Rev. E 70 046131
• [25] Meloni S, Arenas A and Moreno Y 2009 Proc. Natl. Acad. Sci. USA 106 16897
• [26] Boguñá M and Pastor-Satorras R 2002 Phys. Rev. E 66 047104
• [27] Boguñá M, Pastor-Satorras R and Vespignani A 2003 Phys. Rev. Lett. 90 028701
• [28] Wang Y, Chakrabarti D, Wang C and Faloutsos C 2003 Proc. 22nd Int. Symp. on Reliable Distributed Systems (SRDS’03) 25
• [29] Ganesh A, Massoulié L and Towsley D 2005 Proc. of 24th IEEE Conf. on Computer Communications (INFOCOM’05) 2 1455
• [30] Horn R A and Johnson C R 1985 Matrix Analysis (Cambridge: Cambridge University Press)
• [31] Cvetković D, Rowlinson P and Simić S 2010 An Introduction to the Theory of Graph Spectra (Cambridge: Cambridge University Press)
• [32] Barabási A-L and Albert R 1999 Science 286 509
• [33] Barthélemy M, Barrat A, Pastor-Satorras R and Vespignani A 2004 Phys. Rev. Lett. 92 178701
• [34] Hufnagel L, Brockmann D and Geisel T 2004 Proc. Natl. Acad. Sci. USA 101 15124
• [35] Schreiber S J and Lloyd-Smith J O 2009 Am. Nat. 174 490
• [36] Masuda N 2009 New. J. Phys. 11 123018
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters