Gravity model explained by the radiation model on a population landscape.

Understanding the mechanisms behind human mobility patterns is crucial to improve our ability to optimize and predict traffic flows. Two representative mobility models, i.e., radiation and gravity models, have been extensively compared to each other against various empirical data sets, while their fundamental relation is far from being fully understood. In order to study such a relation, we first model the heterogeneous population landscape by generating a fractal geometry of sites and then by assigning to each site a population independently drawn from a power-law distribution. Then the radiation model on this population landscape, which we call the radiation-on-landscape (RoL) model, is compared to the gravity model to derive the distance exponent in the gravity model in terms of the properties of the population landscape, which is confirmed by the numerical simulations. Consequently, we provide a possible explanation for the origin of the distance exponent in terms of the properties of the heterogeneous population landscape, enabling us to better understand mobility patterns constrained by the travel distance.

Introduction

For understanding the mechanisms of human mobility [1-3], optimizing the mobility flows [4], and predicting the dynamics on mobility networks [5-7], a variety of mobility models have been extensively studied [8], such as gravity model [9], intervening opportunities model [10], and radiation model [11]. Among these models, the gravity model has been widely used for predicting the traffic flows between populated areas. The gravity model predicts the traffic flow between an origin and a destination in terms of a simple formula, similar to Newton’s gravity law, using populations of the origin and destination as well as the geographical distance between them [9, 12, 13]. Precisely, the traffic from a site i to another site j is given by where m (m) denotes the population of site i (j) and r is the distance between sites i and j. The value of distance exponent γ is found to range from 0.5 to 3 for several data sets [13]. This original gravity model and its variants have been applied to human mobility and transportation [11-21] ranging from the individual level [22] to the international level [6], and other datasets such as international trade [23], scientific collaboration [24], and mobile phone communication [18, 25] to name a few, mostly due to their simplicity. The gravity models nevertheless have limitations such as the absence of universality regarding the exponent estimation [11].

In order to overcome these limitations of the gravity models, Simini et al. [11] recently suggested the radiation model, similar to the intervening opportunities model, that considers the opportunity for travelers rather than the distance traveled. By employing the radiation and absorption processes of particles, the radiation model describes the mobility patterns without any parameter estimation. Precisely, the traffic from a site i to another site j is given by where T is the outgoing traffic from the site i and s is the total population, except for the sites i and j, within a circle centered at the site i with radius r [11]. The radiation model has several advantages compared to the gravity model such as clear theoretical background, universality due to the absence of parameters to be estimated, and better prediction for long-distance travels, despite some unresolved issues like relatively poor predictability on short-distance travels [17]. Moreover, the radiation model requires additional information on T, in contrast to the gravity model. The variants of the radiation and intervening opportunities models, e.g., a population-weighted opportunities model [26] and a radiation model with an additional scaling exponent [27], have also been studied.

The radiation and gravity models have been compared with each other, often together with other mobility models, in terms of the predictability of mobility patterns observed in various empirical data sets [17, 18, 28]. Here we raise a question: Beyond the comparison, can these radiation and gravity models be more fundamentally connected to each other? The possibility of such connection was briefly argued by Simini et al. [11, 29] such that the surrounding population s was assumed to be proportional to in the case with the uniformly distributed population, and later to be proportional to with the fractal dimension d of the population. These assumptions lead to the asymptotic values of γ = 4 and 2d, respectively. However, the population landscape in reality can be characterized not only by a fractal geometry of populated areas or sites but also by a power-law distribution of the population at each site. In this paper, we first devise a population landscape model characterized both by a fractal dimension d and by the power-law exponent β of the population distribution, and then derive the distance exponent γ as a function of d and β from the radiation model on our population landscapes, which we call the radiation-on-landscape (RoL) model. We also show that the distance exponent can vary according to the population sizes of origin and destination sites. These results shed light on the connection between gravity and radiation models. More importantly, we unveil the origin of the distance exponent in the gravity model in terms of the properties of the heterogeneous population landscape, provided that the radiation model is correct. Therefore we can better understand the mechanism behind the traffic flows constrained by the travel distance.

