key: cord-0685696-xdlk8557
authors: Han, Xiao-Pu; Zhao, Zhi-Dan; Hadzibeganovic, Tarik; Wang, Bing-Hong
title: Epidemic spreading on hierarchical geographical networks with mobile agents
date: 2013-09-13
journal: Commun Nonlinear Sci Numer Simul
DOI: 10.1016/j.cnsns.2013.09.002
sha: e3bf49a9ff91d2428ce8f220a12dc6e8d7eaeff6
doc_id: 685696
cord_uid: xdlk8557

Hierarchical geographical traffic networks are critical for our understanding of scaling laws in human trajectories. Here, we investigate the susceptible-infected epidemic process evolving on hierarchical networks in which agents randomly walk along the edges and establish contacts in network nodes. We employ a metapopulation modeling framework that allows us to explore the contagion spread patterns in relation to multi-scale mobility behaviors. A series of computer simulations revealed that a shifted power-law-like negative relationship between the peak timing of epidemics [Formula: see text] and population density, and a logarithmic positive relationship between [Formula: see text] and the network size, can both be explained by the gradual enlargement of fluctuations in the spreading process. We employ a semi-analytical method to better understand the nature of these relationships and the role of pertinent demographic factors. Additionally, we provide a quantitative discussion of the efficiency of a border screening procedure in delaying epidemic outbreaks on hierarchical networks, yielding a rather limited feasibility of this mitigation strategy but also its non-trivial dependence on population density, infector detectability, and the diversity of the susceptible region. Our results suggest that the interplay between the human spatial dynamics, network topology, and demographic factors can have important consequences for the global spreading and control of infectious diseases. These findings provide novel insights into the combined effects of human mobility and the organization of geographical networks on spreading processes, with important implications for both epidemiological research and health policy.

The complexity of human mobility and interaction patterns has attracted much recent attention in a variety of disciplines, including statistical physics and complex systems science [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] . Ever since the influential work of Brockmann et al. [4] on the bank note dispersal, researchers have been intrigued by the revealed power-law distributions of human travel displacements [3, 5] . Meanwhile, such empirical studies have increasingly benefited from a variety of advanced data collection technologies (e.g. mobile phone, GPS, social-network sites etc. [15] ), reporting the presence of many abnormal properties of

Models that explicitly address the hierarchical organization of geographical networks are capable of mimicking realworld traffic systems in which agents can randomly percolate along the network edges and thereby take part in long-range, inter-layer travels. In one such model [9] , agents can generate not only power-law-like travel displacement distributions, but are also able to display a scaling behavior in the probability density of having traveled a certain distance at a certain time, which is in agreement with recent empirical observations [4] . In addition, it has been shown that inter-event or waiting time distributions also display non-Poisson statistics. We therefore argue here that an implementation of a hierarchical network in our present study can help us understand highly relevant aspects of real epidemics that are otherwise difficult to capture with models that enable only local behavior, i.e. mixing at only one or two scales [22] .

With respect to the infectious disease dynamics, our model is based on the SI spreading process [23, 24] , in which individuals can occupy one of the two possible states: 'susceptible', meaning they are not infected yet, and 'infected', meaning they already have the disease and can spread it to the susceptibles. In the present study, we examine the case of d ¼ 2 dimensional geographical hierarchical network with the following properties:

(i) Structure of the metapopulation network. We use the regular form of the hierarchical network reported in Ref. [9] . In our model, all units, or say subpopulations, are called cities, which can correspond to any kind of human settlements and not only to the real-world cities. These cities are organized in L layers of the network, where the first layer corresponds to the top network layer. A uniform 3-layer structure of the hierarchical geographical network is depicted in Fig. 1 .

