Breast cancer is one of the frequently diagnosed lethal malignancies affecting millions of women worldwide. The risk factors include both genetic and epigenetic abnormalities whereby the latter provides early prognostic biomarkers for breast cancer therapeutics (Chen, 2012; Sun et al., 2017). Furthermore, DNA methylation is the most prominent epigenetic marker that frequently influences the expression or silencing of genes involved in vital cellular activities such as cell proliferation, programmed cell death, and cell differentiation (Chen, 2012). DNA methylation is mediated by DNA methyltransferases (DNMTs), in which the DNMT1 isomer is explicitly responsible for the maintenance of methylation status during replication (S-phase) (Pouliot et al., 2015). Briefly, it transfers the methyl group (−CH3) on cytosine in CpG dinucleotide regions by keeping the CpG islands non-methylated that are preferably located in the proximal promoter region of a gene (Tang et al., 2009). On the contrary, the hypermethylation of promoter CpG islands of TSGs due to the overexpression of DNMT1 is a popular mechanism of gene silencing and a hallmark of cancer (Subramaniam et al., 2014; Pouliot et al., 2015).

DNMT1-regulated hypermethylation of the TSG RUNX3 promoter region is a new insight and an early event in breast cancer predominantly in TNBC (triple-negative breast cancer) (Widschwendter and Jones, 2002; Jung et al., 2007; Jiang et al., 2008; Subramaniam et al., 2009; Shin et al., 2016). Generally, RUNX3 is known to elicit its tumor-suppressive ability through major cancer signaling pathways including TGF-β, Wnt/β-catenin, and KRAS (Chen, 2012). RUNX3 either influences the downstream target of tumor suppressor signaling pathways or acts as an antagonist for oncogenic pathways to exert its antitumor activity. However, despite the evident role of RUNX3 in breast cancer and the relation of DNMT1 with RUNX3 gene, no study has explored the underlying epigenetic molecular events by which DNMT1 silences RUNX3 to promote breast cancer metastasis.

Furthermore, due to the pivotal implication of DNMT1 in various tumors (Li et al., 2019) and the reversible nature of its methylation activity, targeting DNMT1 through small modulators is a promising pharmacological intervention to revive the suppressed TSGs (such as RUNX3). To date, two DNMT inhibitors (i.e., azacytidine and decitabine) have been approved by FDA for the treatment of myeloid malignancies (Gnyszka, Jastrzębski, and Flis 2013). However, the potential side effects associated with these drugs limit their use in high-grade malignancies. Therefore, the development of more potent DNMT1 inhibitors acquiring potential binding pocket features is highly desirable to restore the suppressed TSGs (RUNX3) in cancer therapeutics. Previously, several in vivo and computational studies have targeted DNMT1 for therapeutic purposes (Yoo 2011; Yoo, Kim, and Robertson 2012; Mirza et al., 2013; Krishna et al., 2017); however, its impact on downstream target genes, particularly RUNX3, requires further elucidation.

Therefore, in the present study, a BRN of the DNMT1–RUNX3 signaling was developed to provide an insight into the epigenetic-mediated silencing of RUNX3 leading to cancer development and metastasis. For this purpose, qualitative modeling by the René Thomas formalism was applied in a SMBioNet tool utilizing the existing wet-laboratory data in the form of computation tree logic (CTL) for parameter estimation and verification through the model checking technique (Figure 1). The model trajectories were explored to identify the paths involved in the overexpression of DNMT1, activation of oncogene (c-myc), and suppression of tumor suppressor genes (RUNX3, p21, and p53) leading to cancer invasion or recovery (homeostasis). In addition, DNMT1 was advocated as a potential drug target and favored to be subjected to ML approaches to classify diverse DNMT1 modulators. Decision tree (DT) and neural network classifiers were built utilizing the data set of DNMT1 from the ChEMBL database to identify potential two-dimensional binding features essential for the modulation of DNMT1 activity (Figure 1).

The overall methodology used in the current study is provided in Figure 1.

