Scale-free networks embedded in fractal space

# Scale-free networks embedded in fractal space

K. Yakubo Department of Applied Physics, Hokkaido University, Sapporo 060-8628, Japan    D. Korošak University of Maribor, Institute of Physiology, Faculty of Medicine, Maribor SI-2000, Slovenia University of Maribor, Faculty of Civil Engineering, Maribor SI-2000, Slovenia
July 9, 2019
###### Abstract

The impact of inhomogeneous arrangement of nodes in space on network organization cannot be neglected in most of real-world scale-free networks. Here, we propose a model for a geographical network with nodes embedded in a fractal space in which we can tune the network heterogeneity by varying the strength of the spatial embedding. When the nodes in such networks have power-law distributed intrinsic weights, the networks are scale-free with the degree distribution exponent decreasing with increasing fractal dimension if the spatial embedding is strong enough, while the weakly embedded networks are still scale-free but the degree exponent is equal to regardless of the fractal dimension. We show that this phenomenon is related to the transition from a non-compact to compact phase of the network and that this transition accompanies a drastic change of the network efficiency. We test our analytically derived predictions on the real-world example of networks describing the soil porous architecture.

###### pacs:
preprint: APS/123-QED

## I Introduction

Scale-free organization of networks barabasi2009 (); barabasi1999 (); albert2002 () seems to be the underlying principle common of many complex systems. In real-world networks, examples include the Internet, social networks or communication networks yook2002 (); eubank2004 (); liben-nowell2005 (); gastner2006 (); lambiotte2008 (); bianconi2009 (), the inhomogeneous arrangement of nodes in space strongly impacts the network organization and the linking rules must include the dependence on the distances between nodes xulvi-brunet2002 (); barthelemy2003 (); kaiser2004 (); xulvi-brunet2007 (). However, most of the network models studied so far considered either randomly distributed nodes in metric space dall2002 (); herrmann2003 (); masuda2005 (); morita2006 () or nodes placed on a lattice kleinberg2000 (); warren2002 (); ben-avraham2003 (); rozenfeld2002 (); kosmidis2008 (). The importance of inhomogeneous spatial positions of nodes was emphasized in yook2002 () where it was shown that the fractality of the space is one of the universal parameters that constrain the Internet models describing Internet’s large-scale topology and its observed scale-free character faloutsos1999 (); vazquez2002 (); siganos2003 () at the router and autonomous system level. While preferential attachment barabasi1999 (); albert2002 () seems to be the main underlying mechanism structuring the Internet, the original form of the preferential attachment barabasi1999 () should be altered yook2002 (); vazquez2002 (); shakkottai2010 () to account for the observed spatial and/or functional heterogeneity of the nodes.

Here, we study the structure of networks formed by geographical network model in which the nodes with power-law distributed intrinsic weights (i.e., fitnesses caldarelli2002 (); masuda2004 (); garlaschelli2007 ()) are embedded in a fractal space. We show analytically that the networks produced by such model are scale-free with the degree exponent influenced by the fractal dimension of the embedding space if the spatial embedding is strong enough. By explicitly deriving the degree and the edge-length distribution functions, we classify these networks into non-compact phase with infinite average degree and average edge length and compact phase with finite average degree and average edge length, separated by an intermediate phase characterized by finite average degree and infinite average edge length. It is also shown that the transition between these phases accompanies a drastic change of the network efficiency. Finally, we use our findings in the analysis of the networks describing the soil porous architecture as an example.

## Ii Model

In order to describe an inhomogeneous distribution of nodes in a space of linear dimension , the nodes are isotropically distributed in a fractal manner with the fractal dimension . In the theoretical treatment we employ an artificial boundary conditions for which every node position is statistically equivalent. This is slightly different from usual periodic boundary conditions leading to an anisotropy. The number of nodes is thus given by , where is the -dimensional solid angle and is the density of nodes. Each node has a real and continuous fitness randomly assigned according to the distribution function . If a node pair - satisfies

 F(xi,xj)lmij>Θ , (1)

