A Boolean Model of the Proliferative Role of the lncRNA XIST in Non-Small Cell Lung Cancer Cells

Simple Summary The lncRNA XIST has been verified as an oncogenic gene in non-small cell lung cancer (NSCLC), whereas the tumor suppressors miR-34a, miR-449a, and miR-16 were downregulated in NSCLC. The lncRNA XIST is known to bind to the miR-449a and miR-34a sponges, however, the biological functions linking all these non-coding RNAs (ncRNAs) are still unknown in NSCLC. This study aims to perform dynamic analysis of the gene regulation between these two ncRNAs in NSCLC. Thus, we presented a Boolean model of the plausible network connecting the lncRNA XIST/miR-34a/miR-449a and miR-16 in NSCLC. We observed that miR-449a/miR-34a regulates miR-16 and p21 expression by targeting HDAC1, c-Myc, and the lncRNA XIST. Our results demonstrate that the lncRNA XIST is an attractive target of drug development in NSCLC and that favorable outcomes can be achieved through tumor suppressor miRNAs. Abstract The long non-coding RNA X inactivate-specific transcript (lncRNA XIST) has been verified as an oncogenic gene in non-small cell lung cancer (NSCLC) whose regulatory role is largely unknown. The important tumor suppressors, microRNAs: miR-449a and miR-16 are regulated by lncRNA XIST in NSCLC, these miRNAs share numerous common targets and experimental evidence suggests that they synergistically regulate the cell-fate regulation of NSCLC. LncRNA XIST is known to sponge miR-449a and miR-34a, however, the regulatory network connecting all these non-coding RNAs is still unknown. Here we propose a Boolean regulatory network for the G1/S cell cycle checkpoint in NSCLC contemplating the involvement of these non-coding RNAs. Model verification was conducted by comparison with experimental knowledge from NSCLC showing good agreement. The results suggest that miR-449a regulates miR-16 and p21 activity by targeting HDAC1, c-Myc, and the lncRNA XIST. Furthermore, our circuit perturbation simulations show that five circuits are involved in cell fate determination between senescence and apoptosis. The model thus allows pinpointing the direct cell fate mechanisms of NSCLC. Therefore, our results support that lncRNA XIST is an attractive target of drug development in tumor growth and aggressive proliferation of NSCLC, and promising results can be achieved through tumor suppressor miRNAs.

As we mentioned above, the present study is based on our previous study [17]. Here we added a new layer of complexity that we explain below. The presented logical rules that govern the nodes in our model are suggested based on the biochemical literature demonstrated in Table S1.
DNA damage can activated miR-449a expression as observed by Mao et al. [21]. In this Boolean model, we added miR-449a, lncRNA XIST, YY1 and c-Met. The transcription factor Yin-Yang 1 (YY1) is overexpressed in NSCLC and directly activates lncRNA XIST [22]. YY1 directly inhibits miR-34a expression [23] and miR-449a [24]. YY1 is also a negative regulator of p53 activity [25]. Similarly, c-Met is a partner of the sub-family of receptors tyrosine kinases (RTKs) for hepatocyte growth factor (HGF), which is activated by the Myc proto-oncogene protein (c-Myc) and directly inhibits p53 expression [26], while miR-34a and miR-449a inhibit c-Met. In this way, miR-34a and miR-449a can enhanced p53 expression via knockdown YY1 and/or c-Met [27]. Once p53 is triggered in tumor cells in DNA damage response, it can induce p21 for cell cycle arrest for repair and/or senescence or it can induce PUMA and/or BAX for apoptosis. For more details about multistability predicted by our models [17].

