Scaling theory of fractal complex networks

We show that fractality in complex networks arises from the geometric self-similarity of their built-in hierarchical community-like structure, which is mathematically described by the scale-invariant equation for the masses of the boxes with which we cover the network when determining its box dimension. This approach—grounded in both scaling theory of phase transitions and renormalization group theory—leads to the consistent scaling theory of fractal complex networks, which complements the collection of scaling exponents with several new ones and reveals various relationships between them. We propose the introduction of two classes of exponents: microscopic and macroscopic, characterizing the local structure of fractal complex networks and their global properties, respectively. Interestingly, exponents from both classes are related to each other and only a few of them (three out of seven) are independent, thus bridging the local self-similarity and global scale-invariance in fractal networks. We successfully verify our findings in real networks situated in various fields (information—the World Wide Web, biological—the human brain, and social—scientific collaboration networks) and in several fractal network models.


Revising paradigms of fractal complex networks
It will soon be two decades since it was first shown that some real networks (such as the World Wide Web [WWW] and different biological networks) have fractal properties 1,2 .This means that, when covered with non-overlapping boxes, with the maximum distance between any two nodes in each box less than l B , they exhibit power-law scaling [1][2][3][4][5] : where N B (l B ) is the number of boxes of a given diameter, and d B is the fractal (or box) dimension of the network of size N.Such fractal networks are also said to be self-similar, because their power-law degree distributions, remain invariant under a renormalization scheme 6,7 , according to which a new network emerges from the original one when nodes belonging to the same box in the original network are replaced by one supernode in the renormalized network.In this case, the supernode is connected to another supernode if in the original network there is at least one link between the nodes of the corresponding boxes.
Here, at least two critical remarks can be made.The first remark is that an analogous invariance of the degree distribution with respect to the box-covering renormalization scheme is also observed in networks that do not satisfy Eq. (1) (in this respect, well-known examples are the internet and Barabási-Albert (BA) networks 2,8,9 ).The second remark is that it is not entirely clear, what structural characteristics of fractal networks exhibits geometric self-similarity and remain invariant 10 under the described renormalization.Clearly, the power-law node degree distribution cannot be considered such a characteristic because it is intrinsically invariant under the rescaling of the degree 11 .Its invariance under box-covering renormalization may only suggest the existence of some (presumably) degree-dependent network measure, whose self-similarity under the renormalization procedure could result in the observed invariance of the degree distribution.One argument supporting this statement is that random networks, where the degree distribution is not a power law, can also exhibit fractal properties (in this regard, the best example is the giant component of classical random graphs near the percolation transition).
If the above remarks, indicating an incomplete understanding of fractality in complex networks, are reasonable, pertinent questions would be: What are the real origins and potential consequences of fractality in complex networks?What determines networks' fractal dimension?Indeed, several studies have been published throughout the years that focus on the exploration of the origins of fractality [12][13][14][15][16][17][18][19] .However, these efforts did not lead to consensus.Thus, there is a lack of realistic (and not just deterministic 20,21 , or reflecting the renormalization procedure 2,22,23 ) fractal network models that would allow testing the role of fractality in the context of geometry-involving issues 24 , such as navigability, localization of information sources, prediction of hidden network connections, etc.These are of particular importance when faced with the confirmed fractal properties of different information, biological, and even social networks (see e.g. 25,26 .The goal of this article is to initiate far-reaching changes in this state of affairs. In what follows, we will first argue that the correct scale-dependent network measure, which is self-similar (i.e.geometrically invariant) under the l B -box-covering renormalization procedure, is the normalized mass of the boxµ(L, k) = m(L, k)/⟨m⟩, where m(L, k) is the number of nodes in the box of diameter L ≥ l B and hub degree k, and ⟨m⟩ = N/N B (L) = L d B (1) is the average mass of non-overlapping boxes of this diameter.It should be emphasized here that although the definition of the box is the same throughout the paper, we distinguish between l B -boxes used to renormalize the network and L-boxes (where L ≥ l B ) whose self-similarity we examine.This distinction is crucial to make it easier to understand the main idea of the paper.
Then, we show that one of the consequences of this result is the previously discovered scaling relation between the degree k ′ of the supernode in the renormalized network and the degree k of the hub of the corresponding l B -box in the network before renormalization: , where d k is only one of four scaling exponents that characterize microscopic structure of the fractal complex network and determine its box dimension.We also show that if the fractal complex network has a power-law node degree distribution (which is traditionally referred to as the scale-free property), then the mass box distribution also follows the power-law, and it is invariant under the box renormalization procedure.Furthermore, the characteristic exponents of both distributions are related to the microscopic scaling exponents describing the masses of the boxes, thus bridging local self-similarity and global scale invariance in fractal complex networks.Lastly, we successfully verify our findings in real networks situated in various fields (information -the World Wide Web, biological -the human brain, and social -scientific collaboration networks) and in several fractal network models.

