key: cord-0944715-dcmz4wcn
authors: Gong, Yong-Wang; Song, Yu-Rong; Jiang, Guo-Ping
title: Time-varying human mobility patterns with metapopulation epidemic dynamics
date: 2013-10-01
journal: Physica A
DOI: 10.1016/j.physa.2013.05.028
sha: eb526b1d311cd67f228f3ea1792e32e3dc6dc7c8
doc_id: 944715
cord_uid: dcmz4wcn

In this paper, explicitly considering the influences of an epidemic outbreak on human travel, a time-varying human mobility pattern is introduced to model the time variation of global human travel. The impacts of the pattern on epidemic dynamics in heterogeneous metapopulation networks, wherein each node represents a subpopulation with any number of individuals, are investigated by using a mean-field approach. The results show that the pattern does not alter the epidemic threshold, but can slightly lower the final average density of infected individuals as a whole. More importantly, we also find that the pattern produces different impacts on nodes with different degree, and that there exists a critical degree [Formula: see text]. For nodes with degree smaller than [Formula: see text] , the pattern produces a positive impact on epidemic mitigation; conversely, for nodes with degree larger than [Formula: see text] , the pattern produces a negative impact on epidemic mitigation.

To date, epidemic spreading in complex networks has been studied intensively by many researchers [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] . A key issue in this field is how to find the threshold of an epidemic outbreak, and the remarkable results have been obtained that the epidemic threshold is a positive value in homogeneous networks [1, 2] , but tends to zero in a wide range of SF (scale-free) networks with large size [3, 4] . Most of these studies mainly focused on epidemic dynamics in static networks, where each node represents a single individual and the connections between individuals are fixed.

However, many real social networks are heterogeneous metapopulation networks, wherein the nodes represent social units such as cities and urban areas, and the individuals located in the nodes represent people who can move among cities [14] [15] [16] . In this situation, the spatial structure of population and human mobility pattern are two primary elements that need to be considered in understanding large epidemic dynamics. To do this, recently, Colizza et al. studied the basic RD (reaction-diffusion) process in heterogeneous networks where the reaction (infection) takes place in the same node.

They assumed that the diffusion rates are unitary or null, and investigated the effect of network topology on the epidemic threshold [17] . Later, following this line, more realistic mobility patterns (diffusion processes), such as the traffic-dependent mobility pattern, the population-dependent mobility pattern, bidirectional mobility pattern, and recurrent mobility pattern, have been introduced to study epidemic dynamics in social networks by using the RD process [18] [19] [20] [21] [22] . In addition, more recently, Lund et al. have investigated the effects of city-size heterogeneity on epidemic dynamics in a metapopulation where small villages and big cities exist simultaneously with the basic RD process [23] .

Note that the above studies are all based on the assumption that human mobility rates are unchanged during the epidemic spreading. In fact, people generally have awareness and change their behaviors (including canceling their trips) to respond to an epidemic outbreak, and such behavior changes have been proved to affect the network structure and epidemic dynamics [24] [25] [26] [27] [28] . To the best of our knowledge, so far there has been no report, except Ref. [29] , that has explicitly considered the response of the variation of human travel rates in epidemic dynamics. In Ref. [29] , the authors assumed that the human mobility rate is adjusted according to the local epidemic prevalence in the destination subpopulation. In this case, only the probability of human traveling to the infected subpopulation is affected. However, in some other cases, a lethal epidemic outbreak can cause reduction of travel in global scope, not just in local infected regions, because of people fearing to be infected on trips. For example, in 2003, the spread of SARS (severe acute respiratory syndrome) in several provinces of China led to a decrease of travel all over the country [30] .

In this paper, different from Ref. [29] , we introduce a time-varying human mobility pattern to model the time variation of global human travel due to the emergence of an epidemic. In this pattern, human mobility rates are governed by the global epidemic prevalence and the fraction of people who have awareness of the epidemic spreading. In addition, not just by simulation method, we also theoretically explore the impacts of the mobility pattern on metapopulation epidemic dynamics. Interesting results are found as follows: (i) the time-varying mobility pattern does not alter the epidemic threshold, but it can slightly lower the final average density of infected individuals as a whole; (ii) more importantly, we also find that the pattern can produce different impacts on nodes with different degree, and that there exists a critical degree. For nodes with degree smaller than the critical degree, the pattern produces a positive impact on epidemic mitigation; conversely, for nodes with degree larger than the critical degree, the pattern produces a negative impact on epidemic mitigation. This paper is organized as follows. In Section 2, the metapopulation epidemic model with the time-varying human mobility pattern is suggested. In Section 3, by applying the MF (mean-field) approach, the theoretical analysis of the model is performed. Section 4 gives the simulation results. Conclusions and discussions are presented in Section 5. In addition, Appendix gives the method of calculating the critical node degree in details.

