Identification of Cross-Pathway Connections via Protein-Protein Interactions Linked to Altered States of Metabolic Enzymes in Cervical Cancer

Metabolic reprogramming is one of the emerging hallmarks of cancer cells. Various factors, such as signaling proteins (S), miRNA, and transcription factors (TFs), may play important roles in altering the metabolic status in cancer cells by interacting with metabolic enzymes either directly or via protein-protein interactions (PPIs). Therefore, it is important to understand the coordination among these cellular pathways, which may provide better insight into the molecular mechanism behind metabolic adaptations in cancer cells. In this study, we have designed a cervical cancer-specific supra-interaction network where signaling pathway proteins, TFs, and microRNAs (miRs) are connected to metabolic enzymes via PPIs to investigate novel molecular targets and connections/links/paths regulating the metabolic enzymes. Using publicly available omics data and PPIs, we have developed a Hidden Markov Model (HMM)-based mathematical model yielding 94, 236, and 27 probable links/paths connecting signaling pathway proteins, TFs, and miRNAs to metabolic enzymes, respectively, out of which 83 paths connect to six common metabolic enzymes (RRM2, NDUFA11, ENO2, EZH2, AKR1C2, and TYMS). Signaling proteins (e.g., PPARD, BAD, GNB5, CHECK1, PAK2, PLK1, BRCA1, MAML3, and SPP1), TFs (e.g., KAT2B, ING1, MED1, ZEB1, AR, NCOA2, EGR1, TWIST1, E2F1, ID4, RBL1, ESR1, and HSF2), and miR (e.g., mir-147a, mir-593-5p, mir-138-5p, mir-16-5p, and mir-15b-5p) were found to regulate two key metabolic enzymes, EZH2 and AKR1C2, with altered metabolites (L-lysine and tetrahydrodeoxycorticosterone, THDOC) status in cervical cancer. We believe, the biology-based approach of our system will pave the way for future studies, which could be aimed toward identifying novel signaling, transcriptional, and post-transcriptional regulators of metabolic alterations in cervical cancer.


INTRODUCTION
Cervical cancer is the fourth most frequently occurring cancer and the fourth leading cause of death in women worldwide with an estimate of 5,70,000 cases and 3,11,000 deaths in 2018. Approximately 80-85% of the deaths from cervical cancer occur in lower and middle-income countries compared to highincome countries (1,2). Squamous cell carcinoma (SCC) and adenocarcinomas are the two main types of cervical cancer. Above 90% of patients with cervical cancer belong to SCC (3). The persistent infection with human papillomavirus (HPV), a particularly high-risk type of HPV (mainly HPV16 and HPV18 type), is considered the primary cause of cervical cancer (4)(5)(6). Only HPV16 and HPV18 types are responsible for almost 70% of cases of cervical cancer globally (7). While infection by highrisk HPV is necessary for developing cervical cancer, it alone may not be sufficient. Various studies suggest that the pathogenesis of cervical cancer depends on various other factors acting in concert with disease-associated HPV types (8)(9)(10). Therefore, it is important to understand the molecular mechanism behind the development of cervical cancer.
Metabolic reprogramming is considered one of the emerging hallmarks of cancer cells, and it is essential for cancer cell growth and proliferation to evolve into a more aggressive malignant state (11,12). Understanding the coordination among various biological pathways, such as gene-regulatory, signaling, and metabolic pathways is important and may provide clues into the molecular mechanism of metabolic adaptation in cancer and associated cells. To understand that, one needs to investigate the molecular mechanism by which the impact of signaling, transcriptional, and post-transcriptional aberration is transgressed to metabolic reprogramming. Various studies demonstrated that the metabolic status in cancer cells is regulated by oncogenic changes in signaling pathways (13)(14)(15), transcription factors (TFs) (16)(17)(18), and miRNAs (19)(20)(21). However, these studies are focused either on a single molecule or pathways and may not capture the complex interconnectivity among various biological processes.
In the present study, we have designed a cervical cancerspecific supra-interaction network model incorporating transcriptome data onto a protein-protein interaction (PPIs) network to investigate novel molecular targets and connections regulating the status of metabolic enzymes. We have developed a biology framework of a comprehensive system where signaling (S) pathway proteins, miRNA, and TF-based gene-regulatory modules are connected to metabolic (M) pathway proteins through protein-protein interactors (PPIs; Figure 1). Initially, network topologically IINs, such as a hub, central node (CN), local network perturbing nodes (LNPNs), and global network perturbing nodes (GNPN), were identified in different [S-M, TFmetabolic (TF-M), and miRNA-metabolic] modules of cervical cancer-specific networks using graph theory approach previously reported by our laboratory (38). These IINs may serve as potential diagnostic and/or prognostic biomarkers in cervical cancer. Furthermore, signaling pathway proteins, TFs, and microRNA (miR) to metabolic enzymes interconnecting paths/links [S-PPI-M, TF-target genes (TG)-PPI-M, and miR-TG-PPI-M] were identified in cervical cancer. Publicly available transcriptomic data derived from cervical cancer patients were incorporated into the HMM-based mathematical modeling set-up to weigh and rank the interconnecting link/paths in each module. Additional confidence values based on biological and network topological properties (hub, CN, GNPN, and LNPN) were assigned to each gene/protein/miRNA in the paths/links identified after model implementation to extract out high confident connections/links specific to cervical cancer. In silico validation of these selected genes/proteins/miRNAs and paths has been performed through perturbation analysis, demonstrating the importance of certain genes/proteins/miRNAs forming critical inter-pathway connections. PPI links connecting to key metabolic enzymes, such as RRM2, NDUFA11, ENO2, EZH2, AKR1C2, and TYMS, are identified from signaling proteins (e.g., PPARD, BAD, GNB5, CHECK1, PAK2, PLK1, BRCA1, MAML3, and SPP1), TFs (e.g., KAT2B, ING1, MED1, ZEB1, AR, NCOA2, EGR1, TWIST1, E2F1, ID4, RBL1, ESR1, and HSF2), and miR (e.g., mir-147a, mir-593-5p, mir-138-5p, mir-16-5p, and mir-15b-5p) in cervical cancer scenario. Out of the six metabolic enzymes that are commonly linked by 83 paths/links, EZH2 and AKR1C2 were mapped with deregulated metabolite status. Further, comparative analysis of the identified genes/proteins/miRNAs and the associated molecular pairs and paths in different modules were performed using transcriptomics data obtained from cervical, breast, and ovarian cancer patients. This study led to novel inter-bio-molecular links between signaling, generegulatory components, and metabolic enzymes paving the probable way(s) to identify drug targets to inhibit cervical cancer progression in a more specific manner.

