True scale-free networks hidden by finite size effects

TL;DR

This study investigates whether many real-world networks are truly scale-free or if finite-size effects obscure their scale invariance, finding that many networks do exhibit scale-free properties once finite-size effects are accounted for.

Contribution

The paper demonstrates that many real networks are inherently scale-free, with finite-size effects often hiding this property, challenging previous assumptions based solely on degree distribution analysis.

Findings

Many networks follow finite size scaling without self-tuning.

Biological, technological, and informational networks often exhibit true scale invariance.

Infrastructure and social networks show deviations from scale-free behavior.

Abstract

We analyze about two hundred naturally occurring networks with distinct dynamical origins to formally test whether the commonly assumed hypothesis of an underlying scale-free structure is generally viable. This has recently been questioned on the basis of statistical testing of the validity of power law distributions of network degrees by contrasting real data. Specifically, we analyze by finite-size scaling analysis the datasets of real networks to check whether purported departures from the power law behavior are due to the finiteness of the sample size. In this case, power laws would be recovered in the case of progressively larger cutoffs induced by the size of the sample. We find that a large number of the networks studied follow a finite size scaling hypothesis without any self-tuning. This is the case of biological protein interaction networks, technological computer and hyperlink networks, and informational networks in general. Marked deviations appear in other cases, especially infrastructure and transportation but also social networks. We conclude that underlying scale invariance properties of many naturally occurring networks are extant features often clouded by finite-size effects due to the nature of the sample data.

Tables

Table 1. Table 1: Classification of empirical networks (split into categories). For each category we report the total number of networks and the percentage of SSF, WSF and NSF instances. For detailed results on each network analyzed, see the Supplementary Dataset Table.

	TOTAL	Affiliation	Annotation	Authorship	Biological	Citation	Computer	Hyperlink	Infrastructure	Social	Text
number	185	8	38	15	30	5	13	14	12	39	11
SSF	27%	63%	21%	27%	40%	40%	39%	22%	0%	13%	55%
WSF	23%	12%	24%	20%	30%	0%	38%	21%	17%	18%	27%
NSF	50%	25%	55%	53%	30%	60%	23%	57%	83%	69%	18%

Equations

$$P(k,N)=k^{-\gamma }f(k{N}^d)$$

$$P(k,N)k^{\gamma }=f(k{N}^d).$$

$$P(k,E)k^{\gamma }=f_E(k{E}^{d_E}),$$

$$⟨k^i⟩={\int }_{k_{min}}^{\infty }dk{k}^{i-1}{k}^{-\gamma }f(k{N}^d)\propto N^{-d(i-\gamma )}$$

$$\left\langle k^{i}\right\rangle\Big{/}\left\langle k^{i-1}\right\rangle\propto N^{-d},$$

$$⟨k^i⟩={\int }_{k_{min}}^{\infty }dk{k}^{i-1}{k}^{-\gamma }{f}_E(k{E}^{d_E})\propto E^{-d_E(i-\gamma )}$$

$$\left\langle k^{i}\right\rangle\Big{/}\left\langle k^{i-1}\right\rangle\propto E^{-d_{E}}.$$