In this paper, we consider a total population of N individuals who diffuse in a heterogeneous network with V nodes, wherein each node represents a subpopulation. Each node i contains N i individuals. Two nodes are said to be neighbors with each other given that they share a link in this network. The individuals diffuse from their own node to their neighbor nodes. Here, we extend the RD process in Ref. [17] by incorporating the time variation of human mobility rates to describe the evolution process of the system. The individuals in a node first react and then diffuse at each time step. The classical SIS model is applied to the reaction process inside each node. Schematically, the infection dynamics can be identified by the following set of reaction equations:

The first subequation means that each susceptible (S) individual becomes infected with probability β when it comes into contact with an infected (I) individual. The second indicates that the infected nodes are cured and become susceptible again with probability µ.

At the end of the reaction, individuals begin to move to their neighbor nodes. Similar to Ref. [17] , here we also consider a random walk, but the mobility rates are related to the global epidemic information that can be easily obtained from mass media in today's information society. Let S j (t) and I j (t) denote the number of susceptible and infected individuals in node j at time step t, respectively. The mobility rates are defined as follows:

The quantity d I (t) (d S (t)) is the mobility rate of the infected (susceptible) individuals at time step t. Coefficient ω 1 (ω 2 ) represents the fraction of the infected (susceptible) individuals having awareness of an epidemic spreading. For example, ω 1 = ω 2 = 1 indicates that all individuals have awareness; ω 1 = ω 2 = 0 means that no one has awareness. Clearly, the mobility rates d I (t) and d S (t) are both time varying when ω 1 ̸ = 0 and ω 2 ̸ = 0. At the same time, they will decrease with the increase of ω 1 , ω 2 , and  j I j (t). This reflects the fact that the more people have awareness, and (or) the higher the epidemic prevalence is, the less people travel.

In order to theoretically describe the reaction-diffusion process of the system in heterogeneous networks with large degree fluctuations by applying the MF approach, we assume as in Refs. [17, 18, 21] that the nodes with same degree are statistically equivalent, and define relative densities ρ S,k (t) and ρ I,k (t) denoting the average number of susceptible and infected individuals in these nodes with degree k at time step t, respectively. Therefore, 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. Obviously,

represents the average number of individuals in these nodes with degree k.

According to the working rules of the RD process described in Section 2, we divide each time step into two substeps:

The reactions occur at t → t ′ and the diffusions happen at t ′ → t + 1. Furthermore, we assume that each susceptible individual may react with all the infected individuals in the same node [17] , and then the number

). By using the MF approach, the average number of new infected individuals can be represented as the mathematical expected value of random variable m, i.e., ⟨m⟩ = ρ S,

. At the same time, the number of new recovery individuals is µρ I,k (t).

Here, we define the reaction kernel [17, 18] as

Then, ρ I,k (t ′ ) and ρ S,k (t ′ ) can be expressed as follows:

For mobility dynamics of individuals starting at the end of reaction at time substep t ′ → t + 1, we assume homogeneous mobility along any given connection. From our human mobility pattern (see Eq.

(2)) defined in Section 2, the mobility rates along any given link from a node with degree k to another with degree k ′ will simply be equal to

So we have the following schematic description of the diffusion process:

where the term d I (t)ρ I,k (t ′ ) or d S (t)ρ S,k (t ′ ) identifies the number of infected or susceptible individuals moving out of the class nodes with degree k, and the term k

represents the infected or susceptible individuals arriving from neighbor nodes.

Combining the two subdynamics processes above (i.e., inserting Eq. (5) into Eq. (7)), and then calculating ρ I,k (t + 1) − ρ I,k (t) and ρ S,k (t + 1) − ρ S,k (t), we can easily get the following dynamical reaction rate equations at the MF level:

It is noted that some studies considered reaction and diffusion processes taking place simultaneously. That is, a given particle (individual) either reacts or diffuses to a neighbor node at the same time step (i.e., so-called continuous-time formulation of the reaction-diffusion process). The form of the corresponding dynamical equations can be found in Refs. [31] [32] [33] .

Clearly, Eq. (2) can be rewritten in the following form:

is the average density of infected individuals at time step t and ρ = N V denotes the total average density of all individuals. Moreover, we consider the case of uncorrelated network in which the conditional probability does not depend on the original node, i.e., P(k ′ /k) = k ′ P(k ′ )/⟨k⟩ [34, 35] . Substituting this equation and Eq. (9) into Eq. (8) , and then imposing the conditionsρ I,k (t) = 0 andρ S,k (t) = 0, one can obtain the stationary state equation of Eq. (8), taking the

where Γ =  k P(k)Γ k is the average of Γ k (Γ k has been given by Eq. (4) but with ρ I,k (t), ρ S,k (t) replaced by ρ I,k , ρ S,k respectively) over k. Multiplying Eq. (10a) by P(k) and summing over k, this yields the general result,

Inserting Eq. (11) into Eq. (10), a simple form of Eq. (10) can be written:

The case when ω 1 = ω 2 = 0 was discussed in Ref. [17] . In the following, we analyze epidemic dynamics in networks for ω 1 ̸ = 0 or/and ω 2 ̸ = 0.

First, we analyze the case when ω 1 = 0 (ω 2 ̸ = 0), which means that the mobility rate of infected individuals is unchanged, while the mobility rate of susceptible individuals is changed. In this case, Eq. (12a) takes the following simplified form:

If we insert the stationary form of Eq. (4) (i.e., Γ k = ρ S,k ρ I,k ) and Eq. (13) into Eq. (12b), we can obtain ρ S,k :

On the other hand, we can get from Eqs. (11) and (13) that

By substituting Eq. (14) into Eq. (15) we obtain a self-consistent equation for ρ I , denoted as ρ I = f (ρ I ). It is easy to understand that, if the first derivative value of the function f (ρ I ) is greater than 1 at the origin point, an epidemic will break out. By imposing the condition 

If ω 1 ̸ = 0 and ω 2 = 0, the mobility rate of infected individuals is changed, while the mobility rate of the susceptible individuals is unchanged. Similarly, we can get the same epidemic threshold as for the case ω 1 = 0 (ω 2 ̸ = 0). 

and then the standard quadratic equation form for ρ I,k (regarding ρ I as a parameter) is given by

Clearly, ρ I,k = 0 when ρ I = 0. In contrast, when ρ I ̸ = 0, we can figure out the positive root of Eq. (18), denoted here as

From Eqs. (11), (17) and (19), we have a self-consistent equation for ρ I :

It is not difficult to understand that, if the first derivative value of the function f (ρ I ) is greater than 1 at a point close to the origin, an epidemic will break out. In other words, a non-trivial solution of ρ I = f (ρ I ) exists if the inequality df (ρ I )

is satisfied. Then, we obtain the epidemic threshold:

Here, we point out that the result is obtained in the special case when ω 1 = ω 2 ̸ = 0 and it also holds for the general case ω 1 ̸ = ω 2 , which we will validate later by simulations.

From the above analysis, we can find that the stationary density of infected individuals ρ I (see Eqs. (15) and (20)) and the distribution of infected individuals in nodes with different degree ρ I,k (see Eq. (19) ) are both affected by our pattern.

However, it is also found that the epidemic threshold is independent of time variation of mobility rates, and is identical with the threshold in Ref. [17] . That is to say, the time-varying human mobility pattern does not alter the epidemic threshold. This feature may be traced back to the fact that the mobility rates are prevalence based.

We perform numerical Monte Carlo simulations on uncorrelated heterogeneous networks to validate the theoretical findings in Section 3. The networks are generated with the UCM (uncorrelated configuration model) [36] based on the Molloy-Reed algorithm [37] . Two networks are generated for these simulations: the first network has V = 1000 nodes with degree distribution P(k) ∼ k −3 (2 ≤ k ≤ V 1/2 ) characterized by the first moment ⟨k⟩ = 2.86 and the second moment ⟨k 2 ⟩ = 11.16, respectively. The second network has V = 2000 nodes with degree distribution P(k) ∼ k −2.5 (2 ≤ k ≤ V 1/2 ), and the first and second moments are ⟨k⟩ = 3.6 and ⟨k 2 ⟩ = 23.8, respectively. In order to reduce the influence of random factors on the simulation results, each simulation result shown in the figures is averaged over 20 independent realizations of different initial infected individuals. At the initial time step, all individuals are randomly distributed in the nodes. In Fig. 1 , we report the simulation results for the normalized stationary average infection density as a function of the population density in the network with ω 1 and ω 2 equal to 0 or 1. We can see from ⟨k 2 ⟩ = 0.73 (wherein the parameters are β = 0.1, µ = 0.1, ⟨k⟩ = 2.86, and ⟨k 2 ⟩ = 11.16) regardless of whether mobility rates are changed or not (i.e., ω 1 = ω 2 = 0 or ω 1 ̸ = 0 or (and) ω 2 ̸ = 0). So this validates the theoretical prediction with respect to the epidemic threshold. Fig. 1(b) shows the simulation results obtained for the second network with β = 0.1 and µ = 0.2. Obviously, the simulated values of the threshold are all almost equal to the analytical value ρ c = 1.08 calculated according to the given parameters (i.e., β = 0.1, µ = 0.2, ⟨k⟩ = 3.6, and ⟨k 2 ⟩ = 23.8). So, similarly, this simulation can support the conclusion that the time-varying human mobility pattern cannot alter the epidemic threshold.