We first set the city located in the center of a square field as the top layer (i.e. the first layer) of the network, and we evenly divide this field into K regular sub-regions (with a 3 Â 3 arrangement for K ¼ 9). We then set all the centrally located cities of the respective sub-regions to be the 2nd-layer cities. We note here that the 2nd-layer city in the central sub-region is also regarded as the same node of the 1st-layer city due to the overlap of the positions. Here, the sub-regions are called ''2ndlayer sub-regions''. Similarly, each of the 2nd-layer sub-regions are further evenly divided into K 3rd-layer sub-regions with 3rd-layer central cities. This process is iterated until Lth-layer cities are generated. For 2 6 n 6 L, there are K nÀ2 ðK À 1Þ nthlayer cities. The total number of cities in the network (the size of the network) is then defined as S ¼ K LÀ1 .

In the ðL À 1Þth-layer sub-region, each Lth-layer city is connected to the central city of the sub-region. These Lth-layer cities belonging to the same ðL À 1Þth-layer sub-region are also fully connected. Similarly, for 1 < n < L, the nth-layer cities in the same ðn À 1Þth-layer sub-region are also connected with the central city of the ðn À 1Þth-layer sub-region and are also fully connected to each other, as depicted in Fig. 1(b) and (c). The Fig. 1 illustrates an example with L ¼ 3, and K ¼ 9.

(ii) Random-walk of agents. The total number of agents in the model is M. At each time step, we randomly choose qM walkers to move along the edges of the hierarchical network. The value of q corresponds to the rate of mobile agents and is fixed in the model. In reality, the central city (hub) should have a greater attraction relative to smaller peripheral cities, which is represented in our model by a layer-dependent weight, w n ¼ w LÀn , where n denotes the layer and the ratio w P 1 is a free parameter. The probability that the walker will move along an edge is proportional to its weight.

For each move, the move length d is defined as the Euclidean distance from the city of a given movement initialization to a given target city. Here, the minimum distance between the two neighboring lowest-layer cities is set as one unit. For example, the black arrows in Fig. 1(d) illustrate a possible movement trajectory of an agent with the corresponding move lengths d 1 ¼ 2:236; d 2 ¼ 1:414; d 3 ¼ 6:708, and d 4 ¼ 1:0.

As w increases, the long-range travels to higher-layer cities become more frequent. Thus, larger w correlates with the higher heterogeneity or diversity (i.e attractiveness) of a city. The probability P i that an individual will move to a given neighboring city i, is the ratio of the weight of that neighboring city and the sum over the weights of all the connected neighboring cities surrounding the current city populated by the individual. (iii) Contact and infection of agents. At each time step, each agent randomly interacts with another agent in the same occupied city. If only one of the two agents is infected, another agent will then change his state from susceptible to infected with probability k.

In the initial series of simulations, we fix the following parameters: q ¼ 0:01; k ¼ 0:1; K ¼ 9; L ¼ 5, and M ¼ 10 5 (these are also the default parameter settings in all the following discussions, unless a given parameter value is differently specified). Starting from a random initial configuration, we first let the agents randomly move across the network for over 10 3 time steps to initialize the landscape. In the first time step of the spreading process (t ¼ 0), a randomly chosen agent is set as the first infector, while the others are set as susceptibles. Initially, the network takes a form of an 81 Â 81 lattice with L ¼ 5. In later simulations (Section 5), the network size and population density (as well as other parameters) are systematically manipulated to investigate the demographic effects on the spreading process.

As demonstrated in Ref. [9] , the move-length distribution of agents' mobility in the hierarchical geographical network obeys the power-law-like form PðdÞ $ d Àb . Fig. 2 (a) shows that when the value of w increases from 1.0 to 2.0, the powerlaw exponent, b, monotonically decreases from 2.0 to 0.74, which is in agreement with a whole range of empirical observations on human mobility patterns [3] [4] [5] . We also see that after a sufficient number of iterations, the majority of agents resides in higher-layer cities, which again depends on the value of the parameter w. Thus, as w increases, the long-range jumps to higher-layer cities become more frequent ( Fig. 2(b) ). Interestingly enough, we observe that intervals between the two consecutive contacts established among pairs of agents in our model also obey a power-law-like distribution with an almost unchangeable exponent À1:35 and slight fluctuations Hierarchical structure of the geographical network used in the presented study (K ¼ 9 and L ¼ 3), (a) shows the 3-dimensional representation of the planar hierarchical organization and the edges between the cities in different layers, (b) and (c) display the full-connected structure in the two sub-graphs of the network, (d) shows the corresponding distribution of cities in different layers of the network on the plane. The blue, red, and green circles are the cities in the 1st-, 2nd-, and 3rd-layer, respectively. The blue and light pink arrows denote the corresponding parts in different panels. The single topmost node is the first or highest layer, n ¼ 1. The black arrows show the possible movements of an agent. The panel (e) zooms in an individual node, showing an invasion of an infector (blue agent) and the local spreading within a city, whereas the grey arrow denotes the interaction and contagion spreading between the infector and a previously healthy agent (now infected). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.) ( Fig. 2(c) ). This observation is very close to the temporal properties of natural, real-world social contacts [25] that are clearly different from the behaviors generated by a traditional Lévy-walk system [26, 27] .

