p53/E2F1/miR-25 axis regulates apoptosis induction in glioblastoma cells: a qualitative model

p53 is an anti-cancer protein for inducing cell cycle arrest and apoptosis. In glioblastoma multiform (GBM), p53 is able to induce apoptosis via inhibition of its negative regulator Mdm2. Experimental studies have shown that microRNA-25 (miR-25) can repress Mdm2 expression and, in turn, stabilizes p53 to induce the G1/S checkpoint in GBM cells. miR-25 expression is regulated by the cell cycle inducer molecule E2F1, which has been reported to induce apoptosis when overexpressed in GBM. However, the way p53, E2F1, and miR-25 coordinately regulate apoptosis induction in GBM is still obscure in the literature. In this study, we propose a logical model contemplating the regulatory influence of miR-25 and its regulator E2F1 on cell fate decision. Through in silico results for the wild-type case, we observed that miR-25 may stabilize p53 expression through Mdm2 inhibition inducing a G1/S checkpoint arrest or apoptosis in cells overexpressing E2F1. The predicted probabilities of our model are in good agreement with published experimental data. Moreover, we show that miR-25-induced p53 stabilization might contribute to apoptosis induction in GBM cells. These findings highlight some unrecognized mechanisms that may guide to alternative ideas for GBM therapeutic strategies.