Geometric self-similarity
In classical fractals 27 , which reproduce themselves at different space scales, self-similarity manifests itself in the scale-invariant equation 10 , which describes how the mass m(L) of the system changes with its linear size L: where b > 0. In theoretical physics, this type of equation is, for example, encountered in the theory of critical phenomena 11,28 .Mathematically, this equation defines a homogeneous function.Its solution is simply a power law: which, in the case of fractals, determines their fractal dimension, d f = ln µ/ ln b, and leads to the well-known scaling relation 29 : Moving forward, to address the problem of geometric self-similarity in complex networks, we first argue that Eq. ( 1) can be treated as a special case of Eq. ( 5).Then, building on this observation, we assume that Eq. ( 5) is also a special case of a more general equation in which the masses of the system and its parts, which are further identified with the number of nodes in the network and the number of nodes in different L-boxes extracted from this network, respectively, do not only depend on the diameter of the examined set of nodes (i.e. the entire network or a box) but also on the degree of the best-connected node in this set.This assumption leads us to the consistent scaling theory of fractal complex networks.
To grasp the relation between Eqs. ( 1) and ( 5), it is enough to analyse the meaning of Eq. ( 5), which can be interpreted in two ways.More directly, it states that if one considers a smaller part of the system, let's say of size L ′ = bL (with b < 1), then m(L ′ ), as compared to m(L), is decreased by a factor µ(b) = b d f , which only depends on b.However, this equation also applies to the masses of the system on two different scales, or resolutions, which, from a formal point of view, can be treated as two stages of some renormalization procedure applied to that system.(A network-based illustration of these two interpretation schemes is shown in Figs. 1 and 2(a-c).)Accordingly, to make Eq.( 5) more operationalizable, it can be rewritten as: where the notation with the apostrophe is introduced to indicate the relation between the mass of the system before renormalization, m(L), to its mass after renormalization, m ′ (L ′ ).Now, it is easy to see that Eq. ( 1) is indeed a special case of Eq. ( 6), with: , where L and L ′ stand for diameters of the network before and after renormalization, respectively.Figure 1.Schematic illustration of the idea of geometric self-similarity in complex networks on the example of the fractal model of nested BA networks (for the definition of the model, see "Methods" section).Part a) of the figure shows that the network can be subdivided into parts-boxes of a given diameter-each of which is (at least approximately) a reduced-size copy of the entire network.In the top picture shown, one such box, marked in red, is extracted from the original network and treated as a new network (shown below).It is divided again into new smaller boxes, some of which are marked with different colours.Both macroscopic and microscopic characteristics of this new network (represented by green squares in Fig. 2) are similar to those of the original network (indicated by navy circles in Fig. 2).Part b) of this figure illustrates renormalization procedure applied to the same network as in part a.The top original network is divided into boxes of a fixed diameter, some of which are marked with different colours.In the new network after renormalization (shown below), these boxes are replaced by nodes with the corresponding colours.Again, the macroscopic and microscopic characteristics of the network after renormalization (represented by red triangles in Fig. 2) are similar to those of the original network.
In what follows, to extend the concept of geometric self-similarity to fractal complex networks, we assume that Eq. ( 6) can be rewritten in the form: where d B is the box dimension of fractal networks, whereas m(L, k) and m ′ (L ′ , k ′ ) stand for the number of nodes and supernodes in the same box, before and after its renormalization with boxes of diameter l B < L, respectively.In other words, in Eq. ( 7), m ′ (L ′ , k ′ ) is equal to the number of l B -boxes used to cover the initial box of mass m(L, k).As indicated in this equation, during renormalization, when l B -boxes are replaced with supernodes, not only the mass of the initial box changes, but also its diameter (from L to L ′ ) and the degree of its hub (from k to k ′ , where k ′ is the degree of the best-connected l B -box within the initial L-box).Now, since Eq. ( 7), like Eqs. ( 5) and ( 6), defines a generalized homogeneous function 28 of the form: after its substitution into (7), we obtain several scaling relations characterizing fractal networks.The first relation poses: where d L = 1 is a direct consequence of the applied renormalization procedure, assuming perfect tiling of the network with boxes of diameter l B each 2 .The second relation has the form of the long-confirmed empirical relation 1 , but in the context of Eq. ( 7), which applies to boxes of any diameter L ≥ l B , the range of its applicability is much wider than previously thought, which in view of our approach was limited to the case of L = l B and L ′ = 1.Finally, taken together Eqs. ( 7)- (10) give the following scaling relation: which is one of the most important results of this article.According to Eq. ( 11), the box dimension d B of fractal networks is only determined by the scaling exponents characterizing the microscopic structure of the network at the box level.In particular, as follows from Eq. ( 7), beside the method of determining the box dimension, which involves counting non-overlapping boxes (see Fig. 2(a)), d B can also be obtained by subsequent renormalizations of a single L-box with smaller boxes of a given diameter l B < L (see Fig. 2(d)).The decreasing sequence of renormalized masses of the box obtained in this way: m, m ′ , m ′′ , . . ., m (i) , m (i+1) . . ., when presented on a graph in the form of points (m (i) , m (i+1) ), can be used to determine the coefficient l −d B B in Eq. (7).Repeating this procedure for different values of l B , a set of points (l B , l −d B B ) can be obtained which, when fitted with straight line on a double logarithmic scale, gives the same value of d B , as that resulting from the classical method based on Eq. ( 1).
The form of Eq. ( 11) is also very suggestive.It is the sum of two components, each of which is the product of scaling exponents relating to specific quantities characterizing the mass of the box before and after renormalization.In classical fractals,

