Contact graphs of disk packings as a model of spatial planar networks

Spatially constrained planar networks are frequently encountered in real-life systems. In this paper, based on a space-filling disk packing we propose a minimal model for spatial maximal planar networks, which is similar to but different from the model for Apollonian networks [J. S. Andrade, Jr. et al., Phys. Rev. Lett. {\bf 94}, 018702 (2005)]. We present an exhaustive analysis of various properties of our model, and obtain the analytic solutions for most of the features, including degree distribution, clustering coefficient, average path length, and degree correlations. The model recovers some striking generic characteristics observed in most real networks. To address the robustness of the relevant network properties, we compare the structural features between the investigated network and the Apollonian networks. We show that topological properties of the two networks are encoded in the way of disk packing. We argue that spatial constrains of nodes are relevant to the structure of the networks.


I. INTRODUCTION
The past decade has witnessed a great deal of activity devoted to complex networks by the scientific community, since many systems in the real world can be described and characterized by complex networks [1,2]. Prompted by the computerization of data acquisition and the increased computing power of computers, researchers have done a lot of empirical studies on diverse real networked systems, unveiling the presence of some generic properties of various natural and manmade networks: power-law degree distribution P (k) ∼ k −γ with characteristic exponent γ in the range between 2 and 3 [3], small-world effect including large clustering coefficient and small average path length (APL) [4], and degree correlations [5]. These findings are important for our understanding of the real-life systems, since they strongly affect almost all aspects of various dynamical processes taking place on networks [6,7,8].
In order to reproduce or explain the above-mentioned striking common features of real-life systems, there has been a concerted effort in the last few years within the physical circle and elsewhere to develop network models to uncover and understand the complexity of real systems [1,2]. In addition to the seminal Watts-Strogatz's (WS) small-world network model [4] and Barabási-Albert's (BA) scale-free network model [3], a huge variety of models and mechanisms have been proposed to mimic real-world systems, including initial attractiveness [9], aging and cost [10], fitness model [11], duplication [12], weight or traffic driven evolution [13,14], accelerating growth [15,16], coevolution [17], visibility graph [18], to name but a few. For reviews, see Refs. [1,2,6,7]. Although significant progress has been made in the field of network modeling and has led to a significant improvement in our understanding of complex systems, it is still a fundamental task and of current interest to construct models mimicking real networks and reproducing their generic properties from different angles.
Most previous network models concentrated on topological aspects and ignored the geographical effects. In these work, the node position has no special signification, which is reasonable for some networks that can be considered as lying in an abstract space, such as scientific collaboration network [19], metabolic network [20] and so forth. However, a plethora of real networks have well-defined node positions, including the Internet [21], power grid [22], airline networks [23], to name only a few. This class of networks is often referred to as geographical or spatial networks, which is a promising kind of networks, since its geography has a pivotal influence on the dynamics running on the networks, such as robustness [24], cascading breakdown [25], synchronization [26], disease spreading [27], among others. In addition to the geography, some real-world geographical networks are also planar. Typical examples are street networks [28], electronic circuit networks [29], ant-trail networks [30], and neural networks [31]. Recently several models [27,32,33,34,35,36,37] have been presented that take the the geographical position into consideration. However, models for spatial planar networks are much less, with the exception of space-filling networks-Apollonian networks [38] that have received a great amount of attention [39,40,41,42,43].
In the paper, inspired by the disk packing and the construction method of Apollonian networks [38], we propose a spatial planar network, where nodes (located at the center of disks) correspond to the disks and edges represent contact relation. The suggested network belongs to a deterministically growing class of self-similar graphs. On the basis of the recursive construction and self-similar structure of the network, we determine analytically many relevant topological properties, such as degree distribution, clustering coefficient, average path length, as well as degree correlations. The obtained results show that the network is simultaneously scale-free, small-world, and disassortative, as observed in a variety of real networks. Except for the fact that all nodes of the network have separate spatial positions that encode the network structure, we also show that the network is not only planar, but also maximally planar. Note that although Apollonian networks are also derived from a recursive disk packing, a main subject of the current paper is to compare the topological features of the presented network and Apollonian networks, with an aim to address the effect of the way of disk packings (thus the spatial constrains of nodes) on the network structure, which is important to the theory of complex networks.

