{'fraction_non_alphanumeric': 0.03175944825601784, 'fraction_numerical': 0.013033782629416155, 'mean_word_length': 4.8345050102303135, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Human Mobility: Models and Applications

Recent years have witnessed an explosion of extensive geolocated datasets related to human movement, enabling scientists to quantitatively study individual and collective mobility patterns, and to generate models that can capture and reproduce the spatiotemporal structures and regularities in human trajectories. The study of human mobility is especially important for applications such as estimating migratory flows, traffic forecasting, urban planning, and epidemic modeling. In this survey, we review the approaches developed to reproduce various mobility patterns, with the main focus on recent developments. This review can be used both as an introduction to the fundamental modeling principles of human mobility, and as a collection of technical methods applicable to specific mobility-related problems. The review organizes the subject by differentiating between individual and population mobility and also between short-range and long-range mobility. Throughout the text the description of the theory is intertwined with real-world applications.


Motivation and History
While the term mobility has multiple connotations, in the context of this review it refers to the movement of human beings (individuals as well as groups) in space and time and thus implicitly refers to human mobility. Indeed, from the migration of Homo sapiens out of Africa around 70,000 years ago, through the European discovery of the "New World", to the existence of expatriate populations in contemporary times, the existence of human beings has always been inextricably linked with their movement. While earlier movement patterns were primarily driven by factors such as climate change, inhospitable landscapes, conflict and food scarcity, in modern times, socio-economic factors such as wage imbalance, differences in welfare and living conditions and globalization play an increasing role.
From hunters and gatherers of prehistoric times to present-day commuters of large metropolitan areas, humans are bound to move on a daily basis in order to earn their livelihood. However, daily trips are also undertaken to perform social and leisure activities. The temporal and spatial scales of these trips are much shorter than those of migratory flows, and they are often characterized by the regularities and periodicities that mark human lives. The daily movement of a large and growing population has important impact on the lives of the individuals and the environmental conditions. Studies conducted in Europe and the United States found that the average household spending on transportation is between 15 and 25 percent of the total expenditures, making transportation the second largest expenditure category after housing. Transportation is also the second source of greenhouse gas emissions to the atmosphere. From these few examples it should be clear that mobility has an enormous impact on human societies and an accurate quantitative description of human mobility is of fundamental importance to understand the processes related to human movement and their impact on the community and the environment. A quantitative theory of human mobility ought to be able to provide answers to relevant questions, such as, what determines the decision to start a trip? which factors determine the choice of the destination? to what extent is human movement predictable, and what is the intrinsic degree of randomness? is it possible to find general rules or laws to explain empirical patterns and regularities exhibited by travels in many diverse countries, such as the distribution of commuting times and distances, and the degree of predictability of future whereabouts?
While the study of human mobility currently spans several disciplines, arguably, geography was the first discipline to analyze mobility data and put forward corresponding theories to describe travel patterns. Indeed, the pioneers of quantitative and theoretical geography in the 50's defined "geography as (the scientific discipline of) spatial interaction" [1]. Early quantitative studies of the movements of people and vehicles were held in large US metropolitan areas [2] (see the first chapters of [3] for details), and initial studies of human travel were of scales ranging from international migrations [4,5,6,7] to journey-to-work commuting [8,9,10,11,12]. Indeed, the elucidation and understanding of these patterns was motivated by its relation to several real-world applications such as traffic forecasting [3,13,14], urban planning [15,10], internal security [16,17] and epidemic modeling [18,19,20] to name but a few.
The first systematic analysis of the concept of distance as a constraint to movement was proposed in the 19 th century: in his 1965 review [21] Gunnar Olsson cites Henry C. Carey's Principles of Social Science (1858) as the first work to explicitly make the observation about the amount of interaction between two cities being proportional to their population size and inversely proportional to the intervening distance. Few decades later, the geographer Ernst Ravenstein further developed and popularized the idea in a seminal work where he formulated his laws of migration [22]. Further refinements on this theme were made in the 1940's by the American sociologist Samuel Stouffer [5,23] in his law of intervening opportunities, and by the American linguist George Kingsley Zipf [6]. Zipf's formulation led to what is now conventionally referred to as the gravity law. The increasing availability of datasets on population movements at various levels of granularity, coupled with the quantitative revolution in geography [24,25,26], led to the introduction of more sophisticated mathematical methods such as hidden Markov models and diffusion processes. Gender and socio-economic factors behind population movement were further analyzed thanks to richer datasets resulting from surveys and interviews [27,28,29], and through theories of labor economics focusing on wage differential between locations [30]. Thus through the 20 th century, contributions to the theories of human movement were made chiefly in geography, sociology, and economics, while the scale at which this was primarily studied was at the population level.
To provide context for what is to follow, we briefly describe a selection of influential historical contributions, keeping in mind "Stigler's law of eponymy" which states: "No scientific discovery is named after its original discoverer" [31]). However, if we restrict ourselves to the recognition of distance as a primary factor in determining movement and interactions between places, then it is reasonable to start with the work of Ernst Georg Ravenstein [22].
Laws of Migration. Ravenstein was a German-English geographer who made important contributions to cartography as well as providing one of the first rough estimates of a "maximal" global population based on resource consumption. He was also one of the first to attempt an explanation and prediction of migration patterns within and between countries. Considering the effect of distance as well as the type of migrant (male or female, old or young) as primary factors, he posited the following seven laws: 1. Most migrants only travel short distances, and "currents of migrations" are in the direction of the great centers of commerce and industry given that these can absorb the migrants. 2. The process of absorption occurs in the following manner: inhabitants of the areas immediately surrounding a rapidly growing town flock to it, thus leaving gaps in the rural areas that are filled by migrants of more remote districts, creating migration flows that reach to "the most remote corner of the kingdom". 3. The process of dispersion is inverse to that of absorption, and exhibits similar features. 4. Each main current of migration produces a compensating counter-current. 5. Migrants traveling long distances generally go by preference to one of the great centers of commerce or industry. 6. The natives of towns are less migratory than those of the rural parts of the country. 7. Females are more migratory than males.
Ravenstein added another two laws in 1889 [32]: 8. Towns grow more by immigration than by natural increase. 9. The volume of migration increases as transport improves and industry grows.
While the laws are non-quantitative and observational in character, Ravenstein correctly identified socio-economic factors as well as distance-constraints to be the essential ingredients behind modern population movement. Consequently, his laws stimulated an enormous volume of work, and although they have been refined and adjusted over the years, the essential ingredients of his formulation remain relevant even today.
Law of Intervening Opportunities. One of the most important refinements was made in the 1940's by the American sociologist Samuel Stouffer [5]. Roughly speaking, Stouffer was looking to expand upon Ravenstein's observations regarding migrants moving shorter distances and flocking to commercial centers. To account for this, he proposed that the number of people going a given distance is directly proportional to the number of opportunities at that distance and inversely proportional to the number of intervening opportunities. In other words, trips between two locations are driven primarily by relative accessibility of socio-economic opportunities that lie between those two locations. In this context, opportunity is defined as a potential destination for the termination of a traveler's journey, whereas an intervening opportunity is one that the traveler rejects in favor of continuing on. In Stouffer's original formulation, this can be mathematically expressed as dy(r) dr where y(r) is the cumulative number of migrants that move a distance r from their original location, and x(r) is the cumulative number of intervening opportunities. Assuming that x itself is a continuous function x(r) of distance, then the expression above can be integrated to yield y(r) = log x(r) + C, where C is some constant denoting the number of opportunities at the origin location. Thus, the relationship between mobility and distance is indirect and is established only through an auxiliary dependence via the intervening opportunities: the higher the number of intervening opportunities between two locations at distance r, the smaller the number of migrants that would travel that distance. This may explain why rural migrants may flock to urban centers over large distances, whereas those already in commercial centers are comparatively more stationary.
Distance-Decay and the Gravity Law. Around the same time as Stouffer, the Harvard Philologist, George Kingsley Zipf, was expanding upon his famous observation of the rank-frequency dependence in linguistics; the eponymous Zipf's law, where the frequency of a word ranked z-in terms of usage-has the statistical dependence f z ∼ 1/z [33]. Zipf found that this relation was expandable to other realms of society, specifically the size of cities [34], where the occurrence of a city with population P and consequent rank z also follows the relation Broadly speaking, Zipf's argument for this relation was due to the tension between two competing factors. The first, which he called Force of Diversification relates to the likelihood of populations living close to the source of raw materials (commodities) in order to minimize the cost of transportation to production centers. The second effect, referred to as the Force of Unification, is the tendency of populations to aggregate in urban centers due to the minimization of work required to transport finished products to consumers. While the former leads to the formation of multiple urban centers (given that the commodity sources are not localized in one part of a country) with smaller populations, the latter has the competing effect of urban agglomeration in a few centers of large population. Assuming some kind of equilibrium between these quasi-forces, the rank-frequency relation of Eq. (3) with α = 1 naturally follows [35]. Deviations from equilibrium, where one force dominates over the other, then leads to a change in the exponent. Carrying the argument further, under somewhat unrealistic assumptions of equitable share of national income as well as urban centers being self-sufficient (i.e. production and consumption at equal levels), the share of any center i in the total flow of goods is proportional to its population P i . Therefore, the flow of goods between two centers is proportional to the product of their populations. Finally, if one would like to minimize the cost and work associated with transportation of goods, this flow must be inversely proportional to the distance between centers. Putting all this together we arrive at the relation where w i j represents the flow of goods between two population centers i, j and r i j is the distance between these two centers. Zipf tested his theory on both freight and population movement data and got good qualitative agreement. This formulation is of course quite different form Stouffer's explanation as the effect of distance is quite explicit in Eq. (4). Indeed, the form of Eq. (4) is such that it leads to a distancedecay effect suppressing long-range movement. Furthermore, the flow of populations and goods is seen to be as a result of some dynamic equilibrium between the cost of transportation, manufacture, and distribution of goods and services. The references to "forces" as well as its functional dependence on distance eventually led to Zipf's formulation being dubbed the Gravity law in analogy to Newtonian mechanics. The common theme connecting Ravenstein, Stouffer and Zipf of course is the geographical distance, though its functional effect on movement is quite different in the Intervening Opportunities and Gravity Law models. Nevertheless, both these models were quite influential on subsequent work, setting off two major strains of parallel research as well as attempts at unification.
Time geography and ICT Data. Measuring, understanding and forecasting the displacements of individuals in space and time has long been part of the program of quantitative and theoretical geography, a branch of geography "born" academically in the 1950's [25,26]. While the first efforts in capturing human displacements focused on the aggregated levels of flows between spatial units, some also focused on individual trajectories. Torsten Hägerstrand, a Swedish geographer, laid down in the early 1950's the basis of time geography, and brought a number of conceptual and graphical tools to formalize the trajectories of individuals through space and time. His seminal work [36] remains famous for its proposal to represent individual trajectories in a cube (also known as the "space-time aquarium"), in which the horizontal plane represents the geographical space, while the vertical axis represents time (as depicted in Fig. 1). Hägerstrand proposed a number of graphical conventions (a "notation system" in his own words) to picture the constraints imposed by social life on individual daily trajectories. He also provided means to represent the co-presence and synchronization of several individuals in space, and more generally a set of (essentially graphical) conventions useful to represent the structure and behavior of individual human mobility.
Time geography was naturally invoked (and somewhat rediscovered) in the 1990's, when the modeling of human mobility shifted towards individual-based simulation (micro/multi-agent/agent-based simulation) [37]. This can be understood in the twin contexts of increased computing power and the development of more expressive programming frameworks, allowing for semantically richer and more ambitious models of human dynamics. However, while the models increased in complexity, they were somewhat artificial; relevant data of comparable complexity and resolution was not available for their calibration. Indeed, while the models progressed, the data lagged behind and the best one had to work with was longitudinal survey data collected since the 1970's.
The beginning of the 21 st century saw the introduction and subsequent widespread adoption of mobile phones, as well as the pervasive usage of Global Positioning System (GPS). This led to an exponential increase in data-generation on human move-  [36]. The geographical space is represented by the 2D plan, while time is figured by the vertical axis. (Left) The two curves represent the daily space-time trajectories of two individuals living in the same neighborhood and working in the same place. (Right) The geographical footprints continuously and passively produced by individuals through the use of their ICT devices allow to approximate their trajectories. While these re-constructed trajectories are partial and contain errors that might mislead the understanding of underlying trajectories, they are nonetheless more precise nowadays than they were 10 years ago, and produced by a constantly growing number of individuals worldwide. ment. Coupled with further progress in computing power and sophisticated datamining methods, it enabled to capture the movement not just of populations at finer levels of spatial granularity, but potentially of individuals. In particular, the large volume and frequency of Call Detail Records (CDRs) from mobile phones (see section 2.3) enabled the analysis of human movement at very fine temporal scales. Thus statistical information on mobility became available, on the scale of hours to decades; from the individual, through communities to the level of country-wide populations.
This review predominantly focuses on these later developments. As we will see throughout the text, these new sources of data have rejuvenated the scientific interest for human mobility, opening the door to new questions and measures, as well as enhancing and validating the insight gained from studying traditional data sources. For reasons which are probably technical, historical, and societal, research teams in physics and computer science have broadly "invested" in these new georeferenced (meta)data resulting from human activity [38]. This is evidenced by a quick search on the online website arXiv 1 for physics papers that include the terms "Human Mobility" or "Mobility Patterns". As seen in Fig. 2, since about 2004-2005, the number of such papers displays an almost exponential increase. One must keep in mind that this is probably an incomplete sample, and is underestimating the volume of research in this field. It is likely that one would see an even more dramatic trend if a more comprehensive list of publications were to be included. While this growing interest of the Physics community in human mobility has multiple reasons, it is somewhat to be expected given the traditional interest of statistical physicists in studying emergent collective properties of systems of many interacting particles and in describing the dynamics of randomly moving particles undergoing (anomalous) diffusion. The availability of data at multiple spatio-temporal scales is helping uncover robust statistical patterns as well as aiding the development of phenomenological theories of human mobility, that are well suited to the tools, methods and paradigms of statistical physics.

Scope and Limitations
The recent pace of developments and the volume of work published are such that a comprehensive review of the findings appeared necessary. Moreover, the study of human mobility is a highly interdisciplinary endeavor and progress has occurred in parallel across different academic communities, sometimes with one unaware of the others' works. This is reflected in different terminology of common metrics and similar results being cast in different and seemingly unrelated contexts. Consequently with this review, we also aim at (i) bringing disparate communities together by unifying the findings in a common context, and (ii) providing new researchers in the field with an accessible starting point and a minimal set of tools, metrics, concepts and models.
One must note that humans, of course, do not move in a vacuum. Just as the diffusion of molecules in materials is mediated by the structural and thermodynamic properties of the material, similarly, spatio-temporal patterns of human movement are necessarily shaped by spatial constraints and limitations of geography. Examples of this are the topography of urban centers or the pattern of roads and transportation infrastructure, the properties of which are studied under the aegis of spatial networks [39]. Yet, this survey is not about Spatial Networks, nor do we give details on the so-called "science of cities" (e.g., city forms, properties of urban transport networks, urban scaling laws, etc. [40]); these two aspects will be considered only when made necessary for the reader to understand the concepts related to human mobility.

Organization
We begin in Sec. 2 with a comprehensive (but by no means exhaustive) list of data sources used for empirical studies. In Sec. 3 we introduce key metrics, measures, spatio-temporal scales, as well as an overview of the fundamental physics behind the study of human mobility. In Sec. 4 we describe the state-of-the-art families of models (both generative and phenomenological) that best describe the empirical observations of human mobility. The models are categorized according to scale, starting from the level of individuals, through to the collective level of populations, and finally a mixture of the two incorporating the concept of (inter) modality. Continuing with the theme of scales, in Sec. 5 we describe some selected applications of the framework and models ranging from intra-urban movement to flows between countries and continents including the case of epidemic spreading, transportation systems, and a brief digression on new results on virtual mobility (web browsing). Finally, we conclude in Sec. 6 with challenges and future directions for the field. We also added Appendix A where we provide descriptions of some basic tools and algorithms (agent-based modeling, random walkers, etc.) that may be of use to researchers making initial forays into the field.

Data Sources
A natural starting point is to describe the nature of empirical data which has been used in mobility research. Indeed, empirical data has been vital to both aggregate and individual mobility investigation, by providing means of parameter calibration as well as model validation. In this section, we outline the main sources available for mobility research and the relevant information that can be extracted from them. The data sources are presented in (rough) chronological order of their historical availability and consequent use in research.

Census Data & Surveys
Census data is collected in periodical national surveys in which householders are asked questions about the socio-demographic and economic status of the household members. Of particular interest in terms of mobility are questions related to the location of the workplace, or the place of current and previous residence. Collectively, this information can then be used to estimate commuting flows or internal migration flows within a country.
Different countries introduced national censuses at different times, and the type of information collected also varies across countries and time periods. In the United Kingdom, for example, the first national census was done in 1841 and contained the names, ages, sexes, occupations, and places of birth of each individual living in a household. In 1921, the location of workplaces was added to the questions asked, and since 1961 the location of previous residence was also included. In the United States, the first population census was taken in 1790, while information on workplace and transportation activities was incorporated in 1940 [41]. Nowadays, national censuses are held in most countries typically every ten years. Below we outline some main sources of census and survey data along with the methodology used to extract commuting and migration flows.

Census Data
The data on commuting trips between United States' counties is available from the United States Census Bureau 2 . These county-to-county work flow files are available for each state as well as at the aggregated level across the country. Furthermore, files are available in either county of residence format-containing work destinations of people who reside in each county-or county of work format-which contain the county of residence of people working in each county. The inclusion of a commuter in the data is preconditioned upon them working during the week leading up to the census. Due to the questionnaire design, any individual not satisfying this condition is precluded from entering a workplace location. A worker is defined as someone who has spent at least 1 hour in paid employment or 15 hours in a voluntary role. Any individual required to provide their workplace location was asked further questions about their mode of transport and travel time.
From this information, the census Bureau compiles and publishes aggregated flows corresponding to the residence and workplace location of everyone classed as a worker. In order to validate models and determine statistics such as trip distance and relationship to the origin, destination and intervening population density, further information on the spatial distribution of counties and socio-demographic data is required. This supporting data is available online 3 . Among other things, such census data is invaluable for validating generative models of commuter flow and migration. For instance, Census data from the 2000s (consisting of county level work flow data, corresponding to 3,141 counties and 34,116,820 commuters within the United States) was used to evaluate the so-called radiation model (described in sec. 4.2.2) by Simini et al. [42] Commuting and migration flows data for the United Kingdom is also available online 4 . The information is available at different spatial resolutions: ranging from the geographical area covered by a single postcode, to areas covered by the 7,201 electoral wards. Files contain origin and destination flows divided into subsets of age and gender allowing for a deeper socio-demographic analysis if required. There have been many different countries whose commuting data has been analyzed at the level of municipalities (or smaller administrative units such as zipcodes or electoral tracks) including France 5 and most European countries. In most cases, when not directly available on an open data portal, commuting data is available to researchers on request from national statistics institutes.

Tax Revenue Data
An alternative source for estimating aggregated migration flows is the Statistics of Income Division (SOI) of the Internal Revenue Service (IRS) in the United States. The IRS maintains records of all individual income tax forms filed in each year. Using individual tax return files, it is possible to determine who has, or has not, moved residence or workplace locations in the intervening fiscal year. To do this, first, coded returns for the current filing year are matched to coded returns filed during the prior year. The mailing addresses on the two returns are then compared to one another. If the two are identical, the return is labeled a non-migrant. Any relevant information change during the fiscal year results in the return labeled as that of a migrant. US migration data from 1992-1993 to 2006-2007 is available online 6 . Alongside tax returns, estimated migrations flows are also included in each Census 7 , based on yearly surveys carried out by the American Community Survey (ACS).

