 Research
 Open Access
 Published:
Investigation on changes of modularity and robustness by edgeremoval mutations in signaling networks
BMC Systems Biology volume 11, Article number: 125 (2017)
Abstract
Background
Biological networks consisting of molecular components and interactions are represented by a graph model. There have been some studies based on that model to analyze a relationship between structural characteristics and dynamical behaviors in signaling network. However, little attention has been paid to changes of modularity and robustness in mutant networks.
Results
In this paper, we investigated the changes of modularity and robustness by edgeremoval mutations in three signaling networks. We first observed that both the modularity and robustness increased on average in the mutant network by the edgeremoval mutations. However, the modularity change was negatively correlated with the robustness change. This implies that it is unlikely that both the modularity and the robustness values simultaneously increase by the edgeremoval mutations. Another interesting finding is that the modularity change was positively correlated with the degree, the number of feedback loops, and the edge betweenness of the removed edges whereas the robustness change was negatively correlated with them. We note that these results were consistently observed in randomly structure networks. Additionally, we identified two groups of genes which are incident to the highlymodularityincreasing and the highlyrobustnessdecreasing edges with respect to the edgeremoval mutations, respectively, and observed that they are likely to be central by forming a connected component of a considerably large size. The geneontology enrichment of each of these gene groups was significantly different from the rest of genes. Finally, we showed that the highlyrobustnessdecreasing edges can be promising edgetic drugtargets, which validates the usefulness of our analysis.
Conclusions
Taken together, the analysis of changes of robustness and modularity against edgeremoval mutations can be useful to unravel novel dynamical characteristics underlying in signaling networks.
Background
Robustness and modularity are key properties to understand complex dynamics in largescale biological networks. The former means the capability of a network to maintain functioning against external and internal perturbations [1], and the latter describes the divisibility of a network into clusters [2]. The robust dynamics [3,4,5] and the modularized structures [6,7,8] have been ubiquitously observed through various biological examples. It is also notable that these properties can be changed by structural mutations because they are highly dependent on the network structure. For example, a few studies showed that the modularity is greatly changed by the removal of hubs [9] or by stabilizing events in protein–protein interaction networks. Some other studies also proved that the robustness is considerably changeable according to a variety of mutations [10,11,12,13]. Additionally, there were some previous studies to investigate a relation between the robustness and the modularity. For example, it was shown that the modularized structure of bone networks improves the robustness compared to a regular network of the same size [14]. Some other studies observed that both the robustness and the modularity characteristics could be emergently improved through a network evolution process [15, 16]. Moreover, there were some studies to explicitly examine linear correlations between the robustness and the modularity over differently structured networks [17,18,19]. In metabolic networks, the robustness against the mutant concentrations of metabolites or the mutant expression of enzymes has increased or decreased, respectively, as the modularity increases [17]. On the other hand, the robustness against a gene state perturbation was negatively correlated with the modularity in signaling networks [18, 19]. Although these previous studies found interesting relations between the robustness and the modularity, there are some issues needed to be investigated as follows. The first issue is that there is little known knowledge about changes of the modularity and the robustness. In particular, there was no intensive study about the relationship of the changes of the modularity and the robustness by structural mutations. We note that the previous studies [17,18,19] focused on the robustness and the modularity over networks with very different structures, whereas this study focuses on the changes of the robustness and the modularity over mutant networks with a slight structural modification. This means that the findings in the previous studies do not necessarily hold in our analysis. Another interesting issue is whether some wellknown motifs are relevant to the changes of the modularity and the robustness or not. In fact, some previous studies have shown that network motifs such as feedback loops (FBLs) and feedforward loops (FFLs) ubiquitously found in various biological networks can affect the robustness [13, 20]. For instance, it was reported that more positive and less negative FBLs are observed in robust networks [21]. Another study showed that coherent coupling of FBLs is a design principle of a robust signaling network [22]. It was also reported that coherent FFLs strengthen the robustness against updaterule perturbations [13]. To our best knowledge, even there was no reported motif which is relevant to the modularity property. Taken together, there is little known about motifs which indicate the changes of the modularity, the robustness, or both. The last issue is that there was no previous study to compare sets of nodes or interactions which efficiently control the changes of the modularity and the robustness. This can be impressive because the result can be used to identify functionally important nodes or interactions such as drug targets.
In this work, we tried to investigate the changes of the modularity and the robustness by edgeremoval mutations in signaling networks. Through intensive simulations using a Boolean network model [23, 24], we first found that both the modularity and the robustness increased on average against edgeremoval mutations, but the change of modularity is negatively correlated with the change of robustness. More intriguingly, the modularity change was positively correlated with the degree, the number of FBLs, and the edge betweenness of removed edges, whereas the robustness change was negatively correlated with them. Additionally, we found that these findings are consistently conserved in the random networks. Moreover, we identified two groups of genes which are incident to the highlymodularityincreasing and the highlyrobustnessdecreasing edges against the edgeremoval mutations, respectively, and observed that they are likely to be central by forming a considerably large connected component. The geneontology enrichment of each of the gene groups was clearly different from the rest of genes. Finally, we found that the highlyrobustnessdecreasing edges can be promising edgetic drugtargets. Taken together, the analysis of the changes of the robustness and the modularity against the edgeremoval mutations can be useful to reveal novel dynamical characteristics of signaling networks.
Methods
Network modularity
In this study, we examined the modularity by using the method in a previous study [25] and it has been widely used in many previous studies [18, 19, 26, 27]. Given a directed graph G(V, A) where V and A denote a set of nodes and a set of directed edges, respectively, we consider a partition P = {V _{1}, V _{2}, …, V _{ M }} of V(i.e. V _{ i } ∩ V _{ j } = ∅ for all i ≠ j, and \( \bigcup \limits_{i=1}^M{V}_i=V \)). The modularity of P is defined as \( M(P)=\sum \limits_{i=1}^M\left(\frac{\omega_{V_i{V}_i}}{\omega }\frac{\omega_{V_i}^{in}{\omega}_{V_i}^{out}}{\omega^2}\right) \), where \( {\omega}_{V_i{V}_i} \) is the number of directed edges whose both end nodes belong to V _{ i }, \( {\omega}_{V_i}^{out} \) and \( {\omega}_{V_i}^{in} \) are the numbers of directed edges whose starting or ending node only, respectively, belongs to V _{ i }, and ω is the total number of directed edges in the network. Then the network modularity can be defined as M(G) = max_{ P } M(P). Since it is difficult to optimize the partition, we computed the averaged modularity value of 30 trials of partitions optimized by an existing algorithm [28].
Boolean network dynamics
In this study, we employed a Boolean network model introduced in previous studies [29, 30] to investigate the complex dynamics of biological networks. In a Boolean network of a directed graph G(V, A), V and A denote a set of Boolean variables and a set of ordered pairs of the Boolean variables called directed edges, respectively. A state of each v _{ i } ∈ V is 1 or 0 which represents on or off state of the gene, respectively. Then a state of a network G is defined as a vector of the states of all nodes. A directed edge (v _{ j }, v _{ i }) ∈ A has a positive (activating) or negative (inhibiting) relationship from v _{ j } to v _{ i }. Here, we used a nested canalyzing function (NCF) model [31] (see Additional file 1: Supporting Text section for details), which can represent a variety of canalyzing rules in real molecular interactions [32] can be generated by using the NCF model. Additionally, NCFs properly fit the experimental data gained from literature [31], and can also describe logical interaction rules extracted from gene expression experiments [32, 33]. In this study, each NCF is randomly generated by specifying all I _{ m } s and O _{ m } s between 0 and 1 uniformly at random.
Let G(V, A) a Boolean network with a list of updaterules F = {f _{1}, f _{2}, …, f _{ N }}. Every initial state converges to an attractor which can describe diverse network dynamics such as multistability, homeostasis, and oscillation [34, 35]. Let α(s, G, F) the attractor which the initial state s converged. The network is considered as robust against a perturbation at v _{ i } if the attractor is conserved and we herein considered an updaterule mutation which describes a scenario that F is changed to \( {F}_i^{\prime }=\left\{{f}_1,\dots, {f}_i^{\prime },\dots, {f}_N\right\} \), where \( {f}_i^{\prime } \) means that every canalyzing and canalyzed values were flipped (i.e., all I _{ m } and O _{ m } are changed into 1 − I _{ m } and 1 − O _{ m }, respectively). This updaterule mutation may represent a deleterious change in the function of a protein or gene [36], and have been used in a previous study [13]. Then the network robustness γ(G) is defined as follows:
where S is a set of initial states (i.e. S = 2^{N}), and I(∙) is a function which outputs 1 or 0 if the condition is met or not, respectively. Because ∣S∣ is a very large number, we used a sample subset \( \overset{\sim }{S}\subseteq S \) with \( \left\overset{\sim }{S}\right=2N \) instead of S to calculate γ(G).
Changes of modularity and robustness by edgeremoval mutations
This study focuses on how the modularity and the robustness of a network are changed by edgeremoval mutations. Let m _{1} and r _{1} be the modularity and the robustness of the wildtype network, respectively. Given a removal rate parameter n (%), the mutant network is constructed by simultaneously removing approximately n percent of a total number of edges from the wildtype network. Then let m _{2} and r _{2} be the modularity and the robustness of the mutant network. We defined the changes of the modularity and the robustness by the edgeremoval mutations as (m _{2} − m _{1}) and (r _{2} − r _{1}), respectively. An illustrative example of the notion about the changes of modularity and robustness by edgeremoval mutations is shown in Fig. 1.
Signaling network datasets
To investigate real signaling networks, we used three datasets of signaling networks: a TLGL survival network (TLGL) [37] consisting of 60 genes and 142 interactions, a signal transduction network in fibroblasts (STF) [38] consisting of 139 genes and 557 interactions, and a HIV1 interaction network in Tcell (HIV1) [39] consisting of 138 genes and 368 interactions collected by manually curating signaling pathways from cellcollective (www.cellcollective.org) [40].
Generation of interactionshuffled random networks
We need to extensively simulate randomly structure networks to verify that the new findings in real networks are generally conserved. In this study, we employed a shuffling model to generate random networks [10, 18]. Given a reference network, it rewires some edges in a way that indegree and outdegree of every node are conserved. Accordingly, the structure of the generated random network is considerably similar to that of the original network.
Edgebased structural properties
A previous study has reported that there exists a relationship between a structural property with respect to genes or interactions and the global stability in biological networks [41]. In this regard, we investigated the relations of the following edgebased structural characteristics to the changes of the modularity and the robustness.

Degree of a node means the number of links incident to the node in a graph. On the other hand, the degree of an interaction (DEG) means the sum of the degrees of both end nodes of the edge.

An FBL is a circular chain where nodes are not revisited except the starting and the ending nodes [42]. It plays an important role in controlling the dynamical behaviors of cellular signaling networks. Specifically, v _{0} → v _{1} → v _{2} → … → v _{ L − 1} → v _{ L } is an FBL of length L (≥1) if there exist links from v _{ i − 1} to v _{ i } (i = 1, 2, …, L) with v _{0} = v _{ L } and v _{ j } ≠ v _{ k } for j, k ∈ {0, 1, …, L − 1} and j ≠ k. The number of FBLs of a link e denoted by NuFBL(e) means the number of different FBLs involving e.

Edge Betweenness (EBEW) is defined as the number of shortest paths between pairs of nodes that run along an edge [2], similar to Betweenness of a node. EBEW has been used as an important edgebased centrality measure in a previous study [43].
Software for statistical tests
In this study, IBM SPSS statistics [44] was used to conduct all statistical tests.
Results and Discussion
Relationship between changes of modularity and robustness by edgeremoval mutations
We first investigated the changes of the modularity and the robustness by edgeremoval mutations in three real networks TLGL, STF, and HIV1 (see Methods), and the results are shown in Fig. 2. (TLGL) and Fig. S1S2 (STF and HIV1, respectively) in Additional file 1. In this study, we computed the average changes of the modularity and the robustness values over 5000 trials of edgeremoval mutations. In addition, we varied the removal rate, which denotes the percentage of the number of removed edges over the total number of edges, from 1% to 5%. We first tested whether the average changes are significantly positive using onesample ttest. We note that the average changes were normally distributed, as assessed by KolmogorovSmirnov’s test (see Fig. S3S5 in Additional file 1 for details) and there were no or very few significant outliers, as assessed by a boxplot inspection (see Fig. S6S8 in Additional file 1 for details). As shown in Fig. 2(a), we observed that both average changes were positive for all removal rates, which means that the modularity and the robustness values were increased by edgeremoval (All Pvalues <0.0001; see Additional file 1: Figure S1 (a) and Figure S2 (a) for the results of STF and HIV1 networks, respectively).
In addition, the increase of the robustness was positively related to the removal rate. To examine the relationship between the changes of modularity and robustness values, we scattered them in the cases that the removal rate is 1% (Fig. 2(b)) and 2% (Fig. 2(c)). Intriguingly, there was a negative correlation between the modularity change and the robustness change, and this was consistently observed in the cases of larger removal rates (see Additional file 1: Figure S9) and the other networks (see Additional file 1: Figure S1 (b)(c)) and Fig. S2 (b)(c)). Figure 2 (d) shows the trend of the correlation coefficient values between the changes of modularity and robustness values against the removal rate, and we observed that they were significantly negative irrespective of the removal rates (see Additional file 1: Figure S1 (d) and S2 (d) for the results of STF and HIV1 networks, respectively). Actually, the negative relationship between the modularity and the robustness in signaling networks was observed in our previous studies [18, 19]. However, it should be noted that the previous finding does not imply any relation between the changes of the modularity and the robustness by edgeremoval mutations. To further examine if the negative relationship we found is a general property in randomly structured networks, we generated three sets of 100 random networks shuffled from TLGL, STF, and HIV1 (see Methods), and could observe consistent (see Additional file 1: Figure S10). This implies that such the negative relation between the changes of the modularity and the robustness can be regarded as a general principle conserved in randomly structured networks.
Structural characteristics to affect the changes of the modularity and the robustness
We showed that the changes of the modularity and the robustness are correlated when a network is subject to edgeremoval mutations. To reveal structural characteristics to affect the changes of the modularity and the robustness, we investigated the correlations of each of the changes of the modularity and the robustness with each of three edgebased structural properties, DEG, NuFBL and EBEW (see Methods for the definitions) in TLGL signaling network (Fig. 3; see Additional file 1: Figure S11 and S12 for the results of STF and HIV1, respectively). In Fig. 3, average DEG, NuFBL, and EBEW values of the removed edges over 5000 trials with 1% of the removal rate were examined. Intriguingly, we found that the change of the modularity is positively correlated with the average DEG, EBEW and NuFBL of removed edges (The correlation coefficients in Fig. 3 (a)(c) were 0.24708, 0.13786, and 0.11720, respectively, with all pvalues < 0.001). That is to say, removing edges with a higher degree, EBEW, or NuFBL is more likely to increase the network modularity. These results can be relevant to previous results. For example, the edges with high betweenness values are most likely to lie between subgraphs [45], and thus removing those edges could make a network more separately or more modularized. We also found that the change of the robustness is negatively correlated with the average DEG, EBEW and NuFBL of the removed edges (The correlation coefficients in Fig. 3(d)(f) were −0.21738, −0.14694, and −0.10537, respectively, with all pvalues < 0.0001). In other words, removing edges with a higher degree, EBEW, or NuFBL is more likely to decrease the network robustness. These results can be compared with some previously known results regarding nodebased mutations. For example, some studies reported that a node involving more FBLs is likely to be sensitive against nodebased mutations. To show that these results hold in random networks, we generated three sets of 100 random Boolean networks each of which was shuffled from TLGL, STF, and HIV1 networks, respectively. Through extensive simulations with the removal rate of 1%, we could observe consistent results (see Additional file 1: Fig. S13S15, All Pvalues < 0.0001 using ttest). In other words, the degree, the edgebetweenness and the number of FBLs were positively correlated with the change of the modularity whereas they were negatively correlated in the random networks. It means that those structural characteristics might be a vital factor in controlling both the changes of the modularity and the robustness.
Topological distribution of highly modularityincreasing and robustnessdecreasing edges by removal mutations
In the previous subsection, it was shown that the change of the modularity is positively correlated with the degree, the edge betweenness, and the number of involved FBLs with respect to the removed edges whereas the change of robustness is negatively correlated with them. From these results, we hypothesized that the edges whose removal will increase the modularity or decrease the robustness tend to be centrally located in signaling networks. To validate this hypothesis, we first specified “Highlymodularityincreasing” (HighMI) and “Highlyrobustnessdecreasing” (HighRD) sets of edges as follows: We examined the changes of the modularity and the robustness over 5000 trials of edgeremoval mutations with 1% removal rate, and collected topK set of edges among them in an increasing (resp. decreasing) order of the change of the modularity (resp. the robustness). Considering the distributions of the change of the modularity (resp. robustness), K was chosen to 20, 20, and 18 (resp., 31, 18, and 16) for TLGL, STF, and HIV1 networks, respectively. Then HighMI (resp. HighRD) denotes the union of the edges each of which was included in the modularityincreasing (resp. robustnessdecreasing) topK edges. Accordingly, we identified HighMI (HighRD) groups consisting of 22, 79, and 42 edges (resp. 30, 69, and 33 edges) in TLGL, STF, and HIV1 networks, respectively. Furthermore, we defined HighMIincident (HighRDincident) group which is a set of genes incident to an edge in the HighMI (resp. HighRD) edge group, and found the number of genes in the HighMIincident (resp. HighRDincident) were 29, 81, and 59 (resp. 33, 72, and 48) in TLGL, STF, and HIV1 networks, respectively. The topological distributions of HighMI and HighRD edge sets in TLGL, STF, and HIV1 networks are shown in Fig. 4 and Figure S16S17 in Additional file 1, respectively. As expected, it was observed that the edges in HighMI and HighRD groups are likely to be located at the centre of the signaling network. In order to more clarify this observation, we compared nodebased centrality values between each set of HighMIincident and HighRDincident groups and the set of rest genes. Specifically, we computed average degree, nodebased betweenness [43], stress [46], closeness [47], and the number of involved FBLs [22] for each group of nodes (Fig. 5). As depicted in the figure, we found that genes of HighMIincident and HighRDincident groups showed higher degree, nodebased betweenness, stress, closeness, and the number of FBLs than the rest of genes (Only three cases among 30 comparisons did not show significant differences.) In other words, the genes incident to the interactions whose greatly increase the modularity or decrease the robustness tends to be central in the signaling network. Additionally, we visualized the connectedness of edges of HighMI and HighRD groups by projecting them into a subnetwork from TLGL (see Fig. 4(c) and (d), respectively), STF (see Additional file 1: Fig. S16(c) and (d), respectively), and HIV1 (see Additional file 1: Figure S17(c) and (d), respectively) networks. As shown in the figures, every subnetwork forms a single connected component. This implies that the highly modularityincreasing or robustnessdecreasing edges with respect to edgeremoval mutations are closely located in signaling networks.
Gene ontology analysis of a set of genes incident to highlymodularityincreasing or highlyrobustnessdecreasing edges
We conducted Gene Ontotlogy (GO) enrichment analysis (The Gene Ontotlogy Consortium, 2008) using ClueGO tool [48] to investigate the locational and functional characteristics of sets of HighMIincident and HighRDincident genes. The results are shown in Table 1 (see Additional file 1: Table S1S2). Some GO terms such as protein tyrosine kinase and peptidase activity are more highly observed in HighMIincident and HighRDincident groups. The former is an enzyme which transfers a phosphate group from adenosine triphosphate to a protein in a cell, and the latter is catalysis of the hydrolysis of a peptide bond. In addition, HighMIincident and HighRDincident gene groups showed a greater fraction of response function terms. Regulation of adaptive immune response is any process that modulates the frequency, rate, or extent of an adaptive immune response regarding to robustness change. Furthermore, HighMIincident group showed a greater portion of vital binding functions. For example, a protein phosphatase is an enzyme that removes a phosphate group from the phosphorylated amino acid residue of its substrate protein, and its binding function is interacting selectively and noncovalently with any protein phosphatase. On the other hand, HighRDincident group showed a greater fraction related to signaling pathway. For instance, necroptosis is a programmed form of necrosis, or inflammatory cell death, and its signaling pathway is a series of molecular signals which triggers the necroptotic death of a cell. Taken together, significantly different functions between HighMIincident/HighRDincident groups of genes and the rest of genes can be characterized.
Edgebased drug discovery
We performed a case study to show an application for edgebased drug discovery. For every interaction in HighRD group, we examined the inclusion frequency of the interaction in topK edge sets ranked by a decreasing order of the robustness change among 5000 trials of edgeremoval mutations with 1% removal rate. We found that (JAK → STAT3), (IP3R1 → Ca), and (gp41 → CD28) showed the highest frequency in the TLGL, STF, and HIV1 networks, respectively. We hypothesized that these edges can be candidates of edgetic drugtargets, because they most frequently caused the highest decreasing robustness through removal mutations. To validate this, we surveyed some recent experimental studies. Regarding (JAK → STAT3) interaction of TLGL network (Fig. 6), it was shown that the interaction is associated with oncogenesis, proliferation, survival, metastasis, angiogenesis, and immune evasion in gastrointestinal cancers [49, 50]. For example, a colorectal cancer might be developed by dysregulation of the interleukin (IL)6mediated JAK → STAT3 pathway, and therefore strategies targeting the IL6/JAK/STAT3 pathway have emerged as attractive options to treat colorectal cancer [51]. Next, the (IP3R1 → Ca) interaction of STF network (see Additional file 1: Fig. S18) played an important role of dynamical relationship between IP3R1 and PI3K, which are the most influential components associated with drug resistance [52]. Systemic analysis of these components and their upstream components has resulted in identifying novel combinations of drug targets. In HIV1 network, (gp41 → CD28) was found to be the highest frequency interaction (see Additional file 1: Figure S19), but there was no relevant experimental study to support it. However, we could find biological evidence related to the second highest frequency interaction, (PI3K → PIP3). It is included in PI3K/Akt/mTOR pathway [53], which is known to be frequently activated in ovarian cancer. Therefore, inhibitors targeting this pathway can be evaluated as treatment strategies for ovarian cancer, either monotherapy or in combination with cytotoxic agents [54]. Another interesting point was the feedback loops involved with those interactions. We found a large number of feedback loops were related to (JAK → STAT3), (IP3R1 → Ca), and (PI3K → PIP3) in TLGL, STF, HIV1, respectively (The numbers were 70, 286, and 872, respectively). Considering that the number of involved FBLs was shown to be associated with the functional importance of a node or an interaction, it implies that the found interactions can be promising drugtargets.
Conclusions
There have been many computational studies about the network robustness and modularity, whereas there are few studies on investigating the modularity change and the robustness change. Through extensive simulations, we found that both the modularity and the robustness increased on average in mutant networks by edgeremoval mutations in this study. However, it was interesting that the changes of the modularity and the robustness were negatively correlated. Another interesting finding is that the changes of the modularity and the robustness are positively and negatively, respectively, correlated with each of the degree, the number of FBLs, and the edge betweenness of removed edges. These results were consistently observed in randomly structure networks. Additionally, we identified two sets of genes which are incident to the highlymodularityincreasing and the highlyrobustnessdecreasing edges, respectively, and observed that they are likely to be central by forming a large connected component. These two gene sets were enriched with different GO terms and the investigation on the reason why such GO terms are related to modularity and robustness will be a future study. Finally, we found that the highlyrobustnessdecreasing edge can be considered for promising edgebased drugtargets. Taken together, our results in this study can be useful to unravel novel dynamical characteristics of signaling networks.
Abbreviations
 DEG:

Degree
 EBEW:

Edge Betweenness
 FBLs:

Feedback loops
 FFLs:

Feedforward loops
 GO:

Gene ontology
References
 1.
Kitano H. Biological robustness. Nat Rev Genet. 2004;5(11):826–37.
 2.
Girvan M, Newman MEJ. Community structure in social and biological networks. Proc Natl Acad Sci. 2002;99(12):7821–6.
 3.
Ingolia NT. Topology and robustness in the drosophila segment polarity network. PLoS Biol. 2004;2(6):e123.
 4.
Yi TM, Huang Y, Simon MI, Doyle J. Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc Natl Acad Sci. 2000;97(9):4649–53.
 5.
Little JW, Shepley DP, Wert DW. Robustness of a gene regulatory circuit. EMBO J. 1999;18(15):4299–307.
 6.
Kreimer A, Borenstein E, Gophna U, Ruppin E. The evolution of modularity in bacterial metabolic networks. Proc Natl Acad Sci. 2008;105(19):6976–81.
 7.
Lin YS, Hsu WL, Hwang JK, Li WH. Proportion of solventexposed amino acids in a protein and rate of protein evolution. Mol Biol Evol. 2007;24(4):1005–11.
 8.
von Dassow G, Munro E. Modularity in animal development and evolution: elements of a conceptual framework for EvoDevo. J Exp Zool. 1999;285(4):307–25.
 9.
Han JDJ, Bertin N, Hao T, Goldberg DS, Berriz GF, Zhang LV, Dupuy D, Walhout AJM, Cusick ME, Roth FP, et al. Evidence for dynamically organized modularity in the yeast proteinprotein interaction network. Nature. 2004;430(6995):88–93.
 10.
Trinh HC, Kwon YK. Edgebased sensitivity analysis of signaling networks by using Boolean dynamics. Bioinformatics. 2016;32(17):i763–71.
 11.
Paroni A, Graudenzi A, Caravagna G, Damiani C, Mauri G, Antoniotti M. CABeRNET: a Cytoscape app for augmented Boolean models of gene regulatory NETworks. BMC Bioinformatics. 2016;17(1):64.
 12.
Kaneko K. Evolution of robustness to noise and mutation in gene expression dynamics. PLoS One. 2007;2(5):e434.
 13.
Le DH, Kwon YK. A coherent feedforward loop design principle to sustain robustness of biological networks. Bioinformatics. 2013;29(5):630–7.
 14.
Viana MP, Tanck E, Beletti ME. Costa LdF: modularity and robustness of bone networks. Mol BioSyst. 2009;5(3):255–61.
 15.
Variano EA, McCoy JH, Lipson H. Networks, dynamics, and modularity. Phys Rev Lett. 2004;92(18):188701.
 16.
Hintze A, Adami C. Evolution of complex modular biological networks. PLoS Comput Biol. 2008;4(2):e23.
 17.
Holme P. Metabolic robustness and network modularity: a model study. PLoS One. 2011;6(2):e16605.
 18.
Truong CD, Tran TD, Kwon YK. MORO: a Cytoscape app for relationship analysis between modularity and robustness in largescale biological networks. BMC Syst Biol. 2016;10(Suppl 4):122.
 19.
Tran TD, Kwon YK. The relationship between modularity and robustness in signalling networks. J R Soc Interface. 2013;10(88):20130771.
 20.
Kim JR, Yoon Y, Cho KH. Coupled feedback loops form dynamic motifs of cellular networks. Biophys J. 2008;94(2):359–65.
 21.
Kwon YK, Cho KH. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics. 2008;24(7):987–94.
 22.
Kwon YK, Cho KH. Coherent coupling of feedback loops: a design principle of cell signaling networks. Bioinformatics. 2008;24(17):1926–32.
 23.
Kauffman S. A proposal for using the ensemble approach to understand genetic regulatory networks. J Theor Biol. 2004;230(4):581–90.
 24.
Graudenzi A, Serra R, Villani M, Colacci A, Kauffman SA. Robustness analysis of a Boolean model of gene regulatory network with memory. J Comput Biol. 2011;18(4):559–77.
 25.
Leicht EA, Newman MEJ. Community structure in directed networks. Phys Rev Lett. 2008;100(11):118703.
 26.
Fortunato S. Community detection in graphs. Phys Rep. 2010;486(3):75–174.
 27.
Mucha PJ, Richardson T, Macon K, Porter MA, Onnela JP. Community structure in timedependent, multiscale, and multiplex networks. Science. 2010;328(5980):876.
 28.
Noack A. Modularity clustering is forcedirected layout. Phys Rev E. 2009;79(2):026102.
 29.
Campbell C, Albert R. Stabilization of perturbed Boolean network attractors through compensatory interactions. BMC Syst Biol. 2014;8(1):53.
 30.
Steinway SN, Biggs MB, Loughran TP, Papin JA, Albert R. Inference of network dynamics and metabolic interactions in the gut microbiome. PLoS Comput Biol. 2015;11(6):e1004338.
 31.
Kauffman S, Peterson C, Samuelsson B, Troein C. Random Boolean network models and the yeast transcriptional network. Proc Natl Acad Sci. 2003;100(25):14796–9.
 32.
Harris SE, Sawhill BK, Wuensche A, Kauffman S. A model of transcriptional regulatory networks based on biases in the observed regulation rules. Complexity. 2002;7(4):23–40.
 33.
Naldi A, Carneiro J, Chaouiya C, Thieffry D. Diversity and plasticity of Th cell types predicted from regulatory network Modelling. PLoS Comput Biol. 2010;6(9):e1000912.
 34.
Bhalla US, Ram PT, Iyengar R, Kinase Phosphatase MAP. As a locus of flexibility in a mitogenactivated protein kinase signaling network. Science. 2002;297(5583):1018.
 35.
Pomerening JR, Sontag ED, Ferrell JE. Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2. Nat Cell Biol. 2003;5(4):346–51.
 36.
Ng PC, Henikoff S. SIFT: predicting amino acid changes that affect protein function. Nucleic Acids Res. 2003;31(13):3812–4.
 37.
Saadatpour A, Wang RS, Liao A, Liu X, Loughran TP, Albert I, Albert R. Dynamical and structural analysis of a T cell survival network identifies novel candidate therapeutic targets for large granular lymphocyte leukemia. PLoS Comput Biol. 2011;7(11):e1002267.
 38.
Hirabayashi T, Murayama T, Shimizu T. Regulatory mechanism and physiological role of cytosolic phospholipase A_{2}. Biol Pharm Bull. 2004;27(8):1168–73.
 39.
Oyeyemi OJ, Davies O, Robertson DL, Schwartz JM. A logical model of HIV1 interactions with the Tcell activation signalling pathway. Bioinformatics. 2015;31(7):1075–83.
 40.
Helikar T, Kowal B, McClenathan S, Bruckner M, Rowley T, Madrahimov A, Wicks B, Shrestha M, Limbu K, Rogers JA. The cell collective: toward an open and collaborative approach to systems biology. BMC Syst Biol. 2012;6(1):96.
 41.
Kaiser M, Hilgetag CC. Edge vulnerability in neural and metabolic networks. Biol Cybern. 2004;90(5):311–7.
 42.
Ananthasubramaniam B, Herzel H. Positive feedback promotes oscillations in negative feedback loops. PLoS One. 2014;9(8):e104761.
 43.
Freeman LC, Set A. Of measures of centrality based on Betweenness. Sociometry. 1977;40(1):35–41.
 44.
IBM SPSS Statistics. https://www.ibm.com/products/spssstatistics. Accessed 11 Sept 2017.
 45.
Yoon J, Blumer A, Lee K. An algorithm for modularity analysis of directed and weighted biological networks based on edgebetweenness centrality. Bioinformatics. 2006;22(24):3106–8.
 46.
Shimbel A. Structural parameters of communication networks. The bulletin of mathematical biophysics. 1953;15(4):501–7.
 47.
Wuchty S, Stadler PF. Centers of complex networks. J Theor Biol. 2003;223(1):45–53.
 48.
Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman WH, Pagès F, Trajanoski Z, Galon J. ClueGO: a Cytoscape plugin to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3.
 49.
Nikolaou K, Sarris M, Talianidis I. Molecular pathways: the complex roles of inflammation pathways in the development and treatment of liver cancer. Clin Cancer Res. 2013;19(11):2810.
 50.
Bournazou E, Bromberg J. Targeting the tumor microenvironment: JAKSTAT3 signaling. JAKSTAT. 2013;2(2):e23828.
 51.
Wang SW, Sun YM. The IL6/JAK/STAT3 pathway: potential therapeutic strategies in treating colorectal cancer (review). Int J Oncol. 2014;44(4):1032–40.
 52.
Puniya BL, Allen L, Hochfelder C, Majumder M, Helikar T. Systems perturbation analysis of a largescale signal transduction model reveals potentially influential candidates for cancer therapeutics. Frontiers in Bioengineering and Biotechnology. 2016;4:10.
 53.
Slomovitz BM, Coleman RL. The PI3K/AKT/mTOR pathway as a therapeutic target in endometrial cancer. Clin Cancer Res. 2012;18(21):5856–64.
 54.
Mabuchi S, Kuroda H, Takahashi R, Sasano T. The PI3K/AKT/mTOR pathway as a therapeutic target in ovarian cancer. Gynecol Oncol. 2015;137(1):173–9.
Acknowledgements
This work was supported by the 2017 Research Fund of University of Ulsan.
Funding
Publication charges for this work were funded by the 2017 Research Fund of University of Ulsan.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information files.
About this supplement
This article has been published as part of BMC Systems Biology Volume 11 Supplement 7, 2017: 16th International Conference on Bioinformatics (InCoB 2017): Systems Biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume11supplement6.
Author information
Affiliations
Contributions
YKK conceived the study, CDT performed simulations, CDT and YKK designed the analysis, CDT and YKK wrote and revised the manuscript. Both authors read and approved the final manuscript.
Corresponding author
Correspondence to YungKeun Kwon.
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1: Figure S1.
Analysis of the changes of the modularity and the robustness by edgeremoval mutations in STF signaling network. Figure S2. Analysis of the changes of the modularity and the robustness by edgeremoval mutations in HIV1 signaling network. Figure S3. Analysis of normal distributions of averages of modularity changes and robustness changes in TLGL network. Figure S4. Analysis of normal distributions of averages of modularity changes and robustness changes in STF network. Figure S5. Analysis of normal distributions of averages of modularity changes and robustness changes in HIV1 network. Figure S6. Analysis of outliers of averages of modularity changes and robustness changes in TLGL network. Figure S7. Analysis of outliers of averages of modularity changes and robustness changes in STF network. Figure S8. Analysis of outliers of averages of modularity changes and robustness changes in HIV1 network. Figure S9. Relationship between the changes of the modularity and the robustness in TLGL signaling network. Figure S10. Relationship between the changes of the modularity and the robustness in random networks. Figure S11. Relationship of each of the changes of the modularity and the robustness with the structural properties in STF signaling network. Figure S12. Relationship of each of the changes of the modularity and the robustness with the structural properties in HIV1 signaling network. Figure S13. Relationship of each of the changes of the modularity and the robustness with the structural properties in random networks shuffled from TLGL network. Figure S14. Relationship of each of the changes of the modularity and the robustness with the structural properties in random networks shuffled from STF network. Figure S15. Relationship of each of the changes of the modularity and the robustness with the structural properties in random networks shuffled from HIV1 network. Figure S16. Topological distributions of HighMI/HighRD edges and their incident nodes in STF signaling network. Figure S17. Topological distributions of HighMI/HighRD edges and their incident nodes in HIV1 signaling network. Figure S18. Edgeremoval analysis for edgetic drug discovery in STF signaling network. Figure S19. Edgeremoval analysis for edgetic drug discovery in HIV1 signaling network. Table S1. GO analysis results between HighMIincident/HighRDincident group and the rest of genes in STF network. Table S2. GO analysis results between HighMIincident/HighRDincident group and the rest of genes in HIV1 network. (PDF 3490 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Published
DOI
Keywords
 Boolean dynamics
 Edgeremoval mutations
 Robustness
 Modularity
 Feedback loops
 Centrality
 Geneontology
 Drugtargets