MATERIALS AND METHODS
Messenger Ribonucleic Acid (mRNA) and Micro Ribonucleic Acid (miRNA) Expression Datasets mRNA and miRNA expression datasets of cervical, breast cancer, and ovarian cancer patients were extracted from the Gene Expression Omnibus (GEO) database (39) to study the gene and miRNA expression profiles. For the mRNA expression datasets, similar microarray (Affymetrix microarrays) platforms were used for each cancer type to minimize the undesirable variations that occurred due to different microarray platforms. Further, the cancer samples in each dataset were considered in a 1:1 ratio with normal samples to avoid sample heterogeneity (as shown in Table 1).

Differential Expression Analysis
The differential expression analysis of each dataset was performed separately using the GEO2R web tool (40)(41)(42) available at the GEO database. Genes having log 2 FC ≥ +1.5 and log 2 FC ≤ −1.5 were considered as upregulated and downregulated genes, respectively. The genes with log 2 FC values between −1.5 and +1.5 were considered only as neutrally expressed genes (EG). Benjamini and Hochberg's method (43) was used to control the false discovery rate. Genes with an adjusted p ≤ 0.05 were considered significant. For ovarian cancer datasets, the log 2 FC and p-value threshold of ± 2 and ≤0.01, respectively, were considered for the selection of upregulated, downregulated, and EG. For differential miRNA expression adjusted p ≤ 0.05 and log 2 FC ≥ +1.0 and log 2 FC ≤ −1.0 were considered as thresholds for the identification of upregulated and downregulated miRNAs, respectively.

Construction of Human Protein-Protein Interaction Network
The HPPIN was constructed by extracting experimentally verified (confidence score ≥0.7) HPPIs available in the STRING v11.0 (44) database. The resulting network (proteins as nodes and edges demark interaction) consisted of 5,048 proteins and 18,044 interactions, respectively.

Construction of Cancer-Specific Protein-Protein Interaction Network
Differentially expressed genes (dEXP) and EG from each cervical cancer dataset (GSE9750, GSE63514, and GSE52904; Table 1) were mapped onto the HPPIN to construct a cervical CC-PPIN. The interactions were considered up to the second level (i.e., interactors of interactors). In the first level, interactions mediated by only deregulated genes were considered where their interactors could be either deregulated or neutrally expressed. The resulting network consisted of 2,240 proteins interconnected via 5,452 edges. Similarly, breast and ovarian CC-PPINs were constructed using the corresponding transcriptomic datasets ( Table 1).

