A general model of hierarchical fractal scale-free networks

We propose a general model of unweighted and undirected networks having the scale-free property and fractal nature. Unlike the existing models of fractal scale-free networks (FSFNs), the present model can systematically and widely change the network structure. In this model, an FSFN is iteratively formed by replacing each edge in the previous generation network with a small graph called a generator. The choice of generators enables us to control the scale-free property, fractality, and other structural properties of hierarchical FSFNs. We calculate theoretically various characteristic quantities of networks, such as the exponent of the power-law degree distribution, fractal dimension, average clustering coefficient, global clustering coefficient, and joint probability describing the nearest-neighbor degree correlation. As an example of analyses of phenomena occurring on FSFNs, we also present the critical point and critical exponents of the bond-percolation transition on infinite FSFNs, which is related to the robustness of networks against edge removal. By comparing the percolation critical points of FSFNs whose structural properties are the same as each other except for the clustering nature, we clarify the effect of the clustering on the robustness of FSFNs. As demonstrated by this example, the present model makes it possible to elucidate how a specific structural property influences a phenomenon occurring on FSFNs by varying systematically the structures of FSFNs. Finally, we extend our model for deterministic FSFNs to a model of non-deterministic ones by introducing asymmetric generators and reexamine all characteristic quantities and the percolation problem for such non-deterministic FSFNs.


