key: cord-0843801-ayzbh7bw
authors: Gong, Yong-Wang; Song, Yu-Rong; Jiang, Guo-Ping
title: Epidemic spreading in metapopulation networks with heterogeneous infection rates
date: 2014-12-15
journal: Physica A
DOI: 10.1016/j.physa.2014.08.056
sha: c6d8ebcd8dfac7fdfa27ca40d2d4f44f24bcd514
doc_id: 843801
cord_uid: ayzbh7bw

In this paper, we study epidemic spreading in metapopulation networks wherein each node represents a subpopulation symbolizing a city or an urban area and links connecting nodes correspond to the human traveling routes among cities. Differently from previous studies, we introduce a heterogeneous infection rate to characterize the effect of nodes’ local properties, such as population density, individual health habits, and social conditions, on epidemic infectivity. By means of a mean-field approach and Monte Carlo simulations, we explore how the heterogeneity of the infection rate affects the epidemic dynamics, and find that large fluctuations of the infection rate have a profound impact on the epidemic threshold as well as the temporal behavior of the prevalence above the epidemic threshold. This work can refine our understanding of epidemic spreading in metapopulation networks with the effect of nodes’ local properties.

Epidemic spreading is one of the major threats to human health. The recent outbreaks of human epidemics such as SARS (severe acute respiratory syndrome) and H1N1 influenza have caused the death of many people [1, 2] . Analogously, the propagation of malware on computer networks and cell-phone networks can cause system damage, information leakage, and even economic loss [3] [4] [5] [6] [7] [8] . So, controlling epidemic spreading effectively has become an important issue today. To this end, it is necessary for us to gain a deep understanding of epidemic dynamics in the real-world environments.

In the last decade, epidemic spreading in complex networks has been extensively explored by researchers from a number of different disciplines such as theoretical biology, statistical physics, mathematics, and computer science [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] . In most studies, the substrate network is represented as a static network where each node stands for a single individual and edges between individuals are fixed. The remarkable result is that the epidemic threshold is a positive value in homogeneous networks such as regular lattices and random networks, but tends to zero in a wide range of scale-free (SF) networks [9] . On the other hand, some other studies investigated epidemic spreading in the dynamical network model with the assumption that all mobile agents have the same action radius and move randomly on a plane satisfying periodic boundary conditions [19] [20] [21] [22] [23] [24] . The mobility of individuals has been proved to have a significant impact on the epidemic behavior. For example, Buscarino et al. studied disease spreading in a system of random walkers which additionally perform long-distance jumps, and found that a small percentage of jumps in the agent motion are sufficient to produce a large drop in the epidemic threshold [21] .

Recently, the metapopulation epidemic models based on the reaction-diffusion (RD) process have been introduced to investigate epidemic spreading in a spatially structured social network [25] [26] [27] [28] [29] [30] . In these models, the social network is characterized by the metapopulation network where each node represents a city (subpopulation) consisting of any number of individuals and links denote transportation routes among cities. The individuals' reaction (infection) process occurs within nodes and the diffusion (mobility) process happens among nodes. In the diffusion process, some infected individuals as carriers of a disease may infect previously unaffected nodes. The network topology and diffusion pattern are two key factors in understanding such an epidemic spreading [25, 26, 29] . More recently, people have begun to pay much attention to how the nodes' local properties, such as heterogeneous length of stay of hosts' movements [31] , city-size heterogeneity [32] , local population structure [33] and local epidemic curing rate [34] , affect the metapopulation epidemic dynamics.

In previous metapopulation epidemic dynamics studies, the infection rate has been defined as a same constant for each node. In reality, at different subpopulations, human disease infectivity may possess discrepancies because of social, geographical or cultural aspects. It has been shown that the location-specific factor induces the substantial variation of disease incidence between populations [35, 36] , and the impact of the epidemic in European countries is highly different because of the marked differences in the sociodemographic structure of European populations [37] . In the present work, we introduce a heterogeneous infection rate to characterize the effect of the distinctive nodes' local properties, such as population density, individual health habit, and social conditions, on epidemic infectivity, and study its impact on the metapopulation epidemic dynamics. By using a mean-field (MF) approach and Monte Carlo simulations, it is shown that the heterogeneity of the infection rate has a significant impact on the epidemic threshold as well as the temporal behavior of the prevalence above the epidemic threshold. Note that some related works studied epidemic dynamics on Erdös-Rényi (ER) networks with other heterogeneous infection rates across links that are proportional to the intensity of interactions between people [38, 39] .

