Identifying and modeling the structural discontinuities of human interactions

The idea of a hierarchical spatial organization of society lies at the core of seminal theories in human geography that have strongly influenced our understanding of social organization. Along the same line, the recent availability of large-scale human mobility and communication data has offered novel quantitative insights hinting at a strong geographical confinement of human interactions within neighboring regions, extending to local levels within countries. However, models of human interaction largely ignore this effect. Here, we analyze several country-wide networks of telephone calls - both, mobile and landline - and in either case uncover a systematic decrease of communication induced by borders which we identify as the missing variable in state-of-the-art models. Using this empirical evidence, we propose an alternative modeling framework that naturally stylizes the damping effect of borders. We show that this new notion substantially improves the predictive power of widely used interaction models. This increases our ability to understand, model and predict social activities and to plan the development of infrastructures across multiple scales.

The idea of a hierarchical spatial organization of society lies at the core of seminal theories in human geography that have strongly influenced our understanding of social organization.In the same line, the recent availability of large-scale human mobility and communication data has offered novel quantitative insights hinting at a strong geographical confinement of human interactions within neighboring regions, extending to local levels within countries.However, models of human interaction largely ignore this effect.Here, we analyze several country-wide networks of telephone calls and uncover a systematic decrease of communication induced by borders which we identify as the missing variable in state-of-the-art models.Using this empirical evidence, we propose an alternative modeling framework that naturally stylize the damping effect of borders.We show that this new notion substantially improves the predictive power of widely used interaction models, thus increasing our ability to predict social activities and to plan the development of infrastructures across multiple scales.
Globalization has led us to believe that our world is becoming borderless and deterritorialized.The rise of novel information technologies has even prompted the forecast of the "death of distance" [5].However, even a most basic organization of society requires categories, compartments and borders to maintain order [6].Confinement of human interactions to limited spatial areas is the key message of the long-standing hypothesis of Central Place Theory (CPT) [1,2] which posits the existence of regular spatial patterns in regional human organization.In short, CPT assumes the existence of a "hierarchy" of regions that aims to explain the number, size and locations of human settlements with spatio-economic arguments.Despite its highly simplifying geometric assumptions (Supplementary Information), empirical evidence for CPT's main prerequisite of systematically limited human interactions has been collected in a number of recent studies on massive interaction networks which have indeed observed a substantial impact of political or socio-economic boundaries on human interactions [3,4,[7][8][9][10][11].Typically, if we construct regions by clustering those locations that have strong interactions with each other, we divide countries into contiguous geographical regions with separating boundaries often following surprisingly close existing administrative boundaries.These findings point towards spatial regularities in human organization, prompting us to ask: is there an underlying, rigorously quantifiable, principle behind these patterns?If so, can we exploit this principle to develop better models of human interaction?