As in an SI model, the spreading velocity is represented by the growth of the total infection rate RðtÞ ¼ M I ðtÞ=M and the rate of new infected cases IðtÞ ¼ RðtÞ À Rðt À 1Þ, where M I ðtÞ is the number of infectors at time step t. As shown in Fig. 3 (a) and (b), the spreading velocity is monotonically accelerated with the increasing value of the layer-dependent weight w. When the ratio of layer-dependent weight w is held constant (Fig. 3(d) and (e)), we find that the contagion spreads first from its source location to the closest higher-layer city (directly linked to the contagion source), and then diffuses to the surrounding same-layer cities of the initially affected higher-layer city via available long-range links. Cities in the lower layers of the network are the last ones to become infected. These results suggest that higher-layer cities are at a higher risk for fastgrowing epidemic, and long-range travelers may play the crucial role in this spreading process. Moreover, it is certainly not by chance that hub-like nodes with great infectivity in epidemic contact networks are actually dubbed the superspreaders in the epidemiological literature [28] .

The layer-by-layer spreading observed in our simulations implies that the spreading process in the hierarchical network is sensitive to the source location of the epidemic onset. Setting the first infector as static in an L 0 -th layer city at t ¼ 0, as shown in Fig. 3 (f) and (g), the disease transmission becomes faster when the source is located in higher-layer cities (i.e. when L 0 is small), and especially at larger w. This confirms previous findings suggesting that it is of utmost importance to set the disease monitoring focus on the big cities and transport hubs with greater attraction for long-range travelers.

To compare the dynamics of disease transmission emerging on our hierarchical geographical network with that generated by the Lévy flight model, we additionally simulated the spreading process with Lévy-walking agents in same-sized, 81 Â 81 lattices. To enable comparison, the move-length distribution originating from the Levy-flight model was the same as the one found in the hierarchical network structure (see Appendix A for details). However, relative to the hierarchical network scenario, the movements of agents in the Lévy-flight scenario were not limited to those along the edges and contact intervals did not show scaling properties.

As depicted in Fig. 4 (a) and (b), the spreading velocity in the hierarchical network structure is somewhat slower than that observed in the scenario with the Lévy-flight model (when w is close to 1:0). When the value of w > 1, the growth rate of infections in the hierarchical structure is higher in the early stages of the spreading process (relative to the scenario with Lévy-flight agents), but slower in the later stages, such that the overall relaxation time from the start of the epidemic to the full-infected state is longer than that found in the Lévy-flight model, as illustrated in Fig. 4 (c) and (d).

Previous research on the role of socio-demographic factors in disease spreading yielded many conflicting results [29] [30] [31] . In fact, relatively little is known about the relationship between the peak timing of epidemics and socio-demographic factors such as population density and city size (e.g. epidemic spreading differences between rural and urban areas). In particular, the effects of population size and density have not been comprehensively studied in the context of hierarchical spreading behaviors, driven by large, hub-like population centers [30] .

In Fig. 5 , we see that population size M significantly affects both the height and timing of IðtÞ. Moreover, we find that the population size affects only the peak timing of epidemics in lower-layer cities (see Fig. 5 (c)), while this effect vanishes in higher-layer cities, resulting thus in rather synchronized epidemic peaks across a variety of population sizes. Furthermore, we observe that the population size effect becomes obvious only in sufficiently large networks (see the inset in Fig. 5(c) ).

