How the individual human mobility spatio-temporally shapes the disease transmission dynamics

Human mobility plays a crucial role in the temporal and spatial spreading of infectious diseases. During the past few decades, researchers have been extensively investigating how human mobility affects the propagation of diseases. However, the mechanism of human mobility shaping the spread of epidemics is still elusive. Here we examined the impact of human mobility on the infectious disease spread by developing the individual-based SEIR model that incorporates a model of human mobility. We considered the spread of human influenza in two contrasting countries, namely, Belgium and Martinique, as case studies, to assess the specific roles of human mobility on infection propagation. We found that our model can provide a geo-temporal spreading pattern of the epidemics that cannot be captured by a traditional homogenous epidemic model. The disease has a tendency to jump to high populated urban areas before spreading to more rural areas and then subsequently spread to all neighboring locations. This heterogeneous spread of the infection can be captured by the time of the first arrival of the infection (Tfi)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(T_{fi} )$$\end{document}, which relates to the landscape of the human mobility characterized by the relative attractiveness. These findings can provide insights to better understand and forecast the disease spreading.

Over the past several decades, outbreaks of emerging infectious diseases have been occurring at an increasing rate 1 . They have a considerable impact not only on public health and health care but also on various forms of social and economic well-being 2 . However, despite the increased attention on infectious diseases, there is still a lack of comprehensive knowledge of mechanisms controlling the spatiotemporal propagation of infectious diseases. Traditional homogenous epidemic models, e.g., a compartmental epidemic model, have long been helpful in understanding the foundation of transmission dynamics of infectious diseases. These homogenous epidemic models assume all human individuals to be well-mixed and ignore heterogeneity that may be resulted from, for example, ages, locations, and contact patterns of human individuals. Although these homogeneous models are useful in certain situations, in order to gain more accurate modelling results, incorporating factors describing spatiotemporal heterogeneity in disease transmission dynamics into the epidemic models may be necessary.
In reality, the transmission dynamics of real-world infectious diseases are governed by several spatiotemporal heterogeneous factors, including the seasonality of diseases 3 , the contact networks of the host population 4 , population heterogeneity 5 , and human mobility 6 . Human mobility, in particular, undeniably plays a crucial role in the temporal and spatial transmission dynamics of infectious diseases. During the past few decades, the inclusion of the human mobility mechanism in the epidemic models has been considered in many studies for modeling of the geographical spread of diseases [7][8][9][10][11][12][13][14] . Human traveling not only has a significant effect on the epidemic spreading [15][16][17] , it also remarkably influences the speed of disease spreading 7,8,15 . Moreover, understanding the roles of human mobility on disease spreading is notably useful for designing effective disease control strategies, such as vaccination or travel restriction [18][19][20][21][22] .
Due to the approachability of human travelling data (e.g., check-ins from locations, GPS navigators, or mobile-phone records), scientists have discovered several properties and characteristics of human mobility [23][24][25][26][27][28][29] . A myriad of human mobility models has been developed to understand the fundamental mechanism hidden behind these findings 25,[30][31][32] . These human mobility models are becoming more and more accurate, capable Scientific RepoRtS | (2020) 10:11325 | https://doi.org/10.1038/s41598-020-68230-9 www.nature.com/scientificreports/ of reproducing human mobility either at an individual level 25 or at collective information 31 . Although several research works have studied the effects of human mobility on the spread of diseases [33][34][35][36] , incorporating a realistic individual human mobility model into an epidemic model is still needed. Integrating an epidemic model with an individual human mobility model is necessary for investigating how an individual infectious person geotemporally spread the disease.
In this article, we address the question of how individual human mobility spatiotemporally shapes the spread of communicable infections. To this end, we integrate a classical SEIR (susceptible, exposed, infectious and recovered) individual-based model with an individual human mobility model proposed by Yan et al. 37 , both at the individual and population level, where the transition probabilities between localities are governed by both a gravity-like part and a memory part. To follow the spread of the infection, we have introduced an indicator, "the time of the first arrival of the infection ( T fi ) ", so that an infection which started at an origin point reaches another given point for the first time. To characterize the landscape of human mobility, we also defined the "relative attractiveness (RA)" based on a transition matrix between localities. Now, our objective of investigating how individual human mobility shapes the spread amounts to determining or constructing the relationship between T fi and the explanatory variable, RA. Indeed, the quantity T fi can not only elucidate the direction of epidemics driven by human mobility but also provide the velocity of the spreading by calculating its spatial derivative in the RA space. As an illustrative application, we considered the spread of human influenza in two contrasting countries, namely, Belgium and Martinique, as case studies, to assess the specific roles of human mobility on infection propagation. Interestingly, we found in both cases that T fi follows a truncated power law as a function of RA, indicating a one-to-one relationship between spatiotemporal infection spread and host mobility landscape.