Construction of Cancer-Specific Transcriptional Regulatory Network
The transcriptional regulatory network was constructed by collating deregulated TF-TG interaction and associated protein-protein interactions. Experimentally verified strong evidenced TF-TG interactions were retrieved from Human Transcriptional Regulation Interactions database (HTRIdb) (45) and Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRUSST) (46) databases and were merged generating 22,480 interactions among 697 TFs and 12,407 TG. The combined deregulated genes (differentially

Construction of Cancer-Specific Post-transcriptional Regulatory Network
The post-transcriptional regulatory network was constructed by collating miRNA-TG interactions and protein-protein interactions. The experimentally verified miRNA-TG interactions were retrieved from mirTarbase (47) and Tarbase (48) databases and merged, which consisted of 8,407 miRNA-TG interactions forming among 743 miRNA and 2,891 TG.
Subsequently, the same procedure described for the construction of the TF-TG-PPIN network was followed for the construction of the miRNA-TG-PPIN network. The miRNA-TG-PPIN using cervical cancer-specific miRNA transcriptomics data consisted of 1,017 edges connecting 1,718 nodes. Likewise, breast and ovarian TF-TG-PPINs consisting of 2,668, 2,657 nodes and 6,197, 5,702 interactions, respectively, were also constructed.

Characterization of Cancer-Specific Networks
The cancer-specific networks, PPIN, TF-TG-PPIN, and miRNA-TG-PPIN, described above were compared with the respective random networks of the same number of interactions. Ten random networks were generated by the NetworkX program (49) against each cervical cancer-specific network described above, and the degree distribution of each network was compared with the respective random networks. The degree distribution was calculated using the following formula: Where the degree distribution of network P(k) signifies the fraction of the node with degree k. For the network with a node size of N, n k nodes will have the degree k.

Identification of IINs
Topologically IINs (genes/proteins/miRNA) of the constructed cancer-specific networks described above were identified by utilizing procedures based on graph theory methods described earlier in Bhattacharyya and Chakrabarti (38). Identification of important interacting genes/proteins/miRNAs in the network is based on some independent network properties, such as hub (highly connected nodes in the network), CNs of the network, GNPN, and LNPN. The nodes (genes/proteins/miRNAs) that were identified as topologically important in at least two categories (Hub, CN, LNPN, and GNPN) were considered as IINs.

Over-representation Analysis
Kyoto encyclopedia of genes and genomes (KEGG) pathwaybased ORA was performed with deregulated genes extracted from the mRNA expression datasets used and IINs (except miRNAs) identified in each of the regulatory networks described above using "protein-coding gene set" as the reference gene set in WebGestalt (50) web tool. The top 20 pathway categories were ranked based on significant false detection rate (FDR) calculated using Benjamini and Hochberg procedure (43) and enrichment ratio. Additionally, Gene Ontology (GO)-based molecular functions and online mendelian inheritance in man (OMIM)-based disease pathway over-representation analyses were also performed for the deregulated genes/proteins in cervical cancer.

Construction of S-M Enzyme Cross-Connecting Paths and Network
A signaling-metabolic inter-connection network was constructed using 23 signaling pathway (cancer-specific) genes/proteins and all the metabolic pathways (85 pathways) genes/proteins. Signaling and metabolic gene/protein datasets were created by extracting all the genes/proteins from the KEGG (51, 52) database. All possible unique connections (maximum three proteins involved in between) to a metabolic pathway protein (M) were established through PPIs (up to the second level), considering a signaling pathway protein (S) as a starting point in the HPPIN. Four different types of linking paths were established where signaling proteins were connected to metabolic pathway proteins either directly (S-M) or via one (S-P-M), two (S-P-P-M), or three (S-P-P-P-M) PPIs, respectively. NetworkX program (49) was used to construct all possible signaling to metabolic interconnecting paths. These paths/connections were converted into a network to construct a signaling-metabolic interaction network (SMIN).

Construction of TF to Metabolic Enzyme Cross-Connecting Paths and Network
Connections between TFs and metabolic pathway genes/proteins were established through TF-TG interactions and PPIs of TF-TG (up to the second level) considering TFs as a source. In this study, five different types of the path were established where TFs were connected to metabolic pathway proteins either directly (TF-TG/M), or TFs were connected to their TG, and their TG were connected to metabolic proteins directly (TF-TG-M), or through one (TF-TG-P-M), two (TF-TG-P-P-M), and three (TF-TG-P-P-P-M) PPIs of TG, respectively. The resulting paths/links were converted into a network to construct a TF-metabolic interaction network (TFMIN).

Construction of miRNA to Metabolic Enzyme Cross-Connecting Paths and Network
MicroRNAs to metabolic pathway proteins interconnecting all possible paths were established using the NetworkX program (49). The miRs to metabolic pathway proteins (M) interconnecting paths were established using miRNA-TG interactions and PPIs (up to the second level) of miRNA TG considering miRNA as the source. The resulting paths were of five types viz; miRNA-TG/M, miR-TG-M, miR-TG/P-P-M, miR-TG/P-P-P-M, and miR-TG/P-P-P-P-M paths. The resulting paths were converted into a network to form a miR-metabolic interconnecting network (miRMIN).

Contextualization of Regulatory Molecules (Signaling Pathway Proteins, TF, and miRNA) to Metabolic Enzyme Cross-Connecting Paths and Network
The deregulated (upregulated and downregulated) and neutrally EG and miRNAs identified from the cervical, breast, and ovarian cancer patients specific transcriptomic datasets were mapped onto all possible paths/connections mentioned above to filter cancer-specific regulatory molecules (signaling pathway proteins, TF, and miRNA) to metabolic enzymes cross-connecting paths. For the signaling to metabolic interconnecting paths/connections/links, the paths having deregulated (upregulated and downregulated) genes/proteins at the terminals and deregulated or EG/proteins in middle were filtered out and converted into a network to form cervical, breast, and ovarian cancer-specific signaling-metabolic interconnecting subnetwork (CC-SMIN, BC-SMIN, and OC-SMIN, respectively).
Similarly, to construct the cancer-specific TF and miRNA to metabolic interconnecting sub-networks (CC/BC/OC-TFMIN and CC/BC/OC-miRMIN, respectively), the paths having deregulated genes/miRNA at the terminal and their target and deregulated or EG/proteins in the middle were considered. The respective resulting paths were converted into a network to form cancer-specific (CC/BC/OC-TFMIN and CC/BC/OC-miRMIN) sub-networks.

Calculation of Edge Weight
The edge weights in each sub-networks (CC/BC/OC-SMIN, CC/BC/OC-TFMIN, and CC/BC/OC-miRMIN) were defined in terms of local entropy using an in-house program (37). The probability of interactions of a gene/protein with its interactors in the sub-network was determined by using the principle of mass action to define the local entropy of a gene/protein. The calculation of the interaction probabilities is based on the assumption that two proteins known to interact will have a higher probability of interaction when they are highly expressed. The normalized expression values of sub-network genes in the samples of cancer patients used in this study were utilized to calculate the interaction probabilities.

Calculation of Node Weight and Effect-On-Node
To incorporate the importance and impact of the interactors of a particular node in the sub-networks, the node-weight (W i ) of every node i was defined based on its biological properties [signaling cross-talk (SC) protein and rate-limiting enzyme (RLE)], differentially expressed gene (dEXP), and network topological properties [hub, CP, local network perturbing protein (LNPP), and global network perturbing proteins (GNPP)] using the following formula.
Effect of interactors on a node in the sub-network was defined as effect-on-node (effs) depending on the node weight of its interactors up to the second level.
Where k is the degree of node s, n i is the degree of node i, n j is the degree of node j, and w i and w j are the weights of nodes i and j.

Identification of Significant Regulatory Molecules (S/TF/miR) to a Metabolic Enzyme (M) Interconnecting Pairs and Path
To understand the information flow starting from a regulatory molecule (S/TF/miRNA) to metabolic enzyme (M), cancerspecific cross-connecting paths mentioned above were scored by implementing an HMM-based mathematical model established in our laboratory earlier (37). In this study, two separate models were used to identify the significant S/TF/miR-M pairs (the source; signaling pathway protein/TFs/miR; and destination: metabolic pathway protein) and S/TF/miR-M interconnecting paths. Model 1 was applied to identify the S/TF/miR-M pairs. Model 2 was applied to identify the S/TF/miR-M interconnecting paths between the S/TF/miR-M pairs selected after Model 1. Edge weight and node weight of genes/proteins/miRNAs involved in the S/TF/miR-M path were used to calculate the path scores. The path score of each S/TF/miR to M linking path calculated by Model 1 and Model 2 was converted into a statistical z-score to identify paths deviating from the mean. A z-score (Z) ≥ 1 filter was applied to select the significant S/TF/miR-M pairs. Paths having path score ≥ 80% of the highest path score for every S/TF/miR-M pair were considered as significant S/TF/miR-M interconnecting paths from Model 2. All the networks were visualized and represented by using Cytoscape (53). The signaling, TF, and miR to metabolic pathway connections were represented as the Circos plot (54).

In-silico Perturbation Analysis
In-silico perturbation analysis was performed for each signaling pathway protein, TF, and miRNA-based gene regulatory module to identify the paths that change significantly upon removal of a node (protein/TF/miRNA). To identify the key nodes in the final paths/networks of every module, each of the nodes present in the paths/network having z-score ≥ 1 was removed individually from the HPPIN, and the path score was recalculated for the resulting paths/network by using the HMM Models 1 and 2. The perturbation score was calculated by using the average path score before and after perturbation as below.
Perturbation score = Path score ′ − Path score Where, Path score and Path score ′ are average path scores before and after perturbation, respectively.
The difference of average path score (before vs. after perturbation) for each perturbed node was converted into a zscore (Z), and the nodes for which z-scores deviated from the mean as −1 ≥ Z ≥ 1 were selected as effective or key nodes in significant paths/network.

Metabolomics Data Collection and Integration Into Cancer-Specific Paths
The deregulated metabolites in cervical, breast, and ovarian cancer patient were extracted from literature (55)(56)(57). The cervical cancer metabolomic dataset consisted of 55 downregulated and seven upregulated metabolites. Twenty-one downregulated and 41 upregulated metabolites were found in breast cancer whereas the ovarian cancer metabolomic dataset consisted of 46 downregulated and 116 upregulated metabolites. The metabolic genes corresponding to these metabolites were obtained from the Human metabolome database (HMDB) (58). The deregulated metabolites were mapped to the paths obtained after model implementation.

Survival Analysis
Kaplan-Meier (KM) plotter software (59,60) was used to perform the overall survival (OS) analysis of the constituent genes/miRNAs of the identified cross-pathway paths/links. We used a KM plotter using survival and expression data of 307 cervical cancer patients obtained from the TCGA dataset (project ID: TCGA-CESC; phs000178). To estimate the survival prognostic value of a specific gene/miRNA, the patient samples were divided into high-and low-expression cohorts according to the median expression of the given gene/miRNA, and KM plots were created. Additionally, the hazard ratio (HR) and the log-rank p-value were calculated. The survival estimate of a gene/miRNA with a p-value < 0.05 was considered to be statistically significant.

Drug/Chemotherapy Response Analysis
The receiver operator characteristic (ROC) plotter (61) was used to predict the utility of the genes as predictive biomarkers with respect to drug/chemotherapy response. ROC plotter is capable to link gene expression and response to therapy using transcriptome-level data of 3,104 breast cancer patients.

mRNA and miRNA Expressions in Cervical Cancer Patients
The individual differential expression analysis of three different cervical cancer mRNA expression datasets (GSE9752, GSE63514, and GSE52904) leads to the identification of several upregulated, downregulated, and neutrally EG in cervical cancer (Figure 2, Table 1). However, a comparison of the gene expression patterns across datasets showed that only 54 upregulated, 59 downregulated, and 421 neutrally EG were found to be overlapped (Figures 2A-C). The differential expression analysis of miRNA resulted in four upregulated, 12 downregulated, and 110 neutrally expressed miRNAs in GSE30656 and 15 upregulated, 20 downregulated, and 103 neutrally expressed miRNAs in the GSE81137 dataset (Figures 2D-F).
Over-representation analysis-based enrichment for cellular pathways, molecular functions, and biological processes was performed using a merged list of deregulated (upregulated and downregulated) genes of all the three mRNA datasets. Supplementary Figures 1A,B show the top 20 most enriched pathway filters based on FDR for proteins encoded by upregulated and downregulated genes, respectively. Most highly enriched pathways were found to be DNA replication (p = 4.42E −13 ) and arachidonic acid metabolism (p = 2.76E −08 ) for the proteins encoded by upregulated and downregulated genes, respectively. GO-based molecular function and biological process ORA was also performed using deregulated genes in cervical cancer. The most significantly enriched molecular functions were found to be collagen binding (p = 4.36E −07 ) and oxidoreductase activity (p = 4.35E −04 ) for proteins encoded by upregulated and downregulated genes, respectively. However, microtubule cytoskeleton organization involved in mitosis (p = 2.13E −14 ) and peptide cross-linking (p = 8.81E −11 ) were highly enriched biological processes, for upregulated gene-encoded proteins and downregulated gene-encoded proteins, respectively (Supplementary Figure 1).

Construction and Characterization of Cervical Cancer-Specific Networks
The cervical CC-PPIN, TF-TG-PPIN, and miR-TG-PPIN were constructed by mapping differentially expressed genes, miRNA, and neutrally EG from each cervical cancer dataset (see Methods; Supplementary Figures 2A, 3A, 4A, respectively). The resulting CC-PPIN, TF-TG-PPIN, and miR-TG-PPIN were validated by comparing them with corresponding 10 random networks of the same size. The degree distributions of CC-PPIN, TF-TG-PPIN, and miR-TG-PPIN networks followed the Power law and were considered to have scale-free organization. However, the degree distribution of 10 corresponding random networks showed binomial distribution (Supplementary Figures 2B, 3B, 4B).

IINs in Cervical Cancer-Specific Networks
Network topologically important nodes, such as hubs (highly connected nodes in the network), CNs, LNPNs, and GNPNs in a scale-free network could play important roles in maintaining the network integrity and function. Various topologically important interacting genes/proteins/miRNAs (hubs, CN, GNPN, and LNPN) in cervical cancer-specific networks (CC-PPIN, TF-TG-PPIN, and miR-TG-PPIN) were identified by implementing a graph theory-based method described earlier by Bhattacharyya and Chakrabarti (38). A total of 165, 167, 96, and 67 nodes/proteins in CC-PPIN (Supplementary Figure 2C, Supplementary Table 1), 62, 36, 60, and 80 nodes/genes/proteins in TF-TG-PPIN (Supplementary Figure 3C, Supplementary Table 2), and 45, 45, 30, and 24 nodes/genes/proteins/miRNAs in miR-TG-PPIN (Supplementary Figure 4C, Supplementary Table 3) were identified as hubs, CNs, GNPNs, and LNPNs, respectively. The nodes/genes/proteins/miRNAs common in at least any two of the four categories (hubs, CN, GNPN, and LNPN) were considered as IINs in cervical cancer networks. When the IINs from each network were compared, 30 IINs were found to be common in CC-PPIN and TF-TG-PPIN, 38 IINs were common in CC-PPIN, and miR-TG-PPIN, and 17 IINs were shared by TF-TG-PPIN and miR-TG-PPIN. However, 17 IINs were found to be common in all three cervical cancer networks (Supplementary Figure 5A,  Supplementary Table 4).
Over-representation analysis-based pathway enrichment was performed using IINs in the above described cervical cancerspecific regulatory networks. Top three enriched pathways were found to be DNA replication (p = 1.30E −10 ), basal TFs (p  Figure 5B). The common pathways were cell cycle, DNA replication, nucleotide excision repair, mismatch repair, prostate cancer, herpes simplex infection, oocyte meiosis, and viral carcinogenesis pathways. To identify the potential disease-specific paths/pairs, each node and each edge of the CC-PPIN network were weighted based on their biological properties, differential expression status, and network topological properties (hubs, CN, GNPN, and LNPN) of CC-PPIN. Local signaling entropy (S i ) was integrated to understand the system-level network property. The significance of each node (gene/protein) in the cancer-specific network was estimated in the form of effect-on-node (effs) based on SC, RLE, dEXP (differentially expressed genes), hub, CN, GNPN, and LNPN, respectively. To identify the probable and significant paths of information flow from signaling pathway to metabolic pathway in cervical cancer cells, local signaling entropy (S i ) and effect-on-node (effs) properties were incorporated as node weights. The edge weight of every two interacting nodes of CC-PPIN was defined as the probability of interaction using their normalized expression value in cervical cancer patient samples. Node weight and edge weight were integrated into HMM-based mathematical models (Models 1 and 2) to identify S-M linking pairs and paths. Model 1 was applied to identify the S-M pairs (the source signaling pathway protein and destination metabolic pathway protein). Model 2 was applied to identify the S-M interconnecting paths between the S-M pairs selected after Model 1. The path score of each S-M linking path calculated by Models 1 and 2 was converted into a statistical z-score to identify paths deviating from the mean. A z-score ≥ 1 filter was applied to select the significant S-M pairs. Using these filtering criteria, we identified 81 S-M pairs and 94 S-PPI-M paths in cervical cancer. The selected paths were converted into a network to form a significant CC-PPIN network which consisted of 152 interactions forming among 135 genes/proteins (Table 2, Figure 3A).

S-M Enzyme Cross-Connecting Paths and Network in Cervical Cancer
Mapping pathway information to the terminal nodes (source signaling protein and destination metabolic enzyme) showed that the Ras signaling pathway had maximum connections (50) to all the six metabolic pathways followed by cell cycle (44), Map-kinase (44), epidermal growth factor receptor (EGFR) signaling pathway (42), and p53 signaling pathway (34), respectively. However, among metabolic pathways, nucleotide metabolism had maximum connections (87) with signaling pathways, followed by amino acid (63), energy (58), xenobiotics (44), carbohydrate (42), and lipid metabolism (11), respectively. 1:1 interconnections between signaling and metabolic pathways showed that the cell cycle had maximum connections with nucleotide metabolism (18 connections; Figure 3B, Supplementary Table 5).  Mapping the deregulated metabolites in cervical cancer to significant S-M paths yielded 12 S-PPI-M interconnections/paths where four metabolites [L-lysine, oxoglutaric acid, tetrahydrodeoxycorticosterone (THDOC), and pyruvic acid] were regulated by eight signaling pathway proteins (BAD, CHEK1, GNB5, MAML3, MAP3K1, PAK2, PPARD, and SPP1). The metabolic enzymes connecting these four metabolites were enhancers of zeste homolog 2 (EZH2), procollagen lysine hydroxylase and glycosyltransferase LH3 (PLOD3), aldo-keto reductase family 1 member C2 (AKR1C2), and 3-mercaptopyruvate sulfurtransferase (MPST). L-lysine is the substrate of both EZH2 and PLOD3. Whereas, THDOC and pyruvic acid are the products of metabolic enzymes AKR1C2 and MPST, respectively. Hence, these paths/connections showed the correlated status of the metabolic enzymes and the corresponding metabolites ( Figure 3C). Each node and each edge in the CC-TFMIN were weighted to identify the potential significant paths in the cervical cancer-specific network. Node weights and edge weights were incorporated into HMM-based mathematical models (Models 1 and 2) to identify TF-M pairs and TF-PPI-M paths forming between them (refer to Methods). Model 1 resulted in the identification of 172 significant (z ≥ 1) TF-M pairs and Model 2 resulted in 236 significant TF-PPI-M paths connecting the TF-M pairs obtained after Model 1. The significant TF-PPI-M paths in cervical cancer were converted into a network to form significant CC-TFMIN that consisted of 179 interactions among 141 genes/proteins ( Figure 4A, Table 3).

TF to Metabolic Enzyme Cross-Connecting Paths and Network in Cervical Cancer
The metabolic pathway information was mapped onto the terminal metabolic enzymes of the significant TF-PPI-M paths to identify metabolic pathways that are highly connected to specific TFs. Nucleotide metabolism yielded maximum connections followed by energy metabolism, amino acid metabolism, carbohydrate metabolism, lipid metabolism, and xenobiotics biodegradation and metabolism (Figure 4B,  Supplementary Table 6). Mapping of deregulated metabolites onto the terminal metabolic enzymes resulted in 28 paths where four metabolites were deregulated in cervical cancer. The deregulated metabolites were L-lysine, oxoglutaric acid, pyruvic acid, and THDOC. Llysine was found to be linked with AR, ESR1, ZEB1, NCOA2, HSF2, RBL1, ID4, E2F1, EGR1, MED1, TWIST1, KAT2B, ING1, and PGR TFs via 19 different paths. Oxoglutaric acid was linked with PGR, EGR1, and ESR1 via three paths. MED1 was found to be regulating (probably) pyruvic acid whereas THDOC was found to be linked with AR, TWIST1, and ING1 in four different paths ( Figure 4C).

miR to Metabolic Enzyme Cross-Connecting Paths and Network in Cervical Cancer
Similar to SMIN and TFMIN, all possible paths/links were established considering miRNA as a source node and metabolic enzymes as a destination by collating miRNA-TG interactions and HPPIN. The resulted paths/links consisted of 577 direct (miR-TG/M) paths, 1,145 via their TG (miR-TG/P-M), 26,330 via one PPI (miR-TG/P-P-M), 826,207 via two PPI (miR-TG/P-P-P-M), and 33,934,931 via three PPI (miR-TG/P-P-P-P-M) paths formed by 577, 1,128,9,904,44,271, and 111,730 miR-M pairs, respectively. The deregulated miRNA, genes, and neutrally EG were mapped onto all possible miR-PPI-M paths to filter the cervical cancer-specific miRNA to metabolic enzymes interconnections. The filtered paths/links consisted of 13 miR-TG/M, 1 miR-TG/P-M, 7 miR-TG/P-P-M, 95 miR-TG/P-P-P-M, and 1,851 miR-TG/P-P-P-P-M paths formed between 13, 1, 7, 38 and 149 miR-M pairs, respectively. The resulted paths were converted into a network to form a cervical cancer-specific miRmetabolic enzyme interaction network (CC-miRMIN) which consisted of 309 nodes and 952 interactions among them ( Table 4).
The nodes and edges in the CC-miRMIN were weighted based on their biological properties, differential expression status, and network topological properties in cervical cancer-specific  Figure 5A).
Mapping metabolic pathway information to the terminal metabolic enzymes in significant miR-PPI-M paths showed that amino acid metabolism was highly regulated by miRNAs, followed by nucleotide metabolism, xenobiotics biodegradation and metabolism, energy metabolism, carbohydrate metabolism, and lipid metabolism (Figure 5B).