II. CONSTRUCTION OF THE NETWORK
Before defining the network, we first introduce a disk packing [44,45], which is a variation of the celebrated Apollonian packing [46]. The introduced disk packing, shown in figure 1, is constructed as follows. We start with three mutually touching disks, the interstice of which is a curvilinear triangle, and we denote this initial configuration by generation t = 0. Then in the first generation t = 1, three smaller mutually touching disks are added to fill the interstice of the initial configuration, each of which touches two adjacent sides of the initial curvilinear triangle, giving rise to seven new smaller curvilinear triangles. For subsequent generations we indefinitely repeat the packing process for all the new curvilinear triangles. In the limit of infinite t generations, we obtain a disk packing.
The above-mentioned disk packing can be used as a basis for a network, which is the research object of current paper. The translation from the disk packing to network generation is quite straightforward. Let the nodes (vertices) of the network correspond to the disks and make two nodes connected if the corresponding disks are in contact. Alternatively, one can also connect the centers of the touching disks by lines to obtain the network. Figure 2 shows the network corresponding to the disk packing in figure 1.
In the construction process of the disk packing, for each interstice at arbitrary generation, once we add three disks to fill it, seven new interstices are created that will be filled in the next iteration. When building the network, it is equivalent to say that for each group of three new nodes added, seven new triangles are generated in the network, into each of which a group of three nodes will be inserted in the next iteration. According to this, we can introduce an iterative algorithm to create the network, denoted by F t after t generation evolutions. The iterative algorithm for the network is as follows: For t = 0, F 0 consists of three nodes forming a triangle. Then, we add three nodes into the original triangle. These three new nodes are linked to each other shaping a new triangle, and both ends of each edge of the new triangle are connected to a node of the original triangle. Thus we get F 1 , see figure 3. For t ≥ 1, F t is obtained from F t−1 . For each of the existing triangles of F t−1 that has never generated nodes before, we call it an active triangle. We replace each of the existing active triangles of F t−1 by the connected cluster on the right hand side of figure 3 to obtain F t . The growing process is repeated until the network reaches a desired order (node number of network). Figure 2 shows the network growing process for the first two steps.
Next we compute the order and size (number of all edges) of the network F t . Let L v (t), L e (t) and L ∆ (t) be the number of nodes, edges and active triangles created at step t, respectively. By construction (see also figure 3), each active triangle in F t−1 will be replaced by seven active triangles in F t . Thus, it is not difficult to find the following relation: Note that each active triangle in F t−1 will lead to an addition of three new nodes and nine new edges at step t, then one can easily obtain the following relations: and L e (t) = 9 L ∆ (t − 1) = 9 × 7 t−1 for arbitrary t > 0. From these results, we can compute the order and size of the network. The total number of vertices N t and edges E t present at step t is and respectively. So for large t, the average degree k t = 2Et Nt is approximately 6, which shows the network is sparse as most real systems.
From equations. (1) and (2), we have E t = 3N t − 6. In addition, by the very construction of the network, it is obvious that arbitrary two edges in the network never cross each other. Thus our network is a maximal planar network (or graph) [47].

III. STRUCTURAL PROPERTIES OF THE NETWORK
In this section, we study the statistical properties of the network, in terms of degree distribution, clustering coefficient, average path length, and degree correlations.