Quantifying the inhibitory effect of borders on human interaction
To quantify the hypothesized effect of hierarchical organization on human interactions, we first define consistent nested regional partitions by recursively applying a community detection algorithm to country-wide phone call networks from the United Kingdom, Portugal, France, Ivory Coast and an anonymous country, Country X (Methods).Partitions resulting from this algorithm reflect the communities defined by underlying social interactions, and, contrary to official administrative boundaries, are independent of country-specific historical or political contexts [4] (Supplementary Information).The resulting partition consists of three levels, L 1 , L 2 , and L 3 , that have a natural interpretation: the whole country is divided into L 1 -level regions (regional scale), which are divided into L 2 -level regions (county scale) which in turn split into L 3 -level regions (city scale) composed of several "elementary" locations (cell phone tower or exchange area), Fig. 1.The number of levels is not imposed, but for all countries the process naturally stops subdividing regions at the city scale.To our surprise, we find that, although no spatial constraints are applied, communities consist of contiguous locations at all levels (an observation that had previously been reported just for L 1 -level regions [3,4]), and are strikingly similar to administrative regions as highlighted by comparison with random partitions (Supplementary Table 1).
Several insights that we first derive from these hierarchical partitions of empirical networks are in line with CPT.The L 3 regions have typical spatial extension of a town with its neighborhood [12] (between 15 km to 23 km, depending on the country); similarly, L 2 and L 1 conform to the scale of districts and regions (Fig. 2a).The distribution P (n) of the number n of L 3 communities inside a L 2 community is strongly peaked (around 6) providing quantitative confirmation to the main hypothesis of regular spatial organization of CPT (Supplementary Information).Figure 2f shows that distribution, as well as the distribution of the number of L 2 communities within an L 1 community (#L2/#L1), for the UK.We observe similar peaks in all other countries ( Figs.S1-S4).Following the idea that borders inhibit human interaction, we introduce the notion of hierarchical distance to characterize their impact on communication flows (Fig. 1).Two locations i and j are at a hierarchical distance h ij = 1 if they are in the same L 3 region, h ij = 2 if they are in different L 3 regions, but in the same L 2 region, h ij = 3 if they are in different L 2 regions, but in the same L 1 region and h ij = 4 if they are in different L 1 regions.In other words, the hierarchical distance corresponds to the number of different types of borders separating two locations.This metric only contains limited information about the spatial structure of the regions and is only partly correlated with geographic distance: two locations that are close in terms of geographic distance can still be situated in two distinct L 1 regions and hence far from each other in terms of hierarchical distance.Thus, the hierarchical distance is not a mere discretization of geographical distance, but encodes a qualitatively different, socio-economic notion of distance.To understand the impact of borders on human interaction on each hierarchical level, we define and measure the following damping parameters = j:hij =h T ij is the total duration of calls from location i to all locations at hierarchical distance h.Defining the weight of node i, w i = j T ij , as the total duration of calls originating from node i (including self-loops), W (h) i = j:hij =h w j is the total duration of calls originating from locations at a hierarchical distance h from location i.The ratio measures the relative strength of communication between location i and the locations at a hierarchical distance h from it.In particular, this ratio corresponds to the amount of communication sent to all locations at hierarchical distance h per unit of communication produced there.The damping value q (h) i hence measures the relative importance of locations at hierarchical distance h + 1 compared to those at hierarchical distance h from i.For example, q (h) i = 1 for all h means that communication from i is independent of the hierarchical distance, because there is no damping in the amount of communication sent per unit of communication produced as the hierarchical distance increases.Figure 2k shows the distributions of the damping values in the UK, all well peaked around a strikingly similar mean value which does not substantially change with the hierarchy level (h = 1, 2, 3), Table I.Similar observations are made for all studied countries (Figs.S1-S4).This finding -signing a structural discontinuity of human interactions -and its consequences on modeling, see below, is our main discovery.It means that the damping effect of a boundary is approximately the same irrespective of the level h and origin location i, i.e. q (h) i q.If the probability for two people who live in the same L 3 region to communicate is p 0 , it will be qp 0 for people living in different L 3 regions but in the same L 2 region, q 2 p 0 if they live in the same L 1 but different L 2 region and q 3 p 0 if they live in different L 1 regions, Fig. 3b.

Why and how standard models fail
Using a standard quality of fit statistic (Methods) and comparing distributions of high level per low level regions, we tested to which extent the most widely used models, namely gravity [13] and radiation [14], commit a systematic bias by failing to account for the observed boundary effects.To this end, we compute the communication networks predicted by these models as well as the corresponding partitions resulting from the community detection algorithm (Methods).As previously demonstrated elsewhere [14], the gravity model strongly underestimates and fails to predict high-range flows, i.e. flows between locations where the number of calls is high (Fig. 4a andFigs.S5a to S8a).This certainly explains why the gravity model generates less and larger L 1 regions whose subdivisions do not follow the narrow distributions observed in the data (Fig. 2b and 2g).The damping value predicted by the gravity model is otherwise well peaked, although its average values vary from one h level to another (Table I).In contrast, the radiation model overestimates long-range flows (Fig. 4b), resulting in more and smaller L 1 regions (Fig. 2c).Therefore, the distribution of L 3 within L 2 regions is, although well-peaked, shifted to the left (Fig. 2h).The distribution of damping values in the radiation model is moreover strongly spread out (Fig. 2m and Table I), and does not reproduce the existence of a single typical damping parameter.Similar systematic biases of gravity and radiation models become evident if we measure the probability P dist (d) of a call between locations at distance d (Fig. S9).

Accounting for strong border effects with the Hierarchy model
The two most commonly used models thus fail to reproduce the boundary effect.By design, a model taking into account the observed hierarchical structures by assuming a constant damping value q, would overcome this issue (Fig. 3b).Consider the minimal model in the stylized form T ij ∝ N ij q hij , where N ij represents the potential pairs of contacts between two distinct locations i and j and q hij the probability for two people from these locations to communicate.This model would implement highly discretized hierarchical distances instead of considering a continuum of geographical distances.Similarly to the gravity model, N ij can be taken as proportional to the weight w i and w j of both origin and destination locations.We therefore propose a simple hierarchy model that predicts an interaction strength as where 0 < q < 1 is a parameter to be determined and C i are local normalization factors ensuring w Hier i = w i .This normalization also ensures that the damping values are constants, q Hier i = q (Supplementary Information).The bestfit values of q are very close to the observed values (Table I and Table S2) and robust to small variation (Fig. S10).These values slightly depend on the country, varying between 0.10 and 0.25, reflecting differences in the structural properties of the studied networks.The hierarchy model reproduces almost perfectly the nested structure of regions (Figs.2d and 2i).To our surprise, the hierarchy model also outperforms the state-of-the-art models in terms of goodness of fit measures (Table II and Extended Table 3).In particular, it estimates high-range flows with a greater accuracy than the radiation or gravity models, as can be seen on the top right corners of Fig. 4c and Figs.S5c to S8c, where the markers are typically closer to the equality line than in state-of-the-art models.

DISCUSSION
Focusing on flows between locations at specific hierarchical distance, the goodness of prediction of the different models is informative to understand why the hierarchy model outperforms the others.The radiation model overestimates the flows at h = 1, 2, the corresponding markers being above the equality line in Fig. 4f and Figs.S5f to S8f resulting in an overall overestimation quantified by values of R h=1,2 (Methods) greater than 1 (Table III), and underestimates those at h = 3, 4. On the contrary, the gravity model underestimates the flows at h = 1, 2 and overestimates those at h = 3, 4 (Table III, Fig. 4e, Figs.S5e to S8e).The hierarchy model produces more balanced predictions (R h closer to 1) and thus outperforms existing models.
The hierarchy model requires the knowledge of the communication flows in order to determine the three hierarchical levels each location belongs to.However, it can also be applied in the absence of communication data, using the administrative boundaries and a general damping value q = 0.2 which matches robustly all countries (Extended Data Fig. 10).The resulting hierarchy-admin model based on this administrative partition, is parameter-free and yet it provides similar or sometimes better estimates than the gravity model in terms of communication flow (Fig. 4k, Figs.S5k to S8k: see in particular the case of high-range flow in Portugal and Ivory Coast) or benchmark measures (Table II and Extended Data Table 3).We also tested different constraint conditions and deterrence functions f in the hierarchy model T Hier ij = C i w α i w β j f (h ij ).We compared them to multiple variations of the gravity and radiation models which are widely outperformed by hierarchy models (Supplementary Information).
In summary, we first defined communication flows-induced boundaries by applying standard community detection methods on large-scale human interaction networks and found that these networks have a nested structure reflecting historic, socio-political borders which can be related to the structure predicted by CPT.We introduced the notion of damping parameter, representing the normalized ratio of interactions between locations at different hierarchical distances, to quantify the inhibiting effect of boundaries.Surprisingly, the distributions of damping parameters are well-peaked and largely independent of the hierarchical level, revealing a structural discontinuity effect in each studied country.We further showed that current models of human interaction, based only on population and/or geographical distance, cannot correctly reproduce the characteristic hierarchical structure of interaction networks.We proposed a simple model based on the discrete hierarchical distance that outperforms the state-of-the-art models of human interaction in a number of different countries, demonstrating its general applicability and emphasizing the impact of the borders on human interactions.The development of more sophisticated models combining both geographic and socio-political information will further boost our ability to understand and reproduce the structure of social systems.

Telephone call data
We consider several country-wide data sets of telephone calls, including the four European countries of the UK, France, Portugal and an anonymized Country X, and Ivory Coast.All data sets comprise mobile phone data with the exception of landline calls in the UK.Data was provided by single phone providers with possibly heterogeneous coverage over the respective countries -we have no information on local market shares and on resulting possible inhomogeneities in spatial coverage.Specific details of the different datasets are provided in Extended Data Table 1, all of them gathering millions of users making billions of calls during time frames ranging from 1 to 15 months.We construct interaction networks between different locations of a country based on the aggregated duration of calls having origin in the first and destination in the second location.This process generates weighted directed networks in which loop edges from locations to themselves are also considered, and where the link weight T ij between a location i and location j is defined as the total duration (or, in case of Country X, total number) of calls from location i to location j.The nodes of the network are the locations, corresponding to exchange areas or cell towers areas as reported in Extended Data Table 2.In all datasets, the users are attached to the actual locations where the calls occur (i.e.not necessarily their residential locations).

Network partitioning
A recently developed algorithm for community detection, referred to as "Combo" [3,4,8] is applied to the extracted communication networks to detect communities of highly connected locations.The method follows a standard modularity optimization approach [15,16], scoring the edges of the networks according to their relative strength compared to a null-model based on the weight of the nodes they connect and aiming at maximizing the cumulative score inside the communities.Given a partition of the nodes in a set of clusters c i , the modularity score Q is given by where T ij is the weight of the link between node i and node j, w i = j T ij is the weight of node i and W = i w i /2 is the total weight of the network.While the outcome of partitioning networks is in general not qualitatively dependent on the particular algorithm used, the Combo algorithm has the ability to consistently provide the best results in terms of modularity score compared to other algorithms [17].The modularity optimization approach yields communities whose size and properties are only based on the informations of the links' weights.See [18] for a more explicit interpretation of the modularity, its properties and limits.
Applying the Combo algorithm yields a first partition of the network into communities, further referred to as "level 1" or "L 1 " partition.To obtain the substructure of these communities, we iteratively apply the Combo algorithm on each L 1 community, thus creating a"level 2" or "L 2 " community partition, and and then again on each L 2 community, thus creating a "level 3" or "L 3 " community partition.We find that most of the L 1 and L 2 communities display a clear substructure with high values of internal modularity scores, typically around 0.4 and 0.7 (Supplementary Table I).The resulting communities consists in geographically cohesive regions, which can seem surprising since the algorithm uses only the networks topology and no geographical information, such as the distance between the nodes (Supplementary).This cohesiveness is also linked to the spatial scale of the studied network: we would not expect any contiguous communities if that analysis was done at a city scale, where the movements and communications of individuals are more evenly distributed in space.

Interactions models and goodness measures
The radiation model is a parameter-free model recently introduced in the context of migration patterns [14].Given the geographic distance d ij between two locations i and j, the model predicts that the flow of individuals moves T ij between these two distinct locations will depend on the population at the origin, the population at the destination and on the population s ij within the circle of radius d ij centered on the origin location i. Applied to our case (using the total communication w i at location i as a proxy for its population), the radiation model is written as where s ij = k, 0<d ik <dij w k is the total amount of communication originating from locations at a distance shorter than d ij from location i and C i is a normalization factor ensuring that the predicted total activity of each node is the same than the actual one, i.e. j =i T Radiation ij = j =i T ij .The model is otherwise parameter-free.The gravity model is one of the oldest models describing human mobility and interaction, formulated in analogy to Newton's law of gravity.The classical form predicts that the interaction strength between two distinct locations varies with the distance between them according to a power law: where C is a global normalization constant ensuring that i, j =i T Gravity ij = i, j =i T ij and α, β, and γ are parameters to fit.
We also computed the generalized version of the radiation model proposed in [19], as well as different versions of the gravity and hierarchy models, comparing the results obtained using a power law or exponential deterrence function (Supplementary Information).All parameters in these models were estimated through a regression analysis minimizing the deviance E [20], a measure based on the log-likelihood of model compared to a sturated model that can be interpreted as a generalization of the residual sum of squares R 2 .
We also quantify the fits between communication networks and models through different benchmark measures, namely the Dice distance D, the Sorensen distance S, and the cosine distance C defined by: These three benchmark measures cover most families of distance measures [21], which allows us to ensure that our findings are stable with respect to the distance measure used.They all vary between 0 and 1 and the lower they are, the more similar the model is to the original data.Finally, we also computed the correlation corr between each model and the data defined by which is a measure of similarity varying between -1 and 1 (the closer to 1, the higher the similarity).

Over-and underestimation measure
In order to determine whether a given subset of links are over-or underestimated by the models, we define for any given set E of links, the following ratio: where we use the notation A for the data and M for the model.Values of R E larger (resp.smaller) than 1 hence correspond to an overestimation (resp.underestimation) of the model.The measure R E provides an aggregated • The authors declare that they have no competing financial interests.Correspondence and requests for materials should be addressed to S.G. (email: sebgrauwin@gmail.com).Partitioning of a country based on telephone call networks.Hierarchical distances between two locations are defined through three regional levels -either administrative ones or those found by applying iterative community detection on human interaction networks.Two distinct locations are at a hierarchical distance h = 1 if they are in the same L3 region, h = 2 if they are in different L3 regions but in the same L2 region, h = 3 if they are in different L2 regions but in the same L1 region and h = 4 if they are in different L1 regions.Note that a higher hierarchical distance does not necessarily correspond to higher geographical distance.In the classic gravity model, the probability p that two people communicate is a continuous (e.g.exponential) function of the distance between them.b, In our hierarchy model, that probability is a discontinuous function induced by the assumption of a constant damping value q independent of the point of origin and the hierarchy level h.In both cases, the left panel shows in grayscale the probability of communication from a given point in space in a schematic country, partitioned in three regional levels with the same color coding as Fig 1.
The link between the borders and the deterrence function is clearly apparent in the second case.Gray markers are scatter plots for each pair of locations.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.S1) Properties of the data sets.Country-wide telephone data sets are provided by single telephone operators, covering different time frames, with different numbers of phones, calls, total call durations and on various spatial resolutions.The abbreviations bn. and m. stand for billion and million, respectively.Resolution numbers are given as approximate values.These locations constitute the nodes of the corresponding telephone call networks, while the sum of durations of calls between locations span their weighted links.The last columns report the percentage of non-zero links between pairs of nodes in the extracted network and whether to not that network is directed.The durations of calls are unknown in the case of Country X.All datasets corresponds to mobile phone network except for UK, where the dataset corresponds to a landline network.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.S2) Hierarchical properties of spatial organization from human interactions in France.a-e, Maps of L1 communities in telephone call networks detected from data and from various interaction models.Black lines correspond to the administrative partitioning of the 22 NUTS1 regions of France, colored areas to regions detected by a community detection algorithm applied to (a) the data, and to the (b) gravity, (c) radiation, (d) hierarchy, and (e) administrative models.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.a-e, Maps of L1 communities in telephone call networks detected from data and from various interaction models.Black lines correspond to the administrative partitioning of the 19 administrative regions of Ivory Coast, colored areas to regions detected by a community detection algorithm applied to (a) the data, and to the (b) gravity, (c) radiation, (d) hierarchy, and (e) administrative models.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.administrative models.Gray markers are scatter plots for each pair of locations.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.The proportion P dist (r) of communication occurring between two locations at a distance r (in km) from each other, is measured in the data and in the models.The radiation model is characterized by a lower than actual proportion of communication between distant (more than 100km) locations up to two orders of magnitude.It also presents a higher than actual proportion of communication between close (less than 10km with a peak between 0.5 and 5 km depending on the country) locations up to one order of magnitude.The gravity model presents in all countries a higher than actual proportion of communication between very close locations (100m-1km).In general, it also overestimates the low-range fluxes by 2 to 3 orders of magnitude and slightly underestimates the top-range fluxes.The hierarchy model fits almost perfectly at low distances (less than 10km).Depending on the country, it only deviates slightly from the data at top-range fluxes or estimates them properly.The fit of the hierarchy-admin model depends strongly on the country, but qualitatively comparable to the gravity and hierarchy models.The dashed lines and circle markers show the 'optimal' values of q (reported in Supplementary Table 2) minimizing the deviance E. In all countries, either for the hierarchy or hierarchy-admin model, these optimal values of q are also close to those minimizing the Dice (D), Sorensen (S) and Cosine (C) distances and maximizing the correlation (corr) between the data and the models.These optimal values are also stable, in the sense that close values of q (roughly between 0.1 and 0.3) still provide benchmark measures close to their optimum.The dotted lines in the case of the Hierarchy-Admin model indicate the damping value q = 0.2 matching robustly all countries.

