Quantitative evaluation and reversion analysis of the attractor landscapes of an intracellular regulatory network for colorectal cancer

Background Cancer reversion, converting the phenotypes of a cancer cell into those of a normal cell, has been sporadically observed throughout history. However, no systematic analysis has been attempted so far. Results To investigate this from a systems biological perspective, we have constructed a logical network model of colorectal tumorigenesis by integrating key regulatory molecules and their interactions from previous experimental data. We identified molecular targets that can reverse cancerous cellular states to a normal state by systematically perturbing each molecular activity in the network and evaluating the resulting changes of the attractor landscape with respect to uncontrolled proliferation, EMT, and stemness. Intriguingly, many of the identified targets were well in accord with previous studies. We further revealed that the identified targets constitute stable network motifs that contribute to enhancing the robustness of attractors in cancerous cellular states against diverse regulatory signals. Conclusions The proposed framework for systems analysis is applicable to the study of tumorigenesis and reversion of other types of cancer. Electronic supplementary material The online version of this article (doi:10.1186/s12918-017-0424-2) contains supplementary material, which is available to authorized users.

which is closely related to DNA repair and drug resistance [3,17]. In summary, we have integrated important subnetworks related to colorectal cancer and thereby established an essential intracellular regulatory network model of colorectal cancer by including critical marker genes that represent each subnetwork.
It is indeed true that SLUG and Snail are well-known regulators of EMT. On the other hand, they are also known as marker molecules highly associated with the stemness of colorectal cancer [3]. Recently, overexpression of Snail was shown to induce not only EMT but also a cancer stem cell-like phenotype in human colorectal cancer cells [18]. Moreover, Snail is known to be required for the expression of putative markers of stem cells, and can induce the expression of stemness-promoting genes, such as Nanog, KLF4, and TCF-4 [19]. Therefore, we can consider that SLUG and Snail might play critical roles in determining both EMT and stemness phenotypes in colorectal cancer. In our network model, the SLUG and Snail signals were transmitted to E-cadherin and MMP which are also known to be critical regulators of EMT in colorectal cancer. Hence, we chose E-cadherin and MMP as key marker genes of the EMT module to determine the functional phenotype of EMT, and assigned SLUG and Snail to the stemness module and set them as key marker molecules.
In this study, we aimed to construct a simple but essential regulatory network model of a colorectal cancer cell based on extensive literature survey and data-driven network inference. There are indeed several important regulators not included in our network model, such as TGFβ, SMAD4 or DCC. However, we have integrated the mutation and CNA profiles of these regulators when we construct the network model, so their regulatory effects were implicitly reflected in our network model. For instance, for a cell line with the genetic alteration of TGFβ, we fixed the expression levels of c-Myc and p21 as these are the target nodes of TGFβ [20,21]. As a result, we could successfully reproduce the expression patterns and malignancy ranks of the cancer organoids, as shown in Fig. 3A [22]. Likewise, to reflect the effect of SMAD4, we controlled the expression levels of p21 and Snail [23,24], which are regulated by SMAD4, and thereby we could successfully find out the cancer reversion phenomena (Fig. 3D). In summary, we tried 5 to construct a simple but essential regulatory network model with minimal nodes and to reflect implicitly the effects of some other missing regulators to the network wiring and regulatory logics. In this way, we have successfully recapitulated both the dynamics of the regulatory network of a colorectal cancer cell and previous experimental observations on cancer reversion phenomena. 6

