Research Article

Percolation analysis for constructing arobust modular topology based on abinary-dynamics modelInternational Journal of Distributed
Sensor Networks
2017, Vol. 13(4)Ó The Author(s) 2017
DOI: 10.1177/1550147717701141
journals.sagepub.com/home/ijdsn

Shinya Toyonaga1, Daichi Kominami2 and Masayuki Murata1

Abstract
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 forinvestigating percolation dynamics on a modular network. Evaluation by simulation showed that graphs in which thenumber 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.
Keywords
Percolation, analysis, modular structure, robust connectivity, graph theory

Date received: 30 May 2016; accepted: 2 March 2017
Academic Editor: Valerio Freschi

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 cooperatewith each other. Therefore, mechanisms for sharing
networks as infrastructure are quite important. As acrucial technique for sharing substrates of networks,
virtualization of wireless sensor networks (WSNs) has
been attracting a great deal of attention.1 One way ofvirtualization 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

1Graduate School of Information Science and Technology, Osaka
University, Suita, Japan
2
Graduate School of Economics, Osaka University, Suita, Japan
Corresponding author:
Shinya Toyonaga, Graduate School of Information Science andTechnology, Osaka University, Suita 565-0871, Osaka, Japan.
Email: s-toyonaga@ist.osaka-u.ac.jp

Creative Commons CC-BY: This article is distributed under the terms of the Creative Commons Attribution 4.0 License(http://www.creativecommons.org/licenses/by/4.0/) which permits any use, reproduction and distribution of the work without
further permission provided the original work is attributed as specified on the SAGE and Open Access pages (http://www.uk.sagepub.com/aboutus/
openaccess.htm).

2
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 resultsin fragmentation.2 In sensor networks, the robustnessof the connectivity is a crucial issue because the robust
structure of the networks leads to low overhead ofmaintaining 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 toinvestigate 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 adaptivityof 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 robustnessof 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 andshow 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 abinary-dynamics model4 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 makea policy for embedding a new link. Existing studies canbe applied only for estimating robustness of given
graph ensembles. However, our proposal enables toinvestigate 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 theprobability distribution does not change can make theproblem 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 work3 is composedof constructing intra-module topology and adding
inter-module links to connect them. Therefore, in thisarticle, we focus on the link additions for improving
our method proposed in our previous work.3
Through simulation evaluations, our analyticalresults are in good agreement with numericalInternational Journal of Distributed Sensor Networks
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 issmall because the target of our approach is average
properties of random graph ensembles. For the similarreason, we show that the result of analysis cannot completely capture the result of a percolation process on agraph 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 ofthe giant component of a random graph ensemble canbe derived from a probability distribution as given by adegree 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 notconsidered changes in graphs, such as link additions.
A generating function approach has been proposed
for estimating robust connectivity of random graphensembles.5 In this method, the targeted ensemble ofrandom graphs is defined by a generating function
G(x), which represents a probability distribution, suchas the degree distribution, using an auxiliary variable x.
A generating function for the probability distributionof 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 thesize of the giant component is related to the generating
function H(x), there is no need to calculate H(x) toobtain the size of the giant component (although it ispossible to do so). In the simpler and more straightforward way, it is sufficient to use the generating functions
for the degree distribution (G0 (x)) and for the excess
degree distribution (G1 (x)) to calculate the size of the
giant component. In this research area, many researchers have studied percolation processes on various typesof graphs, such as random graphs,6 networks of networks,7 multiplex networks,8 and interdependent networks.9 However, the generating function approach iscomplex because many auxiliary variables and generating functions are necessary, differing according to thecomplexity of the targeted graphs.
Another important analytical method is a binarydynamics model for evaluating percolation and otherdynamic processes.4 This method is relatively simple. In
this method, the probability distribution of links is used

Toyonaga et al.

3

to obtain the percolation behavior of the network. The
probability distribution of links represents modular
structure, degree–degree correlation within modules, anddegree–degree correlation between modules. The type of
the percolation model can be configured by changing
the definition of a response function. The detail of thebinary-dynamics model is shown in the next section. The
binary-dynamics model can be applied to broad classesof 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 wherea probability distribution is given, and the problem of
how to embed links for constructing a robust topology
is not examined. Therefore, we propose an analyticalmethod that takes into account changes in the probability distribution due to addition of inter-module links.

where rik (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 atstep (n + 1) is then given byrik (n + 1) = rik (0) + (1  rik (0))

m=0 m
m
i
i
3 (qk (n)) (1  qk (n))km Fi (m, k)

ð3Þ

From the above, the ratio of the number of active
nodes to all nodes at step n, denoted as r(n), is given by
P Pi,k,ik0 0r(n) =

X i0 , k, k 0
i

k

0

P

Pi,k,ik 0

i, i0 , k, k 0

ri (n)

ð4Þk

P Pi,k,ik0 0

Methodwhere ri (n) =Binary-dynamics modelPqik (n) =X i0 , k 0kBefore 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 (i0 , k 0 ) node by f(i, k), (i0 , k 0 )g. The probability
0
distribution of links is then defined by tensor ½Pi,k,ik 0  inwhich each element represents the probability that a randomly chosen link is an f(i, k), (i0 , k 0 )g type link.
In the binary-dynamics model, each node takes oneof two states: active or inactive. An inactive node ofwhich a neighboring node is active changes its status toactive stochastically. The dynamics of binary state ofnodes can be regarded as a percolation process. The
probability that an inactive (i, k) node of which mneighboring nodes are active changes status to active isdefined by Fi (m, k). Fi (m, k) is called response function.
Note that we can change the way of percolation byconfiguring only Fi (m, k).
In the binary-dynamics model, when the probability
that an (i, k) node is active at time step n is denoted byqik (n), the probability that a neighbor node of an inactive (i, k) node is active at time step n is given by

i0 , k 000Pi,k,ik 0 qik 0 (n)Pi0 , k 00

Pi,k,ik 0

ð1Þ

Then, qik (n) is given byqik (n + 1) = rik (0) + (1  rik (0))

k1 
X
k1
qik (n))k1m Fi (m, k) ð2Þ
3 (qik (n))m (1  mm=0qik (0) = rik (0)


k 
X
k

k

P Pi,k,ik0 0
i0 , k, k 0

rik (n)ð5Þk

In percolation, the ratio of the active nodes to allnodes 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 ofactive nodes when the
dynamics converge.
0
It is true that Pi,k,ik 0 for the newly created network canbe easily calculated, given the adjacency matrix of the
network0 after inter-module links addition. However,the Pi,k,ik 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 resultin 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 ofinter-module links
To investigate differences in robustness depending on
the connection patterns between modules using thebinary-dynamics model, we analyze site percolation
when the connection patterns between modules arechanged and the probability distribution of links withineach module is given. However, we need to consider the4International Journal of Distributed Sensor Networks

change in the probability distribution of links withineach 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 modulei0 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 ½Previ,a,ib  in which each element represents the probability that a randomly chosenedge is present in the network and connects a degree-anode 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 0modules asi, i
the target value, we define tensor ½Targeta,b  in which
each element represents the probability that a randomly
chosen link is an f(i, a), (i0 , b)g type link. Because tensori, i0½Previ,a,ib  is for intra-module links and tensor ½Targeta,b
is for inter-module links, the conditions of equation (6)are satisfied
0Previ,a,ib = 0, for i 6¼ i0
0

i, i
0
Targeta,b = 0, for i = i

ð6Þ

We define Tintra as the number of links within modules and Tinter as the number of links between modules.
Then, we configure the tensors so as to satisfy the following equations
X

0

i, i
Targeta,b=1

i, i0 , a, b

i, a, b

X

X

Previ,a,ib +Previ,a,ib :

X

0

i, i
Targeta,b = Tintra : Tinterð7Þ

i, i0 , a, b

i, a, b

We consider the conditions specified by equations(6) and (7) to calculate the probability distribution oflinks after addition of inter-module links, denoted by
0
0
i, i0
tensor ½Subi,a,ib , using ½Previ,a,ib  and ½Targeta,
b .
i, i0

Therefore, Suba, b is given by(0
Subi,a,ib =

0

0

0

Previ,a,ib + DPrevi,a,ib
i, i0
Targeta,b

where i = i0
otherwise

ð8Þ

where DPrevi,a,ib denotes the amount of change in theprobability distribution of links, which
is what we need00
to calculate. Then, we use Subi,a,ib instead of Pi,k,ik 0 in
equations (1), (4), and (5) for analysis.
In this article, we consider the case that the numberof modules is two and the number of types of intermodule links is one. When the type of inter-module
links is f(i, d), (i0 , d 0 )g, it gives the following equation

Figure 1. Example of inter-module link addition.

i, i0Targeta,b=

TinterTintra + Tinterwhere a = d, b = d 0
otherwise

0ð9Þ

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,
Subi,k,ib is smaller than or equal to Previ,k,ib , and Subi,a,ik issmaller than or equal to Previ,a,ik . Simultaneously, Subi,d,ib
is larger than or equal to Previ,d,ib and Subi,a,id is larger
than or equal to Previ,a,id . These relations can0 be holdi, ibecause we first configure Previ,a,ib and Targeta,b so as tosatisfy equations (7). This configuration enables us toassume that some amount of the probability of Previ,a,ik(Previ,k,ib ) moves to Previ,a,id (Previ,d,ib ) 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 shownin Figure 1 and add an inter-module link to a degree-1
1
1, 2
1, 1
node. In this example, Prev1,
a, b , Targeta, b , and Suba, b are
shown in equations (10), (11), and (12), respectively. b0
in equation (11) is some constant value. Then, we needto derive equation (12) from equations (10) and (11)b = 1, 2, 301
0 0 4 a=1
ð10Þ
1 B
C
1
Prev1,@0 0 2Aa=2a, b =
26
4 2 0 a=3(
1
where a = 2, b = b0
1, 2
2, 1
Targeta,26
b = Targetb, a =0otherwiseð11Þ
b = 1, 2, 3
1
Sub1,
a, b =

0

0

1 B
@0
26
3

1

a=1
C
0 3Aa=2
3 0 a=3
0 3

ð12Þ

Toyonaga et al.

5

Next, we define the magnitude of influence by link
addition. Here, the magnitude of influence means theratio 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 bylink addition to a degree-1 node because the number oflinks changing their type is proportional to the degreeof 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 byaddition of one link to a degree-k node because thenumber of links changing their type is proportional tothe inverse of the number of links added to one node.
Therefore, in the case when the number of createdinter-module links is fixed and these inter-module links
are connected only to nodes with a fixed degree k suchthat their degree is raised to d, we can define the magnitude of influence by adding (d  k) links to somedegree-k nodes, denoted by Inf (d, k), as follows
Inf (d, k) = H(d  k)ki, i0
Targetd,d0(d  k)ð13Þ

where H is the Heaviside step function. As the networks
considered here have only two modules and also as
just
i, i0
repone fixed value of d and d 0 is considered, Targetd,
0dresents only one scalar quantity.
Moreover, when we add some links to a degree-knode, thePprobability that degree of its neighbor is k 0 isi, iPrevk, k 0 = k 0 Previ,k,ik 0 . When we consider the above, theamount of change in the probability distribution oflinks, denoted by DPrevi,a,ib (d, k), can be calculated by
the following equations
Previ,k,ib
DPrevi,a,ib (d, k) =  da, k P
i, i  Inf (d, k)
k 0 Prevk, k 0Previ,a,ik
 db, k P
i, i  Inf (d, k)
k 0 Prevk, k 0Previ,k,ib+ da, d Pi, i  Inf (d, k)
k 0 Prevk, k 0

Figure 2. Example of link addition in which different
probability distributions arise according to the connected node.

shown in equations (15) and (16), respectively. b0 inequation (16) is some constant value. In this case,f(1, 1), (1, 3)g link changes its type to f(1, 2), (1, 3)g andf(1, 3), (1, 1)g link changes its type to f(1, 3), (1, 2)gwith probability 2=3. Also, f(1, 1), (1, 2)g link changes
its type to f(1, 2), (1, 2)g and f(1, 2), (1, 1)g link changes
1
its type to f(1, 2), (1, 2)g with probability 1=3. Sub1,a, b isshown in equation (17) for the former case and in equation (18) for the latter case. We average them and canobtain the expected probability distribution shown in1equation (19). When we calculate DPrev1,
a, b (2, 1) for thisexample, the result is consistent with equation (19)b = 1, 2, 3,

0

1

a=1
ð15Þ
1 B
C
@1 0 1Aa=218
2 1 0 a=3
1
0
1, 2
2, 1
18 where a = 2, b = b
Targeta,b = Targetb, a =0otherwise
0

1

2

1
Prev1,a, b =

ð16Þ
ð14Þ

Previ,a,ik+ db, d Pi, i  Inf (d, k)
k 0 Prevk, k 0

where d is the Kronecker delta function.
1
When we calculate DPrev1,
a, b (2, 1) for the example
1
1, 1
shown in Figure 1, DPrev1,
1, 3 (2, 1) = DPrev3, 1 (2, 1)
1
1, 1=  1=26, DPrev1,
2, 3 (2, 1) = DPrev3, 2 (2, 1) = 1=26, andthe 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
1
1, 2
added to a degree-1 node. Prev1,
a, b and Targeta, b are

b = 1, 2, 3
1
Sub1,
a, b =

0

1 B
@1
18
1

b = 1, 2, 3
1
Sub1,
a, b =

0

0

0

1 B
@0
18
2

b = 1, 2, 3

0

0
1 B
1
Sub1,
=
@2
a, b
54
4

1

2

a=1
C
2Aa=2
0 a=3

0

2

1
0

1

ð17Þ

1

2
1

a=1
C
1Aa=2
0 a=3

2
2

1
4 a=1
C
5Aa=2

5

0

a=3

ð18Þ

ð19Þ

6

International Journal of Distributed Sensor Networks
Mki (d) is the number of (i, k) nodes selected as boundary
nodes (i.e. endpoint nodes of inter-module links) andconnected by (d  k) inter-module links. This configuration enables our method to derive equation (22) forthe example shown in Figure 3.
To calculate Mki (d), we calculate the total number ofinter-module links between (i, d) nodes and(i0 , d 0 ) nodes
i, i0
after the link addition, denoted by Md, d 0 . This is given
as follows
i, i0
i0 , i(Targetd,1d 0 + Targetd 0 , d )i, i0Md, d 0 = Nzprev
P
02
Previ,a,ib
i, i0 , a, b0

= Nzprev
Figure 3. Example of multiple links addition in which different
probability distributions arise according to the connected node.

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 fromf(i, k), (i, k)g to f(i, d), (i, d)g even though it can
describe the change from f(i, k), (i, b)g to f(i, d), (i, b)g
and from f(i, a), (i, k)g to f(i, a), (i, d)g. 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,
1
1, 2
Prev1,a, b and Targeta, b are shown in equations (20) and(21), respectively. b0 in equation (21) is some constantvalue. The expected probability distribution after intermodule link addition is shown in equation (22). When1
we calculate DPrev1,
a, b (3, 2) for this example, equation
1(23) is derived and Sub1,
3, 3 remains zerob = 2, 3
1
Prev1,a, b =

1
20


(

1, 2
2, 1
Targeta,b = Targetb, a =b = 2, 3
1
Sub1,a, b =

1
60

0

0

0



a=2

ð20Þ

a=3

2
where a = 3, b = b0200otherwise
ð21Þ
4

8

8

4



a=2

1
20



where N is the total number of nodes and zprev is theaverage degree of the graphPbefore adding inter-module
links. For this, zprev equals i, k kpik where pik is the ratioof the number of (i, k) nodes to the total number ofnodes before the inter-module link additions and canbe calculated as follows
P Previ,k,ib
pik =

0

4

4

0

a=2

k

b

P

0

i, i0 , a, b

Previ,a,ibað25Þ

Note that (1=2) Nzprev describes the total
numberofi, i0i0 , iintra-module
links
and
(Target0 + Targetd 0 , d )=
d,
d
P
i, i
i, a, b Preva, b describes the ratio of the number ofinter-module links that connect degree-d nodes fromone module to degree-d 0 nodes of another module tothe number of intra-module links. The second equality
in equation (24) is satisfied because our target is anundirected graph. When we add all inter-module links
to some degree-k nodes and add (d  k) inter-module
links to each of them, Mki (d) can be calculated by thefollowing equation0

Mki (d) =Md,i, id 0
dk

ð26Þ

Therefore, when the expected total number of nodes
is given, Mki (d) can be calculated.

ð22Þ

a=3


ð24Þ

i, i0 , a, bConstraints of input variables0b = 2, 3
1
Sub1,a, b =



8i, i
Targetd,d0P0Previ,a,ib

ð23Þ

a=3

This problem can be solved by replacing Inf (d, k)with Inf (d, k)=Mki (d) and applying equation (14) Mki (d)times with recalculating ½Previ,a,ib  each time, wherei, i
In our approach, ½Targetk,
k 0  needs to satisfy some constraints. First, the number of endpoint (i, k) nodes ofinter-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 ) andthe number of endpoint (i0 , k 0 ) nodes of inter-module
links must be more than (d  k) because we assume that

Toyonaga et al.

7

multiple links between a pair of nodes are not allowed.
Therefore, the constraints can be described as follows(d 0  k 0 )  Mki (d)  Npik
i0
k0

ð27Þ

0

(d  k)  M (d 0 )  Npik 0

We can rewrite the above into a constraint oni, i0
Targetd,d 0 as follows
P
00 i, a, b(d  k)(d  k )Previ,a,ib

N zprev

0i, i Targetd,d0PPrevi,a,ib


0
i,
a,
b
 min (d  k)pik , (d 0  k 0 )pik 0
zprev

ð28Þ

Because the lower bound depends on N, there arecases 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 configurationi, i0
model according to ½Previ,a,ib , ½Targeta,b , and N.
The method for constructing a graph based oni, i0½Previ,a,ib , ½Targeta,b , and N is listed as follows:1.

2.

Constructing a graph within each module(a) Calculating a degree distribution of each
module from ½Previ,a,ib  and assigningdegree to each node;(b) Calculating the number of each type of
intra-module links from ½Previ,a,ib  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 ratioof Tinter to Tintra ;(b) Determining a set of the candidates forthe boundary nodes at random according
to the connection pattern of inter-module
links;(c) Selecting a pair of boundary nodes ofeach 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 thatthe 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 combinationsof d, d 0 , k, and k 0
b = 2,

3,4, 5, 61
0
4 6 8 10 12 a = 2CB
B 6 9 12 15 18 C a = 3
C
1 B
i, i
B 8 12 16 20 24 C a = 4
Preva, b =
C
B
812 B
C
@ 10 15 20 25 30 A a = 5

ð29Þ

12 18 24 30 36 a = 6
 6
where a = d, b = d 0
1, 2
2, 1
Targeta, b = Targetb, a = 812
0otherwise
ð30Þ

In the percolation process, we use two modes ofnode removal: random failure and targeted attack. In
random failure mode, a removal node is selected uniformly at random, while in targeted attack mode, aremoval 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 isdefined as follows

Fi (m, k) =0
Qik

where m = 0
otherwise

ð31Þ

where Qik is the occupation probability of (i, k) nodes.
When we assume (1  p) of all nodes have failed, then
Qik is defined by equation (32) for random failure mode4
and by equation (33) for targeted attack mode
Qik = p

Qik =

8
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
:

where

1
p

k 1 P
P
i

P ii

l=1

fk

i

fl

ð32Þ
k P
P

fil  p

1 i
l =
kP
1 P

fil \p
where l =k1 i
PP i
\
fl
l=1 i

0

otherwise

ð33Þ

8
where fil is the ratio of the number of (i, l) nodes to thetotal number of nodes after the inter-module
addiPlinks
0
k1 P
i
tion and can be calculated from Subi,a,ib .
i fll=1
describes the ratio of the number of nodes that havelower degree than k to the total number of nodes.
Therefore, the first line of equation (33) means that allnodes 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.
Here, we consider that analytical results are obtainedby iterative calculation of equations (1), (2), and (4).
Analysis of percolation in random failure mode. The resultsof 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

Figure 4. Analytical results of percolation analysis in random
failure mode: (a) d = 5, d0 = 5 and (b) d = 7, d0 = 7.

International Journal of Distributed Sensor Networksof d = 5 and d 0 = 5 are shown in Figure 4(a) and theresults for the case of d = 7 and d 0 = 7 are shown in
Figure 4(b). Each legend shows the connection patternofinter-module
links
and
formatted
asi-k-(d-k) i0 -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 sizeof 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 shownin Figure 4, we construct a graph with the number ofnodes set at 1000 and investigate the site percolation
process on it in random failure mode. The number oftrials 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 theboundary nodes are removed, the graph in which thenumber of boundary nodes is large has a slightly morerobust connectivity.
To compare the graphs in which the number ofboundary nodes is the same, we show the results of thesite 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 ofboundary 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 islarge has a slightly more robust connectivity.
Analysis of percolation in targeted attack mode. The resultsof 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 andthere 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 agraph within a module occurs before the removal of allinter-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Toyonaga et al.

9

Figure 5. Numerical results of the site percolation process on graphs in random failure mode with k and k0 fixed: (a) k = 2, k0 = 2,(b) k = 3, k0 = 3, (c) k = 4, k0 = 4, and (d) k = 5, k0 = 5.

Figure 6. Numerical results of the site percolation process on graphs in random failure mode with (d  k) and (d0  k0 ) fixed: (a)d  k = 1, d0  k0 = 1, (b) d  k = 2, d0  k0 = 2, (c) d  k = 3, d0  k0 = 3, and (d) d  k = 4, d0  k0 = 4.

10

Figure 7. Analytical results of percolation analysis in targeted
attack mode: (a) d = 5, d0 = 5 and (b) d = 7, d0 = 7.

phase transition indicates the division of the graph intomodules 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 Mki (d) means thegraph has more robust connectivity for small (1  p)values because modules are connected until removal ofall (i, d) or (i0 , 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(Mk1 (7) + Mk20 (7)) means the graph has more robust
connectivity for small (1  p) values. Therefore, thegraph with larger (1=(d  k) + 1=(d  k 0 )) has morerobust 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 Mki (d)values. For such (1  p) values, parameters of theInternational Journal of Distributed Sensor Networks

Figure 8. Numerical results of the site percolation process ongraph in targeted attack mode: (a) d = 5, d0 = 5 and (b)d = 7, d0 = 7.

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 Qik .
To evaluate validity of our analytical results showed
in Figure 7, we construct a graph with the number ofnodes set at 1000 and investigate the site percolation
process on it in targeted attack mode. The number oftrials 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 ineach method is slightly different, the order of robustness of each connection pattern is the same. From theseresults, analytical results are in good agreement withnumerical results and the graph in which the number ofboundary nodes is large has more robust connectivity

Toyonaga et al.

11

for small (1  p) values when we use the targeted attack
mode in which a removal node is selected in order ofdecreasing degree. After the fragmentation of the network, the network robustness is not significantly different for networks with different Mki (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 theresults 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 methodof 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 modeand Figure 9(b) shows the results in targeted attack
mode. Both show that difference between analyticaland numerical results becomes larger as the number ofnodes becomes smaller. This means that it is hard to
apply our proposal to a graph in which the number ofnodes is small because our approach derives average
properties of random graph ensembles. However, sucha 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 ofnodes 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 intomodules until all degree-7 nodes are removed. In the
numerical simulations, however, the graph can bedivided 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 thegraph. Therefore, the numerical results show a phasetransition 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 numberof degree-7 nodes is small and the graph fragments
easily.
The analytical results cannot completely capture anyof the numerical results in targeted attack mode becauseof the reasons described above. The result for 10,000
nodes is almost the same as for 20,000 nodes, which

Figure 9. Comparison of the analytical and numerical results:(a) random failure mode and (b) targeted attack mode.

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 thatfails in the numerical simulation, the way that nodes failis also an important factor which determines whether
analytic results completely capture numerical results ornot.

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.

12International Journal of Distributed Sensor Networks

Table 1. Parameters for constructing a topology within amodule.
ModelParameter

Value

ER
BA

p(ER)
m0(BA)
m(BA)
m(WS)
b(WS)

0.012
7
6
6
0.01

WS

ER model: Erdó´s–Rényi model; BA model: Barabási–Albert model; WS
model: Watts–Strogatz model.

Simulation settings. In this part, we investigate robustnessof a graph in which modules are constructed by Erdó´ s–
Rényi (ER) model, Barabási–Albert (BA) model, andWatts–Strogatz (WS) model,10 respectively. We use theparameters 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 amodule 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 oflinks within modules. We compare the numerical results
for the site percolation process on a graph constructedby modules and inter-module link, with the analyticalresults using the probability distribution obtained fromexisting modules.
Results using ER model. The numerical and analyticalresults 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 andthe 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, thegraph 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 analyticalresults in targeted attack mode are shown in Figure 11.
The results in Figure 11(b) are consistent with theresults shown above in terms of the first sudden decayof 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

Figure 10. The results in random failure mode (ER model):(a) numerical and (b) analytical results.

nodes are removed, and degree-6 nodes are crucial forthe connectivity of the intra-module because of theirhigh degrees. This leads to the vulnerable connectivityof intra-module network after the fragmentation. In
most of the results of Figure 11, after removal of theboundary nodes, the internal connectivity of module 1
is larger than that of module 2. Therefore, in suchcases, k is more dominant than k 0 after the fragmentation. Figure 11(a) shows that the order of robustness ofeach connection pattern is the same as the results ofanalysis although the giant component sizes are notcompletely the same.
Results using BA model. The numerical and analyticalresults 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 andthe maximum degree within a module is 76.

Toyonaga et al.

Figure 11. The results in targeted attack mode (ER model):(a) numerical and (b) analytical results.

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 graphwithin a module, the numerical and analytical results intargeted 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 arealmost the same, the graph in which the number ofboundary nodes is large has a slightly higher robust
connectivity before the fragmentation of the network
into two modules, in both analytical and numerical13

Figure 12. The results in random failure mode (BA model):(a) numerical and (b) analytical results.

simulations. After the fragmentation of the network,
the ranking of the sizes of the giant componentschanges. The reason for this is the same as that mentioned for the results of Figure 11.
Results using WS model. The numerical and analyticalresults 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 eachtopology is 100. The minimum degree within a moduleis 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 theresults of analysis. This is because the graph constructed by the WS model has a ring-shaped structure14International Journal of Distributed Sensor Networks

Figure 13. The results in targeted attack mode (BA model):(a) numerical and (b) analytical results.

Figure 14. The results in random failure mode (WS model):(a) numerical and (b) analytical results.

and can fragment easily in numerical simulations, while
the theory assumes random intra-module connections.
When we use the WS model for constructing a graphwithin a module, the numerical and analytical results intargeted attack mode are shown in Figure 15. From theresults shown in Figure 15(a) and (b), the results ofanalysis 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 theresults of analysis because of the same reason discussed
in random failure mode.
From these results, our proposal can capture theorder 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 whena 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 model4 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 therobustness 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 showedthat graphs have robust connectivity in terms of thepoint at which the network fragments into two modules

Toyonaga et al.

15(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 thepoint of fragmentation of the network into modules
when we fix the degree of the boundary nodes after thelink addition.
However, because we do not consider the structuralproperties of a graph within a module, a difference ofanalytical and numerical results occurs when a graphwithin a module has a special structure. It is our future
work how to resolve this problem.
For further investigation, we want to consider additionof 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 theconditional probability of the probability distribution oflinks 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 withrespect to the research, authorship, and/or publication of thisarticle.

Funding

Figure 15. The results in targeted attack mode (WS model):(a) numerical and (b) analytical results.

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 validityof the analysis results, we evaluate the percolation process on a graph constructed by a configuration modeland 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 thenumerical results, the results are in agreement qualitatively. Moreover, we investigate the applicable range ofour 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 agiven internal structure. The results show that our proposal can explain the order of the ranking of robustness

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

References
1. Khan I, Belqasmi F, Glitho R, et al. Wireless sensor network virtualization: a survey. IEEE Commun Surv Tut
2015; 18: 553–576.
2. Hong S, Lv C, Zhao T, et al. Cascading failure analysis
and restoration strategy in an interdependent network. J
Phys A: Math Theor 2016; 49(19): 195101, http://stacks.iop.org/1751-8121/49/i=19/a=195101
3. Toyonaga S, Kominami D and Murata M. Brain-inspired
method for constructing a robust virtual wireless sensor
network. In: Proceedings of the 10th international conference on computing and network communications (CoCoNet
2015), Thiruvananthapuram, India, 16–19 December
2015, pp.1–7. IEEE.
4. Melnik S, Porter MA, Mucha PJ, et al. Dynamics on
modular networks with heterogeneous correlations.
Chaos 2014; 24(2): 023106, http://www.ncbi.nlm.
nih.gov/pubmed/24985420
5. Newman ME. The structure and function of complex
networks. SIAM Rev 2003; 45(2): 167–256. http://dx.doi.org/10.1137/S003614450342480

16
6. Newman ME, Strogatz SH and Watts DJ. Randomgraphs with arbitrary degree distributions and theirapplications. Phys Rev E 2001; 64(2 Pt 2): 026118.
7. Leicht EA and D’Souza RM. Percolation on interacting
networks. ArXiv e-prints 2009; 2: 5, http://arxiv.org/
abs/0907.0894
8. Min B, Yi SD, Lee KM, et al. Network robustness ofmultiplex networks with interlayer degree correlations.
Phys Rev E 2014; 89(4): 1–9.

International Journal of Distributed Sensor Networks
9. Dong G, Gao J, Du R, et al. Robustness of network of
networks under targeted attack. Phys Rev E 2013; 87(5):
1–11.
10. Chopra M and Madan M. Network analysis by using
various models of the online social media networks. Int J
Adv Res Comput Sci 2015; 6(1): 111–116.