In-silico Perturbation of Nodes in the Final Weighted Paths/Network
In-silico perturbation analysis was performed to identify the paths that change significantly upon removal of a node (protein/TF/miRNA). To identify the key nodes in the final paths/networks of every module, each of the nodes present in the paths/network having Z ≥ 1 was removed individually from the HPPIN, and the path score was recalculated for the resulting paths/network by using HMM Models 1 and 2. Accordingly, new significant pairs and paths were identified based on Models 1 and 2, respectively. The difference of average path score (before vs. after perturbation) for each perturbed node was converted into a TF-TG/M, TFs connected to metabolic pathway proteins; TF-TG-M, TFs connected to TG and TG connected to metabolic proteins directly; TF-TG-P-M, TFs connected to TG and TG connected to metabolic proteins directly through one PPI of TG; TF-TG-P-P-M, TFs connected to TG and TG connected to metabolic proteins directly two PPIs of TG; TF-TG-P-P-P-M, TFs connected to TG and TG connected to metabolic proteins directly three PPIs of TG. z-score and the nodes for which z-scores deviated from the mean as −1≥ Z ≥ 1 were selected as effective or key nodes in significant paths/networks. Sixteen nodes/proteins (CDC5L, PAK2,  CHECK1, NDUFA9, MCM, POLA1, PIK3CA, PIK3R1,  PDGFRA, LUC7L3, SERPINE1, VTN, ZRANB2, TYMS, CD14, and NDUFA11) in the signaling module (Figure 2A), 16 nodes/proteins (CDKN2A, POLA1, CDC45, CCND1, MCM5, CDC7, RRM2, MCM3, DHFR, AKR1C3, AKR1C2, AKR1C1, RRM1, TYMS, AR, and E2F1) in TF-based gene regulatory module (Figure 4A), and nine nodes/miRNA/proteins (CDK2, MCM3, mir-147a, AKR1C2, CDC5L, PLK1, mir-593-5p, TYMS, and mir-196a-5p) in miRNA-based gene regulatory module ( Figure 5A) were identified as an effector or key nodes in the significant paths/networks.