Results
Individual human mobility-driven disease transmission dynamics. A geo-temporal spreading pattern of an epidemic in Belgium and Martinique is illustrated in Fig. 1. At the beginning of the simulation, there is only one symptomatic infectious individual at the most populated location (red dot in Fig. 1), while all other individuals are susceptible to the disease. Each individual stays at a particular location for a certain period of time, ( T ), following the waiting time distribution with the average waiting time � T� = 5 days. We found that at the early time of the epidemic, the disease is likely to reach urban areas with high population densities and after an intermediate time (e.g., at day 50), it spreads to rural areas with lower population densities. Lastly, at the late time, the disease is gradually transmitted to all neighboring locations and reaches the boundary of the country. However, the epidemic in Belgium seems more dispersed than in Martinique, which might be due to the fact that the population distribution in Belgium is more heterogeneous than in Martinique.
To have an idea of how fast the disease spatially spreads, we computed the radial-averaged time of the first arrival ( T r ) of infection in each location (Fig. 2a) and the radial speed of spread at radius r centered at the first infected location (Fig. 2b). At the early time of the epidemic, the disease spreads locally and slowly, spending time around 50 days to reach 15 km in Belgium (Blue line) and around 40 days to reach 10 km in Martinique (red line). In contrast, after some periods of time, the disease spread faster and faster due to the jumping of infectious individuals to distant locations (Fig. 2b). Note, however, that, although the disease spatially spreads faster in Belgium, the epidemic curve showing the percentage of infectious individuals in Belgium increases slower than that in Martinique (Fig. 2c). This might be because Belgium is bigger, and the population distribution in Belgium is more heterogeneous than that in Martinique. Fluctuations in T r indicate that different locations at similar r have different mobility, causing the variety of the arrival time of infectious individuals (Fig. 2a). Moreover, when considering the effect of the averaged waiting time on the radial speed of disease spreading, we found that increasing the averaged waiting time affects the speed of the spatial spreading of disease in Belgium more than in Martinique (see Figure S4).
In addition, the proposed model allows us to explore the movement behaviors of infectious individuals by tracking their movement trajectories. We elucidated the disease dispersal by computing the maximum distance of disease transmission, the time-averaged mean square displacement (MSD) of infectious individuals, and the number of infected locations (Fig. 3) for both Belgium and Martinique. However, since tracking the movement trajectories of the vast population requires a large amount of computer memory; therefore, the location and movement tracking of infectious individuals in Belgium was performed only in the case of the population size is reduced by 100 times. We have checked that after scaling down the population density in Belgium, there was no location with a zero population. The effects of scaling down the population density are discussed in the supplementary information. Also, to investigate the effects of waiting time on the dynamics of disease transmission, the model was simulated using a different value of waiting time ( T = 1, 3, 5, 10, and 30 days). The maximum distance of disease transmission was measured from the radius between the first infected location and the farthest infected location at a certain time (Fig. 3a,d). We found that the disease initially spreads locally around the first infected location and then spreads rapidly and widely to other locations until it reaches the locations at the boundary with the maximum distance of 243 km and 42 km for Belgium and Martinique, respectively. Figure 3b, e show the MSD, δ 2 (�) , of infectious individuals as a function of lag time (Δ). Note that the MSD of infectious individuals could be measured only during the first 20 days after they become infectious. This is because the probability that an infectious individual will get recover after time Δ since they become infectious decays as 1 − exp (−γ �) , where γ = 1/(3 days); therefore, less than 0.2% of infectious individuals are infectious longer than 20 days. Finally, according to Fig. 3c, f, the increasing number of infected locations is consistent with the increase of the maximum distance of disease transmission shown in Fig. 3a, d. However, due to a higher heterogeneity in the population distribution in Belgium, in the case of the averaged waiting time of 30 days, at day 150, only 79% of locations in Belgium are infected while almost 100% of locations in Martinique are infected (see also the epidemic curves in Figures S11a, S12a). www.nature.com/scientificreports/  www.nature.com/scientificreports/ The relative attractiveness (RA) was introduced to characterize the human mobility landscape. The RA maps for Belgium and Martinique are shown in Fig. 4a, b, respectively. We found that the RA is proportional to the population size in each location ( Figure S7a, b), and the distribution of RA in Martinique seems more homogeneous than in Belgium ( Figure S7c). The RA was classified into four groups with different colors; 0 < RA < 0.82, 0.82 < RA < 1.22, 1.22 < RA < 2, and RA > 2 representing low, medium, high, and very high relative attractiveness of locations, respectively. As can be seen in the RA maps, the locations with very high RA are more scattered in Belgium than in Martinique, and the locations with low RA are in a greater proportion in Belgium than in Martinique. Interestingly, for both Belgium and Martinique, there are approximately six locations in which the RA is comparatively very high ( Figure S7d). Moreover, we found that the average RA of those six locations is three times greater in Belgium than in Martinique. Therefore, in Belgium, individuals have a tendency to be trapped in only a few particular locations in contrast to Martinique, where individuals are likely to visit more various locations.  www.nature.com/scientificreports/ To examine how human mobility shapes the disease spreading, RA is compared with the contour of the first arrival time ( T fi ) of disease occurrence in each location, as shown in Fig. 4c for Belgium and Fig. 4d for www.nature.com/scientificreports/ Martinique. The contours show that the disease frequently jumps to an area with high RA and then spread into a lower RA area. Interestingly, the relationship between RA and T fi (Fig. 4e) follows a truncated power law as, where t 1 represents the arrival time that disease reaches on "flat lands" defined as locations where RA = 1, ν and RA cut−off describes the distribution of T fi for each location. Due to a more homogeneity of RA in Martinique, the disease is likely to spatially spread faster than in Belgium (Fig. 2b). Flat lands, especially, are reached by the disease around 39 days in Martinique and around 55 days in Belgium after the first infectious individual was introduced (see the fitted equations in Fig. 4e).
In addition, we also investigated the effect of waiting time and the population density on the relation between RA and T fi . As, in our model, the human mobility and the epidemiological dynamics are deeply integrated, changing the population density will simultaneously affect both of them. We found that both factors significantly change the relationship in Belgium (Figures S8, S11) but not in Martinique (Figures S9, S12). All fit parameter values are shown in Figures S10 and S13. Besides RA, to understand how infectious individuals move, we measured the infectious-individual visitation frequency of a location by counting the number of times that infectious individuals visit a location. We found that the pattern of the visitation frequency map ( Figure S5) is similar to the pattern of the RA map. Moreover, we calculated the probability distribution of the number of locations visited by infectious individuals, P(S), and the probability distribution of the number of step lengths among the visited locations of infectious individuals shown in Figure S6.