Similar effects of population density were observed in our model (Fig. 6) , where the average population density hm s i ¼ M=S. However, one important difference is that the population density plays a role in epidemic timing even in a wider range of network layers (however, the effect is again stronger in the lower layers; see Fig. 6 (a), revealing a shiftedpower-function-like relationship). Next, as we were wondering why the spreading velocity negatively correlates with the population density, we found that this particular relationship emerges from the gradual fluctuation enlargement effect on the export time of infectors (see Appendix B for details). When the population density is large enough, then the population size in each city and the amount of travelers moving between different cities is also large, which consequentially results in much smoother contagion spreading patterns across cities.

However, the case with a small population density is clearly different. Here, it is very likely that in the long term, no travelers will move from small-sized to larger cities simply because of the random selection of travelers from a small-sized population pool. Thus, a given small-sized city cannot immediately export infectors even though it is fully infected. This phenomenon will significantly delay the spreading process, which in turn will be gradually enlarged by the hierarchical organization of the system, ultimately resulting in the negative correlation between the spreading velocity and population density. This further implies that the actual population distribution deeply impacts the spreading process: The hierarchical system with many small and loose settlements will result in slower spreading velocity and higher unpredictability. The detailed discussions of this phenomenon based on the mean field approach can be found in the Appendix B.

In addition, we discovered that the peak timing of epidemics scales logarithmically with network size (and thus linearly with L), showing a regularity across several 'urban systems' in our model, i.e. across a range of population densities (see the inset in Fig. 6 ) and across different time-scales. Our semi-analytical approach revealed that this logarithmic growth results from the gradual fluctuation enlargement effect on the export time of infectors (see again the Appendix B).

Finally, our simulations revealed that the density effect is weakened for the condition with a large geographical diversity (large w), as shown in Fig. 6(b) . The reason is that, in this condition, many agents gather in the higher layer cities where the spreading is mainly driven by the local infections within a city, and thus the gradual enlargement effect no longer plays an important role.

One of the most prominent features of the hierarchical network structure is that a few higher-layer cities are usually connected by several long-range edges. Since higher-layer cities are also facing higher risks for fast-growing infectious outbreaks, it is the question how to select an optimal countermeasure that would effectively delay the timing of local outbreaks in higher-layer cities in order to prevent further spreading to less endangered, peripheral areas.

In the influenza pandemic of 1918, strict border control measures delayed the disease invasion for more than one year in several Pacific countries [32] . This entry screening policy has also played a major role in the combat against the Severe Acute Respiratory Syndrome (SARS) in 2003. However, the full extent to which border control measures actually contribute to realworld epidemic preventions is still rather unknown. While some reports [33] [34] [35] [36] provided rather pessimistic estimations on the effectiveness of border control, recent evaluations of the entry screening policies adopted by different nations during the 2009 H1N1 pandemic [37] demonstrated that long-range travel control significantly delayed disease transmission, relative to countries which did not employ a strategy for screening of incoming passengers. Due to such conflicting reports, more studies on the relevance of border control measures are needed.

Although the time delay of an epidemic outbreak (such as that resulting from entry screening of long-range travelers) has been explored in a few recent studies [36, 37] , the impacts of various intervention strategies on the height and timing of peak infections, especially in the context of more realistic epidemic scenarios, have remained rather understudied. In this section, we will therefore investigate the efficiency of long-range travel screening as implemented in on our present model. To this effect, we first set the threshold d c for the move length d of agents that will be subject to a border screening procedure.