Text S2. Boolean operator logic construction and weighted sum logic conversion
The states of biological nodes such as expression levels of genes or the activities of proteins can be represented by two quantitatively distinguishable states without any intermediate states. A Boolean network model is constructed solely with nodes that are either 'ON' or 'OFF' state in order to mimic the biological interactions effectively. Boolean network is beneficial for a network that contains well-defined structure of the nodes without having to know kinetic parameters. The three Boolean operators can regulate its nodes: AND, OR, and NOT. If a biological molecule acts inside a cell, it means that the molecule exists and be activated. If a biological molecule exists inside a cell, it means the molecule is frequently transcribed or stable. Thus, we can denote the state of a specific node in time step t + 1 with the Boolean operator logic equation of input node states as below: Where + represents a molecular activity state of ℎ node at time step t + 1 , represents input signals for transcriptional activation, represents input signals for stabilization and represents input signals for functional activation. The transcriptional activation is a biological interaction related with gene expression level changes such as transcription factor-target gene(TG-TF) interaction. The stabilization is a biological interactions related with molecular stability such as ubiquitination, degradation or polymer formation. The activation is a biological interactions related with molecular activity such as phosphorylation or methylation. As a result, we have constructed 34 Boolean operator logic equations based on the action mechanisms of each molecule in our network.
Although the Boolean operator logic is appropriate for reflecting biological action mechanisms to construct the update logic, it is hard to project the mutation and copy number variation data differently. Also, both defining the node state without any input node condition and adjusting the logic based on the experimental data are difficult in Boolean operator logic. The reason of these limitation is that there is no value where can adjust the node state quantitatively because the Boolean operator logic is a 7 qualitative equation. However, all of these problems can be solved when we use a weighted sum logic, a logic which defines both the basal node states and node relationships quantitatively as weights. Thus, not only projecting biological action mechanisms but also making easier to quantitatively adjust the relations, The next state + is then determined by the next input.
The basic algorithms of this conversion are (i) define the basal levels based on the prior biological knowledge, (ii) generate a standard truth  If the increase of normal-like score during the double-node perturbation is larger, similar or smaller than the sum of the normal-like scores from the individually perturbed nodes A and B, then each incident can be classified as synergistic, additive or antagonistic, respectively. The range between increased normal-like scores has been calculated as 20% in average after performing double node perturbation in various pairs.
As a result, we are able to separate the synergistic pairs from various combination of two nodes perturbations and analyze the specific mechanism of them which have lead us to reveal the functional motifs. the method of evaluating network robustness suggested in our study uses the average basin size of five major attractors rather than one single largest attractor, our method can be better in measuring the robustness of the whole network. Note also that we used both basin size and number of attractors to measure the network robustness, which enables to determine the ranking of network robustness at each tumorigenesis step. Using this concept of network robustness, we were able to trace the network robustness during tumorigenesis and could finally use the information to find out a novel strategy for cancer reversion. combinations for four or more link deletion due to computational limitation. After all, we have calculated the average increase of normal-like score for every case. We found that our network model is reasonably robust with respect to the increase of normal-like score for about six random link deletion (see Fig. S5A), indicating that our results remain robust even if the network structure is changed in that range. In addition, to check the dynamical robustness of our network model, we have added random noise from 0 to 15% of each element value of the connectivity matrix and the basal level column vector. After each simulation, we have calculated the average increase of normal-like score for all different levels of random noise. We found that our network model again is reasonably robust against up to 10% noise level (see Fig. S5B), indicating that the dynamical property of our network model is robust against such perturbations.  [27]. As a result, the cancer reversion for the perturbation of the suggested 11 cancer reversion targets is found to be statistically significant (Fig. S6), which indicates that the increase of the normal-like score is large enough to confirm successful cancer reversion by controlling the suggested targets. The one or two stars on top of each bar of the graph represent the p-value below 0.05 or 0.001 respectively. Blue, red and green bar in the figure represent three different types of perturbations, activation, inhibition, and restoration, respectively. The increase of the normal-like score for most of the reversion targets was found to be statistically significant compared to the average increase of the normal-like score when we perturb non-reversion targets as indicated by an empty bar with solid line according to each perturbation type. and p21(SMAD4), which were known to induce colorectal cancer when sequentially mutated. As a result, no significant synergistic effect was observed among these pairs, as shown in Fig. R4. Although Bax and Caspase-3 are likely to be synergistic for every pair, they are not synergistic in terms of cancer reversion since they continuously induce apoptosis. Except for them, almost every synergistic pairs are included in the stable motifs shown in Fig. 5A. The pairs among APC, P53, KRAS and SMAD4 are not significantly the EGF signal transmission motif in Fig. 4C, the synergy of cancer reversion can emerge when multiple nodes of a single signaling pathway were perturbed simultaneously. Thus, adjacent nodes can have a higher chance to cause such a synergistic effect than distant nodes. However, we found that four of the frequently mutated nodes in colorectal cancer are located in different modules and are relatively far apart in our network model, which might be one of the reasons why no significant synergistic effect was observed among them. The blue color indicates a synergistic pair while the red color indicates an antagonistic pair.
The pairs among APC, p53, KRAS, and p21 (SMAD4) were denoted by red squares. The darker the color, the more synergistic or antagonistic the pair is.