Effects of restriction on symptomatic infectious individual movement.
Our model allows us to investigate disease control strategies at the individual human level. In the present work, the model was employed to investigate the effectiveness of restriction of infectious individual movements in Belgium. To reduce the time required for simulating the model, the simulations presented in this section were performed only in the case of the population size of Belgium is reduced by 100 times. The effects of scaling down the population size are discussed in the supplementary information. In practice, only symptomatic infectious individuals may be able to be identified, hence in our model, only (some) symptomatic infectious individuals are banned from traveling between locations while asymptomatic infectious individuals are free to travel according to the human mobility mechanism. Note that although the travel-restricted symptomatic infectious individuals cannot transmit the infection to other locations, they can transmit the disease within the locations that they are staying while they are infectious. The time course of the numbers of infectious individuals at a different percentage of symptomatic infectious individuals restricted to travel to other locations, η , are shown in Fig. 5a. The numbers of infectious individuals decrease significantly when the restriction percentage increases. For 100% travel restriction (green curve), it could delay the epidemic peak from day 89 to day 101 and could also reduce the epidemic size from 65.22 to 57.85%. Figure 5b shows a change in the relation between RA and T fi , in which all fit parameter values are shown in Figure S14. The disease spreads slower when a higher number of symptomatic infectious individuals is banned from travelling; the arrival time of disease reaching on flat lands, t 1 , is 78, 86, 88, 90, and 97 days for η = 0% , 20% , 50% , 70% , and 100% , respectively. All contours of T fi with different η are shown in Fig. 5c-g.