$$\{\begin{array}{l}{x}_{nj}=k_j{n}^d\\ y_{nj}=P(k_j,n)k_j^{\gamma }\end{array}$$

$$S=\frac{1}{3|M|}\sum _{(n,j)\in M}\frac{(y_{nj}-Y_{nj}{)}^2}{d{y}_{nj}^2+d{Y}_{nj}^2},$$

$$Y_{nj}=\frac{W_{xx}{W}_y-W_x{W}_{xy}}{\eta }+x_{nj}\frac{W{W}_{xy}-W_x{W}_y}{\eta }$$

$$d{Y}_{nj}^2=\frac{1}{\eta }(W_{xx}-2x_{nj}{W}_x+x_{nj}^{2}W)$$

Full text

network form and function, degree distribution, power laws, finite size scaling, statistical physics

Networks play a vital role in the development of predictive models of physical, biological, and social collective phenomena [1, 2, 3]. A quite remarkable feature of many real networks is that they are believed to be approximately scale-free: the fraction of nodes with $k$ incident links (the degree) follows a power law $p(k)\propto k^{-\lambda}$ for sufficiently large value of $k$ [4, 5]. The value of the exponent $\lambda$ as well as deviations from power law scaling provides invaluable information on the mechanisms underlying the formation of the network such as small degree saturation, variations in the local fitness to compete for links, and high degree cut-offs owing to the finite size of the network. Indeed real networks are not infinitely large and the largest degree of any network cannot be larger than the number of nodes. Finite size scaling [6, 7, 8, 9, 10, 11, 12], firstly developed in the field of critical phenomena and renormalization group, is a useful tool for analyzing deviations from pure power law behavior as due to finite size effects. Here we show that despite the essential differences between networks and critical phenomena, finite size scaling provides a powerful framework for analyzing the scale-free nature of empirical networks.

The search of ubiquitous emergent properties occurring in several different systems and transcending the specific system details is a recurrent theme in statistical physics and complexity science [13]. Indeed the presence and the type of such “universal” law gives insights on the driving processes or on the characteristic properties of the observed system. Notably, complex systems have the propensity to display “power law” like relationship in many diverse observables (such as event sizes and centrality distribution, to name a few). In particular the power law shape of the degree distribution, which is the hallmark of scale-free networks, leads to important emergent attributes such as self-similarity in the network topology, robustness to random failures and fragility to targeted attacks. Notably scale invariance extends far beyond the degree distribution, affecting many other quantities as weighted degree, betweenness [14] and degree-degree distance [15].

In the last decade the existence of such power laws in complex networks (but also in other areas [16], e.g., law in language [17]) has been questioned [18]. A reason of the shift in such conclusion is in the availability of larger (and new) datasets, and especially in improved statistical methods. Recently, Broido and Clauset[19] fitted a power law model to the degree distribution of a variety of empirical networks and suggested that scale-free networks are rare. Voitalovet al.[20] rebutted that scale-free networks are not as rare if deviations from pure power law behavior are permitted in the small degree regime. The different conclusions may depend on very fine but critical assumptions at the basis of the statistical test for the power law hypothesis. Moreover, a crucial point that is typically ignored but represents the condition for the proper use of maximum likelihood methods is the independence of the empirical observations [21]. In this work we tackle the problem of detecting power laws in networks from a different perspective, based on the the machinery of finite size scaling.

Statistical physics of critical phenomena teaches us that a system at criticality exhibits power law singularities of physical quantities such as, for example, the compressibility, the specific heat, the density difference between the liquid and vapor, as well as the latent heat. Water at its critical point exhibits fluctuations at all scales between the molecular length scale and the size of the container, which could be macroscopically large. Moreover, one finds thoroughly mixed droplets of water and bubbles of gas. Indeed, any large part of the system looks like the whole – the system is self-similar. The length scale of these droplets and bubbles extends from the molecular scale up to the correlation length, which is a measure of the size of the largest droplet or bubble. The divergence of the correlation length in the vicinity of a phase transition at the thermodynamic limit thus suggests that properties near the critical point can be accurately described within an effective theory involving only long-range collective fluctuations of the system. However, both in experiments and in numerical simulations, the infinite size limit cannot be reached and thus one observe deviations from the predicted thermodynamics limit behavior. The finite size scaling (FSS) ansatz has been developed precisely to infer the singular behavior (i.e., the exponents determining the universality classes) of the physical properties of a system in the thermodynamic limit, having only information on the system properties at finite sizes.

FSS has yet a more general validity and does not require the existence of a phase transition or an evolution process. Indeed, even though it was initially used to study finite systems near the critical point of the corresponding infinite system, FSS can be actually applied to describe structures that are self-similar when observed in a certain range of scales. As an example, we consider a Cantor set where we stop the procedure to divide intervals in three parts and removing the middle one at a scale $s_{0}=3^{-m}$. This corresponds to a fractal structure on scales between $s_{0}$ and 1, and to a non-fractal structure on scale smaller than $s_{0}$. If we measure the total length, $L(s)$, of the set with a stick of length $s=3^{-n}$ we find $L(s)=s^{1-D}F(s/s_{0})$ where $F(x)=1$ when $x>1$ whereas $F(x)=x^{1-D}$ when $x<1$ and $D=\log_{3}2$ is the Hausdorff–Besicovitch (or fractal) dimension of the Cantor set. Another illustration of FSS analysis is given by the truncated geometrical series $S(x,N)=\sum_{0}^{N-1}x^{n}$. When $x$ is close to 1 it is easy to see that $S(x,N)=t^{-1}F(tN)$, where $t=1-x$ and $F(z)=1-e^{-z}$. As a matter of fact, the FSS approach has been used to test scale invariance (and self similarity) also for non-critical systems such as (just to mention some very famous examples) polymers in confined geometries [23] and interfaces [24, 25]. In view of the above, FSS can also implemented on the well-established models of scale-free networks-like the Barabási-Albert model where the scale free behavior is not an emergent property at a critical point. Whether or not the same hypotheses holds for real world network does not undermine the possibility of applying FSS to them.

Employing the FSS machinery to test whether empirical networks display scale-free behavior in their degree distribution is not straightforward though. Unlike for physical systems, representations of a network at different scales are typically not available. Thus in order to test whether a network shows a power law distribution of the degree, we construct by hand smaller sized representations of it in an unbiased manner. We then use the characteristics of the large original network as well as the derived sub-networks to test the scale-free hypothesis. Figure 1 shows an illustration of this procedure for a snapshot of the structure of the Internet at the level of autonomous systems [22]. Subsection A provides a brief summary of finite size scaling applied to network topology. Subsection B presents an independent method of determining whether networks are scale-free based on analyses of the size dependence of the ratio of moments of the degree distributions. Subsection C provides information on the sampling scheme used to build sub-networks and on the region selected for the scaling analysis.

In the Results section we test the scale-free hypothesis (intended as the power law behavior in the degree distributions) on around two hundred large empirical networks (those considered in [19] and [20]). Remarkably, we find that such a venerable hypothesis cannot be rejected for many (but not all) networks. Moreover the two scaling exponents for such networks satisfy an additional scaling relationship, which derives from the shape of the degree cross-over in scale-free networks. We benchmark our results against the quality measure of the well-known scale-free graph introduced by Barabási and Albert [4]. Further we show that finite size scaling allows discerning pure power laws from log-normal and Weibull distributions. In conclusion, our results support the claim that scale invariance is indeed a feature of many real networks, with finite size effects accounting for quantifiable deviations.

A. Finite Size Scaling of networks

A scale-free network is postulated to have a degree distribution $p(k)\propto k^{-\lambda}$ beyond some lower degree cut-off $k_{min}$. For an infinitely sized network, since $k_{min}\geq 1$, the exponent $\lambda>1$ in order for $p(k)$ to be normalizable. In what follows, we will consider the cumulative distribution $P(k)=\int_{k}^{\infty}p(q)dq\propto k^{-\gamma}$ where $\gamma=\lambda-1>0$.

Networks are of course not infinitely large. In a network comprising $N$ nodes, $k$ can be at most equal to $N-1$. This is the intrinsic limit on $k$ given by the network size. Thus it is plausible that, below some $k_{c}$ (cross-over value), the degree distribution follows a power law behavior as would be expected for an infinite network but falls more rapidly beyond $k_{c}$. The finite size scaling hypothesis states that

[TABLE]

where $d<0$. The remarkable simplifying feature of the scaling hypothesis is that $P$ is not an arbitrary function of the two variables $k$ and $N$ but rather $k$ and $N$ combine in a non-trivial manner to create a composite variable. The behavior of the system is fully defined by the two exponents, $\gamma$ and $d$, and the scaling function $f$. The exponent $d<0$ so that, for an infinite size network ($N\to\infty$), the argument of $f$ approaches zero. A pure power law decay of $P(k,N)$ with $k$ for very large $N$ requires that $f(x)\to$ constant as $x\to 0$. The additional normalization condition is $f(x)\to 0$ sufficiently fast when $x\to 1$. The finite size effects are quantified by the behavior of the function $f$ as its argument increases, e.g., when $k\gtrsim k_{c}$. For a network with a finite number of nodes, the degree distribution does not follow a pure power law but is modified by the function $f$ (see also [26] for a discussion of finiteness in the context of growing network models).

A powerful way of assessing whether a network is scale invariant is to confirm the validity of the scaling hypothesis and determine the two exponents and the scaling function $f$ by using the collapse plot technique. One may recast Eq. (1) as

[TABLE]

Then the path forward is simple. For networks belonging to the same class but with different $N$, one optimally selects two fitting parameters $\gamma$ and $d$ by seeking to collapse plots of $P(k,N)k^{\gamma}$ versus $kN^{d}$ for different $N$ on top of each other [27]. The fidelity of the collapse plot provides a measure of self-similarity and scale-free behavior, the optimal parameters are the desired exponents, and the collapsed curve is a plot of the scaling function.

We start out with a single representation of an empirical network with $N$ nodes. For purposes of the scaling collapse plot, we seek additional representative networks of smaller sizes. In order to accomplish this, we obtained the mean degree distributions of multiple sub-networks of sizes $\frac{N}{4}$, $\frac{N}{2}$ and $\frac{3N}{4}$, which were then collapsed on to each other and the original network to create a master curve. The quality $S$ of the collapse plot is then measured as the mean square distance of the data from the master curve in units of standard errors. $S$ is thus like a reduced $\chi^{2}$ test, and should be around one if the data really collapse to a single curve and much larger otherwise [28].

Note that as a measure of the size of a network (or sub-network), one may use the number of nodes $N$ or alternatively the number of links $E$. The scaling function in this case reads as follows:

[TABLE]

where the exponent $\gamma$ is the same as before and the exponent $d_{E}<0$ ought to be equal to the previously introduced exponent $d$ for networks satisfying the finite size scaling hypothesis (see next section).

B. Ratio of moments test

A simple alternative and independent test of the scale-free hypothesis is to study the size dependence of the ratio between the $i$-th and the $(i-1)$-th moments of $k$, for various $i$. The $i$-th moment $\left\langle k^{i}\right\rangle$ is defined to be

[TABLE]

provided $i>\gamma$. Instead if $i\leq\gamma$, $\left\langle k^{i}\right\rangle$ converges to a constant value for $N\to\infty$. Therefore when $i-1>\gamma$,

[TABLE]

independently of $i$. Thus, for a scale-free network, a log-log plot of the ratio of consecutive moments versus $N$ is a straight line with slope $-d$. Likewise

[TABLE]

when $i>\gamma$, otherwise $\left\langle k^{i}\right\rangle$ goes to a constant for $E\to\infty$. Therefore when $i-1>\gamma$,

[TABLE]

The exponents $d$ and $d_{E}$ are not independent for scale-free networks. On the one hand, equations (4) and (6) imply $E\propto N^{d/d_{E}}$. On the other, in general $\left\langle k\right\rangle\propto E/N\propto N^{d/d_{E}-1}$. Due to the above equations $\left\langle k\right\rangle$ is constant for scale-free networks with $\gamma>1$, implying that $d=d_{E}$. Thus the difference between $d$ and $d_{E}$ values (that we statistically assess through their $Z$-score) provides an independent quality measure of the scale-free attributes of a network.

C. Sub-sampling and scaling region

In order to generate a sub-network of a given size $n<N$, we pick $n$ nodes at random among the $N$ nodes of the original network, removing all the other nodes and the links originating from them. It is well known that the sub-sampling procedure modifies the shape of the degree distribution of the network. In particular, sub-networks of scale-free networks are not scale-free because of deviations at low $k$ values [29] (this happens independently of the sampling scheme adopted [30]). The problem of the left tail of the distribution however applies more generally, because deviations from the scale-free behavior at low degrees are rather common in empirical and network models. Therefore we perform the scaling analysis described in subsections A and B only for $k\geq k_{min}$, where the lower bound of the scaling region $k_{min}$ is chosen such that the empirical distribution of the original network and its best power law fit (with exponent $\Gamma$, computed with the maximum-likelihood method of Clauset, Shalizi and Newman [18], see Methods) are as similar as possible above $k_{min}$ [31]. In the Supplementary Information we show that this allows us to get rid of any deviations induced by the sub-sampling scheme. However, when the empirical distribution of the network deviates substantially from a power law over its entire domain, then the estimated $k_{min}$ can become very large and may even diverge. In these cases the number of nodes $n^{*}$ of the (sub-)network with $k\geq k_{min}$ becomes very small or vanishing, yielding an unstable or undefined collapse. We thus use $n^{*}\geq\ln N$ as a condition on as the minimum number of nodes in each (sub-)network for the feasibility of the scaling analysis.

Results

To sum up, two independent statistical tests of the scale-free attributes of a network explained in subsections A and B are the quality of the collapse $S$ (i.e., the reduced $\chi^{2}$ between data and master curve) and the compatibility of $d$ and $d_{E}$ (measured through their $Z$-score). Figure S1 in the Supplementary Information outlines the flow of the analysis. In line with Broido & Clauset [19] and Voitalov et al. [20], we use these tests to define a classification for the degree distribution of empirical networks:

• •

SSF (strong scale-free) if $S\leq 1$ and $Z_{dd_{E}}\leq 1$,

• •

WSF (weak scale-free) if $S\leq 3$ and $Z_{dd_{E}}\leq 3$,

• •

NSF (non scale-free) otherwise or when $n^{*}<\ln N$ for the original network or any of its sub-networks.

Note the nestedness of the classification, for which a SSF network is also WSF.

Power law and Poisson distribution

We start analyzing the reference cases of Barabási-Albert [4] and Erdős-Rényi [33] models whose behavior is known. In the former case $p(k)\sim k^{-3}$, whereas, in the latter case $p(k)\sim\mbox{Poisson}_{\bar{k}}(k)$. Figure 2 shows that for a realization of the Barabási-Albert graph the degree distributions of the (sub-)networks result in a collapse of very high quality. The power law exponent $\gamma$ yielding the best collapse is consistent with the value $\Gamma$ obtained by maximum-likelihood fitting the degree distribution of the mother network with a power law [18]. Additionally, the moments ratio are indeed parallel lines, with compatible slopes $d$ and $d_{E}$. A more robust statistics is obtained by analysing $1000$ realizations of the Barabási-Albert model (Figure 3). Within this sample, $98\%$ of the networks are classified as SSF while $2\%$ as WSF. The estimated scaling exponents are all consistent with each others among the different realizations.

For the Erdős-Rényi model the estimated $k_{min}$ for the degree distribution is so large that it is not possible to have (sub-)networks with number of nodes $n^{*}\geq\ln N$ (in principle, for this network, the $k_{min}$ estimated from the KS test should be larger than the largest degree of the network). As such, the Erdős-Rényi graph is classified as NSF. We obtained the same outcome in an ensemble of 1000 realization of this network model.

Alternative fat tail distributions

While the power law is the only distribution featuring scale invariance, there are other distributions characterized by a fat right tail that can resemble a power law in finite systems. Hence determining which of these distribution better fits empirical network data is often a nontrivial task. In particular the classical approach based on $p$-values computed from a Kolmogorov-Smirnov test (see Methods) is able to rule out some competing hypothesis but not to confirm one [18]. Moreover, the hypothesis testing approach may fail when applied to regularly varying distributions [20]. It is therefore meaningful to put our finite size scaling approach to the test of alternative fat tail distributions. Here we consider the representative cases of the log-normal and Weibull distributions. The log-normal distribution $p(\ln k)=\mbox{Normal}(\mu,\sigma)$ is characterized by parameters $\mu$ and $\sigma$, respectively the mean and standard deviation of the variable’s natural logarithm. For large values of $\sigma$ this distribution is highly skewed and features a fat tail for large $k$ values. The Weibull distribution $p(k)=(h/l^{h})k^{h-1}\exp\left[-(k/l))^{h}\right]$ is characterized by parameters $h$ (shape) and $l$ (scale). The fat tail in this case appears for $h\to 0$. We use the Viger-Latapy algorithm [34] to generate networks with these degree distributions.