Partition overlap measures
We use three classical measures of clustering similarity to quantify partition overlaps, i.e. of how well two different partitions of the same set of locations match: Rand's criterion R [? ], Jaccard index J [? ] and the Fowlkes and Mallows index F [? ].Consider two partitions C and C of a set of n nodes and note • n 11 the number of pairs of nodes in the same community both in C and C ; • n 01 the number of pairs of nodes not in the same community in C but in the same community in C ; • n 10 the number of pairs of nodes in the same community in C but not in the same community in C ; • n 00 the number of pairs of nodes in different communities both in C and C .
All types of pairs being taken into account, we have n 11 + n 01 + n 10 + n 00 = n(n − 1)/2.The R, J and F indices are then defined by: which are different ways of quantifying how well the partitions match pairs of nodes.The values of R, J , F lie between 0 and 1, and values close to 1 indicate a perfect match between the two compared partitions.However, even for the case of two completely unrelated clusterings, all indices are in general strictly larger than zero, more so for R [? ].Therefore, to have a baseline, we calculated the average indices over 10000 random shufflings of the nodes's clusters, denoted by Rr , Jr and Fr (see Table V).

Results
The partitions obtain for levels L 1 , L 2 and L 3 correspond in general to geographically cohesive regions, and are rather similar to administrative regions in number and size.As previously noted in [? ??] (from which the partitions resulting from the Combo algorithm for UK and Portugal are reproduced), these results may come as a surprise regarding that the modularity approach of the Combo algorithm has no spatial constraint nor does it impose any restriction on the number of communities.
The partition overlap measures given in Table V, together with p-values assessing the statistical significance of the partition with respect to a null model in which communities of cell towers were randomly reshuffled, quantitatively confirm the similarity between the administrative partitions (we choose to refer to the european Nomenclature of Territorial Units for Statistics -or NUTS -standard for the european countries) and the Combo partition.For example, the L 1 UK partitioning shows values of R = 0.954 with a baseline Rr = 0.825, J = 0.618 with a baseline Jr = 0.049, and F = 0.769 with a baseline Fr = 0.096 -with all significance measures (p < 10 −4 ) indicating a good match between the administrative and the Combo partitions.Going a step further and comparing the L 2 and L 3 communities with NUTS regions of corresponding levels, the match between the Combo and administrative partitions is still good, as indicated by significantly higher than average overlap measures.The L 2 UK partitioning hence shows values of R = 0.972 with a baseline Rr = 0.957, J = 0.222 with a baseline Jr = 0.007, and F = 0.439 with a baseline Fr = 0.018 and the L 3 UK partitioning shows values of R = 0.988 with a baseline Rr = 0.986, J = 0.099 with a baseline Jr = 0.001, and F = 0.292 with a baseline Fr = 0.004.Note that while all these values stay significant, the differences between the Combo / administrative overlap values and the random / administrative overlap values decrease when we look at more fine-grained partitions.This has to be expected since the deviations between the Combo and administrative partition can only increase with the level of partition as the deviations at a given level automatically impact the successive levels.Similar results can be drawn from the other countries, see Table V.
Table V also display the modularity scores Q * Combo and Q * of f of the Combo and administrative partitions at the different levels.To be consistent with the procedure of the Combo algorithm, which builds the level n + 1 subpartition by decomposing each community of the level n partition, these modularity scores indicate in case of L 2 (resp.L 3 ) partitions the average modularity score of each subpartition defined with respect to a subnetwork inside each corresponding L 1 (resp.L 2 ) community.Our measures show that the Combo partition always has a better modularity score than the administrative partition.The modularity scores also decrease when we look at higher level partitions, indicating that the L 1 communities are the most relevant.