then these two nodes are connected by an edge, where is the Euclidean distance between the nodes and , is a real parameter quantifying the strength of the spatial embedding (namely, the strength of the geographical effect), is a threshold value, and is a function relating the two fitnesses to the connectivity. If , we have a conventional fitness model which provides a scale-free network for a variety of combinations of forms of and caldarelli2002 (); masuda2004 (). Here, we concentrate on the case of

 F(x,y)=xy , (2)

and

 s(x)=s0x−α , (3)

where (and ) are in the range of and . From the normalization condition of , we have

 s0=(α−1)xα−1min . (4)

Similar geographical network models have been studied under the assumption that nodes are homogeneously distributed in a Euclidean space masuda2005 (); morita2006 (), while here we investigate how the inhomogeneity (fractality) of spatial distribution of nodes affects properties of the network.

## Iii Degree distribution function

First, we calculate the degree distribution function of the network formed by the above geographical algorithm. Let be the number of nodes connected to the node and included in a thin spherical shell of the radius (width ) centered at the position of the node . Since the distance between the node and a node in the shell is , the connectivity condition Eq. (1) tells us that nodes with can connect to the node . Thus, the number of connected nodes is

 ki(l)dl=n(l)dl⋅∫∞Θlm/xis(x)dx , (5)

where is the number of nodes in this spherical shell. In this expression, we assume that , which is equivalent to where

 lmin(xi)=(xminxiΘ)1/m . (6)

Nodes within the distance from the node can connect to the node . Thus, using Eq. (4), we have

 ki(l)= ρΩ(xminΘ)α−1xα−1ilD−1−m(α−1) , l>lmin(xi) , (7a) ki(l)= ρΩlD−1 , l≤lmin(xi) . (7b)

Integrating with respect to over we obtain the total number of nodes connected to the node , or the degree of the node given by

 ki = ρΩxα−1minLD−m(α−1)[D−m(α−1)]Θα−1xα−1i (8) +ρΩxD/mminΘD/m[1D−1D−m(α−1)]xD/mi.

It should be noted that geometrical coefficients in two terms of Eq. (8) coming from the volume integration are not very meaningful for realistic systems because of our artificial boundary conditions. In the case of , the first term of Eq. (8) dominates the second term when is sufficiently large and we have . The degree distribution function calculated from the relation is then proportional to independently of , , and . If however, , the second term of Eq. (8) becomes much larger than the first term. The relation in this case leads to . Therefore, the degree distribution function of the present geographical network model is given by

 P(k)∝ k−2 , m≤mc0 , (9a) P(k)∝ k−m(α−1)/D+1 , m>mc0  , (9b)

where

 mc0=Dα−1 . (10)

The distribution function obeys power-law forms in both cases of and , and the degree exponent is

 γ=⎧⎨⎩2, m≤mc0 ,mD(α−1)+1, m>mc0 . (11)

The geographical inhomogeneity of the node distribution does not affect the degree distribution of the network for a weak geographical effect (small ), while depends on for a strong effect (large ).

In the above argument, we assumed that the system size is always larger than for any , because we are interested in the thermodynamic case (). This is, however, not obvious, because [and then ] can also diverge in the thermodynamic limit under a constant density . Let us consider carefully this condition . In a finite system with nodes, the fitness is truncated at a finite value. The maximum fitness is given by . Using the fitness distribution Eq. (3) with Eq. (4), the quantity is given by

 xmax=N1/(α−1)xmin . (12)

Thus, the length can be as large as , where

 lmin(xmax)=⎡⎣(ρΩ)1/(α−1)LD/(α−1)x2minΘD1/(α−1)⎤⎦1/m. (13)

From the above expression, the condition is equivalent to , where

 Θ0=(ρΩD)1α−1x2minLDα−1−m . (14)