Figure 4 shows the scaling analysis for a realization of a network with log-normal $p(k)$ and for another realization with Weibull $p(k)$. In both cases we observe that the quality of the collapse is poor and that the moment ratios are not parallel lines. Therefore both networks are classified as NSF. Moreover, $S$ as a function of $\gamma$ does not show any minimum in the region around $\Gamma$ (the minimum does exist, but is located elsewhere). This means that the exponent estimated by finite size scaling $\gamma$ and that obtained from maximum likelihood power law fitting $\Gamma$ are substantially different: the outcome of the scaling analysis is not consistent in this case. However, the result depends much on the choices of parameters characterizing the distribution. Indeed Figure 5 shows that the percentage of networks classified NSF decreases by increasing $\sigma$ in the log-normal case, as well as by decreasing $h$ in the Weibull case – up to a point where the variance of the distributions becomes so large that the scaling analysis can hardly distinguish these distributions from power laws at finite $N$. For these cases, the value of $\gamma$ that minimizes $S$ is indeed compatible with $\Gamma$.

Real world networks

At last we move to real network data. We consider a large set of empirical networks taken from the Index of Complex Networks (ICON) as well as from the Koblenz Network Collection (KONECT). These are the datasets used by Broido & Clauset [19] and Voitalov et al. [20]. See the Methods section for a discussion on how we built the dataset. Overall, we have networks belonging to ten different categories: biological (PPI), social (i.e., friendship and communication), affiliation, authorship (including co-authorship), citation, text (i.e., lexical), annotation (i.e., feature, folksonomy, rating), hyperlink, computer, infrastructure. Figure 6 shows results of the finite size scaling analysis for selected network instances, whereas, Figure 7 and Table 1 summarize results of the scaling analysis for all the networks considered. The main outcomes of the analysis are the following.