Discussion and conclusions
Of several factors, human mobility is an important factor affecting the spread of epidemics. Many researchers have been extensively investigating the role of human mobility on the spatial spread of influenza-like illnesses (ILIs) and other emerging infectious diseases 6,9,10,13,18,[38][39][40] . However, how the infection or the transmission chain would spread and spread with what associated speeds due to human mobility remains to be a challenging question. Therefore, to deeply understand the impact of human mobility on the epidemics, we exploited a classical SEIR individual-based model that integrates with a model of individual human mobility as proposed by Yan et al. 37 where human mobility is described, both at the individual and population level, with transition probabilities between localities governed by both a gravity-like part and a memory part. Yan et al. 37 have shown that, in the long-time regime, the mobility patterns produced by their model are stable and robust, their model is therefore capable of simulating the disease transmission dynamics in a long duration of time.
Our approach can provide a geo-temporal spreading pattern of the epidemics that cannot be captured by a traditional homogenous epidemic model. The model predicts that the disease has a tendency to jump from (high populated) urban to urban areas before spreading to more rural areas and then subsequently spread to all neighboring locations (Fig. 1); a finding supported by other modelling studies 12,41 . In addition, this heterogeneity level of the epidemic spreading can be captured well by the time of the first arrival of the infection ( T fi ) (Fig. 4c, d) which is the time when the first infection is found in each location after the beginning of epidemics introduced at an origin point (the most populated location). This indicator was also used in previous work 33,42 for comparing the epidemic outcomes between countries. Since the distribution of T fi in each location is heterogeneous, the radial-averaged of time of the first arrival of the infection ( T r ) is used to represent how disease propagates through the landscape (Fig. 2). At the early phase of epidemics, the disease spreads slowly to reach neighbor locations due to a small number of infectious individuals, but at a certain time, it spreads rapidly and widely within few days to reach distant locations.
In order to gain a better understanding of how human mobility shapes the disease spreading, we used the "relative attractiveness (RA)" based on only the transition matrix between localities to characterize the landscape of human mobility. We found that the RA map highlights highly attractive areas versus the others, the locations with a higher (lower) RA are more (less) attractive and thus having a higher (lower) tendency to be visited by infectious individuals (Fig. 4a, b). In addition, the T fi map can provide information on how disease spreads in both temporal and spatial scales. By integrating the information from the RA and T fi maps, our analysis indicates that the disease is more likely to reach the highly attractive areas faster than the low attractive areas. Interestingly, the relation between RA and T fi follows the truncated power equation. The fitting parameter t 1 can also be Scientific RepoRtS | (2020) 10:11325 | https://doi.org/10.1038/s41598-020-68230-9 www.nature.com/scientificreports/ www.nature.com/scientificreports/ used to illustrate the time of the first arrival of the infection spreading into 'flat lands' , defined as locations with RA = 1. At the individual level, we tracked the movement trajectories of infectious individuals and recorded the visitation frequency of infectious individuals in each location ( Figure S5). We found that the number of times that infectious individuals visit a location is proportional to the value of RA of that location. The locations with higher RA are more likely to be visited by infectious individuals than those with a lower RA. Our model also allows us to investigate intervention strategies at the individual level. We showed that the travel restriction of symptomatic infectious individuals provides modest results; it can decrease the epidemic size and delay the epidemic peak date (Fig. 5). The simulated incidence profile clearly shows a slowdown in the growth of the number of new infections when the travel restriction percentage increases. Our model predicts a 1-2 weeks delay of the epidemic peak when the travel restriction measure is implemented. These results are similar to those of a previous study, which examined the effect of both border and internal travel restrictions on an influenza pandemic in the United States 22 . In addition, we found that the individual waiting time also affects the epidemic dynamics ( Figures S11-S13). A longer averaged waiting time of human mobility can delay the time of the first arrival of the infection in each location, extend the epidemic peak, and reduce the epidemic size.
The epidemic model presented here also has some limitations. Some parameters in the model, such as the waiting time exponent and the cutoff time in the waiting time distribution, and the memory parameter, are very difficult to measure and therefore need to be estimated. The values of these parameters depend on human mobility behavior and hence can vary in different countries 43 . These may make it difficult to apply the proposed model in a different setting. However, in order to investigate how the individual human mobility spatio-temporally shapes the disease transmission dynamics in different geographical structures, in this study, these human mobility parameters were kept the same when simulating the disease transmission in both Belgium and Martinique. Besides, since our model keeps track of each individual in the population, it can be time-consuming and expensive when simulating the model in a country with a very large population. However, we have performed a sensitivity analysis for the change in the population density ( Figures S8-S10), we found that when the population size increases, the epidemic peaks are relatively shifted to the right. This might be due to the fact that when the population size increases, there will be more susceptible individuals available in the system; it, therefore, takes longer for the epidemic to reach its peak. However, the percentages of the number of infected individuals are not much different when the population size increases, especially in Martinique. Hence, these findings could help us investigate the roles of human mobility at the country-level solely by scaling down the population size in the simulations. For the sake of simplicity, we also neglected the human migration between countries. Furthermore, the model does not consider the heterogeneity in the population, such as ages, households, and social contacts; as well as specific locations of, for example, schools, homes, or workplaces. All of these factors may affect disease transmission 4,6,[44][45][46] . Finally, we did not consider a change in human behavior under the threat of the epidemic; for instance, in reality, during an epidemic, people may avoid travelling to crowded locations.
Due to the explicit representation of human individuals, our model can allow us to understand the mechanism hidden in human mobility that affects the epidemics both at population and individual levels. Future work could employ the model to investigate the characteristics of the areas or the behaviors of an infectious individual which accelerates disease transmission or a super-spreader. Moreover, it can be used to explore other intervention strategies for containing or mitigating epidemics such as quarantine of infectious individuals and targeted vaccination. We believe that the proposed model would be beneficial for the modeling of disease spreading and imperative for the preparedness plans.