The paper is organized as follows. In Section 2, the metapopulation epidemic model with a heterogeneous infection rate is presented based on the RD process. In Section 3, we analyze theoretically the impact of the heterogeneous infection rate on the epidemic threshold by using a MF approach in detail. Section 4 performs extensive numerical simulations to validate the theoretical findings in Section 3, and additionally studies on the temporary behavior of the prevalence above the epidemic threshold.

Since most real-world networks have a heterogeneous topological structure, here we consider a heterogeneous metapopulation network with V nodes and a power-law degree distribution (i.e., p(k) ∼ k −r ), in which there exists a total population of N individuals. The RD process on this network is described as follows. First, reaction (infection) occurring within each node i follows the dynamics of the classical susceptible-infected-susceptible (SIS) model with two parameters β i (a local infection rate related to the properties of node i) and µ (a global constant recovery rate). The reaction process can be schematically represented by

The first sub-equation means that each susceptible (S) individual becomes infected with probability β i if it meets an infected (I) individual, and the second sub-equation indicates that the infected individuals are cured and become susceptible again with probability µ. Within node i, let S i (t) and I i (t) denote the number of infected and susceptible individuals at time step t, respectively. Thus the number of new infected individuals within node i at time step t can be represented as β i S i (t)I i (t) when each susceptible individual may react with all the infected individuals in the same node [25] .

For the diffusion process which happens at the end of reaction at each time step, as in Refs. [25, 32, 33] , we take a simplistic view and use a stochastic migration model (i.e., each individual selects randomly one of the neighbor nodes as its moving destination node) with coefficients 0 ≤ D S ≤ 1 and 0 ≤ D I ≤ 1 that denote the diffusion rates of the susceptible and infected individuals, respectively. This implies that at each time step susceptible and infected individuals in each node with degree k will move into one of their neighbor nodes with probabilities D S /k and D I /k, respectively. The schematic representation of the reaction and diffusion process is shown in Fig. 1 , where k i (or k j ) is the degree of node i (or node j).

In our model, how to reasonably define the local infection rate β i in each node i is an important matter. We assume statistical equivalence for subpopulations of the same degree, which has been successfully applied to many dynamical processes on complex networks [25, 26, 30, 31] ; thus it is reasonable for us to provide a same local infection rate to the same degree nodes (i.e., β i = β j if k i = k j ). Here, the local infection rate for nodes with degree k is defined as where ⟨k α ⟩ =  k k α P(k). In Eq. (2), the parameter β is a constant, and the exponent α is a tunable variable by which we can characterize different heterogeneity levels of the infection rate. Obviously, when α ̸ = 0, Eq. (2) represents a heterogeneous infection rate at the subpopulation level.

Considering that population density, health habits, and social conditions may have a different impact on the infection rate [37, 40] , we further introduce two infection regimes based on Eq. (2): negative-correlation infection regime (NIR) (i.e., the case of α < 0) and positive-correlation infection regime (PIR) (i.e., the case of α > 0). For NIR, the infection rate is negatively related to the node degree. That is to say, the infection rate in nodes with higher degree is smaller than that in nodes with lower degree. This may be understood as follows. Individuals in well-connected central cities (i.e., nodes with higher degree) generally have better health habits, such as washing hands and wearing masks, and have more opportunities for vaccinations than those in peripheral locations (i.e., nodes with lower degree) when an epidemic breaks out. Conversely, for PIR, the infection rate is positively related to the node degree. That is, central cities have a larger infection rate compared with peripheral locations. This can be explained by the fact that central cities generally have a higher population density leading to closer distance and stronger interaction among individuals [41] . This in some sense means that a susceptible individual can be infected with higher probability when it contacts with an infected individual. In the special case of α = 0, we have a homogeneous infection rate β k = β that has been considered in the previous studies [25] [26] [27] [28] [29] [30] [31] [32] [33] .