Interpretation
The similarity between the regions emerging from the communication network through the Combo procedure and the administrative boundaries can be interpreted as a natural evidence towards the latter's validity [? ?].In view of the deviation between the two partitions, Combo partitioning appears to be more aligned with human interactions, as measured by the modularity score, suggesting that the border of administrative regions sometimes deviates from the underlying reality of interactions.Most interestingly, the partitioning created by our approach provides a unified hierarchical framework to compare the geographical structure of human interactions in different countries, which present an important alternative to the administrative boundaries whose shape and number depend substantially on the historical and political context of each country as well as particular local regional delineation policy.The "Network" column indicates the levels of the administrative partition and of the community partition that are compared [? ].
Columns N of f and N combo respectively refer to the number of NUTS regions and the number of communities found by the community detection algorithm at the considered level.Columns Q * of f and Q * Combo indicate the average modularity score of all the administrative or Combo sub-partitions at the considered level.Rr, Jr and Fr give the baselines for Rands criterion R, the Jaccard index J and the Fowlkes-Mallows index F .The closer R, J or F to 1, the better the overlap of the detected communities with the administrative regions.The baseline values of the similarity indices are computed on 10000 random shuffling of the nodes's clusters.For the 3 used measures, none of these random shufflings is more similar to the administrative partition than the partition found by the community detection algorithm: the statistical significances of all the similarities have a p-value < 10 −4 .

