Effects of diffusion rates on epidemic spreads in metapopulation networks

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.


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 [1,2,5,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 [2,[7][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.

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 N nodes, the N by N adjacency matrix is denoted by A; A ij = 1 when node i and node j are adjacent, and A ij = 0 otherwise. A is symmetric and we set A ii = 0 (1 ≤ i ≤ N).
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 Nρ 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 i occurs at a rate of βρ S,i ρ I,i , where ρ S,i and ρ I,i are the numbers of susceptible and infected individuals at node i, respectively. This assumption implies all-to-all interaction within each metapopulation.
At the same time, individuals perform a simple random walk. In infinitesimal time ∆t, a susceptible (infected) individual at node i with degree k i moves to one of its neighboring nodes with equal probability ≈ D S ∆t/k i (≈ D I ∆t/k i ), where D S (D I ) is the diffusion rate for the susceptible (infected) individual and k i is the degree of node i.

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 [2,[7][8][9][10] and solutions are known for uncorrelated random networks [13,[18][19][20], we are concerned with the effect of the diffusion rate on the epidemic threshold in general heterogeneous networks.
The master equations are given by To derive the epidemic threshold, we consider the steady state in which ρ I,i (1 ≤ i ≤ N) is infinitesimally small. In this situation, Eq. (1) represents the master equation for the continuoustime pure diffusion of susceptible individuals because βρ S,i ρ I,i ≈ 0 and µρ I,i ≈ 0. Given D S > 0, the steady state for the number of susceptible individuals is given by Eq. (4) has ρ * I,i = 0 (1 ≤ i ≤ N) as the disease-free solution. The destabilization of the diseasefree 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 β c . 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 D I = ∞, Eq. (4) implies where ρ I ≡ N i=1 ρ I,i /N ≈ 0 is the total density of the infected individual. By substituting Eq. (5) in Eq. (4) and taking the summation over i, we obtain which reproduces the threshold β c = µ k 2 /ρ k 2 known for uncorrelated networks [13]. When D I = ∞, heterogeneous degree distributions decrease the epidemic threshold in arbitrary metapopulation networks.
When D I = 0, Eq. (4) represents N decoupled dynamics, and β c is equal to the epidemic threshold at an isolated node i that hosts the largest number of individuals among the N nodes.
We obtain where k max is the degree of node i [20]. When β > β c , the nodes that satisfy k i > k µ/(βρ) can eventually have infected individuals, whereas the other nodes can not. In heterogeneous networks, we obtain Therefore, the value of β c for D I = 0 is smaller than that of β c for D I = ∞. Although this fact apparently indicates that the endemic state occurs more easily for D I = 0 than for D I = ∞, this statement is somewhat misleading because D I = 0 implies that the infection does not travel across metapopulations. When D I = 0, the infection is confined in the metapopulations that contain index patients.
Assume that D S and D I are arbitrary. We develop a method to calculate the epidemic threshold for an arbitrary structure of networks and diffusion rates. The Jacobian J = (J ij ) of the dynamics represented by Eq. (4) around the disease-free solution is given by where δ ij is the Kronecker delta. When the real parts of all the N eigenvalues of J 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].
J is isospectral to where diag(d 1 , . . . , d N ) denotes the diagonal matrix with the ith diagonal element equal to d i . Because is symmetric, all the eigenvalues of J ′ and hence those of J are real. Let λ max denote the maximum eigenvalue of J.
Fix β and express λ max using the Rayleigh quotient: where ⊤ denotes the transpose. Suppose that the maximum of Eq. (11) is realized by x =x, i / k , where λ max is the largest eigenvalue before β is replaced by β + ∆β. Therefore, λ max monotonically increases with β; this guarantees that the minimum value of β satisfying λ max = 0 is equal to β c . The condition λ max = 0 is equivalent to With regard to the value of β c , Eq. (12) is equivalent to Equation (13)  where is equal to β c . For arbitrary networks, we can reliably compute the minimum eigenvalue of M in O(N 3 ) time by combining the Cholesky decomposition, which is applicable to symmetric matrices, and the inverse power method.
For uncorrelated networks, substituting A ij = k i k j /( k N) in Eq. (14) yields The inverse of Eq. (15) is given by Because all the eigenvalues of M are positive, all the eigenvalues of M −1 are also positive.
Therefore, we can rapidly calculate β c by applying the standard power method to M −1 .

Epidemic threshold increases with the diffusion rate
In this section, we show that β c monotonically increases with D I in arbitrary heterogeneous networks. To this end, we decompose M as where and L is a Laplacian matrix such that its minimum eigenvalue is 0. All the other eigenvalues of L 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].

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 µ = 1 and population density ρ = 50 without loss of generality and vary the infection rate β. We first examine the SIS dynamics on scale-free networks generated by the Barabási-Albert (BA) model [32] with N = 200 and mean degree 6 (i.e., 3 links are added per new node). To generate a network, N − 3 = 197 nodes are sequentially added to a triangle according to the preferential attachment. The degree distribution of the network is the power law p(k) ∝ k −3 [32]. For a generated scale-free network, in Fig. 1(a), we compare β c obtained theoretically and numerically for a range of D. The theoretical value of β c is equal to the minimum eigenvalue of matrix M given by Eq. (14). We obtain a numerical estimate of β c 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 D, we generate 400 realizations of the scale-free network and run the SIS dynamics 5 times. The mean fraction of infected individuals ρ * I is calculated on the basis of 2000 runs. The dependence of ρ * I , normalized by ρ, on β and D are shown in Fig. 1(b). In agreement with Fig. 1(a), the epidemic threshold increases with the diffusion rate.
However, when D = 0, ρ * I 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 D = 0.1, the infection reaches multiple metapopulations above the epidemic threshold, but only to a limited extent because of a small diffusion rate. When D = 1, above the epidemic threshold, the infection seems to spread to more metapopulations and individuals than when D = 0.1. A visible fraction of infected individuals (ρ * I /ρ ≈ 0.03) survives between β ≈ 0.4 and β ≈ 0.6, for which larger diffusion rates do not allow the endemic state. When D = 20, ρ * I 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 D = 20 is nonetheless larger than that for D = 0.1 and D = 1, which is consistent with Fig. 1(a).
To quantify the localization of infected individuals, in particular when an appropriately small D (i.e., D = 1) yields relatively many infected individuals (i.e., 0.4 ≤ β ≤ 0.6), we plot the inverse participation ratio defined by Y 2 = N i=1 (ρ I,i /ρ I ) 2 (see [33] for an exemplary use of Y 2 for epidemic dynamics on networks). When Y 2 is of the order of unity, the infected individuals are localized in a small number of metapopulations. When Y 2 = O(1/N), the infected individuals have spread to O(N) metapopulations. In Fig. 1(c), Y 2 are plotted for the values of D used in Fig. 1(b). For D = 0, Y 2 = 1 irrespective of β, as expected. For D = 0.1, Y 2 is relatively large for β ≈ 0.5. For D = 1, Y 2 in this range of β becomes smaller to approach the values for D = 20 (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 [15,16,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 ρ * R 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 D facilitates spreads of infection to many metapopulations, we measure Y 2 . For the SIR dynamics, we need to modify the definition of Y 2 because runs without any secondary infections lead to Y 2 = 1. Runs with a minor fraction of infected individuals also yield a spuriously large Y 2 .
Therefore, in the calculation of Y 2 , we exclude the runs in which less than 1% of the individuals are eventually infected. The dependence of Y 2 on β and D 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 N = 332 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 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 N = 200 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 β c is independent of D for the regular random graph. As shown in Fig. 4(a), the numerically obtained β c is largely independent of D and close to the theoretical value, albeit there is some disagreement for small D. For a range of β, the normalized ρ * I and ρ * R for the SIS and SIR models are shown in Figs. 4(b) and 4(c), respectively. Because β c little depends on D, a large diffusion rate results in a large fraction of infected individuals for any β.

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.