• •

Figure 7(a): the scaling exponents $d$ and $d_{E}$ obtained from the moment ratio test are compatible in most of the cases.

• •

Figure 7(b): the value of $\gamma$ computed from finite size scaling is often in good agreement with $\Gamma$ obtained from the maximum likelihood power law fit of the degree distribution [18].

• •

Figure 7(c): the exponents $\gamma$ and $d$ of the scaling function are not independent but satisfy a universal relation $d\simeq-(\gamma+1)^{-1}$, which derives from the nature of the degree cross-over in scale-free networks – namely the maximum degree for which the power law behaviour holds. According to Eq. (1), this is the value $k_{c}$ for which the scaling function $f(x)\to 0$ (graphically speaking, when the master curve $P(k)k^{\gamma}$ falls down), corresponding to $x\gtrsim 1$ whence $k_{c}\sim N^{-d}$. The analysis presented in Figure 7(c) suggests that $k_{c}\sim N^{1/(\gamma+1)}$, and in agreement with theoretical results we find that also the maximum degree of the network $k_{max}$ scales in the same way (see Supplementary Information). However this scaling behavior is somehow different from the $k_{c}\sim N^{1/\gamma}$ as predicted by hand-waving argument [39, 40, 41], likely due to inner correlations in the networks which modify the value of the cross-over [40].