If , the quantity goes to zero in the thermodynamic limit and any finite satisfies , namely, . Thus, the degree distribution function is given by Eq. (9b) in the thermodynamic limit with . On the contrary, if , is always less than because diverges as . In this case, there must be nodes satisfying for which . Since the condition implies that any node in the whole system can connect to the node , the degree of such a node is independent of . Thus, these nodes give an additional -functional contribution [] to the degree distribution given by Eq. (9a) for nodes with . Let us estimate the magnitude of this -functional contribution. It is proportional to the number of nodes having . The quantity is given by and can be written as

 n0=N(Lξ)−m(α−1) , (15)

where is the node-pair distance defined by

 ξ=lmin(xmin)=⎛⎝x2minΘ⎞⎠1/m , (16)

below which the two nodes are connected independently of the fitness. The properly normalized -functional part of is then presented by . Since is always positive, the -functional contribution vanishes in the thermodynamic limit. Therefore, Eq. (9b) is valid both for the cases of and in an infinite system. For a finite , however, can be chosen to be less than independently of and the -functional contribution remains finite. Thus, the degree distribution function of a finite system with must have the -functional correction term, that is,

 P(k)=p0k−γ+(Lξ)−m(α−1)δ(k−N+1) , (17)

for both and . Here, is a normalization constant and the exponent is given by Eq. (11). It should be noted that Eq. (17) holds for a finite but large because the first term of Eq. (17) is valid in the large limit.

The quantity is a characteristic value of peculiar to a finite system with a fixed density . There are two other characteristic values of for a finite . One is below which every node can connect to all other nodes. Network constructed under becomes the complete graph. Obviously, is given by

 Θmin=x2minLm . (18)

Another characteristic is above which no node can connect to any other nodes. Network with is a set of isolated nodes. We have

 Θmax=x2maxΔlm , (19)

where is the minimum edge length. Using Eq. (12), is written as

 Θmax=(ρΩD)2α−1+mDx2minL2Dα−1 . (20)

In a finite system we always assume that satisfies the condition leading to non-trivial networks. It is not necessary to consider this condition in an infinite system because and vanishes and diverges, respectively, in the thermodynamic limit.

## Iv Relation between ⟨k⟩ and Θ

From Eq. (9b), it is clear that the average degree diverges for because of and it remains finite for in the thermodynamic limit. The average degree in a finite system is, however, always finite and depends on . A large restricts a connection of a node pair and leads a small . It is important to know the relation between and for a finite but large system. We calculate the average degree by

 ⟨k⟩=∫xmaxxminkis(x)dx , (21)

instead of , using the derived expression (not asymptotic form) of for a finite system.

Let us consider Eq. (21) separately for and for . In the case of , the length can be larger than , which implies . Thus, the integral of Eq. (21) is separated into two regions

 ⟨k⟩=∫ΘLm/xminxmink(x)s(x)dx+(N−1)∫xmaxΘLm/xmins(x)dx , (22)

where is given by Eq. (8) regarding as a continuous variable . The coefficient in the second term comes from the fact that nodes satisfying connect to all nodes. Using Eq. (8) and approximating by , we have

 ⟨k⟩ = Xρ⎛⎝x2minΘ⎞⎠α−1LD−m(α−1)log(YΘLmxmin2) (23) +Zρ⎛⎝x2minΘ⎞⎠D/m ,

where

 X = Ω(α−1)D−m(α−1) , (24) Y = (25) Z = Ωm2(α−1)2D[D−m(α−1)]2 . (26)

It should be again emphasized that these geometrical quantities , , and resulting from the volume integration depend strongly on the boundary conditions and are not very meaningful. We can evaluate the asymptotic behavior of for large by using Eq. (23). In the case of , the first term of Eq. (23) obviously dominates the second term. Then, ignoring unimportant geographical coefficient, behaves asymptotically as

 ⟨k⟩∼Θ1−α(logΘ+c) . (27)

where is a constant depending on the boundary condition. On the other hand, a careful treatment is required for . It seems that the second term of Eq. (23) dominates the first term for in the thermodynamic limit. However, we should note that must be infinitesimal to satisfy the condition in this calculation because for goes to zero for so it is not obvious which term is dominating in Eq. (23). In order to find the dominant term, we evaluate the lower bounds of these terms by replacing by and using (14). Then, the lower bounds of the first and the second terms are proportional to and respectively. This suggests that the second term dominates the first term for large because the exponent is positive for . We have then

 ⟨k⟩∼Θ−D/m , (28)

