Percolation analysis for constructing a robust modular topology based on a binary-dynamics model

In the context of Internet of Things, virtualization of wireless sensor networks is a crucial technology for sharing sensors as infrastructure. In our previous work, we proposed a brain-inspired method for constructing a robust and adaptive virtual wireless sensor network topology and showed that the method of constructing links between modules has crucial effect on robustness and adaptivity of the constructed virtual wireless sensor network topology. However, the best way of constructing a robust and adaptive virtual wireless sensor network topology is still unclear. Therefore, in this article, we use an analytical approach and propose a method for clarifying robustness of a topology according to the method of constructing links between modules. We add a new tool to a binary-dynamics model which is an analytical method for investigating percolation dynamics on a modular network. Evaluation by simulation showed that graphs in which the number of nodes selected as endpoint nodes of inter-module links and the degrees of the endpoint nodes before the link addition are large have robust connectivity in terms of the point of fragmentation of the network into modules when we fix the degree of the endpoint nodes after the link addition. After the point, the internal structure of modules may matter more. We additionally investigate an applicable range of our proposed method.


Introduction
In the context of Internet of Things, various devices, such as sensor nodes or actuator nodes, are deployed all over the world and each subset of them compose a local network. When we consider user's various service demands, these networks are required to cooperate with each other. Therefore, mechanisms for sharing networks as infrastructure are quite important. As a crucial technique for sharing substrates of networks, virtualization of wireless sensor networks (WSNs) has been attracting a great deal of attention. 1 One way of virtualization of WSN is to construct a logically connected overlay network for each application. In a virtualization scenario, multiple sensors in multiple WSNs can be used as shared infrastructure, with some sensors integrated for running each application. The virtualization of WSNs improves manageability and flexibility.
When we integrate local networks, modular structure emerges. A module consists of a group of nodes connected densely by a large number of intra-module links and modules are connected sparsely by a few number of inter-module links. Such modular structure can 1 be seen in many real networks such as Internet, social networks, or biological networks. An artificial network integrated by interdependent systems is highly vulnerable to targeted attacks or cascading failures and results in fragmentation. 2 In sensor networks, the robustness of the connectivity is a crucial issue because the robust structure of the networks leads to low overhead of maintaining stable services running over them. The study of the way of connecting sensor networks to earn the robustness thus leads to improve the cooperation of deployed sensor devices. Therefore, it is important to investigate effect of modular connection pattern on efficiency of a network, especially robust connectivity.
In our previous work, we proposed a brain-inspired method for constructing a robust and adaptive virtual wireless sensor network (VWSN) topology. 3 We showed that the method of constructing links between modules has crucial effect on robustness and adaptivity of the resulting VWSN topology. 3 However, the best way of constructing a robust and adaptive VWSN topology, and particularly of constructing links between modules, is still unclear. Toward clarifying this, we use an analytical method to study robustness of a topology according to the method of constructing links between modules. Note that we regard a network topology as an undirected graph in analytical theory.
In this article, we use an analytical approach and show the way of connecting modules so that a constructed network has the most robust connectivity. We propose a method to investigate percolation dynamics on a modular network, especially graph ensembles after addition of inter-module links. In consideration of addition of inter-module links, we add a new tool to a binary-dynamics model 4 which is an analytical method for estimating robustness of modular networks.
Our main contribution is that our analytical method considers the link addition and can be applied to make a policy for embedding a new link. Existing studies can be applied only for estimating robustness of given graph ensembles. However, our proposal enables to investigate percolation behavior according to different embedding patterns of inter-module links when a probability distribution of intra-module graph ensembles and that of inter-module link ensembles are given independently. Note that rewiring strategy in which the probability distribution does not change can make the problem for the analytical theory simple. However, link additions are more general than rewiring in the actual environment. Also, our virtual topology construction method proposed in our previous work 3 is composed of constructing intra-module topology and adding inter-module links to connect them. Therefore, in this article, we focus on the link additions for improving our method proposed in our previous work. 3 Through simulation evaluations, our analytical results are in good agreement with numerical simulations in a configuration model network. Additionally, we show that it is hard to apply our proposal to the graph in which the number of nodes is small because the target of our approach is average properties of random graph ensembles. For the similar reason, we show that the result of analysis cannot completely capture the result of a percolation process on a graph having special structural properties, such as ringshaped structure.