4/18
in which the mass of the box depends only on its linear size L, this sum has only one component.For this reason, in classical fractals, the fractal dimension can be determined by one of two methods: the box-covering method or the cluster-growing method, which are equivalent to each other.However, this is not the case of fractal complex network 1,30 , where α, playing the role of the spreading dimension, only describes how the mass of the box, Eqs. ( 8), varies with its diameter: where b > 0. In a similar vein, the second addend in (11), which is further called the mass exponent (in analogy to the degree exponent, d k ), only characterizes, how the local network density (understood as the number of nodes per local area of dimeter L = L ′ ) changes as a result of renormalization: Finally, an observation of great importance for the scaling theory of fractal complex networks (see "Scale-free property" section) is that when the box mass m(L, k) (8) is divided by the average mass ⟨m⟩ = N/N B (L) = L d B (1) one gets the normalized mass: which turns out to be the invariant of the l B -renormalization procedure, since Indeed, the normalized box mass ( 15) is a local network measure that behaves the same regardless of the scale of observation, as the following scaling relations clearly describe: and A proper perspective on the meaning of these two relations is gained when comparing them with the corresponding relations for classical fractals, namely Eqs. ( 5) and (6).From this perspective, the scaling exponent d m appears to be the self-similarity dimension of fractal complex networks, which, remarkably, is different from the box dimension d B .