A. Degree distribution
When a new node i is added to the graph at step t i (t i ≥ 1), it has a degree of 4. Let L ∆ (i, t) be the number of active triangles at step t that will create new nodes connected to node i at step t + 1. Then at step t i , L ∆ (i, t i ) = 4. From the iterative generation process of the network, one can see that at any step each two new neighbors of i generate three new active triangles involving i, and one of its existing active triangle is deactivated simultaneously. We define k i (t) as the degree of node i at time t, then the relation between k i (t) and L ∆ (i, t) satisfies: (3) It should be mentioned that the initial three vertices created at step 0 have a little different evolution process from other ones. We can easily obtain: L ∆ (0, t) = 3 t and k i (t) = 3 t + 1. Thus, at step t, the degree of the initial three nodes is less than that of those three nodes born at step 1 but larger than that for those nodes emerging at other steps.
Equation (4) shows that the degree spectrum of the network is discrete. It follows that the cumulative degree distribution [6] is given by which is valid for all t i ≥ 2.
When t is large enough, one can obtain So the degree distribution follows a power law form with the exponent γ = 1 + ln 7 ln 3 .

B. Clustering coefficient
The clustering coefficient [4] C i of node i is defined as the ratio between the number of edges e i that actually exist among the k i neighbors of node i and its maximum possible value, k i (k i − 1)/2, i.e., The clustering coefficient of the whole network is the average of C ′ i s over all nodes in the network. For our network, the analytical expression of clustering coefficient C(k) for a single node with degree k can be derived exactly. When a node enters the system, both k i and e i are 4. In the following iterations, each of its active triangles increases both k i and e i by 2 and 3, respectively. Thus, e i is equal to 4 + 3 2 (k i − 4) for all nodes at all steps. So one can see that there exists a one-to-one correspondence between the degree of a node and its clustering. For a node of degree k, we have In the limit of large k, C(k) is inversely proportional to degree k. The same scaling of C(k) ∼ k −1 has also been observed in several real-life networks [48]. Using equation (8), we can obtain the clustering C t of the network at step t: where the sum runs over all the nodes and D r is the degree of the nodes created at step r, which is given by equation (4). In the infinite network order limit (N t → ∞), equation (9) converges to a nonzero value C = 0.603, as shown in figure 4. Therefore, the average clustering coefficient of the network is very high.

C. Average path length
We represent all the shortest path lengths of network F t as a matrix in which the entry d ij is the distance between node i and j that is the length of a shortest path joining i and j. A measure of the typical separation between two nodes in F t is given by the average path length (also called average distance) d t defined as the mean of distances over all pairs of nodes: where denotes the sum of the distances between two nodes over all couples. Note that in equation (11), for a pair of nodes i and j (i = j), we only count d ij or d ji , not both. (η = 1, 2, · · · , 7), which are connected to one another at the six edge nodes, i.e., A, B, C, X, Y , and Z.

Recursive equation for total distances
We continue by exhibiting the procedure of determining the total distance and present the recurrence formula, which allows us to obtain D t+1 of the t + 1 generation from D t of the t generation. The network F t under consideration has a self-similar structure that allows one to calculate D t analytically [49,50,51]. As shown in figure 5, network F t+1 may be constructed by joining at six edge nodes (i.e., A, B, C, X, Y , and Z) seven copies of F t that are labeled as F t . According to the second construction method, the total distance D t+1 satisfies the recursion relation where ∆ t is the sum over all shortest path length whose endpoints are not in the same F (η) t branch. The last term -9 on the right-hand side of equation (12) compensates for the overcounting of certain paths: the shortest path d AX between A and X, with length 1, is included in both F (1) t and F (2) t ; similarly, the shortest paths d AY , d BX , d BZ , d CY , d CZ , d XY , d XZ , and d Y Z are all computed twice. To determine D t , all that is left is to calculate ∆ t .

