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.


I. INTRODUCTION
For understanding the mechanisms of human mobility [1][2][3], optimizing the mobility flows [4], and predicting the dynamics on mobility networks [5][6][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 i (m j ) denotes the population of site i (j) and r ij is the distance between sites i and j.The value of distance exponent γ is found to range from 0.5 to 3 for some data sets [13].This original gravity model and its variants have been applied to human mobility and transportation [11][12][13][14][15][16][17][18][19][20], international trade [21], scientific collaboration [22], and mobile phone communication [18,23] to name a few, mostly due to their simplicity.On the other hand, the gravity models have limitations such as unclear theoretical background for the distance exponent and the absence of universality regarding the exponent estimation.
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 i is the outgoing traffic from the site i and s ij is the total population within a circle centered at the site i with radius r ij .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 longdistance travels, despite some unresolved issues like relatively poor predictability on short-distance travels [17].
The variants of the radiation and intervening opportunities models, e.g., a population-weighted opportunities model [24] and a radiation model with an additional scaling exponent [25], 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,26].Here we raise a question: Beyond the comparison, can these radiation and gravity models be connected to each other?The possibility of such connection was briefly discussed in Refs.[11,27], where the surrounding population s ij was assumed to be proportional to r 2 ij , i.e., the case with the uniformly distributed population in space.This leads to the asymptotic value of γ = 4.However, the population landscape in reality is well-known to be highly heterogeneous rather than being uniform.Such a heterogeneous population landscape can be characterized by a fractal geometry of populated areas, or sites, and a power-law distribution of the population at each site.In this paper, we will show that the distance exponent γ can be derived from the radiation model on a heterogeneous population landscape characterized by a fractal dimension and by the power-law exponent of the population distribution.By doing so, we can connect the radiation and gravity models, and more importantly, we can understand the origin of the distance exponent in the gravity model in terms of the properties of the heterogeneous population landscape.Therefore we can better understand the mechanism behind the traffic flows constrained by the travel distance.
The fractal geometry and the power-law distribution of populations are well-known characteristics of human settlement, e.g., cities and towns.Since the monumental paper on the fractal geometry by Mandelbrot [28], this concept has been applied to various topics, including the fractal structures of cities [29][30][31] and their growth patterns [32][33][34].In particular, the fractal dimension of cities is found to range from 1.3 to 1.9 [30,31].Next, 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 [35][36][37].The power-law exponent of the population distribution is found to have the value ranging from 1.7 to 3 [9,35,37,38].Despite the ongoing debate on whether populations are characterized by a power-law or a log-normal distribution [39][40][41], the power-law distribution of populations can be still a reasonable assumption for model studies.
Our paper is organized as follows: In Sec.II we introduce our setting for the heterogeneous population landscape.In Sec.III the radiation model on the heterogeneous population landscape is related to the gravity model to derive the distance exponent in terms of the fractal dimension of the landscape and the power-law exponent of the population distribution, which is also numerically confirmed.Finally, we conclude our work in Sec.IV.

