Epidemic spreading on hierarchical geographical networks with mobile agents.

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 τ 0 and population density, and a logarithmic positive relationship between τ 0 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.

Introduction

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 human motion, including ultraslow diffusion [3], [4] and strikingly high predictability [6]. It is these abnormal properties that actually shape the difference between the natural human movements and the so-called Lévy flights [7], [16].

Several studies have been conducted so far to address the underlying mechanisms of the discovered abnormal properties of human spatial behavior. For the urban and intercity travels, Song et al. investigated the impact of exploration and preferential return on mobility patterns [8]. For the long-range travels, Han et al. [9] demonstrated that hierarchical geographical organization of traffic systems plays a crucial role in the emergence of scaling in human trajectories.

An important, but less studied issue in this context, is the impact of the abnormal features of human motion on the dynamical evolution of spreading processes. For example, based on the mobile phone user data published in Ref. [3], Wang et al. analyzed the spreading of a mobile-phone-virus outbreak [10]. Ni and Weng [11] investigated the effects of heterogeneous spatial properties in metapopulation networks, and more recently, Belik et al. [12], [13] reported on the dynamics of a spreading process under bidirectional mobility, that was further investigated by Balcan and Vespignani [14]. All of these studies converge on the finding that spatially constrained mobility networks play a crucial role in shaping the key features of spreading processes.

While the focus of most previous studies has largely been on the local temporal development of diseases and epidemics [17], [18], their geographical, multi-scale circulation is less well understood [19], [20]. In the present paper, we propose a metapopulation model [21] that allows us to explore the contagion spread patterns in relation to multi-scale mobility behaviors and the organization of traffic systems embedded in realistic geography. Thus, we assume here that besides local, close-contact-driven and smaller outbreaks, it is the movement of individuals among the many different scales of a large geographical network that is actually essential to the spread of a global epidemic.

A previously reported model [9] demonstrated that a number of realistic mobility patterns can naturally be generated by random walks on a hierarchical geographical network, mimicking the properties of real-world traffic systems. Here, we directly implement a hierarchical geographical network with randomly-walking agents into the susceptible-infected (SI) epidemic model, combining thereby the study of epidemic dynamics with human-like mobility behaviors and the real-world-like organization of traffic systems (Section 2).

We first investigate the elementary properties of the generalized SI spreading process on hierarchical geographical networks and then compare the obtained epidemic patterns against the scenario generated with the Lévy flight model (Section 3; details of the employed Lévy flight model are further presented in Appendix A). We then systematically investigate how population density and city size affect the epidemic curves in different layers of a hierarchically organized geographical network (Section 4). A semi-analytical method is further introduced to better understand the importance of pertinent demographic factors (Appendix B). Finally, we study the effects of a border control measure on the timing and height of epidemic peak which is one of the main aspects of most mitigation strategies (Section 5).

The model

Models that explicitly address the hierarchical organization of geographical networks are capable of mimicking real-world 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 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 .

Fig. 1

Hierarchical structure of the geographical network used in the presented study ( and ), (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, . 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.)

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 arrangement for ). 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 “2nd-layer 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 , there are nth-layer cities. The total number of cities in the network (the size of the network) is then defined as .

In the 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 th-layer sub-region are also fully connected. Similarly, for , the nth-layer cities in the same th-layer sub-region are also connected with the central city of the 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 , and .

(ii) Random-walk of agents. The total number of agents in the model is M. At each time step, we randomly choose walkers to move along the edges of the hierarchical network. The value of 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, , where n denotes the layer and the ratio 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 , and .

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 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 .

In the initial series of simulations, we fix the following parameters: , and (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 time steps to initialize the landscape. In the first time step of the spreading process (), 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 lattice with . 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.

Simulation results

As demonstrated in Ref. [9], the move-length distribution of agents’ mobility in the hierarchical geographical network obeys the power-law-like form . Fig. 2 (a) shows that when the value of w increases from 1.0 to 2.0, the power-law exponent, , 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)).