Scale-free property
At this point, we would like to emphasize the lack, in our considerations so far, of scale-free node degree distributions, whose invariance due to the renormalization procedure is considered an attribute of fractal networks [1][2][3] .Interestingly, this lack clearly shows the otherwise obvious fact that fractal networks may not have the scale-free property.Nevertheless, when they reveal the property, then both the node degree distribution, P(k) ∼ k −γ , and the box mass distribution, P(m) ∼ m −δ , are invariant under the box-covering renormalization procedure, with their invariance being a consequence of the already discussed geometric self-similarity of boxes and the scale-free property of the distribution of normalized masses P(µ) ∼ µ −δ , from which P(m) inherits its characteristic exponent δ (see Fig. 2(b,c)).
To show this, let us assume that P(µ) is scale-free: where, by writing P(µ; L) instead of P(µ), we emphasize that all boxes in the network have the same diameter L. To clarify, this distribution refers to the normalized masses of non-overlapping boxes of diameter L used to cover the network.Invariance of this distributions in networks after l B -renormalization is due to the property ( 16) which implies that P(µ; L) = P(µ ′ ; L ′ ), i.e.
Now, having the relationship between µ, m and k (15) and using it together with Eq. (19) in the balance equations between the corresponding distributions, i.e.P(µ)dµ = P(m)dm and P(µ)dµ = P(k)dk, it is easy to show that and where the characteristic exponent γ is given by: The invariance of these distributions in networks after l B -renormalization is obvious due to Eq. ( 20).
The above reasoning shows that the geometric self-similarity of the boxes ( 16)-( 18) and the scale-free distribution of their normalized masses (19) themselves guarantee the invariance of P(k) and P(m) under renormalization.Another consequence of these two assumptions (i.e.self-similarity and scale-freeness), which is not obvious, although it may seem so at first glance, is independence of P(µ; L) (19) from the diameter of the boxes L (of course, the same applies to P ′ (µ ′ ; L ′ )).In general, this feature can be shown to be true by comparing the numbers of boxes having the same hub nodes when the network is covered with boxes of different diameters.Because the diameter and degree of the hub determine the mass of the box, such a comparison comes down to comparing the number of boxes with the given normalized masses: where the relationship between the considered masses is determined by Eq. ( 17).Making the appropriate substitutions in this equation, i.e. 1), ( 15) and (19), not only do we confirm that the distribution P(µ; bL) is scale-free regardless of b ( 19), but we also obtain a new relation between the scaling exponents: Interestingly, using Eqs.( 13) and ( 23), the above relation can be easily transformed into the well-known relation 1 At this point, a natural question to ask is: How is it possible that the relation (26) was originally derived without having to refer to the self-similarity of the boxes?In fact, as we show below, self-similarity cannot be ignored, and the derivation described in Ref. 1 takes it into account, albeit implicitly.
More precisely, in the original reasoning leading to Eq. ( 26), one starts with the following equation: where is the number of nodes with k (respectively, k ′ ) links in the network before (after) renormalization.Then, substitutions are made in this equation: (1) stand for the number of nodes in the network before and after renormalization, respectively.These substitutions lead to the following density balance equation: from which Eq. ( 26) is obtained under the assumptions that both node degree distributions are scale-free with the same scaling exponent, i.e.P(k) ∼ k −γ and P ′ (k ′ ) ∼ k ′−γ , and that Eq. ( 10) is met between k and k ′ .It should be emphasized, however, that what underlies validity of these assumptions is the geometric self-similarity of the boxes and the scale-free distribution of their masses.Furthermore, this derivation itself is a special case of more general considerations, in which the starting point is the below equation: whose logic is similar to that behind Eq. ( 24).To explain, the left-hand side of Eq. ( 29) represents the number of boxes with diameter L and hub degree k in the network before renormalization, while the right-hand side is the number of boxes with the same diameter L and hubs of degree k ′ in the network after l B -renormalization.The numbers of these boxes must match because hubs of degree k ′ in the network after renormalization arise from those l B -boxes in the network before renormalization that contained hubs of degree k.The relation between the masses of the considered boxes is given by Eq. ( 18).Interestingly, for arbitrary value of L, the scaling analysis of Eq. ( 29) leads to the scaling relation (25).However, when L = 1 is assumed, then Eq. ( 29) can transformed to Eq. (27).In particular, the left hand side of Eq. ( 29) becomes: where we one by one used Eqs.( 1), ( 19), (15), and (23).Similar transformations applied to the right-hand side of Eq. ( 29) result in: what was to be shown.

From microscopic to macroscopic scaling exponents in real and model-based fractal networks
All scaling exponents discussed in this article, which describe fractal complex networks, can be divided into two groups.The first group refers to the macroscopic characteristics of the network (d B , γ, and δ ), and the second group includes the exponents that characterize the network structure at the microscopic level (d k , d m , α and β ).Interestingly, exponents from both groups are related to each other and, as in the scaling theory of critical phenomena, only a few of them, three to be exact, are independent.The choice of the three fundamental exponents depends on the focus of the study.Here, to validate our results in real and model-based fractal networks, we take the easier to measure macroscopic exponents as independent.This choice results in the following set of test relations, cf.Eqs. ( 25) and ( 26): and, cf.Eqs. ( 11) and ( 13): of which only the relation for d k (30) has been verified in real 1 and model 2 networks, and the results of the validation of relations (31) are summarized below.
The real networks analyzed in this paper come from various fields and represent information, social, and biological networks.We analyzed: 1) a sample of the WWW with nodes corresponding to web pages and links standing for hyperlinks 31 ; 2) a coauthorship network (DBLP), where nodes are scientists and edges are placed between two scientists if they have co-authored a paper 32,33 ; 3) a functional brain network, which reflects the correlation between the activity of different areas in the human brain 34,35 .In addition to real networks, we have also analyzed several fractal network models, including our own network model, which is based on nested BA networks 36 , the Song-Havlin-Makse (SHM) model 2 and (u,v)-flowers 20 .Detailed information on all these networks (real and synthetic) can be found in "Methods" section.