Survival Analysis of the Genes/miRNAs of Identified Paths/Links in Cervical Cancer
Potential prognostic values of the genes and miRNAs of the signaling, transcriptional and post-transcriptional crossconnecting paths/links to metabolic enzymes in cervical cancer patients were explored by evaluating the correlation and OS. A total of 53 genes and 10 miRNAs were found to be significantly associated with the OS in the log-rank test with a p < 0.05 (Supplementary Table 12). Mapping these genes and miRNAs onto the corresponding cross-connecting paths/links yielded 16, 34, 20, and 9 paths/links to have 1-25, 26-50, 51-75, and 76-100% of their component as a prognostic marker in cervical cancer patients ( Figure 7A). Almost all the final selected paths/links (79 out of 83) possess at least one node (gene/miRNA), whose expression is significantly associated with cervical cancer patients' survival. In total, 38% (30 out of 79) of the selected paths/links have more than 50% nodes to be significant (p < 0.05) prognosis marker ( Figure 7A). Further, we checked the status of these prognostic markers in different types of paths/links. Most of the two-component (2C) paths were found to have 100% of their component as a prognostic marker. Three-component (3C) paths were found to have 51-75% of their component nodes as a prognostic marker. Similarly, significantly higher numbers of longer or higher component paths (e.g., 4, 5, and 6C, respectively) also possess more than 25% of their nodes as a prognostic marker ( Figure 7B).