II. MODELING POPULATION LANDSCAPE
We describe how to generate the heterogeneous population landscape characterized by a fractal geometry and a power-law distribution of populations.For this, we first generate a set of sites in a two-dimensional space with a fractal dimension d f .Then we assign to each site i the population m i independently drawn from P (m) ∼ m −β with an exponent β > 1, which will be called the population exponent.Note that the geometry of the sites can be implemented irrespective of the functional form of the population distribution.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 [42], 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 one site in the first layer and a circle centered at this site with radius R. Within this circle, η > 1 sites are randomly created in the second 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 η L−1 sites in the Lth layer.In our work, we consider the set of sites only in the last layer to find its fractal dimension as [43] Once the set of N = η L−1 sites with a fractal geometry is generated, N values are independently drawn from a population distribution P (m) to be assigned to each site.
As for the population distribution we adopt the powerlaw distribution with the population exponent β: where C is a normalization constant.We set m min = 100 and m max = 10 7 .Figure 1(b) shows an example of the generated population landscape in the two-dimensional space using d f = 1.8, η = 2, L = 14, and β = 3, where the value of λ is determined by Eq. ( 3).The height in the vertical axis indicates the population assigned to each site.
Finally, although there exist many other modeling approaches for generating heterogeneous population landscapes [32,[44][45][46], we have adopted the Soneira-Peebles model for the fractal geometry, mostly because the implementation of this model is efficient and scalable.

III. CONNECTING RADIATION AND GRAVITY MODELS A. Scaling behavior of surrounding populations
The connection between radiation and gravity models can be made by the observation that the surrounding population s ij of the radiation model in Eq. ( 2) must be correlated with the distance r ij of the gravity model in Eq. ( 1).We remind that the surrounding population s ij is defined as the total population within a circle centered at the site i with radius r ij .As a limiting case, if the population is distributed uniformly in space, one can get s ∝ πr 2 [11].However, the relation between s ij and r ij can be nontrivial in general, particularly for the heterogeneously distributed population modeled in Sec.II.Precisely, let us denote by Λ ij the set of sites within a circle centered at the site i with radius r ij .Here Λ ij does not include sites i and j, and the number of sites in Λ ij is denoted by n ij .Then the surrounding population is written as Here m l s are independently drawn from the population distribution P (m) ∼ m −β in Eq. ( 4).For a sufficiently small value of β, the summation in Eq. ( 5) is dominated by the maximum value, i.e., max{m l } l∈Λij , which is proportional to n 1/(β−1) ij [47,48].On the other hand, for a sufficiently large value of β, all terms in Eq. ( 5) contribute to the summation, leading to s ij m l n ij with m l denoting the average of m l .Finally, by considering the fractal geometry of the sites, we can relate n ij to r ij : In sum, we obtain the scaling relation between s ij and r ij , by introducing the scaling exponent α as with where the crossover value for β has been set as 2 considering the continuity between two regimes of β.
We numerically test the scaling relation between α, d f , and β in Eq. ( 8), using our model for the heterogeneous population landscapes that we have described in Sec.II.For this, we generate the population landscapes using the same parameter set of d f = 1.8, η = 2, R = 1, and L = 14, then assign to the sites the populations drawn from P (m) ∼ m −β for several values of β.After logbinning the measured s ij and r ij for all possible pairs of i and j, we estimate the scaling exponent, denoted by α sim .For example, as shown in Fig. 2, the value of α sim turns out to be comparable to the theoretical value of α from Eq. ( 8): α sim = 2.23(2) and α = 2.25 for β = 1.8, while α sim = 1.74(2) and α = 1.8 for β = 3.Note that the scaling exponent α relating s ij and r ij was empirically estimated for each county in the United States to find that the distribution of α is peaked around at α ≈ 2.53, see Fig. S1 of Ref. [11].

B. Derivation of the distance exponent
Now we connect the radiation model on the heterogeneous population landscape to the gravity model.For this, we begin with a radiation model for finite systems [17]: where M is the total population in the system, i.e., M = i m i .p ij denotes the probability of an individual from a site i to travel to another site j, satisfying j p ij = 1.By assuming that m i M , m i s ij , and m j s ij , one can get the approximated equation: By plugging s ij in Eq. ( 7) into the above equation, we get the normalized travel probability as a function of r ij : As a result, the travel probability from a site i to another site j normalized by their populations is expressed as a power function of the distance r ij , similar to the gravity model in Eq. ( 1): By comparing the distance dependence of the radiation and gravity models, we finally obtain the distance exponent γ as a function of the fractal dimension d f and the population exponent β: Our theoretical expectation in Eq. ( 13) implies that the underlying configuration of the population landscape can account for the scaling between mobility flows and distances 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 distance dependence of the traffic flows.Let us consider a jobseeking situation as in the original radiation model [11].Since the number of cities is proportional to r d f , a higherdimensional geometry with a larger d f 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 f 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 each 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 γ.On the other hand, 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  13) is denoted by black dashed lines, and the estimated γsim from the simulations by red circles with error bars depicting the interval of regression coefficients.The 95% confidence interval for γ(d f,sim , β), denoted as violet shade, is obtained from the estimated fractal dimension d f,sim of the generated population landscapes.

C. Numerical validation
By performing numerical simulations, we validate the scaling relation in Eq. ( 13) for the distance exponent γ as a function of the fractal dimension d f and the population exponent β.We use the same heterogeneous population landscapes generated in Subsec.III A. On these population landscapes, we directly calculate the travel probability p ij in Eq. ( 9) for all possible pairs of sites i and j [49].Then we plot the normalized travel probability, i.e., p ij /(m i m j ), against the distance r ij .Then, after log-binning the scatter plot of (r ij , p ij /(m i m j )) for all possible pairs of i and j, we estimate the scaling exponent, denoted by γ sim .The value of γ sim turns out to be comparable to the theoretical value of γ from Eq. ( 13), e.g., we find γ sim = 4.50(2) and γ = 4.5 for β = 1.8, while γ sim = 3.37(2) and γ = 3.6 for β = 3, see Fig. 3.
We remark that the values of γ sim have been estimated in the scaling regime of 10 −2 < r ij < 1 for both values of β as depicted in Fig. 3.That is, we find the lower and the upper cutoff scales for the scaling regime.The lower cutoff can be related to the smallest length scale that is R/λ L = Rη −L/d f ≈ 10 −2 for the parameter values used.The upper cutoff is related to the largest length scale, which is trivially R = 1.
Next, we test the reliability of our results in Eq. ( 13) for the wider range of parameter values, i.e., 1.5 ≤ β ≤ 3 for several cases with (d f , η, L) = (1.5, 2, 14), (1.8, 2, 14), and (1.8, 5, 6), respectively, as shown in Fig. 4. The estimated values of γ sim from the numerical simulations turn out to be comparable with the theoretical values of γ(d f , β) in Eq. ( 13), but with some systematic deviations in the vicinity of β = 1.5 as well as around β = 2.The deviations in the vicinity of β = 1.5 can be partly understood as the generic finite-size effect regarding the powerlaw distribution.That is, the finite number of samples drawn from the power-law distribution might result in the natural cutoff in the distribution of the samples especially when the original distribution has a too heavy tail.In addition, the existence of sites with extremely large populations can undermine our assumptions, e.g., m i M , for deriving the scaling relation.The deviations of γ sim observed around β = 2 might be due to the deviations of the fractal dimension of the generated population landscape from the theoretical value in Eq. ( 3).In order to test this hypothesis, we measure the fractal dimension of the generated population landscape and denote it by d f,sim .Based on the standard deviation of d f,sim , we obtain the 95% confidence interval of γ(d f,sim , β), which is depicted by shaded regions in Fig. 4. The deviations of the fractal dimension are found to account for the deviations of γ sim only to some extent.