In order to further validate the finding with respect to the epidemic threshold in a more general condition that 0 < ω 1 < 1 and 0 < ω 2 < 1, we report the simulated epidemic threshold ρ c as a function of parameters ω 1 and ω 2 in Fig. 2 . It can be found that the epidemic thresholds still approximate to 7.3 for the first network (see Fig. 2(a) ) and 1.08 for the second network (see Fig. 2(b) ), respectively, although they have some slight fluctuations. This further proves that the time-varying mobility pattern does not alter the epidemic threshold, regardless of whether ω 1 is equal to ω 2 or not.

To emphasize the impacts of the time variation of human mobility rates on epidemic behavior when ρ > ρ c , we compare the simulated results in two typical cases of ω 1 = ω 2 = 0 (not considering time variation of human mobility rates) and ω 1 = ω 2 = 1 (time-varying human mobility rates). In the following, for convenience, we denote the former case as the Case I pattern and the latter as the Case II pattern.

In Fig. 3 , we show the normalized stationary average infection density ρ I /ρ and relative difference in ρ I /ρ (see the insert) versus the total density of individuals ρ when ρ > ρ c . It can be seen that the stationary infection density in the Case II pattern is always less than that of the Case I pattern irrespective of the ρ value. This means that the time-varying mobility pattern can slightly mitigate the global epidemic spreading as a whole. Moreover, we also find that, in a certain range of ρ (e.g., 5 ≤ ρ ≤ 15 in Fig. 1(a) and 2 ≤ ρ ≤ 18 in Fig. 1(b) ), the value of ρ I /ρ in the Case II pattern is significantly lower than that in the Case I pattern, whereas such difference is very small if parameter ρ does not lie in this range. The main reason for the occurrence of this phenomenon can be explained as follows. On the one hand, when the value of ρ is too small, each node has on average very few individuals. In such a case, a mobility rate change can hardly alter the total average number of contacts. On the other hand, the larger the density of individuals ρ is, the more individuals move among the nodes, and the faster the infected individuals move to most nodes. Once most nodes have infected individuals, decreasing the mobility rate cannot effectively affect the epidemic prevalence in the stationary state.

Furthermore, we show the distribution of ρ S,k with respect to degree k in Fig. 4 . It is easy to see from Fig. 4 (a)-(d) that there always exists a common point at which the ρ S,k values are equal for the two patterns, regardless of the ρ value. Here, we denote the nodes' degree corresponding to the common point as the critical degree k c . For nodes with degree smaller/larger than k c , the relative density of susceptible individuals ρ S,k of the Case II pattern is higher/lower than that of Case I pattern. Moreover, it is noticed that the equation ρ S,k + ρ I,k = kρ/⟨k⟩ is always satisfied for these two cases. So we can derive that the relative density of infected individuals ρ I,k of the Case II pattern is lower/higher than that of the Case I pattern for nodes with degree smaller/larger than ρ c . This means that k c is a critical threshold with respect to the node degree separating the positive impacts from the negative impacts on epidemic mitigation. This phenomenon can be understood as follows. The nodes with higher degree have more individuals, and this is convenient for epidemic spreading inside those nodes. So the infection densities in the stationary state are higher than those in nodes with lower degree nodes. Compared with the Case I pattern, the Case II pattern can decrease the number of individuals moving among nodes. This will restrain, in some sense, the percolation of infected individuals from high infection density nodes to low ones.

In the following, we try to obtain the theoretical value of the threshold k c . In fact, from Eqs. (11) and (12), we can get the expressions of ρ S,k for the Case I pattern and the Case II pattern, respectively. Letting the two expressions be equal, one can derive k c as (see the Appendix for details)