Methods
Population data. The spread of infection was studied in two contrasted areas, Belgium and Martinique, to investigate different human mobility landscapes. Belgium, a European country, covers an area of 30,688 km 2 with a population of 11,250,658 inhabitants, and Martinique is a French territory of 1,128 km 2 with 408,596 inhabitants. The geographical areas of Belgium and Martinique were divided into 2,252 and 86 grids with the 2.5-arcminute resolution (5 km at the equator), respectively. The grid structures and population sizes within each grid cell (locality) for both Belgium and Martinique were obtained from the socioeconomic data and applications center (SEDAC) 43 .

The SEIR epidemic model with human mobility.
To investigate the roles of human mobility on the spreading of epidemics, we consider the spread of human influenza over an area comprising of N localities or locations, each of which with a population m i such that the total population of the area is, M = N i=1 m i . Two processes are considered: (1) individual human mobility between distinct localities using the model adapted from the model proposed by Yan et al. 37 (Fig. 6a) and (2) within each locality, the transmission of infection between individuals is described using a compartmental SEIR model (Fig. 6b).
Human mobility. The locations in the study area are represented in Fig. 6a by circles of size proportional to their population size m i . In this model, the probability that an individual initially at location i at time t travels to location j at time t + T is given by G ji (t, t + �T) = P i (�T) × p ij , where P i (�T) is the distribution of waiting times at the initial location before moving and p ij is the transition probability. The distribution of waiting times is given by 23 where α i is the waiting time exponent and τ i is the cutoff time of the location i . For simplicity, we assumed in all simulations that P i (�T) is independent of the location and, therefore, identical for all locations www.nature.com/scientificreports/ ( P i (�T) = P(�T) ). Details of the method for estimating the waiting time parameters, α and τ , are provided in the supplement information. The transition probability for traveling from location i to location j is governed by two elements. The first element is a gravity-like contribution, where m j is the population of destination location j and W ji is the total population within the circular region centered at location j having the radius equal to the distance between locations i and j as illustrated by the dashed circle in Fig. 6a. In this model, the transition matrix B ij is characterized by peaks of the points more likely to be visited at short and long distances from the origin (jump transitions) than it would be in the case of the classical gravity-like model. This heterogeneity of the transition matrix would be more marked when the distribution of population densities is heterogeneous. The second part is a memory effect where r j is the index indicating the order rth of the visit of location j along the movement trajectory (shown at the bottom of Fig. 6a), with r = 1 for the home location, and λ is a constant parameter representing the strength of the memory effect. Therefore, the transition probability that an individual travels from location i to locations j is, Note that A j = 1 when location j has never been visited (gray circle) before. It follows that the probability for an individual exploring new locations is simply, p ij = m j /W ji . Full description of the model parameters and their meaning can be found in 37 .
SEIR individual-based model. The population in each location (enlarged circle in Fig. 6a) is classified into four epidemiological classes: susceptible (S), exposed (E), infectious (subdivided into asymptomatic (I A ) and symptomatic infectious (I S )), and recovered (R) (Fig. 6b) 4 . Within each location, susceptible individuals become is shown at the bottom. The transition probability for an individual moving from location i to location j is governed by two elements, namely, a gravity-like part ( B ij = m j /W ji ) and a memory effect A j = 1 + /r j . If location j has never been visited before (grey circle), its A j value is assumed to be unity, A j = 1 ; thus, it will be visited with p ij = m j /W ji . For the previously visited location j, it can be visited again with p ij = m j W ji 1 + r j , where r j is an index indicating that location j is the rth visited location. At a certain time, each location contains a number of human individuals who will stay at the location for different periods of time drawn from the waiting time distribution P(�T) 23 . The figure was adapted from the model proposed by Yan et al. 37 . (b) Based on the individual's disease transmission dynamics, individuals are categorized into five compartments, namely, the susceptible ( S ), exposed ( E ), asymptomatic infectious ( I A ), symptomatic infectious ( I s ), and recovered ( R ) compartments. Individuals in the exposed compartment move to the asymptomatic infectious compartment with the probability p a . www.nature.com/scientificreports/ infected, and progress to the exposed class, following contacts with a symptomatic or an asymptomatic infectious individual staying at the same location as described by the force of infection F i given by where m i is the population size at location i, β is the disease transmission rate, and ω is a scaling parameter that takes into account the reduced infectiousness of asymptomatic infectious individuals. It is worth noting that, since influenza viruses spread from person to person primarily through large-particle respiratory droplet transmission, which requires close contact between source and recipient persons. In our study, we, therefore, assumed that the influenza transmission is the frequency-dependent transmission. Such an assumption was also used to predict the spread of influenza in several studies 6,20,38,41,47,48 . After being infected, exposed individuals become either asymptomatic (with a probability p a ) or symptomatic infectious (with a probability 1 − p a ) at a rate ε , which is inversely proportional to the incubation period. Both symptomatic and asymptomatic infectious individuals recover from the disease at a rate, γ, which is inversely proportional to the infectious period, and have lifelong immunity. The disease transmission rate was calculated using the following relationship 41 : The descriptions and values of all parameters used in the model are summarized in Table 1.