Definition of crossing distance
In order to compute ∆ t , we classify the nodes in F t+1 into two categories: the six edge nodes (such as A, B, C, X, Y , and Z in figure 5) are called connecting nodes, while the other nodes are named non-connecting nodes. Thus ∆ t , named the crossing distance, can be obtained by summing the following path length that are not included in the distance of node pairs in F (η) t : length of the shortest paths between non-connecting nodes, length of the shortest paths between connecting and non-connecting nodes, and length of the shortest paths between connecting nodes (i.e., d AZ , d BY , and d CZ ).
Denote ∆ α,β t as the sum of all shortest paths between non-connecting nodes, whose end-points are in F be the set of non-connecting nodes in F t+1 , which belong to F (η) t . Then the total sum ∆ t is given by By symmetry, equation (13) can be simplified as ∆ t = 9 ∆ 1,2 t + 9 ∆ 1,3 t + 3 ∆ 1,4 t + 21 To calculate the crossing distance ∆ 1,2 t , ∆ 1,3 t , ∆ 1,4 t , and j∈Ω 4 t d Aj , we classify interior nodes in network F t+1 into seven different parts according to their shortest path lengths to each of the three peripheral nodes (i.e. A, B, C). Notice that nodes A, B, C themselves are not partitioned into any of the seven parts represented as P 1 , P 2 , P 3 , P 4 , P 5 , P 6 , and P 7 , respectively. The classification of nodes is shown in figure 6. For any interior node v, we denote the shortest path lengths from v to A, B, C as a, b, and c, respectively. By construction, a, b, c can differ by at most 1 since vertices A, B, C are adjacent. Then the classification function class(v) of node v is defined to be for a < b = c, P 2 for b < a = c, P 3 for c < a = b, P 4 for a = c < b, P 5 for a = b < c, P 6 for b = c < a, P 7 for a = b = c.
It should be mentioned that the definition of node classification is recursive. For instance, class P 1 and P 4 in F (1) t belong to class P 1 in F t+1 , class P 3 and P 5 in F (1) t belong to class P 2 in F t+1 , class P 2 , P 6 , and P 7 in F (1) t belong to class P 5 in F t+1 . Since the three nodes A, B, and C are symmetrical, in the network we have the following equivalent relations from the viewpoint of class cardinality: classes P 1 , P 2 , and P 3 are equivalent to one another, and it is the same with classes P 4 , P 5 , and P 6 . We denote the number of nodes in network F t that belong to class P 1 as N t,P1 , the number of nodes in class P 2 as N t,P2 , and so on. By symmetry, we have N t,P1 = N t,P2 = N t,P3 and N t,P4 = N t,P5 = N t,P6 . Therefore in the following computation we will only consider N t,P1 , N t,P4 , and N t,P7 . It is easy to conclude that Considering the self-similar structure of the network, we can easily know that at time t + 1, the quantities N t+1,P1 , N t+1,P4 , and N t+1,P7 evolve according to the following recursive equations where we have used the equivalent relations N t,P1 = N t,P2 = N t,P3 and N t,P4 = N t,P5 = N t,P6 . With the initial condition N 2,P1 = 4, N 2,P4 = 2, and N 2,P7 = 6, we can solve the recursive equation (17) to obtain For a node v in network F t+1 , we are also interested in the smallest value of the shortest path length from v to any of the three peripheral nodes A, B, and C. We denote the shortest distance as f v , which can be defined to be min(a, b, c).

Calculation of crossing distances
Having obtained the quantities N t,Pi and d t,Pi (i = 1, 2, · · · , 7), we now begin to determine the crossing distance ∆ 1,2 of ∆ 1,2 Proceeding similarly, we obtain and With the obtained results for δ i t , we have Analogously, we find and Substituting equations (28), (29), (30), and (31) into equation (14), we the final expression for cross distances ∆ t ,

