Time-varying human mobility patterns with metapopulation epidemic dynamics.

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 k c . For nodes with degree smaller than k c , the pattern produces a positive impact on epidemic mitigation; conversely, for nodes with degree larger than k c , the pattern produces a negative impact on epidemic mitigation.

Introduction

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.

Model description

In this paper, we consider a total population of individuals who diffuse in a heterogeneous network with nodes, wherein each node represents a subpopulation. Each node contains 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 individual becomes infected with probability when it comes into contact with an infected 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 and denote the number of susceptible and infected individuals in node at time step , respectively. The mobility rates are defined as follows: The quantity is the mobility rate of the infected (susceptible) individuals at time step . Coefficient represents the fraction of the infected (susceptible) individuals having awareness of an epidemic spreading. For example, indicates that all individuals have awareness; means that no one has awareness. Clearly, the mobility rates and are both time varying when and . At the same time, they will decrease with the increase of , and . This reflects the fact that the more people have awareness, and (or) the higher the epidemic prevalence is, the less people travel.

Theoretical analysis

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 and denoting the average number of susceptible and infected individuals in these nodes with degree at time step , respectively. Therefore, we have where is the number of nodes with degree and sums run over all nodes having degree equal to . Obviously, represents the average number of individuals in these nodes with degree .

According to the working rules of the RD process described in Section  2, we divide each time step into two substeps: and . The reactions occur at and the diffusions happen at . Furthermore, we assume that each susceptible individual may react with all the infected individuals in the same node  [17], and then the number of new infected individuals (random variable) in nodes with degree at time substep will follow the binomial distribution . By using the MF approach, the average number of new infected individuals can be represented as the mathematical expected value of random variable , i.e., , which is approximately equal to . At the same time, the number of new recovery individuals is . Here, we define the reaction kernel  [17], [18] as Then, and can be expressed as follows:

For mobility dynamics of individuals starting at the end of reaction at time substep , 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 to another with degree will simply be equal to So we have the following schematic description of the diffusion process: where the term or identifies the number of infected or susceptible individuals moving out of the class nodes with degree , and the term or 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 and , 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: where is the average density of infected individuals at time step and 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.,    [34], [35]. Substituting this equation and Eq. (9) into Eq. (8), and then imposing the conditions and , one can obtain the stationary state equation of Eq. (8), taking the form where is the average of ( has been given by Eq. (4) but with replaced by respectively) over . Multiplying Eq. (10a) by and summing over , this yields the general result,

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

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

Case 1: or

First, we analyze the case when , 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., ) and Eq. (13) into Eq. (12b), we can obtain :

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

By substituting Eq. (14) into Eq. (15) we obtain a self-consistent equation for , denoted as . It is easy to understand that, if the first derivative value of the function is greater than 1 at the origin point, an epidemic will break out. By imposing the condition , we can derive the epidemic threshold with respect to the total density of individuals in the network:

If and , 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   .

Case 2: and

In this case, the mobility rates of infected and susceptible individuals are both changed with the epidemic prevalence. From Eq. (12) combined with , we have and then the standard quadratic equation form for (regarding as a parameter) is given by with

Clearly, when . In contrast, when , we can figure out the positive root of Eq. (18), denoted here as From Eqs. (11), (17), (19), we have a self-consistent equation for :

It is not difficult to understand that, if the first derivative value of the function is greater than 1 at a point close to the origin, an epidemic will break out. In other words, a non-trivial solution of exists if the inequality is satisfied. Then, we obtain the epidemic threshold:

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

From the above analysis, we can find that the stationary density of infected individuals (see Eqs. (15), (20)) and the distribution of infected individuals in nodes with different degree (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.

Simulation results

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 nodes with degree distribution characterized by the first moment and the second moment , respectively. The second network has nodes with degree distribution , and the first and second moments are and , 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 and equal to 0 or 1. We can see from Fig. 1(a) (simulation results obtained for the first network) that the simulated values of the threshold all agree approximately with the analytical value (wherein the parameters are , and ) regardless of whether mobility rates are changed or not (i.e.,  or or (and) ). 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 and . Obviously, the simulated values of the threshold are all almost equal to the analytical value calculated according to the given parameters (i.e.,  , and ). So, similarly, this simulation can support the conclusion that the time-varying human mobility pattern cannot alter the epidemic threshold.

Fig. 1

The normalized stationary average infection density versus the total density of individuals with and equal to 0 or 1. The red arrows indicate the theoretical value of the epidemic threshold . (a) The first network with and . (b) The second network with and .

In order to further validate the finding with respect to the epidemic threshold in a more general condition that and , we report the simulated epidemic threshold as a function of parameters and 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 is equal to or not.

Fig. 2

The epidemic threshold (i.e., the critical value of the average population density) as a function of and . The color indicates the simulated values of . The red arrows indicate the theoretical value of . (a) The first network with and . (b) The second network with and .

To emphasize the impacts of the time variation of human mobility rates on epidemic behavior when , we compare the simulated results in two typical cases of (not considering time variation of human mobility rates) and (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 and relative difference in (see the insert) versus the total density of individuals when . 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., in Fig. 1(a) and in Fig. 1(b)), the value of 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.

Fig. 3

The impacts of mobility patterns on the infection density. The figure compares the normalized stationary average infection density between the Case I pattern and the Case II pattern versus the total density of individuals when . The insert shows the relative difference in of the two patterns, defined as with and denoting the values of in the Case I pattern and the Case II pattern, respectively. (a) The first network with and . (b) The second network with and .

Furthermore, we show the distribution of with respect to degree in Fig. 4 . It is easy to see from Fig. 4(a)–(d) that there always exists a common point at which the 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 . For nodes with degree smaller/larger than , the relative density of susceptible individuals of the Case II pattern is higher/lower than that of Case I pattern. Moreover, it is noticed that the equation is always satisfied for these two cases. So we can derive that the relative density of infected individuals of the Case II pattern is lower/higher than that of the Case I pattern for nodes with degree smaller/larger than . This means that 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.

Fig. 4

Stationary density as a function of degree in the Case I pattern or the Case II pattern for different . The simulations are based on the first network with and . (a), (b), (c), and (d) denote the cases , and 20 respectively.

In the following, we try to obtain the theoretical value of the threshold . In fact, from Eqs. (11), (11) and (12), we can get the expressions of for the Case I pattern and the Case II pattern, respectively. Letting the two expressions be equal, one can derive as (see the Appendix for details) where and are normalized average infected densities for the Case I pattern and the Case II pattern, respectively. From Fig. 2, we know that is always less than and tends to when is sufficiently large, which leads to a decrease of the value. In the extreme case that , the minimum value of is derived: . So, we can get approximately by using when is sufficiently large . That is,

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

In particular, we have for homogeneous networks; thus 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.

Conclusions and discussions

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.