Here, to be clear, in the real world, there are unlikely too large fluctuations of the infection rate between different cities (nodes) even if they have a great difference in local properties. So, similar to Refs. [42, 43] , the heterogeneity exponent α is limited to the interval −1 ≤ α ≤ 1 when epidemic dynamics is discussed in our study.

The theoretical analysis of epidemic dynamics in heterogeneous networks is generally a difficult task due to the large degree fluctuations. To this end, we define the relative densities ρ S,k (t) and ρ I,k (t) as in Ref. [25] denoting the average number of susceptible and infected individuals in these nodes with degree k at time step t, respectively. So we have

where V k is the number of nodes with degree k and sums run over all nodes j having degree k j equal to k. Based on the description of the model in Section 2 and following Refs. [25, 30, 32] , we formulate the dynamical reaction rate equations for ρ I,k (t) and ρ S,k (t) at the MF level where P(k ′ /k) is the conditional probability that a node with degree k connects to another node with degree k ′ , and Γ k (t) is the reaction kernel given by

The first term on the right-hand side of Eq. (4a) represents the variation of the number of the infected individuals due to local reaction within the nodes. The second and the third terms of Eq. (4a) denote the number of outgoing and incoming infected individuals due to migration including the new infected individuals generated by the reaction term Γ k (t), respectively. Similarly, we can understand the meaning of each term of Eq. (4b) for the susceptible individuals.

In the following, we will consider the case of uncorrelated networks in which the conditional probability does not depend on the originating node (i.e., P(k ′ /k) = k ′ P(k ′ )/ ⟨k⟩ [44] ). By substituting this equation into Eq. (4) and imposing ∂ρ I,k (t)/∂t = ∂ρ S,k (t)/∂t = 0, we have the stationary state equation of Eq. (4) that is independent of time, taking the form

where ρ S =  k P(k)ρ S,k and ρ I =  k P(k)ρ I,k denote the average densities of susceptible and infected individuals in the network at the stationary state, respectively, and

If let ρ represent the total average density of the population in the network, it will satisfy that

Multiplying Eq. (6a) by P(k) and summing over k, we have the general result

Inserting Eq. (9) into Eq. (6), a simple form of Eq. (6) can be obtained

This is a general form of the RD process for any diffusion rates D I and D S . One notes that a thorough analysis of Eq. (10) is not easy. To focus on how the heterogeneity of the infection rate influences the epidemic dynamics and for the sake of simplicity, we will only consider the following two special mobility patterns: (i) the susceptible individuals cannot move (D S = 0) and the infected individuals always move (D I = 1); (ii) both the susceptible and infected individuals always move (D S = D I = 1).

In this case, the susceptible individuals cannot move and the infected individuals instead always move to their neighbor nodes. From Eq. (10), we have

Substituting the definition of Γ k (it has been given by Eq. (5) but with ρ I,k (t) and ρ S,k (t) replaced by ρ I,k and ρ S,k , respectively) into the second sub-equation of Eq. (11), one gets

By replacing β k in Eq. (12) with Eq. (2), then multiplying the resulting equation by P(k) and summing over k, this will The second sub-equation of Eq. (13) is obtained by using Eq. (8) (i.e., ρ = ρ I + ρ S ). Clearly, the condition for the presence of epidemic phase is thus given by ρ I > 0. By setting ρ I = 0 we will obtain the epidemic threshold ρ c with respect to the total average density of the population in the network, represented as

From Eq. (14), it can be seen that ρ c is closely related to the heterogeneity exponent α, and its values are equal for the same absolute value of α. This implies that in the case of D S = 0 and D I = 1, the epidemic threshold is profoundly impacted by the heterogeneity of the infection rate. For the special case of α = 0 (i.e., the homogeneous infection rate), the epidemic threshold is simplified as ρ c = µ/β which is identical with that in Ref. [25] .