for and . At , both terms in Eq. (23) should be considered.

Next, we treat the case of . Here is always less than , and the degree is given by Eq. (8) for any in the range of . The average degree is then simply presented by

 ⟨k⟩ = ∫xmaxxmink(x)s(x)dx (29) = ρXα−1⎛⎝x2minΘ⎞⎠α−1LD−m(α−1)log(WLD) −ρZ[D−m(α−1)]⎛⎝x2minΘ⎞⎠D/m(WLD)Dm(α−1)−1 +ρZ[D−m(α−1)]⎛⎝x2minΘ⎞⎠D/m ,

where and we used Eq. (8) for and Eq. (12). For , it is easy to understand that the third term dominates other two terms for large . So we have

 ⟨k⟩∼Θ−D/m . (30)

In the case of , the infinitely large must be considered when we find the dominant term of Eq. (29), because is larger than and diverges for in the thermodynamic limit. As in the case of and , replacing in Eq. (29) by , the dependence of the upper bounds of these terms points to the dominant term. Since the upper bounds of the first, second, and third terms are proportional to , , and , respectively, the first term dominates the second and third terms because of . Therefore, the average degree is asymptotically given by

 ⟨k⟩∼Θ1−α . (31)

This relation differs from Eq. (27) by a logarithmic correction.

In summary, the relation between and is given by

 ⟨k⟩∼ Θ1−α(logΘ+c) , mmc0 (32b)

for and

 ⟨k⟩∼ Θ1−α , mmc0 . (33b)

for . At , is related to through Eq. (23) for and through Eq. (29) for , because every term contributes equally to even in the thermodynamic limit. A similar result to Eq. (33b) has been obtained by morita2006 () where the nodes were uniformly distributed in a -dimensional space with the -max norm and is fixed at . The asymptotic dependence of for a fixed can be also evaluated from Eqs. (23) and (29). In the case of and enough large (then, ), the dominant first term of Eq. (23) gives . For (and then ), we have . These dependences are consistent with those calculated by by using Eqs. (17) and (9b) for and , respectively, and taking into account the dependence of in Eq. (17).

## V Edge-length distribution function

In this section, we derive the edge-length distribution function of our geographical networks embedded in a fractal space. To this end, we regard given by Eq. (7b) as a continuous function of the fitness and the edge length . The average number of edges, , of length from a given node is obtained by averaging over the fitness , i.e.,

 k(l)=∫xmaxxmins(x)k(x,l)dx. (34)

Equation (7b) expresses the forms of by separating two cases and for a fixed . Corresponding to this classification, in Eq. (34) for a fixed has different forms for and , respectively. Thus, the integral of Eq. (34) is calculated as

 k(l) = ρΩ(xminΘ)α−1lD−1−m(α−1)∫xlxmins(x)xα−1dx (35) +ρΩlD−1∫∞xls(x)dx ,

where . Here, is assumed to be larger than (namely, ) and the upper cut-off of the integral is, as an approximation, extended to infinity. Using Eqs. (3), (4), (7b), and (16), is expressed by

 k(l)=ρΩlD−1(lξ)−m(α−1)[1+m(α−1)log(lξ)] . (36)

The probability distribution function is given by

 R(l)=k(l)∫Lk(l′)dl′ , (37)

where the integration in the denominator is done over the whole range of . Neglecting the normalization constant, the edge-length distribution is

 R(l)∝lD−1(lξ)−m(α−1)[1+m(α−1)log(lξ)] . (38)

We should remark that for goes to infinity as . Since the length does not exceed the size in a finite system, the distribution is actually truncated at . In the thermodynamic limit, however, we must consider the infinitesimal normalization constant coming from the denominator of Eq. (37).

In the case of (namely, ), in Eq. (34) is given by Eq. (7b) for any in the integration range . Thus, is given by , namely,

 k(l)=ρΩlD−1 , (39)