More specifically, all agents with the move length d who at a given time step satisfy the condition d P d c will be screened. Here we define the cost of screening C as the average number of agents traveling along the edges longer than the threshold length d c : C ¼ hm c i=M. Obviously, the cost C correlates negatively with d c : When d c is large, only a few higher-layer cities will be screened. On the other hand, small d c implies that much more cities will be involved in the screening procedure. Since there are more long-range travelers for larger w, the cost C also increases with the value of w. Fig. 7(a) shows the relationship between C and d c for different values of w (in the inset). The bumps visible in Fig. 7 (a) are due to the bumpy move length distribution PðdÞ (shown in Fig. 2(a) ). Under this border screening regime, each infected individual participating in a long-range travel (with move length not smaller than d c ) will subsequently be quarantined with probability p, or it will keep its infectivity with probability 1 À p; thus, p can serve here as the measure of efficiency of border control.

We now focus our discussion on the peak timing ratio between the situation with (s) and without (s 0 ) border control. The higher the ratio s=s 0 the longer the relative delay of the peak timing of epidemic spreading. When compared against the scenario without border control, our simulations show that the efficiency of the screening of long-range travelers is deeply affected by several socio-demographic factors such as city size and population density. First, we find that border control of long-range travelers would be a more efficient measure when the diversity of cities is rather small (i.e. small w). As shown in Fig. 7(b) , the ratio s=s 0 rapidly decreases with the growth of w. As expected, the border control of long-range travelers was found to be more efficient for small w, which corresponds to more homogeneous cities and a lower influx of passengers from distant regions. As mentioned previously, when w is large, only a few higher-layer cities can manage to collect the largest fractions of the simulated population, such that a swift local outbreak in higher-layer cities has much more contribution to the overall spread of infections in the rest of the network. Moreover, since larger w implies a more costly border control for the same threshold d c (Fig. 7(a) ), the feasibility of this control measure for the condition with the larger city diversity is further reduced. Additionally, the relationship between the ratio s=s 0 and the cost of border screening C is logarithmic for w ¼ 1 (the inset in Fig. 7) , indicating that the improvement of the border screening efficiency would be rather limited if one would only enlarge the screening range. Thus, the above discussed results regarding the efficiency of border control are rather pessimistic, however, the observed sensitivity of the epidemic delay to the probability p of isolating an infected agent shows another picture. The growth of the ratio s=s 0 with p is quite fast and surprisingly shows a super-exponential form (Fig. 7(c) ), suggesting that the efficiency of the proposed border control mechanism can largely be improved if the screening procedure is more dedicated to increasing the probability of infector detection and isolation.

Moreover, in this context, we further investigated the impact of population density on the efficiency of border control. Fig. 7 (c) shows that border screening is more efficient for the condition with a lower average population density. The specific relationship between the epidemic delay and population density is plotted in the inset of Fig. 7(c) , and is shift-power-law like for the condition with smaller w. When w is large, the effect of population density is weakened, because in this case, the rise of local infections within a city weakens the gradual enlargement effect. Thus, we find here that besides w, the population density also has an additional impact on the outcome of a control measure, which would be more efficient in hierarchical systems with many small and loose settlements.

The hierarchical structure of our epidemic model contains two relevant aspects. Firstly, our hierarchical geographical network mimics the real-world traffic systems, and secondly, as indicated in Ref. [9] , random walks in such hierarchical networks can naturally generate scaling mobility patterns, which are quantitatively in agreement with recent empirical studies on human mobility. Consequentially, the geographical hierarchical network presented in our model not only provides a novel framework for modeling the interaction among various subpopulations, but it also introduces the effects of uniquely human patterns of motion on the contagion spread processes.

Taken together, our simple metapopulation model was able to reproduce some of the key features that are typically observed in more complex natural systems with spreading processes. We first investigated contagion spread patterns on hierarchical networks and then compared the observed behavior against the scenario generated by the Lévy-flight model. In the early stages of disease transmission, higher-layer cities are at a higher risk for a fast-growing outbreak, which is mainly caused by the infected long-range travelers. This in turn may necessitate top-down procedures such as exit and entry screenings at international borders [38] which were quantitatively studied in our paper.

However, we note here that larger cities and transport hubs with greater attraction for long-range travelers are not necessarily the key players in the spreading process. For example, Kitsak et al. [39] have shown that quite differently from the common belief, the most central or the most connected nodes in a network are not the best spreaders of a contagion. Instead, it is those individuals positioned within the core of a given network who are actually responsible for the largest-scale spreading of a disease [39] .