In this case, both the susceptible and infected individuals always move to their neighbor nodes. Then from Eq. (10) we have In addition, from Eqs. (7) and (9), one obtains

Inserting Eqs. (2) and (15) into Eq. (16) , and then after a simple algebraic calculation, it is easy to get

So, the condition for the presence of epidemic phase can be given by ρ I > 0. This implies that the epidemic threshold is

Similarly, we can see from Eq. (18) that the epidemic threshold ρ c is strongly impacted by the heterogeneity exponent α. For the special case of α = 0, we recover ρ c = µ⟨k⟩ 2 β⟨k 2 ⟩ [25] .

In order to understand more intuitively the impact of the heterogeneity of the infection rate on the epidemic threshold, the analytical curves of the epidemic threshold ρ c are plotted in Fig. 2 . In Fig. 2(a) , the two curves completely fall onto the same line and they are rising with the absolute value of α increasing.

This means that, in the Case 1 mobility pattern, PIR and NIR infection regimes have the same impact on the epidemic threshold (i.e., increasing the epidemic threshold by the same extent), and the higher the heterogeneity of the infection rate is (i.e., the larger the absolute value of α is), the larger the epidemic threshold is. By contrast, for the Case 2 mobility pattern, PIR and NIR infection regimes have a different impact on the epidemic threshold as shown in Fig. 2(b) . That is, the PIR infection regime decreases the epidemic threshold and the NIR infection regime increases it. In addition, such an impact becomes stronger with the increase of the heterogeneity of the infection rate. These are the main results in our study.

In this section, we perform an extensive set of Monte Carlo simulations to validate the theoretical predictions in Section 3. Here we first carry out simulations in a configuration network generated by an uncorrelated configuration model (UCM) [45] .

The configuration network has V = 2000 nodes with the power-law degree distribution p(k) ∼ k −2.5 (2 ≤ k ≤ 42) and average degree ⟨k⟩ = 3.9.

The dynamics proceeds in parallel and considers discrete time steps representing the unitary time scale τ of the process.

The reaction and diffusion rates are therefore converted into probabilities at each time step and the system is updated according to the following rules [25] : in each node with degree k, each infected individual is cured and becomes susceptible again with probability µτ . At the same time, each susceptible individual becomes an infected individual with probability 1−(1−β k τ ) n I , where n I is the total number of the infected individuals in the node. After all nodes have been updated for the reaction, we simulate the diffusion process: each individual moves into a randomly chosen neighbor node with probabilities D S τ and D I τ for the susceptible and infected individuals, respectively. In our simulation, the dynamics parameters are µ = 0.05, β = 0.05 and τ = 1 if not otherwise specified. −1) is the same as that in the case of α = 0.5 (or α = 1). We also find that NIR and PIR infection regimes increase the values of the epidemic threshold ρ c by the same extent, and the larger the absolute value of α is, the larger the simulated epidemic threshold is. These simulation results match well the analytical prediction (see Eq. (14)).

For the Case 2 mobility pattern (see Fig. 3(b) ), the curve with α = 0 separates the curves with α > 0 (PIR) from those with α < 0 (NIR). This validates that the PIR infection regime decreases the epidemic threshold and NIR increases it. Here for α equal to 0, 0. The simulation results (see Fig. 3 (a)) indicate that for the Case 1 mobility pattern and the same absolute value of α, PIR and NIR infection regimes have the same impact on the epidemic threshold as well as the stationary infected density. In the following, we want to know in such a case how these two infection regimes affect the temporal behavior of the prevalence above the epidemic threshold. To this aim, in Fig. 4 , we plot the time evolution of the normalized epidemic prevalence ρ I (t)/ρ and the corresponding spread speed v(t) (the inset) for different α values in the configuration network with the average population density ρ = 3, where v(t) is here defined as Refs. [46, 47] v(t) = ∂(ρ I (t)/ρ) ∂t ≈ ρ I (t)/ρ − ρ I (t − 1)/ρ.