S, TF, and miR Cross-Talks in Breast and Ovarian Cancers
To investigate whether the cross-pathway links are specific to cervical cancer, we identified such paths from two other female-specific cancers, such as breast and ovarian cancers. As mentioned in the Methods, paths originating from S, TF, and miR connecting metabolic enzymes (M) were identified in breast and ovarian cancers using the cancer-specific transcriptomics data mapping followed by implementation of HMM-based mathematical models. , and posttranscriptional (miR-M) regulatory links identified from cervical, breast, and ovarian cancer networks were compared to estimate the common and specific regulators and regulatory links (Figure 8). Interestingly, a very little overlap of regulatory paths and pairs was observed among the three types of cancers (Figures 8A-F). In total, 32% of terminal signaling proteins and 43% of terminal metabolic enzymes forming CC-specific S-M enzymes paths were found to be common with that extracted from breast and ovarian cancer networks (Figures 8G,J). Similarly, 46% of terminal TFs and 30% of terminal metabolic enzymes forming CC-specific TF-metabolic enzyme paths were found to be common with those extracted from breast and ovarian cancer networks (Figures 8H,K). Three metabolic enzymes (EZH2, ENO2, and RRM2) were found to be commonly regulated by S-M, TF-M, and miR in all three cancers (Figure 8).
The metabolic enzyme EZH2 was connected to deregulated metabolites L-lysine in cervical, breast, and ovarian cancer ( Figure 6B, Supplementary Figure 6). Twenty-five out of the 27 paths/links connected to metabolic enzyme EZH2 in the breast cancer network (Supplementary Figure 6A) possesses at least one gene, whose expression is significantly associated with drug/chemotherapy response. Seventeen out of the