In the present model, the links between the nodes of the investigated hierarchical network were undirected. However, the directedness of links [40, 41] was recently reported to play an important role in the spreading process [42] , resulting in behaviors that are distinct from those observed on undirected networks [43, 42] . Generally, studies of epidemic spreading on directed networks are relatively rare [44] , even though directed flow is known as an essential characteristic of many asymmetric real-world processes such as contagion spread in trade networks [45] , in which the availability of a given network link does not guarantee the reciprocal connection. Future investigations of epidemics on hierarchical networks should not ignore these issues. For example, in a hierarchically organized system, some starting peripheral nodes could have only out-links, while the majority of others could allow for mobility in both directions; however, some inter-layer transitions or quarantined areas could be specified with in-links only.

For geographical regions with greater variations in socio-demographic factors (e.g. larger amount of city layers) and sufficiently large sizes, we found that population density and city size become good predictors of the epidemic timing. Moreover, we discovered a shifted power-law-like negative relationship between the peak timing of epidemics s 0 and population density hm s i, and a logarithmic positive relationship between s 0 and the network size S. This may further indicate that just as many other natural phenomena, also epidemic spread patterns may have some universal properties that are shared by all world cities and sustained across different cultures and times [46] . In addition, our mean-field approximation showed a gradual enlargement effect in the spreading process on a hierarchical network, and revealed that the smooth-spreading assumption cannot always provide an appropriate description of epidemic spreading in complex networks, especially not in hierarchically organized systems.

Finally, the efficiency of the employed border control strategy and its cost are quantitatively studied. At first, the reduced screening efficiency for situations with large city diversity and the logarithmic relationship between the screening efficiency and cost of border control confirm the previously reported pessimistic outcomes of travel restrictions (see e.g. [35] for travel control measures related to the global spread of 2009 H1N1 pandemic, exploring mobility restrictions that differ with respect to their magnitude and timing). On the other hand, the observed super-exponential growth of s=s 0 with p is an indication of a possible improvement of this control measure e.g. by means of enhancing the detectability of infectors at each entry during the screening procedure.

However, we point out here that the ideal intervention, especially in the early stages of an outbreak, may depend on the actual aim(s) of an intervention [47] , as some intervention goals may have a multitude of roughly equally effective strategies, while others may require only one specific action. This, in turn, indicates the relevance of more precise descriptions of policy objectives when planning a given mitigation strategy. Our study finds that the effectiveness of border control is sensitive to the specific pattern of cities and the associated population distribution. This in turn suggests that the conclusions of previous studies related to special regions [33] [34] [35] could be carefully promoted to other situations. Generally speaking, the underdeveloped and developing countries usually have more homogeneous cities, and some countries with low population density typically have loose population configurations. Accordingly, border control measures may be of higher relevance for arresting epidemics in these countries.

In summary, based on the metapopulation modeling framework, we have introduced a novel disease transmission model, derived from the combined study of real-world human mobility patterns and the multi-scale spread of infectious diseases, which can be employed to investigate the effects of early epidemic countermeasures on the height and timing of peak infections in hierarchically organized, complex geographical networks. Our results suggest that the interplay between the dynamics of human spatial behavior and network topology can have important consequences for the global spreading of infectious diseases. We therefore hope that our results can motivate further model extensions that should provide broader insights into the combined effects of human mobility and the organization of geographical networks on spreading processes, with relevant outcomes for both epidemiological studies and health policy. m I ðt þ 1Þ À m I ðtÞ ¼ k hm s i ðhm s i À m I ðtÞÞm I ðtÞ ð B:1Þ

We consider a temporal mean field assumption that ignores the randomness of the timing of infection events within a city; the left-hand side of Eq. (B.1) can be approximately replaced by the differential form dm I ðtÞ dt . According to the initial condition m I ð0Þ ¼ 1, we can obtain the approximate function of the m I ðtÞ, m I ðtÞ % hm s ie kt e kt þ hm s i À 1

: ðB:2Þ