Results

Modeling heterogeneous population landscapes

As for the properties of heterogeneous population landscapes, we consider the fractal geometry of cities and the power-law distribution of their populations, both of which are well-known characteristics of human settlement. On the one hand, the fractal geometry suggested by Mandelbrot [30] has been applied to the landscapes of human settlements in several states of the United States of America [31] and over the world [32]: The fractal dimension in those datasets is found to range from 1.4 to 1.9. The fractality has also been studied regarding the inner structures of cities [33-35] and their growth patterns [36-39]. On the other hand, the power-law distribution of urban populations was presented in the classic paper by Zipf [9] as well as in a number of recent studies [40-44]. The power-law exponent of the population distribution of cities is found to have the value ranging from 1.7 to 3 [9, 40, 42, 45]. Despite the ongoing debate on whether populations are characterized by a power-law or a log-normal distribution [43, 46, 47], the power-law distribution would be still a reasonable assumption for model studies.

For modeling the heterogeneous population landscape, we first generate a set of sites in a two-dimensional space with a fractal dimension d. Then we assign to each site i the population m independently drawn from P(m) ∼ m− with an exponent β, which will be called the population exponent. Note that the geometry of the sites can be implemented irrespective of the functional form of P(m). In our work, we focus on the case in which the location and population of each site are fully uncorrelated with each other.

In order to generate a fractal geometry of sites, we employ the Soneira-Peebles model [48], originally developed for simulating the self-similar galaxy distribution. The model on the two-dimensional space iteratively locates sites within each circle centered at the site in the previous layer whose radius is decreasing as the layer deepens, see Fig 1(a). Precisely, we consider a circle centered at the origin with radius R. Within this circle, η > 1 sites are randomly placed in the first layer and each of these sites is assigned a circle with radius R/λ, with λ > 1 denoting the contraction factor between layers. The same process is repeated until the depth of the layer reaches L, eventually leaving us with η sites in the Lth layer. Here L denotes the number of layers. In our work, we consider the set of sites only in the last layer to find its fractal dimension as [49] Once the set of N = η sites with a fractal geometry is generated, we draw N independent values from a population distribution P(m) to randomly assign them to the sites. As for the population distribution we adopt the power-law distribution with the population exponent β > 1: where m0 is the lower bound of the population. We set m0 = 100 to scale the population to a realistic size. Fig 1(b) shows an example of the generated population landscape in the two-dimensional space using η = 2, λ = 21/1.5 (i.e., d = 1.5), L = 13, and β = 3. The height in the vertical axis indicates the population assigned to each site. Although there exist many other modeling approaches for generating heterogeneous population landscapes [36, 39, 50, 51], we have adopted the Soneira-Peebles model for the fractal geometry, mostly because the implementation of this model is efficient and scalable.

Fig 1

Modeling heterogeneous population landscapes.

(a) Schematic diagram of the Soneira-Peebles model in the two-dimensional space with η = 3. The number in each circular symbol denotes the layer which it belongs to. The sites at the first layer (blue) are randomly placed within the circle with radius R. Similarly, the sites at the second layer (red) are randomly placed within the circles with radius R/λ. (b) An example of the generated population landscape using the Soneira-Peebles model with η = 2, λ = 21/1.5 (i.e., d = 1.5), and L = 13, and a population distribution P(m) ∼ m− with the population exponent β = 3. The height in the vertical axis represents the normalized value of the population assigned to each site.

Connecting the radiation-on-landscape model to the gravity model

The connection between radiation and gravity models can be made by the observation that the surrounding population s of the radiation model in Eq (2) might be correlated with the distance r of the gravity model in Eq (1). The relation between s and r can be analytically derived in our population landscape model. Using this relation, the radiation model in our population landscape, i.e., the radiation-on-landscape (RoL) model, can be described by Eq (2) but in terms of r. By expanding the RoL model with respect to r, one can derive the distance exponent γ as a function of the fractal dimension d and the population exponent β of population landscapes.