A schematic model depicting the stepwise process of replication and methylation of TSG RUNX3 by DNMT1 and its subsequent influence on various cancer signaling pathways was established using the literature-driven information. Qualitative modeling was then performed to unfold the cellular events involved in the epigenetic, that is, DNMT1-mediated silencing of RUNX3 that leads to cancer invasion. The comprehensive literature-driven pathway of DNMT1–RUNX3 is elaborated in the Supplementary Material Section S3, Figure 2. Briefly, a series of stepwise events at the replication fork led by macroprotein complexes (DNMT1, UHRF1, HAUSP, Tip60, HDAC1, and PCNA) ensures the methylation of the newly synthesized daughter strand during the S-phase (Figure 2). As a result, the nascent RUNX3 strand acquires the normal status of methylation, whereas the promoter region remains hypomethylated to facilitate active transcription of RUNX3 gene. Subsequently, RUNX3 (transcription factor) functions as a tumor suppressor protein to combat cancer initiation and metastasis through various signaling pathways. The TSG RUNX3 elicits an antitumor activity by regulating the transcription of target genes (p21, c-myc, etc.) of the major cancer signaling cascades, such as transforming growth factor-beta (TGF-β), Wnt/β-catenin, and mitogenic KRAS pathway as explained in Figure 2.

From literature-driven pathways described in Figure 2, a qualitative biological regulatory network was generated based on preferred entities. DNMT1, RUNX3, p21, c-myc, p53, and MDM2 nodes were selected due to their well-established functionality in the biological system (Supplementary Material Section S3). As a result, the BRN (Figure 3A) composed of total six nodes and nine interactions exhibiting all the significant activation and inhibition relationships was obtained. Initially, negative feedback loops were observed from an inhibitory set of genes required by the system to generate the stable states. According to the interaction graph (Figure 3A), RUNX3 transactivates p21 to maintain the concentration of DNMT1 through a negative feedback loop. In addition, RUNX3 inhibits the onset of oncogene c-myc that might upregulate the expression level of DNMT1 in a positive manner. Interestingly, p53 also prevents upregulation of DNMT1 through the activation of p21 and inhibition of oncogene c-myc. Furthermore, the proposed network characterized oscillatory behavior of two regulatory loops involving DNMT1 and RUNX3: 1) a negative feedback loop between p21 and DNMT1 through RUNX3 and 2) a positive feedback loop between c-myc and DNMT1 through RUNX3 (Figure 3A). It is notable that the representation of the BRN as a weighted, directed graph (Figure 3A) was obtained using GENOTECH through the implementation of discrete formalism. Therefore, the interaction graph utilized the qualitative data only, that is, type of interaction (activation +, or inhibition −) among the nodes (genes), and the threshold value for each interaction. However, modeling the dynamic behavior of such complex systems that include both positive and negative loops requires the computation of accurate logical parameters. These parameters were generated in the SMBioNet software using the known experimental observations encoded as CTL formulas reported in Figure 3B. The first CTL observation searches for a state with high expression of DNMT1 and oncogenes leading to tumor invasion. The later CTL formula seeks for the stable state or homeostasis that has normal expression level of tumor suppressor genes (Figure 3B). These observations were coded as input in SMBioNet to generate the preferred set of logical parameters.

TNBC is the most aggressive subtype of breast cancer that lacks obvious treatment options due to the availability of limited information about definite biomarkers. Previously, aberrant epigenetic modifications have been implicated in breast cancer, highlighting RUNX3 as a promising prognostic biomarker (Wang et al., 2014). Moreover, several former studies have reported the depleted level of RUNX3 in breast cancer cell lines predominantly due to local hypermethylation at the proximal promoter region of TSG RUNX3, which is an early event in carcinogenesis (Chen, 2012; Lau et al., 2006; Chen et al., 2016). Therefore, the aberrant DNA methylation pattern is known as the hallmark of cancer epigenetics characterized by DNMT1, also known as the maintenance methyltransferase. So far, the individual role of RUNX3 as a TSG in breast cancer and DNMT1 as the methylation maintenance protein has been reported in various scientific studies (Jiang et al., 2008; Chen et al., 2016; Lau et al., 2006; Subramaniam et al., 2009). However, the underlined epigenetic-mediated RUNX3-shared signaling that instigates and disseminates the tumorigenesis remains largely elusive. This study provides an insight into the integrated DNMT1–RUNX3 signaling by taking into account various cancer-related significant upstream and downstream regulators (such as p21, c-myc, p53, and MDM2) and presents the DNMT1–RUNX3 signaling cascades as one consolidated network. To the best of authors’ knowledge, this study is one of its kind to illustrate and model the epigenetic-mediated silencing of TSG RUNX3 through the René Thomas framework modeling that has largely been a common practice of systems biology to investigate the dynamics of biological networks (Saeed et al., 2018).