The probability that the source city exports an infector to another city for the first time at some time point t is proportional to the current total number of infectors m I ðtÞ, and we thus obtain this probability immediately from:

where P tÀ1 t 0 ¼0 X 1 ðt 0 Þ represents the probability that the city exports an infector before the time point t. We find it difficult to get the analytical formula of X 1 ðtÞ. From Eqs. (B.2) and (B.3), the numerical result of X 1 ðtÞ is shown in Fig. B.8(a) , which is unimodally varying with time. For smaller population densities m s , due to the slower decay of X 1 ðtÞ (Fig. B.8(a) ), the width of the peak is wider and the average value of the export time of infectors hti is slightly larger, which is mainly caused by stronger fluctuations of export time at smaller m s .

After the export of an infector from the source city, the second connected city will become infected. Since the spreading in the second city starts from the time at which the export from the source city has been carried out, and since the local spreading process also obeys Eq. (B.1), the probability X 2 ðtÞ that the second infected city further exports an infector for the first time at a time point t depends on X 1 ðtÞ,

where X 1 ðt À t 0 Þ and X 1 ðt 0 Þ respectively represent the probability that the initially infected city (the source city) exports an infector for the very first time at time t À t 0 , and that the second city exports an infector for the first time after t 0 import steps.

Since X 1 ðtÞ is unimodal, X 2 ðtÞ is also having a single-peak, and due to the slower decay of X 1 ðtÞ for small m s and the corresponding larger average export time, the peak-time of X 2 ðtÞ will be delayed and the peak-time difference for different m s will also be enlarged.

Similarly, for each additional cross-city spreading, a new local spreading process is added, and the X for the following infected cities in the spreading chain can be written directly where n s is the rank of the city in the spreading chain of cities. The difference of the peak-time will be enlarged along with the growth of n s . Defining the peak time difference of X 2 ðtÞ and X 1 ðtÞ as Ds, and since the enlargement of this difference is mainly generated by the local spreading process, from Eqs. (B.4) and (B.5), the peak-time of X ns ðtÞ would be semi-linearly correlated with the rank n s in the spreading chain.

The numerical solutions of Eqs. (B.2)-(B.5) confirmed the cascading enlargement effect of the fluctuations of the interregional spreading process (see Fig. B.8 ). In the hierarchical structure, the length of the spreading chain is generally proportional to L; for example, in our model, if the source city is in the bottom layer, the lowest possible minimum length of the spreading chain is 2L À 1, so the peak time s p for n s ¼ 2L À 1 is generally proportional to s 0 . As shown in the inset of Fig. B.8(e) , the linear dependence of the peak time s p of X ns ðtÞ vs n s is observed in the numerical solution. Because the size of the network is M LÀ1 , the logarithmic relationship s 0 vs. L (the inset in Fig. 6(a) ) is easily understandable from the linear form of s p . Also, as shown in Fig. B.8(e) , the curve of s p vs. hm s i also shows a shifted-power-function-like form, which is generally in agreement with the simulation results of our model.

Similarly, the fluctuations in the local contagion spreading also showed the cascading enlargement effect, but they were ignored in the above mean field discussions. This is why the fitting slope in the inset of Fig. B .8(e) is smaller than that in the inset of Fig. 6(a) . Above all, both the shifted-power-function-like negative relationship between s 0 and the population density hm s i and the logarithmic positive relationship between s 0 and S can be explained by the gradual enlargement of the fluctuations in the spreading process. This analysis additionally implies that the smooth-spreading assumption cannot always provide a meaningful description of the spreading process, in particular in hierarchically organized networks. Furthermore, the observed enlargement of the width of the curve of X with larger n s and smaller hm s i implies that the spreading process has less predictability and more variability in the hierarchical organization [48] , especially in networks with more layers or many small and loose settlements.

A universal model for mobility and migration patterns

A review of power laws in real life phenomena

Understanding individual human mobility patterns

The scaling laws of human travel

Characterizing the human mobility patterns in a large street network

Limits of predictability in human mobility

A Lévy flight for light

Modelling the scaling properties of human mobility