Scaling behavior of surrounding population

We first remind that the surrounding population s is defined as the total population, except for the sites i and j, within a circle centered at the site i with radius r. Let us denote by Λ the set of sites, except for i and j, within a circle centered at the site i with radius r, and the number of sites in Λ is denoted by n. In a d-dimensional space, one can write as with a coefficient c. The surrounding population is written as where m denotes the population of the kth populated site in Λ, such that . As all ms are statistically independent of each other, one can relate m with its rank k using P(m) as [52] From Eq (4) we have where we note that β > 1, leading to Therefore one gets where we have used Eq (5). When β > 2, the term of dominates s for large r, while the term of does for β < 2. Therefore, we obtain the scaling relation between s and r for large r: with

Expansion of the RoL model

The relation between s and n in Eq (9), together with the relation between n and r in Eq (5), allows us to rewrite the radiation model in terms of r, i.e., the RoL model. From Eq (2) we define the travel probability as and the rescaled travel probability as For the expansion, we consider three cases: (i) m, m ≪ s, (ii) m ≪ s ≪ m, and (iii) m ≫ s.

(i) If m, m ≪ s, the rescaled travel probability is expanded as Here we find the leading term of from Eq (11), leading to This scaling form of the distance dependence enables us to compare our RoL model to the gravity model in Eq (1): By comparing the distance dependence of the RoL and gravity models, we obtain the distance exponent γ as a function of the fractal dimension d and the population exponent β: Note that the result of γ = 2d has been suggested in a previous work [29].

(ii) If m ≪ s ≪ m, one gets From the leading term of , we obtain

(iii) Finally, if m ≫ s, one has irrespective of m. Since the leading term is independent of r, we have However, the subleading terms are still functions of r, leading to a weak distant-dependent behavior of the rescaled travel probability.

From the above analysis, it is remarkable to find how the distance exponent γ can vary according to the population sizes of origin and destination sites, i.e., m and m, respectively. This strongly implies that a given data set does not necessarily have to be characterized by the single value of the distance exponent. In reality, travelers from small towns may have different reasons for selecting their destinations, hence different travel distances, than those from big cities; the population size of the destination can also affect the traveling behaviors.

We provide an intuitive explanation for our results in Eqs (18) and (20). We remind that in the gravity model, the distance exponent γ plays a role of spatial cost in determining the traffic flows because the larger γ leads to the stronger dependence of the traffic flows on the distance. Let us consider a job-seeking situation as in the original radiation model [11]. Since the number of cities is proportional to , a higher-dimensional geometry with a larger d would provide more opportunities in the same range of r from the origin. It implies that a job-seeker can find a job at a closer city and does not need to travel farther in a higher-d space, leading to a larger γ. Dependency of γ on the heterogeneity of the population distribution can also be understood with the job-seeking example. In the original radiation model, a place with the larger population provides more opportunities, and a job-seeker finds a job at the closest city providing the better opportunity than the origin. For example, let us consider a homogeneous case with 10 medium-sized cities with two workplaces per city, which can be contrasted to a heterogeneous case with one extremely large city with 11 workplaces and nine small cities with one workplace per each. Then the job seekers in the homogeneous case tend to travel to any other cities providing a little better opportunities, implying a smaller γ. In contrast, the job seekers in the heterogeneous case tend to travel only to the extremely large city and do not have to travel farther than that city, implying a larger γ. Since the smaller β implies a more heterogeneous population distribution, one can relate the smaller β to the larger γ, closing our arguments for Eqs (18) and (20).

Numerical validation

We numerically test the analytic results using the heterogeneous population landscapes described in Fig 1. We generate 100 different population landscapes with the same parameter set of η = 2, λ = 21/1.5 (i.e., d = 1.5), R = 1, and L = 13, then assign to the sites the populations drawn from P(m) ∝ m− in Eq (4). We also set the upper bound of m as 107. Once the population landscapes are generated, one can calculate for every pair of sites i and j the distance r, the number of sites for the surrounding population n, the surrounding population s, and the travel probability p using the following Eq (23) for the finite system. The travel probability for the finite system [17] is given by where M ≡ ∑ m denotes the total population in the system. As almost all ms are much smaller than M, our analytic results remain valid.