Fig. 2

The move-length distributions of agents in the hierarchical geographical network, each of which lasted for time steps (a), and the average number of static agents in each network node as a function of the layer rank (b), in both plots for the ratio of layer-dependent weights , and , respectively, (c) shows the interval time distribution between the two consecutive contacts established among pairs of agents for different w. All results are averaged over 100 independent runs.

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 and slight fluctuations (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 and the rate of new infected cases , where 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. Fig. 3(c) depicts the relationship between the peak timing of epidemics (i.e. the timing of the maximum ) and w, in which the visible bumps are induced by the fluctuations in the move length distribution .

Fig. 3

Spreading patterns for different values of w in the absence of an intervention. The top panels show the growth of the total infection rate (a), the rate of new infections as a function of time (b), and the peak timing of epidemics vs. the ratio of layer-dependent weight w (c). Panels (d) and (e) display the rate of new infections vs. time for different layers L and the ratio of layer-dependent weight and , respectively. (f) and (g) show the spreading patterns () for different infection source layers and w, with (f) for and (g) for . All displayed results are averaged over 100 independent simulation runs.

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 fast-growing 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 -th layer city at , 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 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, 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 ). When the value of , 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).

Fig. 4

Contagion spread patterns for the hierarchical network (black curves) for different w, and the corresponding epidemic scenario generated with the Lévy-flight model on lattices (red dashed curves), both in the absence of an intervention. (a) and (c) show the growth of total infection rate , and (b) and (d) are the rate of new infections . All results are the averages over 100 independent simulation runs. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Population size and population density effects

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 . 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)).

Fig. 5

The rate of new infections for different L and fixed M (a), and different M and fixed L (b). Panel (c) is the dependence of the peak time on different L and M (without border control); the inset shows the peak time as a function of the size of the network S for different M. Here and other parameters are specified as in the default settings. All results are the averages over 100 independent simulation runs.

Similar effects of population density were observed in our model (Fig. 6 ), where the average population density . 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 shifted-power-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.

Fig. 6

Panel (a) shows the dependence of on and different L when ; the inset in (a) shows the dependence between the peak time and the size of the network S for different ; The dashed lines in panel (a) are the fitting shifted-power functions: (black), (red), and (blue). The dashed lines in the inset represent the fitting logarithmic functions: (black), (red), and (blue). Panel (b) shows the dependence of on different w and when . The dashed lines in panel (b) are the fitting shifted-power functions: (black), (red), and (blue). All other parameters are set as in the default settings and all displayed results are the averages over 100 independent simulation runs. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

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.

Entry screening of long-range travelers

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 real-world 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 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 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 . Obviously, the cost C correlates negatively with : When is large, only a few higher-layer cities will be screened. On the other hand, small 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 for different values of w (in the inset). The bumps visible in Fig. 7(a) are due to the bumpy move length distribution (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 ) will subsequently be quarantined with probability p, or it will keep its infectivity with probability ; thus, p can serve here as the measure of efficiency of border control.

Fig. 7

Disease spreading patterns in the hierarchical network () under long-range travel control, (a) shows the cost of border screening C vs. the ratio of layer-dependent weights w for the critical move length and 27, and C vs. for different values of w (the inset), (b) shows the ratio between the peak time of for the situations with () and without () entry screening of long-range travelers vs. w for different when the efficiency of border control , and the inset of (b) shows the ratio vs. the cost C for different w in semi-log plot, where the total number of agents . (c) shows the ratio as a function of p for and and and 30 in a semi-log plot, where the dashed pink curve shows the fitting function , and the inset of (c) depicts in a log–log plot the ratio vs. the average population density for different values of w when , where the dashed line shows the fitting function + 0.99. All results are the averages over independent simulation runs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We now focus our discussion on the peak timing ratio between the situation with () and without () border control. The higher the ratio 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 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 (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 and the cost of border screening C is logarithmic for (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 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.

Discussion

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 and population density , and a logarithmic positive relationship between 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 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.