Natural human mobility patterns and spatial spread of infectious diseases

We investigate a model for spatial epidemics explicitly taking into account bi-directional movements between base and destination locations on individual mobility networks. We provide a systematic analysis of generic dynamical features of the model on regular and complex metapopulation network topologies and show that significant dynamical differences exist to ordinary reaction-diffusion and effective force of infection models. On a lattice we calculate an expression for the velocity of the propagating epidemic front and find that in contrast to the diffusive systems, our model predicts a saturation of the velocity with increasing traveling rate. Furthermore, we show that a fully stochastic system exhibits a novel threshold for attack ratio of an outbreak absent in diffusion and force of infection models. These insights not only capture natural features of human mobility relevant for the geographical epidemic spread, they may serve as a starting point for modeling important dynamical processes in human and animal epidemiology, population ecology, biology and evolution.

The geographic spread of emergent infectious diseases, epitomized by the 2009 H1N1 outbreak and subsequent pandemic [1], the worldwide spread of SARS in 2003 [2,3], and recurrent outbreaks of influenza epidemics [4][5][6], is determined by a combination of disease relevant human interactions and mobility across multiple spatial scales [7]. While infectious contacts yield local outbreaks and proliferation of a disease in single populations, multiscale human mobility is responsible for spatial propagation [8]. Therefore, a prominent lineage of mathematical models has evolved that is based on reaction-diffusion dynamics [9,10] in which the combination of local exponential growth and diffusive dispersal captures qualitative aspects of observed dynamics. Typically these systems exhibit constant velocity epidemic wave fronts. A related class of phenomenological models is based on the concept of an effective force of infection across distance and thus does not require explicit modeling of dispersal [11,12].
Although more sophisticated models [2,4,13,14] have been developed to describe the dynamics of recent emergent epidemics such as the H1N1 pandemic or SARS, taking into account long distance travel and multiscale mobility networks [15], the majority of these models are still based on the interplay of local reaction kinetics and diffusion processes on networks of metapopulations. The key assumptions of diffusive transport [ Fig. 1(a)] are that (a) individuals behave identically, (b) movements are stochastic, (c) spatial increments are local and as consequence individuals eventually visit every location in the system. Although it is intuitively clear that these assumptions are idealizations and in conflict with everyday experience, the difficulty is how to refine them when data on mobility is lacking, insufficient, or incomplete. Fortunately, a series of recent studies [16][17][18][19] substantially advanced our knowledge on multiscale human mobility. An important discovery that emerged from these studies are individual mobility networks, i.e., individuals typically only visit a limited number of places frequently, predominantly performing commutes between home and work locations and possibly a few other locations. Consequently, individuals or groups of individuals exhibit spatially constrained movement patterns despite their potentially high mobility rate. It has remained elusive how this novel and important empirical insight on individual mobility networks can be reconciled with epidemiological models, to what extent it may impact spatial disease dynamics, and how it may alter spreading scenarios and predictions promoted by ordinary reactiondiffusion models, in which mobile hosts can reach every location in the system.
In this paper we demonstrate how individual mobility networks can be incorporated into a class of models for spatial disease dynamics, and address the questions of whether and how significantly the existence of spatially constrained individual mobility networks impacts on key features of disease dynamics. To this end we investigate a model that explicitly accounts for the bidirectional mobility of individuals between their unique base location (e.g. their home) and a small set of other locations [ Fig. 1(b)]. The entire population is therefore represented by a set of overlapping individual mobility networks associated with each base location, see, e.g., [20][21][22]. We focus on the analysis of epidemics on regular lattices and complex metapopulation networks and systematically compare the dynamics to reaction-diffusion systems as well as the force of infection models. In conflict with reactiondiffusion models, in which front velocities increase with travel rates unboundedly, we show that our model predicts a saturation of wave front velocities, a direct consequence of the rank of locations in individual mobility patterns. This suggests that estimates for propagation speeds may have been considerably overestimated in the past. Furthermore, we analyze a fully stochastic model to show that both, in regular lattices as well as complex metapopulation networks, the global outbreak of a disease is determined by a novel type of threshold for the attack ratio of the disease which depends on the characteristic time spent at distant locations. Finally, we show that in the limit of low and high travel rates our model agrees with reaction-diffusion models and the direct force of the infection class of models, respectively.
Consider a system of M populations labeled m and assume that in each an epidemic outbreak can be described by a compartmental SIR model, i.e., where S m , I m , and R m label and quantify susceptible, infected, and recovered individuals of population m. Infections and recovery events occur at rates and , respectively, with > . The number of individuals in population m is given by N m ¼ S m þ I m þ R m , which, however, is conserved only in the statistical sense at equilibrium. This system of reactions yields the mean-field descriptions @ t I ¼ IS=N À I and @ t S ¼ ÀIS=N. A natural and plausible extension to a system of M coupled populations is diffusive dispersal among those populations defined by hopping rates w nm > 0 from population m to n, which yields a metapopulation reaction-diffusion system, i.e., for infected and susceptible individuals, @ t I n ¼ S n I n =N s n À I n þ X m ðw nm I m À w mn I n Þ; where N s n is the number of individuals in population n in diffusive equilibrium. Travel rates ! mn are usually estimated by gravitylike laws [23]. N s n is determined by detailed balance, i.e., N s n =N s m ¼ w nm =w mn and the conservation of the number of individuals in the metapopulation N ¼ P m N s m . These types of models have been employed in numerous recent studies [2,14,24,25]. The tight connection to spatially continuous reaction-diffusion systems is revealed for the special case of a linear grid of populations placed at regular intervals l at positions x m ¼ ml and mobility between neighboring populations, i.e., w nm ¼ ! nÀ1;m þ ! nþ1;m which yields (3) is related to a classic form of a Fisher equation which has been deployed in mathematical epidemiology [9,10]. For sufficiently localized initial conditions this systems exhibits traveling waves with speed Note that this velocity increases monotonically with the mobility rate !.
In order to account for individual mobility patterns we propose the following generalization: individuals possess two indices n and k, characterizing their current and their base location, respectively. The dispersal dynamics is governed by a Markov process where X stands for each infectious state S, I, and R. Equation (5) implies that individuals with base location k possess a unique mobility rate w k nm that determines how they travel between locations n and m. This yields a generalization of (2) where I k n and S k n are the number of infected and susceptible individuals of type k that are currently located in population n. N n is the total number of individuals at n, i.e., N n ¼ P k ðI k n þ S k n þ R k n Þ. If ! k nm are k independent, we recover the reaction-diffusion case described above. In the following we consider the case of overlapping star-shaped networks corresponding to commuting between base and destination locations only. This implies that ! k nm ¼ 0 if either k Þ n or k Þ m. We further assume the constant return rate w k mk ¼ w À for all k and m. Realistically we have ! n kn =! À ( 1, implying individuals belonging to n remain at their base most of the time. If mobility rates are large compared to the rates associated with the infection and recovery, i.e., ! k mk , ! À ) , , the detailed balance condition is fulfilled for infected and susceptible individuals separately and the last term in Eq. (6) vanishes in which case the above model is equivalent to the remote-force-of-infection model-see also [11,21].
In general, it is difficult, if not impossible, to extract the dynamic differences between the systems defined by Eqs. (2) and (6) for an arbitrary metapopulation. Yet, it is crucial to understand the key dynamic differences that emerge when individual mobility patterns replace ordinary diffusion. We therefore investigate the impact of individual mobility patterns in the instructive one-dimensional lattice case. We assume that only infected individuals can travel and that recovery is absent (SI model), i.e., ¼ 0 (relaxing these two restrictions does not change the main results but clarifies the analysis). We denote the number of infected individuals with base location n that are located on site n (i.e. their base) and sites n AE 1 by I n n and I AE n , respectively. At rates ! þ and ! À individuals leave and return to their base. In order to obtain a spatially continuous description we approximate I AE nAE1 ! I AE ðx AE lÞ % I AE AE lrI AE þ l 2 2 ÁI AE and introduce relative concentrations u n ¼ I n =N, v n ¼ ðI þ n þ I À n Þ=N leading to where D ¼ l 2 =2. The function uðx; tÞ is the fraction of individuals based on and located at x whereas vðx; tÞ is the fraction of individuals based on but not located at x. Note that we discarded here a third equation for w ¼ ðI þ À I À Þ=N independent from (7) with a solution vanishing at long times. This does not change the main result (8). The nontrivial steady state is given by One of the key questions is if the system (7) exhibits stable propagating waves as solutions and how their velocity depends on system parameters. The ansatz uðx; tÞ 1 uðx À ctÞ, and likewise for v, leads to a stable solution with a velocity given by Fixing ! À , this implies that c ! 0 as ! þ ! 0. On the other hand, if ! À is small, the system is entirely determined by the forward rate ! þ . For a systematic comparison to the ordinary reaction-diffusion system, we have to establish a relation between ! AE and ! in Eq. (4). For the sake of simplicity we consider the symmetric case, i.e., ! þ ¼ Numerical simulations of a fully stochastic system nicely agree with the analytical predictions (Fig. 2). The essential feature of cð!Þ is its saturation as ! ! 1 whereas c $ ffiffiffiffi ! p for reaction-diffusion systems. This effect is a direct consequence of the spatial constraints imposed by individual mobility patterns, i.e., increasing the mobility rate ! yields a higher frequency of travel events but is restricted to the set of locations connected to the base location and thus is universal and also holds for more complex metapopulation topologies [26]. For high rates, the deviation between ordinary reaction-diffusion and bidirectional mobility increases without bounds. This is an important insight as the expression for wave front speeds for reaction-diffusion have been used in the past to estimate wave front speeds from mobility rates and vice versa. Our results indicate that if constrained mobility patterns and bidirectional movement patterns are at work these estimates must be treated with particular care. Note that in the limit ! ! 1 Eq. (7) reduces to the heuristic equation in Ref. [27].
A key question in epidemiological contexts concerns conditions under which an epidemic outbreak propagates or wanes. Outbreak dynamics are usually triggered by crossing thresholds inherent in the system's dynamics [28]. For example the basic reproduction number R 0 of a disease, i.e., the expected number of secondary cases caused by a single infected individual in a susceptible population represents such a threshold [8]. In the single population SIR dynamics introduced above R 0 ¼ =. If R 0 > 1, an outbreak occurs, otherwise the epidemic dies out. In the extended metapopulation system with diffusive host movements a second threshold, known as the global invasion threshold [24], is induced by the flux of individuals between populations, which must be sufficiently high for an epidemic to spread throughout the metapopulation.
An interesting and important feature of the model we propose is the existence of another invasion threshold which is determined by the return rate ! À or equivalently by the time an individual spends outside the base location. This finding is illustrated in Fig. 3, which depicts the attack ratio as a function of ! À for a bidirectional SIR epidemic on a lattice. The attack ratio is defined as the fraction of the overall population affected during an epidemic. We fixed the entire interlocation flux at a value above the global invasion threshold ensuring a global outbreak in the reaction-diffusion system. This result is universal for the bidirectional system and is independent of the particular topology. Insets (a) and (b) of Fig. 3 display simulation results for uncorrelated scale-free and Erdős-Rényi networks. For low return rates the attack ratio is close to unity. However, with growing values of the return rate, the attack ratio decreases steeply and vanishes, reflecting the absence of the global outbreak. The regime of high return rates corresponds to small dwelling times on distant locations. Infected individuals do not have enough time to transfer the disease to susceptibles in unaffected locations before they return. Note that for a wide range of local population sizes (50-2000) the attack ratio collapses onto a universal curve $ ð! À =NÞ which suggests that the basic mechanism of the invasion threshold and its functional dependence on return rate ! À can be understood theoretically.
To estimate the observed invasion threshold on a onedimensional lattice for R 0 * 1 we calculated the number of infected individuals that originate from an affected location and can seed an outbreak in an unaffected one [28]. These seeders are approximately given by the total number of individuals entering the new location per unit time, i.e., N! multiplied by the typical time they remain active in the destination location. This time is given by the inverse rate r of exiting the seeders class. Since two processes (returning to the base location and recovery) can trigger this exit, r ¼ þ ! À . Therefore, % N!=ð þ ! À Þ. In the tree approximation the threshold condition for an epidemic on a network with an average degree " k is given by ð " k À 1ÞðR 0 À 1Þ > 1 [24]. Thus, for a one-dimensional lattice the threshold condition reads where we kept the total interlocation flux ! ¼ 2! þ ! À = ð2! þ þ ! À Þ constant (this relation and positivity of travel rates impose the restrictions: ! À > ! and ! þ < !=2) and used ! þ ( ! À . Moreover, as extensive simulations show, the global invasion threshold in terms of the flux rate ! exists in bidirectional systems only for low return rates ! À . With increasing return rate the disease fails to propagate globally, even for constant interlocation flux !. This effect, observed in lattice topologies and paradigmatic more realistic network topologies, is a fundamental consequence of birectional mobility patterns. This effect is absent in reaction-diffusion systems fully determined by the total interlocation flux. This property implies that the invasion of a front is limited only by return rates and not as is typically assumed by the overall particle flux !. In the context of disease dynamics this implies that more efficient containment strategies could be potentially devised that do not target the overall mobility but rather, modify mobility patterns in an asymmetric way. Note that the transition takes place only in a system with a finite number of agents per site; in the system with an infinite number of individuals this effect would disappear, as also confirmed in Fig. 3. Indeed from Eq. (10) for ( ! À the scaling $ ð N! ! À Þ follows-with an increasing number of individuals per site N, the threshold value of the return rate ! À increases. Recently, an unprecedented amount of detailed information on human activity became available, requiring the revision of established models and the development of new more sophisticated ones. We investigate an epidemiological model explicitly incorporating such important properties of human mobility as individual mobility networks and frequent bidirectional movements between home and Results of stochastic simulations. Attack ratio ð! À =NÞ as a function of the return travel rate ! À for SIR epidemics on a lattice with 10 2 sites with N agents/site. Epidemic parameters: ¼ 1, ¼ 0:1. Insets: (a) ð! À =NÞ for a scale-free network ( ¼ 1:5) of 10 3 nodes populated uniformly with hNi ¼ 250=site, with k min ¼ 5 and k max ¼ 50; (b) ð! À =NÞ for an Erdős-Rényi network with 500 nodes and " k ¼ 10. Results were averaged over 50 realizations. Interlocation flux rate was h!i ¼ 1.
destination locations. It manifests surprising dynamical features such as the existence of bounded front velocity and a novel propagation threshold, which are universal for any metapopulation topology. As more detailed data continue to become available our approach is a promising foundation for further research.