where η 1 and η 2 are normalized average infected densities ρ I /ρ for the Case I pattern and the Case II pattern, respectively. From Fig. 2 , we know that η 2 is always less than η 1 and tends to η 1 when ρ is sufficiently large, which leads to a decrease of the k c value. In the extreme case that η 1 = η 2 , the minimum value of k c is derived: k min c = ⟨k 2 ⟩ ⟨k⟩ . So, we can get k c approximately by using k min c when ρ is sufficiently large (η 2 ≈ η 1 ). That is,

From the given network parameters, we have the theoretical value of k c = k min c = 3.9 for the first network. It can be seen from Fig. 4 that the simulated values of k c are all greater than 3.9. Moreover, these simulated values of k c approach 3.9 when ρ increases. This agrees well with the theoretical results.

In particular, we have ⟨k 2 ⟩ = ⟨k⟩ 2 for homogeneous networks; thus k c ≈ ⟨k⟩ from Eq. (23) . This means that our pattern has the same effect on the final infected density for all nodes in such networks. This deduction agrees with our intuitive judgment. 

When an epidemic breaks out, people usually make some response, such as reducing the amount of travel to lower the risk of being infected. In this paper, we focus on how the response of reducing mobility rates affects the metapopulation epidemic dynamics. We introduce a time-varying human mobility pattern in which human mobility rates are dynamically changed according to the global prevalence. The impact of the pattern on metapopulation epidemic dynamics is studied. It is found that the time variation of human mobility rates does not alter the epidemic threshold, but it can lower the final average density of infected individuals as a whole. More importantly, we have also derived a critical value with respect to the node degree. For nodes with degree smaller than the critical value, our pattern produces a positive effect on epidemic mitigation. Conversely, for nodes with degree larger than the critical value, it produces a negative effect on epidemic mitigation.

Our present study has considered the time variation factors in the pattern and have found some interesting results on epidemic dynamics. These results enable us to better understand epidemic dynamical behavior in heterogeneous social networks. Certainly, more complicated time-varying mobility patterns can be defined based on other more realistic mobility patterns such as recurrent mobility or traffic-dependent mobility. This deserves further investigation in future work.

For convenience, here, we let ρ 0 . and ρ 1 . denote the related quantities in the cases ω 1 = ω 2 = 0 and ω 1 = ω 2 = 1, respectively.

In the case ω 1 = ω 2 = 0, from Eq. (12) or results in Ref. [17] , we have where η 1 = ρ 0 I /ρ and η 2 = ρ 1 I /ρ.

The mathematical theory of infectious diseases

Epidemic spreading in scale-free networks

Epidemic dynamics in finite size scale-free networks

Behavior of susceptible-infected-susceptible epidemics on heterogeneous networks with saturation

Epidemic spreading on uncorrelated heterogenous networks with non-uniform transmission

Epidemic spreading on heterogeneous networks with identical infectivity

Epidemic dynamics on complex networks

An SIS model with infective medium on complex networks

Imperfect targeted immunization in scale-free networks

Thresholds for epidemic spreading in networks

Modelling of discrete-time SIS models with awareness interactions on degree-uncorrelated networks

A second-sound based, hyperbolic SIR model for high-diffusivity spread

The architecture of complex weighted networks

The worldwide air transportation network: anomalous centrality, community structure, and cities' global roles

Influence of dynamical condensation on epidemic spreading in scale-free networks

Reaction-diffusion processes and metapopulation models in heterogeneous networks

Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations

Natural human mobility patterns and spatial spread of infectious diseases

Modeling epidemics dynamics on heterogenous networks

Phase transitions in contagion processes mediated by recurrent mobility patterns

Heterogeneous length of stay of hosts movements and spatial epidemic spread

Effects of city-size heterogeneity on epidemic spreading in a metapopulation: a reaction-diffusion approach

Epidemic dynamics on an adaptive network

Epidemic reemergence in adaptive complex networks

The impacts of subsidy policies on vaccination decisions in contact networks

Modeling the influence of information on the coevolution of contact networks and the dynamics of infectious diseases

The impact of awareness on epidemic spreading in networks

Modeling human mobility responses to the large-scale spreading of infectious diseases

Forecast and control of epidemics in a globalized world

Continuous-time formulation of reaction-diffusion processes on heterogeneous metapopulations

Analysis and monte carlo simulations of a model for the spread of infectious diseases in heterogeneous metapopulations

Effects of diffusion rates on epidemic spreads in metapopulation networks

Dynamical and correlation properties of the internet

Evolution of Networks: From Biological Nets to the Internet and WWW

Generation of uncorrelated random scale-free networks

A critical point for random graphs with a given degree sequence

The authors would like to thank the editors and the reviewers for their constructive comments and valuable suggestions. Appendix. The method of calculating the critical degree k c