Surrounding population

The results of (r, n) for all possible pairs of sites i and j are depicted as a heat map in Fig 2(a), from which we estimate the fractal dimension and the coefficient in Eq (5). Here the scaling behavior is observed in the intermediate range of r. The lower bound of this range is related to the smallest length scale, i.e., R/λ ≈ 10−2 for the parameter values used, while the upper bound is related to the largest length scale, which is trivially R = 1.

Fig 2

Properties of the heterogeneous population landscapes.

(a) Numerical validation of the scaling relation between n and r in Eq (5) for the fractal geometry of sites generated using Soneira-Peebles model with η = 2, λ = 21/1.5 (i.e., d = 1.5), R = 1, and L = 13, averaged over 100 different landscapes. (b, c) Numerical validation of the analytic relation between s and r in Eq (9), together with Eq (5), on the same fractal geometry of sites as in (a), but also with P(m) ∼ m− for the values of β = 1.5 (b) and of β = 3 (c), respectively. For each gray-colored heat map, the darker color implies more pairs of sites around the point (r, n) or (r, s). The log-binned curve (red circles) of the heat map is compared to the corresponding equation (light blue curve).

Similarly, from the results of (r, s) for all possible pairs of sites i and j, we get the heat map for a few values of β, as shown in Fig 2(b) and 2(c). When log-binned, it gives the curve of s as a function of r, which turns out to be comparable to the analytic result in Eq (9) when using estimated values of and for both cases with β < 2 and β > 2. Accordingly, the scaling relation between α, d, and β in Eq (12) is also validated.

Rescaled travel probability

Next, we test the validity of the expanded forms of rescaled travel probabilities in Eqs (15), (19), and (21), by comparing them to the numerical results on the generated population landscapes using Eq (23). In particular, for studying the effects of origin and destination populations on the scaling behavior of the rescaled travel probability, the sites are decomposed into 10 groups according to their population sizes, denoted by G for v = 1, ⋯, 10. Then all pairs of origin and destination sites can be decomposed into 100 groups of pairs, such that G = {(i, j)|i ∈ G and j ∈ G} for v, w = 1, ⋯, 10. For each group of pairs, say G, we calculate the rescaled travel probabilities for all pairs in G using p in Eq (23) to obtain a heat map for (not shown). By log-binning the heat map, one gets the curve of the rescaled travel probability as a function of r for each G, as shown in Fig 3. We find that these numerical results are in good agreement with the expanded forms of rescaled travel probabilities for m, m ≪ s in Eq (15), for m ≪ s ≪ m in Eq (19), and for m ≫ s in Eq (21), respectively. Accordingly, the scaling relations between γ and α, i.e., the scaling relations between γ, d, and β in Eqs (18), (20), and (22) are also validated. This implies that the distance exponent γ can vary according to the population sizes of origin and destination sites, even in the same population landscape. Here we remark that a recent empirical study showed that the origin and destination populations affect the travel patterns, whereas the distance exponent has been assumed to be the same irrespective of the populations [53].

Fig 3

Numerical validation of the expanded forms of the radiation-on-landscape (RoL) model.

The expanded forms of the rescaled travel probabilities in Eqs (15), (19), and (21) (solid curves) are tested against the numerical results using p in Eq (23) calculated on the same population landscapes used in Fig 2 (symbols), for the values of β = 1.5 (top) and 3 (bottom), respectively. For clearer visualization, we show only the curves of the rescaled travel probability for the group of the smallest ms (left) and those for the group of the largest ms (right), but for all groups of m in each panel. The analytic results of the distance exponent γ in Eqs (18), (20), and (22) are also plotted by black dashed lines for comparison.