Introduction
Many of the complex systems around us and in various research fields of science and technology can be described by networks [1][2][3][4][5]. Since nodes and edges, the constituents of networks, represent a wide variety of objects and interactions, respectively, Euclidean distances are not always defined for networks. The absence of the Euclidean distance eliminates the limitation of the number of edges connecting to a node, namely the degree k of a node, and thus allows a large fluctuation of k. In fact, degree distributions P(k) of many real-world networks obey power-law functions for large k, i.e. P(k) / k −γ with an exponent γ [6]. We define the shortest path distance to be the minimum number of edges connecting two nodes even for a network where the Euclidean distance is not defined. We can quantify the linear distance over a network by the average shortest-path distance hli or the network diameter L defined as the largest shortest-path distance. If the diameter L of a network G scales with the number of nodes N as L / log N, G is referred to as a small-world network [7]. In contrast, if the relation L / N 1=D f holds, G is said to be a fractal network with the fractal dimension D f [8]. Due to the smallworld nature of many real-world networks, a lot of structural models of small-world and scalefree networks have been proposed [6,[9][10][11][12][13][14][15][16], and various phenomena or dynamics on them have been extensively studied [17][18][19][20][21]. Actual scale-free networks, however, often possess fractal structures in shorter length scales than their network diameters or the average shortestpath distances. World Wide Web, protein interaction networks, and actor networks are known to be examples of real-world fractal scale-free networks [8,22,23]. Nevertheless, there is less research on fractal scale-free networks (FSFNs) than on small-world scale-free networks.
In particular, the lack of a structural model of FSFNs that can freely control the exponent γ, fractal dimension D f , and other structural features delays the study of phenomena or dynamics on FSFNs.
There are two representative models of FSFNs, the (u, v)-flower model [24] and the Song-Havlin-Makse (SHM) model [25]. The (u, v)-flower model constructs hierarchically a highly cycle-rich FSFN and can vary the network structure to some extent by adjusting the parameters u and v (u � v � 2). The clustering coefficient is, however, always zero independently of u and v. Meanwhile, the SHM model forms a tree-like FSFN and can change the fractal dimension D f (and then γ followed by the change in D f ) by controlling the parameter z characterizing the growth rate of nodes. Since the network has a tree structure, the clustering coefficient is also zero for any value of the parameter z in the SHM model. In addition to these two structural models, several deterministic models for FSFNs have been proposed so far [26][27][28][29][30][31][32]. These are, however, classified into derivatives of the (u, v)-flower model or the SHM model, or synthetic models in which the network structure cannot be freely controlled.
In this paper, we propose a general model of hierarchical FSFNs which can change widely and freely the structural features of obtained networks. This model constructs an FSFN hierarchically by replacing each edge in the previous generation network iteratively with a small graph called a generator. By choosing the generator appropriately, it becomes possible to control the scale-free nature, fractality, clustering property, and the nearest-neighbor degree correlation. We actually calculate analytically the exponent γ describing the scale-free property, fractal dimension D f , average clustering coefficient C, global clustering coefficient C 4 , and joint probability P(k, k 0 ) defining the nearest-neighbor degree correlation for networks formed by the present model. Furthermore, the bond-percolation problem on constructed FSFNs is investigated as an example of analyses of phenomena occurring on them. Specifically, we analytically derive the percolation critical point and critical exponents for infinitely large FSFNs, by using structural information of the generator. In addition, a model of non-deterministic FSFNs is provided by introducing asymmetric generators, and various statistical properties of resulting networks are examined theoretically.

Model
We first prepare a small connected graph G (called generator hereafter) in which two particular nodes are specified as root nodes. As shown in Fig 1A, a network in the t-th generation, G t , is constructed by replacing every edge in G tÀ 1 with the generator G so that the terminal nodes of the edge coincide with the root nodes of G. This procedure is an inversion operation of the renormalization transformation that replaces small subgraphs with superedges [24,33]. Although the initial network G 0 can be arbitrarily chosen, we fix, in this work, G 0 to be a single edge connecting two nodes for simplicity. In order to obtain a deterministic fractal scale-free network (FSFN), the generator G must satisfy the following three conditions: (1) The degree of the root node is no less than 2.
(2) The shortest-path distance between the two root nodes is 2 or longer.
(3) The two root nodes are symmetric to each other in G.
The first and second conditions guarantee the scale-free property and the fractal nature of the constructed network, respectively. If the first condition is violated and the degree of the root node is one, our model produces a non-scale free network with an exponentially damped degree distribution. Such networks can also be formed by existing models [34,35]. In the case of the violation of the second condition, namely, when the two root nodes are directly connected by an edge, our model gives hierarchical scale-free networks with the small-world property such as a network modeled by [36]. It should be emphasized that even if the conditions (1) and (2) are violated, all the analytical arguments below still hold, except for those in the sections "Scale-free property", "Fractal property", and "Degree correlation". The third condition is for the model to be deterministic. The meaning of "symmetric" in the third condition is the following. If a network constructed by removing one root node and its edges from G has the same topology as a network formed by removing another root node and its edges from G, these root nodes are called symmetric.
The most distinct advantage of this model is that a generator G can be chosen arbitrarily (within the three conditions) and thus a wide variety of hierarchical FSFNs are produced by this model, whereas previous models for FSFNs based on inverse renormalization procedure either fix G or limit the structural freedom of G to a narrow range. Therefore, the present model is able to construct FSFNs proposed by previous models [24][25][26][27][28][29][30]. For example, the (u, v)-flower [24] is reproduced by choosing a generator G of a cycle with u + v nodes as shown in Fig 1B. If we start with a star graph with 4 leaves as the initial graph G 0 instead of an edge and employ a generator consisting of two connected star graphs with z leaves as shown in Fig 1C, we obtain the Song-Havlin-Makse (SHM) model [25]. Furthermore, by using a generator with the different number of edges connecting the root nodes in Fig 1C, we can generate the fractal scale-free tree proposed by [29,30].

Properties of fractal scale-free networks
Any characteristics concerning a constructed network G t in the t-th generation are completely determined by the nature of the generator G. Many of quantities defining G t can be analytically calculated because of the simplicity of the model. We present here several indices of G t by using quantities describing the generator G.
Numbers of edges and nodes. The most fundamental quantities are the numbers of nodes and edges in G t . The number of edges M t in the t-th generation network G t is m gen times larger than the number of edges in G tÀ 1 , namely, where m gen is the number of edges in the generator G. This relation with M 0 = 1 immediately leads to The number of nodes N t in G t is the sum of the number of nodes N t−1 in the previous generation network G tÀ 1 and the number of newly added nodes in the replacement of edges in G tÀ 1 . It is convenient for counting newly added nodes and for later discussion to define remaining nodes which are the nodes in the generator G other than the root nodes. For each edge in G tÀ 1 , the edge replacement with G introduces n rem (= n gen − 2) nodes, where n gen and n rem are the numbers of the entire nodes and remaining nodes in G, respectively. The number of newly added nodes in G t is thus given by M t−1 n rem and N t is expressed by N t = N t−1 + M t−1 n rem . Solving this recurrence relation with N 0 = 2 with the aid of Eq (2), the number of nodes in G t is given by For t � 1, this equation leads to the approximate relation The quantities M t and N t are not influenced by the structure of the generator G but are determined only by the numbers of nodes and edges in G. Many other indices characterizing G t , however, depend on the topology of G as shown below. Let us consider the number of nodes N t (k) of degree k in G t . A node of degree k in G tÀ 1 has the degree κk in G t , where κ is the degree of the root node in G. Therefore, the number of nodes of degree k in G t is the sum of the number of nodes of degree k/κ in G tÀ 1 only if k/κ is an integer and the number of degree k nodes which are newly added in the edge replacement operation. Thus, for t � 1, we have the relation, where N t (k/κ) = 0 if k/κ is non-integer, the summation is taken over the n rem remaining nodes in G, k rem n is the degree of the n-th remaining node, and δ k,k 0 is the Kronecker delta.

PLOS ONE
Considering that N 0 (k) = 2δ k,1 for G 0 of a single edge, Eq (5) gives We can calculate N t (k) if κ, m gen , n rem , and k rem n characterizing the generator G are given. Using the above expression of N t (k), the average degree hki t and the average squared degree hk 2 i t of G t are evaluated. Since it is obvious that the average degree must be given by hki t = 2M t /N t , Eqs (2) and (3) lead to hki t ¼ 2m t gen ðm gen À 1Þ 2ðm gen À 1Þ þ n rem ðm t gen À 1Þ : This can be confirmed by calculating ∑ k kN t (k)/N t with the aid of Eq (6). In the limit of infinite size N t (t ! 1), the average degree is then given by The average squared degree hk 2 i t = ∑ k k 2 N t (k)/N t is obtained by using Eq (6) as Taking the summation over t 0 in Eq (9) and using Eq (3) for N t , hk 2 i t is expressed as hk 2 i t ¼ ðm gen À 1Þ½2k 2t ðm gen À k 2 Þ þ K rem 2 ðm t gen À k 2t Þ� ðm gen À k 2 Þ½2ðm gen À 1Þ þ n rem ðm t gen À 1Þ� for m gen 6 ¼ k 2 k 2ðtÀ 1Þ ðk 2 À 1Þð2k 2 þ tK rem 2 Þ 2ðk 2 À 1Þ þ n rem ðk 2t À 1Þ for m gen ¼ k 2 For t ! 1, this quantity becomes As will be shown in the next section, the condition for the convergence of hk 2 i 1 is equivalent to γ > 3, where γ is the exponent of the degree distribution P(k) / k −γ for k � 1. Scale-free property. Let us consider the asymptotic behavior of the degree distribution P (k) of G t for k � 1 and t � 1. The degree of a node in the (t − 1)-th generation network G tÀ 1 is multiplied by κ in G t . Therefore, the number of nodes with degree k in G tÀ 1 is identical to the number of nodes with degree κk in G t if k is larger than the maximum degree of remaining nodes in G, that is, This can be directly confirmed by Eq (5). Since it is natural to suppose that the degree distribution for a large generation converges to a specific functional form P(k), the above equation is written, in a continuum approximation for k, as for large t and k. Using Eq (4), this relation leads to and we have the solution of this functional equation as where This implies that a constructed network possesses the scale-free property and the exponent of the power-law degree distribution is determined by the degree κ of the root node and the number of edges m gen in the generator. If the condition (1) for the generator G is violated, the exponent γ diverges and the network G t does not exhibit the scale-free property. The exponent γ always satisfies because m gen � 2κ for any generator. It is easy to check that the exponent γ for the (u, v)flower [24], the SHM model [25], and their derivatives [26][27][28][29][30] can be reproduced by Eq (17). It should be emphasized that Eq (16) does not mean N t (k) / k −γ for large t and k. This is because degrees of nodes in G t can actually take exponentially discretized values, such as k, κk, κ 2 k, . . ., whereas P(k) is defined in the continuum approximation for k. The number of nodes of degree k is then given by N t ðkÞ ¼ N t R kk k Pðk 0 Þdk 0 , which provides the asymptotic behavior of N t (k) as where for a large generation. Finally, we consider the condition for the convergence (or divergence) of hk 2 i 1 as shown in Eq (12). If m gen > κ 2 , log m gen /log κ is larger than 2, namely, γ given by Eq (17) is large than 3. Therefore, hk 2 i 1 becomes finite if γ > 3, while it diverges for γ � 3. This is quite reasonable because of the form of Eq (16).
Fractal property. Fractality of a network can be examined by the relation between the number of nodes in the network and the network diameter, the maximum value of the shortest-path distance, according to the definition of fractal networks mentioned in the Introduction section. Let L t−1 be the diameter of G tÀ 1 . Then, each of L t−1 edges between two nodes separated by the diameter is replaced with the generator G in the network G t and the shortestpath distance between these two nodes becomes λL t−1 , where λ is the shortest-path distance between the two root nodes in G. If the generation t is large enough, the distance λL t−1 is almost the longest shortest-path distance in G t . More precisely, the diameter of G t is given by λL t−1 + L 0 , where L 0 is a constant. If λ gives the diameter of G, we have L 0 = 0. Even in the case that the diameter of G is longer than λ and L 0 is finite, however, λL t−1 is much larger than L 0 for t � 1 and we can ignore the term of L 0 . Therefore, the diameter L t of G t is expressed as On the other hand, the numbers of nodes N t and N t−1 are related by Eq (4) for t � 1. Eqs (4) and (21) implies that if the network diameter is λ times larger, the number of nodes becomes m gen times larger. This leads to the relation between N t and L t as where This result shows that the network G t formed by our model exhibits the fractal nature with the fractal dimension D f given by Eq (23). The fractal dimensions of the (u, v)-flower [24] and the SHM model [25] are reproduced by Eq (23). If the condition (2) for the generator G is violated, namely if the two root nodes are directly connected in G, the fractal dimension D f diverges and G t becomes a small-world network.
Clustering property. The clustering property of a network is often characterized by two types of quantities, namely the average clustering coefficient [7] given by 2D n k n ðk n À 1Þ ; ð24Þ and the global clustering coefficient [37,38] defined as where N is the network size, k n is the degree of node n, and Δ n is the number of triangles including the node n. The quantity C is the average of the local clustering coefficient C(n) = 2Δ n /[k n (k n − 1)], and C 4 is the ratio of three times the total number of triangles in the network (∑ n Δ n ) to the total number of connected triplets of nodes [∑ n k n (k n − 1)/2]. We can calculate analytically these two clustering coefficients for FSFNs formed by our model. At first, we derive the average clustering coefficient C t of G t in the generation t. It should be noted that all triangles in G t are produced in the edge replacement procedure to form where C rem is the average clustering coefficients of the remaining nodes in the generator G, i.e., The N t−1 nodes left in G t are those inherited from G tÀ 1 . The degree of a node whose degree in G tÀ 1 is k becomes κk in G t , and this node contributes to kΔ R triangles in G t , where Δ R is the number of triangles in the generator which include one of the root nodes. Therefore, the average clustering coefficient C t of the network G t is written as Substituting Eq (6) for N t−1 (k), we can express C t by quantities characterizing the generator. In the limit of t ! 1, the average clustering coefficient becomes This result implies that if there exists even one triangle in the generator, the average clustering coefficient of G t will be finite, no matter if the size of G t is finite or infinite. Next, we consider the global clustering coefficient C 4 defined by Eq (25). The denominator ∑ n k n (k n − 1) is obviously expressed as N t (hk 2 i t − hki t ) for G t . The quantity ∑ n Δ n in the numerator is equal to three times the total number of triangles in the network. Since all triangles in G t are produced in the procedure of replacing M t−1 edges in G tÀ 1 with the generator G, the total number of triangles in G t is given by M t−1 Δ gen , where Δ gen is the number of triangles included in G. Thus, the global clustering coefficient C D t for G t is given by where hki t and hk 2 i t are presented by Eqs (7) and (11), respectively. In contrast to C 1 , even if for m gen > κ 2 is finite and is, by means of Eqs (8) and (12), expressed as where hki gen and hk 2 i gen are the average degree and average squared degree of the generator G, respectively. Degree correlation. The degree correlation between nodes adjacent to each other is described by the joint probability P(k, k 0 ) that one terminal node of a randomly chosen edge has degree k and the other terminal node has degree k 0 . This probability is expressed by using the number M(k, k 0 ) of (k, k 0 )-edges, i.e. edges connecting nodes with degrees k and k 0 , as where M is the total number of edges in the network. We can derive the joint probability P t (k, k 0 ) for G t by counting the number M t (k, k 0 ) of (k, k 0 )-edges in G t . Since all edges in G t are yielded in the replacement of M t−1 edges in G tÀ 1 with the generator G, M t−1 m rem (k, k 0 ) edges in G t contribute to M t (k, k 0 ), where m rem (k, k 0 ) is the number of (k, k 0 )-edges connecting remaining nodes in G. In addition, M t (k, k 0 ) includes the contribution from edges between nodes inherited from G tÀ 1 and their neighboring nodes. Considering these contributions, where μ(k) is the number of nodes with degree k adjacent to one of the root nodes in G. The quantity N t−1 (k/κ) represents the number of nodes inherited from G tÀ 1 whose degrees become k in G t . Note that N t−1 (k/κ) is zero for non-integer k/κ and is related to N t (k) by Eq (5). Substituting Eq (32) into Eq (31), we obtain The first term of Eq (33) gives the contribution from edges between remaining nodes in G accompanying the edge replacement and the second term comes from edges between inherited nodes from the previous generation and their neighbors. All quantities on the right-hand side of Eq (33) are determined from the structure of G. It is thus possible to evaluate analytically the nearest-neighbor degree correlation in the FSFN formed by a given generator. If a network has no degree correlation, the joint probability must be given by a product of a function of k and a function of k 0 . Since P t (k, k 0 ) presented by Eq (33) cannot be expressed in the form of such a product, there exists some sort of nearest-neighbor degree correlation in G t . The type of the degree correlation is revealed by various measures, such as the assortativity [39], Spearman's degree rank correlation coefficient [40], and average degree of the nearest-neighbors of nodes with degree k [41]. These measures are calculated by the joint probability P t (k, k 0 ).

Percolation problem
A wide choice of generators allows us to prepare a variety of FSFNs. Many indices characterizing structural features of constructed FSFNs are analytically evaluated as shown in the previous section. It becomes possible, by employing these FSFNs, to investigate systematically how various phenomena and dynamics occurring on FSFNs are influenced by the characteristics of networks. The percolation transition is one of the most fundamental phenomena on complex networks and is deeply related to the robustness of networks to failure of nodes or edges [42][43][44][45] or the spread of disease [21,[46][47][48]. Although percolation processes in small-world networks have been extensively studied so far, as reviewed by [19], our understanding of the percolation problem for FSFNs is still limited [49][50][51][52][53]. For example, the relation between the percolation threshold and a specific structural feature of FSFN, such as clustering property, has not yet been systematically studied. In this section, we discuss the percolation transition in FSFNs formed by our model. We concentrate here on the bond-percolation process corresponding to edge failure, because the site-percolation problem reflecting node failure is much more complicated. We calculate analytically the percolation transition point and some critical exponents for the bond-percolation problem in which edges are randomly removed from an infinite FSFN G 1 with probability 1 − p. For this purpose, it is convenient to define the renormalized root nodes (RRNs) of G t as two nodes corresponding to the root nodes of the renormalized network (in the sense of the edge renormalization) of G t by G tÀ 1 , which takes the same topology as the generator G. Namely, RRNs of G t are the oldest two nodes in G t . An example of RRNs is illustrated by red nodes in Fig 2B or 2C, which will be referred to in the section "Order parameter exponent".
Critical point. For an infinitely large FSFN (t ! 1), the shortest-path distance between the RRNs diverges. Thus, if the RRNs of G t are still connected to each other in a networkG t ðpÞ which is formed by removing edges randomly with probability 1 − p from G t , the network G t ðpÞ is considered to be percolated. Since the network G t is composed of m gen pieces of G tÀ 1 , the probability R t (p) of the RRNs of G t being connected to each other inG t ðpÞ is related to the probability R t−1 (p) of the RRNs of G tÀ 1 being connected inG tÀ 1 ðpÞ as where π(p) is the probability that the two root nodes of the generator G are connected to each other in a network where edges are randomly removed from G with probability 1 − p. Eq (34) has an unstable fixed point at p = p c , i.e. R t (p c ) = R t−1 (p c ) = p c . The probability p c gives the percolation critical point, because the percolation probability R 1 (p c ) is finite. Therefore, the critical point p c is presented by the non-trivial and meaningful solution of the equation, We can determine the functional form of π(p) for a given generator G. The probability π(p) is expressed as where s m is the number of subgraphs of G with m edges in which the root nodes are connected. The above summation starts from m = λ, because s m = 0 for m < λ. It is easy to count the number s m by finding numerically such subgraphs from all subgraphs of G, because the total number of subgraphs of G is only 2 m gen with m gen not very large usually. The function π(p) is in general an m gen -th degree polynomial in p, but the polynomial degree can be reduced if there exist edges in G that do not contribute to paths connecting the root nodes. In such a case, the degree of π(p) becomes equal to the number of edges in G˘, where G˘is the core subgraph of G consisting only of all edges that contribute to paths between the root nodes. The coefficient Correlation length exponent. The correlation length (in the sense of the shortest-path distance) of a percolation networkG 1 ðpÞ near the critical point p c behaves as where ν is the critical exponent for the correlation length and ξ 0 is a constant. When we renormalize the substrate network G 1 by the generator G, the edge occupation probability in the renormalized network G 0 1 is given by π(p). Thus, the correlation length of the renormalized percolation networkG 0 1 ½pðpÞ� is the same as ξ, as expressed by The coefficient x 0 0 is λ times larger than ξ 0 , i.e., x 0 0 ¼ lx 0 , because the root nodes in G are separated by λ. Eqs (37) and (38) then lead to Therefore, taking the limit p ! p c , the correlation length exponent is given by where π 0 (p) is the first derivative of π(p). As well as the critical probability p c , the exponent ν also depends only on the structure of the core subgraph G˘of G, because λ is determined by G˘.
The correlation volume N ξ is the average number of nodes within the radius ξ from a node, which is presented by N x ¼ N x0 jp À p c j Àñ , whereñ is the correlation volume exponent and N ξ0 is a constant. Similarly to the argument of ξ, the correlation volume of the renormalized percolation networkG 0 1 ½pðpÞ� is the same as N ξ ofG 1 ðpÞ, namely N x ¼ N 0 x0 jpðpÞ À p c j Àñ . The coefficient N 0 x0 is given by N 0 x0 ¼ m gen N x0 because of Eq (4). These relations give us the exponentñ as From Eqs (40), (41), and (23), we have the simple relatioñ which is derived also from N x / x D f . Order parameter exponent. Let us consider the critical exponent β for the order parameter P 1 , i.e., the probability of a randomly chosen node belonging to the giant connected component inG 1 ðpÞ. We can calculate β by extending the argument for the (u, v)-flower [49] to our general model. The order parameter P 1 is the limiting value of P t (p, N t ) for N t ! 1, where P t (p, N t ) is the probability that a node in the t-th generation FSFN G t with N t nodes belongs to the largest component inG t ðpÞ. According to the finite-size scaling theory [54], the quantity P t (p, N t ) for a large generation t and near the critical point p c must have the form of

PLOS ONE
where F(x) is a scaling function. At the critical point p = p c , we then have the relation here we used Eq (4) for the last equation. Therefore, Eq (41) leads to the expression for the exponent β as where ω c is the ratio of P t (p c , N t ) to P t−1 (p c , N t−1 ), namely, Since the probability π(p) is presented by Eq (36) for a given generator G, we can calculate analytically the exponent β from Eq (45) if ω c is obtained for G.
In order to calculate the ratio ω c , we introduce two probabilities S t and T t . The quantity S t is the probability that a randomly chosen node Q is connected to one of the RRNs of G t in the percolation networkG t ðpÞ, and T t is the probability that the chosen node Q is connected to both RRNs. The probability P t (p, N t ) for a large t is then given by because the probability of the node Q being at a finite distance from either of the RRNs is almost zero for t ! 1. As illustrated in Fig 2, the probabilities S t and T t can be expressed as functions of S t−1 , T t−1 , and the probability R t−1 (p) of the two RRNs of G tÀ 1 being connected to each other. Since these functions are linear with respect to S t−1 and T t−1 , we can express (S t , T t ) as where W is a two-by-two matrix whose matrix element w ij with i, j = 1 or 2 is a function of R t−1 (p). For a large enough generation t, the largest eigenvalue ω of the matrix W gives the ratio (S t + T t )/(S t−1 + T t−1 ) and thus P t (p, N t )/P t−1 (p, N t−1 ) from Eq (47), because the vector (S t , T t ) T becomes proportional to the eigenvector belonging to ω. Therefore, the ratio ω c defined by Eq (46) is simply the largest eigenvalue of W at p = p c and for t ! 1. Let us determine the functional forms of the matrix elements w ij [R t−1 (p)] from the structure of the generator G. As seen from the form of Eq (48), the element w 11 is, for example, the conditional probability that a randomly chosen node Q is connected to one of the two RRNs of G t under the condition that the node Q is connected to one of the RRNs of the subgraph G tÀ 1 including Q. Various connection patterns contribute to this conditional probability w 11 . A situation illustrated in Fig 2B also contributes to w 11 . Since the probability of this situation occurring is R 2 tÀ 1 ½1 À R tÀ 1 � 2 S tÀ 1 , this connection pattern is incorporated into w 11 as a term of R 2 tÀ 1 ½1 À R tÀ 1 � 2 . By the same token, w 12 contains a contribution from a situation shown in Fig 2C. As demonstrated by these examples, the matrix element w ij is expressed as a polynomial consisting of terms of R m tÀ 1 ½1 À R tÀ 1 � m gen À mÀ 1 with 0 � m � m gen − 1. Since R t (p c ) becomes equal to p c for t ! 1, the matrix element w ij in the thermodynamic limit at p = p c is presented by where c ij (m) is the number of connection patterns, with m connected subgraphsG tÀ 1 ðp c Þ, that a randomly chosen node Q is connected to RRNs of G t . The prefactor 1/m gen for j = 2 is the probability that the node Q is included in a specific subgraph G tÀ 1 , and 1/2m gen for j = 1 is the probability that the node Q is in a subgraph G tÀ 1 and one of the RRNs of G tÀ 1 is chosen as the connection point of Q. In order to calculate the coefficient c ij (m) from the structure of G, we consider subgraphs G 0 (e 0 ) formed by removing single edges e 0 from G. The edge e 0 corresponds to the subgraphG tÀ 1 ðp c Þ including the node Q. The coefficient c ij (m) is the number of subgraphs of G 0 (e 0 ) for any e 0 with m edges, in which j(= 1 or 2) terminal nodes of e 0 are connected to i(= 1 or 2) root nodes of G under the condition that only for j = 2 the terminal nodes of e 0 are considered to be connected directly to each other even though the edge e 0 is absent in G 0 (e 0 ). Because of the small size of G, the numbers of these subgraphs can be counted numerically, as in the case of the evaluation of s m in Eq (36). We can eventually obtain the order parameter exponent β from Eq (45) by calculating the largest eigenvalue ω c of the matrix W whose elements are given by Eq (49).

Examples.
As examples of the above general argument, let us demonstrate the structural features of FSFNs formed by two kinds of generators, and compare the critical properties of percolation on these networks. We employ two generators G A and G B shown in Fig 3A and 3B, respectively. These generators have the same number of nodes (n gen = 6), number of edges (m gen = 8), degree of the root node (κ = 2), shortest-path distance between the root nodes (λ = 3), and degrees of remaining nodes (k rem n ¼ 3 for any n). The similarity found in the generators G A and G B leads to similar structural features of the FSFNs constructed by them. In fact, the FSFNs G A 1 and G B 1 formed by G A and G B , respectively, in the infinite generation possess the same average degree hki 1 = 7/2 calculated by Eq (8), same second moment hk 2 i 1 = 63/4 by Eq (12), same scale-free property γ = 4 by Eq (17), and same fractal dimension D f = log 8/ log 3 = 1.893 by Eq (23). In addition, the joint probabilities P t (k, k 0 ) given by Eq (33), which describe the nearest-neighbor degree correlation, are also the same for G A t and G B t for any t, because the nearest-neighbor degree correlations in G A and G B are the same. As a result, the Spearman's degree rank correlation coefficient becomes % = −21/64 for both FSFNs in the infinite generation.
Clustering properties of these two FSFNs are, however, different from each other. Since there is no triangle in G B , the average clustering coefficient C t and the global clustering coefficient C 4 t are zero for any generation FSFN formed by G B . On the contrary, the generator G A has two triangles and the local clustering coefficients of all nodes in G A are finite. Therefore, both of the two kinds of clustering coefficients C t and C 4 t take finite values. The average clustering coefficient of G A 1 is C 1 = 0.31486 which is obtained from Eq (28) with Δ R = 1, C rem = 4/3, and k rem n ¼ 3 for any n. We can also calculate the global clustering coefficient of G A 1 as C 4 1 ¼ 3=14 from Eq (30) (49) for the matrix elements w ij are also obtained by counting numerically the numbers of subgraphs of the generators satisfying required conditions. The largest eigenvalues of the matrix W are then computed as ω c = 0.9649 for G A and 0.9344 for G B . These quantities characterizing the generators, such as π(p), ω c , m gen , and λ, give the percolation critical point p c and critical exponents for bond percolation on G A 1 and G B 1 as shown in Table 1. The validity of these values, as well as structural measures, has been confirmed by numerical calculations. Although most of the structural features of G A 1 and G B 1 are the same, except for the clustering property, G B 1 is more robust than G A 1 against edge elimination and the percolation transitions on these FSFNs belong to different universality classes with different critical exponents. Universality classes of the percolation transition in complex networks are reviewed in [18]. As can be seen from these examples, our generalized model enables us to examine systematically the relationship between certain structural properties of FSFNs, such as the clustering property, and phenomena or dynamics on them.

Asymmetric generator
Up to here, a generator must satisfy the conditions (1), (2), and (3) as mentioned in the Model section. The condition (3) guarantees that the resulting FSFN has a deterministic structure. We consider in this section the case that the condition (3) is violated, namely, the generator is Table 1. Values of the percolation critical point p c and critical exponents ν,ñ, and β for FSFNs formed by the generators G A and G

PLOS ONE
asymmetric. An asymmetric generator G is a network whose subgraph obtained by removing one root node and its edges from G has a different topology from a network obtained by removing another root node and its edges from G. Since the two root nodes of an asymmetric generator are not equivalent, there are two ways to replace an edge with the generator as illustrated by Fig 4. Here, we assume that the way of edge replacement is randomly chosen with the probability 1/2. This stochasticity makes a final network G t non-deterministic, but G t is still fractal and has the scale-free property in a statistical sense, as will become clear in the discussion below. We show how the various results obtained for symmetric generators are modified by the asymmetry of generators. The numbers of nodes N t and edges M t are given by Eqs (2) and (3), respectively, as in the case of symmetric generators, because these quantities are not influenced by the symmetry of generators. The number of nodes of degree k is, however, different from Eq (6). The probability that k 1 edges from a node of degree k are replaced with the generator G in one way and the remaining (k − k 1 ) edges are replaced with G in another way is given by k ing the probability that k 1 edges from a node of degree k 0 in the (t − 1)-th generation network G tÀ 1 are replaced with G in one specific way, the number of nodes N t (k) of degree k in G t is presented by where κ 1 and κ 2 are the degrees of the two root nodes of G. If κ 1 = κ 2 , the above relation becomes equivalent to Eq (5). We can obtain N t (k) by solving the recurrence Eq (50) with the initial condition N 1 ðkÞ ¼ d kk 1 þ d kk 2 þ P n rem n¼1 d k;k rem n . The moments hki t and hk 2 i t can be calculated from this N t (k). The average degree hki t is given by Eq (7) [or (8) for t ! 1], because N t and M t are unchanged from those for symmetric generators. Computing hk 2 i t = ∑ k k 2 N t (k)/N t , the second moment in the infinite generation limit is expressed as
In order to examine the scale-free property of G t , let us consider the asymptotic behavior of N t (k) for large values of t and k. Neglecting stochastic fluctuations in edge replacements taken into account in Eq (50), one can find the asymptotic form of N t (k). In an average sense, k/2 edges from a node of degree k in G tÀ 1 are multiplied by κ 1 and the remaining k/2 edges are multiplied by κ 2 in G t . Therefore, we have This expression is the same as Eq (13) with � k instead of κ. Thus, according to the argument deriving Eq (17) from Eq (13) in the section "Scale-free property", the degree distribution P(k) obeys asymptotically P(k) / k −γ , where the exponent γ is given by This implies that a network G t formed by an asymmetric generator also exhibits the scale-free property. It is interesting to note that G t is scale-free as long as � k is greater than 1, even if the degree of one of the root nodes is unity, regardless of the condition (1) for generators. Here, we should emphasize that the relation between N t (k) and P(k) depends on the values of κ 1 and κ 2 . If κ 1 = κ 2 holds, degrees in G t take exponentially discretized values, such as k, � kk, � k 2 k, . . ., as in the case of symmetric generators, and N t (k) is proportional to kP(k), i.e. N t (k) / k −γ 0 with g 0 ¼ log m gen =log � k. On the other hand, degrees in G t take nearly uniform values if κ 1 6 ¼ κ 2 . In this case, the exponent γ 0 becomes equal to γ, because N t (k) / P(k). Therefore, the exponent γ 0 describing the asymptotic behavior of N t (k) is presented by The fractal property of a network G t is not influenced by the symmetry of the generator. This is because the relation between the diameter L t−1 of the (t − 1)-th generation network G tÀ 1 and L t for G t is still expressed by Eq (21) even for the asymmetric generator, and the number of nodes N t also remains as given by Eq (3). Therefore, G t formed by an asymmetric generator keeps the fractal property with the fractal dimension given by Eq (23).
In order to obtain the expression of the average clustering coefficient of the FSFN formed by an asymmetric generator, we need to distinguish the numbers of triangles Δ R1 and Δ R2 including the root nodes R 1 and R 2 , respectively. The probability that k 1 Δ R1 triangles arise from the replacement of k 1 edges from a node of degree k with the generator G and (k − k 1 )Δ R2 triangles arise from the replacement of the remaining (k − k 1 ) edges is k Considering this probability, the average clustering coefficient is presented by where h(k, k 1 ) = κ 1 k 1 + κ 2 (k − k 1 ) and N t−1 (k) is given by Eq (50). This corresponds to Eq (27) for symmetric generators. Since the total numbers of triangles and connected triplets in the network G t are not affected by the symmetry of the generator, the global clustering coefficient C D t is presented by Eq (29). In the thermodynamic limit (t ! 1), C D 1 ¼ 0 for m gen � � k 2 and C D 1 is given by Eq (30) with � k instead of κ if m gen > � k 2 . In the symmetric generator case, only nodes with degree k/κ in G tÀ 1 can become nodes with degree k in G t . The joint probability P t (k, k 0 ) is then provided by Eq (33). For FSFNs by asymmetric generators, however, all nodes with degree k 00 in G tÀ 1 can be nodes with degree k in G t if k 00 satisfies h(k 00 , k 1 ) = k for arbitrary integer k 1 2 [0, k 00 ]. Taking into account the probability of choosing k 1 edges from k 00 edges, we can write the joint probability P t (k, k 0 ) as where J k 00 k 1 ðk; k 0 Þ ¼ ½k 1 m 1 ðkÞ þ ðk 00 À k 1 Þm 2 ðkÞ�d hðk 00 ;k 1 Þ;k 0 and μ 1(2) (k) is the number of nodes with degree k adjacent to the root node R 1 (2) in G.
The percolation process on our FSFN depends only on how the two root nodes are connected to each other in the generator G, irrespective of the symmetry of G. Thus, the critical point and critical exponents of the bond-percolation transition on FSFNs formed by asymmetric generators are also obtained by the same argument as that for symmetric generators. Although the renormalization of G t by G tÀ 1 is not uniquely determined in the asymmetric generator case because of the structural fluctuation in these networks, G t can be renormalized into the network with the same topology as the generator G by using any of the various realizations of G tÀ 1 . The RRNs of G t or its subgraph G tÀ 1 by which G t is renormalized are also defined as the nodes corresponding to the root nodes of the renormalized network. Eventually, the critical point p c is given by the solution of Eq (35), and the critical exponents ν,ñ, and β are provided by Eqs (40), (41), and (45), respectively. The only difference from the case of the symmetric generator is that the numbers of subgraphs s m and c ij (m) are calculated for the asymmetric generator G.
As an example, let us consider the FSFN G C t formed by the asymmetric generator G C illustrated by Fig 5A. Basic measures characterizing G C are m gen = 7, n rem = 3, κ 1 = 3, and κ 2 = 2. Since κ 1 6 ¼ κ 2 , the degree distribution P(k) of G C t is simply obtained by N t (k)/N t and is analytically calculated by using Eqs (3) and (50). The solid line in Fig 5B represents the theoretical result of P(k) for the 6th generation FSFN G C 6 . Symbols showing the numerical result agree well with this theoretical curve. The peaks equally spaced on the logarithmic k-axis reflect the nested appearance of the binomial distribution in N t (k) given by Eq (50). Eq (55) gives the power-law exponent describing P(k) as γ = 3.1237. The envelope of P(k) surely exhibits this power-law behavior as shown in Fig 5B. The average degree and average squared degree of G C 1 are hki 1 = 4 and hk 2 i 1 = 220/3, respectively. The fractal dimension of G C 1 is D f = log 7/log 2 � 2.8074. Since the generator contains triangles and m gen > � k 2 , the network G C 1 is clustered in both senses of C 1 and C D 1 , as is clear from the fact that C 1 = 0.3952 and C D 1 ¼ 9=182 � 0:04945. The degree correlation of G C t can be evaluated by the joint probability P t (k, k 0 ) which is computed by Eq (58). We can calculate the assortativity r t [39] and Spearman's degree rank correlation coefficient % t [40] from P t (k, k 0 ). Although r t is zero for t ! 1 because γ < 4 and then hk 3 i 1 = 1, the Spearman's correlation coefficient % t is finite for any t. In the thermodynamic limit, this coefficient is calculated as % 1 = −0.5221, which indicates that G C 1 exhibits disassortative degree correlations between neighboring nodes. Percolation properties of G C 1 can be also clarified by adapting the argument in the section "Percolation problem" for symmetric generators. The critical point and critical exponents for the bond-percolation process on G C 1 are calculated as p c = 0.4473, ν = 1.2561,ñ ¼ 3:5262, and β = 0.2761 from Eqs (35), (40), (41), and (45), respectively. The validity of these results is demonstrated in Fig 5C that is a scaling plot depicting P t ðp; N t ÞN b=ñ t as a function of N 1=ñ t jp À p c j for G C t in various generations t. According to Eq (43), the fact that these plots for different N t fall on the same curve, as shown in Fig 5C, implies that the calculated p c ,ñ, and β are correct values. It should be emphasized that data points for each size in Fig 5C are obtained from a single realization of G C t and P t (p, N t ) is almost independent of samples for large t. This is because the global connectivity of G t is irrelevant to the way of edge replacements.

Conclusion and discussion
In this work, we have proposed a general structural model of fractal scale-free networks (FSFNs) and calculated analytically various measures characterizing structures of constructed https://doi.org/10.1371/journal.pone.0264589.g005

PLOS ONE
networks. As an example of analyses of phenomena occurring on our FSFNs, the percolation problem on infinite FSFNs has been studied. Using the present model, one can provide a wide variety of deterministic and non-deterministic FSFNs which include those formed by existing models and examine systematically the influence of a specific structural property on a phenomenon on FSFNs.
To construct an FSFN, we first prepare a small graph called a generator G in which two particular nodes are specified as root nodes. The degrees of the root nodes must be no less than 2, and the shortest-path distance between the two root nodes have to be 2 or longer. An FSFN in the t-th generation, G t , is formed by replacing every edge in the previous generation network G tÀ 1 with the generator G iteratively so that the terminal nodes of the edge coincide with the root nodes of G. If the generator G is symmetric with respect to the root nodes, the constructed network G t is deterministic, and vice versa. The most distinct advantage of this model is that a generator G can be chosen arbitrarily and this enables us to control the scale-free property, fractality, and other structural properties of FSFNs. Using topological information of G, we have analytically presented various indices and quantities that describe the structure of the FSFN. The obtained analytical expressions ensure that these quantities can be changed independently by varying the structure of the generator G.
We have also studied the bond-percolation problem on infinite FSFNs built by our model and computed analytically the critical point p c and various critical exponents. Furthermore, the effect of the clustering property on the percolation transition has been examined by comparing the critical points of FSFNs whose structural properties are the same as each other except for the clustering coefficient. As demonstrated by this example, the present model makes it possible to elucidate how a specific structural property influences a phenomenon occurring on FSFNs by varying systematically the structures of FSFNs.
The present model builds an FSFN by replacing every edge with a single specific generator. This model can be extended to a model in which two or more generators are employed. Networks formed by such an extended model will keep the fractal and scale-free properties. In a model with two generators, for example, an edge is replaced with a generator G 1 with the probability p or with another generator G 2 with the probability 1 − p. It is easy to show that the exponent γ is presented by Eq (17) with hm gen i and hκi instead of m gen and κ, respectively, where hm gen i is the mean number of edges and hκi is the mean degree of the root nodes of the multiple generators. In the two-generator model, these mean quantities are simply given by hm gen i = pm gen 1 + (1 − p)m gen 2 and hκi = pκ 1 + (1 − p)κ 2 , where m gen 1 (2) is the number of edges in G 1(2) and κ 1(2) is the degree of the root node of G 1 (2) . If the generators are asymmetric, κ 1(2) is the average degree of the two root nodes in G 1 (2) . The fractal dimension D f is also written as Eq (23) with hm gen i and hλi = pλ 1 + (1 − p)λ 2 , where λ 1(2) is the shortest-path distance between the root nodes in G 1 (2) . Other measures characterizing the constructed network G t are computed in similar ways to the calculations for a single asymmetric generator. Since these measures are continuous functions of the probability that a generator is adopted for an edge replacement, we can control more freely and finely the structural properties of G t by adjusting the adoption probability. The idea of constructing a network by means of mixed or probabilistic edge replacements with two kinds of small graphs has already been considered in the SHM model and the extension of the (u, v)-flower model, though the obtained network is not fractal [25,55,56]. The above extended model can be regarded as a generalization of this idea. The extension to multiple generators does not just provide a highly controllable mathematical model. The multi-generator model could be relevant to the formation mechanism of realworld FSFNs. As seen in the growth process of the World Wide Web or trading networks, many real networks grow by replacing their constituent elements with small motifs or hierarchical combinations of them. The multi-generator model suggests that networks become fractal and scale-free if the replacing procedure satisfies some conditions. Therefore, the present model and its extensions open up avenues for a systematic understanding of phenomena occurring on FSFNs and for the elucidation of formation mechanisms of real-world FSFNs.