Central Place Theory
Classical theories of urban planning have traditionally suggested economical and geographic laws to systematically determine the arrangement of towns and cities.In particular, the central place theory (CPT) developed by Christaller [? ?] seeks to explain the number, size and location of human settlements in an urban system.The basic assumption of the theory is that settlements function as 'central places' providing services to surrounding areas.The hierarchy of the cities is based on the range of goods and services they provide.Low order goods and services (groceries, bakeries, post offices) are present in all places, including small and large centers.Higher order goods and services (jewellery, large malls, universities) are only present in large centers, which are less numerous.These centers are supported by a large population, including its own and those of the surrounding smaller centers.The lowest settlements should form an hexagonal lattice, this being the most efficient regular pattern to serve areas without any overlap (in terms of radial distance and perimeter for a fixed area).Settlements of higher order (villages, cities) should then be regularly spaced on an hexagonal pattern of higher radius, with their centers placed on centers of hexagons of the lowest order.
Christaller defined K-hexagonal landscapes as arrangements where each higher order settlement is supported by K − 1 lower order settlements and itself.Christaller and later Lösch [? ] both developed arguments over which value of K is adapted to describe different situations.For example, K = 3 is suited for sharing local goods (marketing principle), K = 4 is suited for reducing cost of transport (traffic principle) while K = 7 -a case where each satellite depends only on one center -is suited for political stability (administrative principle).Christaller conceived these models as hierarchical, with all higher order places in the hexagon surrounded by lower order places to explain not only local but regional economics and spatialization of urban centers, the value of K possibly changing from level to level.
Filling the gap between CPT and reality, several distortions to the original model have been introduced over time to account for inhomogeneous population densities, resource locations, or specialization of cities.However, nowadays CPT is not a part of modern regional science and it has been criticized for being a static theory, not explaining how central places emerge and develop [? ].It has also been shown that while CPT is ideally suited to describe agricultural areas, it is less relevant to describe industrial areas which are highly diversified in nature.But despite its imperfections, CPT still remains one of the strongest economic theories for understanding the spatial organization of the society as the hierarchy of urban centers.It has been applied in regional and urban economies, describing the location of trade and service activity, and for describing consumer market-oriented manufacturing.

