Abstract
Cancer heterogeneity may reflect differential dynamical outcomes of the regulatory network encompassing biomolecules at both transcriptional and post-transcriptional levels. In other words, differential gene-expression profiles may correspond to different stable steady states of a mathematical model for simulation of biomolecular networks. To test this hypothesis, we simplified a regulatory network that is important for soft-tissue sarcoma metastasis and heterogeneity, comprising of transcription factors, micro-RNAs, and signaling components of the NOTCH pathway. We then used a Boolean network model to simulate the dynamics of this network, and particularly investigated the consequences of differential miRNA degradation modes. We found that efficient miRNA degradation is crucial for sustaining a homogenous and healthy phenotype, while defective miRNA degradation may lead to multiple stable steady states and ultimately to carcinogenesis and heterogeneity.
Heterogeneity is one of the key features of cancer. Even for cancer cells of the same tissue, their gene expression profiles show great heterogeneity. We hypothesize that these heterogeneous expression profiles correspond to different steady states of the biomolecular network underlying cancer development. In normal cells, the network dynamics converge to only a few steady states (e.g. states A, B and C), and only one steady state (e.g. state A) is actually viable. For example, if state A corresponds to the viable phenotype, then cells of states B and C would possess some kind of defect and have to undergo apoptosis. As a result, normal cells manifest great homogeneity. In comparison, cancer signaling networks may possess many more steady states (e.g. 100 steady states). Although some of these correspond to non-viable phenotypes, there are still many steady states (e.g. states A, B, C, D, E, F) corresponding to viable and robust phenotypes, such as, for instance, six abnormal gene-expression patterns that govern carcinogenesis. The tumor, thus consists of cells of heterogeneous gene expression: a cell can adopt any one of the six gene-expression profiles.
In the present article, we use soft-tissue sarcomas (STS) as an example to test the above hypotheses through computational modeling of the cancer signaling network. STS are heterogeneous cancers of mesenchymal tissues, comprising more than 50 histological subtypes derived from muscle, fat, blood vessels, nerves, tendons and the synovia. Some 40% to 50% of patients with STS will develop lethal metastases (1). Soft-tissue sarcoma is notoriously heterogeneous–even cells of the same subtype of STS have significantly different gene expression (2, 3). After examining a variety of STS gene expression profiles, Samantarrai et al. found that the NOTCH signaling pathway is markedly aberrant in STS cells (4). The aberration extends to regulators of NOTCH signaling, including transcription factors (TF) such as Finkel-Biskis-Jinkins murine osteosarcoma viral oncogene homolog (FOS), v-rel avian reticuloendotheliosis viral oncogene homolog A (RELA), tumor protein p53 (TP53), v-myc avian myelocytomatosis viral oncogene homolog (MYC) and micro-RNAs (miRNAs) such as miR-20b, miR-21, miR-29b, miR-34a. This greater NOTCH signaling pathway was illustrated elsewhere [Figure 3 of (4)]. To highlight these regulatory steps (TF → gene, TF → miRNA, miRNA → gene), we simplified the greater NOTCH signaling network into a core network comprising of 12 nodes and 26 edges (Figure 1). We refer to this simplified version as the miRNA-TF-NOTCH network throughout the present article. We then used a Boolean network model to simulate the dynamical changes of the miRNA-TF-NOTCH network, namely the ON/OFF switches of all the 12 nodes through a succession of discrete time steps. For this network, we found that the dynamical changes always converge to a steady state. Because the steady states are not unique, the actual steady state at which the dynamics finally arrive will depend on the initial state at which the simulation began. We use the total number of steady states to measure the degree of heterogeneity.
Given the importance of miRNA degradation in STS development and metastasis, we investigated how the modes of degradation influence network dynamics, especially with regard to the number of steady states. We considered the two following modes of degradation: Mode 1: miRNAs undergo efficient self-degradation. Mode 2: miRNAs are stable by default; they may degrade only after interacting with other molecules (5-8). The network dynamics differs under the two modes of miRNA degradation. We used computer simulation to compare the two dynamics, which revealed their possible correspondence with physiological or various pathological phenotypes.
Materials and Methods
The Boolean network model. Given the complexity of biomolecular networks, mathematical modeling based on traditional differential equations is difficult to use to obtain the system's steady state in an acceptable time. We thus use logic-based (Boolean) approaches to model biomolecular networks here. A Boolean network consists of N kinds of interacting molecules, each of which is modeled as being either ON (active) or OFF (inactive). At any given time, these ON/OFF combinations constitute a network state. Over time, the system dynamically changes from one network state to another, depending on the interactions between the molecules. Thus, from a given state at the start, there is a well-defined sequence of network states that finally converge to a fixed network state (known as the attractor in dynamical systems theory).
The dynamics of a Boolean network model (determining the next state from the current state) can be described as follows (9, 10): where rji represents an inhibitory (red) arrow from node j to node i; gji represents a stimulatory (green) edge from node j to node i; “addition” represents the Boolean operator OR; “multiplication” represents the Boolean operator AND; the bar on a variable represents the Boolean operator NOT. Note that rii represents self-degradation of node i. Indeed, a molecule does not inhibit itself, thus rii represents degradation instead of inhibition. An edge is defined as active at time t only when its source node is active at time t. For example, gji=1, sj(t)=1, sj(t')=0 means that there is a green edge from node j to node i, and the edge is active at time t and inactive at time t'.
Figure 1 depicts an example: r28=1 and g28=0 because there is a red arrow pointing from node 2 to node 8; r59=0 and g59=1 because there is a green arrow pointing from node 5 to node 9. Figure 1 does not provide information of degradation, thus here, rii remains as a variable for all i=1, 2, …, 12.
miRNA-TF-NOTCH regulatory network. By removing those peripheral nodes having low connectivity with the major nodes, the greater NOTCH signaling network described in (4) is simplified into a core network comprising 12 nodes and 26 edges (Figure 1). This simplified network is called the miRNA-TF-NOTCH network. In the network, the four purple nodes shown in Figure 1 represent the four miRNAs whose degradation behaviors were the focus of the present study. We considered the two modes of miRNA degradation: (i) Mode 1: Active degradation. In terms of our discrete time dynamical model, if an miRNA is ON (i.e. present at a high concentration) at time step t, then it switches off automatically (i.e. the concentration becomes low) at time step t + 1, unless some other node(s) in the network stimulate the miRNA at time step t. If a miRNA (node i) is active at time t [i.e. si(t)=1], then it automatically becomes inactive at time t + 1 [i.e. si(t + 1)=0], unless there is an active green edge pointing to it. As far as the computational model is concerned, this means rii=1 for i=1, 2, 3, and 4. This actually removes the term Eq. 1 for nodes i=1, 2, 3 and 4 only.
(ii) Mode 2: Degradation not automatic. Only if one of its targeted nodes is active at time t, does the miRNA degrade. Under this condition, rii(t) (for i=1, 2, 3, and 4) become time-dependent. For example, node 1 in Figure 1 (miR-20b) has one targeted node (node 6, namely signal transducer and activator of transcription 3 (STAT3)); thereby r11(t)=s6(t). That is, if node 6 is active at time t [s6(t)=1], then node 1 degrades [r11(t)=1]; if node 6 is inactive at time t [s6(t)=0], then node 1 does not degrade [r11(t)=0]. Node 2 (miR-21) has two targeted nodes (nodes 6 and 8); thereby r22(t)=s6(t) + s8(t). Note that addition means OR in our model; thus node 2 degrades when either node 6 or node 8 is active.
Both of these two modes were considered in our simulations. Starting from an arbitrary initial network-state (e.g. the t=0 row of Table I), one obtains the next network state (e.g. the t=1 row of Table I) by Eq. 1, according to the network connections shown in Figure 1 and one of the two degradation modes. The subsequent network states are obtained in the same way. When Eq. 1 no longer produces new network states, a steady state has been reached. For example, t=5 row of Table I is identical with the t=4 row of Table I, indicating that t=5 row is a steady state. One then starts from another initial network state to obtain a new state transition trajectory.
Figure 2 gives an overview of all state transitions, where each blue dot represents one network state. Because the network has 12 nodes and each node has only two states, there are 212=4096 states, namely 4096 blue dots in Figure 2. Among the 4096 states, there are only nine steady states, namely A, B, C, … I. The collection of all the network states that converge to a steady state is called the basin of that steady state; and the number of network states in a basin is called the basin size. For example, the basin of steady state A is illustrated in the upper-left corner of Figure 2. Because there are 2764 states in the basin, the basin size of A is 2764. More intuitively, we say that the steady state A comprises 2764 states.
Computational facilities and methods. The computations were performed on the computer cluster of the Department of Biology, South University of Science and Technology of China, consisting of five AMAX ServMax servers (https://www.amax.com/default.asp). The software used to implement the mathematical model (i.e. Eq. 1) was Mathematica 10.3.0.0. Because the model is deterministic, no probability and statistics are involved in the present study. For example, the number of iterations is completely determined by Eq. 1 and the network connections in Figure 1. In Figure 2, each arrow corresponds to one transition of a network state, namely one iteration of Eq. 1.
Results
Simulation of network dynamics with miRNA degradation mode 1. Table I gives an example sequence obtained with degradation mode 1, where the t=0 row gives an initial network-state (0 1 0 0 0 0 0 0 0 1 1 0), namely the molecules miR-21, Yin Yang 1 (YY1), specificity protein 1 transcription factor (SP1) are active and molecules miR-20b, miR-29b, miR-34a, RELA, STAT3, STAT4, MYC, TP53 and FOS are inactive. Note that the t=4 and t=5 rows are identical, which implies that the network state has reached a steady state, after four time steps of molecular interactions dictated by the network structure in Figure 1. The steady state is (0 1 0 1 1 0 1 0 1 0 0 1). This corresponds to the stable phenotype, namely when molecules RELA, miR-21, STAT4, TP53, miR-34a and FOS are active (or highly expressed) and molecules STAT3, MYC, miR-29b, YY1, SP1 and miR-20b are inactive.
Note that Table I only represents the trajectory starting from one particular initial network state. If the initial state changes, then the trajectory changes, but the final steady state may or may not change. Because there are in total 212=4096 states, there are 4096 trajectories and that shown in Table I is only one of them. Figure 2 gives a full view of the 4096 states (represented by the blue dots) and the direction of transitions among them (the orange arrows). These 4096 states converge to nine steady states A, B, … I. Steady state A is the state described above (0 1 0 1 1 0 1 0 1 0 0 1). It is actually the steady state with the largest basin size, comprising 2764 states.
Table II presents the identities of all nine steady states, as well as their basin sizes.
Simulation of network dynamics with miRNA degradation mode 2. With degradation mode 2, the state transition landscape is shown in Figure 3, which has 58 steady states, which is many more than in the degradation mode 1. Table III presents nine steady states whose basin sizes are the largest among the total 58 steady states using this mode. Even the steady state with the largest basin size comprises only 995 states, far fewer than 2764, the largest basin size found with degradation mode 1.
Mapping steady states to physiological/pathological phenotypes. Figure 4 presents the histogram of basin size for both mode 1 (red) and mode 2 (green). Mode 1 is characterized by one super stable steady state, namely the one comprising 2764 states. This super stable steady state may correspond to the physiological (healthy) phenotype, for a deviation from the steady state can easily be corrected after at most four steps of adjustments. One may argue that there is still a total chance of (4096-2764)/4096=33% that the deviated state finally returns to one of the other eight states (see Figure 2). Although that is possible, the other eight steady states may correspond to nonviable phenotypes, with which the cells will undergo apoptosis. In this event, at the macroscopic level is healthy tissue comprising of homogeneous cells, namely cells expressing miR-21, miR-34a, RELA, STAT4, TP53 and FOS, but with inactive miR-20b, miR-29b, STAT3, MYC, YY1 and SP1.
In comparison, mode 2 of miRNA degradation is characterized by the co-existence of multiple steady states. There are in total 58 steady states, and many of them are of comparable basin sizes. As a consequence, even the tallest green bar apparent in Figure 4 is much shorter than the middle red bar (which represents 2764, the largest basin size under mode 1). This multiplicity of stable steady states may correspond to cancer heterogeneity. The two most stable steady states, namely the one with basin size 995 and that with basin size 983, are equally stable and may correspond respectively to two dominant cancer sub-types. The actual number of viable cancer phenotypes may well be much larger than two. Taking a conservative assumption that only one in 10 of these 58 steady states are viable, there are about six pathological phenotypes. The great heterogeneity is certainly a characteristic of cancer in general and STS in particular.
Our analysis revealed that the mode of miRNA degradation is crucial to the correct regulation of the NOTCH signaling pathway and ultimately to a person's well-being. Micro RNAs are small non-coding RNA molecules of about 22 nucleotides found in plants, animals, and some viruses, which function in RNA silencing and post-transcriptional regulation of gene expression. According to our Boolean network analysis, miRNA should undergo efficient degradation (mode 1) to correctly regulate TF-NOTCH signaling. If miRNA somehow becomes resistant to degradation, then the network state would be stuck at one of many stable states that may correspond to various cancer phenotypes. Clinically this may manifest as different gene-expression profiles of different cells with the same macroscopic tumor type.
Discussion
Soft-tissue sarcomas are tumors of mesenchymal tissues with a high degree of heterogeneity, comprising of more than 50 histological subtypes derived from muscle, fat, blood vessels, nerves, tendons and the lining of the joints. STS also possess a high metastatic potential (including primary site escape, local invasion, intravasation, extravasation, colonization, expansion). They account for more than 90% of cancer-associated mortality and are thus of extreme clinical relevance. Both aspects of malignancy (i.e. heterogeneity and metastasis) are related to de-regulation at the transcriptional and post-transcriptional levels.
As we have shown, a complex network of interactions between biomolecules such as TFs, mRNAs, miRNAs and protein kinases that ultimately account for the different physiological or pathological phenotypes. Computational systems biology thus can provide a better understanding of cancer biology. It can also serve for the planning of individualized cancer therapy and determination of prognostic data. Samantarrai et al. constructed a moderately large miRNA-TF network that modulates NOTCH signaling pathway through curated regulations, with a focus on the mechanisms of metastasis (4).
We focused more on the heterogeneity aspect of STS complexity. We first simplified the network constructed by Samantarrai et al. into a minimal network comprising of 12 nodes and 26 edges. This simplified network allowed us to use Boolean network modeling to reveal systems level properties of miRNA-TF regulation of NOTCH signaling pathway in STS with satisfactory computational efficiency. Study of a larger network by Boolean modeling is possible but computationally costly. One key assumption of our Boolean network modeling is that the number of steady states correlates strongly with the degree of heterogeneity, namely the number of subtypes of a macroscopic STS tumor.
Because miRNA degradation appears to be important to STS development, we studied how the two modes of miRNA degradation affect the number of steady states. We found that with the first mode (fast degradation), the system usually converges to a single steady state, which may correspond to homogeneity of normal tissue. With the second mode (conditional degradation), the system has many stable steady states, which may correspond to heterogeneous cancer phenotypes.
We, therefore, propose that under normal conditions, the miRNA action is transient, and efficient degradation of miRNA is crucial to the cell's normal functioning. Given the power of miRNAs in gene silencing and regulation, miRNAs are double-edged swords – their actions are necessary but prolonged actions would do harm to the cell. Once their tasks have been completed, miRNAs have to be degraded rapidly to avoid any ill effects. Indeed, we have shown that defective miRNA degradation generates many stable steady states. Since the normal phenotype usually corresponds to only one steady state, the other steady states may well correspond to pathological phenotypes such as cancer. These discoveries may stimulate studies on the mapping between our theoretical steady states and existing experimental data, such as differential profiles of gene expression or protein concentration within the same tumor. Such intratumoral heterogeneity has been well documented in recent years (11-13). Our results would also provide valuable insights into future experimental investigations.
This work also demonstrates the power of Boolean modeling. For complex interactions such as miRNA–TF network of STS, the strength of detailed interaction between a pair of molecules is far less important than the topology of the network. This allows molecular concentrations to be condensed simply to one of two molecular states, ON or OFF. Moreover, interactions are modeled as either stimulatory or inhibitory. Such abstractions allow an efficient analysis of complex networks with possibly some reduction of precision. In comparison, mathematical models based on differential equations are usually difficult for studying large-scale networks due to their computational complexity.
Acknowledgements
This work was partly supported by the National Natural Science Foundation of China (no. 61471186), Startup Funds from Education Ministry of China, and Shenzhen Municipal Research Funds (JCYJ20140417105816347), Fundamental Research Grant of South University of Science and Technology of China (FRG-SUSTC1501A-38).
- Received November 24, 2015.
- Revision received December 16, 2015.
- Accepted December 29, 2015.
- Copyright© 2016 International Institute of Anticancer Research (Dr. John G. Delinassios), All rights reserved