Moreover, the model checking technique was adopted to build the interaction graph rendered as a state graph (Figure 6B) by utilizing the literature-driven information of the DNMT1–RUNX3 signaling (Figure 2). The similar approach of the model checking has been previously applied in different former studies such as parameter estimation through formal modeling (Ahmad et al., 2012; Bibi et al., 2017), tail resorption in tadpole metamorphosis (Khalis et al., 2009), and immunity control in bacteriophage lambda (Richard, Comet, and Bernot 2006). We utilized this technique in the SMBioNet tool to exhaustively explore the model space of the DNMT1–RUNX3 BRN for the estimation of precise parameters by encoding the existing wet-laboratory data in the form of computation tree logic (CTL) (Figure 3B). Once the computational verification with laboratory data is completed, the resultant model parameters and all its trajectories were further analyzed to understand how systematic evolution of the DNMT1–RUNX3 system takes place with time, which might lead the system to invade cancer metastasis (Figure 6B).

As a result, two important behaviors of the DNMT1–RUNX3 system were highlighted that shows the successive genetic events of the recovery cycle (Figure 5) and a bifurcation state leading to cancer invasion (Figure 6C). The trajectory (right side of Figure 6C) from our model plotted as a state graph articulates that the onset of oncogene c-myc infects the system and anticipates the upregulation of DNMT1 by downregulating the RUNX3 expression level. It generates a feedback loop where the consistent onset of oncogene c-myc accompanied by the persistent suppression of TSG RUNX3 contributes toward the overproduction of DNMT1. These findings are in agreement with the previous experimental studies, which reported the transcriptional upregulation of DNMT1 in TNBC due to the amplified c-myc expression level, emphasizing the importance of DNMT1 and c-myc relation in tumorigenesis (Wu et al., 2021).

On the contrary, the activation of RUNX3 (left side of Figure 6C) and subsequent TSGs (such as p21, p53) regulate the adequate level of DNMT1 to acquire homeostasis. In addition, network analysis emphasized the tumor-suppressive role of RUNX3 by demonstrating the higher probability of the infected system to recover to normal homeostasis in the presence of RUNX3 (left side of Figure 6C) than its likelihood to invade the tumorigenesis. This is in line with the former findings that have conversely related the reactivation of RUNX3 with a reduced potential of metastasis and invasiveness in breast cancer cells (Widschwendter and Jones, 2002). Although computational models cannot replace experiments, they are a step to demonstrate whether or not a proposed mechanism is sufficient to produce an observed phenomenon or an underlying assumption on the basis of their mathematical framework.

Moreover, the findings of our qualitative modeling suggested DNMT1 as a critical hub regulator of various oncogenic and tumor suppressor proteins that determine the ultimate status of the system. These outcomes advocate DNMT1 as a potential drug target in epigenetic cancer signaling through logic-based temporal evolution. Notably, the hypermethylation of TSGs by DNMT1 is a reversible process; therefore, RUNX3 expression could be restored using demethylating compounds. Consequently, designing specific small modulators that inhibit DNMT1 activity in order to restore the TSG RUNX3 could represent new clinical avenues for breast cancer therapeutics. Several pioneer studies have discussed the probable restoration of RUNX3 and a reduced carcinogenic potential in cancer cell lines when treated with demethylating drugs (Lau et al., 2006; Jung, Park, Young Kim, et al., 2007).