We remark that the number of pairs of highly populated sites is in general much lower than those of other cases, so that the corresponding curves of the rescaled travel probability tend to be more fluctuating or even apparently missing, e.g., in the case with the groups of large m and m for β = 3 in Fig 3(d). Except for this case, we generically observe clear scaling behaviors of the rescaled travel probability in the intermediate range of r. In addition, the curves are found to saturate to a constant for sufficiently small r, whereas for sufficiently large r, these curves converge to eventually approach the lower bound of the rescaled travel probability, . These findings can be explained by Eq (23): On the one hand, for sufficiently small r, s becomes negligible as there would be only few or even no sites in the surrounding area between i and j. Thus, the rescaled travel probability becomes independent of r as . On the other hand, if r becomes sufficiently large, s approaches the total population M, irrespective of m and m. This is why all curves converge and eventually approach the lower bound of the rescaled travel probability as .

Finally, we discuss the generic behavior of the distance exponent according to the population sizes of origin and destination sites. We first point out that the scaling relations in Eqs (18), (20), and (22) have been derived in the limiting cases of m and m. Therefore, these results cannot be simply applied to the scaling behaviors observed for the cases with intermediate ranges of m and m. For these cases, one can estimate the apparent distance exponent γ based on the assumption of the simple scaling form as for each group of pairs G. It is found that the value of γ appears to be continuously varying according to the origin population m for the smallest and the largest groups of m, as depicted in Fig 4. For example, when m ≪ s, the value of γ is ≈ 2α for m ≪ s, and then it continuously decreases as m increases. We also find the clear dependency of γ on the destination population m for a given m.

Fig 4

Behaviors of the apparent distance exponent γ according to the origin and destination populations.

We estimate the values of the apparent distance exponent γ, defined by Eq (24), from the numerical curves of the rescaled travel probability shown in Fig 3 for the values of β = 1.5 (a) and 3 (b), respectively. These values (symbols) are compared to analytic values for the limiting cases, i.e., γ = 2α for m, m ≪ s and γ = α for m ≪ s ≪ m, which are plotted by black dashed lines with gray shadows denoting γ ± σ. Here σ is determined using the standard deviation of the estimated .

Conclusion

Although two representative mobility models, i.e., gravity and radiation models, have been compared to each other against the empirical traffic data sets [17, 18, 28], the more fundamental connection between these models has been far from being fully understood. In order to study such a connection in a realistic population landscape, we first model the heterogeneous population landscape by assuming a fractal geometry of sites and the population at each site following a power-law distribution. Then the radiation model on such population landscapes, namely, the radiation-on-landscape (RoL) model, can be written in terms of the distance between two sites. By expanding the rescaled travel probability in the RoL model and comparing it to the gravity model, we derive the distance exponent in the gravity model as a function of the fractal dimension and the population exponent of the population distribution. We also find that this distance exponent can vary according to the population sizes of origin and destination sites. These analytic expectations are confirmed by numerical simulations on our population landscapes. Consequently, we could connect two representative mobility models, and more importantly, the origin of the distance exponent could be related to the properties of the heterogeneous population landscape as well as the population sizes of origin and destination sites. Therefore we can better understand the mechanism behind the traffic flows constrained by the travel distance. In particular, the effects of the populations of origins and destinations on the distance exponent can be empirically studied as a future work.

In our work we have assumed that the location and population of each site are fully uncorrelated with each other, while there might be some correlations between them in reality. One can study the effects of spatial correlations, e.g., by the positively correlated populations at close sites, on the traffic flows and their characteristic distance exponent. In addition, as for the functional form of the population distribution, one can adopt other functional forms than the power law, such as the log-normal distribution given by Gibrat’s law [46].

Finally, we remark that the mass term m in many mobility models has been used to denote the population at the site, while it can be interpreted as other sources of attraction of sites, e.g., each site’s traffic volume [17], economic indicator [23], communication volume [18], and citations [24]. Indeed, the diverse values of distance exponent have been observed according to the modes of transportation, geographical regions, and granularities [13]. Considering our above findings on the mass dependency of the traffic flows, it is of crucial importance to empirically and theoretically relate various observables attributed to the site for better understanding of the human mobility.