Local Travel Surveys
An important source of data used to construct trips, and therefore flows between two locations, is local surveys. It provides an advantage over censuses as it is possible to collect a more detailed data set that includes information such as the purpose of the trip and mode of transport used. However, the increase in accuracy, resolution, and variety of meta-data leads to a corresponding sacrifice in scale. As such, surveys typically have a much lower number of respondents compared to a census. Additionally, these typically cover smaller areas such as a single city (or indeed neighborhoods within cities), as opposed to the level of a state, or an entire country.
A representative example is provided by the household Travel Tracker Survey for the Chicago metropolitan area, conducted between 2007 and 2008, and carried out by the Chicago Metropolitan Agency for Planning 8 . This survey contains information on travel activities of household members, among other details. Information was collected from a total of 10,552 households over a one or two-day period, and households were required to provide a detailed travel inventory for each member, including details such as trip purpose, transport mode, departure and arrival times, and public transport information such as boarding location, distance to final destination, and fare paid. A similar survey is available for the Los Angeles region 9 and was analyzed by Liang et al. [43], where 46,000 trips between 2,017 zones within the Los Angeles county were extracted. The use of small scale surveys is often augmented by other forms of data (for example GPS tracks, see Sec. 2.4) in order to have accurate records of an individual's position combined with an annotated description of the purposes of each trip [43,44].
While historically most prevalent, census data provides us only with coarse-grained patterns of human migration and movement. Surveys, on the other hand, provide far more detail but on a much smaller scale, thus limiting their use in terms of statistical validation of models. While the census reveals little information about where people spend their time, or where and how they travel, at the survey level, data reliability is restricted by scale and self-reporting errors. Thus both sources lack the ability to provide a dynamic picture of human mobility [45].

Dollar Bills
An unusual (and perhaps not immediately apparent) source of mobility data is related to the tracking of currency. An archetypal example of this are online bill-tracking websites, which are designed to monitor the dispersal of individual bank notes worldwide. The movement of a bank note between geographical locations occurs when it is carried by an individual, therefore data collected from such websites may be used to infer trajectories and patterns of human mobility. Brockmann et al. [46] obtained trajectories of 464,670 dollar bills from the website www.wheresgeorge.com, allowing for an analysis of bank note dispersal in the United States (excluding Hawaii and Alaska). The data set contains entries corresponding to time-stamped reports of the geographic location of individual bills, numbering around 250 million as of 2017. All bills registered on the site are initially entered by a user which involves the input of their  Figure 1 | Dispersal of bank notes and humans on geographical scales. a, Relative logarithmic densities of population (c P ¼ log r P /kr P l), report (c R ¼ log r R /kr R l) and initial entry (c IE ¼ log r IE /kr IE l) as functions of geographical coordinates. Colour-code shows densities relative to the nationwide averages (3,109 counties) of kr P l ¼ 95.15, kr R l ¼ 0.34 and kr IE l ¼ 0.15 individuals, reports and initial entries per km 2 , respectively. b, Trajectories of bank notes originating from four different places. City names indicate initial location, symbols secondary report locations. Lines represent short-time trajectories with travelling time T , 14 days. Lines are omitted for the long-time trajectories (initial entry in Omaha) with T . 100 days. The inset depicts a close-up view of the New York area. Pie charts indicate the relative number of secondary reports coarsely sorted by distance. The fractions of secondary reports that occurred at the initial entry location (dark), at short (0 , r , 50 km), intermediate (50 , r , 800 km) and long (r . 800 km) distances are ordered by increasing brightness of hue. The total number of initial entries are N ¼ 2,055 (Omaha), N ¼ 524 (Seattle), N ¼ 231 (New York), N ¼ 381 (Jacksonville). c, The short-time dispersal kernel. The measured probability density function P(r) of traversing a distance r in less than T ¼ 4 days is depicted in blue symbo It is computed from an ensemble of 20,540 short-time displacements. T dashed black line indicates a power law P(r),r 2(1 þ b) with an exponen b ¼ 0.59. The inset shows P(r) for three classes of initial entry location (black triangles for metropolitan areas, diamonds for cities of intermed size, circles for small towns). Their decay is consistent with the measur exponent b ¼ 0.59 (dashed line). d, The relative proportion P 0 (t) of secondary reports within a short radius (r 0 ¼ 20 km) of the initial entr location as a function of time. Blue squares show P 0 (t) averaged over 25, initial entry locations. Black triangles, diamonds, and circles show P 0 (t) the same classes as c. All curves decrease asymptotically as t 2h with an exponent h ¼ 0.60^0.03 indicated by the blue dashed line. Ordinary diffusion in two dimensions predicts an exponent h ¼ 1.0 (black dashe line). Lévy flight dispersal with an exponent b ¼ 0.6 as suggested by b predicts an even steeper decrease, h ¼ 3.33 (red dashed line).
4 Fig. 4: Trajectories of bank notes originating from four different locations. Tags indicate initial, symbols secondary report locations. Lines represent short time trajectories with traveling time T < 14 days. Lines are omitted for the long time trajectories (initial entry: Omaha) with T > 100 days. The inset depicts a close-up of the New York area. Pie charts indicate the relative number of secondary reports coarsely sorted by distance. The fractions of secondary reports that occurred at the initial entry location (dark), at short (0 < r < 50km), intermediate (50 < r < 800km) and long (r > 800km) distances are ordered by increasing brightness of hue. The total number of initial entries are N = 2, 055 (Omaha), N = 524 (Seattle), N = 231 (New York), N = 381 (Jacksonville). Figure from [46]. ZIP code and a serial number printed on the bill. While the majority of the bills are not entered more than once, around 11% of the bills are reported multiple times, with 3-5 entries per bill being fairly common (Any more than 5 entries is considered rare). Once a bill has been registered, any future hits allow for time and distance between reports to be determined and recorded in the database (see Fig. 4).
While a convenient and large-scale source of mobility data, the use of currency to infer patterns in human mobility is problematic. For instance, bank note dispersals do not contain information on the number of individuals that have carried a given note during two instances of measurement (i.e., when it appears in the database). Consequently, the trajectories of bank notes are likely a convolution of the mobility of several individuals; coupled with the relatively few samples per dollar bill, this makes it a rather inaccurate and problematic source for inferring individual mobility patterns. Nevertheless, similar to census data, currency tracking provides a convenient, and relatively easily accessible picture of the coarse-grained patterns of human mobility.