Models' variations
In parallel with the models presented in the main text, we also tested different variations of the gravity, radiation and hierarchy models to predict human interactivity.In the following section, we present the definition of these variations and standard benchmark measures quantifying how well these models fit the data.As for the models presented in the main text, we always use the total amount of communication w i originating from a location i as a proxy for its population.

Definitions
Radiation model versions The radiation model is a parameter-free model recently introduced in the context of migration patterns [? ].Given the distance d ij between two locations i and j, the model predicts that the flux of individual moves T ij between those two locations will depend on the population at the origin, the population at the destination and on the population s ij within the circle of radius d ij centered on the origin location i. Applied to our case, the radiation model is written as where s ij = k, 0<d ik <dij w k is the total amount of communication originating from locations at a distance shorter than d ij from location i and C i is a normalization factor ensuring that the predicted total activity of each node is the same as the actual one, i.e. j =i T Rad ij = j =i T ij .The model is otherwise parameter-free.We also applied the generalized version of the radiation model proposed in [? ], introducing a parameter λ which can be interpreted in our case as a fraction of individuals people will not consider as potential contacts because of lack of information about them of other reasons.This version of the radiation model can be written as where C i is again a normalization factor ensuring that j =i T GenRad ij = j =i T ij .Notice that the case λ = 0 corresponds to the original radiation model.

Gravity models
The gravity model is one of the oldest models describing human mobility and interaction, formulated in analogy to Newton's law of gravity.Here interaction strength or mobility fluxes between a source i and a destination j are originally proposed to be related to a function of the distance d ij between the two locations and the product of the (powers of) population at the source an at the destination, T Grav ij ∼ w α i w β j f (d ij ).We tested the two classical forms of the gravity model using a power or an exponential law as deterrence function: where in both cases C is a global normalization constant ensuring that i, j =i T Grav ij = i, j =i T ij and α, β, γ, d 0 are parameters to fit.For example, taking the logarithm on both sides of Eq. ( 4), we obtain log T GravP L ij = log C + α log w i + β log w j + γ log d ij .We then use, as is customary, log T GravP L ij to estimate the different parameters through a regression analysis [? ].We selected the set of parameters that minimized the deviance E (see definition in main text, Methods section).We also computed versions of these models where the population exponents are fixed, i.e. α = β = 1.
At this point, we observe an often overseen difference between the gravity model relying on a single global normalization factor C and the radiation model relying on local normalization factors {C i } -one for each location.One could argue that this difference gives an advantage to the radiation model.For the sake of comparison, we then also tested locally constrained gravity models that can be written as: where in each case C i is a normalization factor ensuring that the predicted total activity of each node is the same as the actual one.This type of constrained model was already presented in [? ] and more recently in [? ].
Hierarchical models The general idea behind the hierarchical model is to simply replace the actual distance used in the gravity models by a hierarchical distance.In its most generic definition, the hierarchical model predicts an interaction strength between location i and j to be written as T Hier ij = C i w α i w β j f (h ij ), where h ij is the hierarchical distance between these two locations based on the Combo partition or any other partition and f is a deterrence function.As it is done for the gravity model, one could a priori choose any deterrence function.We tested different simple forms of hierarchy models, using the hierarchical distances {h ij } provided either by the Combo partition or the administrative partition.Hierarchy models using a global normalization framework -as the usual gravity models do -read: where C is a global normalization constant, and α, β, γ and q are parameters which we fit by minimizing the deviance E, see below.We also computed versions of these models where the population exponents are fixed, i.e. α = β = 1.The locally constrained versions of models given in Eqs. ( 8) and ( 9) with fixed population exponent read: T HierEXP loc ij = C i w i w j q hij (10) It turns out that the model given by Eq (10) is naturally related to the notion of damping parameter.Indeed, assuming β = 1 in the generic functional form T Hier ij = C i w α i w β j f (h ij ) and taking into account that the normalization factor C i ensures w Hier an equation which immediately implies that choosing an exponential deterrence function f (h) = q h will ensure a constant damping parameter with respect to the locations and hierarchy levels q Hier(h) i = q for all i, h.
Hierarchical-Admin models The hierarchical models rely on the notion of hierarchical distances between locations, which depend on a given partition.For cases when the communication network or it's partitioning are not known, we can as well defined a model based on Administrative partitions of the countries.In the following, we refer to hierarchy models based on Administrative partition as 'hierarchy-admin' models.