• •

No particular relation between quality of collapse $S$ and estimated exponent $\gamma$ is found, nor any clusterization of networks amenable to categories within the plane defined by these two variables (see Supplementary Information). However this result is obtained when the different network categories are well balanced in the dataset, because networks that are very similar tend instead to cluster together. This is for instance the case of protein interaction networks belonging to different species. In order to remove this artificial clustering effect, we have not considered in our dataset these (and other) cases of very similar networks nor repetitions of the same network (see Supplementary Information). This is the main reason why our dataset is apparently smaller than that used by Broido & Clauset [19].

• •

Overall, as shown in Table 1, the 185 networks of our dataset are classified as strong scale-free (SSF) in the 27% of cases, weak scale-free (WSF) for the 23% and non-scale-free (NSF) for 50%. This classification however does vary substantially among the different network categories. On the one hand, biological networks are very often classified at least as WSF. The same happens for computer and hyperlink networks, with outliers respectively given by the Gnutella peer-to-peer file sharing network (that has the same character of a social networks [42]) and by some hyperlink networks restricted to specific domains. Citation and text networks are few in our analysis, but are often scale-free. On the other hand, infrastructure networks (i.e., road and flights network) are rarely scale-free (with the notable exception of Air traffic control systems), possibly because of the heavy cost of establishing a connection. Between these two extremes, there are the social and other kinds of networks (see for instance the well-known discussion of the Facebook case presented in [43, 44], and that of other information sharing social network presented in [45]).