and then,

 R(l)∝lD−1 . (40)

It should be noted that the distribution depends on the fractal dimension independently of while the degree distribution does not depend on for (weak geographical effect region).

Two expressions Eqs. (38) and (40) of for and must coincide at . This condition concludes that the proportionality coefficients for Eqs. (38) and (40) are identical. We can derive the common coefficient for an infinite system from the normalization condition of given by

 1 = C∫ξ0lD−1dl (41) +C∫∞ξlD−1(lξ)−β[1+βlog(lξ)]dl ,

where . When , this equation provides a finite coefficient expressed by

 C=DξD(1−mc0m)2 , (42)

while is infinitesimal for as mentioned above. It should be noted that the coefficient does not depend on the boundary condition because the boundary-condition dependent factors in the numerator and the denominator in Eq. (37) are canceled out.

From the above argument, we can immediately derive the probability of two nodes with distance to be connected. This probability is given by the ratio of the number of connected nodes to nodes located at distances from a given node, namely, . Since as mentioned below Eq. (5) and is given by Eq. (36) or (39) for or respectively, we have

 g(l)=⎧⎪⎨⎪⎩(lξ)−m(α−1)[1+m(α−1)log(lξ)], l>ξ ,1, l≤ξ . (43)

This expression can be alternatively derived directly from the meaning of ,

 g(l)=∫∞xmins(x)dx∫xy/lm>Θs(y)dy . (44)

Considering that two nodes are always connected if the fitness of one node exceeds , we can separate the above integration into two parts as

 g(l)=∫∞xls(x)dx+∫xlxmins(x)dx∫∞Θlm/xs(y)dy , (45)

for (). For , the integral range of the second integral in Eq. (44) is spread over the whole region of , then . These equations again lead Eq. (43) if we use Eq. (3). We should note that the probability given by Eq. (43) does not depends on the fractal dimension for any value of .

## Vi Numerical confirmations

The above analytical results are numerically verified in this section. At first, we confirm the behavior of the degree distribution function for geographical networks on the Sierpinski node set (Sierpinski geographical networks). Nodes are located on vertices of the Sierpinski gasket with the fractal dimension and connected by edges according to the condition Eq. (1) with Eqs. (2) and (3). In our numerical calculations in this work, the lower cut-off is set to be unity. A typical network on the 6th generation Sierpinski gasket () is depicted in Fig. 1. There exist nodes (hubs) possessing a large number of edges. Numerically calculated degree distribution functions for two networks with different and are presented in Fig. 2. The scale-free property of networks formed by our algorithm is clearly shown in this figure. The scale-free exponents calculated numerically for many Sierpinski geographical networks with different combinations of and are plotted in Fig. 3 as a function of . In this figure, values of for networks with nodes distributed homogeneously in two-dimensional Euclidean space (2d-geographical networks) are also plotted. We should remark that our theoretical arguments are valid for a geographical network with homogeneously distributed nodes which is a special case with the Euclidean dimension instead of the fractal dimension . In fact, Eq. (11) with reproduces the results by morita2006 (); masuda2005 (). We see that all numerical results collapse onto the theoretical line given by Eq. (11). These results strongly support the theoretical predictions on presented in Sec. III.

Next, we confirm numerically the relation between and for 2d-geographical networks. We have to evaluate defined by Eq. (14) to distinguish four regions with respect to and corresponding to Eqs. (32b) and (33b). In addition, must satisfy the condition , where and are given by Eqs. (18) and (20), respectively. Rewriting Eqs. (14) and (20) by and , we can estimate and without treating the boundary-condition dependent geometrical factor . Figures 4(a) and 4(b) show numerically calculated as a function of for . In these calculations, the exponent characterizing the strength of the geographical effect are chosen to be , and nodes are distributed in a square of size . The fitness is allocated to each node according to the distribution function given by Eq. (3) with . The square system has the periodic boundary conditions in the and directions. Since , , , and from these parameters, the conditions and are satisfied for Figs. 4(a) and 4(b), respectively, and for both figures. Numerically calculated s are well described by solid lines representing Eqs. (32a) and (33a), where the constant in Eq. (32a) and prefactors are suitably chosen. Results for are presented in Figs. 4(c) and 4(d) which show the - relation for and , respectively. As in the cases of Figs. 4(a) and 4(b), nodes are distributed in a square of size . Parameters and are chosen as and so that the condition is satisfied. In order to realize the condition in Fig. 4(c), we treated a large network with , which has , , and . On the contrary, the number of nodes for Fig. 4(d) is , in which is satisfied with , , and . Similarly to Figs. 4(a) and 4(b), numerical results agree well with the theoretical results shown by solid lines.