Mobile Phone Records
Perhaps the most important, game-changing data of the last decade for inferring human mobility are the mobile phone Call Detail Records (CDRs). Most cell phone companies maintain detailed records of customer information related to their calls and text messages. The information contained in these CDRs include the time of the Each dot corresponds to a mobile phone tower, and each time a user makes a call, the closest tower that routes the call is recorded, pinpointing the users approximate location. The gray lines represent the Voronoi lattice, approximating each towers area of reception. The colored lines represent the recorded movement of the user between the towers. Figure from [47].
call/message as well as the company's cell tower routing the communication. Knowing the exact geographic coordinates of the cell tower makes it possible to estimate the approximate location of the user. As mobile phones are (typically) personal devices and are mostly carried by a single person, the corresponding location data can be used to infer a single individual's movements over the recorded period (an example is shown in Fig. 5). Unlike the census, surveys, or currency tracking, CDRs allow for the characterization of individual mobility patterns, and -depending on the datasetat an unprecedented spatio-temporal resolution. Furthermore, due to the tremendous trend of global cell-phone adoption, one can conduct multi-scale studies ranging from the level of neighborhoods to the country-wide level (including international travel and movement across borders).
CDR data concern both phone calls and SMS exchange and always include the following information: time stamp, caller ID, recipient ID, call duration and antenna (cellular tower) code. Due to privacy concerns, customer identifiers are anonymized before the data is made available to scientists for analysis. Quite often, additional clean up of the data is required to make the data more analyzable; for instance, one typically needs to restrict the dataset to users who visit a minimum threshold number of towers, as well as make calls at a high enough frequency to maintain statistical reliability for each user.
One of the first work using CDR's was performed by Gonzaléz et al. [48] on two anonymized datasets provided by a major mobile operator in a large European country. The first set contained the locations of 100,000 randomly selected mobile users whose positions were recorded over a period of 6 months, each time a call or text was received or sent. The second set (collected for validation purposes) corresponded to the locations of 206 mobile users, recorded every two hours over a period of a week. From this data, traces of mobility patterns such as the distribution of displacements between the locations of consecutive calls (jump-lengths) were measured. The temporal period between consecutive calls was also measured to determine a distribution of wait times (i.e., characteristic time spent by users in a given location). The knowledge of users' consecutive locations also enabled the calculation of so-called return time distributions, measuring the probability of a user to return to a given location within a given time. Other more complex metrics, including the radius of gyration were also extracted. (For details on metrics see Sec. 3.1).
The quality of mobility information extracted from mobile phone data depends on the frequency at which the location of the user is recorded. For example, the larger data set used in [48] displayed an irregular call pattern; many calls were placed over a short time period, followed by long periods of inactivity. Consequently, the data displayed a highly inconsistent temporal frequency, which may confound the analysis of mobility patterns. The smaller-scale validation dataset was therefore collected to account for these irregularities. Other problematic issues also arise, such as the accuracy of recorded locations. In general, it is the position of the cell tower routing the call that is recorded and not the exact location of the individual carrying the phone. Since there is significant variability in the areas covered by cell towers-with coverage ranging from tens of meters in the densest neighborhoods of urban areas, to a few kilometers in rural areas-a recorded user in a rural area may transmit all communication via a single tower during their daily routine, while moving the exact same area as that of an urban user who may ping several towers. Thus, in the former case, despite possibly significant movement on part of the user, for purposes of analysis they are considered stationary, which is obviously misleading.
We must note that obtaining CDRs for research purposes poses significant challenges. Since they are not typically publicly available, one has to directly approach a mobile phone operator. As the data contains sensitive information [49], any data collected by a specific group of researchers may not be shared with other groups, making reproducibility problematic. Yet, recent initiatives have seen some phone companies release large scale CDR data within the context of "challenges" among researchers. The purpose of making this data available is to provide researchers some material allowing them to extract useful information to address major societal challenges. Such recent datasets are those that have been provided by Orange in the framework of their Data 4 Development (D4D) challenge. These contain four anonymized tables for mobile phone users in the Ivory Coast (Côte d'Ivoire) and Senegal. For example, in the case of Ivory Coast, data contained (i) hourly antenna-to-antenna traffic; (ii) individual trajectories of 50,000 users over two weeks with antenna location information; (iii) individual trajectories of 500,000 users over an entire year with sub-prefecture location information; (iv) communication graphs for 5,000 users [50]. In particular (iii) contains information regarding individual trajectories for the entire observation period at the cost of reduced spatial resolution; position is recorded as a sub-prefecture rather than antenna location. The raw data consists of a user ID, date and time stamp of the communication and the sub-prefecture number  which can be linked to a file containing the latitude and longitude of the center of each sub-prefecture [51].
Telecom Italia's "Big Data Challenge" is another example of a large scale CDR data set readily available to researchers. Alongside anonymized call records, addi- tional metadata related to weather, news, Twitter activity, and electricity data from two areas in Italy (Milan and Trentino) was provided [52]. Data from this challenge has been used to investigate the relationship between individual daily trips and personal constraints [53], estimate urban population from call volume [54] or determine the relationship between urban communication and happiness [55]. The increasing availability of such large-scale, anonymized datasets for research purposes is an encouraging trend for future research on human mobility. For more detailed information on CDR sources and some associated results (including other research topics than human mobility), see the review by Blondel et al. [56].

GPS Data
The greatest level of accuracy on movement trajectories is provided by Global Positioning System (GPS) data, which is capable of providing precise information on any location covered by at least four GPS Satellites. GPS trackers are units which receive signals from GPS satellites and compute the device's position at regular intervals. This technology allows researchers to trace the movement of individuals with a high degree of accuracy and temporal frequency, thus providing a rich source of data that can be analyzed and directly mapped to human mobility patterns.
Rhee et al. [58] used GPS trackers with a position accuracy of less than three meters to record the movements of over 150 individuals across two University campuses, a state fair (US), Disney World, and New York City. The position of the individuals carrying the trackers was recorded every 10 seconds while they were in these locations resulting in daily traces that had on average a temporal scale of 9-10 hours. Additionally, Microsoft has made GPS trajectories from 182 individuals freely available online via their GeoLife project 10 . The dataset consists of 17,621 trajectories recorded by GPS trackers and GPS enabled mobile devices, and has a high spatial and temporal scale for over 90% of individuals with positions recorded every 1-5 seconds or 5-10 meters for time periods that span several years. Since its release, GeoLife has been used in studies of both human mobility [59,60] and social interactions [61]. For a review of potential applications see [62].
Transmitters attached to vehicles [63,64] are another source of GPS data. These are useful for studying questions related to urban traffic: monitoring, prediction, and prevention of congestion [64]. In Italy, 2% of privately owned cars have a GPS system; each vehicle has assigned a unique identifying number, thus enabling the tracing of anonymous individual cars' trajectories with high precision. A GPS signal is transmitted approximately every 2km, as well as when the engine is switched on and off [63]. Each signal is converted to a data point containing information on position, velocity and distance covered. While the data suffers from inaccuracy due to variables such as satellite coverage and position precision, the temporal precision is of high quality, as is both instantaneous velocity and distance covered. From such GPS vehicle traces, path length distributions can be measured, which corresponds to the distribution of distances traveled in single trips, where a trip is defined as travel that occurs between instances of the engine being started and stopped. Activity distributions, which are equivalent to waiting time distributions in other studies, can also be determined from the elapsed time between consecutive trips. A further advantage of GPS traces over CDRs is that the locations are recorded to a much higher degree of accuracy and at a constant frequency. Yet, there are drawbacks: GPS transmitters do not work indoors and generally rely on batteries as their source of power. As a result, traces can contain periods during which no signal is received or may be terminated prematurely. Furthermore, GPS datasets typically feature a smaller number of individual users, up to several thousand, in comparison to mobile phone data which can provide information on the mobility of millions of users.

Online Data
Another valuable source of location information are Online Social Network (OSN) or Location-Based Social Network (LBSN) services, that attract hundreds of millions of users worldwide. Since the introduction of GPS and wifi chips in smartphones, social network providers have been able to collect valuable data on both the social connections as well as precise geographical locations of their users. Indeed, services such as Twitter, Facebook, Foursquare and Flickr collect geotagged data every time a user enables localization for the content being posted (e.g., checking-in at a restaurant with friends); this is associated with geographic coordinates, a time stamp, and additional information relative to the location itself or the content published by the user. The mobility profiles of the users can be obtained from the list and number of locations other than their home which they have visited over the period of the study (see Fig. 7 for an example). These profiles can range from intra-urban routine trips to worldwide travel [65,66]. From the geographical locations, by inferring a user's home location, their radius of gyration can be determined using the geographical distance between the home location and reported locations. The jump length distribution can be determined by defining a jump as the geographical distance between consecutive reports for a single user and frequency of reports can be used to assess temporal patterns in human mobility [66]. Other measures such as visitation frequency and predictability are also obtainable [67].
Geotagged data can also be used to assess social interactions between individuals. Scellato et al. [68] analyzed three online location-based social networks: Foursquare, Brightkite, and Gowalla, finding their exhibition of common properties found in most other real world complex networks, such as fat tailed degree distributions, high clustering coefficient, and short average path length between nodes, confirming the smallworld nature of LBSN's [69]. Furthermore, the probability of a social link between two users was found to exhibit a distance dependence, decreasing with geographical separation between the users. Compared to CDR's, data collected from OSN's and LBSN's, have the advantage of having more contextual information associated with the geographical positions and the users, enabling the study of mobility in a broader context. For example, data collected from Twitter users has been applied to topics such as social networks [70,71,72], evolution of moods [73,74,75], and crisis management [76,77,78] among other features.
Yet, as in all other data sources, limitations persist. Besides requiring extensive clean up (such as removing users with suspiciously high activity, unreasonably fast movement between two consecutive check-ins or incomplete data [66]), mobility and interaction data originating from online social networks needs further validation in order to be considered truly representative of the general population. Sloan et al. [79] found that Twitter data originating from the UK has significant demographic differences compared to the wider population data from UK census. Furthermore, there appear to be demographic differences between the users who enable geotagging and those who do not (only ≈ 3% of the users were found to have geotagging enabled).

Metrics, Physics and Scales
In this section, we will discuss some of the fundamental metrics used to characterize mobility as well as the associated spatio-temporal scales at which they are relevant. We will then move onto a discussion of some of the physics associated with mobility, including the relations between distance, time and velocity. We end the section with a discussion of energy arguments and interpolation of spatio-temporal scales through the lens of multimodality.

General Metrics 3.1.1. Jump Lengths
A key factor in modeling human mobility is the distance an individual travels in a given time period. Measures of distance are often dependent on the source of data being used, and the terms flight length, jump length, displacement, and trip refer to different distance measures that may be extracted from data. For example, early measures of mobility leveraged information about the spatial trajectories of bank-notes (one can think of this as an aggregation of many individual trips as a given banknote necessarily changes many hands). Conventionally, the distance between two instances of the appearance of a banknote in the measured data is termed a jump length [46]. Later, higher resolution measurements of movement were provided by Call Detail Records (CDRs). Here one could reliably measure the location (and corresponding displacement) of individuals based on placement of mobile phone towers that are pinged when one makes a call. Associated with the displacement between two successive calls is a time interval [48], that provides an estimate of the stay of an individual in a location (waiting-times), although it is not possible to detect the position of the user between two consecutive calls. Even higher resolution spatio-temporal data available from Global Positioning System (GPS) allows one to make more stringent definitions of displacement. In addition to the standard definition of jump-lengths [58], one can now define a displacement between stops, i.e. the displacement of individuals between two locations, given that they spent a minimum amount of time per location [63]. Indeed, the use of the terms "stop", "trip" or "displacement" reflects human behavioral tendencies that motivate people to go from point A to point B. It is therefore important to verify that a "stop" corresponds to an actual behavior and that it is not artificially generated by oversampled or undersampled signals [80].
Regardless of how one chooses to constrain it, the jump length, typically denoted as ∆r, is defined as the euclidean distance between r(t) and r(t + dt) corresponding to locations recorded at intervals t and t + dt. Of particular interest is the distribution of ∆r within a population; how likely is it that a random member of the population will travel a distance r from their origin location in a time dt? In order to measure this, when modeling human mobility it is common to consider the probability distribution function (PDF) of jump lengths, P(∆r). This may be defined as the probability of finding a displacement ∆r in a short time step dt.
Brockmann et al. [46] determined P(∆r) for the trajectories of dollar bills and the distribution of jump lengths was observed to follow a power law, P(∆r) ∼ ∆r −(1+β) , with exponent β = 0.59 ± 0.02, independent of the size (in population) of the entry point of a bank note (see Fig. 8). The power-law behavior of P(∆r) has also been observed in the trajectories of mobile phone users [48,81]. The empirical distribution of displacements obtained from mobile phone data analyzed in [48] is well approximated by a truncated power-law distribution: where r 0 and κ represent cutoffs at small and large values of ∆r. The value of the observed scaling exponent in Eq. (5) is β = 0.75 ± 0.15, not far from β = 0.59 ± 0.02 obtained from bank note dispersal, suggesting that the two distributions may capture the same fundamental mechanism driving human mobility patterns. Zhao et al. [82] measured the distribution of jump lengths using two GPS data sets that, along with location and time stamps, also included transportation mode. From this, the jump length ∆r was determined as the longest straight line distance between two locations without a change of direction. Wait times, ∆t, were taken to be the time spent in a particular transportation mode. A single trip may be made up of several flights, each of which have a corresponding transportation mode taken directly from the dataset. The transportation modes were grouped into four categories; Walk/Run, Car/Bus/Taxi, Subway/Train and Bike. The distribution of jump lengths is found to be log-normal if each transportation mode is considered individually, whereas P(∆r) for all modes combined follows a truncated power-law with β = 0.55 and 0.39 for the two datasets. The distribution of wait times, P(∆t), was found to be exponential, corresponding to a large number of Walk/Run flights that connect the use of other transportation modes within a single trip and therefore have only a short duration. This result provides further insight into studies using on bank note and mobile phone data [46,48,81] in that the power-law behavior of the distribution of jump lengths is the result of combining several distinct log-normal distributions, corresponding to the different transportation modes available, with exponentially distributed wait times. GPS devices installed in cars transmit the location of the car to a high degree of accuracy every few seconds when the car is in motion. The analysis of such data from over 150,000 cars in Italy over a period of one month [64] found that the probability density function, P(∆r) had two different regimes; an exponential distribution up to a characteristic distance of ∼20 km corresponding to inter-city travels, and a power-law with β = 1.53 corresponding to intra-city movement. Here, ∆r corresponds to the true length of a single trip. Trips are distinguished by records that are greater than 20 mins apart; as the GPS device only transmits when the car is moving, using a threshold of 20 mins allows for small stops, such as waiting at traffic lights or re-fueling, to be accounted for and included within a single trip. The exponent, β, is significantly higher than those obtained for bank note braic ail in atial ion 15 which clearly shows that dispersal cannot be described by a pure Lévy flight. The measured decay is even slower than the decay expected from ordinary two-dimensional diffusion (h ¼ 1, dashed black line). Therefore, we conclude that the slow decay in P 0 (t) reflects the effect  dispersal [46] and mobile phone data [48,81]. The authors suggest this is due to the limitation in the GPS data caused by restrictions on the size of the area covered (500km in length), and the physical limitation on the distance an individual will drive in a single trip.

Mean Square Displacement (MSD)
Research into the distribution of jump lengths suggests that individual trajectories can be described by Lévy flights, a family of models associated with random walks (described in detail in Sec. 4.1). In this context, a common measure of the potential exploration of area by an individual is the Mean Square Displacement (MSD): Here r 0 is a vector marking the origin of the individual relative to some reference point, or in other words, the location at which a particular trajectory starts, while r(t) mea-  In particular, for a random walk (ordinary diffusion) ν = 1. Yet a CTRW type process is not entirely a realistic representation of how humans actually move. In particular, a random walker will tend to drift away from the origin of its trajectory rather rapidly, the longer the elapsed time. Whereas, as one would imagine, empirical measurements indicate that people have a tendency (on average) to return home on a daily basis [48]. Furthermore, the means of transportation constrain the maximum jump length, therefore it has been observed that the PDF's P(∆t) and P(∆r) have exponential cutoffs; this implies that their variances are finite and hence in the long-time limit the scaling of the MSD will asymptotically converge to that of Brownian motion, MSD(t) ∼ t 1/2 .
To account for this discrepancy, refinements were made in [81], where analysis was restricted to high-resolution data with a combination of mobile phone records that provided locations each time a call was made, along with location data recorded by mobile services at hourly intervals. The newer measurements revealed that human mobility apparently follows an ultra-slow diffusive process characterized by a slower than logarithmic growth of the MSD with time, MSD(t) ∼ (log log(t)) 2/β (see Fig. 9).
Yet, the picture is typically more complex. Indeed, a common problem arising when calculating the MSD, is that the location of an individual's origin is often ambigu-ous [58]. One method of overcoming this issue is to take the average of MSD values measured by varying the origin among all locations that an individual visits [83,84]. Such an approach was taken in [58] where they analyzed the GPS traces of individuals in a number of different location, including a university campus (KAIST), New York City, Disney World, and a US state fair. The averaging procedure applied to the MSD revealed two types of temporal behavior; up to a certain time (∼ 30 minutes) the participants moved super-diffusively (as expected in a CTRW process), whereas for periods greater than 30 minutes, the individuals moved sub-diffusively. It is here that the notion of scales becomes important. The authors explain this dichotomous behavior by showing that at short time scales, the distribution of jump lengths follow a power-law distribution, whereas at longer time scales, jump lengths are more homogeneous, following a Gaussian distribution. While the combination of power law waiting times and jump lengths lead to super-diffusive behavior at shorter spatio-temporal scales, the Gaussian distribution of jump lengths (while maintaining a power law behavior for the waiting times) leads to mobility that appears sub-diffusive [85] at larger spatio-temporal scales.

Radius of Gyration
The sub-diffusive behavior of humans at certain scales suggests that they tend to move a characteristic distance away from their starting locations. This distance can be quantified by the so-called, radius of gyration, r g , defined as the root mean square distance of a set of points from a given axis. Traditional applications of this measure are rooted in physics (where it is related to the mass moment of inertia) and engineering (distribution of cross sectional area in a column). A general formulation of r g is given by: where r i are the coordinates of the N individual points (or N measurements of location) and r 0 is the position vector of the center of mass of the set of points, r cm = N i=1 r i /N. When applied to human mobility, the radius of gyration can be used to characterize the typical distance of an individual from the center of mass of their trajectory. The radius of gyration is not equivalent to the square root of the Mean Squared Displacement, defined in Eq. 7, although the two quantities are related because the former is a central second moment while the latter is a second moment about the origin.
Gonzalez et al. [48] used mobile phone call records to determine the radius of gyration for two sets of users; one set with high temporal resolution (position recorded every two hours over one week), and one set with lower resolution on a longer (6 month) scale. The radius of gyration was calculated for each user, a, up to a time t using the recorded positions, r (a) i , i = 1..., N (a) . The value of r 0 was taken as the center of mass of each users trajectory; r a cm = N (a) i=1 r (a) i /N (a) . The distribution of r g , P(r g ), was found to follow a truncated power law of the form: where r (0) 0 represents the minimum radius of gyration cutoff due to the spatial sensitivity erns. a, Week-long trajectory of 40 st individuals travel only over short ver hundreds of kilometres. b, The e different phone towers are shown as grey marks the approximate reception d by us records only the identity of the e can not identify the position of a user of the user shown in b is constructed which the user visited a total of 12 Among these, the user is found on 96 erred locations; the frequency of visits for each location is shown as a vertical bar. The circle represents the radius of gyration centred in the trajectory's centre of mass. c, Probability density function P(Dr) of travel distances obtained for the two studied data sets D 1 and D 2 . The solid line indicates a truncated power law for which the parameters are provided in the text (see equation (1)). d, The distribution P(r g ) of the radius of gyration measured for the users, where r g (T) was measured after T 5 6 months of observation. The solid line represents a similar truncated power-law fit (see equation (2)). The dotted, dashed and dot-dashed curves show P(r g ) obtained from the standard null models (RW, LF and TLF, respectively), where for the TLF we used the same step size distribution as the one measured for the mobile phone users. of the data, and κ r represents an upper cutoff mostly due to the finite size of the study area (see Fig. 10).

Nature Publishing Group ©2008
Indeed, the measured value of the scaling exponent β r = 0.65 ± 0.15 indicates a significant degree of heterogeneity in the travel habits of the observed population. Measuring the conditional jump length distribution, P(∆r|r g ), revealed that users with small r g travel mostly over small distances, whereas those with large r g tend to display a combination of many small and a few larger jump sizes. Once one accounts for this heterogeneity in travel habits by rescaling the distribution with respect to r g , it leads to a collapse of the data onto a single curve thus, Here β c = 0.2 ± 0.1 is a "universal" scaling exponent, and F(x) is a scaling function that is constant for x < 1 and rapidly decreasing for x 1.
It is important to note that the distributions described above are not independent but are related by the equality P(∆r) = ∞ r 0 P(∆r|r g )P(r g )dr g , where P(∆r) is the jumplength distribution introduced previously. If P(∆r) has a power-law scaling exponent β, then we have β = β c + β r , which seems to be in good agreement with empirically measured values. This form of scaling corresponds to a type of random walk called a Lévy flight (Sec. 4.1.2), and the results suggest that this may be the behavior of individuals up to their associated characteristic distance r g (and beyond which saturation effects take over).
A way to uncover the saturation effect is to examine the time evolution of the radius of gyration, i.e. r g (t). In [48], it was found that the average r g of mobile phone return probability for each individual F pt (t) (first passage time probability) 21,22 , defined as the probability that a user returns to the position where he/she was first observed after t hours (Fig. 2c). For a two-dimensional random walk, F pt (t) should follow ,1/(t ln 2 (t)) (ref. 21). In contrast, we found that the return probability is characterized by several peaks at 24 h, 48 h and 72 h, capturing a strong tendency of humans to return to locations they visited before, describing the recurrence and temporal periodicity inherent to human mobility 23,24 . remaining time in 5 to 50 places, visited w Therefore, the observed logarithmic satur the high degree of regularity in the daily tra captured by the high return probabilities frequented locations (Fig. 2d).
An important quantity for modelling h the probability density function W a (x, y) t given position (x, y). As it is evident from F travel in different regions, yet each user defined area, defined by home and workp be found most of the time. We can comp ferent users by diagonalizing each trajecto ing the probability of finding a user in a giv the user's intrinsic reference frame (see Su for the details). A striking feature of W (x, anisotropy in this intrinsic reference fram in Fig. 3a); we find that the larger an indi nounced is this anisotropy. To quantify t anisotropy ratio S ; s y /s x , where s x and deviation of the trajectory measured in th frame (see Supplementary Information). W monotonically with r g (Fig. 3c), being S*r {g g for g < 0.12. Given the small valu other functional forms may offer an equa nistic models are required to identify if th law or only a reasonable approximation to To compare the trajectories of differen individual anisotropies, rescaling each respective s x and s y . The rescaledW W À ( Fig. 3b) is similar for groups of users w r g , that is, after the anisotropy and the r g all individuals seem to follow the same u lity distribution. This is particularly evid show the cross section ofW W x=s x ,0 ð Þ users, finding that apart from the noise i indistinguishable.
Taken together, our results suggest that t in bank note measurements capture a conv heterogeneity shown in equation (2) and users. Individuals display significant regu to a few highly frequented locations, suc regularity does not apply to the bank no the trajectory of its current owner; that i humans do not.
The fact that individual trajectories same r g -independent two-dimensional W W x=s x ,y s y À Á suggests that key statistic vidual trajectories are largely indisting Therefore, our results establish the basi agent-based models, requiring us to place tional with the population density of a giv user an r g taken from the observed P(r g predicted anisotropic rescaling, combined W W x,y ð Þ, the shape of which is provided a we can obtain the likelihood of finding a u the known correlations between spatial p our results could help quantify the role of ment and evolution [25][26][27][28][29] and improve our u processes 8,30 . AEr g (t)ae versus time for mobile phone users separated into three groups according to their final r g (T), where T 5 6 months. The black curves correspond to the analytical predictions for the random walk models, increasing with time as AEr g (t)ae | LF,TLF , t 3/2 1 b (solid curve) and AEr g (t)ae | RW , t 0.5 (dotted curve). The dashed curves corresponding to a logarithmic fit of the form A 1 B ln(t), where A and B are time-independent coefficients that depend on r g . b, Probability density function of individual travel distances P(Dr | r g ) for users with r g 5 4, 10, 40, 100 and 200 km. As the inset shows, each group displays a quite different P(Dr | r g ) distribution. After rescaling the distance and the distribution with r g (main panel), the different curves collapse. The solid line (power law) is shown as a guide to the eye. c, Return probability distribution, F pt (t). The prominent peaks capture the tendency of humans to return regularly to the locations they visited before, in contrast with the smooth asymptotic behaviour ,1/(t ln(t) 2 ) (solid line) predicted for random walks.  Main After rescaling the distance and the distribution with r g , the different curves collapse to a power law (solid line). (c) Return probability distribution, F pt (t). The prominent peaks capture the tendency of humans to return regularly to the locations they visited before, in contrast with the smooth asymptotic behavior ∼ 1/(tln(t) 2 ) (solid line) predicted for random walks. (d) A Zipf plot showing the frequency of visiting different locations. The symbols correspond to users that have been observed to visit n L = 5, 10, 30, and 50 different locations. Denoting with L the rank of the location listed in the order of the visit frequency, the data are well approximated by R(L) ∼ L −1 . The inset is the same plot in linear scale, illustrating that 40% of the time individuals are found at their first two preferred locations; bars indicate the standard error. Figure from [48]. users displays a logarithmic increase with time, r g (t) ∼ A + B ln t (see Fig. 11), in contrast to a pure random walk where the scaling is of the form r g (t) ∼ t θ . The observed saturation in r g can be attributed to the regularity in travel patterns of the individuals, and specifically to the high probability for an individual to return to a few highly-frequented locations. As pointed out before, the true picture is a mixture of behaviors. Supporting the conclusions made in [58] while measuring the MSD, Zhao et al. [86] argue that humans engage in a mixture of superdiffusive and subdiffusive | Population heterogeneity and ultraslow growth of the radius of gyration. a, The P(r g ) distribution of the radius of gyration r g for mobile-phone users at different moments of time. The straight line, shown as a guide to the eye, represents a power-law decay with the exp 1 + ↵ ⇡ 1.55. b, The time evolution of the radius of gyration r g (t) for the mobile-phone users, where r g (t) is normalized by its maximal value T = 1 year. The orange curve is a guide to the eye, following (log t) 2 . c, The distribution of the radius of gyration P(r g ) for the individual-mobi ↵ = 0.75, = 0.6, = 0.2 and ⇢ = 0.1, the values found to be of direct relevance to human mobility. d, The time evolution of the radius of g the individual-mobility model with the same parameters as in c, where r g (t) is normalized by its value r g (T) at T = 10 6 . The green curve is a eye, following (log t) 2 .
model that assigns a quenched visitation variable to each site. Our approach not only avoids the need to parameterize such a model, but also achieves the trajectory selection dynamically, through its self-quenching character.
We also note that the model is designed to capture the long-term spatial and temporal scaling patterns; thus, in its present from it does not reproduce the short-term temporal order and correlations potentially present in individual mobility. Our choice to focus the asymptotic properties is driven not only by theoretical arguments (we aim to reproduce the universal and not the transient patterns), but also by practical considerations: many human-mobility-driven processes, from epidemic spreading to city planning, are driven by the asymptotic characteristics of human mobility. To achieve a as well potential correlations in spatial mobility (that B is between locations A and C in space, the likely orde will be A ! B ! C or C ! B ! A). These correla constrain the human trajectories, being partly respo high degree of predictability characterizing individ patterns 39 . Finally, we note that the dynamical quen recurrent behaviour are unique to human trajectorie restrict banknote diffusion or foraging behaviour 13, our model represents an improvement over the CTR models, as it is adapted to capture the specifics of hum By reproducing the basic scaling laws character trajectories, the present model offers a conceptual fr has the flexibility to absorb future extensions, potentia behaviors (depending on temporal scale), and the latter specifically can be explained as a consequence of the saturation of r g . This was later reproduced by Song et al. [81] and demonstrated to be the consequence of the fat-tailed distribution of the jump lengths (Fig. 12).
The radius of gyration is a measure to characterize the typical distance traveled by an individual and it depends both on the mutual distance of the locations visited and on the time spent (or the total number of visits) in each location. However, the radius of gyration does not allow us to quantify the relevance of each location in determining an individual's characteristic mobility. Indeed, an individual who spends a majority of time in their most visited locations, e.g. home and work, will have a large r g if these two locations happen to be quite far from each other. Conversely, even if these most visited locations are close to each other, large values of r g may be reported if the individual happens to visit a number of distant locations.
To disentangle these effects, one can study the influence of frequency of location visits on the characteristic distance traveled, by studying the k-radius of gyration, r (k) g , defined as the radius of gyration computed over an individual's k most frequently visited locations L 1 , . . . , L k (see next section) [87]. One can express this modified form as where N k is the sum of the visits to the k most frequented locations, r (k) cm is the center from the population (see Supplementary Note 4, Supplementary  Fig. 2). The three methods produce similar results, indicating that the two classes are clearly separated and well-defined.
Consequently, in the following we method to split the population into k-The ratio s k ¼ r ðkÞ g =r g measures the recurrent mobility on her overall mobi higher is the weight of the top k locatio individual. Figure 4 shows the probab ratio for different k. We observe two s k ¼ 0 corresponds to k-explorers, who smaller than the total r g ; the peak at k-returners, individuals whose r ðkÞ g is v Note that only for a few individual suggesting that for the great majority most frequented locations are on av of mass than their remaining less Supplementary Note 9, Supplementary the k-explorers gradually become explorers and returners peaks to respectively. The population reaches and k-explorers for k ¼ 4 for GSM. In k, we always have more k-returners th reason is that GPS data only contai vehicles, hence missing long dist frequented by a particular individual, These trips increase the total r g w Neglecting these trips results in a lower total r g , increasing the chance to class Returners and explorers are also ch spatial distribution of the visited locat representative examples of individua two-returners and two-explorers with profiles, the visited locations tend to g The mobility networks of returners and explorers for k ¼ 2. Nodes (circles) indicate the geographic locations visited by the individual, and each link denotes a travel observed between two locations. When the total r g is small, the two most important locations (red and blue) are close to each other for both two-explorers and two-returners. As the total radius increases the behaviour of two-returners and two-explorers starts to differ; for two-returners, the two most important locations move away from each other; for two-explorers, they stay close and other clusters of locations emerge far from the centre of mass (the grey cross). The mobility networks of returners and explorers for k = 2. Nodes (circles) indicate the geographic locations visited by the individual, and each link denotes a travel observed between two locations. When the total r g is small, the two most important locations (red and blue) are close to each other for both twoexplorers and two-returners. As the total radius increases the behavior of two-returners and two-explorers starts to differ; for returners, the two most important locations move away from each other; for explorers, they stay close and other clusters of locations emerge far from the center of mass (the gray cross). Figure from [87].

2-Explorers
of mass computed on those locations and n j is the number of visits to the j-th most visited location. As an example, if r (2) g r g , then the characteristic traveled distance is dominated by the two most frequented locations, whereas if r (2) g r g than the two most frequented locations do not offer an accurate characterization of the individual's travel pattern, requiring us to consider more locations.
The analysis of GSM and the GPS data revealed that there exist two distinct classes of individuals, returners and explorers. The characteristic distance traveled by kreturners is dominated by their recurrent movement between a few preferred locations, and their radius of gyration is well approximated by their k-radius of gyration for k ≥ 2. In contrast, k-explorers have a tendency to wander between a varying number of different locations and their k-radius of gyration is very small compared to their overall r g (see Fig. 13).

Most Frequented Locations and Motifs
When considering patterns in human mobility, particularly movements within a single day or week, it is essential to distinguish between locations based upon their importance. As already mentioned, people have a tendency to return home on a daily basis and therefore most daily and weekly trajectories will start and finish at the same location. One method of quantifying the importance of a location is the use of ranks; the most visited location (likely home or work) would have rank 1, a school or local shop may have rank 2 or 3, etc. In [48], the rank of a location was determined for each individual mobile user by the number of times their position was recorded in the vicinity of the cell tower covering that location. It was found that visitation frequency follows a Zipf law, that is the probability of finding a user at a location of rank L is approximately P(L) ∼ 1/L.
Another method of distinguishing between locations is to construct each individual's mobility pattern as a network. Schneider et al. [44] used data from both mobile phone users and travel survey respondents to construct weekday mobility networks for each individual. These profiles consisted of nodes to represent locations visited and directed edges to represent trips between the locations. Every daily trip network started and ended at the home location of each user, which was determined as the location they were most often recorded at between 3:00am and 3:30am. For mobile phone users, only days where calls were made, and therefore location was recorded, during 8 or more 30min time slots were included. Daily networks were constructed for weekdays only in order to identify patterns in mobility during a typical day.
It was found that ∼ 90% of the recorded trips made by all users can be described with just 17 daily networks. These 17 trip patterns can be described as motifs; a subnetwork within a complex network [88]. In this case, a daily network was considered to be a motif if it occurred more than 0.5% in the datasets. The result that the daily patterns of human mobility can be constructed by just 17 trip networks suggests that these motifs represent the underlying regularities that exist in our daily movements and are therefore useful for the accurate modeling and simulation of humans' mobility patterns.

Origin-Destination Matrices
The Origin-Destination (OD) matrix T is the standard object in aggregated mobility studies and transport planning. It provides an estimate of the number of individuals traveling between locations in a given area, over a given period of time. More precisely, an OD matrix is a n × m matrix where n is the number of different "Origin" zones, m is the number of "Destination" zones, and T i j is the number of people traveling from zone i to zone j. More commonly, an area under study is partitioned to an equal number of origin and destination points and therefore, n = m. The size n of the OD matrix then depends on the spatial scale/resolution at which the data has been collected. Traditionally, zones are administrative units, whose size may vary from census and electoral units to entire municipalities, departments or states, depending on the question that motivated the development of the OD matrix. Obviously, the maximal spatial resolution of an OD matrix depends on the data source.
An OD matrix can be empirically derived, being estimated from household or roadside travel surveys, traffic counts, and more recently from individual digital footprints (see below). It can also be the output of a model, like in the classic 4-step transport model widely used in urban planning 11 . OD matrices and the 4-steps framework have been used in transport planning since at least the middle of the 20 th century, and there exists an enormous literature that covers its origins, uses, models, synthetic indicators, and so on [89].
OD matrices are of specific importance to models of aggregate flow, that is models of human mobility at the population level, rather than individual level. Examples of such models are gravity and intervening opportunities models (Sec. 4). In empirical OD matrices, T i j indicates the number of travelers from i to j measured during observation time window, whereas if the OD matrix is a model's estimate or prediction T i j usually indicates the expected (average) number of travelers between the two locations. The OD matrix diagonal elements are usually 0 to indicate that people only move between distinct locations. In some studies the diagonal elements are larger than zero and indicate the number of non-traveling individuals in the origin location (or those who travel within the origin location). The fraction of travelers from the origin i to all other locations can be calculated as In this case, the entry p i j may represent the probability of an individual located at i to select location j as their destination over all other possible locations.
Until the 2000's, the traditional approaches to estimate OD matrices consisted in relying on travel surveys or counting. The drawback of such approaches are related to their cost to setup, frequency of updates, and that they cover a limited sample of individuals (a few hundreds or thousands of households) [90]. Consequently, to circumvent these difficulties, a large body of literature has dealt with extracting OD matrices from recently available, individual digital footprints [91,92,93,94,95,90,96,97,98]. Early work on this topic has already been summarized in a dedicated review a few years ago [99]. Here we don't go deeply into the details of the numerous methods that have been proposed, and give only the general picture for the case of average daily journey-to-work commuting.
Starting from time-sequences of successive locations, the goal is to identify important places ("anchor points") for each individual, such as the residence and workplace. The heuristic underlying method is straightforward: due to the circadian regular daily rhythms of human activity, one can make the reasonable assumption that for most individuals, the most frequent location of mobile devices will be very near the residence during non-working hours of weekdays, and during weekends. Similarly, during working hours on weekdays, it is likely that for most individuals their devices will be located just next to their workplace or location of primary activity (which include students, unemployed, retired, and employed individuals with a non-fixed workplace) [20,96,97,98]. In some cases validation of methodology was conducted by asking a small group of volunteers to complete a survey, whose results were then compared to the key locations extracted from the recorded activity of their mobile devices [93].
Alexander et al. [97] proposed a method for constructing an OD matrix from mobile phone records. Such a method allows for data at the individual level to be aggregated in a way that lends itself to analysis at the population level. Implementing the proposed method, in order to infer a user's origin and destination, clustered locations are extracted from the mobile phone record. By clustering locations in close proximity to one another into a single origin or destination, noise from the data, such as inexact triangulation, is eliminated thus allowing for a more accurate analysis. The clustered locations are then assigned a type: work, home, or other. The assignment takes into account the time of day a user is observed there, the duration of their stay, and the day of the week. By the nature of mobile phone records, location is only recorded when a user makes or receives a call, therefore the arrival times and duration of stay at a given location may not reflect the true times. To correct for this, surveys on trips may be used to assign true arrival times and duration from a probability distribution derived from the survey data. Then, trips can be constructed for every user between each location at which they are observed. In order to ensure that the final OD matrix is representative, users that do not satisfy certain conditions, such as on average of at least one visit to the home location per day, may be removed from the data set. The daily trips are then combined for all users and locations allowing for the aggregate flow between all possible locations to be determined and an origin destination matrix to be constructed.
Lenormand et al. [96] compared the journey-to-work OD matrices extracted from three different data sources (census, mobile phone CDR data, and Twitter data) in Madrid and Barcelona, two of the largest European urban areas. Projecting residential and workplace locations identified in each dataset on regular square grids, in order to have a common spatial resolution for each dataset, they measured the correlation between the OD matrices estimated through each data source, as a function of the level of spatial aggregation. They found good correlation for square grids of size 2km (ρ > 0.9), and very good correlation between the various sources when aggregating the individual data sources (Twitter, CDR data) on the map of administrative units (municipalities) used by the travel surveys (ρ > 0.99) (here ρ is the Pearson correlation coefficient). Using a slightly different method that aggregates nearby cell towers when identifying important locations, Alexander et al. [97] estimated OD matrices of daily average commuting trips from CDR data in the Boston area, and then performed a sensitivity analysis of the correlation between this OD and the OD issued from transport surveys at different levels of spatial aggregation. Tizzoni et al. [20] also compared OD estimated from census and mobile phone data at different scales of spatial aggregation. They found lower correlation values, but applied a simpler method to estimate the OD from mobile phone data, when the aforementioned papers applied more stringent filters to remove outliers, fusion phone antennas close to one another, etc. Several algorithms in the literature are reviewed in [98], and further successive refinements and improvements in the methodology have been made since [95,90,100,97].

Physics of mobility 3.2.1. Distance, travel time, and effective speed
In essence, the metrics described above, deal with the concepts of distances traveled and the time elapsed within and between journeys. Another important feature, which has not been discussed, is the notion of speed. To understand the relation between these quantities, it is important to note that mobility occurs over a large variety of distances. Framing it in the context of urban agglomerations, intra-urban mobility occurs over distances typically in the range of 1-10 kms, while inter-urban displacements occur generally over a wider range of the order of 100 kms (or more for large countries such as the US), and inter-country and intercontinental trips cover distances typically of the order of 1000 kms. While to first order, one may think of time spent on a given trip being roughly proportional to the distance traveled, one has to take into account the mode of transportation, which itself depends on the distance. For shortrange travel, slow transportation modes with many stops are utilized (such as walking or public transportation), and for longer distances, one typically takes fast trains or planes with comparatively fewer stops. Thus the relation between distance, traveltimes, and therefore speed, is fairly nuanced, as can be seen in Fig. 14, where we show the relation between apparent speed (i.e. the geodesic distance divided by the travel time) and the travel distance [101].
In general, a trip can consist of multiple connections and may in fact be multimodal (multiple components of travel) with the corresponding waiting-time distributions associated with the modes of transportation [102]. Also, it is observed that the apparent speed v increases with travel distance according to the following power-law functional form [101] v ∼ r β , where β ≈ 0.5. This particular dependence is primarily due to a combination of the hierarchical structure of transportation systems [103] and the fact that waiting-times (parking, take-off, landing, etc) decrease in proportion to trip distance. For example, it has been observed for public transportation as well as personal cars [103] that a plot of the average trip velocity as a function of the trip duration reveals an effective acceleration (Fig. 15) of the form This effective acceleration is a direct consequence of the hierarchical organization of roads: the longer the trip, the more likely that higher velocity roads such as highways are used, thus trips can be decomposed into a ascending cascade to faster roads, followed by a descending cascade when approaching the target location. Eq. (12) implies that the distance evolves with trip duration as from which the relation v ∼ √ r directly follows [101].

Travel time budget
In the context of intra-urban mobility, Marchetti [105], based on ideas developed by Zahavi [106], proposed the idea of a travel time budget of about one hour per day irrespective of location. This implies that with improvements in transportation technology and corresponding increase in speed, a greater amount of distance is covered within the budgeted time, thus allowing for urban sprawl. This assumption can be reformulated as the rational locator hypothesis [107] which posits that individuals maintain approximately a constant journey-to-work travel times by adjusting their home and workplace. This was revisited in [107] on data for travel times in Washington DC for 1968, 1988, and Twin Cities for 1990 and 2000. The results show that there is a strong dependence on geography: for the Washington DC greater urban area, travel times are relatively stable, while the results for Twin Cities show a marked increase of the commute time. This has also been confirmed in another study [104] that discussed the impact of congestion on mobility patterns. In particular, it has been measured for US cities that typical travel delay due to congestion increases with the population as ∆τ ∼ P 1+δ (14) where δ ≈ 0.3, indicating that average commuting time increases with population, due to the (nonlinear) congestion effects, in sharp contrast with the travel time budget hypothesis (See Fig. 16).

Energy arguments
Mobility, ultimately, is about energy. Indeed, moving an object of a certain mass at a certain speed requires a given amount of energy. Consequently, one would expect to apply energy concepts to understand human travel behavior [108]. First, it has been observed that the average travel times for different modes of transportation are inversely proportional to the energy consumption rates measured for the respective physical activities. Second, when daily travel-time distributions of different transport modes such as walking, cycling, bus, or car travel are appropriately scaled, they turn out to have a universal functional relationship (see Fig. 17).
This corresponds to a canonical-like energy distribution (with exceptions for short trips) hinting at a law of constant average energy consumption related to daily traveling. The argument, first proposed by Kolbl et al. [108], goes as follows. If we define the energy E i spent per transportation mode i, the average energy consumption per day E is constant and independent of the mode of transportation. The corresponding entropy is given by and the constraints on the energy distribution read  which leads to the canonical distribution However, the data indicates that the probability to expend small energy is vanishing, and this was accounted for by introducing a cutoff of the form e −αE/E i . This cutoff term is reminiscent of the so-called "Simonson effect" [109], where short trips are not very likely to be taken with a given mode if the energy spent is much less than αE. The complete energy distribution is then given by This finding, highlights the importance of physical concepts in understanding mobility and more generally social phenomena. In particular, it contains only physical variables such as travel times and energies, that are both measurable, in contrast with utilities introduced in classical choice modeling that describe preferences that are not measurable. Yet, as the authors point out, there are several issues that need to be accounted for, including measurement errors as well as multimodal trips that combine different transportation modes.

Interpolation of scales: the importance of multimodality
As discussed previously, mobility occurs over multiple spatio-temporal scales (multimodal) and thus a comprehensive picture of human mobility requires an accounting of the effects of multimodality. Indeed, the transition between different modes of transport necessarily implies a temporal cost that increases in proportion to the complexity of a trip (number of modes, or indeed level of mode). For example, in the UK, a rough estimate shows that roughly 23% of travel time is accounted for by connections made between modes in trips [102].
As showed in Fig. 18, the average structure of travel time varies with trip length and between cities. With the exception of London, short trips are composed of longer waiting times than riding times (i.e more time is spent being stationary than on the move), where the waiting times are mostly intra-layer in nature, due to bus-bus interchanges. If the transportation network is particularly multi-modal, inter-layer waiting times and walking times start to play a significant role for longer values of .
The effect of multiple scales on movement can be unpacked by defining two types of trajectories. One is the (ideal) quickest path that connects two points and neglects delays due to inter-modal connections (walking and waiting times). For this ideal path to exist, one would need perfect synchronization between modes (bus schedules aligning perfectly with train schedules for example). The other (more realistic) type of trajectories are so-called time-respecting paths that are the quickest paths accounting for inter-modal effects such as arrival, departure and connection times.
For OD pairs i and j, denoting the ideal path time as τ m (i, j) and the corresponding time-respecting path as τ t (i, j), the delay due to the lack of synchronization can be measured by Averaging this ratio over all OD pairs in a city, one obtains a characteristic delay δ.
As an example, for cities in the UK the variation of δ with trip length l is unimodal, reaching its maximum δ max for short trips and then decreasing with distance according to where ν ≈ 0.5. This surprising collapse (see Fig. 19) suggests an underlying process describing the accumulation of waiting and walking times along time-respecting paths. It is clear for a given trip that the number of modes and the frequency of their use play an important role in determining δ. To measure their effect a natural quantity is the average number of stop events per unit time where C α is the number of stop events in mode α and ∆t α the duration spent in that mode.
The quantity Ω can be thought of as a measure of the efficiency of transportation modes in terms of synchronization. The connection betweenδ and Ω appears to follow a power law [102] where µ ≈ 0.3 ± 0.1. The trend is due to the fact that larger values of Ω implies larger frequency and thus a better synchronization between modes. The observed small value of µ is however bad news in terms of efficiency: Decreasing inefficiency by a factor of 2 requires a 10-fold increase in Ω meaning essentially that the trip would need to be unimodal. Thus as the results show, for large cities with multi-modal transportation system, the temporal cost due to interchange between layers plays a key role in the statistics of trips. Indeed, a challenge is posed due to the entanglement between the temporal and multilayer aspects of the system.

General Mobility Models
Models of human mobility can be aimed at reproducing individual mobility patterns or general population flows. In both cases one must necessarily take into account the characteristic spatial and temporal scales of the mobility process, which can vary from hundreds of meters to thousands of kilometers and from hours to years. For this reason, each of these cases have been tackled with distinct modeling frameworks. Individual mobility is subject to a certain level of uncertainty associated with free will and arbitrariness in the actions of individuals, leading to a degree of stochasticity in trip patterns. Consequently, minimal models borrow concepts and methods from random walks and Brownian motion. However, several studies highlighted that individual trajectories are far from random, possessing a high degree of regularity and predictability, which can be exploited to predict an individual's future whereabouts and to construct realistic generative models of individual mobility. At the level of population flows, models describe the aggregate mobility of many individuals and aim to reproduce Origin-Destination (OD) matrices by estimating the average number of travelers between any two spatial regions (e.g. municipalities) per unit time (i.e. daily in the case of commuting flows, yearly in the case of migration flows). Most modeling approaches derive the mobility flows as a function of a few relevant variables of the regions considered, such as mutual distances, population levels, GDP per capita, etc. The next sections review the state-of-the-art mobility models starting from individuallike approaches, moving on to flows at the population level and concluding by the interpolation between these scales (inter-modality).

Individual-Level (Random walks)
A random walk is mathematically defined as a path formed by successive discrete random steps although they can be described in the continuum limit as well (Sec. 4.1.3). The simplest version, however, deals with spatial displacements ∆X i that are taken at discrete times t i . If x 0 = 0 corresponds to the initial position of the walker at time t i = 0 (we can arbitrarily fix the origin in the initial position of the walker), then the position after N steps is given by the random variable where each displacement ∆X i is a random variable extracted from a probability distribution f (∆x), and draws are assumed to be statistically independent. The probability distribution f (∆x) can be used to determine the probability density function (PDF), P(x, t) for the process to be at position x at time t which completely characterizes the nature of the walk and the related spatial and temporal measures. For example, the mean square displacement (Sec. 3.1.2) corresponds to the second moment MSD(t) = X(t) 2 , and the n-th moment is obtained from where brackets indicate ensemble averages over multiple realizations of walks.
Of particular interest when analyzing models of individual mobility is the scaling of the square root of the mean squared displacement (RMSD), R(t) = √ MSD(t), with time t; a relationship that characterizes the speed of displacement from the origin with time. The scaling of the MSD can be used to categorize the type of diffusive motion of the random walker. Ordinary Brownian motion (Sec. 4.1.1) has a MSD that scales linearly with time (for any spatial dimension) therefore on average we expect that after a time t the distance of the walker from the origin is proportional to the square root of the elapsed time; R(t) ∼ t 1/2 . Taking this definition of diffusion, random walks that have displacement growing at a slower rate than t 1/2 are said to be sub-diffusive. In contrast, if displacement grows at a rate faster than t 1/2 , the random walk is classed as super-diffusive. The models of random walks that can lead to these different behaviors are outlined below.
As an example, the most basic form of a random walk is the discrete symmetric random walk in 1D, for which ∆X i = ±1 with equal probability. For this simple case after an elapsed time time t = N (where N is the number of steps taken by the walker), the first two moments are X N = 0 and X 2 N = N.

Brownian Motion
Brownian motion is a class of random walk originally developed to describe the motion of a particle suspended in a fluid (liquid or gas) [110]. Such a particle undergoes many rapid instantaneous collisions with much smaller particles in the medium, resulting in a trajectory characterized by a series of irregular and random displacements. Mathematically, a 1-dimensional Brownian motion is a random walk in the space of real numbers R with independent and normally distributed increments where the probability to observe a displacement of magnitude X from the origin location after a time t is Gaussian distributed with mean zero and variance proportional to t. Brownian motion can be defined as a limit of the discrete symmetric random walk. Let us assume that the particle can take steps of length 1/ √ k to its left or right with equal probability, and that after time t the particle has taken N = t k steps. For a given k, the displacement of the particle at time t is given by and taking the limit k → ∞ leads to Brownian motion as a consequence of the Central Limit Theorem (CLT). In fact, in the large k limit X k (t) tends to X(t), whose PDF is where µ = ∆X and σ 2 = ∆X 2 are the mean and variance of the random walk displacements. For this case, the first moment is zero and σ 2 = 1, so the MSD is simply t and the scaling of the RMSD is R(t) ∼ t 1/2 , corresponding to ordinary diffusion. Using a similar limit, it can be shown that a d-dimensional random walk converges to a d-dimensional Brownian motion, for which the probability that the walker is found at distance x from the initial position is P( 4 t , and the RMSD always scales as the square root of t, R(t) ∼ t 1/2 .

Lévy Flight
Unlike Brownian motion, there is a class of random walks called Lévy Flights, for which one cannot use the CLT. A Lévy Flight is composed of a series of small displacements, interspersed occasionally by a very large displacement. It is formally defined as the sum of independent identically distributed random variables whose PDF for a single jump has a divergent second moment due to a long-tailed distribution of the form with 0 < β < 2. If the displacement after k steps of size 1/k 1/β is defined as X k (t), then the random variable Z N defined as the re-scaled sum of N = t k independent random variables distributed as Eq. (27), takes the form The rescaled variable satisfies a generalization of the CLT, namely the Lévy-Khintchin theorem which states that the PDF of Z N in the limit N → ∞ is a so called α-stable (Lévy) distribution. These distributions do not have a closed formed in real space; instead the characteristic function can be written in Fourier space, and for this particular example, the tail would show the same power law behavior as Eq. (27). Changing variables from Z N to X k via Eq. (28) we obtain The MSD is equal to the second moment and hence RMSD for a Lévy flight scales super-diffusively: R(t) ∼ t 1/β . The jumps in Lévy flights can be considerable, but occurs within the same time step as a short one (which is unrealistic), and thus Lévy flights can be seen as only a rough approximation to actual human trajectories.

CTRW
The random walk models discussed thus far have been discrete in time. In each time interval, an instance of a jump is dictated by the corresponding jump-length distribution. A continuous time random walk (CTRW) is a random walk in which the number of jumps made in a time interval dt is also a random variable or equivalently, the time elapsed between jumps (wait-time ∆T ) is also a random variable.
If the PDF of jump-lengths is f (∆x) and that for wait-times is φ(∆t), and these are independent, the CTRW consists of pairwise random and independent events with ∆X and ∆T drawn from the joint PDF, P(∆x, ∆t) = f (∆x)φ(∆t), which denotes the probability that a jump of length ∆x is taken after a time ∆t. According to this model, after N steps, the total displacement, X N , and the total elapsed time T N are given by: The PDF of the process, P(x, t) can be Fourier-Laplace transformed to give: whereφ(u) andf (k) are the Laplace and Fourier transforms of φ(∆t) and f (∆x) respectively. Taking the inverse transform we get This expression for P(x, t) can be analyzed according to the asymptotic behavior of distributions φ(∆t) and f (∆x). In particular four types of models are obtained, depending on whether none, either, or both distributions have heavy tails. The properties of these models and their application to human mobility are discussed below. Ordinary Diffusion: If the expectation value of φ(∆t) and the variance of f (∆x) are both finite, i.e 1 −φ(u) ∼ τ u andf (k) ∼ 1 − (σ k) 2 then from Eq. (32) and Eq. (33), asymptotically: where D is a diffusion constant. Therefore a CTRW with well defined jump-length and waiting-time distributions, is asymptotically (with time) equivalent to Brownian Motion.
Lévy Flights: If f (∆x) ∼ ∆x −(1+β) (0 < β < 2), and φ(∆t) has a finite variance, a Lévy flight is recovered. Following the same procedure as Sect. 4.1.2, the PDF of the process is given by where G(x/t 1/β ) is a scaling (limiting) function. Therefore if the distribution of jump lengths is a power-law and the wait-times are well defined, R(t) ∼ t 1/β and a continuous time random walker will follow a super-diffusive path, equivalent to a Lévy Flight. Fractional Brownian Motion: Conversely, if jump-lengths with a finite variance are combined with a power-law distribution of wait-times φ(∆t) ∼ ∆t −(1+α) , (0 < α < 2), the PDF is with H being a non-Gaussian limiting function. In this case, the effect of the distribution of waiting times is to slow down the random walk. Here R(t) ∼ t α/2 , consequently the walk is sub-diffusive for α < 1 and super-diffusive for 1 < α < 2.
Ambivalent Processes: The fourth variant occurs when both f (∆x) and φ(∆t) are heavy-tailed. In this case, R(t) ∼ t α/β , and hence the nature of the diffusive behavior is fully specified by α and β. For β < 2 α, the CTRW is super-diffusive and for β > 2 α, it is sub-diffusive. If β = 2 α, the random walk converges to ordinary diffusion/Brownian motion, despite the diverging moments of the respective distributions. A schematic of these limiting cases is provided in Fig. 20. Given the sensitivity of the models to the parameters of the wait times and jump lengths, the importance of their accurate measurement from data cannot be over-stated. In any event, the ambivalent process model is most often used for describing the mobility of individuals. Analysis of various data sources (GPS, CDRs, Dollar bills) has found that both jump length distributions and the distribution of wait times display power-law behavior [46,48,86,81]. The parameter ranges measured from data correspond to α estimated from empirical data range from 0.42 ≤ α ≤ 0.8 [86,81] and values of 0.31 ≤ β ≤ 0.75 [86,48].

Preferential Return
An important aspect of human behavior, missing from the models discussed thus far, is the tendency of individuals to return to one or more locations on a daily basis (socalled preferential return). In a CTRW process, the number of distinct sites S visited by a random walker in time t is given by with µ = α for a CTRW and µ = 1 for a Lévy Flight [111]. The analysis of CDR's (Sec. 2.3), revealed a power-law distribution of jump-size with exponent α = 0.8, while µ was independently measured to be µ = 0.6 ± 0.02, considerably less than the theoretical prediction. Indeed for random walks, the probability of an individual to visit any distinct location becomes asymptotically uniform, whereas analysis of visitation patterns [48] from data suggests that the rank-frequency visitation frequency follows a Zipf's law: where k corresponds to the rank of location according to frequency of visit. This implies the distribution of visitation frequencies follows P( f ) ∼ f −(1+1/ζ) . Furthermore the scaling of the MSD in the CTRW model suggests that an individual will asymptotically drift away from the origin (home) in contrast for what is known from daily experience. To account for the human trait of returning to locations, Song et al. [81] include two extensions to the CTRW model: exploration and preferential return.
Exploration is defined as the probability for an individual to move to a previously unvisited location according to where ρ and γ are parameters of the model and S increases by one for each new visit. Preferential return on the other hand is the probability to return to a previously visited location and is thus complementary to exploration: P ret = 1 − ρ S −γ (see Fig. 21 for a schematic). The probability to visit a previous location i, denoted Π i , is determined by the number of previous visits to i, i.e., Π i = f i where f i is the visitation frequency. The parameters of the model are restricted to 0 < ρ < 1 and γ ≥ 0. Considering the total number of steps, N, an individual makes in time t, and noting that dS /dN = P new , the number of distinct locations visited in time t is given by S ∼ N 1/(1+γ) . For a power-law waiting time distribution φ(∆t), we have t ∼ n 1/α and hence the exponent in Eq. (37)is thus leading to slower rate of exploration than predicted by CTRW processes. Furthermore, the number of visits to location i at step N, m i (N), increases according to where propose that the main reason for this discrepancy is that two generic mechanisms, exploration and preferential return, both unique to human mobility, are missing from the traditional random-walk (Lévy-flight or CTRW) models. (i) Exploration: random-walk models assume that the next diffusive step is independent of the previously visited locations. In contrast, the scaling law (1) indicates that the tendency to explore additional locations decreases with time. Indeed, the longer we observe a person's trajectory, the harder it is to find locations in the vicinity of their home/workplace that they have not yet visited. (ii) Preferential return: in contrast with the random-walk-based models for which the visitation probability is random and uniform in space, humans show significant propensity to return to the locations they visited frequently before, such as their home or workplace.
In what follows we present an individual-mobility model that incorporates ingredients (1) and (2), showing that they are sufficient to explain anomalies A-C. The model, intended to describe the trajectory of an individual, assumes that at time t = 0 the individual is at some preferred location (see Fig. 2). After a waiting time 1t chosen from the P(1t ) distribution, the individual will change his/her location. We assume that the individual has two choices. (i) Exploration: with probability the individual moves to a new location (different from the S locations he/she visited before). The distance 1r that he/she covers during this exploratory jump is chosen from the P(1r) distribution and his/her direction is selected to be random. As the individual an assumption known as preferentia advantage in network and social scienc

Model predictions
The individual-mobility model has tw 0, both of which control the user location during his/her next move ver visited location. The numerical values be determined later from the empirical To explain the origin of the anom individual-mobility model the probabi to a new location is proportional to predicting S ⇠ n 1/(1+ ) , where n is th moves the individual had up to time t . distribution P(1t ) ⇠ |1t | 1 the tim jumps n as t ⇠ n 1/ (Supplementary Sec follows (1) with the exponent To verify the validity of this prediction model, in Fig. 1d we calculated S(t values of , finding that the asympto follows equation (5). Therefore, we with the empirical data.
To account for anomaly B we no visits to location i, increases as dm 5 i = f i = m i /6 i m i (n) is the probabili i during step n. When > 0, in t probability to explore a new location is return visits; thus, asymptotically we As 6 i m i (n) = n, we obtain m i (n) = jump during which location i was fir m i (n i ) = 1. Owing to preferential retur visited, the more it is visited later. Thu i coincides with the order in which Note that for = 0 and ⇢ = 1 eq indicating, as expected, that the visita in the CTRW limit.
To test the validity of prediction (6) different values of and ⇢ in the indivi that the numerically observed scaling with equation (6). The model also sug law is rooted mainly in the preferent calculate f k for a 'democratic' model in . Combined with the fact that f i is proportional to m i (N) this results in the relation f k ∼ k −ζ with ζ = 1 + γ, thereby accounting for the relation in Eq. (38).
The MSD is related to S via where c is a constant. Given the measured range of the parameters, this leads to three asymptotic regimes: MSD ∼ (log t) 2/β for ζ < 1; MSD ∼ (log log t) 2/β for ζ = 1; and MSD → X max for ζ > 1 where X max denotes the saturation point of the MSD. Measured values suggest that ζ = 1.2 ± 0.1 indicating a saturation of movement at long times [81]. This regime corresponds to an individual's motion being dominated by their most visited location and is more in line with expected human behavior.

Recency
The concept of recency was introduced to solve discrepancies that emerges under the standard preferential return assumptions. Specifically, the fact that the earlier a location is discovered, the more visits it will receive, leading to a cumulative advantage precluding people from changing preference of location (unlike what is observed in data). Barbosa et al. [112] proposed a model in which the exploration phase in human movements also considers recently-visited locations and not solely frequently-visited locations. Two rank variables K f and K s are defined to characterize, respectively, the frequency and recency of a given location in the context of an individual's trajectories. More precisely, the rank variables can be described as: K s is the recency-based rank. A location with K s = 1 at time t means that it was the previous visited location. K s = 2 means that such location was the second-mostrecent location visited up to time t and so on.
K f is the frequency-based rank. A location with K f = 1 at time t means that it was the most visited location up to that point in time. Similarly, a location with K f = 2 is the second-most-visited location up to time t, and so on.
The model can be described as follows: first, a population of N agents is initialized and scattered randomly over a discrete lattice with L × L cells, each one representing a possible location. The initial position of each agent is accounted as its first visit. At each time step agents can visit a new location with probability similar to the preferential return model, Eq. (39).
The return phase happens with probability P ret analogous to Eq. (39), with the caveat that return jumps selects location i from frequently visited ones with probability α and recently visited locations with probability 1 − α thus: where k s (i) is the recency-based rank and k f (i) is the frequency-based rank of the location i. When α = 1 the preferential return model is recovered. The empirical measurements for the two datasets considered had α = 0.1, and ν = 1.6. γ = 0.6 was kept the same as [81]. As seen from Fig. 22 the preferential return model does not capture the broader distribution of p(K f , K s ) for recently-visited locations, an effect captured by the recencybased refinements. The primary differences can be seen for the distribution of the recency rank K s seen in log-linear scale.

Social-based models
It is natural to assume that two individuals who have a social interaction, such as friends, family or colleagues, do not always move independently [113,114,115]. Occasionally, they will share full trips, destinations or origins. The trips can be also synchronized if the objective is to meet somewhere or go back home after a meeting. Furthermore and closing the loop, the social network of an individual typically reflects the geography of their life with tighter connections with people spatially closer or at least in clusters related to the places in which the person has previously resided. These correlations have been observed in several publications. For example, a relation between distance and "online friendships" was described in 2005 by Liben-Nowell et al. [116] and was later confirmed using surveys [117,118], social networks [119,120] and (c) Distribution of the frequency ranks, the preferential return mechanisms (labelled EPR) exhibits a powerlaw distribution. The activation of the recency mechanism does not affect the frequency rank distribution. (d) K s distribution, the EPR mechanism does not capture the power-law behavior observed on the empirical data. Figure from [112]. mobile phone records [121,122,123]. Indeed this feature formed the basis of a generative model describing the behavior of population level aggregate economic indicators across urban systems [124]. Possible methods for exploiting this intuitive observation towards improving forecasting of individual movement has been explored in [125]. For example, nonlinear time series prediction methods based on the delay embedding theorem by Takens [126] can be used to forecast an individual's future mobility given a detailed history of past movements and the assumption of a certain degree of determinism in mobility patterns. Furthermore, improvements in prediction accuracy by about one or two orders of magnitude has been observed when including the time series of past movements of members of the social network (individual's acquaintances) Conversely, online [127,128,129] social links can be inferred from co-occurrences of individuals in space and time. In particular, the probability that two users have a friendship link on the Flickr social network grows significantly with the number of distinct geographical locations that they both visited within a given time threshold, where the spatial and temporal information are extracted from the geo-tagged photos they published. The probability of being friends is higher if the size of the regions and the temporal range between the two observations decrease. This result can be explained by simple mobility models exploiting the fact that friends are more likely to spend time together in the same place, often live close to each other and that the jump-size distribution decays as a power law. Similarly, the location of a person may be predicted by those of their near contacts [130]. Two individuals with similar mobility patterns are generally in close proximity in the who-calls-whom social network of mobile phone users [131]. Indeed, strong correlations exist between various classical topological measures of the proximity of two mobile phone users in the social network, such as Adamic-Adar, Jaccard and Katz, and various measures of spatio-temporal proximity, such as co-location rates and spatial cosine similarity. Combining co-location information with information on the relative positions of users in the social network it is possible to predict the formation of new social ties with higher precision than using only information derived from the proximity in the social network.
The idea of interconnecting mobility and social interactions has been considered in different theoretical contexts. Detailed models have been proposed in the area of transportation to take into account the relation between social network and transport demand. This includes microsimulation of transport systems [113,114,115] in which agents may have communication between them and have common objectives [132]. Social groups and their sizes can impact transport demand [133], daily schedules, social relations and face-to-face interactions [134,135,136,137,138]. A simpler framework has been introduced from a Physics perspective, as for example formulated by [139]. The system is composed of N agents within a radius r moving in a 2D square of lateral size L and periodic boundary conditions. The starting position of the agents and their directions of movement are randomly selected. All the agents have the same initial speed v 0 and they travel in straight lines until they "crash". Every crash produces a new social tie between the pair of agents involved. After every crash, the velocities of the two agents are updated randomly selecting new directions of movement and increasing their speeds according to the expression: where v i (t) is the velocity of agent i at time t immediately after the crash, c is a constant and k i (t) is the number of connections of i in the social network. The number of agents N is constant, but after a certain period of time old agents are removed and substituted with new agents with speed v o and no social links, allowing the system to reach a stationary state. The characteristics of the emerging social network are analyzed. The degree distribution depends on the constant c and on the lifetime of the agents in the system. For small c, the distribution is well approximated by a Poissonian distribution, but deviates from this shape with increasing c and average degree k . Beyond the distribution, the model can reproduce other properties of real social networks such as the degree-degree correlations, the size of the large connected cluster and the abundance of cliques. More realistic mobility models in the context of generated social networks were explored in [140,141]. In particular the former proposed a model inspired by the methods described in Sec. 4.1.4 with important variations in the main ingredients (See Fig. 23). The initial conditions were set such that the agents are placed according to measured population density in the area considered. Then at every step, each agent decides to visit a social contact, with probability p v or, otherwise moves according to a Lévy-like flight with probability 1 − p v , with the destination chosen in proportion to the population density in that area. This precludes agents from moving to unphysical locations such as water-bodies, geographic barriers and the like. After movement, the agent generates a new directed social link with probability p in the neighborhood of her present position and with probability p c with a random agent. The parameter p c represents social relations that emerge online, inspired by the observations made by [116]. Note that the jump distribution can be calibrated from measurements made in [81], and the tie formation parameter p can be sampled at different scales from measurements made in [128]. This leaves only two free parameters in the model: the probability of visiting friends p v and that of creating new links regardless of the distance p c . These parameters were adjusted to reproduce an error function containing topological and geographical properties of a set of online social networks such as Twitter, Gowalla and Brightkite. With these parameters, the model reproduces the degree distributions found in [139] as well as the spatial dependence of link distribution, the probability of reciprocal ties, the overlap of the social environment, the density of triangles and the disparity of the triangles in terms of distances of their vertices. The model is also amenable to a mean-field analytical treatment that demonstrates the necessity of accounting for social effects on mobility (visiting friends for example) to obtain realistic reproductions of the behavior seen in social networks. The variant proposed by [141] instead focuses on behavior seen at shorter timescales such as intra-day mobility, for which period the social network can be approximated as static. Like the preferential return model, at each step an agent decides to return to a previously visited location with probability 1 − ρ S γ or explore a new one with ρ S γ , but there is an additional "social pressure" component. This occurs with probability α where one of the agent's contacts are selected at random with probability proportional to the co-similarity in the location visiting profiles. Consequently, the new location is chosen from the list of the contact's previously visited locations. In the absence of social pressure, (probability 1 − α) the behavior is the same as Sec. 4.1.4. The introduction of the social pressure component does a better job at reproducing the re-visitation profiles measured in CDR data than models without this component.

Population-Level
Mobility information gathered at the individual level can be aggregated to study the flows of individuals traveling from one region to another at different spatio-temporal scales. These flows can be organized in the framework of Origin-Destination (OD) matrices (see Sec. 3.1.5). Such a format, with all possible combinations of origins and destinations for trips, is easily transformed into a directed weighted network in which nodes denote locations (for example counties or municipalities) and link weights correspond to the flow of travelers between the two locations. As discussed earlier in the review, OD matrices can be empirically estimated from transportation surveys, traffic counts or individuals' geolocated ICT data. An OD matrix, or its corresponding network, provides useful information on the travel demand between the origin and destination areas, representing a valuable asset widely studied and used in Geography, Transportation research and Urban Planning. Therefore, being able to obtain accurate estimation of OD matrices is crucial for both modeling and applications, and this problem has attracted the interest of researchers and decision makers for decades.
Note that individual mobility patterns are straightforwardly aggregated into flows, however the inverse problem, i.e the disaggregation of flows, is typically not possible. Instead, one has to resort to determine dependences between mobility flows and a limited set of "static" attributes of the locations that would allow predictions on how changes in these attributes can potentially influence future travel demand. To this end, various spatial interaction models have been proposed to predict flows of individuals based on a small number of key local attributes. Considering a region of interest divided into n locations, the purpose of these models is to estimate the number of trips T i j from location i to location j from the socio-economic characteristics of the populations of i and j, and their spatial distribution. Models of spatial flows have been traditionally developed starting from the principle of entropy maximization subjected to various constraints. In the strongest version, the constraints are the number of people leaving and entering each location, but in softer ones there can be other proxy variables to represent the demand and attraction of the trips' origins and destinations such as population levels, the finite amount of resources for travels, utility functions to describe individual choices over competing alternatives, etc [142,143,144,145,146].
In the simplest version, the objective is to infer flows from the product of two types of variables: one type that depends on an attribute of each single location (e.g. the population), and the other type that depends on a quantity relating a pair of locations (e.g. the distance or travel time). The differences between the various models consist of the choice of variables considered, and the specific functional forms in which these variables enter. Over the years, two main schools of thought have emerged. The first one assumes that the number of trips between two locations is a decreasing function of their distance, giving rise to the so-called gravity models [148,6]; the second variant postulates that the number of intervening opportunities, defined as the number of potential destinations between two locations, determines the mobility flow between them, and models that share this assumption are called Intervening Opportunities models [5]. An illustrative example that highlights the fundamental differences between distance and intervening opportunities is shown in Fig. 24. In addition to providing the mathematical framework to model human mobility flows [39,149,150], these models have found successful applications in estimating other spatial flows, including cargo shipping volume [151] and social interactions from inter-city phone calls [122,152].

Gravity models
Preceded conceptually by the work of H.C. Carey in 1858 about land use [148] and the retailing models of W. J. Reilly in 1931 [153], George K. Zipf proposed in 1946 an equation to calculate mobility flows inspired by Newton's law of gravitation [6] 12 . In his work, Zipf highlights the importance of the distance for human migration patterns where the magnitude T i j of a migratory flow between two communities i and j can be approximated by where P i and P j are the respective populations and r i j the distance between i and j.
The basic assumptions of this model are that the number of trips leaving i is proportional to its population, the attractivity of j is also proportional to P j , and finally, that there is a cost effect in terms of distance traveled. These ideas can be generalized assuming a relation of the type: where K is a constant, the masses m i and m j relate to the number of trips leaving i or the ones attracted by j, and f (r i j ), called a "friction factor" or "deterrence function", is a decreasing function of distance. As in the original version of the model, the masses usually are functions of population (not necessarily linear); common functional forms used in the literature are m i = P i α or m j = P j γ [89]. Unlike the original version, however, other variables, such as gdp-per-capita, may factor in the definition of the masses [154,155]. The distance function f (r i j ) is commonly modeled with a powerlaw or an exponential form, although more complicated functions can be considered, such as a combination of the two, Indeed, the optimal form of the function may change according to the purpose of the trips, the spatial granularity of the locations, and the transportation mode [39]. For example, in the case of commuting flows, the value of the exponent is highly correlated with the scale [156,157] as can be seen in Fig. 25. In other applications, the distance may not be the appropriate variable to quantify the cost of travel between two locations, and other variables such as travel time or economic (i.e. monetary) cost of a trip may offer better characterizations. The ability to estimate, even if as a crude approximation, trip-flows, and consequently, traffic demand between two different locations as a function of their local properties, has made the gravity model widely popular in transport planning [158,89], in studies of geography [159] and spatial economics [160,161]. It has been also used in situations where the knowledge of mobility flows is essential but the empirical data is not available as in the case of the characterization and modeling of epidemic spreading patterns [162,163,164,165,155]. Despite its widespread use and historical popularity, one must keep in mind that the gravity model is a gross simplification of travel flows, and in many cases, falls far short of capturing actual empirical observations [42,166,157]. Furthermore, the model requires the estimation of a number of free parameters, making it rather sensitive to fluctuations or incompleteness in data [150,42].
Constrained gravity models. Some of the limitations apparent in the gravity model, can be resolved via certain constrained versions. For example, one may hold the number of people originating from a location i to be a known quantity O i , and the gravity model is then used to estimate the destination, constituting a so-called singlyconstrained gravity model of the form, As can be seen, in this formulation, the proportionality constants K i depend on the location of the origin and its distance to the other places considered. One can go further and fix now also the total number of travelers arriving at a destination j as D j = i T i j , leading to a doubly-constrained gravity model. For each Origin-Destination pair, the flow is calculated as where there are now two flavors of proportionality constants, usually calibrated with an Iterative Proportional Fitting procedure [167]. The use of singly-, doubly-or non-constrained models depends on the amount of information available and on the pursued objective. If the aim is to approximate the mobility flows and transport demand from indirect socio-economic variables of different geographical areas, then one employs non-constrained models. On the other hand, if out-going or in-going flows are empirically measured quantities, and the objective is to estimate the elements of the OD matrix T i j , then one employs constrained models.
Maximum entropy derivation of gravity models. While the gravity model may seem rather ad hoc in terms of the form of Eq. (46), an argument in favor of this functional form was provided by Alan Wilson [168] in the framework of classical transportation theory combined with entropy maximization. Essentially, the argument posits that in the absence of detailed information, out of all possible variants of OD matrices, the most probable ones are those that can be obtained with the highest number of trip configurations under opportune global and local constraint satisfaction. The objective is thus to find the set of flows {T i j } that maximizes the number of possible configurations of trips associated with it, respecting the possible constraints.
Let Ω({T i j }) be the number of distinct arrangements of individuals (configurations) that give rise to the set of flows {T i j }, corresponding to the number of ways in which T 11 individuals can be selected from the total number of travelers T = i j T i j ; T 12 from the remaining T − T 11 and so on and so forth. Then we have that One then determines the maximum of Ω using Lagrange multipliers, subject to the constraints: where C i j is the cost of travel from location i to location j. The first two constraints ensure that trips originating and terminating in every location are equal to their observed values, whereas the third constraint fixes a total cost for all trips, C. In the limit of a large number of trips, T , the configuration that maximizes Ω is by far the most probable, and hence dominates over all other configurations. The resulting doublyconstrained gravity model takes the form where the values of K i and L j are set in order to fulfill the first two constraints in Eq. (52). The Lagrange multiplier β appearing in Eq. (53) controls the effect of cost on flows, and its value is empirically calibrated. It is worth noting that one can introduce a power law distance decay f (r i j ) = r −β i j by considering a cost function of the form C i j ∝ ln(r i j ). These arguments have been further developed in a series of references [145,146] that also includes the distinguishable nature of the travelers and therefore necessarily modifying the statistics.
Using gravity models to estimate mobility flows. In order to calibrate and apply gravity models to estimate mobility flows within a given region the following procedure is usually adopted. The best fit values of the parameters are determined using an optimization algorithm that either minimizes some error function between the model's estimates and the observed data [169], or maximizes the likelihood function of the observed data given the model's parameters [170]. Generalized Linear Models (GLM) [171] are a generalization of linear regression that are usually applied to fit the parameters of globally and singly constrained gravity models. GLM methods are more adapt than Ordinary Linear Regression (OLM) as it allows for the use of a wider and more realistic range of probabilistic models to capture fluctuations in flow estimates.
Both the model training and parameter calibration steps can be employed either on the entire dataset or on subsamples, particularly if the aim is to validate the performance of the model in relation to the remaining set.

Intervening opportunities models
Along with the gravity model, one of the first attempts to provide a conceptual and formal model of human mobility was introduced in 1940 by Stouffer [5]. Contradicting a long tradition in the social sciences -that distance is the central factor in determining mobility -Stouffer proposed a conceptual framework in which distance and mobility are not directly related. Instead, Stouffer suggested that what plays the key role in determining migration is the number of intervening opportunities or the cumulative number of opportunities between the origin and the destination. In the paper, Stouffer does not provide a precise definition for "opportunities", leaving it to be defined depending on the social phenomena under investigation.
The basic idea behind the intervening opportunities (IO) model is that the decision to make a trip is not explicitly related to the distance between origin and destination, but to the relative accessibility of opportunities for satisfying the objective of the trip. An opportunity is a destination that a trip-maker considers as a possible termination point for their journey, and an intervening opportunity is a location that is closer to the trip maker than the final destination but is rejected by the trip maker. The law of intervening opportunities as proposed by Stouffer in 1940 states "The number of persons going a given distance is directly proportional to the number of opportunities at that distance and inversely proportional to the number of intervening opportunities". Stouffer used this model to estimate migration patterns between services and residences [5]. The traditional form of the intervening opportunities model is usually given by Schneider's version of Stouffer's original model [172]. Schneider's hypothesis states that "The probability that a trip ends in a given location is equal to the probability that this location offers an acceptable opportunity times the probability that an acceptable opportunity in another location closer to the origin of the trips has not been chosen". More formally, the flow T i j from the origin location i to the j-th location ranked by travel cost from i is given by: Here, O i is the total number of trips originating from i and the second term represents the probability that one of these trips ends in location j. The denominator is a normalization factor ensuring that the probabilities sum to 1 (i.e. j T i j = O i ). This probability depends on V i j , the cumulative number of opportunities up to the j-th location ranked by travel cost from the origin location i (and n is the total number of locations in the region considered). Usually, the population, m j , or the total number of arrivals, D j = i T i j , are assumed to be proportional to the number of "real opportunities" in location j. The value of the parameter L can be seen as the constant probability of accepting an opportunity destination. As in the case of the gravity model, the value of L is adjusted in order to obtain simulated flows as close as possible to observed data. Several studies and variants of intervening opportunities models have been developed on this concept [173,174,175,159,176,177]. The gravity and the intervening opportunities models have been compared several times during the second half of the twentieth century, showing that generally both models perform comparably [178,179,180,181]. In fact one can cast the intervening opportunities model as a special variant of the gravity model with the friction factor f (r i j ) replaced by a function of the number of opportunities between the two locations, f (V i j ) = e −L V i j−1 [182,181]. Consequently, the intervening opportunities model can also be derived from entropy maximization methods as discussed earlier in the context of the constrained gravity model. Indeed, some authors have proposed hybrid gravity-opportunities models taking into account both the effect of distance and the number of opportunities between locations [183,184].
However, unlike the gravity approach and despite relative good performances, the intervening opportunities model has suffered a growing loss of popularity. This is mainly due to the lack of research effort into the implementation and calibration of the model, attributable to the fact that the theoretical and practical advantages of opportunities models over the gravity ones are not overwhelming [89].

The radiation model
The concept of intervening opportunities and Schneider's hypothesis has recently triggered a renewed interest thanks to the recently proposed radiation model [42]. The radiation model assumes that the choice of a traveler's destination consists of two steps. First, each opportunity in every location is assigned a fitness represented by a number, ersal model for mobility and migration ns 2,3 , Marta C. González 4 , Amos Maritan 2 & Albert-László Barabási 1,5,6 contemporary form in 1946 (ref. 1), but with roots he eighteenth century 2 , the gravity law 1,3,4 is the prerk with which to predict population movement 3,5,6 , lume 7 and inter-city phone calls 8,9 , as well as bilateral een nations 10 . Despite its widespread use, it relies on eters that vary from region to region and suffers lytic inconsistencies. Here we introduce a stochastic g local mobility decisions that helps us analytically ng and mobility fluxes that require as input only he population distribution. The resulting radiation obility patterns in good agreement with mobility tterns observed in a wide range of phenomena, from tion patterns to communication volume between . Given its parameter-free nature, the model can be where we lack previous mobility measurements, proving the predictive accuracy of most of the cted by mobility and transport processes [11][12][13][14][15][16][17][18][19][20][21][22][23] . th Newton's law of gravity, the gravity law assumes of individuals T ij that move between locations i and j roportional to some power of the population of the estination (n j ) locations, and decays with the distance as re adjustable exponents and the deterrence function fit the empirical data. Occasionally T ij is interpreted as te of individuals travelling from i to j, or an effective the two locations 24 . Despite its widespread use, the otable limitations: e, we lack a rigorous derivation of (1). Whereas zation 25 leads to (1) with a 5 b 5 1 , it fails to offer rm of f(r). o, lacking theoretical guidance, practitioners use a ce functions (power law or exponential) and up to to fit the empirical data 5,7,8,11,14 . ee, as (1) requires previous traffic data to fit the para-, it is unable to predict mobility in regions where we raffic data, areas of major interest in modelling of s. r, the gravity law has systematic predictive discrep- Fig. 1a we highlight two pairs of counties with similar ination populations and comparable distance, so the flux between them should be the same. Yet, the upplementary Information) documents an order of rence between the two fluxes: only 6 individuals n the two Alabama counties, whereas 44 do in Utah. , equation (1) predicts that the number of commuters t limit as we increase the destination population n j , yet  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 model (b, c). b, The definition of the radiation model: an individual (for example, living in Saratoga County, New York) applies for jobs in all counties and collects potential employment offers. The number of job opportunities in each county (j) is n j /n jobs , chosen to be proportional to the resident population n j . Each offer's attractiveness (benefit) is represented by a random variable with distribution p(z), the numbers placed in each county representing the best offer among the n j /n jobs trials in that area. Each county is marked in green (red) if its best offer is better (lower) than the best offer in the home county (here z 5 10). c, An individual accepts the closest job that offers better benefits than his home county. In the shown configuration the individual will commute to Oneida County, New York, the closest county whose benefit z 5 13 exceeds the home county benefit z 5 10. This process is repeated for each potential commuter, choosing new benefit variables z in each case. Each offers attractiveness (benefit) is represented by a random variable with distribution p(z), the numbers placed in each county representing the best offer among the jobs in that area. Each county is marked in green (red) if its best offer is better (lower) than the best offer in the home county. (c) An individual accepts the closest job that offers better benefits than his home county. Figure from [42]. z, chosen from some distribution p(z), whose value represents the quality of the opportunity for the traveler. Second, the traveler ranks all opportunities according to their distances from the origin location and chooses the closest opportunity with a fitness higher than the traveler's fitness threshold, which is another random number extracted from the fitness distribution p(z) (see Fig. 26). As a result, the average number of travelers from location i to location j, T i j , takes the form: Here again, the destination of the O i trips originating in i is sampled from a distribution of probabilities that a trip originating in i ends in location j. This probability depends on the number of opportunities at the origin m i , at the destination m j and the number of opportunities s i j in a circle of radius r i j centered in i (excluding the source and des- tination). This conditional probability needs to be normalized so that the probability that a trip originating in the region of interest ends in this region is equal to 1. In case of a finite system it is possible to show that this is equal to 1 − m i M where M = i m i is the total number of opportunities [166]. In the original version of the radiation model, the number of opportunities is approximated by the population, but the total inflows D j to each destination can also be used [156,166,157]. The great advantage of the radiation model compared with other spatial interaction models is the absence of a parameter to calibrate with observed data. In particular, the flows defined in Eq. (55) are independent of the fitness distribution p(z). However, this advantage represents also a limitation since the model does not seem to be very robust to changes in the spatial scale [156,166,43,157]. To overcome this drawback, a radiation model with opportunities' selection [185] and an extended radiation model [186] have been proposed. In this extended version, the conditional probability to perform a trip between two locations according to the spatial distribution of opportunities is derived under the survival analysis framework introducing a parameter α to control the effect of the number of opportunities between the source and the destination on the location choice. The addition of this scaling parameter seems to greatly improve the performance of the model, and similar to the gravity model, its value can be inferred from the scale of the study region according to the homogeneity of the opportunities' spatial distribution [186,157].
An interpretation of the radiation model in the context of job search and, consequently, for the formation of commuting flow networks has been given in [187]. The basic components are individuals residing within a demarcated geographical area, who are seeking employment. Job locations are uniformly distributed in space and characterized by a fitness parameter z, which itself is drawn from a distribution that includes aspects such as wage and worker-convenience. Individuals have a certain tolerance level, z * , and will search in a progressively increasing radius from their residence, with the search terminating at the first instance of the condition z > z * . Using extreme value statistics, it is possible to determine the distance distribution, P(r), between residence and job places, which is independent of the distribution of z (for non-pathological cases), and was calculated to be of the form where ρ is the density of jobs. The model satisfactorily reproduces some of the empirical features of the travel distance distributions in commuting data as can be seen in Fig. 27.

Comparison between models
The models described so far provide theoretical expectations for the OD matrices in terms of the flow values T i j between geographical regions i and j. Before employing these for practical applications, it is essential to validate the expected (calculated) flows against empirical evidence, that may be limited in the sense of being restricted in space compared to regions of interest, or indeed to shorter time windows than the temporal period of interest. By evaluating the performance of each model, one selects the one that most closely matches the empirical data (if available). The accuracy, in general, will depend on the spatio-temporal scale as well as on available (meta) information and the extent of missing data. It is crucial that comparison between models is done on an equal footing, i.e using the same extent of available data and in equal spatiotemporal resolutions. As discussed, the model output is the value of T i j , subject to constraints, which range from Eq. (46), where only the total sum of the flows is imposed via the parameter K, the singly constrained framework models based on intervening opportunities or the closest opportunity case, which is production constrained, in the sense that outflows O i are fixed a priori and must be preserved for every location, and finally to the doubly-constrained models, where both the outflows O i and the inflows D i in every region i are preserved.
In general, all of the described models share some common assumptions: 1. The trip distribution P i j that generate the flows T i j is independent of the trip production O i . 2. In both the unconstrained and singly-constrained models, choice of travel destinations are statistically independent (i.e there are no memory effects). 3. Flows are estimated as a product of variables related to opportunities and distance (i.e. variables are "separable").
Under these assumptions, it is possible to draw statistical laws that govern the distribution of travels P i j , and from thereon, build models with varying levels of constraints. This is important to adapt the models to the available information and to fairly compare between them. One such comparison was carried out in [157], where commuting data was obtained from the census offices of England and Wales (E&A), France (FRA), Italy (ITA), Mexico (MEX), Spain (SPA), USA and then in shorter spatial scales, from the cities of London (LON) and Paris (PAR). Several variants of the gravity, intervening opportunities and radiation models were considered with multiple levels of constraints. The predicted flows T i j were compared with the empirical data using two related metrics: The common part of commuters (CPC) for all location pairs i, j with positive flows both in the empirical data set and in the model prediction, defined as Here T = i, j T e i j is the total number of commuters, T m i j the model prediction for the flow, and T e i j the empirical value. The CPC is one, if the agreement is perfect and zero if there is no overlap between data and model. The ability of the models to recover the topological structure of the original network was assessed through the second metric, termed the common part of links (CPL) and defined as where 1 T i, j >0 is the Heaviside function. The CPL measures the proportion of links in common between the simulated and the observed networks, it is null if there is no link in common and one if both networks are topologically equivalent. The results in terms of the CPC are shown in Fig. 28. For this choice of datasets, the gravity model is seen to moderately outperform the other models. One must note of course, that this might change across different datasets.

Intermodality
A complete description of mobility must take into account the mutlimodal structure of transport and the transitions between the different modes as alluded to in Sec. 3.3. A good framework for carrying out such an analysis has been developed in recent years in the context of multilayer or multiplex networks [188,189]. These are networks in which the nodes can be present in one or multiple layers, and each layer by itself contains a set of links or node-node interactions; if a large fraction of nodes in the network are present in all apparent layers, then it is termed multiplex. On the other hand if one finds different flavors of nodes restricted within their own putative layers, then the network is dubbed a multilayered one. Indeed, each layer of the network constitutes a separate network in its own right depending on the context of its interactions. For example, in social networks, layers may correspond to different groupings such as coworkers, family, different cohorts of friends and so on, all of which have withinlayer interactions. Information diffusion in the network, however, occurs through each layer and across the full multilayer structure. Translating this to the case of mobility and transportation networks, layers may correspond to different transportation modes (road, subway, airline), while connections between layers constitute the interchanges between these modes. Within this context, the standard way to construct multilayer networks is to associate locations to nodes and flows, and frequency of travel (or simply unweighted simple connections) to links between different transportation modes (see, for example Fig. 29).
Traditionally, single-mode transport networks have been studied using mono-layers. Examples include metro networks in Beijing [191], Boston [192], London [193,194], Montreal [195], Paris [193] and Seoul [196]; bus networks in China [197,198] and Poland [199]; train networks in Boston, Vienna [200] and in India [201]; and combined public transportation networks in 14 cities around the world [202]. In most of these cases, the transport network was generated using different perspectives. One such perspective is the L-space representation, where nodes correspond to stations in the transportation map, and are linked together if they are consecutive stops in a transportation line. Similarly, the network can be built in the P-space framework, if nodes are stations linked together if they correspond to stops in the same line (irrespective of sequence). This latter representation corresponds to a one-mode projection of a bipartite network [69] consisting of stations and lines. (For a graphic explanation of this notation see [202].) In addition to geographical aspects, the link weights can also Fig. 29: An example of a multilayered network approach to mobility. Public transport networks in the London metro area, separated into multiple layers consisting of the bus, subway and rail networks. In this particular instance, one can see that the London bus network is the most "used" layer. Figure from [190].
incorporate information on the number of passengers in the same route, among other features.
The fact that a single transport network can be represented in multiple ways, naturally paved the way for the use of multilayer structures [203,204]. In this case, each layer contains the expression of the network in a different space, and the multilayer allows for a deeper study of the network accessibility for travelers. As mentioned before, a multilayer network can be also employed to condense the transport information coming from different modes. An example for this can be seen in Fig. 29 which shows elements of the transportation network in London, combining the bus, train and metro networks, whereas additional elements including air and ferry can also be incorporated [190]. The more intricate topology of multilayer networks necessarily affects the spreading processes taking place on the system (diffusion and navigation) [205,206,189] that are modeled using the methods discussed in Sec. 4.1 and 4.2. This leads to new features such as congestion [207], change in the navigability of the network subject to disruptions [208,209], produces new dynamic regimes that are absent in monolayers [210] and abrupt reactions to minimal topological changes [211].
A particular type of multilayer representation of transport networks is worth a separate discussion. In this case, nodes still correspond to locations, and links are assigned according to one of the space-frameworks (L or P), but the weight of the links now account for the temporal duration of trips between the nodes. This network configuration, which some authors refer to as a PT -space, appear in each of the layers, where each layer usually corresponds to a single transport mode. With such a representation, the optimal time for a traveler traversing between a given OD pair, can be calculated using optimal path algorithms across the multilayer structure, containing the full transport network and taking into account intermodality. These types of multilayers have been studied, for example, in British [102] and French [212] cities. In the British case, one detects a shift in preference for different transportation modes as well as the overlap between temporally and spatially optimal paths as a function of distance. In the French case, a few particular (temporally optimal) connections were singled out, so-called efficient connections, and their occurrence (and hence influence) in the general mobility modes were studied across a number of cities. Using a similar framework, with two layers corresponding to road and metro networks in the cities of London and New York, congestion effects as a function of varying the speed of metro lines was studied [213]. Finally, it is worth mentioning the point raised in [214], where the authors discuss the differences between building the network by considering a layer as a single transport mode or as a single line. The latter description provide a deeper view of the inefficiencies due to line transfer of the passengers, while the former is well-suited to study topological aspects and geographical coverage.

Selected Domains of Application
Here we discuss some selected applications of the frameworks, concepts, models and datasets introduced thus far. We organize this section by scale, ranging from single scale (mostly dealing with transportation modes), multi-scale applications with a focus on movement in cities and epidemic spreading, and finally end with some new developments related to virtual scales, i.e mobility patterns seen in online activity.

Pedestrian movement
Large concentrations of people such as those occurring during festivals, sport events or religious ceremonies can present a challenge in terms of service and transport demand as well as public safety. Avalanches in stadiums, processions and even subway stations have occurred in recent decades leading in high casualties [215]. In the face of this, strategies need to be devised for ensuring prompt evacuation in cases of emergency when considering the configuration of new buildings, public spaces and transport systems. Consequently, understanding the empirical patterns governing pedestrian movement and being able to build corresponding models, has the potential to benefit a wide range of disciplines including architecture, urban planning and public service stakeholders among others. Given the broad range of interest, a large amount of research has been devoted to this sphere, over the past few decades. Indeed, several reviews exist on this topic [215,216,217,218,219,220], so instead of providing a detailed retrospect, we provide here an overview of the most recent developments.
Modeling frameworks on pedestrian movement can be broadly classified according to the approach used to consider individual pedestrians. For example, when the aim is crowd control, then usually the number of individuals considered is quite large, and the metrics of interest are crowd velocities or pressures. In such a case, one can define fields at every point in space and time to characterize the state of the system. These fields usually include the density of pedestrians ρ( x, t), the local velocity v( x, t) or the flow across an exit, J( x, t). In analogy with fluid dynamics, one sees a continuity equation of the form, stemming from the conservation number of people in the crowd. The motion of pedestrians can be assumed to minimize the time spent in the system. This condition is included via the introduction of a field φ( x, t) coupled with the density in such a way that the local velocity shows a relation where the gradient of φ determines the direction of motion and f (ρ) is a simple function of the density ensuring that the pedestrians speed decreases as the density increases. Possible choices are a linear relation, f (ρ) = ρ − ρ max , keeping ρ < ρ max or a non linear one, f (ρ) = (ρ − ρ max ) 2 , among others [221,222,223]. The connection with φ is established through the function f (ρ) thus, leading to the so-called Hughes model. In a scenario with individuals exiting a room, it produces two regimes in the pedestrian mobility depending on the global density, one with free flow and another with congestion. Furthermore, it turns out that shock-waves and moving congestion fronts emerge due to the narrowing of the crowd close to the door, leading to intermittent exit flow, avalanches and stop-and-go effects [222]. An alternative approach is to use agent based modeling where the state of every individual (agent) i is simulated in detail. This framework is especially applicable to simulations of heterogeneous populations in which the agents may show different features. The most popular model falling into this category is the so-called social force model proposed by Helbing and Molnár [225,226]. The idea behind this model is to apply the analogues of Newton's laws of motion to each pedestrian i. If the agent's velocity is v i , the location r i at time t, then the acceleration can be expressed as The first term on the right hand side accounts for the tendency of the individuals to walk at a certain desired speed v 0 i . It is usually written as which ensures the recovery of the desired speed in a characteristic reaction time τ. The second term f a i is an attraction force depending on the agent's position introduced to guarantee that walkers are compelled to a certain target or direction. The third force, f rb i , is repulsive depending on the agent's location and takes into account the effects of the boundaries and other obstacles on the agent's trajectory. The next term represents The social force model reproduces several phenomena observed in empirical circumstances. For example, when two groups of pedestrians moving in opposite directions meet in a corridor, unidirectional stripes are formed, or when they meet in a door, the flow displays intermittence due to clogging [226] (Fig. 30). Variations of this model have been used to analyze escape patterns from a building in low visibility conditions [224,215], to assess evacuation plans in skyscrapers [228,229], ships [230] and aircraft [231]. In relation to crowd control, they have also been used to analyze how to mitigate avalanches in the Hajj religious gathering in Saudi Arabia [227,232], in the Notting Hill Carnival in London [233] or, along with mobility models for other transport media, to study evacuation of cities and populated areas in case of natural disasters such as flooding [234,235,236]. Additionally, these models provide the basic tool to model scenarios such as the formation of mosh and circle pits as a mark of collective motion in heavy music concerts [237], as well as the risk of disease propagation in massive gatherings [238].
In the last decade, there has been a large volume of research related to the constituent terms of Eq. (62), including empirical calibration of the different coefficients as well as the introduction of a cognitive basis for the different forces [217,232]. In high density regimes, such as those occurring in videos taken of the disaster of the Jamarat Bridge in the Hajj pilgrimage of 2006, the trajectories of single pedestrians become turbulent [227,232] (Fig. 31). The original version of the social force model is unable to reproduce this turbulent regime and consequently modifications to the force terms in high density regimes, were introduced. One of these modifications is the centrifugal force model, so-called because of its similarity with the classical centrifugal force in rotating bodies, and involves a modification of the mutual repulsion term f r i j thus, Here m i is the mass of pedestrian i, K i j is a coefficient enforcing the influence only of pedestrians in front of i on its movements and e i j is the unit vector pointing from i to j [239,240] . This modification has the effect of introducing both laminar as well as turbulent regimes, a feature missing in earlier versions of the social force model. Further refinements also included the interactions of individuals based on specific boundaries and shapes of obstacles [241].
A second modification to the social force model deals with the cognitive capacities and limitations of pedestrians and how they can affect the force terms of the model [242,243]. Examples include analysis of the time scales regulating pedestrian interactions [244], the addition of a field to take into account the (partial) global perception of individuals and the environment under consideration [245,246], or the role of leaders [247]. The capacity of the pedestrians to foresee the movements of others and the dangers of obstacles have been modeled within this framework in [242,248]. The model has also been adapted to agents with an anisotropic perception of their environment [249] and its instabilities in a context of stop-and-go oscillations have been studied in [250,251]. A systematic stability analysis of its numerical solutions has been provided in [252] together with a version of the model that notably increases its computational performance.
In addition to this, there exist a family of models that lie between the macroscopic and microscopic approaches, characterized by discretization of space into cells, and the specification of rules on how agents navigate between these cells. Such cellular automata models have the advantage of computational scaleability, as for example, in [253] where simulations of crowds of size 10 5 , were conducted with running time scaling linearly with the number of individuals. The price to pay for this increased computational performance is the loss of local details within each cell. These models have been used to study clogging in exits, and due to obstacles by adding a "friction term" between agents [231,254], pedestrian flow through subsequent bottlenecks [255] and the formation of stripes in junctions with agents moving in different directions [256]. Some modifications to these models include, for example, behavioral factors of the agents such as path selection based on the time necessary to reach the destination [257], heterogeneous agents [258] and long-range spatial awareness [259].
While some of the models described thus far have been inspired by empirical data obtained from crowds in different contexts, challenges remain in actually monitoring crowds to get accurate and precise measurements on individual movement (see [219] for a recent review on methods used for this). One way to calibrate models is to take a controlled approach, i.e conduct experiments on a smaller scale, in which groups of people are directed to move under specific constraints. This controlled environment can provide insight into questions such as how pedestrians walk in social groups [263], the quantification of pressure and density in different types of bottlenecks and junctions [264], the perception of personal space [265], the clogging effects and instabilities in narrow areas [266,267] and the asymmetric flow produced by stairs [268]. One factor that is common to such experiment-based research is the generation and analysis of "fundamental diagrams" [269] (Fig. 32), that represent the relationship between the average velocity v and the density ρ of pedestrians, as well the flux through a surface J (an exit) among other parameters. In the free flow regime, with no obstacles and low density, pedestrians tend to move at an optimal velocity. As the density increases, the velocity decreases until the crowd stops. Similarly, the flux through an exit also decreases with increasing pedestrian concentrations. The relations between these  [260].On the right, experiments in Germany with two set of people: students and soldiers, figure originally taken from [261]. The composite figure comes from [262].
quantities are typically non-linear and are of great relevance for the design of buildings and urban spaces. Due to this, a large effort has been invested in the analysis of these diagrams. The results suggest while qualitatively the functional forms of v(ρ) and J(ρ) have common features across different settings, quantitative differences emerge based on extraneous factors. Some of these are cultural perception of personal space (a comparison between India and Germany was performed in [260]), the presence of T-junctions [270], the individual's perception of open versus closed spaces [271] and the bidirectional flow of agents [272].

Air Transportation
On the opposite end of the spectrum from pedestrian mobility (which is necessarily short-range) is long-range transportation over thousands of kilometers, which is primarily represented by air transportation. Air transport is of paramount importance for modern global connectivity, playing a central role in the world economy and in the interchange of people, ideas and, unfortunately, also the propagation of diseases. While the study of such systems has a long tradition in the engineering community, primarily related to aeronautics and air traffic management [274,275], recent years has provided further impetus to this line of research with the development of tools based on network science [69]. In this setting, initial studies were conducted from a static perspective, where airports are the nodes of the network, with a link between a pair of airports if a direct flight exists between them, discarding any directional information as well as temporal dynamics. The basic statistical features of such networks were analyzed, for instance, in [276,277,273] for the World Airport Network (WAN), in [278] for the Chinese network and in [279] for Brazil. Figure 33 shows the location of the high-degree nodes (maximum traffic) and those with the highest betweenness centrality (bottlenecks) in the WAN. Statistical distributions of the degree (number of destinations of the airports), the traffic (number of passengers per route or per airport) and betweenness were found to be typically heavy-tailed. Furthermore, it was found that the networks can be divided into modules with a certain level of self-contained traffic [273], although this division is not seasonably stable when the focus is set on smaller geographical scales such as the US network [280]. The average number of passengers per route w i j was determined to be a nonlinear function of the traffic in airports of the form where k i is the traffic at the origin, k j at the destination and θ ≈ 1/2 [277]. An explanation of these measured topological properties has been proposed by simple models based on cumulative advantage [276,277]. The robustness of these networks to disruptions has been considered in terms of a unipartite framework [281] as well as multiplex-like approaches [282], where the network of each airline forms a layer and the airports appearing in the different layers connect the structure. A short review of network analysis on air transportation can be found in [283].
Of course the static topological framework does not take into account the dynamics occurring on the network. Indeed a key component of air transportation is that goods and people arrive in their destination in time. In particular logistical failures can propagate in the network where travel schedules are not met, with the delay propagating to effect a significant fraction of the system and paralyzing airport operations. This idea was first explored in the context of air transportation in [285], where the concept of delay multipliers was introduced. This can be seen as a combination of a delayed flight not only affecting the departure of other flights in the same airport, but also affecting connecting flights down the transportation chain, leading to the formation of avalanches or "delay-trees" [284] as shown in Fig. 34. Direct costs of flight delays in Europe amount to more than one billion euros [286], while in the US direct and indirect costs go beyond 40 billion dollars [287]. Delays directly affect airlines since they increase operational costs, but as an indirect factor they also bring associated reputational costs [288,289]. Passengers, on the other hand, undergo a direct loss of time, which can be further exasperated by lost connections leading to missed business and leisure opportunities. Finally, efforts to recover delays airborne usually imply excess fuel consumption and larger CO 2 emissions.
Given these challenges, it is important to characterize the sources of initial (primary) delays [290,291] and devise mitigating strategies. Unfortunately this is a highly complex problem with the involvement of numerous (and sometimes unrelated) factors such as weather conditions, technical failures in the aircraft, aircraft rotations, crew and passenger connections, organizational issues in airports, and the mechanisms through which these delays are transferred and amplified in subsequent oper- ations [292,285,293]. Airline schedules typically include a buffer time to absorb these delays [294], but there are frequent scenarios when this is not nearly enough. The Central Office for Delay Analysis of Eurocontrol (CODA) in Europe and its counterpart in the US, the Bureau of Transportation Statistics of the Department of Transportation, release monthly reports on flight performance including statistics on the major causes of delays and the airport affected. Based on these reports, studies have been conducted in single hubs such as Newark [295] as well as parts of the European [275,296] and US [289,297,284] networks. There have been primarily two approaches considered to model such effects. One class of models probe the network structure to search for weak points, analyzing jamming and congestion phenomena, especially trying to isolate those airports that cause cascades of delays that span a finite-fraction of the network [299,300,301,302]. This approach has also been extended to multiplex networks of airports and airlines also attempting to identify carriers that are most responsible for delays [282]. In addition to these more-stylized approaches, there exist other families of models whose objectives are to provide realistic predictions that can be implemented in real-time. Such models are primarily agent-based and incorporate fine-grained details of air operations while sacrificing analytical tractability. Indeed, this has been the main approach adopted in the area of transport engineering related to Air Transport Management (ATM). There are hundreds of contributions on this line mainly presented in sectorial conferences such as the USA/Europe ATM Seminar, the Air Transport Research Society (ATRS) Conference, ICRAT, the SESAR Innovation Days, etc, whose proceedings are accessible online as well as publications in journals such as Journal of Air Transport Man-   agement (JATM). The advances in this area cover topics such as the optimization of runway management, the effects of bad weather, sector congestion, how the advent of drone technology affects commercial air transportation, among many other minute details. While a comprehensive listing of all the findings would require a review by itself, we mention a few that are directly related to delay propagation in networks. Micro-modeling of the connection between air-transport logistics and delays was first explored in [303,304]. Later a similar framework was used to assess the economic impact of a major disruption in a European hub [305]. Also in Europe, within the umbrella of the WPE of the Joint SESAR Undertaking, several projects as NEWO, TREE and POEM have studied the modeling of delay propagation from different perspectives: the first two from the network and airline managers' and the latter from the passenger's point of view.
In the same context, data-driven models have been developed to reproduce the delay propagation patterns in the US [306]. These models were validated against empirical performance data [307], showing an acceptable level of accuracy and precision in the prediction of delayed flights, identification of airports displaying major departuredelays, network-wide congestion and generalized bad weather conditions [298]. The models are built with an agent-based approach applied at the level of aircraft, with delays appearing at any point in the operations due to technical, airport organizational or weather-related issues incorporated as initial conditions to the model. The airport management is simulated with a system of queues that depend on the specific network considered (US or Europe). One of the primary finding of such analysis is that the mechanism with the strongest potential to effect delay is passenger and crew connections (Fig. 35).

Sea Networks
Another example of long-range transportation are sea transportation networks, which have received comparatively less attention than their air traffic counterparts. While his-torically, sea transportation was the main mode of long-range human mobility, in recent times they mostly involve movement of cargo. Still, for sake of completeness we provide here a brief description on studies conducted on such networks. Typically such networks are studies as representing the ports as nodes, with links connecting them whenever major cargo liners travel from one to other [308]. The links typically are given weights that correspond to the number of ships (trips) along the route. The resulting network spans the globe, connecting different continents as shown in Fig. 36. Similar to that seen in air transportation networks, statistical distributions of typical topological properties are heavy-tailed, with only a few number of ports accounting for a majority of global traffic, with the majority playing a peripheral role [151]. The centrality of ports was studied in [309], where it was found that Singapore, Antwerp, Bushan and Rotterdam top the list in terms of both traffic and load. A similar relationship between connectivity and traffic as in Eq. (65) was noted with an exponent of 1.46 ± 0.1 [151]. Furthermore, the network can be segmented by the type of ship (containers, oil tankers, etc) in terms of a multiplex network framework [310]. A systematic analysis on the similarities between the world airport network and the boat network was performed in [311], including a study of the resilience of both networks. From a modeling perspective certain aspects of the flows can be reproduced with a gravity model (Sec. 4.2.1), although there are some strong limitations. An important additional application of this network, besides the transport of goods, is the study of the invasion of species brought by carriers from different parts of the world [312,313].

Intra urban mobility
One of the more important and much studied applications of mobility is that seen within cities. Indeed, there are good reasons for this. One is just the availability of more data. With rapid urbanization, an increasing part of the global population is living in urban areas. Large cities are early adopters of new technologies have large populations and consequently more mobile phone users and high densities of of phone antennas. Consequently the spatial resolution of CDR data (Sec. 2.3) collected in cities is order of magnitude higher when compared to rural/non-urban areas. The same considerations hold for the data generated from social networking applications, such as Twitter or Foursquare. Indeed, the usage of these applications has been shown to have socially and spatially biases; they are more likely to be used by urban residents, educated individuals, young people and middle age employed adults [314].
The other of course is the unique set of challenges that abound from large agglomerations of people resulting in high population density. One of these is congestion that has an associated energy cost as measured for example by C02 emissions. As a large part of human movement occurs in individual vehicles, accurate measurements of transportation flows are needed to infer how congestion affects carbon emissions. For example, it is known that CO2 emissions mainly depend on the time spent traveling, and in places of high density such as cities where congestion appears, traffic jams are responsible for carbon emissions that scale super-linearly with the population size [315,104]. Yet another reason for studying intra urban mobility is related to that fact that cities typically have high levels of socioeconomic inequality among their residents, Fig. 37: Statistics of the co-location times for "familiar strangers". In (a), typical contact network. In (b), inter-encounter distribution times. In (c), joint probability distribution of inter-encounter times. In (d), distribution time between encounters for groups. Figure from [316].
with residential segregation that is amplified by other forms of segregation, perceived in mobility practices and facilities that may differ depending on the neighborhood of residence. Additionally, the multi-modal nature of transport networks in urban areas (which include privately owned vehicles, buses, metros, taxis, tramways, bike sharing systems, and pedestrians) gives rise to specific questions related to the navigability of these multi-layer networks. For example, the combination of the regularity of daily circadian rhythms and smart transport card technology provides data that allows for the study of intriguing phenomena, such as familiar strangers [316], and can be used to uncover the properties of "hidden" temporal networks of the city, formed by individuals that temporarily co-locate in the same spaces, but do not know each other (see Fig. 37).
There is an important body of mixed literature, quantitative and qualitative, on intra-urban mobility that relies upon longitudinal travel surveys, started more than 60 years ago. Covering it is beyond the scope of this review, and possible starting points for those interested may be found in [317,318,319], which focus on North-American metropolises. Ref. [320] was one of the first papers to explicitly discuss the possibilities offered by individual digital footprints for urban analysis. It underlined that such data could be used to gain understanding of a number of phenomena associated to intra-urban mobility. Following this, a series of papers highlighted possibilities of real-time visualization and monitoring of displacements in cities [321,322,94]. The first studies that utilized metadata produced by handheld devices carried by individuals focused on a single city, and essentially consisted of visualizations and simple aggregated measures, providing information on the density of individuals using a device in various areas of the city and at various times of the day. Metadata provided by mobile phone operators contains information on the nationality of the phone number, which is also common on social networking user profiles. Consequently one of the first uses of such databases of individual geotagged data was to model the movements of foreign visitors in cities. [323] studied the location and mobility patterns of tourists in Rome, and subsequent studies highlighted the specific travel patterns of visitors, when compared to the movements of city residents, notably in Paris [324,325], Jakarta and Singapore [326].
Another popular case study are ones which focus on the influence an area of residence has on intra-urban mobility patterns. Most of the first intra-urban studies using ICT data in these categories have focused on the organization of mobility in one or two large cities. One of the very first in this category was [93], which used mobile phone CDR data recorded during a period D of 62 consecutive days in New York and Los Angeles. They compared aggregated mobility statistics; the individuals' mean daily range -defined for each individual as the mean of the maximal distance traveled each day, d∈D max i j (d(i, j))/|D| -and the individuals' maximal daily range -defined as max D (max i j (d(i, j))) (see Fig. 38). They showed that residents of Los Angeles globally commute over further distances, but also highlighted heterogeneities; some residents in Manhattan commute farther than their counterparts in Los Angeles and the neighborhood of residence has an influence on the average journey-to-work commuting distance. Combining mobile phone data and odometer data (a distance measuring device for bikes and automobiles), [327] performed similar comparative measures between the two data sources, highlighting the influence of a number of built environmental factors on the travel distances of Boston residents. [328] went further by combining mobile phone data and census data which provided information on the proportion of racial groups in each of the cities' census tracks. They showed that residential segregation was amplified by a mobility segregation in the urban areas of Boston and Los Angeles.
Focusing on a single city, some authors have investigated the relation between morphology and mobility patterns: [329] used smart cards data collected over several weeks in the London subway to provide measures linking the polycentric spatial structure of London and the organization of mobility patterns in the city. Relying upon mobile phone data collected during 2 months in 31 Spanish urban areas, [330] showed how the spatial organization of hotspots and the shape of cities evolve during the course of a typical weekday, as residents move inside the cities, highlighting different typical city forms (see Figure 39). Thanks to a large Foursquare dataset that records GPS positioning data, [65] compared the distribution of jump lengths in 34 different cities worldwide. Fitting the distribution in each city, they conclude that at the city scale, human displacements cannot be properly fitted by a power-law distribution, a result in agreement with several other case studies held for particular cities, which found that the distribution of travel distances inside an urban area is more properly fitted by an exponential distribution P(∆r) ∝ exp(−∆r/r 0 ). They propose an individual model derived from Stouffer's intervening opportunities, that states that the important parameter for explaining destination choice in cities is not the distance by itself, but rather the density of possible alternative destinations in a perimeter of distance d. Given a set U of places in a city, the probability of moving from place u ∈ U to a place v ∈ U is formally defined as where rank u (v) ∝ | {w : d(u, w) < d(u, v)} |, which can be seen as an urban counterpart to the radiation model proposed in the same year by [42] (see section 4). One of the most debated questions in the recent literature on human mobility has been the shape of the statistical distributions of the distances traveled by individuals P(∆r), and their travel times at various spatial and temporal scales and in different environments. Providing a neat characterization of these distributions is of extreme importance, and a prerequisite to modeling. Indeed, the distribution of travel distances P(∆r) is a key ingredient that any reasonable model should be able to reproduce and explain. [331] contains a very clear recap of the efforts of the last decade on this question, and a table summarizing the parameters values of the fitting distributions that have been proposed. These distributions have been characterized thanks to various sources of individual data (dollar bills, mobile phones, private car GPS, taxis, longitudinal travel surveys, and geolocated tweets, among others). In most cases, the fit of P(∆r) is consistent with the general form of a truncated power law distribution with κ = 0 corresponding to the case of a non-truncated power-law, and β = 0 corresponding to an exponential distribution of scale κ. The majority of papers propose truncated power law fits for P(∆r), and argue that the distribution of displacements follows an underlying Lévy-flight process.
The studies compared in [331] provided very different results for the parameters values β, r 0 and κ of the general form above, and consequently they proposed different underlying processes for explaining the shape of this distribution. To settle the discrepancies, a convincing modeling framework, able to reproduce a number of stylized facts, and P(∆r) is required, and more research and cross-checks should be performed. Data accessibility is always an issue in this type of research, so one workaround passes is repeating the same measures on many different datasets, to enable a systematic comparison of the empirical results as done by [331].
Apart from this, some of the many aspects studied on intra-urban mobility relate to the typical distance scales of movement [63,332]; measures of travel duration [332]; effects of population size on the structure and dynamics of mobility [104,333]; the effect of city topology on the displacements of its inhabitants (in terms of travel time, distance, and spatial organization of flows) [334] optimal 'strategies' for navigating multimodal transport networks [102], the efficiency of transport itineraries [335]; the level of complexity of transportation maps, which might exceed cognitive limits, and therefore explain sub optimal navigation by city residents [336].

Epidemic Spreading
Perhaps one of the most important applications of studies on human mobility and transportation systems relates to epidemic spreading. While epidemiology as a discipline has several hundred years of tradition, only recently has it been possible to model the global spreading patterns of diseases in a realistic and reasonably accurate way. Of course, every disease has its own peculiarities, and requires tailored modeling with the inclusion of the specific mechanisms associated with it. Consequently, the literature in this area is extensive, with dedicated monographs and a number of reviews [337,338,339,340]. Here we discuss how recent knowledge of human mobility has influenced the latest generation of epidemiological models.
Whilst it has been known for some time that air transportation plays a major role in the global spreading of diseases [342,343], the implementation of air traffic data in epidemic models has only been possible in recent times with the advances in information technologies and computing resources over the past decade or so [344,345,346], with the majority of studies focusing on the spread of influenza. The flu is a paradigmatic example of contact diseases that spreads annually cross the globe, and for which contagion is airborne with a patient clinical process that allows for relatively simple modeling. Indeed, improved surveillance systems have brought new insights on the temporal and spatial patterns of its propagation [163,347].
Of course, new modeling frameworks are not restricted to the flu but can be also extended to other contact diseases such as SARS [348], mosquito-driven infections lile malaria [349] or yellow fever [238], or even Ebola and other emerging health threats [350,351]. In many such cases, the first instance of the disease in a particular location occurs through the arrival of infected humans from the source location via air transportation. The world airport network (WAN) can thus be considered as the backbone through which new emerging diseases arrive where connections with larger flows of passengers marking the most likely pathways for propagation [341] (see Fig. 40). Indeed, within the framework of the WAN (as discussed in Sec. 5.1.2), it is possible to estimate the potential of every airport as a source for a pandemic [352]. Concerning realistic models, two main frameworks have been used in the literature: agent-based modeling and metapopulations. The difference between the two lies in the level of information needed. Agent-based models attempt to generate synthetic population mimicking at an individual (agent) level the realistics aspects related to disease propagation. A detailed knowledge is therefore required of features such as location of residential areas, businesses and schools, leisure activities, population divisions by age, gender and sizes of households as well as individual mobility. These are typically sourced from census-like surveys and then a synthetic population is created assigning each individual to a household, a work/school place and a characteristic mobility profile. The earliest variants of such models were developed in the 1990's within the Los Alamos project Transportation Analysis Simulation System or TRAN-SIMS. The objective of TRANSIMS was to reproduce via agent-based modeling the daily life of Portland, Oregon, including population, mobility, income levels and other details sourced from the US census [354]. This was later adapted to study the evolution of epidemics, called Epidemic Disease Simulation System or EPISIMS, which was used to study different scenarios to prevent the propagation of smallpox [355]. Similar models were applied to study the propagation of a strain of avian flu, H5N1, in South East Asia, where for instance, epidemic spreading was simulated with the source centered in Thailand (see Fig. 41), and including neighboring countries at a range of 100 kms containing about 80 million inhabitants (agents) [353,356]. This was primarily done to assess the effect of several intervention measures (distribution of antivirals and reducing the social contacts) to halt the progress of the disease. Similar models have been applied to even larger populations in the US and UK [357,358,359,360]. The level of granularity present in such models, enables minute analysis such as the effect of closing specific schools, employing different vaccination protocols (say by age) on disease propagation [358,359]. The arrival of the 2009 flu pandemic provided further impetus to such modeling approaches and agent-based models have been now developed for most of Western European countries [361,362,363,364] where a dazzling amount of scenarios have been exhaustively simulated including the risk of spread of a newly created virus strain from a lab following the case of the creation of a new H5N1 influenza strain [365]. One must note that the mobility aspect in such models relates to the movement of the individuals (agents) and is explicitly implemented both at long and short ranges involving multiple models of (simulated) transportation.
While agent-based models provide the best approximation to reality, they are severely limited by the fact that they require highly detailed input information, the data for which in many cases is noisy, unreliable or just simply unavailable. Furthermore, given the sheer complexity of such models, an analytical treatment is nearly impossible. On the other hand, the metapopulation framework, balances granularity with analytic tractability, and has easier to obtain data as input parameters [342,366,367]. In this model, the population of interest is divided into subpopulations, within which contacts are modeled by a fully-mixed mean-field approach, whereas contacts between subpopulations occur through the (measured) mobility networks. One of the earliest such models is the so-called Global Epidemic and Mobility model or GLEaM [368]. GLEaM is a global simulation platform that divides the globe in cells of 15×15 minutes of arc, which are grouped in basins of attraction around the major transportation hubs, which in most cases are airports. The mobility is then simulated at two levels: long range air transport, which is calibrated from the world airport network, and short-range mobility obtained from either census commuting data or estimated using the gravity law [165]. Figure 42 displays the different components of the model.
Most notably, GLEaM was used to study the propagation patterns of the 2009 H1N1 influenza pandemic in real time. The pandemic started in a small population of Mexico early in that year [369], and its unfolding received attention of mass media media as well as the scientific community, which brought unprecedented details on the characteristics of the first imported cases in several countries. This information allowed for the calibration of the infectivity and the seasonality in GLEaM in early May 2009 [368]. Once calibrated, the model predicted the disease propagation patterns: essentially, the arrival and prevalence peak times for the rest of the countries of the world, predictions that were later validated, after the full surveillance reports were available in 2012 [347]. Additionally, the model helped estimate the actual initial number of cases present in Mexico in the early stages of the disease [370]. Ever since, the model has become an important tool for subsequent studies that takes its results as a baseline (see, In both cases, the link weights are plotted using a heatmap code. In (c)-(f), comparison between flows coming from commuting data w D and those produced with a gravity model w M with the aim of extrapolating short-range mobility in zones where the data was not available. Figure from [164].
for example [341,352]) and has been also used to assess the risk of spreading of new emerging health threats [371,351,372].
The flexibility of the multilevel framework inherent in the model allows for the quantitative assessment of the differences between short-range and long-range traveling on the disease propagation. In particular it was determined that travel restrictions via decrease of the air traffic (or temporary quarantines) can at best delay the disease spread by a few days or weeks. Unless one totally shuts down the network, or severely restricts it (which of course carries a massive economic cost) global propagation is almost guaranteed [373,372]. Newer versions of such models factor in facets of human behavior, such as the tendency of people to restrict contacts and travel when made aware of the existence of a dangerous disease [374], as well as heterogeneity in the population (age, gender, socioeconomic indicators) and its effect on travel patterns [375]. Possible future directions regarding these models include the possibility of feeding them with real-time ICT collected data on mobility [376,377,20,378] or hybrid frameworks that merge aspects of agent-based and metapopulation modeling [361].

Web (Online) Mobility
Ostensibly, physical human mobility and virtual navigation are indeed activities of very distinct nature, yet there are a number of reasons why the latter can be cast in the context of the former. Indeed, many daily activities that occurred in specific physical locations, such as work, study, leisure, financial transactions, purchasing, are now increasingly migrating towards being conducted online. Moreover, the world wide web (www) has a humanitarian and social role, in the sense that it helps individuals and populations with limited mobility to have access (albeit virtually) to locations, people and information hitherto unreachable. For instance, the disabled [379,380,381] and the elderly [382,383,384], at an individual level, for whom the www represents a means of overcoming their physical limitations, or even entire populations, whose mobility is restricted for political or geographical reasons, for which the www is their only access to the outside world [385]. In this section we present an overview of the recent empirical evidence on the statistical similarities between virtual and physical navigation.
In the context of this parallel between virtual and physical mobility, we differentiate two classes of activities that although related, have practical and functional distinctions, that is, mobility in virtual spaces (for example, World of Warcraft, Second Life and The Sims) and navigation in virtual spaceless environments (for example, browsing the WWW or navigating through the interface of a device). The difference between these two classes of activities is given by the fact that in the first there is an explicit notion of space such that displacements between different locations come with a temporal cost, often proportional to some distance measure between them [387].
An empirical analysis of the behaviors and trajectories in the virtual space of the online game Pardus suggests a subdiffuse behavior in online mobility [386], similar to that observed in physical mobility [81,48] ( Fig. 43). Additionally, it appears that the return time (in number of discrete steps) of users of the Pardus game follows a power-law distribution with an exponential cut-off (inset of Fig. 43). These similarities between virtual and physical mobility were later empirically corroborated by other similar datasets. A comparative analyses between traces of World of Warcraft (WoW) players and GPS traces [388] revealed that the distributions of jump lengths and waiting times are both approximated by lognormal distributions (Fig. 44).
Furthermore, the analyses of CDRs along with the Web activities of 20,000 mobile phone users [389] suggest that visitation frequencies in both activities have a stretched exponential distribution (see Fig. 45). However, Web traces from mobile phone users tend to overrepresent a particular set of activities such as fact finding, information gathering and communication [390]. More recently, the analyses of complete Web browsing activities of 582 users, has revealed that it is possible to profile Web users based on their exploratory vs. exploitative behaviors in their browsing patterns [391]. More surprisingly, however, is the fact that the distribution of returners and explorers profiles are remarkably similar to the one observed in human mobility [87]. The observed similarities between physical and virtual navigation is beginning to draw increased attention from the scientific community. Of course the analogy can oly be taken so far, as compared to physical movement, virtual activity is relatively recent (on an evolutionary time scale) and much remains to be understood. It is reasonable, however, to believe that online activity is driven by the same psychological mechanism as off-line behaviors, now adapted to new environments or contexts. In fact, such intuition has increasingly found scientific support, both from a theoretical and empirical point of view, such as evidence of foraging behaviors in online shopping [392] and information consumption on the Web [393,394], similar to that observed in human resource acquisition behaviors [395]. (preferential re [12] dynamics, we hypothesize two underlying human activities in cyber exploration and preferential return, a

Conclusions
In this work we discussed the state-of-the art in the field of human mobility. After a short historical description of the field, we focused on current work starting with the type of data utilized in most research projects, followed by the description of metrics and models of individual mobility and population mobility. We illustrated a list of selected applications of human mobility modeling, demonstrating the insights the field can bring to the solution of real-world problems.
It should be clear from the survey that the study, understanding, and modeling of human mobility is of primacy to many areas of society. At a deeper level the current research raises a fundamental question regarding human cognitive limitation. Among the many features described here, one in particular stands out: Why are we likely to return to frequently-visited locations or recently-visited locations? Does it relate to our ability to remember where we have been? Is it a innate characteristic of individuals or something that can be learned (i.e. nature vs. nurture)? The same can be said about regularities whereby human behavior is divided between those that are explorers or returners in their mobility pattern [87]. Indeed, this begs the interesting question, as to whether this observed division is an inherent property of society, say acting as a stabilizing influence? Indeed, this dichotomy is also found in other context, such as social learning, where it has been suggested that a population cooperating to solve a complex task can attain a better performance if it is composed by a mixture of explorers, individuals who try new solutions, and exploiters, individuals that have less propensity to abandon the optimally extant solutions [396]. More generally, the exploration vs exploitation dilemma is a fundamental issue in optimal decision-making and is analyzed in economics (multi-armed bandits) and computer science (reinforcement learning) and finds natural applications in marketing and advertisement, where recommendation algorithms determine whether people are likely to adopt familiar or unfamiliar products. All these problems require an understanding of human dynamics and many of these areas can possibly be studied using approaches inspired from research in human mobility. With the data we continue to collect as a society, and given the methodology developed to explain regularities in human mobility, we now have an additional set of mathematical tools to characterize the regularities in human activities more broadly.
A different structural dimension linked to human mobility relates to the invariant that humans are linked to each other forming social networks. Studies of social networks aim at grasping the multidimensional properties of the structures formed by entities (e.g., persons and organizations) and the ties connecting them together [397]. In fact, the intersection between human mobility and social networks goes beyond the obvious analogy that humans are social and their movement may be related to them being social [398,399,400,401,402]. The field of Social Network research is a huge enterprise in its own right and elucidating the connection between social networks and human movement is a nascent effort that should be explored in more depth.
An interesting research frontier in human mobility in the near future will likely concern the shift from traditional vehicles to autonomous, self-driving, vehicles. The diffusion of autonomous vehicles will transform private and public transportation, with sweeping consequences for society, the economy and the environment. How exactly will these transformations happen and how will mobility habits change in response to these new technologies? How can such vehicles be controlled and how does one plan routes in order to optimize fuel and time consumption, reducing congestion and pollution? Answers to these important questions, may well arise from a fuller understanding of human mobility and behavior at the individual and population (collective) levels.
Finally, it should be clear from this survey that interdisciplinarity and cross-fertilization is a condicio sine qua non to the success of the field. Computer scientists, physicists, social scientists, environmentalists, engineers, government officials and many other professions need to work together to develop and implement solutions to many of the problems society faces today that could benefit from understanding mobility (and dynamics) such as crime, urban planning, energy consumption, social integration, to name but a few. or another, on physical mechanisms and first principle formulations. However, what all those mobility models have in common is the fact that they are based on Monte Carlo methods, in the sense that they rely on generating repeated random samples as an approximation for the behavior of the system (in this case, human displacements).
Computational models can come in different flavors, varying in many aspects such as modeling approach and complexity, and are used to simulate human mobility to validate the analytical models by comparing data generated from the simulation with real world traces. Simulations can be grouped into three major modeling approaches: numerical, particle-based and agent-based models (ABM). Also, there are simulation of discrete events and simulation based on system dynamics [403]. The former consists in conducting entities through processes or queues in order to analyze the general behavior exhibited by a system over time. In the latter, simulations are performed by analyzing the behavior of variables in systems of differential equations [404]. Historically, the first simulations were numerical in nature and although effective, in many cases the results of the simulations of the models did not have good adherence to reality, especially in cases where many abstractions are needed (i.e., intentional disregard of aspects of systems to simplify the model) [405]. The rise in popularity of ABM models, especially to model human mobility and social phenomenas, is due to the fact that they are more suitable when the simplicity of a Monte Carlo method is not sufficient or appropriate to represent more complex entities such as adaptive processes or local interactions. Also ABMs are especially useful when second-order dynamics or emergent phenomena are involved. In the context of human mobility, elements such as past-events memory, social networks and spatial capabilities make ABM the natural modeling approach.
Agent-based models are a class of computational models that emphasizes interactions of entities to assess their effects on a system. The key point is that simple rules and micro-level interactions lead to the emergence of complex macroscopic regularities and organizations that represent the observable outcomes to allow the study of the dynamics of individual behavior. Computational agents, according to Wooldridge and Jennings [406], normally have the following characteristics: Autonomy The agents' operation must take place without the direct control of its actions or internal states; Social behavior Agents must interact with other agents through some kind of standard communication; Reactivity Agents must be able to understand the environment and respond to it; Initiative In addition to reacting to the environment, agents are also able to take the initiative of an action, changing their behavior so as to fulfill a certain purpose.
Based on these aspects, simulations carried out with agent-based models allow for a deep investigation on the system from microlevel behaviors to macrolevel phenomena. Thus, agent-based models, depending on the kind of model to be designed, may show several benefits over other modeling techniques. According to Bonabeau [407], they may be summarized as follows: r ← POWER-LAW(α, x min ) 8: x ← x + r cos θ 9: y ← y + r sin θ 10: n ← n + 1 11: end while 12: end procedure Continuos-Time Random Walk (CTRW). A limitation of the RW models introduced so far is that time intervals between two successive steps are constant (Sec. 4.1.3). In order to simulate a CTRW we can modify Algorithm 3 to include the wait time τ between jumps (Algorithm 4). Notice that the wait time and jump length need not to be drawn from power law distributions, but CTRW with power laws has been shown to be useful to model human mobility [46]. Moreover, the power laws can have different parameters.
Exploration-Preferential Return (EPR). Although there is an extensive literature supporting the idea that human trajectories indeed follow reproducible scaling laws, the predictions of Lévy flight and CTRW models systematically deviate from empirical results. Song et al. [81] proposed two key mechanisms, unique to human trajectories, that were missing from existing models: Exploration and Preferential Return (see Sec. 4.1.4). Another aspect that differs from previous study is the use of power law with exponential cutoff for the jump length and wait time ( Fig. 48(a)). Power law with cutoff does not have a closed-form expression. However, we can generate a random number drawn from a power law with parameters x min , α and an exponential cutoff K with the following process [409]: generate an exponentially distributed random number x = x min − 1 λ ln (1 − r), where r is a uniformly distributed number in the interval [0, 1). Then, accept or reject it with probability p or 1 − p respectively, where p = x Repeating the process until a number x is accepted will guarantee that it will follow the desired distribution. This procedure is illustrated by Algorithm 5. Note that in the pseudo-code x and y are lists, therefore the operator + has the meaning of concatenation of the two operands. Finally, if we define S as the set of previously visited locations, F i the number of times (frequency) that the location i has been visited, and ρ and γ the parameters that control the probability P new of an exploration jump then the EPR random walk is implemented by Algorithm 6. Note that while the x and y coordinates vary continuously, the space is actually divided in discrete "locations", like Voronoi patches or squares of a grid. Therefore, different points in the plane may be part of the same location; to this purpose we use the function CURRENT-LOCATION(x,y) to obtain the location associated with the current coordinates of the agent. r ← POWER-LAW(α, x r min ) 8: x ← x + r cos θ 9: y ← y + r sin θ 10: τ ← POWER-LAW(β, x τ min ) Wait Phase 11: while τ > 0 do 12: if τ >= 1 then 13: WAIT (1) 14: τ n ← n + 1 20: end while 21: end procedure Recency. A further improvement over the EPR model has been introduced by Barbosa et al. [112]. The main observation is that frequently visited locations are also recently visited locations (Fig. 48(b)). However, recently visited locations might not always be frequently visited, because it might be that a location has been visited for the first time (Sec. 4.1.5). The algorithm to implement the recency mobility model is exactly the same as for EPR, but when deciding the for a return jump we can choose either a frequency based return with probability α or a recency based return with probability 1 − α, where α is a parameter extracted from empirical data (Equation (43), Algorithm 7).
Social Models. It is intuitive to understand that the mobility of socially related individuals is correlated; for example friends often hang out and visit places together. Therefore, mobility models that include social network links and proximity based interactions were developed (See Sec. 4.1.6). Two prominent examples are the models described in [140] and [141], which can be seen as variations of the preferential return model (Sec. 4.1.4), where the preferential return is governed by the social strength between the users instead of the visitation frequency of the locations. Given the previously presented building blocks and algorithms, it should be trivial at this point to implement these models, therefore the explicit algorithms are omitted.
In [140], at every step each agent chooses to visit a social contact (i.e., another agent with whom the agent has a link) with probability p v . Otherwise, it chooses a new location with the complementary probability 1 − p v following a Lèvy-like flight. After Algorithm 6 Exploration-Preferential Return (EPR) Random Walk 1: procedure EPR(N) 2: x, y, n, S ← 0 3: while n < N do 4: i ← CURRENT-LOCATION(x,y) 5: if i ∈ S then if RANDOM(0, 1) < P new then Exploration Step 13: θ ← RANDOM(0, 360) 14: r ← POWER-LAW-CUTOFF(α, x r min , K r ) 15: else Frequency Based Return Step 16: r, θ ← FREQUENCY-RETURN() 17:

end if
18: x ← x + r cos θ 19: y ← y + r sin θ 20: τ ← POWER-LAW-CUTOFF(β, x τ min , K τ ) Wait Phase 21: while τ > 0 do 22: if τ >= 1 then n ← n + 1 30: end while 31: end procedure 32: function Frequency-Return 33: tot ← j∈S F j 34: throw ← RANDOM(0,tot) 35: partial ← 0 36: for j ∈ S do 37: partial ← partial + F j 38: if partial < throw then 39: r ← EUCLID-DIST(i,j) 40 k ← ROUND(POWER-LAW(ν)) select the rank from a zipfian law 3: while k ≥ SIZE(S ) do In the recency model, a user can decide to go back to a recently visited location even though it has been visited only once or few times, therefore the visitation frequency is more evenly distributed. Furthermore, the spatial pattern shows several clusters further from the starting location, similarly to what happens with Lévy flights.
the movement, the agent creates a new link with probability p in the "neighborhood", that is the nearby agents, and probability p c with a random agent. An agent can easily keep track of the social links with the other agents by updating a list of the agent's ids it meets. In [141], the first part of the mobility model is identical to the preferential return model, where the agent returns to a previously visited location with probability 1−ρS −γ or visit a new location with probability ρS −γ . However, it adds an additional "social" step. The actual choice of the location to visit, either new or returning, is made based on the social influence with probability α or, with probability 1 − α, just the individual preference (i.e., proportionally to the location visitation frequency). In case the location is chosen according to a social contact, the social contact j is chosen with a probability directly proportional to the current mobility similarity between the two agents. The mobility similarity between two agents i and j is defined as the cosine similarity in the location visiting profiles: where v i and v j represent the total visits made by the agents to each location in the simulation space.

A.2.2. Population Mobility
A simple way to formalize the flow of individuals from one location to another is to divide the area of interest into different zones labeled as i = 1, . . . , n and to count the number of individuals going from location i to location j, ∀i. These numbers T i j are the elements of the so called Origin-Destination (OD) matrix (Sec. 3.1.5) Such matrix defines a directed and weighted network, and in the general case is timedependent. Note that usually n = m and T i j = 0, ∀i = j. The OD matrix is very different from a "segment" measurement, where we can easily count the number of individuals going through a link of the transportation system under consideration (e.g., one of the flights of the airline network, a segment of road, etc.). In contrast, the OD matrix is usually extremely difficult and costly to obtain and measure. There are obviously many factors which control the origin-destination matrix: land use, location of industries and residential areas, accessibility, etc. In order to simulate the OD matrix using the multi-agent modeling framework we need to define a process that drives the movement of the agents from one location to another. Let's consider a discrete-time agent simulation for which at every time tick t we know the elements T i j of the matrix OD(t). T i j can be considered as a frequency count, therefore we can transform the elements of the OD matrix by computing the relative frequency P i j : P i j = T i j i j T i j P i j represents the probability of an agent to move from location i to location j at a specific instant in time t. If we want to run a simulation for an amount of time τ, for each time instant t the simulation algorithm normalizes the matrix OD(t), then loop through each element P i j of the normalized matrix OD and move an agent, MOVE-AGENT(i,j), from i to j with the probability P i j (Algorithm 8). This framework is general and applies to any mobility model that uses origin destination matrices such as the gravity models familiy (Sec. 4.2.1), intervening opportunities, and radiation model. However each model may be different in the way that the elements T i j are computed. It is also interesting to note that if we represent the movement of the agents over e.g., a map layer, then the flux of agents over the connection between two locations should represent and estimate of the intensity of the traffic.