Introduction
p53 pathway is the most important process involved in the interception of cancer initiation (Vousden and Lane 2007). The expression levels of p53 are highly related to cell fate decision governing the balance between survival and death. The levels of p53 are tightly controlled by its feedback loop with Mdm2 ( Bar-Or et al 2000). While Mdm2 negatively regulates p53 through Mdm2-mediated ubiquitination, p53 induces the expression of Mdm2. The inactivation of the p53 is usually associated with the amplification of Mdm2 expression in cancer (Forte et al 2019). Targeting of Mdm2 is an important cancer treatment strategy due to potential reactivation of p53 function via the induction of the p53-responsive phenotypes such as cell cycle arrest and apoptosis. Thus, the study of new mechanisms associated with Mdm2 repression helps to improve strategies of cancer treatment.
Glioblastoma multiforme (GBM) is an aggressive brain tumor that presents high proliferation and dissemination rates (Nørøxe et al 2016), and radiotherapy is one of the current standard treatments (Mann et al 2018). However, the high resistance to radiation of GBM-derived cell lines due to failure of radiation-induced primary apoptosis is one barrier against GBM cell elimination. The attenuation of p53 function by Mdm2 can lead to the resistance of GBM cells to radiation treatment (Afshar et al 2006;Costa et al 2013). In this context, Villalonga-Planells and colleagues showed the p53 activation by Mdm2 inhibition can induce apoptosis in GBM cells via induction of PUMA (Villalonga-Planells et al 2011). In addition, the study of Daniele and colleagues explained that the targeting of AKT can activate p53 expression. They further demonstrated that the downregulation of AKT decreases Mdm2 expression and increases p53 expression (Daniele et al 2015). Thus, Figure 1. Regulatory network for the G1/S checkpoint of GBM cells. Nutlin_3a and E2F1_inductionin dark gray are the inputs of the network. Nodes in light gray represent molecules. Rectangular nodes represent components with two states (active or inactive) and the elliptical node denotes a component with more than two states. Black arrows denote activations and red hammerhead connectors represent inhibitions. these data suggest that the reactivation of p53 by repression of Mdm2 may be a successful therapeutic strategy for GBM treatment.
MicroRNAs (miRNAs) have an important role in gene expression regulation. They are small non-coding RNAs of approximately 20-25 nucleotides in length which function as post-transcriptional gene expression regulators. MiRNAs can play an important role in the molecular dynamics of GBM (Suh et al 2012). Recent reports have shown that one miRNA in particular, microRNA-25 (miR-25), is involved in several biological processes such as proliferation and cell cycle arrest in GBM (Suh et al 2012;Peng et al 2015). miR-25 is normally overexpressed in GBM cells (Peng et al 2015) and its activity has been reported to induce abnormal cell cycle progression for inhibiting the cyclin-dependent kinase inhibitor 1C (also known as p57 or p57Kip2) in GBM cells (Zhang et al 2015). However, Suh and colleagues showed that overexpression of miR-25 can stabilize p53 expression by inhibiting Mdm2 and thus inducing G1/S checkpoint in GBM-derived U87 cell line (Suh et al 2012). These results suggest that miR-25 can reactivate p53 function, being a potential target molecule for GBM treatment.
In addition to p53, other transcription factors (TFs) may present an important role in cell fate decision. Among them, the E2F family is known for regulating cell cycle progression. Specifically, E2F1 has been demonstrated to induce the expression of many genes including miR-25 (Suh et al 2012). In GBM, besides its functional role in cell proliferation, E2F1 overexpression can induce apoptosis via caspase-8 (Shu et al 2000). Moreover, E2F1 can also contribute to cell fate decision inducing p53 activation through target molecules that are able to inhibit Mdm2 expression (Inoue et al 2016). miR-25 has been reported as the link between E2F1 and p53 which performs this function (Suh et al 2012). However, how E2F1 and p53 coordinately regulate apoptosis induction via miR-25 expression in GBM is still unclear.
Motivated by the above considerations, we propose in this paper a logical model of G1/S checkpoint regulation contemplating the regulatory interactions among miR-25, p53, and E2F1 (see figure 1) to clarify the mechanisms of apoptosis induction in GBM. Based on experimental studies, we consider two inputs in the network: nutlin-3a (p53 activator and DNA-damage inducer) and E2F1 induction by adenoviral infection.
Integration of computational and experimental approaches has contributed to the development of new treatment strategies for glioblastoma (Ozdemir-Kaynak et al 2018). Specifically, the development of a molecular regulatory network into a formal framework of a computational-experimental integration process is being recognized as a valuable approach to study cell-fate decisions and many other biological processes (Le Novere

Logical formalism used to construct the regulatory network
The construction of the logical model of a regulatory network (Albert and Thakar 2014) is based on the published biochemical information associated with proteins and miRNAs and their interactions (activatory or inhibitory), which are characterized by a directed graph. The variables values are discrete taking whole number values. A logical function defines the values of each node according to the influence of its regulators through the logical operators AND, OR, and NOT. A model component is considered active (a value greater than 0) due to its regulation level to be sufficient to exert an effect, otherwise it is considered inactive (value 0) (Abou-Jaoudé et al 2016).

Dynamics analysis
The dynamics of a logical model in its state space can be represented by a state transition graph, whose nodes represent states of the network and whose edges represent transitions between these states. The state transition graph represents all the possible trajectories that one initial state can lead to a final state or attractor. Attractors which have no outgoing edges are called steady states, whereas a set of transitions trapped in a fixed group of states compose a cyclic state. Node values in the network can be updated asynchronously or synchronously. In the asynchronous scheme, nodes are selected randomly and then updated, whereas in the synchronous scheme all nodes are updated at the same time. In this work, we used only the asynchronous updating scheme due to its potential to present non-deterministic behavior (Abou-Jaoudé et al 2016). In this study, a Monte Carlo algorithm was used to estimate the attractor reachability using 'exact exit probabilities' and 100 000 runs (Mendes et al 2018). This method allows estimating the reachability probabilities of stable states by avoiding the exploration of states previously visited and appropriately associating new probabilities with state transitions. For more details about the algorithm, see (Mendes et al 2018). Furthermore, the dynamics of a regulatory network is influenced by negative and positive circuits (also known as feedback loops) (Thieffry 2007). Negative circuits can give rise to oscillatory dynamics and positive ones to multi-stationary behavior (Thieffry 2007). These circuits control the dynamics of molecular networks inside cells. To test their effect on the dynamics, the logical approach allows in silico perturbations of nodes, i.e. force them to remain in a fixed state corresponding to loss (LoF) or gain of function (GoF) experiments and also by interaction perturbations (Li and Gao 2019). The perturbations can disable the influence of the functional circuits (multistability or oscillations), which elucidate their specific role in the dynamics.
The model analysis was performed using the tool GINsim 3.0.0b (Naldi et al 2018). The logical rules and official names of the model elements can be found in S1-table (https://stacks.iop.org/JPCOMPLEX/01/035001/ mmedia). The model file is available in S1-file.

Proposed regulatory network of the G1/S checkpoint mechanism in GBM cells
We propose a G1/S cell cycle checkpoint regulatory network for GBM cells contemplating the involvement of miR-25. The proposed network is based on established biochemical information with some reasonable assumptions due to the absence of experimental data. These assumptions will be described in the next paragraph. The network encompasses 22 nodes (proteins and miRNA) and 49 direct interactions among them, and it contemplates two inputs: nutlin-3a and E2F1 induction, representing stimuli (see figure 1). Light gray nodes represent proteins and miRNA, black arrows denote activations, and red hammerhead connectors denote inhibitions.
Nutlin-3a is able to inhibit p53-Mdm2 interaction and activate p53 signaling in cancer cells (Villalonga-Planells et al 2011). Furthermore, this molecule triggers DNA damage response in the cell which, in turn, can help restore p53 function (Verma et al 2010;Costa et al 2013). Both mechanisms were considered in the proposed network. Regarding the molecular processes involved, the activation of the G1/S and G2/M cell cycle checkpoints in response to p53 activation has well-known mechanisms (Sancar et al 2004;Medema and Macůrek, 2012;Gupta et al 2018). DNA double-strand breaks activate the ATM pathway, while single-strand breaks SSB activate the ATR pathway. The official biochemical names of the model elements and the corresponding logical rules can be found in the S1-table. The phosphorylation cascade induced by these pathways leads to the activation of p53 (see review by (Gobbini et al 2013)). ATM and ATR phosphorylate Mdm2 activating p53 expression. p53 can induce Mdm2 which, in turn, degrades p53 forming a negative feedback loop (Bar-Or et al 2000). Activation of p53 can coordinate cell cycle arrest and apoptosis by its phosphorylation at distinct sites (Zhang et al 2011). The two forms of phosphorylated p53 are called p53_A and p53_K, which define the p53 function in the network (Zhang et al 2011). Hence, p53_an induces molecules related to cell cycle arrest as p21 (Zhang et al 2011), while p53_K activates pro-apoptotic proteins as PUMA (Nakano and Vousden 2001). The transition between these forms is controlled by a positive circuit between p53_A/p53_K. For more details on this positive circuit, see our previous publications (Gupta et al 2018;Gupta et al 2020b). P57Kip2 is activated by ATM (Lee 2007). Both cell cycle arrest inducers p21 and p57Kip2 inhibit the activity of the cell cycle promoters cyclin-dependent kinase 4 and 6-cyclin D and cyclin-dependent kinase 2-cyclin E, abrogating cell proliferation and arresting the cell cycle at G1/S (Donjerkovic and Scott 2000). Cdc25A is a phosphatase required to activate the protein complexes: cyclin-dependent kinase 4 and 6-cyclin D and cyclin-dependent kinase 2-cyclin E (Donzelli and Draetta 2003). Activation of the CDK-cyclin complexes phosphorylates the retinoblastoma protein (RB) releasing E2F1 (Sherr and Roberts 1999). RB inhibits Mdm2 to stabilize p53, while Mdm2 inhibits RB, forming a positive circuit between RB/Mdm2 (Yap et al 1999). Activation of E2F1 can induce proliferation or apoptosis according to its level of accumulation in cells as suggested by the studies of Afshar et al 2006 andShu et al 2000. Normal levels induce proliferation, whereas high levels trigger apoptosis. The two forms associated with the accumulation of E2F1 in the model are called E2F1_P and E2F1_K. Therefore, E2F1_P induces molecules associated with proliferation as CDKCycE, while E2F1_K enhances the expression of the pro-apoptotic protein caspase_8. Both E2F1_P and E2F1_K are assumed to activate miR-25 expression. In the model, E2F1_P induces the intermediary level of miR-25 (miR-25:1) found in GBM cells, whereas E2F1_K activates the highest level of miR-25 (miR-25:2). The transition between the two forms of E2F1 in the model is controlled by a double inhibition (a positive circuit) due to their opposite functions in the model; this approach has already been used in other computational models (Steinway et al 2014;Yan et al 2014). ATM and ATR can activate the apoptotic role of E2F1 (Carcagno et al 2009). As such, we assumed that ATM and ATR trigger E2F1_K activation. It was considered in the model that p53 can present an inhibitory role in E2F1 expression, the response of which is observed in GBM cells (Suh et al 2012). The input E2F1_induction corresponds to external stimuli that induce E2F1 expression in the cell.

Coherence between the model and experimental perturbations related to the G1/S checkpoint in GBM
In order to observe the coherence between the model and experimental perturbations related to the G1/S mechanism in GBM, we compared the outcome of the perturbed model to the observed phenotypes of mutants found in the literature (table 1). Although other molecules not included in the model can contribute to the G1/S checkpoint in GBM, the model consistency in terms of its dynamics under perturbations indicates that the G1/S arrest in GBM cells is well characterized by the molecules considered in the network (Annovazzi et al 2015). For further details about model validation, see S1 figure.

p53 can contribute to apoptosis induction in cells overexpressing E2F1 via miR-25
To analyze the dynamics of the model under nutlin-3a and E2F1 stimuli, we performed wild-type simulations (no perturbations are considered) identifying all possible attractors of the dynamics (figure 2). Thus, the approach identified ten steady states which are mainly distinguishable by the state of activation of the inputs nutlin-3a and E2F1_induction (figure 2).
For inactive inputs, we observed a state that can be interpreted as a proliferative phenotype, since all cell cycle promoters, i.e. CDKCycD, CDKCycE, Cdc25A, and E2F1_P, are activated, consistent with experimental observations (Donzelli and Draetta 2003;Zhang et al 2005;Xiao et al 2013). On the other hand, active inputs trigger multistability which is associated with non-deterministic behavior (Silveira and Mombach 2020; Gupta et al 2020a). Nutlin-3a ON induces a bistable switch related to p53 activation as p53 LoF abrogates this state. In the model, G1_Arrest corresponds to the activation of p53_A, p21, p57Kip2, and RB consistent with experimental data (He et al 2005;Lee 2007;Lahav 2008;Biasoli et al 2013). The induction of p53_K and PUMA enhances apoptosis in the model (Villalonga -Planells et al 2011).
For the input condition nutlin_3a OFF and E2F1_induction ON in figure 3, three states were observed (tristable case), which correspond to a proliferative and two apoptotic states. Activation of E2F1_P preserved the proliferative state, which is consistent with experimental data (Xiao et al 2013). In contrast, E2F1_K induction activates caspase_8 and apoptosis, which is also consistent with experimental observations (Shu et al 2000). Also, E2F1_K activated high levels of miR-25 which, in turn, stabilized p53 activity by inhibition of Mdm2. Once p53 is activated, the dynamics is the same observed for nutlin-3a ON. For the cases where the outputs, cell cycle arrest, and apoptosis are simultaneously active, the states were interpreted as apoptotic, since they are induced during cell cycle arrest (Medema and Macůrek, 2012;Gupta et al 2018).
When both inputs (nutlin-3a and E2F1_induction) are active, the model presents multistability composed of four states. Three of which are apoptotic states activated by PUMA and/or caspase_8 in response to p53_K  and/or E2F1_K, respectively. The only state associated with G1_Arrest is induced by p53_A. These results suggest that p53 can contribute to apoptosis induction in GBM cells overexpressing E2F1 via miR-25. The multistability observed in figure 2 is controlled by the functional circuits in the model. By analyzing with GINsim, we observed the functionality of 7 circuits (table 2). To evaluate their specific role in the model dynamics, we perturbed (LoF) the interactions that compose the circuit. We verified that E2F1_K/E2F1_P and p53_A/p53_K are responsible for the multistability of the attractors observed in figure 2 (see S2-figure). For the input condition E2F1_induction ON and nutlin3-a OFF, E2F1_K/E2F1_P triggers bistability through a bistable switch between these components, i.e. E2F1_P ON and E2F1_K OFF switch to E2F1_P OFF and E2F1_K ON. In the state where E2F1_K is ON, p53 is activated due to the activation of miR-25. This mechanism triggers the activation of p53_A/p53_K circuit inducing another bistable switch. Thus, E2F1_P/E2F1_K and p53_A/p53_K circuits trigger a tristability with E2F1_induction ON.
For nutlin_3a ON and E2F1_induction OFF, p53_A/p53_K triggers the bistability, whereas the condition nutlin_3a ON and E2F1_induction ON, both p53_A/p53_K and E2F1_P/E2F1_K induce a multistabiliy involving four states. The remaining circuits are functional during the transient dynamics of the model, i.e.  their functionalities do not affect directly the attractors of the system. Collectively, the results of the circuits indicate that p53 and E2F1 present a pivotal role in cell fate decision of the model.

Reachability probabilities of the multi-stable cases support the role of p53 on E2F1-induced apoptosis via miR-25 activation
To confirm the influence of p53 on E2F1-induced apoptosis in the model, we calculated the probabilities of reachability of the apoptotic states of each multi-stable case shown in figure 2 through a Monte Carlo algorithm with 100 000 runs. We analyzed the cases for nutlin-3a ON (bistability), E2F1_induction ON (tristability), and both ON (multi-stability involving four states). Thus, by using the proliferative state as an initial condition, due to GBM having high proliferation rates (Nørøxe et al 2016), we observed that the model is qualitatively coherent with experimental data ( figure 3(A)). For simulations where the input E2F1_induction is ON, the states of E2F1_P and E2F1_K, which are regulated by the input, were not restricted to a particular value in order to better explore the space state. The histograms of the probabilities were splitted into separate ones according to the apoptotic inducer in each state. As such, the probability of each apoptotic state induced by p53 and/or E2F1 was highlighted, stressing the influence of p53 on apoptosis in this multi-stable case. Thus, the results suggest that the combination of p53 and E2F1 may present a synergistic effect on apoptosis induction in GBM cells.
Probabilities were compared with experiments using the GBM cell line U87 having wild-type p53. The experimental approaches are related to cells treated with nutlin-3a (DNA-damage inducer) (Verma et al 2010;Costa et al 2013), E2F1 infection (Shu et al 2000), and a combination of ionizing radiation (DNA-damage inducer) and E2F1 infection (Shu et al 2000). The concentrations or doses used to treat the cells were 10 μM, multiplicity of E2F1 infection of 100 (the adenoviral infection corresponds to the use of a recombinant adenoviral vector expressing E2F1 in U87 cells), and the combination of a multiplicity of E2F1 infection of 100 and 10 Gy of radiation. Experimental results using IR were compared with nutlin-3a activation due to their ability to induce DNA-damage. Data on a combination of E2F1 induction and nutlin-3a treatment in U87 cells were not found in the literature.
In order to observe whether p53 stabilization via E2F1-induced miR-25 can affect the probabilities of the wild-type case for nutlin-3a ON ( figure 3(B)), we performed the GoF of the highest level of miR-25. This perturbation was performed under this condition to exclude the effect of the E2F1_K/caspase-8 axis (figure 3(B); top-right subnetwork) on apoptosis induction. Thus, we observe that p53 and miR-25 can coregulate cell fate. These results support our findings related to the possible synergistic effect of E2F1 and p53 on the apoptosis induction in GBM cells.

Proposed experimental design
The concept of experimental design is motivated by the studies of Moradimotlagh et al (Moradimotlagh et al 2020) and Daniele and colleagues (Daniele et al 2015), which show that the targeting of Mdm2 can be a promising treatment of GBM. These studies led us to introduce a new experimental design based on our model, which consists of using a combination of E2F1 overexpression and an inhibitor of Mdm2 activity such as nutlin-3a. E2F1 and the Mdm2 inhibitor may act synergistically to stabilize p53 expression rendering the cells more susceptible to apoptosis induction. Hence, the possible outcome of this experiment would be high rates of apoptotic cells induced by E2F1 and p53.

Discussion
Logical formalism was used in this study to investigate the influence of miR-25 and E2F1 on p53 stabilization in GBM cells during nutlin-3a stimulus. The regulatory mechanisms of the G1/S checkpoint-induced p53 involve the inhibition of E2F1, Cdc25A, cyclin E/Cdk2, and cyclin D/Cdk4 (Annovazzi et al 2015). According to experimental studies with GBM cells, p21 and p57Kip2 have an important inhibitory role of these molecules through direct or indirect effect (Tsugu et al 2000;Abbas and Dutta 2009). The upstream regulators of both p21 and p57Kip2 are p53 and ATM, respectively, which are DNA damage transducers. Cancer cells are proliferative due to the enhanced activity of cyclins, E2F1, Cdc25A, and other molecules (Donzelli and Draetta 2003;Zhang et al 2005;Xiao et al 2013). Under an initial condition such as DNA damage, the checkpoint pathway is activated arresting cells at the G1/S checkpoint or inducing apoptosis via p53 stabilization. In GBM cells, miR-25 is able to inhibit the transcriptional activity of p57Kip2 (Zhang et al 2015) and Mdm2 (Suh et al 2012), and its expression is controlled by E2F1 (Brosh et al 2008), in our model miR-25 inhibits p57Kip2 at both levels (miR-25:1 and miR-25:2). Whereas, miR-25 inhibits Mdm2 at high level (miR-25:2), for more details see S1table. The effect of E2F1 on apoptosis was also contemplated by the model as experimental studies report that its overexpression induces caspase-8, a potent apoptosis inducer in GBM cells (Afshar et al 2006). It is important to point out that experimental data report that GBM cells can present resistance to radiation due to the failure of p53 to induce its pro-apoptotic target BAX (Shu et al 2000). However, Villalonga-Planells and colleagues showed the p53 activation via nutlin-3a in GBM cells was able to induce apoptosis via induction of PUMA, another p53-induced pro-apoptotic molecule (Villalonga-Planells et al 2011), which was contemplated in the model. Thus, the proposed network considers the mechanisms related to cycle progression and G1/S checkpoint that are characteristic of GBM cells.
MiR-25 is a well-known negative regulator of Mdm2 and p53 (Sárközy et al 2018). Experimental studies have shown the reduction of p53 concentration caused by miR-25 via direct inhibition of p53 3'UTR. Consequently, miR-25 can act as a cancer promoter molecule. Nevertheless, Suh and colleagues showed that in the U87 cell line miR-25 is only able to repress Mdm2, while p53 levels are not reduced (Suh et al 2012). Thus, the present study may be inaccurate for other cell lines. The presence of wild-type p53 in these cells is essential for the purpose of our model. Moreover, our approach is one version of the mechanisms that can occur in the U87 cells, so this does not exclude the possibility that other proposals might also be coherent with experimental data.
Mdm2 is well known for controlling the balance between survival and death via regulation of p53. Higher levels of Mdm2 maintain the downregulation of p53, whereas lower levels trigger p53 accumulation in the cell allowing the induction of p53-responsive cell fates. As such, the targeting of Mdm2 to activate p53 is a potential strategy for cancer treatment. According to our study, miR-25 might be able to stabilize activated p53 expression inducing G1/S arrest and apoptosis. Indeed, experimental studies have shown that p53 stabilization can induce these phenotypes in GBM cells (Villalonga-Planells et al 2011;Suh et al 2012). It is important to point out that in human cancers with strongly downregulated p53, miR-25 might result in a failure to induce enough p53 activity to inhibit tumor development. Therefore, the combination of miR-25 expression with factors that can help stabilize p53 may render the cells more susceptible to p53-dependent responses.
The activation of miR-25 in both G1/S arrest and apoptosis phenotypes suggests that this molecule may contribute to cell fate decision in GBM. According to our study, overexpression of E2F1 induces miR-25 which stabilizes p53 inducing apoptosis. According to experimental data, E2F1 and p53 are able to induce apoptosis in GBM (Afshar et al 2006;Villalonga-Planells et al 2011). These results suggest that these TFs can contribute to the apoptotic phenotype in GBM. Indeed, studies have shown that E2F1 can trigger p53 accumulation inducing apoptosis in other cell lines (Rogoff et al 2002;Tanaka et al 2005). Thus, the main prediction of the model p53 stabilization by activation of miR-25 induced by E2F1 overexpression can provide a new regulatory role of miR-25 in cells overexpressing E2F1 associated with the apoptosis induction. To test experimentally our findings, our study proposes the analysis of E2F1 overexpression observing its effect on miR-25 and p53 expressions in U87 cells treated with Mdm2 inhibitors such as nutlin-3a. Hence, it might be possible to verify a synergistic effect toward apoptosis which would be useful for a therapeutic strategy of GBM.