Finally, the behavior of the edge-length distribution function is examined also for 2d-geographical networks. Results (dots) shown in Fig. 5 are calculated in the condition of , for which can be properly normalized even in the thermodynamic limit. Solid and dashed curves represent the theoretical prediction given by Eq. (38) and (40) with the common prefactor presented by Eq. (42). It should be emphasized that there is no fitting parameter to obtain the theoretical curve. Numerical results agree quite well with the theoretical prediction. The peak structure and the tail profile slightly deviating from a power law (see the inset of Fig. 5) result from the logarithmic term of Eq. (38). The reason of the slight deviation between numerical results and the theoretical line for (see the inset) is due to a finite defined below Eq. (19).

## Vii Compactness and efficiency

From the functional forms of and , we can immediately find the convergence property of the average degree and the average edge length in the thermodynamic limit. As argued at the end of Sec. IV, the quantity for diverges as and it converges for . On the other hand, the convergence of is governed by Eq. (38) describing in the asymptotic regime. From Eq. (38), the average diverges for and remains finite for , where

 mc1=D+1α−1 . (46)

The behavior of and suggests that the impact of the geographical effect can be classified into three regions. For , both quantities and diverge in the thermodynamic limit. Since a node connects with a huge number of nodes far away we call this region the non-compact phase. In the same sense, the region of is termed the compact phase where both quantities and remain finite and a node connects with only a small number of nodes in its vicinity. In the intermediate phase, i.e., , the average degree converges but diverges. These regions are summarize in Table 1.

Let us consider the relation between the compactness of a network and the geographical efficiency defined by

 e=2N(N−1)∑i>j1rij , (47)

where is the shortest Euclidean distance between two nodes and along a network path. This quantity is a natural extension of the global efficiency defined by

 E=2N(N−1)∑i>j1dij , (48)

where is the shortest network distance between two nodes and introduced in latora2001 (). Both quantities and characterize the efficiency of the information exchange or the flow in the network. If the efficiency is governed mainly by the Euclidean distance along the path rather than the number of steps, the geographical efficiency is more meaningful than , and vice versa.

We calculated numerically these quantities and examined how the compactness is related to the efficiencies and . Solid lines in Fig. 6 represent the efficiencies for Sierpinski geographical networks in three different generations (having the same linear dimension ), while dashed lines show them for 2d-geographical networks with different numbers of nodes in squares of size . In order to eliminate the and dependence of the maximum and minimum values of the efficiencies, we plot the rescaled quantities and defined by and , where () and () are at and , respectively. The rescaled geographical (global) efficiency () rapidly increases (decreases) near as increasing . From the fact that enlarging the system the slope near becomes steeper and the point at the intersection of and approaches , the efficiencies and in the thermodynamic limit are supposed to show step-like forms at . This behavior can be interpreted as follows. In the non-compact phase (), existing short-cut edge (in the topological sense) connects a starting node and a target node by a small number of edges with going back-and-forth around the target node. This back-and-forth motion, however, requires an extra Euclidean distance and leads the low geographical efficiency. On the other hand, in the compact phase, the network structure resembles the structure of a regular lattice (or regular Sierpinski gasket). Although we need lots of edges to connect two distant nodes, the total Euclidean length along the path can be minimized as the geodetic distance between two nodes. Thus, the geographical efficiency becomes large in this region. The above consideration supports that the global efficiency measuring the number of edges to connect nodes is large for and small for . It should be noted that the abrupt change in or occurs at but not at though the transition point is related to the edge length. Since and behave oppositely, it is crucial to clarify which efficiency is more relevant to a given problem by considering how strongly the cost of the flow is influenced by the Euclidean distance. It is also interesting that both efficiencies and are relatively high in the network at at which the competition between order and disorder in the geometrical sense is balanced.