Boolean Methods
The structure of the Boolean model of a regulatory network is predicated on the published biochemical knowledge associated with proteins, miRNAs, and lncRNA and their interactions (activatory or inhibitory), which are described by a directed graph. The variable's values are discrete taking 0 or 1 values. The logical functions are framed supported by the Boolean operators AND, OR, and NOT define the worth of individual nodes in keeping with the influence of its regulators. The Boolean model dynamics in its state-space may be described by a state transition graph, whose nodes define states of the system and whose edges characterize transitions amid these states. The state transition graph serves all the potential trajectories that a single initial state can drive to an ultimate state or attractor. Attractors that haven't any unreserved edges are named steady states or fixed points, while a collection of transitions trapped in an exceedingly fixed grouping of states form a cyclic state.The importance of the nodes within the system will be revised asynchronously or synchronously. Within the asynchronous scenario, nodes are picked randomly so revised, whereas within the synchronous scenario all nodes are updated at equivalent times. During this work, we only operated the asynchronous update strategy for realism because of its ability to provide non-deterministic modes [28]. Furthermore, the dynamics of a regulatory network are manipulated by negative and positive circuits (also referred to as feedback loops). Negative circuits can bring about oscillatory dynamics and positive ones to multi-stable states behavior. These circuits command the dynamics of molecular systems within the cells. To check their consequence on the dynamics, the Boolean technique permits in silico perturbations of nodes, i.e., forcing them to stay very fixed coequal to loss (LoF) or gain of function (GoF) experiments and having by interaction perturbations. Node or interaction perturbations can undermine the consequence of an operational circuit evidencing their precise role within the dynamics.
The GINsim model file is available in the Data File S1.