Discussion

Since the onset of network science, scale invariance of complex networks has been regarded as a universal feature present in real data [46, 47, 48, 49, 50, 18] as well as reproduced in models [4, 51, 52, 32, 53, 54]. Thus the recent claim by Broido & Clauset [19] that scale-free networks are rare created a stir, strengthening previous claims along the same direction [55, 18, 16]. Voitalov et al. [20] replied to these arguments fitting data to generalized power laws, that is, regularly varying distributions $p(k)=l(k)k^{-\lambda}$ (where $l(k)$ is a function that varies slowly at infinity and thus does not affect the power law tail). By allowing deviations from the pure power law distribution at low $k$, they argued that scale-free networks are definitely not rare. Gerlach & Altmann [21] very recently touched on this issue, showing that correlations present in the data can lead to false rejections of statistical laws when using standard maximum-likelihood recipes (in the case of networks, this can be important in the presence of degree-degree correlations).

In this work we go beyond statistical arguments and apply powerful tools from the study of critical phenomena in physics to analyse a wide range of model and empirical networks. Here we have showed that many of these networks spontaneously, without fine-tuning, satisfy the finite size scaling hypothesis, which, in turn, supports the claim that complex networks are inherently scale-free.

While a direct comparison with the results previously discussed would be interesting, the final results would be meaningless, given the differences in the underlying hypothesis of the different models. We showed how different hypothesis can lead to different results. The hypothesis underling our approach, which came from result previously obtained in the field of statistical mechanics and critical phenomena, goes beyond the applications they where initially thought for and it does not need the existence of a critical point. Together with the others, our methodology goes in the bag of tools a researcher can use in order to assess the scale freeness of a network.

Our scaling analysis is based on the extraction of small representations of the networks using a random node selection scheme. Of course, an intrinsic limitation of any rescaling method applied to network data is the impossibility to consider system sizes spanning orders of magnitude. As a further general remark, finding a robust method to rescale (or coarse grain [56, 57]) a network is still an open issue in the literature since networks are not embedded in any Euclidean space. Commonly used approaches lack generality since they are based on the choice of the embedding geometric space [58] or on the average path length [59]. In order to avoid ad hoc assumptions, we decided to follow the simplest (although not necessarily the most accurate) scheme. As shown in the Supplementary Information, by averaging over many extraction of the sub-network we are able to preserve the degree distribution of the original network, that is what we are interested in. Finally note our claims regards the self-similarity of the degree distribution, but we restrain ourselves in making general conclusions about the overall self-similarity of networks – this would involve the study of other quantities such as clustering, average path length and so on [60].

Materials and Methods

Here we report the steps to test the finite size scaling hypothesis of Eq. (2) together with the moments ratio test of Eq. (5). Note that in order to test Eqs. (3) and (7), one uses the number of edges $E$ ($e$) associated with each (sub-)network of size $N$ ($n$), and replaces $d$ with $d_{E}$.