Origin of the scaling law in human mobility: Hierarchy of traffic systems

Understanding the spreading patterns of mobile phone viruses

Impact of travel patterns on epidemic dynamics in heterogeneous spatial metapopulation networks

Natural human mobility patterns and spatial spread of infectious diseases

Recurrent host mobility in spatial epidemics: beyond reaction-diffusion

Phase transitions in contagion processes mediated by recurrent mobility patterns

Digital epidemiology

Random walks and flights over connected graphs and complex networks

Modelling dynamical processes in complex socio-technical systems

The relationship between human behavior and the process of epidemic spreading in a real social network

Forecast and control of epidemics in a globalized world

Modeling the adoption of innovations in the presence of geographic and media influences

Critical paths in a metapopulation model of H1N1: efficiently delaying influenza spreading through flight cancellation

Multiscale, resurgent epidemics in a hierarchical metapopulation model

Dynamical patterns of epidemic outbreaks in complex heterogeneous networks

Behaviors of susceptible-infected epidemics on scale-free networks with identical infectivity

The origin of bursts and heavy tails in human dynamics

Impact of non-Poissonian activity patterns on spreading processes

Impact of human activity patterns on the dynamics of information diffusion

Are SARS superspreaders cloud adults?

The effect of public health measures on the 1918 influenza pandemic in U.S. Cities

The 1918-1919 influenza pandemic in England and Wales: Spatial patterns in transmissibility and mortality impact

Estimation of potential global pandemic influenza mortality on the basis of vital registry data from the 1918-20 pandemic: a quantitative analysis

Border control measures in the influenza pandemic plans of six South Pacific nations: a critical review

Strategies for mitigating an influenza pandemic

Effects of internal border control on spread of pandemic influenza

Human mobility networks, travel restrictions, and the global spread of 2009 H1N1 pandemic

A simple explanation for the low impact of border control as a countermeasure to the spread of an infectious disease

Entry screening to delay local transmission of 2009 pandemic influenza A (H1N1)

Dynamics of an SIQS epidemic model with transport-related infection and exit-entry screenings

Identification of influential spreaders in complex networks

Cooperation in the snowdrift game on directed small-world networks under self-questioning and noisy conditions

Formation control for second-order multi-agent systems with time-varying delays under directed topology

Epidemic spreading in annealed directed networks: susceptible-infected-susceptible model and contact process

Evolution of ethnocentrism on undirected and directed Barabási-Albert networks

Predicting epidemics on directed contact networks

Disease spread in small-size directed trade networks: the role of hierarchical categories

Growth, innovation, scaling, and the pace of life in cities

Mitigation strategies for pandemic influenza A: balancing conflicting policy objectives

Epidemic variability in hierarchical geographical networks with human activity patterns

To ensure comparability with the dynamics of disease transmission on our hierarchical geographical network, the Lévyflight model employed in our simulations had following characteristics:(i) The size of the network (81 Â 81 lattices), the number of agents M, and the algorithms governing the local interactions and the spreading behavior among agents within a city, are identical to those of the hierarchical network scenario. (ii) Differently from the hierarchical network scenario, agents can directly jump between any pair of nodes in the Lévyflight scenario. The move length of each move is here again defined as the Euclidean distance from the starting to the target city, as in the previously discussed hierarchical network structure. (iii) The move-length distribution of the Lévy-flight model was set to accurately correspond to the move-length distribution of the hierarchical network scenario. To this effect, we first create a data set containing a record of a large number of move lengths in the hierarchical network model with given parameter settings. In the Lévy-flight model, for each individual movement (starting e.g. with city A), we then randomly pick a move length (e.g. d ¼ d 0 ) from the hierarchical network dataset and we set it as the length of this move. Next, for the target city of this movement we randomly select one city from the set of cities having exactly the distance d ¼ d 0 from the city A. (iv) As in the hierarchical network scenario, agents cannot move across the boundary of the lattice.

For the source city of the epidemic, assuming the population density is hm s i and the spreading is smooth enough, it is easy to write the equation of the growth of the total number of infectors m I in the city: