 Research
 Open access
 Published:
Effects of ordered mutations on dynamics in signaling networks
BMC Medical Genomics volumeÂ 13, ArticleÂ number:Â 13 (2020)
Abstract
Background
Many previous clinical studies have found that accumulated sequential mutations are statistically related to tumorigenesis. However, they are limited in fully elucidating the significance of the orderedmutation because they did not focus on the network dynamics. Therefore, there is a pressing need to investigate the dynamics characteristics induced by orderedmutations.
Methods
To quantify the orderedmutationinducing dynamics, we defined the mutationsensitivity and the orderspecificity that represent if the network is sensitive against a double knockout mutation and if mutationsensitivity is specific to the mutation order, respectively, using a Boolean network model.
Results
Through intensive investigations, we found that a signaling network is more sensitive when a doublemutation occurs in the direction order inducing a longer path and a smaller number of paths than in the reverse order. In addition, feedback loops involving a gene pair decreased both the mutationsensitivity and the orderspecificity. Next, we investigated relationships of functionally important genes with orderedmutationinducing dynamics. The network is more sensitive to mutations subject to drugtargets, whereas it is less specific to the mutation order. Both the sensitivity and specificity are increased when differentdrugtargeted genes are mutated. Further, we found that tumor suppressors can efficiently suppress the amplification of oncogenes when the former are mutated earlier than the latter.
Conclusion
Taken together, our results help to understand the importance of the order of mutations with respect to the dynamical effects in complex biological systems.
Background
In a tumor cell, DNA damage restoration can have some errors such as chromosome abnormalities or genetic instability that result in a sequence of mutations [1]. The accumulated mutations can cause tumorigenesis or cancer development [2, 3]. Interestingly, this process can be affected by the order of genes subject to mutations. For example, it was observed in patients with Myeloproliferative Neoplasms that JAK2 mutation followed by TET2 mutation influenced the clinical features [3]. In another study [4], it has been shown that the timing of the DNMT3A mutation can affect the phenotypes of myeloid diseases in different ways. It was also found that the mutation of the CSF3R gene arising in the early severe congenital neutropenia stage is crucial to leukemia transformation [5]. Moreover, the mutation order influences the mutagen target size in tumor evolution [6] and results in complications in cancer biology [7].
Based on these observations, many methods attempted to investigate the effects of the order of mutation occurrence on cancer development [8], differences in clinical presentation [3], or response to targeted therapy [9] through biological experiments. In addition, statistical approaches were also developed to estimate the effect of the mutation order [10,11,12]. The understanding of the significance of ordered mutations can be enhanced through the analysis of biological networks. For example, a few previous studies reported that specific sequences of mutations which are efficient in cancer development were related to genegene interaction networks [13,14,15]. In addition, simple network structural characteristics were found to be relevant to the ordered mutations [16, 17]. Other studies also implicitly attempted to relate the dependency of ordered compensatory perturbation with a growth rate of cells in metabolic systems, through the network structural analysis and biological experiments [18,19,20]. Although these previous studies explained the significance of the ordered mutations, the results are limited due to the deficiency of analysis on the network dynamics.
Accordingly, we aim to investigate interesting dynamical characteristics induced by the ordered doublemutations in the signaling networks. To this end, we employed a Boolean network model and defined two measures with respect to orderedmutationinducing dynamics, the mutationsensitivity and the orderspecificity, against a double knockout mutation. The former represents how likely a network state trajectory changes by an ordered mutation whereas the latter indicates the likelihood that the network state trajectories induced by different orders of a mutation are not identical. Through intensive investigations in three real signaling networks, we found that a network is more sensitive when a double knockout mutation occurs in the order to induce a longer path and a smaller number of paths than in the reverse order. In addition, the existence of a feedback loop structure reduced the mutationsensitivity as well as the orderspecificity. Next, we investigated the orderedmutationinducing dynamics of some functionally important genes such as drugtargets, tumor suppressors, and oncogenes. It was interesting that the number of drugtargets subject to mutations was negatively correlated to the mutationsensitivity, whereas the mutation order was more specific in mutations in drugtargets than nondrugtargets. In addition, both the mutationsensitivity and the orderspecificity of samedrugtargets were higher than those of the differentdrugtargets, respectively. Furthermore, we found that tumor suppressors can efficiently suppress the amplification of oncogenes when the former genes are mutated earlier than the latter genes. Taken together, our results enhance the understanding of the dynamical effects of ordered doublemutations in complex biological systems.
Methods
Datasets
In this work, we employed datasets of three molecular interaction networks, a human cancer signaling (HCS) network with 1192 genes and 3102 interactions constructed in previous studies [21, 22] to provide a map of human cancer signaling, another largescale signaling network with 1659 genes and 7964 interactions constructed in previous study [23] which was derived from the Kyoto Encyclopedia of Genes and Genomes database [24] (KEGG) network, and a Tcell large granular lymphocyte survival signaling [16, 25, 26] (TGL) network with 61 genes and 193 interactions about the longterm survival of competent cytotoxic T lymphocytes in humans. Moreover, we retrieved lists of drugtargets, tumor suppressors, and oncogenes from the DrugBank [27], TSGene [28, 29], and ONGene [30] databases, respectively. Accordingly, we found 504, 538, and 20 drugtargets in the HCS, KEGG, and TGL networks, respectively. In addition, we identified 245 tumorsuppressors and 227 oncogenes in HCS, 176 tumorsuppressors and 168 oncogenes in KEGG, and 6 tumorsuppressors and 13 oncogenes in TGL networks (see Additional file 1: Table S1, S2 and S3 for the lists of drugtargets, tumor suppressor genes, and oncogenes in HCS, KEGG, and TGL networks, respectively).
A Boolean network model
To analyze the network dynamics induced by ordered mutations, we applied a Boolean network model, which is the simplest computational model [31,32,33] and has been used to examine complex behaviors of biological networks [34, 35]. A Boolean network is represented by a directed graph G(V,â€‰A) where Vâ€‰=â€‰{v_{1},â€‰v_{2},â€‰â€¦,â€‰v_{N}} is a set of nodes and Aâ€‰âŠ†â€‰Vâ€‰Ã—â€‰V is a set of directed links. Each v_{i}â€‰âˆˆâ€‰V has a value of 1 (on) or 0 (off), which indicates the possible states of the corresponding elements. A directed link (v_{i},â€‰v_{j}) represents a positive (activating) and a negative (inhibiting) relationship from v_{i} to v_{j}. Let v(t) denote the state of node v at timestep t (t is a nonnegative integer). When a state of v_{i} at time tâ€‰+â€‰1 is determined by the values of k_{i} and other nodes \( {v}_{i_1},{v}_{i_2},\dots, {v}_{i_{k_i}} \) with a link to v_{i} at time t, then the update rule of v_{i} is represented by a Boolean function \( {f}_i:{\left\{0,1\right\}}^{k_i}\to \left\{0,1\right\} \). Herein, all nodes are synchronously updated, and we employed a nested canalyzing function (NCF) model [36, 37] to describe an update rule f_{i} as follows:
where I_{m} and O_{m} (mâ€‰=â€‰1,â€‰2,â€‰â€¦,â€‰k_{i}) represent the canalyzing and canalyzed Boolean values, respectively, and O_{def} is set to \( 1{O}_{k_i} \) in general. Unfortunately, it is not easy to infer the canalyzing and the canalyzed values in the real signaling networks, so we specified I_{m} and O_{m} values independently and uniformly at random between 0 and 1. We note that NCFs have been shown to properly fit real biological experimental data [36, 38], and many biological networks were successfully simulated by NCFs [39, 40].
A network state at time t can be denoted by a list of state values of all nodes, v(t)â€‰=â€‰[v_{1}(t),â€‰v_{2}(t),â€‰â€¦,â€‰v_{N}(t)]â€‰âˆˆâ€‰{0,â€‰1}^{N}. Every network state transits to another network state determined by a set of Boolean update functions Fâ€‰=â€‰{f_{1},â€‰f_{2},â€‰â€¦,â€‰f_{N}} which is synchronously updated f_{1}, f_{2}, â€¦f_{N}, and f_{1}, f_{2}, â€¦f_{N} eventually converges to either a fixed point or a limitcycle attractor. The attractor is rigorously defined as follows.
Definition. Let v(0), v(1), â‹¯v(t), â‹¯ be a network state trajectory starting at v(0). Then, the attractor denoted by ã€ˆG,â€‰F,â€‰v(0)ã€‰ is represented by an ordered finite sublist of the trajectory, [v(Ï„),â€‰v(Ï„â€‰+â€‰1),â€‰â€¦,â€‰v(Ï„â€‰+â€‰pâ€‰âˆ’â€‰1)], where Ï„ is the smallest timestep such that v(t)â€‰=â€‰v(tâ€‰+â€‰p) for âˆ€tâ€‰â‰¥â€‰Ï„ with v(i)â€‰â‰ â€‰v(j) for âˆ€â€‰iâ€‰â‰ â€‰jâ€‰âˆˆâ€‰{Ï„,â€‰Ï„â€‰+â€‰1,â€‰â€¦,â€‰Ï„â€‰+â€‰pâ€‰âˆ’â€‰1}. Herein, p is called the attractor length.
To identify an attractor, the network state trajectory is computed by synchronously updating the state values of all nodes until the timestep t is found such that v(t)â€‰=â€‰v(tâ€‰+â€‰p).
Computation of mutationsensitivity and orderspecificity
Given a Boolean network G(V,â€‰A) with a set of nodes Vâ€‰=â€‰{v_{1},â€‰v_{2},â€‰â€¦,â€‰v_{N}} specified by a set of corresponding updaterules Fâ€‰=â€‰{f_{1},â€‰f_{2},â€‰â€¦,â€‰f_{N}}, consider a state trajectory starting from an arbitrary initial state. When a network is subject to a mutation, the trajectory may converge to a different attractor. Then, the network is regarded as sensitive to the mutation. Let Wâ€‰âŠ†â€‰V be a set of genes subject to knockout mutations [41, 42], and we denote by F^{W} as a set of updaterules where every gene in Wâ€‰âŠ†â€‰V is frozen to 0 (off state) in F. In this work, we investigate the effect of the ordered doublemutations on the network dynamics. Let (v_{k},â€‰v_{l}) be an ordered pair of nodes subject to a doublemutation with a time gap T, which means that v_{k} is first mutated at timestep tâ€‰=â€‰0, and then v_{l} is mutated at timestep tâ€‰=â€‰T. In other words, the time gap represents the timestep lag between the occurrences of the first and the second mutation. We can implement it by assuming that \( {F}^{\left\{{v}_k\right\}} \) and \( {F}^{\left\{{v}_k,{v}_l\right\}} \) are effective for 0â€‰â‰¤â€‰tâ€‰<â€‰T and tâ€‰â‰¥â€‰T, respectively. It has been known that the notion of the time gap is important since it can affect the mutation process [14, 43, 44]. Then, when we denote by \( {F}_{\left({v}_k,{v}_l\right)}^{\prime } \) a series of sets of the update rules by (v_{k},â€‰v_{l})ordered doublemutation, we can define the mutationsensitivity as follows:
where S is a set of considered initialstates and I(condition) denotes an indicator function that returns 1 if the condition is true and 0, otherwise. In other words, Î´ represents the probability that a network converges to a different attractor by the double knockout mutation. To quantify the specificity of dynamics with respect to the mutation order, we define the orderspecificity as follows:
In other words, Î” represents the probability of a network converging to different attractors by different mutation orders. Figure 1 shows an illustrative example of the mutationsensitivity and the orderspecificity notions. Let {v_{3},â€‰v_{4}} be a pair of genes subject to mutations and 0100 be a given initial state. Then, a wildtype attractor denoted by Att1 is computed by applying F all the time. Next, we compute Att2 and Att3 to which the network converges against (v_{3},â€‰v_{4}) and (v_{4},â€‰v_{3})ordered mutations, respectively. Note that \( {F}^{\left\{{v}_3\right\}} \) (or \( {F}^{\left\{{v}_4\right\}} \)) and \( {F}^{\left\{{v}_3,{v}_4\right\}} \) applies for 0â€‰â‰¤â€‰tâ€‰<â€‰T and tâ€‰â‰¥â€‰T, respectively, in computing Att2 (resp. Att3). When Att2 (or Att3) is not identical to Att1, the network is regarded as sensitive to the (v_{3},â€‰v_{4})ordered (resp. (v_{4},â€‰v_{3})ordered) mutation. In addition, the network dynamics are specific to the mutation order if Att2 and Att3 are not identical to each other. We note that a pair of genes (v_{k},â€‰v_{l}) with a common child node in the network are excluded from analysis in this study. In case that there exists such a common child node v_{c}, it is probable that the update of v_{c} is differently affected by the (v_{k},â€‰v_{l})ordered and the (v_{l},â€‰v_{k})ordered mutations, because the occurrence order of v_{k} and v_{l} in the NCF to update v_{c} represents the degree of influence on the update of v_{c}. Finally, we note that the network dynamics can depend on the initial network states. Therefore, a total of 1000 initialstates (i.e., Sâ€‰=â€‰1000) were randomly generated to compute the mutationsensitivity and the orderspecificity values in Eqs. (1) and (2) in all simulations of this study.
Structural characteristics of ordered gene pairs
Some structural characteristics of genes are related to network dynamic stability [45]. In this study, we considered the following structural properties to investigate the relationship between ordered mutations.

The path length of (v_{i},â€‰v_{j}) denoted by l(v_{i},â€‰v_{j}) is defined by the number of links included in the shortest path from v_{i} and v_{j}.

The number of paths of (v_{i},â€‰v_{j}) denoted by n(v_{i},â€‰v_{j}) is defined by the number of nonidentical paths from v_{i} and v_{j}.

A feedback loop (FBL) is a sequence of nodes where no node is repeated except the starting and the ending nodes. It has been known that FBLs play an important role in controlling dynamics behaviors of cellular signaling networks [32, 46, 47].
Random network generation
To verify that the results of mutationsensitivity and orderspecificity in the real molecular interaction networks are consistent with randomly structured networks, we generated random networks using the BarabÃ¡si Albert (BA) [48] model which is a kind of network growth model with a preferential attachment scheme.
Parallel computation
For efficient insilico simulations, we basically implemented the program code using PANET [35] which is an analysis tool of the network dynamics analysis tool using the OpenCL library. This enables us to compute a large number of attractors in parallel by assigning each initial random state in Eqs. (1) and (2) to a processing unit of CPUs and/or GPUs.
Statistical analysis
In this paper, we conducted the MannWhitney U test to see if the mutationsensitivity and the orderspecificity are significantly different between any two groups, because they are not normally distributed. The MannWhitney U test combines two groups and ranks them. Then, it calculates a statistic of the difference of the rank sum between two resampled groups. We used MedCalc [49] Statistical Software (version 13.0.6) for the MannWhitney U test.
Results
In this study, we simulated the orderedmutationinducing dynamics of three real biological networks using the Boolean network model (see Methods for details). As explained in Methods section, we note that a total of 1000 initialstates were randomly generated to compute the mutationsensitivity and the orderspecificity values in Eqs. (1) and (2) in all simulations. In addition, we constructed a set of ordered gene pairs to be investigated, Î©, for tractable simulation. It consists of all ordered pairs of genes in the case of TGL network with a small number of nodes (Nâ€‰=â€‰61), whereas 30,000 randomly selected gene pairs in the case of largescale networks of HCS (Nâ€‰=â€‰1192) and KEGG (Nâ€‰=â€‰1659). Considering the different network size, we also set the time gap (T) to 2â€“20 in HCS and KEGG, and 1â€“10 in TGL, respectively.
Distributions of orderedmutationinducing dynamics
To see how frequently the network dynamics are affected by ordered mutations, we examined the accumulative distributions of the mutationsensitivity (Î´) and orderspecificity (Î”) values of examined gene pairs (Î©) in three signaling networks (Fig. 2) in the case of the largest time gap (i.e., Tâ€‰=â€‰20 for HCS and KEGG, and Tâ€‰=â€‰10 for TGL). In the figure, the yaxis value means the cumulative probability of mutationsensitivity or orderspecificity larger or equal to the xaxis value. We observed that the cumulative probabilities of Î´â€‰â‰¥â€‰0.1 in HCS, KEGG, and TGL were 0.56, 0.62, and 0.39, respectively. This implies that it is not rare to observe that the network dynamics are sensitive against the doublemutations. It was also observed that the cumulative probabilities of Î”â€‰â‰¥â€‰0.1 in HCS, KEGG, and TGL were 0.38, 0.54, and 0.32, respectively. We need to note that an orderspecificity of zero can be observed even in the case of gene pairs with nonzero mutationsensitive values according to the definitions. Therefore, the observed distribution of the orderspecificity implies that the mutation order is considerably critical to the network dynamics.
Relation between structural characteristics and orderedmutationinducing dynamics
There have been many previous studies on the relationship between the structural properties and the dynamical behavior in biological networks [45, 50, 51]. Inspired by them, we investigated the relationships between some structural properties and orderedmutationinducing dynamics (Fig. 3). We first classified every ordered gene pair (v_{i},â€‰v_{j}) in Î© into â€˜Shorterpath directionâ€™ and â€˜Longerpath directionâ€™ groups if l(v_{i},â€‰v_{j})â€‰<â€‰l(v_{j},â€‰v_{i}) and l(v_{i},â€‰v_{j})â€‰>â€‰l(v_{j},â€‰v_{i}) (see Methods for the definition), respectively, and forced knockout mutations in the order of v_{i} and v_{j}. We note that the gene pair (v_{i},â€‰v_{j}) which is not bidirectionally connected was excluded from analysis to remove the effect of the connectedness factor on the dynamics. We compared the average mutationsensitivity values between them (Fig. 3(a)(c); all Pvalues using the MannWhitney U test). As shown in the figure, the mutationsensitivity of the â€˜Longerpath directionâ€™ group is significantly higher than that of the â€˜Shorterpath directionâ€™ group in all signaling networks for most time gap parameter values. In other words, the network is more sensitive when the double knockout mutation occurs in the order inducing a longer path than in the reverse order. Next, we classified every ordered gene pair (v_{i},â€‰v_{j}) into â€˜Morepaths directionâ€™ and â€˜Fewerpaths directionâ€™ groups if n(v_{i},â€‰v_{j})â€‰>â€‰n(v_{j},â€‰v_{i}) and n(v_{i},â€‰v_{j})â€‰<â€‰n(v_{j},â€‰v_{i}) (see Methods for the definition), respectively, and forced knockout mutations in the order of v_{i} and v_{j}. We compared the average mutationsensitivity values between them (Fig. 3(d)(f); all Pvalues using the MannWhitney U test). As shown in the figure, the mutationsensitivity of the former group is significantly smaller than that of the latter group in both signaling networks, almost irrespective of the time gap parameter. In other words, the network is more sensitive when the double knockout mutation occurs in the order involving fewer paths than in the reverse order. We note that our previous study showed that the dynamics influence from a gene on another gene is likely to be lessened as the path length increases and the number of paths decreases [50]. Thus, it is interesting that both the â€˜Longerpath directionâ€™ and â€˜Fewerpaths directionâ€™, which showed relatively higher mutationsensitivity values, represent ways to induce a smaller dynamicsinfluence from the first mutated gene on the second mutated gene than the reverse order. Finally, we considered the FBL as another interesting structural property for investigation, because many previous studies have proven the relation of it with the dynamical behavior of biological networks [46, 47, 52]. We classified every ordered gene pair into â€˜FBLâ€™ and â€˜NonFBLâ€™ groups if any gene in the pair is involved in an FBL or not, respectively. Then, we compared the average mutationsensitivity values between them (Fig. 3(g)(i); all Pvalues using the MannWhitney U test). As shown in the figure, the mutationsensitivity of the former group is significantly smaller than that of the latter group in both signaling networks regardless of the time gap parameter. In addition, we further compared the orderspecificity between the two groups (Fig. 3(j)(l)). (Note that it is not feasible to compare the orderspecificity between Longer and Shorterpath direction groups, or More and Fewerpaths direction groups because the relation of ordered gene pair in each group is not symmetric). We found that the orderspecificity of the FBL group was significantly smaller than that of the NonFBL group. This implies that the FBL structure reduced the specificity of the mutation order. Taken together, we can conclude that orderedmutationinducing dynamics are highly related with the structural properties such as the path length, the number of paths, and the FBL. In addition, we investigated the relationship between these structural properties and the orderedmutationinducing dynamics in a numbers of BA random networks (see Methods section) and found consistent results (see Additional file 1: Figure S1). In other words, our findings might be observed in networks with various structures. Finally, the results can be related with a recent study about the occurrence of different cancer types by the mutation order [53]. The authors in that study found that a double mutation in the order of EP300 and TP53 genes was relatively frequent in patients with esophageal and bladder urothelial carcinoma. On the other hand, the mutation in the reverse order was enriched in patients with cervical squamous cell carcinoma and endocervical adenocarcinoma. It is intriguing that the order of EP300 and TP53 belongs to â€˜Shorterpath directionâ€™, â€˜Morepaths directionâ€™, and â€˜FBLâ€™ groups in the HCS network according to our classification, all of which indicated a relatively low mutationsensitivity.
Analysis of drugtarget genes with respect to orderedmutationinducing dynamics
Some previous studies have investigated the characteristics of drugtarget genes through networkbased structural analysis [54, 55], and the findings were useful to understand tumorigenesis in cancer [56]. In this study, we extended it to dynamic analysis by investigating the orderedmutationinducing dynamics of drugtargets in signaling networks. We first specified all genes as â€˜Drugtarget (DT)â€™ and â€˜Nondrugtarget (NonDT)â€™ genes (see Methods and Additional file 1: Table S1, S2 and S3). Then, we classified every ordered gene pair in Î© into four groups: â€˜DT â†’ DTâ€™, â€˜DT â†’ NonDTâ€™, â€˜NonDT â†’ DTâ€™, and â€˜NonDT â†’ NonDTâ€™. After forcing double knockout mutations, we compared the average mutationsensitivity value among them (Fig. 4(a)(c); all Pvalues using the MannWhitney U test). As shown in the figure, the values of â€˜NonDT â†’ NonDTâ€™ and â€˜DT â†’ DTâ€™ groups were highest and lowest, and they were the bounds for the values of other groups. Furthermore, the sensitivity of the â€˜NonDT â†’ DTâ€™ group was significantly higher than that of the â€˜DT â†’ NonDTâ€™ group, for most time gap values. Considering that these two groups are identical to each other except for the order in a gene pair, the result implies that the sensitivity difference was caused by only the mutation order. We further examined the orderspecificity values of DT and NonDT groups (Fig. 4(d)(f); all Pvalues using the MannWhitney U test) and found that the former is larger than the latter. This finding is interesting considering that the mutationsensitivity of â€˜DT â†’ DTâ€™ was smaller than that of â€˜NonDT â†’ NonDTâ€™ in Fig. 4(a)(c). In other words, the network is less sensitive, but the mutation order is more critical when drugtarget genes are mutated than when nondrugtarget genes are mutated. Moreover, we further investigated the â€˜DT â†’ DTâ€™ group by classifying every gene pair in the group into â€˜Same drugâ€™ and â€˜Different drugâ€™ subgroups (see Additional file 1: Table S1, S2 and S3) for cases where both genes of a pair are targeted using the same drug or different drugs, respectively (TGL network was excluded from this analysis because there is no pair of genes belonging to â€˜Same drugâ€™ group). We found that both the mutationsensitivity and orderspecificity of the â€˜Different drugâ€™ group were less than those of the â€˜Same drugâ€™ group (Fig. 5). This implies that the network is more sensitive, and the mutation order is more specific when drugtargets from the samedrug are mutated than when drugtargets from the differentdrug are mutated. Interestingly, this observation can be linked to some previous experimental studies about multiple drug treatments. For example, a specific sequential treatment of roscovitine before doxorubicin is synthetically lethal in breast cancer cell [57] and the treatment order of doubledrugs with the shared targets is significant to the treatment efficiency [58]. In addition, our result implies that the orderedmutationinducing dynamics can be useful to predict a new drugtarget gene which may show relatively lower mutationsensitivity and the higher orderspecificity when it is subject to the ordered mutation with another drugtarget gene together.
Analysis of tumor suppressor and oncogenes with respect to orderedmutationinducing dynamics
It is known that tumor suppressors and oncogenes perform their cellular functions jointly in tumor progressions [59, 60], and tumor suppressors can be considered as therapeutic targets for cancer drugs [61, 62]. In this study, we investigated the orderedmutationinducing dynamics of tumor suppressors and oncogenes in signaling networks. We first specified all genes in the networks as â€˜Tumor suppressor genes (TSG)â€™ and â€˜Oncogenes (OCG)â€™ (see Methods and Additional file 1: Table S1, S2 and S3), and then identified two groups of ordered gene pairs in Î©, â€˜TSG â†’ OCGâ€™ and â€˜OCG â†’ TSGâ€™. For every ordered pair of genes, we computed the mutationsensitivity after forcing double knockout mutations according to the order of gene pair. Then, we compared the average mutationsensitivity value between those two groups (Fig. 6(a)(c); all Pvalues using MannWhitney U test). As shown in the figure, the mutationsensitivity value of the former group was significantly smaller than that of the latter group in all signaling networks, almost irrespective of the time gap. In other words, the network is more sensitive when oncogenes were mutated before tumor suppressors than the reverse order. In addition, we further compared the orderspecificity between two groups, â€˜TSGâ€™ and â€˜OCGâ€™, and found that the orderspecificity values of the former group were smaller than the latter group, almost irrespective of the time gap (Fig. 6(d)(f); all Pvalues were obtained using MannWhitney U test). This finding can be also related to some previous studies on the ordered mutations between oncogenes and tumor suppressor genes. For example, the double mutation in the order of TP53 and NOTCH, which are representative tumorsuppressor and oncogenes, respectively, was frequently observed in early stage of esophageal carcinoma patients [53], whereas the reverseordered mutation is likely to lead to a metastasis progression in mouse experiments [63, 64]. It was also shown that alteration of RAS, which is another oncogene, before loss of P53 formed a malignant tumor with metastatic behavior, but the reverseordered mutation resulted in benign tumors [2, 65].
Discussion
In this study, we defined the mutationsensitivity and the orderspecificity based on a Boolean network model to unravel the effects of ordered mutations on dynamics in signaling networks. It was interesting to observe that some structural properties of signaling networks can be a good indicator to explain the dynamical behavior with respect to orderedmutation experiments. In addition, it was shown that various functionally important genes are related to the orderedmutationinducing dynamic. These results can enhance the understanding of the dynamic effects of ordered doublemutations on complex dynamics of largescale biological systems, which supports the usefulness of our approach. Despite the usefulness of our approach, there are some limitations to be discussed. In this study, we employed the random nested canalyzing function to simulate the Boolean dynamics of the molecular signaling networks. This artificial specification can be a limitation of this study, although some previous studies have proven the usefulness of the model in fitting the update rules from the real biological data [36, 38]. Another concern is the synchronous update scheme, which is less realistic than the asynchronous update scheme. Therefore, a future study will include an approach to more accurately model the update rule inferred from real biological data.
Conclusions
Many previous studies investigated ordered mutations and found statistical relations with cancer development. Recently, these studies were extended to incorporate the analysis of biological networks. However, they are limited in identifying the significance of ordered mutations because they did not focus on analysis of the network dynamics. In this regard, we quantified the orderedmutationinducing dynamics by defining the mutationsensitivity and the orderspecificity measures using a Boolean network model. Specifically, they represent the probability that a network converges to a different attractor by a double knockout mutation, and the probability with which a network converges to different attractors by different mutation orders, respectively. It was not rare to observe both nonzero sensitivity and specificity values in largescale signaling networks. In addition, we examined the relationship between the structural characteristics such as the path length, the number of paths, and the feedback loop with the orderedmutationinducing dynamics in the signaling networks. Interestingly, they showed significant relationships, which implies that such structural properties need to be considered in experimental studies with respect to orderedmutation experiments. Next, we investigated the orderedmutationinducing dynamics of various functionally important genes. The numbers of drugtargets genes were negatively correlated to the mutationsensitivity, whereas the network was more specific to the order of mutations subject to drugtargets genes than the rest genes. In addition, we found that tumor suppressors can efficiently suppress the amplification of oncogenes when the former genes are mutated earlier than the latter genes. Taken together, our results enhance the understanding of the dynamic effects of ordered doublemutations on complex dynamics of largescale biological systems.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its additional files.
Abbreviations
 BA:

BarabÃ¡si Albert
 DT:

Drugtarget
 FBL:

Feedbackloop
 HCS:

Human cancer signaling
 KEGG:

Kyoto Encyclopedia of Genes and Genomes database
 NCF:

Nested canalyzing function
 OCG:

Oncogene
 TGL:

Tcell large granular lymphocyte
 TSG:

Tumor suppressor gene
References
Loeb KR, Loeb LA. Significance of multiple mutations in cancer. Carcinog. 2000;21:379â€“85.
Kent DG, Green AR. Order matters: the order of somatic mutations influences cancer evolution. Cold Spring Harb Perspect Med. 2017;7:a027060.
Ortmann CA, Kent DG, Nangalia J, Silber Y, Wedge DC, Grinfeld J, Baxter EJ, Massie CE, Papaemmanuil E, Menon S, et al. Effect of mutation order on Myeloproliferative neoplasms. N Engl J Med. 2015;372:601â€“12.
Nangalia J, Nice FL, Wedge DC, Godfrey AL, Grinfeld J, Thakker C, Massie CE, Baxter J, Sewell D, Silber Y, et al. DNMT3A mutations occur early or late in patients with myeloproliferative neoplasms and mutation order influences phenotype. Haematol. 2015;100:e438â€“42.
Beekman R, Valkhof MG, Sanders MA, van Strien PMH, Haanstra JR, Broeders L, GeertsmaKleinekoort WM, Veerman AJP, Valk PJM, Verhaak RG, et al. Sequential gain of mutations in severe congenital neutropenia progressing to acute myeloid leukemia. Blood. 2012;119:5071â€“7.
Kimberly L, Toaa A, Libia P, Maya B, George B. Size matters: sequential mutations in tumorigenesis may reflect the stochastic effect of mutagen target sizes. Genes Cancer. 2011;2:927â€“31.
Swanton C. Cancer evolution constrained by mutation order. N Engl J Med. 2015;372:661â€“3.
Drost J, van Jaarsveld RH, Ponsioen B, Zimberlin C, van Boxtel R, Buijs A, Sachs N, Overmeer RM, Offerhaus GJ, Begthel H, et al. Sequential cancer mutations in cultured human intestinal stem cells. Nat. 2015;521:43.
Sun QY, Ding LW, Tan KT, Chien W, Mayakonda A, Lin DC, Loh XY, Xiao JF, Meggendorfer M, Alpermann T, et al. Ordering of mutations in acute myeloid leukemia with partial tandem duplication of MLL (MLLPTD). Leuk. 2016;31:1.
Ascolani G, LiÃ³ P. Modelling the order of driver mutations and metabolic mutations as structures in cancer dynamics. ARXIV. 2017;2:eprint arXiv:1705â€“10862.
Guo J, Guo H, Wang Z. Inferring the temporal order of Cancer gene mutations in individual tumor samples. PLoS One. 2014;9:e89244.
Kang H, Cho KH, Zhang XD, Zeng T, Chen L. Inferring sequential order of somatic mutations during Tumorgenesis based on Markov chain model. IEEE/ACM Trans Comput Biol Bioinform. 2015;12:1094â€“103.
Gerstung M, Eriksson N, Lin J, Vogelstein B, Beerenwinkel N. The temporal order of genetic and pathway alterations in tumorigenesis. PLoS One. 2011;6:e27136.
Lecca P, Casiraghi N, Demichelis F. Defining order and timing of mutations during cancer progression: the TODAG probabilistic graphical model. Front Genet. 2015;6:309.
Misra N, Szczurek E, Vingron M. Inferring the paths of somatic evolution in cancer. Bioinform. 2014;30:2456â€“63.
Cornelius SP, Kath WL, Motter AE. Realistic control of network dynamics. Nat Commun. 2013;4:1942.
Jackson FLC, Niculescu MD, Jackson RT. Conceptual shifts needed to understand the dynamic interactions of genes, environment, epigenetics, social processes, and behavioral choices. Am J Public Health. 2013;103:S33â€“42.
Cornelius SP, Lee JS, Motter AE. Dispensability of Escherichia coliâ€™s latent pathways. Proc Natl Acad Sci. 2011;108:3124â€“9.
Motter AE, Gulbahce N, Almaas E, BarabÃ¡si AL. Predicting synthetic rescues in metabolic networks. Mol Syst Biol. 2008;4.
Wytock TP, Fiebig A, Willett JW, Herrou J, Fergin A, Motter AE, Crosson S. Experimental evolution of diverse Escherichia coli metabolic mutants identifies genetic loci for convergent adaptation of growth rate. PLoS Genet. 2018;14:e1007284.
Cui Q, Ma Y, Jaramillo M, Bari H, Awan A, Yang S, Zhang S, Liu L, Lu M, O'ConnorMcCourt M, et al. A map of human cancer signaling. Mol Syst Biol. 2007;3:152.
Cui Q, Purisima EO, Wang E. Protein evolution on a human signaling network. BMC Syst Biol. 2009;3:21.
Kim JR, Kim J, Kwon YK, Lee HY, HeslopHarrison P, Cho KH. Reduction of Complex Signaling Networks to a Representative Kernel. Sci Signal. 2011;4:ra35.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D457â€“62.
Motter AE, Albert R. Networks in motion. Phys Today. 2012;65:43â€“8.
Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, Albert R, Loughran TP. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci. 2008;105:16308.
Knox C, Law V, Jewison T, Liu P, Ly S, Frolkis A, Pon A, Banco K, Mak C, Neveu V, et al. DrugBank 3.0: a comprehensive resource for â€˜Omicsâ€™ research on drugs. Nucleic Acids Res. 2011;39:D1035â€“41.
Zhao M, Sun J, Zhao Z. TSGene: a web resource for tumor suppressor genes. Nucleic Acids Res. 2013;41:D970â€“6.
Zhao M, Kim P, Mitra R, Zhao J, Zhao Z. TSGene 2.0: an updated literaturebased knowledgebase for tumor suppressor genes. Nucleic Acids Res. 2016;44:D1023â€“31.
Liu Y, Sun J, Zhao M. ONGene: a literaturebased database for human oncogenes. J Genet Genomics. 2017;44:119â€“21.
Helikar T, Konvalina J, Heidel J, Rogers JA. Emergent decisionmaking in biological signal transduction networks. Proc Natl Acad Sci. 2008;105:1913.
Kwon YK. Properties of Boolean dynamics by node classification using feedback loops in a network. BMC Syst Biol. 2016;10:83.
Raeymaekers L. Dynamics of Boolean networks controlled by biologically meaningful functions. J Theor Biol. 2002;218:331â€“41.
Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cellcycle network is robustly designed. Proc Natl Acad Sci U S A. 2004;101:4781.
Trinh HC, Le DH, Kwon YK. PANET: a GPUbased tool for fast parallel analysis of robustness dynamics and feedforward/feedback loop structures in largescale biological networks. PLoS One. 2014;9:e103010.
Kauffman S, Peterson C, Samuelsson B, Troein C. Random Boolean network models and the yeast transcriptional network. Proc Natl Acad Sci. 2003;100:14796â€“9.
Kauffman S, Peterson C, Samuelsson B, Troein C. Genetic networks with canalyzing Boolean rules are always stable. Proc Natl Acad Sci U S A. 2004;101:17102â€“7.
Samal A, Jain S. The regulatory network of E. coli metabolism as a Boolean dynamical system exhibits both homeostasis and flexibility of response. BMC Syst Biol. 2008;2:21.
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:23â€“40.
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:e1000912.
Campbell C, Albert R. Stabilization of perturbed Boolean network attractors through compensatory interactions. BMC Syst Biol. 2014;8:53.
Kwon YK, Kim J, Cho KH. Dynamical robustness against multiple mutations in signaling networks. IEEE/ACM Trans Comput Biol Bioinform. 2016;13:996â€“1002.
Jilkine A, Gutenkunst RN. Effect of dedifferentiation on time to mutation Acquisition in Stem CellDriven Cancers. PLoS Comput Biol. 2014;10:e1003481.
Turajlic S, McGranahan N, Swanton C. Inferring mutational timing and reconstructing tumour evolutionary histories. Biochim Biophys Acta (BBA)  Rev Cancer. 2015;1855:264â€“75.
Klein C, Marino A, Sagot MF, Vieira Milreu P, Brilli M. Structural and dynamical analysis of biological networks. Briefings Funct Genomics. 2012;11:420â€“33.
Thomas R, Thieffry D, Kaufman M. Dynamical behaviour of biological regulatory networksâ€”I. Biological role of feedback loops and practical use of the concept of the loopcharacteristic state. Bull Math Biol. 1995;57:247â€“76.
Ananthasubramaniam B, Herzel H. Positive feedback promotes oscillations in negative feedback loops. PLoS One. 2014;9:e104761.
BarabÃ¡si AL, Albert R. Emergence of scaling in random networks. Sci. 1999;286:509.
Schoonjans F, Zalata A, Depuydt CE, Comhaire FH. MedCalc: a new computer program for medical statistics. Comput Methods Prog Biomed. 1995;48:257â€“62.
Mazaya M, Trinh HC, Kwon YK. Construction and analysis of genegene dynamics influence networks based on a Boolean model. BMC Syst Biol. 2017;11:133.
Prill RJ, Iglesias PA, Levchenko A. Dynamic properties of network motifs contribute to biological network organization. PLoS Biol. 2005;3:e343.
Le DH, Kwon YK. The effects of feedback loops on disease comorbidity in human signaling networks. Bioinformatics. 2011;27:1113â€“20.
Li X. Dynamic changes of driver genesâ€™ mutations across clinical stages in nine cancer types. Cancer Med. 2016;5:1556â€“65.
Kotlyar M, Fortney K, Jurisica I. Networkbased characterization of drugregulated genes, drug targets, and toxicity. Methods. 2012;57:499â€“507.
YÄ±ldÄ±rÄ±m MA, Goh KI, Cusick ME, BarabÃ¡si AL, Vidal M. Drugâ€”target network. Nat Biotechnol. 2007;25:1119.
DurmuÅŸ Tekir S, YalÃ§Ä±n Arga K, Ãœlgen KÃ–. Drug targets for tumorigenesis: insights from structural analysis of EGFR signaling network. J Biomed Inform. 2009;42:228â€“36.
JabbourLeung NA, Chen X, Bui T, Jiang Y, Yang D, Vijayaraghavan S, McArthur MJ, Hunt KK, Keyomarsi K. Sequential combination therapy of CDK inhibition and doxorubicin is synthetically lethal in p53mutant triplenegative breast Cancer. Mol Cancer Ther. 2016;15:593.
Koplev S, Longden J, FerkinghoffBorg J, Blicher BjerregÃ¥rd M, Cox TR, Erler JT, Pedersen JT, Voellmy F, Sommer MOA, Linding R. Dynamic rearrangement of cell states detected by systematic screening of sequential anticancer treatments. Cell Rep. 2017;20:2784â€“91.
Lee EYHP, Muller WJ. Oncogenes and tumor suppressor genes. Cold Spring Harb Perspect Biol. 2010;2.
Zhu K, Liu Q, Zhou Y, Tao C, Zhao Z, Sun J, Xu H. Oncogenes and tumor suppressor genes: comparative genomics and network perspectives. BMC Genomics. 2015;16:S8.
Morris LGT, Chan TA. Therapeutic targeting of tumor suppressor genes. Cancer. 2015;121:1357â€“68.
Li J, Hao D, Wang L, Wang H, Wang Y, Zhao Z, Li P, Deng C, Lj D. Epigenetic targeting drugs potentiate chemotherapeutic effects in solid tumor therapy. Sci Rep. 2017;7:4035.
Chanrion M, Kuperstein I, BarriÃ¨re C, El Marjou F, Cohen D, Vignjevic D, Stimmer L, PaulGilloteaux P, BiÃ¨che I, Tavares SDR, et al. Concomitant Notch activation and p53 deletion trigger epithelialtomesenchymal transition and metastasis in mouse gut. Nat Commun. 2014;5:5005.
Jen KY, Song IY, Banta KL, Wu D, Mao JH, Balmain A. Sequential mutations in Notch1, Fbxw7, and Tp53 in radiationinduced mouse thymic lymphomas. Blood. 2012;119:805.
Herbet M, Salomon A, Feige JJ, Thomas M. Acquisition order of Ras and p53 gene alterations defines distinct adrenocortical tumor phenotypes. PLoS Genet. 2012;8:e1002700.
Acknowledgements
This work was supported by the 2019 Research Fund of University of Ulsan.
About this supplement
This article has been published as part of BMC Medical Genomics, Volume 13 Supplement 4, 2020: 18th International Conference on Bioinformatics. The full contents of the supplement are available at https://bmcmedgenomics.biomedcentral.com/articles/supplements/volume13supplement4.
Funding
Publication charges for this work were funded by the 2019 Research Fund of University of Ulsan. The funding body has not played any roles in the design of the study and collection, analysis and interpretation of data in writing the manuscript.
Author information
Authors and Affiliations
Contributions
YKK conceived the study, MM and THC performed simulations, MM, THC, and YKK designed the analysis, MM and YKK wrote and revised the manuscript. All authors read and approved the final manuscript.
Corresponding author
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.
Additional information
Publisherâ€™s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Additional file 1 Figure S1.
Relations of structural properties with orderedmutationinducing dynamics in BA network. A total of 250 BA random networks with Vâ€‰=â€‰50 and Aâ€‰=â€‰100 were generated. The time gap (T) was set to 1â€“10. (a) Mutationsensitivity result with respect to the shortest path length. All pairs of nodes involving an FBL were classified into â€˜Shorterpath directionâ€™ and â€˜Longerpath directionâ€™ groups according that l(v_{i},â€‰v_{j})â€‰<â€‰l(v_{j},â€‰v_{i}) and l(v_{i},â€‰v_{j})â€‰>â€‰l(v_{j},â€‰v_{i}), respectively. (b) Mutationsensitivity result with respect to the number of paths. All pairs of nodes were classified into â€˜Morepaths directionâ€™ and â€˜Fewerpaths directionâ€™ groups according that n(v_{i},â€‰v_{j})â€‰>â€‰n(v_{j},â€‰v_{i}) and n(v_{i},â€‰v_{j})â€‰<â€‰n(v_{j},â€‰v_{i}), respectively. (c) Mutationsensitivity result with respect to the FBLs. All pairs of nodes were classified into â€˜FBLâ€™ and â€˜NonFBLâ€™ groups, according that any gene of the pair is involved by an FBL or not. (d) Orderspecificity result with respect to the FBLs. All Pvalues were computed using the MannWhitney U test. Table S1. Gene information of HCS consisting 1192 genes, including its association with drugtarget, tumor suppressor, and oncogene. Table S2. Gene information of KEGG consisting 1659 genes, including its association with drugtarget, tumor suppressor, and oncogene. Table S3. Gene information of TGL consisting 61 genes, including its association with drugtarget, tumor suppressor, and oncogene.
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
Cite this article
Mazaya, M., Trinh, HC. & Kwon, YK. Effects of ordered mutations on dynamics in signaling networks. BMC Med Genomics 13 (Suppl 4), 13 (2020). https://doi.org/10.1186/s129200190651z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s129200190651z