Exact result for average path length
With the above-obtained results and recursion relations, we now readily calculate the sum of the shortest path lengths between all pairs of nodes. Inserting equation (32) into equation (12) and using the initial condition D 2 = 717, equation (12) is solved inductively, Substituting equation (33) into equation (10) In the large t limit, d t ∼ t, while the network order N t ∼ 7 t which is obvious from equation (1). Thus, the average path length grows logarithmically with increasing order of the network. We have checked our analytic result provided by equation (34) against numerical calculations for different network order up to t = 8 which corresponds to N 8 = 1 007 772. In all the cases we obtain a complete agreement between our theoretical formula and the results of numerical investigation, see figure 7.
Recently, it has been suggested that for random uncorrelated scale-free networks (SFNs) with degree exponent γ < 3 and network order N , their average distance d(N ) behaves as a double logarithmic scaling with N : d(N ) ∼ ln ln N [52,53]. However, for the deterministic network considered here, in despite of the fact that its degree exponent γ = 1 + ln 7 ln 3 < 3, its average path length scales as a logarithmic scaling with network order, showing a obvious difference from that of the stochastic scale-free counterparts. The logarithmic scaling of d t with N t for our network as well as the Apollonian networks [41] shows that previous relation between APL and the network order obtained for uncorrelated SFNs [52,53] is not valid for disassortative SFNs [52,53], at least for some spatial networks, e.g., the Apollonian networks and the network considered here. This leads us to the conclusion that degree exponent itself does not suffice to characterize the APL of SFNs.

D. Degree correlations
An interesting quantity related to degree correlations is the average degree of the nearest neighbors for nodes with degree k, denoted as k nn (k), which is a function of node degree k [54,55]. When k nn (k) increases with k, it means that nodes have a tendency to connect to nodes with a similar or larger degree. In this case the network is defined as assortative [5]. In contrast, if k nn (k) is decreasing with k, which implies that nodes of large degree are likely to have near neighbors with small degree, then the network is said to be disassortative. If correlations are absent, k nn (k) = const. We can exactly calculate k nn (k) for network F t using equations (3) and (4) to work out how many links are made at a particular step to nodes with a particular degree. By construction, we have the following expression [56,57] for k = 4 × 3 t−ti (t i ≥ 1), where k(t i , t) is the degree of a node i at time t that was born at step t i . Here the first sum on the right-hand side accounts for the links made to nodes with larger degree (i.e. t ′ i < t i ) when the node was generated at t i . The second sum describes the links made to the current smallest degree nodes at each step t ′ i > t i . The last term 2 accounts for the two links connected to two simultaneously emerging nodes. After some algebraic manipulations, we can rewrite equation (35) in term of k to obtain k nn (k) = (3 2 t+1 + 3 t−1 ) 4 k 2−ln 7/ ln 3 2 × 7 t−1 + 8 ln( k 4 ) 3 ln 3 − 10.
For k = 3 t + 1 (t i = 0), we have Therefore, for large t and k, k nn (k) is approximately a power law function of k as k nn (k) ∼ k −ω with ω = 2 − ln 7 ln 3 ≃ 0.229, which shows that the network is disassortative. Note that k nn (k) of the Internet exhibits a similar power-law scaling with exponent ω = 0.5 [54].

IV. CONCLUSION
In summary, motivated by the disk packing and Apollonian networks, we have presented a model for spatial planar networks introducing the influence of geography encoded in the disk packing. According to the construction, we have studied analytically the main structural features of the network. We have shown that the network has a power law distribution with exponent γ = 1 + ln 7 ln 3 , it has a large clustering coefficient 0.603, its APL scales logarithmically with the number of network nodes, and it is disassortative with the average degree of the nearest neighbors for nodes having degree k being roughly a power-law function of k with exponent -0.229.
Note that although both the network considered here and the Apollonian network [38,56] are translated from disk packings, and both networks have qualitatively similar topologies, their structural characteristics are quantitatively different. For example, the exponent of degree distribution is 1 + ln 7 ln 3 for our network, while for Apollonian network it is 1 + ln 3 ln 2 ; the average clustering coefficients for our network and Apollonian network are 0.603 and 0.828, respectively. In addition, the average path length [41] and the degree correlations [56] for both networks are also of quantitative difference. These disparities of the two networks show that the ways of disk packing lead to different spatial constrains of network nodes, which in turn have a significant impact on network properties and thus dynamics running on networks, such as cascading failing [25], random walks [58], and so on. Thus, we can conclude that for spatial networks the positions, where nodes are geographically located, matter greatly and should be incorporated when modeling such networks. Ignoring the geography will lead to miss some important attributes and properties of the systems.