Analysis
Benchmark measures (as defined in main text) of the different models along with their fitted parameters are reported in Tables VI and VII.We make a number of remarkable observations: • In every country and according to all benchmark measures, the generalized radiation model is significantly more appropriate than the original one to describe our communication networks (e.g. in UK, the cosine distance to the data goes from 0.344 in the original radiation model to 0.195 in the generalized radiation model, similarly the correlation to the data goes from 0.656 to 0.805).
• Locally constrained models always perform better compared to their 'global normalization' counterpart.
• When using an exponential deterrence function, the hierarchy models based on our Combo partition are always significantly better than the corresponding gravity models, with respect to the deviance or any other benchmark measure.E.g. in UK, the cosine distance goes from 0.595 in the 'gravity EXP' model to 0.278 in the 'Hierarchy EXP' model, similarly the correlation to the data goes from 0.402 to 0.72.The Hierarchy models using the administrative partition are also slightly better than the corresponding gravity models, except in Country X where we already observed some bias due to the administrative borders (see section II-D above).
• When using a power law deterrence function, the best model can vary with respect to the benchmark measure and the studied country.The 'Gravity PL loc' model is in general the best one.
• As one could expect, the models where the population exponent α and β are not constrained always have a lower deviance E than the corresponding models with α = β = 1 (but be aware that a correct comparison of the performance of the models based on the deviance should also take the number of parameter into account).
According to the other benchmark measures (not used for the fitting), the unconstrained models are also most often better fit than the constrained ones (e.g.their correlation to data is higher 12 times out of 18) and when it's not the case the difference between the fit measures is small.
• In general, the gravity PL model is better than the gravity EXP model.It is the best one regarding all benchmark measures in UK and Country X, and two out of five measures in Portugal.
• The hierarchy EXP models always fit better than the hierarchy PL models, all measures and countries considered.
The 20 tested models can be ranked according to the average of the benchmark measures taken over all countries, see Table VIII.While on average the best model is always the constrained 'Gravity PL loc' model, the 'hierarchy EXP loc' model presented in the main text is a close second (and first in some countries).7 of the top 10 models are versions of the hierarchy model.In particular, the 'hierarchy-admin EXP loc' (presented in main text), the 'hierarchy PL loc' and 'hierarchy PL' based on a power law deterrence function also outperform the state-of-the-art radiation and gravity PL models.
FIG.1.Partitioning of a country based on telephone call networks.Hierarchical distances between two locations are defined through three regional levels -either administrative ones or those found by applying iterative community detection on human interaction networks.Two distinct locations are at a hierarchical distance h = 1 if they are in the same L3 region, h = 2 if they are in different L3 regions but in the same L2 region, h = 3 if they are in different L2 regions but in the same L1 region and h = 4 if they are in different L1 regions.Note that a higher hierarchical distance does not necessarily correspond to higher geographical distance.

FIG. 2 .FIG. 3 .
FIG.2.Hierarchical properties of spatial organization from human interactions.a-e, Maps of L1 communities in telephone call networks detected from data and from various interaction models.Black lines correspond to the administrative partitioning of the 11 NUTS1 regions of UK, colored areas to regions detected by a community detection algorithm applied to (a) the data, and to the (b) gravity, (c) radiation, (d) hierarchy, and (e) administrative models.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.

FIG. 4 .
FIG. 4. Comparison of model predictions.a-d, Comparison of the actual total communication to the predicted communication for each pair of distinct locations, for the (a) gravity, (b) radiation, (c) hierarchy, and (d) administrative models.Gray markers are scatter plots for each pair of locations.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.

FIG. 1 .
FIG. 1. (FigureS1) Hierarchical properties of spatial organization from human interactions in Portugal.a-e, Maps of L1 communities in telephone call networks detected from data and from various interaction models.Black lines correspond to the administrative partitioning of the 5 NUTS1 regions of Portugal, colored areas to regions detected by a community detection algorithm applied to (a) the data, and to the (b) gravity, (c) radiation, (d) hierarchy, and (e) administrative models.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.

FIG. 3 .
FIG. 3. (FigureS3) Hierarchical properties of spatial organization from human interactions in Ivory Coast.a-e, Maps of L1 communities in telephone call networks detected from data and from various interaction models.Black lines correspond to the administrative partitioning of the 19 administrative regions of Ivory Coast, colored areas to regions detected by a community detection algorithm applied to (a) the data, and to the (b) gravity, (c) radiation, (d) hierarchy, and (e) administrative models.All detected regions are cohesive although some of the distinct colors used may appear similar.f-j, Probability distribution of number of subregions by region found in (f) the actual network and (g-j) in each model.The gravity model (g) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (h) strongly overestimates the number of L1 communities.The hierarchy model (i) correctly determines the distributions of sub-communities per community.k-o, Probability distributions of damping values q (h) .The hierarchy model (n) assumes a constant damping value for all levels.

FIG. 4 .FIG. 5 .
FIG. 4. (Figure S4) Hierarchical properties of spatial organization from human interactions in Country X. a-e, Probability distribution of number of subregions by region of Country X found in (a) the actual network and (b-e) in each model.The gravity model (b) underestimates the number of L1 communities but overestimates the numbers of subregions within regions.The radiation model (c) strongly overestimates the number of L1 communities.The hierarchy model (d)correctly determines the distributions of sub-communities per community.f-i, Probability distributions of damping values q (h) .The hierarchy model (h) assumes a constant damping value for all levels.Maps of of L1 communities are not shown as in other countries due to our non-disclosure agreement with the data providers from Country X.

FIG. 6 .FIG. 7 .
FIG. 6. (Figure S6) Comparison of model predictions in France.a-d, Comparison of the actual total communication to the predicted communication for each pair of distinct locations, for the (a) gravity, (b) radiation, (c) hierarchy, and (d)administrative models.Gray markers are scatter plots for each pair of locations.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.

FIG. 8 .
FIG. 8. (FigureS8) Comparison of model predictions in Country X. a-d, Comparison of the actual total communication to the predicted communication for each pair of distinct locations, for the (a) gravity, (b) radiation, (c) hierarchy, and (d) administrative models.Gray markers are scatter plots for each pair of locations.A box is colored green if the equality line y = x lies between the 9th and 91th percentiles in that bin and is red otherwise.Red boxes hence emphasize significant biases of the models.Black circles correspond to the average total communication of the pairs of locations in that bin.e-h, Goodness of prediction with respect to the hierarchical distance h, for the (e) gravity, (f) radiation, (g) hierarchy, and (h) administrative models.Gray markers are scatter plots for each pair of locations.Error bars show the corresponding 9th and 91th percentiles of total communication values.i-l, For each L3 community, comparison of the fractions of activity of model versus data between that L3 community and L3 communities at different hierarchical distances, for the (i) gravity, (j) radiation, (k) hierarchy and (l) administrative models.