To date, only two FDA-approved drugs (azacytidine and decitabine) are available that target DNMT1 in addition to other DNMTs for the treatment of myelodysplastic syndrome (Gnyszka, Jastrzębski, and Flis 2013). However, potential side effects of these drugs necessitate the design of less toxic and more specific inhibitors of DNMT1. Ever since the breakthrough of epigenetics in cancer treatment, several in silico studies such as pharmacophore and QSAR modeling (Yoo, Kim, and Robertson 2012; Phanus-Umporn et al., 2020) have been reported for drug discovery against DNMT1 through small modulators. However, to the best of authors’ knowledge, no model is developed that could learn the diverse features of DNMT1 inhibitors. Therefore, in this study, we adopted machine-learning approaches to identify potential 2D descriptors of DNMT1 modulators by utilizing all the available chemical scaffolds of DNMT1 inhibitors (i.e., natural compounds, synthetic, nucleoside inhibitors, and FDA drugs).

ML classifiers DT and ANN were developed to classify the active and least-active compounds of the DNMT1 data set. All 10 captured features of the best-performing model (DT and ANN) along with feature descriptions are provided in Section 4, Supplementary Table S4. The trained DT identified vsa_acc, h_logP, b_count, radius, PEOE_VSA-5, b_double, SlogP_VSA5, SMR_VSA2, Q_VSA_PNEG, and Kier3 as key descriptors of demethylating compounds targeting DNMT1. However, the overall vDW surface area of hydrogen bond acceptors (vsa_acc) was reflected as the most prominent and distinguishing 2D descriptor by the DT model. This feature has also been previously identified by Hassanzadeh et al. (2017), which further emphasizes the significance of HBA for the activity of DNMT1 modulators. Furthermore, the developed models were able to classify the training set and predict the test set with optimal accuracies (table 1). The final selected trained DT and ANN classifiers showed an optimal classification accuracy of 0.97 and high sensitivity and specificity along with MCC values of 0.93 and 0.92, respectively.

The subsequent screening of the diverse test set from the DT and the ANN model resulted in the classification accuracies of 0.83 and 0.72, respectively. However, the DT classifier outperformed the ANN model showing higher classification accuracy of 0.83 and sensitivity of 0.92. Notably, despite the structural diversity of the DNMT1 data set, the ML classifiers developed in the current study retained an optimal classification accuracy of 97% for the training data and greater than 70% for the test data. The relatively low predictive accuracy of ML models on the test set can be attributed to the diversity of DNMT1 inhibitors as they differ highly in terms of chemical structure, molecular weight, and other pharmacological variables such as clogP and lipophilicity. Therefore, the learned features of the training set might not be sufficient to fully predict and explain the behavior of the test set. Also, the FDA-approved drug (azacytidine) was fairly classified among most active modulators of DNMT1 by our developed model, which strengthens the classification and predictive ability of the ML models.

RUNX3 has been proposed as a potential biomarker in TNBC for an early prognosis, which is known to be downregulated by DNMT1. However, the precise mechanism of epigenetic-mediated silencing of TSG RUNX3, which results in cancer invasion and metastasis, has not yet been explored. This study deals with the formal modeling of the DNMT1–RUNX3 signaling and the development of ML models on a diverse data set of DNMT1 modulators. First, we employed a qualitative modeling approach to provide an insight into epigenetic-inspired RUNX3 signaling. The results revealed that the onset of oncogene c-myc introduces pathogenesis in the system and its consistent activation along with persistent suppression of TSG RUNX3 hyperactivates DNMT1 leading to cancer metastasis. Conversely, the activation of RUNX3 leads the system to acquire normal homeostasis by transactivating other TSGs such as p21 and p53. Moreover, our findings advocate DNMT1 as a potential epigenetic drug target to revive the suppressed TSG RUNX3 in breast cancer therapeutics. Second, predictive ML models (DT and ANN) have been developed that identify some potential 2D descriptors essential to modulate the DNMT1 activity, and the best-performing models have effectively classified the active and least-active inhibitors. The trained (DT and ANN) models acquired 97% classification accuracy on the training data set, and the subsequent screening of the test set through DT and ANN models achieved 83% and 72% predictive accuracy, respectively, emphasizing the optimal efficiency of the developed model. In general, the application of formal methods to unveil the network and model the underline genetic events responsible for DNMT1-inspired TSG RUNX3 silencing along with ML approaches to predict the 2D attributes of hypomethylating compounds (targeting DNMT1) could present new computational avenues for the treatment of breast cancer requiring epigenetic therapy.