## Viii Example

We will demonstrate the above described approach using soil-pore networks as an example of a real-world spatially embedded network. Two approaches to build complex network models of soil-pore organization have recently been developed mooney2009 (); santiago2008 (); cardenas2010 (). Here, we will concentrate on the networks presented in mooney2009 () formed by geographical threshold algorithm as described in Sec. II. We consider a set of pores representing the nodes of the network. The nodes of the network are the centers of the pores and the links between nodes are drawn according to Eq. (1) with Eqs. (2) and (3), where in this case the continuous fitness variable is the size of the pore. Pore sizes and their relative positions of actual soil specimens are obtained from the image analysis of 2d soil X-CT scans mooney2009 ().

An example of the soil-pore network overlain on the 2d soil porous structure image is presented in Fig. 7 (from mooney2009 ()). The pore-size distribution of the soil sample shown in Fig. 7 is analyzed to be of the form with (Fig. 8). The fractal dimension of the soil-pore structure is also determined using box-counting method as shown in Fig. 9. Our theory predicts that the soil-pore network has a power-law degree distribution function and the degree exponent depends on above . Figure 10 shows the degree distributions for three networks based on soil image data constructed with different values of the parameter (, , and ). In this figure, the power-law fits (dashed lines) and calculated exponents of the degree distributions are shown together with theoretically predicted scale-free exponents from Eq. (11). The results clearly demonstrate the agreement of the empirically obtained scale-free exponents and the theoretically predicted ones. This also shows how the fractality of pore spatial distribution is reflected in the network organization. By measuring the scaling exponent of the degree distribution for any and the exponent of the pore-size distribution, the fractal dimension can be determined.

## Ix Conclusions

We have proposed a geographical scale-free network model with the nodes embedded in a fractal space and analytically and numerically studied several network properties. The fractal dimension of the embedding space was found to influence the scale-free exponent as only if the spatial embedding is strong enough (i.e., when ) otherwise . The analyses of the average degree and average edge length revealed that this type of network can exist either in the non-compact, compact, or intermediate phase depending on the importance of the spatial arrangement of nodes. We derived the edge-length distribution functions for our network model and showed that it has a peak-like structure similar to the profile of the shortest-path-distance distribution observed in a large-scale structure of the Internet vazquez2002 (); mahadevan2006 (). It is interesting to apply our approach to modeling the Internet at the autonomous systems level considering the observed long-tailed distribution of autonomous systems sizes lakhina2003 (). The measured degree distribution exponent of the Internet is slightly larger that () vazquez2002 (); siganos2003 () and seems to be decreasing with time siganos2003 (). In our network model this would be an indication of the evolution of the Internet towards the non-compact phase ().

We hope that our work will help in advancing the understanding of the complex systems in which the heterogeneity of intrinsic properties and the spatial arrangement of the elements play an important role.

###### Acknowledgements.
This work was supported in part by a Grant-in-Aid for Scientific Research (No. 22560058) from Japan Society for the Promotion of Science and by grant No. J3-2290 from Slovenian Research Agency. Numerical calculations in this work were performed on the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

## References