Results
The fixed points (also known as stable states) for the wild-type case (WT) and mutants are illustrated in Figure 2 together with experimental results and references. The input of the model is DNA damage which can be ON or OFF. Each line of the figure is a unique fixed point of the network. The first 3 lines correspond to the WT. The first fixed point of the WT is a proliferative state rising from the lack of DNA damage and is defined by activation of cell cycle regulators and inactivation of cell cycle inhibitors: E2F1, Cdk46-CycD, Cdk2-CycE c_Myc, Cdc25A, and lncRNA XIST. The next two fixed points of the WT for when DNA damage ON correlated to the bistability involving senescence and apoptosis. Bistability is a stochastic behavior of the model dynamics where two different fixed points are randomly selected from the same initial condition. The probabilities of these fixed points are not necessarily equal. This bistability is regulated by the p53-A/p53-K circuit [33], which was proposed in previous works [13,17,[33][34][35]. The remaining phenotypes in Figure 2 follow from mutations that will be discussed below. , whereas colored cells such as: purple, green, red, light blue, dark blue and black denote activation (ON). The stable states were pinpointed for different networking scenarios: WT, miR-449a KO, miR-449a E1, miR-34a KO, miR-34a E1, miR-16 KO, miR-16 E1, miR-34a E1 along with miR-449a KO, miR-449a E1 along with miR-34a KO. Cells in dark blue color and black color that describe the model predictions and indicated by the question mark such as mention previously, in this way, miR-34a Vs lncRNA XIST, Negative correlation between miR-34a and lncRNA XIST in DNA damage response, miR-34a E1, and lncRNA XIST E1. Red arrows designate down-regulation, whereas yellow arrows convey up-regulation. The left-most queue exhibits the input DNA Damage and the right-most queues demonstrate the model developments: proliferation, apoptosis, and senescence. Per line symbolizes a unique fixed point be compatible with the input.
In Table 1 we pointed the experiments used as influential references to our study (for more details see Figure S1. The model was fine-tuned to come up with the knowledge provided by each study (see section Methods for more details about the biochemical information and also the network construction). Some studies focused only on the interaction between the lncRNA XIST and one miRNA, while other studies considered more interacting miRNAs with lncRNA XIST (Table 1 and Figure 2). Here we resort to many of those studies. miR-449a E1 Activating Senescence and Apoptosis A549, H460 [14] lncRNA XIST KO Induces Senescence and Apoptosis A549, H460 [14] miR-16 Vs lncRNA XIST in DNA Damage Response Negative correlation A549 [15] miR-16 E1 Triggers senescence and apoptosis A549 [15] lncRNA XIST KO Induces Senescence and Apoptosis A549 [15] First, we interrogated how perturbations of the network connecting the non-coding RNAs affect cell fate (results are presented in Figure 2). To do that, we considered the experiments in Table 1 that focused on the lncRNA XIST and just one or more of the other miRNAs. For example, the investigation of the particular regulatory role of a single miRNA was conducted making an in-silico knockout (KO) of the others. We confirmed that different perturbations of miR-449a and miR-34a induce different phenotypes as observed in experiments using A549 and H460 cells by Luo et al. [38] and He et al. [36]. Other perturbations correspond to predictions, as they were not observed experimentally. For example, the combined perturbation of miR-449a overexpression (E1) with miR-34a KO, or the inverse perturbation, induced only apoptosis. However, the overexpression of both nodes generates a bistable state of senescence and apoptosis. This result is experimentally supported by the fact that both miR-449a and miR-34a directly inhibit the expressions of c-Myc, HDAC1, and Bcl2, which initiates the activation of p16 and p21, leading to cycle arrest. See Figure 2.
GINsim allows us to recognize functional circuits within the network, i.e., circuits that actively sustain the dynamics of network (see Section 2). We found 27 positive and 14 negative circuits (Table 2). More details about the functional circuits of the Boolean model can be found in Table S2. However, we elected only 25 of those circuits containing an utmost of four components, a number of them were already studied experimentally. Only twelve, out of these 25, involving miR-34a, miR-16 and miR-449a and/or lncRNA XIST are unknown. The biochemical interactions explaining these 12 positive circuits are well-known within the literature for NSCLC cells (Table 3), however their performance in the G1/S checkpoint is yet to be determined. Then, as observed experimentally NSCLC, we investigated through in-silico perturbations whether these circuits were important to manage the WT-case. We executed edge perturbations on the circuits (see Table 4). This sort of perturbation is challenging to be obtained experimentally, nonetheless, they are truly valuable to specify circuit significance. First, we select an edge and eliminate it. For example, for the circuit miR-34a/YY1/lncRNA XIST, we removed the interaction between miR-34a and lncRNA XIST. We obtained in this case as output states, senescence and apoptosis and then, removing the interaction between YY1 and lncRNA XIST, we got two apoptotic and one senescence states. And at last, removing the interaction between lncRNA XIST and miR-34a, we obtained senescence and apoptosis. Using this technique, we perturbed all 12 circuits. As can be seen in Table 4, we found that the interaction between YY1 and lncRNA XIST (belonging to the circuits miR-34a/YY1/lncRNA XIST, YY1/lncRNA XIST/miR-34a/c-Myc and miR-449a/c-Myc/YY1/lncRNA XIST) is responsible for generating apoptotic and one senescent state. Whereas, the interaction between Rb and c-Myc in circuits CDK4-6-CycD/RB/c-Myc/miR-16 and CDK2-CycE/RB/c-Myc/miR-16 addresses for two apoptotic, one senescent and one undefined state. In this way, we determined that only five, out of 12 circuits, could produce multistability ( Table 4, highlighted in bold fonts). Table 2. Functional circuits in the network and experimental observations. Circuits for which no experimental data were found are indicated by question marks (?) and highlighted in bold fonts.

Number
Circuits Reference

Discussion
It is well recognized that RNAs interaction networks (lncRNA-miRNA-mRNA) play an essential role in regulatory mechanisms and cancer [1,2]. In this context, lncRNA XIST has been studied in several cancers including NSCLC. It has been observed that lncRNA applies its biological influences on cancer by sponging or indirect interactions with XIST miRNAs. MiR-34a, miR-449a, and miR-16 regulate cell cycle progression by targeting CDK4,6-CyclinD and CDK2-Cyclin E complexes and the lncRNA XIST sponges these miRNAs. In the current study, based on a previously published work [17], we presented a Boolean model of the plausible network connecting lncRNA XIST/miR-34a/miR-449a and miR-16 in NSCLC (see Figure 1).
Lately, Zhang et al. [14] and Zhou et al. [15] demonstrated that lncRNA XIST can regulate proliferation and tumor formation by targeting miR-449a/miR-16 in NSCLC. On the other hand, Liu et al. [16], pointed the role of lncRNA XIST/miR-34a axis in thyroid cancer. They further demonstrated that lncRNA XIST was significantly up-regulated and miR-34a was under expressed in Thyroid cancer [16]. Moreover, lncRNA XIST knockdown inhibits cell proliferation in-vivo and in-vitro [16] in thyroid cancer cell lines. Interestingly, miR-34a is downregulated in NSCLC [36]. However, no other investigation in the literature analyzed the role of the miR-34a/lncRNA XIST axis on NSCLC. Thus, this finding detected in Liu et al. [16] study permits us to predict the mechanisms identified for miR-34a/lncRNA XIST in thyroid cancer cells that might also be observed in NSCLC cells.
Thus, to advance our understanding of the synergistic coordination of cell fate determination in NSCLC cells within a non-coding RNA interaction framework, a mechanistic understanding of this development is needed. Therefore, a Boolean model was constructed that observed in the experimental outcomes regarding the action of these molecules at the G1/S checkpoint in NSCLC cells. In the loss of DNA damage, the model defines only a unique fixed point represented by a proliferative phenotype that is expected. In contrast, in the occurrence of DNA damage, two fixed points (p53-dependent cell fates such as senescence or apoptosis) were obtained. Certainly, in NSCLC both phenotypes are recognized in response to DNA damage. In addition, the Boolean model was investigated through GoF and LoF perturbations of its variables corresponding to the experimental outcomes (see Table 1) from multiple experimental studies [14,15,[36][37][38][39][40] which examined the interaction between lncRNA XIST and miRNAs (miR-449a, miR-16, and miR-34a) in NSCLC cells. Previously, we have shown that miR-34a can regulate miR-16 expression by targeting HDAC1 and c-Myc [17]. Remarkable, we detected that miR-449a can regulate miR-16 as well as p21 expression through HDAC1, c-Myc, and lncRNA XIST. Additionally, we found that miR-34a can control miR-16 and p21 expression via lncRNA XIST knockdown (for more details, see Figures 2 and 3). In addition, in gene regulatory networks (GRN), biological circuits play a vital role. Certainly, GRN is a mixture of various circuits that can interpret the dynamics of a biological system [62,63]. In this regard, we found twelve new positive circuits (see Table 2), which by perturbation analysis suggested that only five of these circuits affect (see Table 4, highlighted in bold fonts) cell fate decisions. We also provided evidence of these circuits' presence according to NSCLC (for more details see Table 3). DNA double-strand break (DSB) can be caused by the radiomimetic chemicals or reactive oxygen species (ROS) or radiation, which triggers autophosphorylation of ATM at serine 1981 by initiating its kinase activity [64]. Downstream phosphorylation on the ATM pathway drives the activation of p53 in response to DNA damage. ATM directly induces Mdm2 phosphorylation, thereby reducing the activity of Mdm2 and resulting in strong induction of p53. Activated p53 induces Mdm2 which inhibits p53 by ubiquitination [49]. Within the model, based on the different phosphorylation events that p53 controls, it is described by two variables: p53-A and p53-K. p53-A represents p53 phosphorylated at serine 15 and serine 20 that induce p21 expression for the cell cycle arrest or senescence, whereas p53-K specifies further phosphorylation at serine 46 that modulates Bax and/or PUMA expression for the apoptosis (see [33] for more information on p53-A/p53-K function). Activated ATM/p53 can induce miR-34a [50,65] and miR-16 [66,67] expression in DNA damage response. On the other hand, miR-449a induced through the DNA damage as suggested by Mao et al. [21]. HDAC1, lncRNA XIST and C-Myc are well-known inhibitors of p21 [68][69][70] and miR-16 [15,71,72]. Interestingly, HDAC1, lncRNA XIST, and C-Myc are directly targeted by the miR-449a [14,21,73] and miR-34a [16,36,74]. Consequently, miR-449a can regulate p21 expression and miR-16 activity through the suppression of HDAC1, lncRNA XIST, and C-Myc. Once p21 expression and/or miR-16 is activated, it can trigger senescence and/or apoptosis in NSCLC at the G1/S checkpoint. Arrowhead edge represents activation. Whereas, a hammer-head edge represents suppression, respectively. Detailed molecular mechanism of DNA damage of the model see our previous publish study [17] and also the latest reviewed article by Huang et al. [75].

Conclusions
In summary, we proposed an interaction network connecting several non-coding RNAs and predicted that miR-449a can control miR-16 and p21 expression through HDAC1, c-Myc, and lncRNA XIST. Our approach may contribute to suggest alternative therapeutic strategies for cancer by targeting non-coding RNAs embedded in checkpoint regulatory networks.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/biology11040480/s1, Figure S1: The experiments used as major references to our study, Data File S1: The GINsim 3.0.0b code of the model. Table S1: Official names of the model elements and logical rules. Table S2: Functional circuits of the Boolean model.