Finite Size Scaling analysis

Given an undirected network of size $N$, our analysis is based on the following steps.

We compute the degree distribution $p(k,N)$ and use the method of Clauset, Shalizi and Newman [18, 31] to estimate the best fitting power law parameters $\Gamma+1$ and $k_{min}$. 2. 2.

We generate an ensemble of 100 sub-networks for each size $n\in\{\frac{N}{4},\frac{N}{2},\frac{3N}{4}\}$. Each sub-sample is obtained by picking $n$ nodes at random from the original network and by deleting all the other nodes and the links incident to them. We then compute the mean degree distribution $p(k,n)$ over each sub-network ensemble. 3. 3.

Both for the original network and for each sub-network, we check whether the (average) number of nodes $n^{*}$ with $k\geq k_{min}$ is larger than $\ln N$. If this condition is not met, we classify the network as non scale-free and the analysis ends. Otherwise, we proceed by removing the region below $k_{min}$ in both $p(k,N)$ and each $p(k,n)$, and renormalize them afterwards. As explained in the main text, this allows us to get rid of deviations at low degrees, including those induced by the sub-sampling (see also the Supplementary Information). 4. 4.

Using the moment ratio test, we determine $d$ (and its associated error) as follows. We compute a given moment ratio $\left\langle k^{i}\right\rangle/\left\langle k^{i-1}\right\rangle$ on each (sub-)network of size $n$, and use least-squares to fit $\ln(\left\langle k^{i}\right\rangle/\left\langle k^{i-1}\right\rangle)$ versus $\ln n$. We then average the resulting fit slope over different choices of the moments (indexed by $i$) to obtain $-d$. Note that since this test is computationally less expensive than the collapse analysis (see below), we use more than four sub-network sizes. In particular we use 20 equally spaced values of $n\in[\frac{N}{4},N]$, for each of which we compute the moments ratio (and associated error used as fit weight) over an ensemble of 100 $n$-sized sub-network built as described above. 5. 5.

For each (sub-)network size $n\in\{\frac{N}{4},\frac{N}{2},\frac{3N}{4},N\}$ we obtain the cumulative degree distribution $P(k,n)$. We then determine the exponents $\gamma$ and $d$ (and their associated errors) that maximizes the quality of the collapse plot (see below). Notably, the scaling exponent $d$ obtained from the collapse is always compatible with that obtained from the moment ratio test. Hence in order to decrease the computational cost of the method, one can in principle vary only $\gamma$ while keeping $d$ fixed at the value obtained from the moment ratios fit.

Quality of collapse

We now describe the procedure for deriving the master curve of the scaling function from the cumulative degree distributions of the various sub-networks, following the steps described in [28, 61]. The key premise is that when these distributions are properly rescaled they can be fitted by a single (master) curve. The quality of the collapse plot is then measured as the distance of the data from the master curve, and the collapse is good if all the rescaled distributions overlap onto each other.

In practice for each (sub-)network size $n\in\{\frac{N}{4},\frac{N}{2},\frac{3N}{4},N\}$ we have the set $\{j\}$ of ordered points for the cumulative degree distribution in the form $\{(k_{j},P(k_{j},n))\}_{j}$. After applying the scaling laws we have:

[TABLE]

so that $x_{nj}$ is the rescaled $j^{th}$ degree in the distribution of the $n$-sized sub-network, and $y_{nj}$ is the rescaled value of such distribution relative to the $j^{th}$ degree. We also assign an error on the latter quantity as $dy_{nj}=dP(k_{j},n)\,k_{j}^{\gamma}$, where $dP(k_{j},n)$ is the Poisson error on the count $P(k_{j},n)$ — see the Supplementary Information.

The master curve $Y$ is the function best fitting all these points. We define the quality of the collapse as

[TABLE]

where $Y_{nj}$ and $dY_{nj}$ are the estimated position and standard error of the master curve at $x_{nj}$, while $M$ is the set of terms of the sum (roughly, the set of points for which the curves for the various $n$ overlap).

For each $x_{nj}$, in order to define $Y_{nj}$ and $dY_{nj}$ we first need to select a set of points $m_{nj}$ as follows. In each of the other sets $n^{\prime}\neq n$, we select (and put in $m_{nj}$) the two points $j^{\prime}$ and $j^{\prime}+1$ that best approximate $x_{nj}$ from below and above, i.e., the two points such that $x_{n^{\prime}j^{\prime}}\leq x_{nj}\leq x_{n^{\prime}(j^{\prime}+1)}$. If this procedure fails to select two points for each $n^{\prime}\neq n$, then $Y_{nj}$ and $dY_{nj}$ are undefined at $x_{nj}$ which thus does not contribute to $S$ (this happens if set $n$ is alone in this region of $x$ and is the master curve by itself). Otherwise, we compute $Y_{nj}$ and $dY_{nj}$ using a linear fit through the selected points in $(n^{\prime},l)\in m_{nj}$, so that $Y_{nj}$ is the value of that straight line at $x_{nj}$ and $dY_{nj}$ is the associated standard error:

[TABLE]

[TABLE]

where $w_{n^{\prime}l}=1/dy_{n^{\prime}l}^{2}$ for the fit weights and $W=\sum_{(n^{\prime}l)\in m_{nj}}w_{n^{\prime}l}$, $W_{x}=\sum_{(n^{\prime}l)\in m_{nj}}w_{n^{\prime}l}x_{n^{\prime}l}$, $W_{y}=\sum_{(n^{\prime}l)\in m_{nj}}w_{n^{\prime}l}y_{n^{\prime}l}$, $W_{xx}=\sum_{(n^{\prime}l)\in m_{nj}}w_{n^{\prime}l}x_{n^{\prime}l}^{2}$, $W_{xy}=\sum_{(n^{\prime}l)\in m_{nj}}w_{n^{\prime}l}x_{n^{\prime}l}y_{n^{\prime}l}$, $\eta=WW_{xx}-W_{x}^{2}$ for the fit parameters.

The quality of the collapse $S$ measures the mean square distance of the sets to the master curve in units of standard errors, analogously to a $\chi^{2}$ test [28]. The number of degrees of freedom can be estimated by noting that each of the $|M|$ points of the sum of $S$ has in turn 3 intrinsic degrees of freedom: $|m|$ points as described above (6 in our case) minus 2 from computing mean and variance of $Y$, minus 1. Hence by using $3|M|$ as normalization factor, $S$ should be around one if the data really collapse to a single curve and much larger otherwise.

We optimize the quality $S$ of the collapse by varying the scaling exponents $\gamma$ in the interval $\Gamma-0.5\leq\gamma\leq\Gamma+0.5$ and $d$ in the interval $d-0.1\leq\gamma\leq d+0.1$. The errors associated with $\gamma$ and $d$ are estimated with a $S+1$ analysis: $\Delta\gamma$ is such that $S(\gamma+\Delta\gamma)=S(\gamma)+1$ and $\Delta d$ is such that $S(d+\Delta d)=S(d)+1$.

Dataset

We extract a collection of real network data from the Index of Complex Networks (ICON) at https://icon.colorado.edu/ as well as the Koblenz Network Collection (KONECT) at http://konect.uni-koblenz.de/. The full list of networks we consider together with detailed results of the finite size scaling analysis are reported in the Supplementary Dataset Table. To define the dataset we select networks (removing duplicates appearing in both ICON and KONECT) according to the following criteria.

First, to allow for a reliable scaling analysis, we only use networks with $N>1000$ and $E>1000$ (for computational reasons, we did not consider networks with more than 50 million links). We then include undirected networks, as well as the undirected version of both directed and bipartite networks. Similarly, we consider binary networks as well as the binarized version of weighted and multi-edge networks. We however ignore networks that are marked as incomplete in the database. Importantly, among database entries that possibly represent the same real-world network we select only one (or at most a few) entry, and consistently we do the same for temporal networks (when there is only one snapshot, we ignore the time stamp of links).

In practice, in KONECT we select only the Wikipedia-related networks in English language. For ICON the implications are more profound. We ignore interactomes of the same species extracted from different experiments, the (almost 100) fungal growth networks, the (more than 100) Norwegian boards of directors graphs, the (more than 100) CAIDA snapshots denoting autonomous system relationships on the Internet, networks of software function for Callgraphs and digital circuits ITC99 and ISCAS89. We consider only one instance of Gnutella peer-to-peer file sharing network, as well as a few instances of the (more than 50) within-college Facebook social networks and of the (about 50) US States road networks. Among the (more than 100) KEGG metabolic networks, we select 17 species trying to balance the different taxonomies.

Thus, in our analysis, we do employ the same data source used by Broido & Clauset [19], but we avoid over-represented network instances. As explained in the main text, this procedure removes the clustering of similar networks shown in Figure 6, and leads to less biased conclusions on the scale-free nature of networks belonging to different categories.

Acknowledgments.

GCi and GCa acknowledge support from the European Project SoBigData++ (GA. 871042). AM acknowledges support from University of Padova through ”Excellence Project 2018” of the Cariparo foundation.

Author Contribution.

GCa conceived the experiment. MS and GCi performed the analyses of the dataset. AM coordinated the activities on finite size scaling analysis. All authors contributed to the interpretation of results, and to the writing of the manuscript.