Table 1.
Values of the scaling exponents for various fractal networks.In the table, N is the number of nodes in the analyzed network, ⟨k⟩ is the average node degree, and d corresponds to the diameter of the network.The empirical values of the scaling exponents were determined by fitting the appropriate scaling relations to real data and results of numerical simulations, for real and model networks, respectively.The theoretical values of microscopic exponents, which are given in brackets, are of two types.For real networks and for the nested BA model, they were calculated on the basis of appropriate scaling relations using empirical values of macroscopic exponents, and for deterministic models of fractal networks, they were calculated on the basis of derived theoretical relations.Table 1 presents the theoretical and empirical values of the scaling exponents of all analyzed networks.The theoretical values, which are given in brackets, are of two types.For the deterministic model-based networks-the SHM model and (u,v)-flowers-their values can be calculated using the appropriate formulas, the details of which are provided in "Methods" section (more precisely: "SHM model" and "(u,v)-Flowers" subsections, respectively).For real networks and for the numerical model of nested BA, the theoretical values of α and β were calculated from Eqs. (31) using the empirical values of the macroscopic exponents.Correspondingly, the empirical values of the scaling exponents were calculated from Figs. 3 and 4 according to the following protocol (the same for each network): First, we determined the box dimension d B of these networks resulting from tiling the network with boxes of different sizes l B , see Figs. 3, 4 (a-c).To this end, we used the algorithm developed by Song et al. 37 , and in the case of deterministic models of fractal networks, shown in Fig. 4, we additionally analysed the tiling consistent with their deterministic construction procedures, finding that they use a much smaller number of boxes than the Song method.We confirmed that the value of d B after renormalization (even multiple times) remains the same as before renormalization, see Fig. 2(a).We then examined the invariance of distributions P(k) and P(µ).The given values of l B refer to the diameter of the boxes that are used to renormalize the network.As already stated, l B = 1 refers to the original network -before renormalization.In the case of P(µ) distributions, the diameters L of the boxes whose mass was studied are also given.With respect to these distributions, the provided values of l B and L should be read as follows: The relevant distribution P(µ) refers to the network that was first renormalized with boxes of diameter l B and then covered with non-overlaying boxes of diameter L. Regarding P(µ), however, due to the low statistical reliability of the data for l B , L > 1, in this paper, we only present data for the largest networks (i.e.WWW and model based networks).It should be noted that in all networks we studied, both distributions are scale-invariant, with well-defined characteristic exponents γ and δ (see Figs. 3, 4 (d-i)).Lastly, having determined the macroscopic scaling exponents: d B , γ, and δ , we were able to calculate the theoretical values of the local exponents-α and β , Eqs. ( 31)-which we used to obtain the adequately rescaled masses of boxes to determine their empirical values (see Figs. 3, 4 (j-l)).In particular, to obtain the empirical value of α, the masses of all the internally connected boxes, obtained during tiling the network with different l B -boxes, were divided by the hub's degree raised to the power of the theoretically obtained β .Such rescaled masses m/k β were then plotted against the actual diameters of the boxes, L < l B , which had been specified individually for each box.A similar procedure was applied to determine the empirical value of β .(For more details see the subsection titled "Numerical calculation of microscopic scaling exponents" in the "Methods" section).
Interestingly, in the case of deterministic fractal network models, only the box-covering method which takes into account the network construction procedure while using a smaller number of boxes than Song's method, leads to microscopic exponents consistent with their theoretical predictions (cf.Fig. 4(k,l) and Tab. 1).In the case of these networks, the poor performance of Song's method is especially visible in the range of small masses of the P(µ) distributions (see subset graphs in Fig. 4(h,i)).We suspect, this latter observation may explain the occurrence of two different scaling behaviours (for small and large µ) in P(µ) of other fractal networks (cf.Figs.3(g-i) and 4(g)).