IV. CONCLUSIONS
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,26], the more fundamental connection between these models has been far from being fully understood.In order to study such a connection in a more realistic population landscape, we first model the heterogeneous population landscape by generating a fractal geometry of sites using the Soneira-Peebles model and then by assigning to each site a population independently drawn from a power-law distribution.Then the radiation model on this population landscape is approximated to be written in terms of the distance between the origin and the destination of travels.This approximated equation is directly compared to the gravity model to derive the distance exponent in the gravity model as a function of the fractal dimension and the population exponent of the population distribution.Numerical simulations show a good agreement with our theoretical expectation.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.Therefore we can better understand the mechanism behind the traffic flows constrained by the travel distance.
So far we have assumed that the location and population of each site are fully uncorrelated with each other.However, in reality there might be some correlations between them.For example, one can consider the spatial correlation such that populations at geographically close sites are positively correlated, which must have an effect on the distance dependence of traffic flows.In an extremely positive-correlated case, the population landscape could show a highly centralized structure, where the population density can be described as a decaying function of the distance from the central site.The effect of the spatial correlation on the distance dependence of traffic flows can be studied as a future work.In addition, as for the functional form of the population distribution, one can adopt other functional forms than the power law, such as the lognormal distribution given by Gibrat's law [39].
Finally, we remark that the mass term m i 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 [21], communication volume [18], and citations [22].Indeed, the diverse values of distance exponent have been observed according to the modes of transportation, geographical regions, and granularities [13].In this respect, it is of crucial importance to empirically and theoretically relate various observables attributed to the site for better understanding of the human mobility.

FIG. 1 .
FIG. 1.(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 second layer are randomly placed within the circle with radius R. Similarly, the sites at the third layer are randomly placed within the circles with radius R/λ.(b) An example of the generated population landscape using the Soneira-Peebles model with d f = 1.8, η = 2, and L = 14, 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.

FIG. 2 .
FIG. 2. Numerical test of the scaling relation in Eq. (8) between the surrounding population sij for a pair of sites i and j and their distance rij for the values of β = 1.8(a) and of β = 3 (b), on the same fractal geometry of sites generated using Soneira-Peebles model with d f = 1.8, η = 2, R = 1, and L = 14.For the gray-colored heat maps, the darker color implies more pairs of sites around the point (rij, sij), obtained from 100 different population landscapes.From the log-binned curve (red circles), we estimate the scaling exponent relating sij to rij as αsim = 2.23(2) for β = 1.8 and αsim = 1.74(2) for β = 3, respectively, which is comparable to the theoretical value of α (black line) from Eq. (8).

FIG. 3 .
FIG. 3. Numerical tests of the scaling relation in Eq. (11) between the normalized travel probability pij/(mimj) for a pair of sites i and j and their distance rij for the values of β = 1.8(a) and of β = 3 (b), on the same population landscapes used in Fig.2.By using the same log-binning method as in Fig.2, we estimate the distance exponent as γsim = 4.50(2) for β = 1.8 and γsim = 3.37(2) for β = 3, respectively, which is comparable to the theoretical value of γ (black line) from Eq. (13).

FIG. 4 .
FIG. 4. Comparison of theoretical and simulation results on the distance exponent γ for a wide range of 1.5 ≤ β ≤ 3 and for the cases with (d f , η, L) = (1.5, 2, 14) (a), (1.8, 2, 14) (b), and (1.8, 5, 6) (c).The theoretical curve of γ(d f , β) in Eq. (13) is denoted by black dashed lines, and the estimated γsim from the simulations by red circles with error bars depicting the interval of regression coefficients.The 95% confidence interval for γ(d f,sim , β), denoted as violet shade, is obtained from the estimated fractal dimension d f,sim of the generated population landscapes.