It can be seen from Fig. 4 that in PIR (α > 0) and NIR (α < 0) infection regimes the time evolution curves of the epidemic prevalence are different for the same |α| ( Fig. 4 (a) |α| = 0.5 and Fig. 4 (b) |α| = 1). More specifically, at the early time steps, the epidemic prevalence in the PIR infection regime grows faster than that in the NIR infection regime. However, as time goes on, it is just the opposite. As a result, after a long time, the two infection regimes make epidemic spreading reach a same equilibrium state. Similarly, in the NIR infection regime, epidemic has a faster speed at the early time steps and reaches the infection peak earlier than that in the NIR infection regime (see the inset of Fig. 4) .

This rationale behind this finding is as follows. According to Eq. (2), in nodes with higher degree, the infection rates are larger in the PIR infection regime than those in the NIR infection regime. In nodes with lower degree, it is just the opposite. On the other hand, as a result of the stochastic migration process, most individuals will gather in nodes with higher degree, and this is convenient for epidemic invasion in those nodes. In such cases, epidemic speed in the PIR infection regime is larger than that in the NIR infection regime at early time steps due to the fact that infection mainly occurs in nodes with higher degree. As time goes on, the epidemic diffuses to nodes with lower degree. From this time, epidemic speed in the NIR infection regime begins to increase because of larger infection rates in those nodes with lower degree.

In order to further confirm the validity of the theoretical results in a real-world setting, we perform simulations in the US air transportation network with 332 nodes and 2126 links [48, 49] . In the network, each node represents an airport and each link denotes a direct flight between two airports. The network link weight is ignored and the cumulative degree distribution is shown in Fig. 5 . The simulation results about the epidemic threshold and the temporal behavior of the prevalence above the epidemic threshold in the US air transportation network are shown in Figs. 6 and 7, respectively. It can be seen from these two figures that the simulation results are also in accord with theoretical results obtained in Section 3. The theoretical and simulated results above are based on the two special mobility patterns. To understand how the heterogeneous infection rate affects the epidemic threshold in a more general mobility pattern, we give the simulation results about the epidemic threshold with D s = 0.5 and D I = 0.8 in the configuration network ( Fig. 8(a) ) and the US airport transportation network (Fig. 8(b) ). It can be seen from Fig. 8 that for a more general mobility pattern, the NIR infection regime increases the epidemic threshold and the PIR infection regime decreases it, which are similar to the results in the special case of D s = D I = 1.

In this paper, we have proposed a heterogeneous infection rate of which the heterogeneity is governed by the exponent α and further denoted the case of α < 0 as the negative-correlation infection regime (NIR) and the case of α > 0 as the positive-correlation infection regime (PIR). By a MF approach and Monte Carlo simulations, the impact of these two infection regimes on epidemic dynamics has been studied in heterogeneous metapopulation networks with individuals diffusing randomly. Interesting results have been found as follows. In the case that the susceptible individuals cannot move and the infected individuals always move (i.e., Case 1 mobility pattern), NIR and PIR infection regimes produce the same impact on the epidemic threshold, i.e., increasing the epidemic threshold by the same extent. However, when both the susceptible and infected individuals always move (i.e., Case 2 mobility pattern), the PIR infection regime decreases the threshold and NIR increases it. Furthermore, we have also provided some simulations of the temporal behavior of the prevalence above the epidemic threshold in the Case 1 mobility pattern. It is shown that, compared with the NIR infection regime, the PIR infection regime favors the epidemic invasion at early time steps, but in the long-time limit, these two infection regimes have the same impact on the final epidemic prevalence.

It is worth stressing that the present theoretical analysis is based on two special mobility patterns and so the analytical solutions are not universal. For the general case, we have only given the simulation results. In spite of this, our study has uncovered an important aspect (i.e., the impact of local properties of the node on epidemic infectivity) that must be considered in modeling and understanding epidemic dynamics in social networks. Certainly, the full theoretical solution of the present model may be obtained by the integration method used in Ref. [50] , which will be our future work.

On the stability of the malware free equilibrium in cell phones networks with spatial dynamics