Perspectives
The origins and consequences of fractality are one of the three main research directions in the geometry of complex networks 24 , next to the hyperbolic geometry of hidden network spaces 38,39 and the geometry induced by dynamic processes in networks [40][41][42] .Although these three geometries, due to the various definitions of distance in each of them, are defined differently, there is no doubt that they must be closely related to each other.While these relationships have yet to be explored, evidence of their existence can be found in our results.
For example, when examining deterministic models of fractal networks (SHM model and (u,v)-flowers, see "SHM model" and "(u,v)-Flowers" subsections, respectively), we noticed that while macroscopic scaling exponents are very stable in the sense that they do not depend on the box-covering method 37,43 , this may not be the case for microscopic exponents.In particular, in the mentioned models, gathering nodes according to their kinship-which is the most optimal, because it corresponds to the smallest number of boxes-gives the values of microscopic exponents closest to their theoretical predictions.Since the degree of kinship can be thought of as a distance in some metric space-the space of kinship-this observation is important.In fact, the fractality of these models may be considered a feature they inherit from their kinship spaces.Here, natural questions arise, such as whether the fractality of real complex networks may result from the properties of hidden (similarity-based) metric spaces 44 .Similar studies on community structure confirm the existence of such a relationship [45][46][47] .The mention of the community structure is not entirely accidental here, because, as the example of the DBLP network shows-in which the removal of weak ties reveals its fractal properties (see also 26,33 )-the fat-tailed community size distribution 48,49 may result from the scale-invariant distributions of box masses observed in (not necessarily tree-like) fractal skeletons 13,14 of these networks.
The second thread that we would like to emphasize concerns the geometry induced by diffusion-like dynamic processes in networks [40][41][42] .In classical fractals, this kind of geometry is closely related to the cluster-growing method of calculating their fractal dimensions, which is actually a way of measuring the distance 27 .In complex networks, establishing an analogous relationship has not been possible so far due to the lack of theoretical foundations distinguishing between the box dimension d B (7) (which can be determined by the box covering method) and the spreading dimensionα (12) (which corresponds to the cluster-growing method).It seems that the scaling theory of fractal complex networks presented in this paper has the potential to break this impasse.This is even more likely since in its general findings, with box masses depending not only on the diameter of the boxes but also on the degree of the best-connected node inside the box, the theory refers to the well-established heterogeneous (degree-based) mean-field theory commonly used to study dynamical processes on complex networks 50 .