Related work
Many researchers have studied percolation processes on various types of graphs using a generating function approach. In this type of approach, the expected size of the giant component of a random graph ensemble can be derived from a probability distribution as given by a degree distribution or a distribution of types of links. We can then estimate robustness of a graph by evaluating percolation transitions of the size of the giant component. However, prior studies have focused on estimating robustness of given graph ensembles and not considered changes in graphs, such as link additions.
A generating function approach has been proposed for estimating robust connectivity of random graph ensembles. 5 In this method, the targeted ensemble of random graphs is defined by a generating function G(x), which represents a probability distribution, such as the degree distribution, using an auxiliary variable x. A generating function for the probability distribution of component sizes, denoted by H(x), can then be calculated from G(x). When a giant component exists, we can calculate the size of the giant component by calculating the ratio of nodes that do not belong to the giant component, from H(x). Note that, for commonly used random networks (configuration model), although the size of the giant component is related to the generating function H(x), there is no need to calculate H(x) to obtain the size of the giant component (although it is possible to do so). In the simpler and more straightforward way, it is sufficient to use the generating functions for the degree distribution (G 0 (x)) and for the excess degree distribution (G 1 (x)) to calculate the size of the giant component. In this research area, many researchers have studied percolation processes on various types of graphs, such as random graphs, 6 networks of networks, 7 multiplex networks, 8 and interdependent networks. 9 However, the generating function approach is complex because many auxiliary variables and generating functions are necessary, differing according to the complexity of the targeted graphs.
Another important analytical method is a binarydynamics model for evaluating percolation and other dynamic processes. 4 This method is relatively simple. In this method, the probability distribution of links is used to obtain the percolation behavior of the network. The probability distribution of links represents modular structure, degree-degree correlation within modules, and degree-degree correlation between modules. The type of the percolation model can be configured by changing the definition of a response function. The detail of the binary-dynamics model is shown in the next section. The binary-dynamics model can be applied to broad classes of graphs by configuring the probability distribution according to the node types or link types.
Almost all existing methods, however, focus on evaluating percolation processes on graph ensembles where a probability distribution is given, and the problem of how to embed links for constructing a robust topology is not examined. Therefore, we propose an analytical method that takes into account changes in the probability distribution due to addition of inter-module links.

Binary-dynamics model
Before we explain our proposal, we explain the binarydynamics model in detail. For convenience of explanation, we denote the type of a degree-k node belonging to module i by (i, k) and the type of a link that connects (i, k) node with (i 0 , k 0 ) node by f(i, k), (i 0 , k 0 )g. The probability distribution of links is then defined by tensor ½P i, i 0 k, k 0 in which each element represents the probability that a randomly chosen link is an f(i, k), (i 0 , k 0 )g type link.
In the binary-dynamics model, each node takes one of two states: active or inactive. An inactive node of which a neighboring node is active changes its status to active stochastically. The dynamics of binary state of nodes can be regarded as a percolation process. The probability that an inactive (i, k) node of which m neighboring nodes are active changes status to active is defined by F i (m, k). F i (m, k) is called response function. Note that we can change the way of percolation by configuring only F i (m, k).
In the binary-dynamics model, when the probability that an (i, k) node is active at time step n is denoted by q i k (n), the probability that a neighbor node of an inactive (i, k) node is active at time step n is given by Then, q i k (n) is given by where r i k (0) is the ratio of the number of active (i, k) nodes to all (i, k) nodes at the initial step. The ratio of the number of active (i, k) nodes to all (i, k) nodes at step (n + 1) is then given by From the above, the ratio of the number of active nodes to all nodes at step n, denoted as r(n), is given by where In percolation, the ratio of the active nodes to all nodes at which the dynamics (i.e. the iterative calculation of equations (1), (2), and (4)) converge describes the size of the giant component. Thus, r(n) for n ! ' describes the size of the giant component consisting of active nodes when the dynamics converge.
It is true that P i, i 0 k, k 0 for the newly created network can be easily calculated, given the adjacency matrix of the network after inter-module links addition. However, the P i, i 0 k, k 0 derived from an adjacent matrix after intermodule links addition denotes only one example, and we only get the robustness of the adjacent matrix. In this strategy, we need to try the whole connection patterns of inter-module links to make topology robust and it is only one result in the situation. This can result in tremendous overhead. In contrast, our method can get the expected robustness of the topology classified according to the strategy of inter-module links addition. Then, our approach narrows the candidates of the link addition strategy to get the robust topology with low overhead when the number of nodes is large. Therefore, our approach can show the policy to make topology robust according to the link addition strategy.

Deriving link probability distribution after addition of inter-module links
To investigate differences in robustness depending on the connection patterns between modules using the binary-dynamics model, we analyze site percolation when the connection patterns between modules are changed and the probability distribution of links within each module is given. However, we need to consider the change in the probability distribution of links within each module according to addition of inter-module links.
In this article, we consider that two previously isolated modules are connected by newly created a fixed number of inter-module links, and that these links connect a number of nodes with a fixed degree k in module i to a number of nodes with a fixed degree k 0 in module i 0 such that the degrees of nodes that receive the intermodule links are raised from k to d or from k 0 to d 0 .
First, we define tensor ½Prev i, i a, b in which each element represents the probability that a randomly chosen edge is present in the network and connects a degree-a node to a degree-b node both located in module i. It is actually the probability distribution of links considering the newly added inter-module links but also neglecting the new degrees of the boundary nodes. Because we use the probability distribution of links between modules as the target value, we define tensor ½Target i, i 0 a, b in which each element represents the probability that a randomly chosen link is an f(i, a), (i 0 , b)g type link. Because tensor ½Prev i, i a, b is for intra-module links and tensor ½Target i, i 0 a, b is for inter-module links, the conditions of equation (6) are satisfied We define T intra as the number of links within modules and T inter as the number of links between modules. Then, we configure the tensors so as to satisfy the fol- We consider the conditions specified by equations (6) and (7) to calculate the probability distribution of links after addition of inter-module links, denoted by where DPrev i, i 0 a, b denotes the amount of change in the probability distribution of links, which is what we need to calculate. Then, we use Sub i, i 0 a, b instead of P i, i 0 k, k 0 in equations (1), (4), and (5) for analysis.
In this article, we consider the case that the number of modules is two and the number of types of intermodule links is one. When the type of inter-module links is f(i, d), (i 0 , d 0 )g, it gives the following equation Focusing on module i, if we consider that the degree resulting from the addition of inter-module links is d, then we add (d À k) links to some (i, k) nodes. Then, is larger than or equal to Prev i, i d, b and Sub i, i a, d is larger than or equal to Prev i, i a, d . These relations can be hold because we first configure Prev i, i a, b and Target i, i 0 a, b so as to satisfy equations (7). This configuration enables us to assume that some amount of the probability of Prev i, i a, k ) when the intermodule links are added in the way mentioned above because some degree-k nodes change to degree-d nodes.
For example, we assume we use the topology shown in Figure 1 and add an inter-module link to a degree-1 node. In this example, Prev 1, 1 a, b , Target 1, 2 a, b , and Sub 1, 1 a, b are shown in equations (10), (11), and (12), respectively. b 0 in equation (11) is some constant value. Then, we need to derive equation (12) from equations (10) and (11) Figure 1. Example of inter-module link addition.
Next, we define the magnitude of influence by link addition. Here, the magnitude of influence means the ratio of the number of intra-module links which change their type because of the change in degree of their endpoint nodes. The magnitude of influence by link addition to a degree-k node is k times as large as that by link addition to a degree-1 node because the number of links changing their type is proportional to the degree of the node that receives the inter-module links. Also, the magnitude of influence by addition of (d À k) links to a degree-k node is 1=(d À k) times as large as that by addition of one link to a degree-k node because the number of links changing their type is proportional to the inverse of the number of links added to one node. Therefore, in the case when the number of created inter-module links is fixed and these inter-module links are connected only to nodes with a fixed degree k such that their degree is raised to d, we can define the magnitude of influence by adding (d À k) links to some degree-k nodes, denoted by Inf (d, k), as follows where H is the Heaviside step function. As the networks considered here have only two modules and also as just one fixed value of d and d 0 is considered, Target i, i 0 d, d 0 represents only one scalar quantity.
Moreover, when we add some links to a degree-k node, the probability that degree of its neighbor is k 0 is Prev i, i k, k 0 = P k 0 Prev i, i k, k 0 . When we consider the above, the amount of change in the probability distribution of links, denoted by DPrev i, i a, b (d, k), can be calculated by the following equations where d is the Kronecker delta function. When we calculate DPrev 1, 1 a, b (2, 1) for the example shown in Figure 1, DPrev 1, 1 1, 3 (2, 1) = DPrev 1, 1 3, 1 (2, 1) = À 1=26, DPrev 1, 1 2, 3 (2, 1) = DPrev 1, 1 3, 2 (2, 1) = 1=26, and the others equal to zero. This is consistent with equation (12). In another example shown in Figure 2, we can get different probability distributions according to the connected node. We assume an inter-module link is added to a degree-1 node. Prev 1, 1 a, b and Target 1, 2 a, b are shown in equations (15) and (16), respectively. b 0 in equation (16) is some constant value. In this case, (17) for the former case and in equation (18) for the latter case. We average them and can obtain the expected probability distribution shown in equation (19). When we calculate DPrev 1, 1 a, b (2, 1) for this example, the result is consistent with equation However, the case of adding links to multiple nodes is not considered in equation (14) because the equation cannot describe the change in link type from Our method cannot be directly applied to the example shown in Figure 3. We assume two inter-module links are added to two degree-2 nodes, respectively. In this example, This problem can be solved by replacing Inf (d, k) with Inf (d, k)=M i k (d) and applying equation (14) M i k (d) times with recalculating ½Prev i, i a, b each time, where M i k (d) is the number of (i, k) nodes selected as boundary nodes (i.e. endpoint nodes of inter-module links) and connected by (d À k) inter-module links. This configuration enables our method to derive equation (22) for the example shown in Figure 3.
To calculate M i k (d), we calculate the total number of inter-module links between (i, d) nodes and (i 0 , d 0 ) nodes after the link addition, denoted by M i, i 0 d, d 0 . This is given as follows where N is the total number of nodes and z prev is the average degree of the graph before adding inter-module links. For this, z prev equals P i, k kp i k where p i k is the ratio of the number of (i, k) nodes to the total number of nodes before the inter-module link additions and can be calculated as follows Note that (1=2) Nz prev describes the total number of intra-module links and (Target i, i 0 d, d 0 + Target i 0 , i d 0 , d )= P i, a, b Prev i, i a, b describes the ratio of the number of inter-module links that connect degree-d nodes from one module to degree-d 0 nodes of another module to the number of intra-module links. The second equality in equation (24) is satisfied because our target is an undirected graph. When we add all inter-module links to some degree-k nodes and add (d À k) inter-module links to each of them, M i k (d) can be calculated by the following equation Therefore, when the expected total number of nodes is given, M i k (d) can be calculated.

Constraints of input variables
In our approach, ½Target i, i 0 k, k 0 needs to satisfy some constraints. First, the number of endpoint (i, k) nodes of inter-module links must be less than the number of (i, k) nodes. Second, the number of endpoint (i, k) nodes of inter-module links must be more than (d 0 À k 0 ) and the number of endpoint (i 0 , k 0 ) nodes of inter-module links must be more than (d À k) because we assume that multiple links between a pair of nodes are not allowed. Therefore, the constraints can be described as follows We can rewrite the above into a constraint on Because the lower bound depends on N, there are cases where the desired graph cannot be constructed if N is small.

Simulation evaluation
We investigate robustness according to the connection patterns of inter-module links using our proposed method. To evaluate validity of our proposal, we compare the analysis results with the numerical results. In numerical simulations, we evaluate the site percolation process on a graph constructed by a configuration model according to ½Prev i, i a, b , ½Target i, i 0 a, b , and N. The method for constructing a graph based on ½Prev i, i a, b , ½Target i, i 0 a, b , and N is listed as follows: 1. Constructing a graph within each module (a) Calculating a degree distribution of each module from ½Prev i, i a, b and assigning degree to each node; (b) Calculating the number of each type of intra-module links from ½Prev i, i a, b and N; (c) Selecting a pair of endpoint nodes of each intra-module link from the candidates at random and connecting them.

Adding inter-module links
(a) Calculating the number of the specified type of inter-module links from the ratio of T inter to T intra ; (b) Determining a set of the candidates for the boundary nodes at random according to the connection pattern of inter-module links; (c) Selecting a pair of boundary nodes of each inter-module link from the candidates at random and connecting them.
In numerical simulations, we assume that inactive nodes are failed nodes.
Percolation on graph ensembles with a given probability distribution Simulation settings. We assume that the number of modules is two and the probability distributions of intramodule and inter-module links are given by equations (29) and (30), respectively. Equation (29) describes that the degree distribution of each module is uniform. In this case, the ratio of the number of intra-module links to inter-module links is 200 to 3. The expected total number of nodes is 1000 when we analyze the percolation process using our proposed method. We analyze the site percolation process with various combinations of d, In the percolation process, we use two modes of node removal: random failure and targeted attack. In random failure mode, a removal node is selected uniformly at random, while in targeted attack mode, a removal node is selected in order of decreasing degree. We then define a response function for each node removal mode in order to analyze the site percolation process using our method. In the binary-dynamics model, 4 the response function for the site percolation is defined as follows where Q i k is the occupation probability of (i, k) nodes. When we assume (1 À p) of all nodes have failed, then Q i k is defined by equation (32) for random failure mode 4 and by equation (33) for targeted attack mode where f i l is the ratio of the number of (i, l) nodes to the total number of nodes after the inter-module links addition and can be calculated from Sub i, i 0 a, b . P kÀ1 l = 1 P i f i l describes the ratio of the number of nodes that have lower degree than k to the total number of nodes. Therefore, the first line of equation (33) means that all nodes that have lower degree relative to the value of p are occupied. The second line means that degree-k nodes are occupied with the probability which equals to the ratio of active degree-k nodes to all degree-k nodes.
The third line means that all nodes that have higher degree relative to the value of p are not occupied.
Analysis of percolation in random failure mode. The results of analysis in random failure mode are shown in Figure 4. We evaluate the percolation with various combinations of d, d 0 , k, and k 0 . The results for the case of d = 5 and d 0 = 5 are shown in Figure 4(a) and the results for the case of d = 7 and d 0 = 7 are shown in Figure 4(b). Each legend shows the connection pattern of inter-module links and formatted as i-k-(d-k) i 0 -k 0 -(d 0 À k 0 ). In each figure, the horizontal axis shows the ratio of the failed nodes, denoted by (1 À p), and the vertical axis shows the size of the giant component. We assume that the ranking of the robustness of networks is equivalent to the ranking of the size of the giant component at each p value.
As shown in Figure 4, there is little difference depending on the connection pattern of inter-module links in random failure mode.
To evaluate validity of our analytical results shown in Figure 4, we construct a graph with the number of nodes set at 1000 and investigate the site percolation process on it in random failure mode. The number of trials is 500. Figure 5 shows the results of the site percolation process in random failure mode with k and k 0 fixed. Each figure shows that there is little difference but the graph in which the number of boundary nodes is small has a slightly more vulnerable structure. Because the graph is divided into modules when the boundary nodes are removed, the graph in which the number of boundary nodes is large has a slightly more robust connectivity.
To compare the graphs in which the number of boundary nodes is the same, we show the results of the site percolation process in random failure mode with (d À k) and (d 0 À k 0 ) fixed in Figure 6. There are small differences in each of the figures even though the number of boundary nodes is the same. This difference arises from differences in the number of neighbors of boundary nodes. Because the graph is divided and isolated when all neighbors of boundary nodes are removed, the graph in which the average number of intra-modular connections of the boundary nodes is large has a slightly more robust connectivity.
Analysis of percolation in targeted attack mode. The results of analysis in targeted attack mode are shown in Figure 7. The results for the case of d = 5 and d 0 = 5 are shown in Figure 7(a) and the results for the case of d = 7 and d 0 = 7 are shown in Figure 7(b).
In Figure 7(a), all results show that the size of the giant component decreases sharply at (1 À p) ' 0:4 and there is little difference between them. Because nodes are removed in order of decreasing degree, all degree-6 and degree-5 nodes are removed when (1 À p) is larger than 0:4. In this evaluation, the fragmentation of a graph within a module occurs before the removal of all inter-module links because d and d 0 are set to 5.
In Figure 7(b), every result shows a phase transition at small (1 À p) value, and the size of the giant component decreases gradually after that. In this evaluation, a   phase transition indicates the division of the graph into modules by removal of all inter-module links. Because the maximum degree is 7 when d and d 0 equal 7, connectivity of the graph makes it vulnerable to targeted attack. In addition, robustness of connectivity is different according to k and k 0 . The larger M i k (d) means the graph has more robust connectivity for small (1 À p) values because modules are connected until removal of all (i, d) or (i 0 , d 0 ) nodes. In this evaluation, because all (1, 7) nodes and all (2, 7) nodes are the boundary nodes and a removal node is selected uniformly and randomly from degree-7 nodes at small (1 À p), the larger (M 1 k (7) + M 2 k 0 (7)) means the graph has more robust connectivity for small (1 À p) values. Therefore, the graph with larger (1=(d À k) + 1=(d À k 0 )) has more robust connectivity in terms of the point at which the network fragments into two modules. But for large (1 À p) values, the network robustness is not significantly different for networks with different M i k (d) values. For such (1 À p) values, parameters of the intra-modular structure determine the slight differences between the sizes of the giant component for different networks exemplified. Note that the targeted attack we use in this evaluation is based on the nodes' degrees. This means these results can be different according to the definition of Q i k . To evaluate validity of our analytical results showed in Figure 7, we construct a graph with the number of nodes set at 1000 and investigate the site percolation process on it in targeted attack mode. The number of trials is 500. Figure 8 shows the results of the site percolation process in targeted attack mode. By comparison with Figure 7, although the giant component size in each method is slightly different, the order of robustness of each connection pattern is the same. From these results, analytical results are in good agreement with numerical results and the graph in which the number of boundary nodes is large has more robust connectivity  for small (1 À p) values when we use the targeted attack mode in which a removal node is selected in order of decreasing degree. After the fragmentation of the network, the network robustness is not significantly different for networks with different M i k (d) values.
Applicable range of our proposed method. When we consider the IoT environment, there may be cases where the number of nodes belonging to each module is large. Given this context, we need to investigate the applicable range of our proposed method for making policies according to the results of analysis.
To show the applicable range of our proposed method, we compare the results of analysis to the results of the percolation process on the graph when the number of nodes is one of 100, 200, 500, 1000, 10,000 and 20,000, respectively. We use the method of 1-6-1 2-6-1 as an example, the number of trials is 100, and the expected number of nodes for analysis is 1000. Figure 9(a) shows the results in random failure mode and Figure 9(b) shows the results in targeted attack mode. Both show that difference between analytical and numerical results becomes larger as the number of nodes becomes smaller. This means that it is hard to apply our proposal to a graph in which the number of nodes is small because our approach derives average properties of random graph ensembles. However, such a small graph can be analyzed by numerical simulations. Therefore, this is not a problem for our proposal whose target is the graphs in which the number of nodes is large.
In Figure 9(b), however, the results of analysis cannot completely capture any of the numerical results. In this evaluation, degree-7 nodes are 6% of the total nodes and half of them belong to module 1 and the rest belong to module 2. In the analysis method, the selection of a removal node is regarded as completely uniform. This means that the graph is not divided into modules until all degree-7 nodes are removed. In the numerical simulations, however, the graph can be divided into modules when half of degree-7 nodes are removed because removal of all degree-7 nodes belonging to module 1 or 2 results in fragmentation of the graph. Therefore, the numerical results show a phase transition at an earlier stage than the analysis. In addition, Figure 9(b) shows that a graph with a small number of nodes becomes vulnerable because the number of degree-7 nodes is small and the graph fragments easily.
The analytical results cannot completely capture any of the numerical results in targeted attack mode because of the reasons described above. The result for 10,000 nodes is almost the same as for 20,000 nodes, which means that the difference between numerical and analytical results cannot get any closer. Because this difference comes from the difference of an order of node that fails in the numerical simulation, the way that nodes fail is also an important factor which determines whether analytic results completely capture numerical results or not.
Percolation on graph ensembles with a probability distribution derived from a given intra-module topology The results shown above are obtained when we first give a probability distribution and investigate the site percolation process on a graph constructed by the configuration model to evaluate validity of our proposed method. When we consider an actual situation, however, we need to show a policy to make a graph robust by connecting multiple existing networks. Simulation settings. In this part, we investigate robustness of a graph in which modules are constructed by Erdo´ś-Re´nyi (ER) model, Baraba´si-Albert (BA) model, and Watts-Strogatz (WS) model, 10 respectively. We use the parameters shown in Table 1 for constructing a module. Using these parameters, the expectation of the average degree in a module is 6 when the number of nodes in a module is 500. The number of modules is two, the number of nodes in a module is 500, and the number of links added between modules is 1% of the total number of links within modules. We compare the numerical results for the site percolation process on a graph constructed by modules and inter-module link, with the analytical results using the probability distribution obtained from existing modules.
Results using ER model. The numerical and analytical results for the random failure mode in a network composed of ER modules are shown in Figure 10. The number of trials for numerical simulations in each topology is 100. The minimum degree within a module is 1 and the maximum degree within a module is 15.
The numerical and analytical results for the case of d = 7 and d 0 = 7 are shown in Figure 10(a) and (b), respectively. All of them are consistent with the results shown above. For the random failure mode, the analytical results depend only slightly on the connection pattern of inter-module links. In numerical results, the graph in which the number of boundary nodes is small has a slightly vulnerable structure.
When we use the ER model for constructing a topology within a module, the numerical and analytical results in targeted attack mode are shown in Figure 11. The results in Figure 11(b) are consistent with the results shown above in terms of the first sudden decay of the size of the giant component. After the fragmentation of the network, the ranking of the size of the giant component changes. This is because the connectivity of the intra-module network remains high after the fragmentation when k (or k 0 ) is set to small value. When we set k and d to 6 and 7, respectively, little number of degree-6 nodes remain after all degree-7 nodes are removed, and degree-6 nodes are crucial for the connectivity of the intra-module because of their high degrees. This leads to the vulnerable connectivity of intra-module network after the fragmentation. In most of the results of Figure 11, after removal of the boundary nodes, the internal connectivity of module 1 is larger than that of module 2. Therefore, in such cases, k is more dominant than k 0 after the fragmentation. Figure 11(a) shows that the order of robustness of each connection pattern is the same as the results of analysis although the giant component sizes are not completely the same.
Results using BA model. The numerical and analytical results for the random failure mode in a network composed of BA modules are shown in Figure 12. The number of trials for numerical simulations in each topology is 100. The minimum degree within a module is 3 and the maximum degree within a module is 76.  The numerical and analytical results for the case of d = 7 and d 0 = 7 are shown in Figure 12(a) and (b), respectively. These results show that there is little difference depending on the connection pattern of intermodule links in random failure mode. Because we use BA model for constructing an intra-module graph, it has a power-law degree distribution. Therefore, difference of robustness will not occur unless we set d and d 0 to greatly high value.
When we use the BA model for constructing a graph within a module, the numerical and analytical results in targeted attack mode are shown in Figure 13. The numerical and analytical results for the case of d = 7 and d 0 = 7 are shown in Figure 13(a) and (b), respectively. Although the sizes of the giant components are almost the same, the graph in which the number of boundary nodes is large has a slightly higher robust connectivity before the fragmentation of the network into two modules, in both analytical and numerical simulations. After the fragmentation of the network, the ranking of the sizes of the giant components changes. The reason for this is the same as that mentioned for the results of Figure 11. Results using WS model. The numerical and analytical results for the random failure mode in a network composed of WS modules are shown in Figure 14. The number of trials for numerical simulations in each topology is 100. The minimum degree within a module is 4 and the maximum degree within a module is 8.
The numerical and analytical results for the case of d = 7 and d 0 = 7 are shown in Figure 14(a) and (b), respectively. Each figure shows that the locus of giant component size in each method is almost same. However, the giant component size in numerical simulation decreases at an earlier step compared with the results of analysis. This is because the graph constructed by the WS model has a ring-shaped structure  and can fragment easily in numerical simulations, while the theory assumes random intra-module connections.
When we use the WS model for constructing a graph within a module, the numerical and analytical results in targeted attack mode are shown in Figure 15. From the results shown in Figure 15(a) and (b), the results of analysis are consistent with the results shown above. A graph constructed by method in which (1=(d Àk) + 1=(d À k 0 )) is large has high robustness for small (1 À p) values. After a graph fragments into modules, however, the giant component size in numerical simulation decreases at an earlier step compared with the results of analysis because of the same reason discussed in random failure mode.
From these results, our proposal can capture the order according to which the networks fragment into two modules. However, because we do not consider the structural properties of a graph within a module, a difference of analytical and numerical results occurs when a graph within a module has a special structure. It is our future work how to resolve this problem.

Conclusion
In this article, we propose a method for estimating robustness of graph ensembles after addition of intermodule links when the probability distribution of links within a module is given. Our proposal is to a binarydynamics model 4 and add a new tool for the model. We investigate robustness according to the connection patterns of inter-module links.
In our simulation evaluation, we compare the robustness of the networks in various connection patterns of inter-module links when we fixed the values d and d 0 that describe the degrees of the boundary nodes after the link addition. Simulation experiments showed that graphs have robust connectivity in terms of the point at which the network fragments into two modules  when the number of nodes selected as boundary nodes and the degrees of the boundary nodes before the link addition are large. After the point, the internal structure of modules may matter more. To evaluate validity of the analysis results, we evaluate the percolation process on a graph constructed by a configuration model and find that the analysis results are in agreement with the numerical results. For the targeted attack mode, although the analytical results do not match well to the numerical results, the results are in agreement qualitatively. Moreover, we investigate the applicable range of our proposed method and showed that the difference between the analysis and numerical results increases as the number of nodes decreases.
To show a policy to make a graph robust by connecting multiple existing networks, we investigate robustness of a graph composed of modules with a given internal structure. The results show that our proposal can explain the order of the ranking of robustness (measured by the removal probability at which the network fragments into two modules) observed in the numerical results. Then, the number of nodes selected as boundary nodes and the degrees of boundary nodes before the link addition should be large in terms of the point of fragmentation of the network into modules when we fix the degree of the boundary nodes after the link addition.
However, because we do not consider the structural properties of a graph within a module, a difference of analytical and numerical results occurs when a graph within a module has a special structure. It is our future work how to resolve this problem.
For further investigation, we want to consider addition of multiple types of inter-module links. When multiple types of inter-module links are added, we cannot ignore the order of addition of the inter-module links because the conditional probability of the probability distribution of links changes after each addition of an inter-module link. This analysis will be realized when we use equation (14) multiple times according to the sequence of the type of new added link. Therefore, our method can also be applied to the graph in which multiple modules exist.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by a Grant-in-Aid for JSPS Fellows JP26Á1639.