FIG. 9 .
FIG. 9. (FigureS9) Comparing the model predictions.The proportion P dist (r) of communication occurring between two locations at a distance r (in km) from each other, is measured in the data and in the models.The radiation model is characterized by a lower than actual proportion of communication between distant (more than 100km) locations up to two orders of magnitude.It also presents a higher than actual proportion of communication between close (less than 10km with a peak between 0.5 and 5 km depending on the country) locations up to one order of magnitude.The gravity model presents in all countries a higher than actual proportion of communication between very close locations (100m-1km).In general, it also overestimates the low-range fluxes by 2 to 3 orders of magnitude and slightly underestimates the top-range fluxes.The hierarchy model fits almost perfectly at low distances (less than 10km).Depending on the country, it only deviates slightly from the data at top-range fluxes or estimates them properly.The fit of the hierarchy-admin model depends strongly on the country, but qualitatively comparable to the gravity and hierarchy models.

FIG. 10 .
FIG. 10. (FigureS10) Stability analysis: benchmark measures of the Hierarchy and Hierarchy-Admin models with varying parameter q.The dashed lines and circle markers show the 'optimal' values of q (reported in Supplementary Table2) minimizing the deviance E. In all countries, either for the hierarchy or hierarchy-admin model, these optimal values of q are also close to those minimizing the Dice (D), Sorensen (S) and Cosine (C) distances and maximizing the correlation (corr) between the data and the models.These optimal values are also stable, in the sense that close values of q (roughly between 0.1 and 0.3) still provide benchmark measures close to their optimum.The dotted lines in the case of the Hierarchy-Admin model indicate the damping value q = 0.2 matching robustly all countries.
P.H. acknowledges support of the German Academic Exchange Service (DAAD) via a postdoctoral fellowship.F.S. acknowledges support of the European FET-Open Project DATASIM (FP7-ICT-270833).S.G., M.S., S.S., and C.R. thank the National Science Foundation, the AT&T Foundation, the MIT SMART program, the Center for Complex Engineering Systems at KACST and MIT, Volkswagen Electronics Research Lab, BBVA, The Coca Cola Company, Ericsson, Expo 2015, Ferrovial, the Regional Municipality of Wood Buffalo and all the members of the MIT Senseable City Lab Consortium for supporting the research.The authors also thank Dashun Wang, Markus Schläpfer, Oleguer Sagarra and Santi Phithakkitnukoon for valuable feedbacks.• All authors designed the research, contributed to the model development and edited the manuscript; Z.S. provided the France dataset and A.-L.B. Country X dataset; C.R. proposed a general idea of the project; S.S. provided the concept of a model based on hierarchical distance and its initial validation; F.S. proposed the concept of damping parameter; S.G., P.H. and M.V. processed and cleaned data; S.G. and S.S. designed the algorithms; S.G., M.S. and M.V. implemented the algorithms; S.G., M.S., F.S. and P.H. analyzed the results; S.G. and M.S. were the lead writers of the manuscript. ACKNOWLEDGMENTS•

TABLE I .
Values of the damping value q for the actual and modeled networks in the UK.
TABLE II.Benchmark measures quantifying the goodness of fit in the UK.The Dice (D), Sorensen (S), Cosine (C) and deviance (E) are four different measures of the distance between the actual and modeled networks.The correlation corr measures a similarity between a model and the data.The parameters of the gravity and hierarchy models were chosen to minimize the value of E. .4560.543 α = 0.65, β = 0.65, γ = −1.46Radiation 1.622 0.624 0.632 0.344 0.656 Hierarchy 0.464 0.233 0.437 0.231 0.768 q = 0.139 Hierarchy-Admin 0.679 0.503 0.527 0.458 0.540 q = 0.2 (imposed)

TABLE III .
Over-/ under-estimation measures of link at specific hierarchical distance in the UK.

TABLE II .
(TableS2) Values of the damping parameter q for the actual and modeled networks in France, Portugal, Country X and Ivory Coast.

TABLE III .
(TableS3) Benchmark measures quantifying the goodness of fit in Portugal, France, Country X and Ivory Coast.The Dice (D), Sorensen (S), Cosine (C) and deviance (E) are four different measures of the distance between the actual and modeled networks.The correlation corr measures a similarity between a model and the data.The parameters of the gravity and hierarchy models were chosen to minimize the value of E.

TABLE IV .
(TableS4) Over-/ under-estimation measures of link at specific hierarchical distance in France, Portugal, Country X and Ivory Coast.

TABLE V .
Overlap between the administrative regions and the community found by the Combo algorithm.

TABLE VI .
Benchmark goodness fit measures The Dice (D), Sorensen (S), Cosine (C) and deviance (E) are four different fit values measuring a distance between the actual and modeled networks.The correlation corr measures a similarity between the model and the data.The different parameters of the gravity and hierarchy models were chosen to minimize the value of E. Stars (*) denote models where we imposed the population exponents α and β to be equal to 1.The rows corresponding to models presented in the main text are highlighted.− λ = 8.65 × 10 −10 UK / gravity EXP 941.1 0.777 0.579 0.595 0.402 α = 0.78, β = 0.78, d0 = 63.8kmUK / hierarchy EXP 547.7 0.278 0.447 0.278 0.720 α = 0.90, β = 0.90, q = 0.138 UK / hierarchy-admin EXP