DISCUSSION
Understanding the molecular mechanisms for cancer progression and subsequent development of potential therapeutics to inhibit this complex disease are difficult from the independent knowledge of ongoing signaling, gene regulatory, and metabolic alterations. Therefore, understanding the intricate coordination of signaling and gene regulatoryinduced proliferation of tumor cells/growth and metabolic processes is very much required. An integrated view of the probable interconnections between oncogenic signalinggene regulatory pathways and the metabolic shift could be one of the better ways to find out possible potential therapeutic targets. Our approach toward the establishment of a cross-pathway metabolic interconnection network is an attempt in that direction.
It is well-established that genetic modifications, altered transcriptional, and post-transcriptional regulations are responsible for mediating the changes in biological processes, which ultimately shape complex pathophysiological situations like cancer. The interconnectivity and regulations are perhaps maintained through the systemic-coordinated interaction of proteins as a complex system, acting as a perfect molecular machine (62)(63)(64). Therefore, the identification of such a precise protein-interaction network responsible for the disease progression is of utmost importance for understanding the disease and potential therapeutic development.
In this study, we have developed a biology framework of a cervical cancer-specific system where signaling (S) pathway proteins, miRNA, and TF-based gene-regulatory modules are connected to metabolic (M) pathway proteins through PPIs. Publicly available transcriptomic data derived from cervical cancer patients were incorporated into a mathematical modeling FIGURE 6 | Signaling pathway proteins, transcription factor, and miR cross-connecting paths/links to metabolic enzymes in cervical cancer. Six metabolic enzymes (A-F) are commonly linked to signaling proteins, transcription factors, and miR. Two metabolic enzymes, EZH2 and AKR1C2 (B,D), are connected to deregulated metabolites in cervical cancer. Terminal signaling pathway proteins, transcription factors, and miRNAs, metabolic enzymes are colored in purple, green, red, and blue. Protein-protein interactors are colored in orange. Gene regulatory edges are represented as black arrows and protein-protein interactions are represented by orange edges. Nodes with an asterisk (*) are key or effector nodes in the significant paths/network.  set-up to weigh and rank the interconnecting link/paths in addition to biological and network topological properties to extract out high confidence inter-pathway connections that are perhaps responsible for facilitating metabolic adaptation in cervical cancer.
In our previous study, we implemented the underlined mathematical model-based approach for the development, test, and validation of S-M interconnecting links using glioblastoma multiform (GBM) cell line-derived transcriptomics and proteomics data (37). Further, in-vitro perturbation of genes/proteins involved in forming a high-score interconnection between S-M pathway proteins showed a significant change in the expression of proteins involved in the metabolic pathway. This validated our model for discovering hitherto unknown connections/involvement between signaling and metabolic genes/proteins. As a natural follow-up study, here we have significantly upgraded the previous model with two entirely new types of connectivity paths linking TF and miRNA-based regulatory mechanisms to altered states of metabolic enzymes. We have used large-scale patient-derived cervical cancer data and implemented additional network topology-based weights to signify the identified cross-pathway links. Further, we have utilized differential metabolite data to extract out paths that correlate with the altered status of the metabolic enzymes that were proposed to be regulated via signaling and regulatory factors. Comparison of the significant paths originated with S, TFs, and miRNAs yielded 88 commonly linked paths connecting six common metabolic enzymes (e.g., RRM2, AKR1C2, ENO2, TYMS, EZH2, and NDUFA11).
The ribonucleotide reductase subunit M2 (RRM2) was found to be significantly upregulated in cervical cancer tissue and is linked to promoting the progression of cervical cancer (65). RRM2 is likely to become a novel potential diagnostic and prognostic biomarker of cervical cancer.
Aldo-keto reductase subfamily 1C2, which plays a major role in regulating the activity of androgens, estrogens, progesterone, and prostaglandins metabolisms, is also implicated with cervical, endometrial, and bladder cancers (66,67). Overexpression of AKR1C2 is found to be a mildly favorable prognostic marker (Supplementary Table 12) but lower expression of NADH:ubiquinone oxidoreductase subunit A11 (NDUFA11) is prognostically unfavorable in cervical cancer (Supplementary Table 12). Enolase 2 (ENO2) and thymidylate synthetase (TYMS) are found to be upregulated in cervical cancer transcription datasets (68,69). Overexpression of enhancer zeste homolog 2 (EZH2) has been linked with proliferation, progression, and prognosis of cervical cancer (70). However, our survival analysis using data of cervical cancer patients from TCGA suggested much lower survival with lower expression of EZH2 (Supplementary Table 12).
Several miRNAs have been identified whose roles have been implicated in cervical cancer progression. Most of the miRNAs except one (miR-593-5p) forming the metabolic pathway PPI links were previously reported to be dysregulated in cervical cancers (71)(72)(73)(74)(75). However, the diverse mechanisms by which these miRNAs could regulate cervical cancer progression were not well-known especially their roles in regulating the metabolic adaptation in cervical cancer cells. Our study provides novel avenues to study the impact of these important miRNAs in the regulation of metabolic reprogramming in cervical cancer. miR-593 plays important role in the regulation of lung, breast, and gastric cancer proliferation (76)(77)(78)(79). Higher expression of miR-593 is found to be unfavorable for the survival of cervical cancer patients (Supplementary Table 12). Hence, its role in cervical cancer especially in its metabolic adaptation is worth investigating further.
Transcription factors are key regulators of cancer proliferation and metastasis. Roles of several transcription factors, such as SOX2, E2F4, E2F1, POU5F1, SMAD3, SMAD2, VDR, ERG, TP53, EWS, c-fos, fra-1, OCT4, KLF4, C-MYC, and NANOG, were established in cervical cancer (80,81). Our study highlights the probable roles of important transcription factors in regulating the metabolic status of cervical cancer cells via modulating the metabolic enzymes. Thirty-one TFs were found to be connected to 30 metabolic enzymes via the TG and their respective PPIs. Fifteen TFs were found to be linked with six metabolic enzymes for which altered metabolite status could be associated (Figure 4).
Ras, cell cycle, MAPK, EGFR, and p53 are among the top five most connected signaling pathways to the metabolic enzymes via PPI interconnectivity (Figure 3). Among the 27 terminal signaling proteins that form significant connections with metabolic enzymes, 9 (CHEK1, MAP3K1, CDKN2C, PAK2, EGFR, FGFR2, PDGFRA, and PTK2) are found to be kinases. TFs and miRNAs are generally regarded as "undruggable, " hence, these regulatory kinases could be ideal candidates for targets of small molecules inhibitors/drugs to check their roles in altering the functional activities of the connected enzymes.
All cervical cancer-based cross-pathway links are provided in Supplementary Tables 8-11. Similarly, an online platform is also created as a separately published work (82) where cervical cancer dataset-specific S-M, TF-metabolic, miRNA-metabolic, and combined paths are made available at http://www.hpppi.iicb. res.in/APODHIN/home.html.
Cervical cancer-based PPI links to metabolic enzymes originated from signaling (S), TF, and miR regulatory molecules were also compared to the same identified from breast and ovarian cancer networks. Comparison of the regulators and regulatory links yielded little overlap among the three cancers (Figure 8) indicating the existence of cancer-specific regulatory mechanisms for probable metabolic alterations. However, some signaling proteins were found to be common regulatory molecules among the three cancers whereas terminal enzymes, such as EZH2, ENO2, RRM2, and TYMS, were found to be commonly regulated in all the cancers by three different regulatory mechanisms (Figure 8).
We understand that our approach is computation heavy and our findings require further in vitro and in vivo experimental validations. Similarly, for effective stratification of the patients, multiple omics data (e.g., transcriptomics, proteomics, and metabolomics) need to be generated from an individual patient. Nevertheless, we believe that our systems biology-based approach of identifying multi-factor signature links connected to the regulation of status and functionalities of metabolic enzymes paves the way for future studies which could be aimed toward identifying novel regulators of metabolic alterations in cervical cancer.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: Gene expression omnibus (GEO).

AUTHOR CONTRIBUTIONS
SC conceptualized the project. KK and SB performed the analysis. KK and SC analyzed the data and drafted the manuscript. All authors contributed to the article and approved the submitted version.