Mobility landscape: the relative attractiveness of locations (RA).
To characterize the human mobility landscape, we calculated the relative attractiveness (RA) of each location of the study area. To this end, we consider the Master equation describing the movements of individuals in the study area as, where ⇀ P is a column vector of probabilities of finding individuals at location i and time t , B t is the transpose of the matrix B (of the gravity-like part of the mobility) with elements B ij = m j /W ji and the matrix K is given as, where N is the total number of locations. The matrix B t − K describes the fluxes of individuals traveling forward and backward between all locations. At the stationary state, the time variation of the probability vector ⇀ P is zero, i.e., B t − K ⇀ P = 0 . It thus follows that the stationary probability is obtained as, U is the (normalized) eigenvector of B t − K associated to the zero eigenvalue. Finally, the relative attractiveness RA i of location i can be defined as a ratio of the stationary probability of finding individuals at location i , P st,i , and the equi-probability, 1/N, of finding individuals at any location; Therefore, locations with RA i > 1 have a higher propensity to be visited by individuals than those with RA i < 1.

Time of the first arrival of the infection.
When it comes to an epidemic, we are mainly interested in knowing how the infection or the transmission chain would spread and spread with what associated speeds. A commonly used approach to analyze an epidemic spread in terms of directions and velocities is the "trend surface analysis (TSA)" [29][30][31][32][33] , to name just one. Unfortunately, such a method cannot be used to predict the profile of the spread and it becomes inappropriate for analyzing the spread of an epidemic due to the heterogeneity of host www.nature.com/scientificreports/ movements. To address how epidemics spread through different locations, we used the time of the first arrival of the infection, T fi (i|j) , which is the average elapsed time that an infection started out at location j reaches location i for the first time. In all simulations reported in this work, we assumed that the first infected location j was the most populated one so as to use T fi (i) instead of T fi (i|j) , for simplicity. Due to different pathways of infection between two locations, the distribution of velocity is heterogeneous. To have an idea of how the expansion of epidemics is, we computed the radial-averaged time of the first arrival of the infection T r for all locations at radius r from the starting point of epidemics. Then the radial speed of expansion v(r) was calculated by Mean square displacement. In the simulations, we tracked the movement trajectories of infectious individuals to understand how infectious individuals move after they got the infection. We recorded the sequences of visited locations, step lengths, and times that infectious individuals spent at each location. To analyze the mobility of infectious individuals, we calculated the time-averaged mean square displacement (MSD), δ 2 (�) , of each infectious individual 49 ; where δ 2 i (�) is MSD of infectious individual i at a different lag time (Δ), T i is a tracking period, a time duration since individual i got infection to the time when he/she recovered, Δ is a lag time ( � < T i ), r i (t) is a spatial displacement of infectious individual i at time t after being infectious. We then averaged the MSD of all infectious individuals, δ 2 (�) . This MSD describes how far, on average, infectious individuals travel during their infectious duration.
Simulation details. In the simulation, individuals in each location initially stay in particular locations for certain periods of time drawn from the truncated power law (Eq. 1) and then move to other locations following the transition probability (Eq. 4). All visited locations and their orders of the visit of each individual are recorded for calculating the memory effect.
Before introducing an infection into the model simulations, human mobility was first pre-equilibrated by simulating the human mobility model until the total human mobility flux of each location reaches a constant, i.e., the number of individuals moving in equals to the number of individuals moving out. Unless stated otherwise, we assumed that the first infectious individual is initially introduced in the most densely populated location. At each time step of one hour, the disease transmission dynamics within each location was simulated using the tau-leaping method [50][51][52] . Besides spreading within a location, the disease can be transmitted spatially to other locations via human mobility in which its speed relies on the waiting time distributions and the transition probability (Eqs. 1,4). In this model, we assumed that the dynamics of human birth and death are much slower than the dynamics of the epidemics and human mobility, and there is no emigration and no immigration; thus, the human population size is constant. The descriptions and values of all parameters used in the model are summarized in Table 1. All simulations were performed using MATLAB software (version R2018b, The MathWorks, Inc) running on personal computers. All maps showing the spatio-temporal disease spreading patterns in Belgium and Martinique were also generated with the use of MATLAB software (version R2018b, The MathWorks, Inc).
In this study, we also examined the effects of restriction on symptomatic infectious individual movement in Belgium with � T� = 5 days. The percentage of symptomatic infectious individuals restricted to travel to other locations, η , was varied from 0%, 20%, 50%, 70%, and 100%. Because infectious individuals partly become asymptomatic with the probability ρ a , therefore, the proportion of travelling infectious individuals, θ , is This means that when all symptomatic infectious individuals are banned from traveling ( η = 100% ), there will be only the mobility of asymptomatic infectious individuals ( θ = ρ a ). On the other hands, when η = 0%, everyone moves normally ( θ = 1).