Real and model-based fractal networks analysed in the paper
The real networks analysed include: • WWW network: The web subset analysed consists of 326 k web pages that are linked if there is a URL link from one page to another 31 .This dataset has been analysed for fractal properties in many other papers (see e.g. 1,2 .It is publicly available in many network repositories (e.g. 52).
• DBLP coauthorship network: DBLP is a digital library of article records published in computer science 32,53 .In this study, we use the 12th version of the dataset (DBLP-Citation-network V12; released in April 2020, which contains information on approximately 4.9 M articles published mostly during the last 20 years).We ourselves processed the raw DBLP data into the form of coauthorship network, from which we extracted the network backbone by imposing a threshold on the minimum number of joint papers (≥ 25) two scientists have.This procedure significantly reduces the size of the studied network (from 2.9 M nodes and 12.5 M links to 2.5 k nodes and 3.2 k edges), but thanks to it the network becomes naturally fractal.
• Human brain networks: The networks are based on functional magnetic resonance imaging (fMRI).The fMRI data consists of temporal series, known as the blood oxygen level dependent (BOLD) signals, from different brain regions.To build brain networks, the correlations C i j between the BOLD signals are calculated and the two nodes (brain regions) are connected if C i j is greater than some threshold value T .In our case we assume T = 0.85.The brain networks analysed here were used in 34,35 and can be found at 54 .
The studied models of fractal networks include: • SHM model: The details of the model are presented in "SHM model" subsection, where local scaling exponents for this model were also derived.
• (u,v)-Flowers: The details of the model are presented in "(u,v)-Flowers" subsection, where local scaling exponents for this model were also derived.
• Nested BA networks: The nested BA network model has three parameters: N -the number of nodes, k max -the degree of the best connected node in the network, and m -the number of edges by which the newly created node connects to the already existing nodes.The network evolution procedure is as follows: 1. First, a BA network with the hub of degree k max is created (that is, the network grows until one of the nodes reaches degree k max ).
2. Then, as long as the size of the network is less than N (see Fig. 5): (a) a node is chosen proportionally to its degree k and it is replaced with a BA subnetwork with the largest node degree k; (b) edges that were connected to the removed node are reconnected to randomly selected nodes of the newly created subnetwork.

Numerical calculation of microscopic scaling exponents
In Figs 3(j-l) and 4(j-l) we presented the microscopic scaling exponents α and β .Their values result from fitting a straight line to the set of points marked with blue circles and red triangles, respectively.These points represent the geometric mean of logarithmically equal sized bins of the original data which were obtained as follows.First, we estimated the range of l B for which dependence of N B on l B is linear in log-log scale.In Figure 6(a), which shows an example analysis of nested BA model, this range is indicated by a gray rectangle.Then, for each l B in this range, we performed box covering, obtaining a triple of values (m, L, k) for each box.Based on the set of such triples and the theoretically calculated value of the α or β exponent, a set of points (k, m/L α ) or (L, m/k β ) was created, respectively.In Fig. 6(d), the later set is represented by yellow circles.Having this raw data, logarithmic binning has been performed and geometric mean has been calculated for each bin.Fig. 6(d), in addition to the blue points denoting the geometric mean, also shows the geometric standard deviation and the fitted line, whose slope corresponds to the α value we are looking.
If we restrict our analysis to one specific value of l B (for example, in Fig. 6(b) we take l B = 32, while in Fig. 6(c) we take l B = 8), we end up with a much smaller set of triples that would not allow for a reliable results.These limitations are particularly severe for small real, DBLP, and brain networks, with sizes N < 3 • 10 3 each.The blue lines in Fig. 6(b) and 6(c) are not a result of fitting but are shown for comparison purposes only.

Microscopic scaling exponents for deterministic fractal network models
In this section, we derive exact formulas for microscopic scaling exponents (α and β ) characterising deterministic fractal network model.

SHM model
In the Song-Havlin-Makse (SHM) model2 , at t = 0, the network starts to grow from two nodes connected by one link.Then, during subsequent, t + 1, time steps, next (t + 1)-generations of the network recursively emerge, in which: s new nodes are attached to the endpoints of each link of the previous t-generation, old links are removed from the network, and new links are created in place of those removed, which connect pairs of offspring-nodes attached to the endpoints of the deleted ones.
As a result of this construction procedure, in successive generations, every node i increases its degree multiplicatively: where is the time that has elapsed since the node appeared in the network for the first time, at time t i , assuming that its initial degree1 was equal to k i (t t ,t i ) = 1.
It is also easy to see that a similar multiplicative dynamics is also shown by the diameter L i and mass M i of the largest box where the i-th node of degree k i (t,t i ) acts as a hub 2 .Specifically, the maximum diameter of such a box is given by: with3 a = 3 and L i (t i ,t i ) = 1, whereas its mass satisfies the following recurrence relation: where4 n = 2s + 1 and M i (t i ,t i ) = 1.At this point, it is worth noting a few remarks regarding Eqs. ( 34) and (35).First, Eq. ( 34), when applied to the largest boxes with the oldest nodes, i.e. those from which the network's evolution began, shows how the diameter of the entire network changes in subsequent generations: Correspondingly, by applying Eq. ( 35) to these boxes, one finds the analogous formula for the total number of nodes in the network: In Ref. 2 , the above recurrence relations, Eqs. ( 36) and ( 37), together with Eq. ( 32), were used to derive exact expressions for the box dimension of the considered fractal network model, cf.Eq. ( 1): and for its degree exponent, cf.Eq. ( 8): thus enabling verification of the scaling relation (17), according to which, in this model, the characteristic exponent of the degree distribution is given by: The second remark regarding Eqs.(34) and (35) it that boxes containing hubs of degrees k i (t,t i ) > 1 may have diameters and masses smaller than L i (t,t i ) and M i (t,t i ), respectively.For example, when the diameter of the i-th box is l i = 1, then the box is confined to the node itself and as a result its mass is equal to m i = 1.Similarly, when l i = a = 3, then the box, apart from the hub itself, also contains all its neighbours, making the mass of the box equal to m i = 1 + k i ≃ k i .More generally, the diameter of the i-th box can be equal to: where with the value of τ affecting its mass, which can be determined from 5 : Now, substituting Eq. ( 32) into ( 43), one gets: Then, using Eq. ( 41) in (44), one obtains the following relation for the mass of the box as a function of its diameter and hub's degree, cf.Eq. ( 9): where the local scaling exponents are given by: and It is easy to see that, together with the previously obtained expressions for d B (38) and d k (39), the obtained above expressions for α and β satisfy the scaling relation: d B = α + d k β , Eq. ( 10), which has been derived in the main text of the paper.

(u,v)-Flowers
In the deterministic fractal network model called (u,v)-flowers 20 , networks start to grow, at t = 0, from two nodes, so-called initial hubs, connected by one link.Then, subsequent (t + 1)-generations of the model are obtained from t-generations by replacing each link by two parallel paths of u > 1 and v ≥ u links long.An essential and not obvious at first glance property of this construction procedure is its equivalence to another procedure in which to obtain (t + 1)-generation one produces w = u + v copies of the previous t-generation and then joins the copies at their initial hubs.
From the second method of constriction, it is easy to see 20 that the number of links in (u,v)-flowers of generation t > 0 is given by: the number of nodes is: and the diameter of the networks grows as: Furthermore, by construction, the networks have only nodes of degree where n = 1, 2, . . .,t, and their node degree distribution is scale-free (2) with the characteristic exponent equal to: It was also shown that the box dimension (1) of (u,v)-flowers is: and their degree exponent (8) is: in accordance with the scaling relation (14).
In what follows, we show that the local scaling exponents, α and β (9), of the model are given by: and respectively, so their values satisfy the scaling relation (10).We first consider the scaling exponent β .From Eq. ( 9), it follows that if the degree of the hub inside the box increases x times, then the mass of the box will increase x β times.Correspondingly, the second method of construction of (u,v)-flowers assumes that in successive generations of these networks, the degrees of the initial hubs double, i.e. x = 2, which is due to the merger of two initial hubs from two copies of the network of the previous generation.Moreover, since the merged copies are identical, the masses of the boxes with the initial hubs also double, i.e. x β = 2. Thus, we come to the conclusion that the masses of the boxes are proportional to the degrees of their hubs, which gives β = 1, i.e.Eq. (56).
To find α, we again consider boxes with the initial hubs of degree k i = 2 t , Eq. ( 51), in networks of generation t > 0. Such boxes can be of various diameters.For example, when the diameter of the box is twice the diameter L (t−1) of the network of

Figure 2 .
Figure 2. Macroscopic characteristics of fractal complex networks: corresponding to (a) the number of boxes -N B (l B ) needed to cover the considered networks as a function of the diameter l B in the box, (b) the node degree distributions -P(k), and (c) the distributions of normalized masses of L-boxes -P(µ), for L = 3.To construct these graphs, a nested BA network of size N ≃ 5 • 10 4 and diameter d = 475 was created (this data are marked with navy circles, see also Tab. 1).To analyse the self-similarity of the network parts, the original network was covered with boxes of diameter l B = 40, and the largest box of size M ≃ 1.4 • 10 3 was extracted as a new network (this data are marked with green squares).To create the renormalized network of size N ′ ≃ 6.1 • 10 3 , the original network was covered with boxes of size l B = 6, and then each of these N B (6) = N ′ boxes was replaced with a supernode (these data are marked with red triangles).In the graph (d), results of repeated l B -renormalizations of a single box of size m ≃ 3.5 • 10 4 and diameter L = 300 are shown, which allow for an alternative determination of the box dimension of the studied fractal network (see the description in the main text of the paper).

Figure 3 .
Figure 3. Scale-invariant and self-similar scaling in real fractal networks.The graphs placed in the same column refer to the same network (i.e.WWW, brain and DBLP, respectively, starting from the left), and those placed in the same row to the same scaling relation.In particular, the following graphs show: (a-c) A log-log plot of N B versus l B revealing the fractal nature of the studied network according to Eq. (1).(d-f) Invariance of the node degree distribution P(k) under the renormalization for different box sizes l B (the case of l B = 1 corresponds to the original network).(g-i) Invariance of the normalized mass box distribution P(µ) (where L represents diameter of the considered boxes).(j-l) Scaling of the masses of boxes according to Eq. (8).(See the description given in the main text.)

Figure 4 .
Figure 4. Scale-invariant and self-similar scaling in model-based fractal networks.As in Fig. 3, the graphs placed in the same column refer to the same model of fractal networks (i.e. the nested BA networks, the SHM model, and (u,v)-flowers, respectively, starting from the left), and those placed in the same row to the same scaling relation.The presentation of data in this figure compared to Fig. 3 differs only in that the graphs relating to deterministic models show two types of points: closed and open.For these models, closed points refer to the box-covering method resulting from their deterministic construction procedure, which uses a significantly smaller number of boxes than the Song's algorithm, whose results correspond to open points (see the main text of the paper for more detailed explanation).Tab. 1 shows the results obtained based on the closed points.

Figure 5 .
Figure 5. Single step of the construction procedure of the nested BA network.First, one node is chosen with probability that is proportional to its degree, in the figure k = 4. Then the node is replaced by the corresponding BA network with the best connected node of the same degree k = 4, as the removed one.Green edges of the removed node are reconnected to randomly selected nodes of the newly created subnetwork.

Figure 6 .
Figure 6.Numerical calculation of microscopic scaling exponents.Detailed description of this figure is given in the text.

Figure 7 .
Figure 7. (2,2)-Flowers of generation t = 5.(a) The network was covered with boxes of diameter l B = 8 according to the degree of kinship of the nodes.(b) The node is covered with boxes of diameter l B = 8 according to the algorithm developed by Song et al.37 .It is clear that in the studied deterministic network, Song's random algorithm performs worse than the covering according to the kinship space.