• (1) A.-L. Barabasi, Science 325, 412 (2009).
• (2) A.-L. Barabasi and R. Albert, Science 286, 509 (1999).
• (3) R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, 47 (2002).
• (4) S.-H. Yook, H. Jeong, and A.-L. Barabàsi, Proc. Natl. Acad. Sci. U.S.A. 99, 13382 (2002).
• (5) S. Eubank, H. Guclu, V. S. A. Kumar, M. V. Marathe, A. Srinivasan, Z. Toroczkai, and N. Wang, Nature 429, 180 (2004).
• (6) D. Liben-Nowell, J. Novak, R. Kumar, P. Raghavan, and A. Tomkins, Proc. Natl. Acad. Sci. U.S.A. 102, 11623 (2005).
• (7) M. T. Gastner and M. E. J. Newman, Eur. Phys. J. B 49, 247 (2006).
• (8) R. Lambiotte, V. D. Blondel, C. de Kerchove, E. Huens, C. Prieur, Z. Smoreda, and P. Van Doorena, Physica A 387, 5317 (2008).
• (9) G. Bianconi, P. Pin, and M. Marsili, Proc. Natl. Acad. Sci. U.S.A. 106, 11433 (2009).
• (10) R. Xulvi-Brunet and I. M. Sokolov, Phys. Rev. E 66, 026118 (2002).
• (11) M. Barthelemy, Europhys. Lett. 63, 915 (2003).
• (12) M. Kaiser and C. C. Hilgetag, Phys. Rev. E 69, 036103 (2004).
• (13) R. Xulvi-Brunet and I. M. Sokolov, Phys. Rev. E 75, 046117 (2007).
• (14) J. Dall and M. Christensen, Phys. Rev. E 66, 016121 (2002).
• (15) C. Herrmann, M. Barthelemy, and P. Provero, Phys. Rev. E 68, 26128 (2003).
• (16) N. Masuda, H. Miwa, and N. Konno, Phys. Rev. E 71, 036108 (2005).
• (17) S. Morita, Phys. Rev. E 73, 035104 (2006).
• (18) J. M. Kleinberg, Nature 406, 845 (2000).
• (19) C. P. Warren, L. M. Sander, and I. M. Sokolov, Phys. Rev. E 66, 056105 (2002).
• (20) A. F. Rozenfeld, R. Cohen, D. ben-Avraham, and S. Havlin, Phys. Rev. Lett. 89, 218701 (2002).
• (21) D. ben-Avraham, A. F. Rozenfeld, R. Cohen, and S. Havlin, Physica A 330, 107 (2003).
• (22) K. Kosmidis, S. Havlin, and A. Bunde, Europhys. Lett. 82, 48005 (2008).
• (23) M. Faloutsos, P. Faloutsos, and C. Faloutsos, Comput. Commun. Rev. 29, 251 (1999).
• (24) A. Vazquez, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. E 65, 066130 (2002).
• (25) G. Siganos, M. Faloutsos, P. Faloutsos, and C. Faloutsos, IEEE-ACM Trans. Netw. 11, 514 (2003).
• (26) S. Shakkottai, M. Fomenkov, R. Koga, D. Krioukov, and K. C. Claffy, Eur. Phys. J. B 74,271 (2010).
• (27) G. Caldarelli, A. Capocci, P. De Los Rios, and M. A. Munoz, Phys. Rev. Lett. 89, 258702 (2002).
• (28) N. Masuda, H. Miwa, and N. Konno, Phys. Rev. E 70, 36124 (2004).
• (29) D. Garlaschelli, A. Capocci, and G. Caldarelli, Nature Physics 3, 813 (2007).
• (30) V. Latora and M. Marchiori, Phys. Rev. Lett. 87, 198701 (2001).
• (31) S. J. Mooney and D. Korošak, Soil Sci. Soc. Am. J. 73, 1094 (2009).
• (32) A. Santiago, J. P. Cardenas, J. C. Losada, R. M. Benito, A. M. Tarquis, and F. Borondo, Nonlin. Processes Geophys. 15, 893 (2008).
• (33) J. P. Cardenas, A. Santiago, A. M. Tarquis, J. C. Losada, F. Borondo, and R. M. Benito, Geoderma, doi:10.1016/j.geoderma.2010.04.024 (in press).
• (34) P. Mahadevan, D. Krioukov, M. Fomenkov, B. Huffaker, X. Dimitropoulos, K. Claffy, and A. Vahdat, Comput. Commun. Rev. 36, 17 (2006).
• (35) A. Lakhina, J. W. Byers, M. Crovella, and I. Matta, IEEE J. Sel. Areas Commun. 21, 934 (2003).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters