Network Biomarkers of Bladder Cancer Based on a Genome-Wide Genetic and Epigenetic Network Derived from Next-Generation Sequencing Data

Epigenetic and microRNA (miRNA) regulation are associated with carcinogenesis and the development of cancer. By using the available omics data, including those from next-generation sequencing (NGS), genome-wide methylation profiling, candidate integrated genetic and epigenetic network (IGEN) analysis, and drug response genome-wide microarray analysis, we constructed an IGEN system based on three coupling regression models that characterize protein-protein interaction networks (PPINs), gene regulatory networks (GRNs), miRNA regulatory networks (MRNs), and epigenetic regulatory networks (ERNs). By applying system identification method and principal genome-wide network projection (PGNP) to IGEN analysis, we identified the core network biomarkers to investigate bladder carcinogenic mechanisms and design multiple drug combinations for treating bladder cancer with minimal side-effects. The progression of DNA repair and cell proliferation in stage 1 bladder cancer ultimately results not only in the derepression of miR-200a and miR-200b but also in the regulation of the TNF pathway to metastasis-related genes or proteins, cell proliferation, and DNA repair in stage 4 bladder cancer. We designed a multiple drug combination comprising gefitinib, estradiol, yohimbine, and fulvestrant for treating stage 1 bladder cancer with minimal side-effects, and another multiple drug combination comprising gefitinib, estradiol, chlorpromazine, and LY294002 for treating stage 4 bladder cancer with minimal side-effects.


Introduction
Bladder cancer is still one of the most common cancers worldwide. Single gene markers have been proposed for improving cancer treatment [1]. However, single gene markers cannot overcome treatment side-effects because the markers are not implicated in genome-wide networks, and the analysis of a genome-wide network is a complicated issue from a systems biology perspective. The rapid development of molecular biology techniques has produced a great deal of high-throughput experimental data, including genome-wide microarray data, genome-wide methylation profiles, nextgeneration sequencing (NGS) data, microRNA (miRNA) profiles, genetic sequences, protein abundance data, and drug response genome-wide microarray data. These kinds of omics data provide an opportunity to design multiple drug combinations for the treatment of bladder cancer by applying the network biomarkers identified by systems biology.
To date, genetic regulation systems, including proteinprotein interaction networks (PPINs) and gene regulatory networks (GRNs), have been applied to analyze the functional mechanisms behind human aging and cancer [2,3]. We now know that epigenetic alterations are much more rapid and adaptive with regard to influencing genome-wide gene expression than genetic changes [4]. Rapid and slow response mechanisms, that is, epigenetic alterations and genetic changes, respectively, coordinate an efficient and robust system. Epigenetic regulation, including DNA methylation and histone modification, results in potentially reversible alterations in gene expression that do not involve permanent changes to the DNA sequence. miRNAs that are influenced by aberrant epigenetic regulation also mediate the regulation of 2 Disease Markers gene expression [5]. It has been found that DNA methylation directly affects the binding affinities of miRNAs, RNA polymerase, and transcription factors (TFs) [6] and indirectly influences protein-protein interactions (PPIs) [7]. Methylation analysis of human genomic DNA in 12 tissues revealed that DNA methylation profiles are tissue-specific [8]. Therefore, omics data and systems biology methods [9][10][11] are required to unravel the mechanisms underlying carcinogenesis from the complex molecular biology and design anticancer drugs for the treatment of bladder cancer.
The Human Genome Project (HGP) has identified 30,000-40,000 genes in human DNA, including miRNAs. The genes, proteins, and their associations, miRNA regulation, and DNA methylation constitute the integrated genetic and epigenetic genome-wide network (IGEN), which coordinates cellular responses. PPIN in human lung cancer [12] and GRN in human aging [13] of the genes with significant expression differences between cancer cells (or aged people) and normal cells (or young people) have been identified for the extraction of the core network biomarkers according to the estimated association abilities between TFs (or upstream proteins) and target genes (or target proteins). Aging is associated with cancer [14]. The association abilities estimated by the network models assume that the binding affinities of TFs (or upstream proteins) to target genes (or proteins) are the same. According to a recent study in primary human somatic and germline cells [6], the impact of the binding affinities of miRNAs, RNA polymerase, and TFs on gene expression is mediated by DNA methylation. According to the available genome-wide methylation profiles and NGS data for bladder cancer in The Cancer Genome Atlas (TCGA), DNA methylation and miRNA regulation can be also characterized by the GRN model to identify the genomewide IGEN. In this study, we identified the IGENs in normal bladder cells and bladder cancer cells and then investigated the impact of epigenetic regulation and miRNA regulation on bladder carcinogenesis by comparing the IGEN in normal bladder cells with that in bladder cancer cells.
Although a genome-wide IGEN can be identified based on well-defined system identification techniques [2,12], the mean by which the core network biomarkers are extracted from the identified genome-wide network is still an important issue. The total association capabilities of a single node can affect the contribution it makes to its neighbors. However, the genome-wide IGEN including transcriptional gene regulations, miRNA regulations, and PPIs constitutes a genomewide network structure. The contribution made by one node to its neighbors is not sufficient to explain its impact on a genome-wide scale network of bladder cells. In this study, we applied a principal genome-wide network projection (PGNP) based on principal component analysis (PCA) to identify core network biomarkers in bladder carcinogenesis, with the objective of extracting the most significant part from a genome-wide network structure. Because the drug response genome-wide microarray data are now available [15], we analyzed the drug response microarray data of the core network biomarkers to design multiple drug combinations with minimal side-effects for bladder cancer treatment. Therefore, the identified core network biomarkers could provide an opportunity to design such drug combinations for bladder cancer treatment. Furthermore, it has been reported that aging (over 45 years old) and smoking are two major risk factors for bladder carcinogenesis [16]. Therefore, we used the core network biomarkers to elucidate the cellular mechanisms by which aging and smoking elevate bladder cancer risk through epigenetic regulation, miRNA regulation, and signaling pathways.
According to the strategy shown in the flowchart (Figure 1), we integrated omics data, including genome-wide methylation profiles, NGS expression data, miRNA profiles in TCGA, drug response genome-wide microarray data in the Connectivity Map (CMAP) [15], drug-gene interaction data in the Drug Gene Interaction Database (DGIdb) [17], miRNA-target gene association data in TargetScan [18], PPIs in BioGRID, transcription regulations in the Human Transcriptional Regulation Interactions database (HTRIdb) [19], the Integrated Transcription Factor Platform (ITFP) [20], and the TRANSFAC [21], biological processes and pathways in a gene ontology (GO) database, the National Center for Biotechnology Information (NCBI) Entrez Gene database, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database [22]. We used miRNA-target gene association data, PPIs, and transcription regulations to build the candidate IGEN for general molecular mechanisms. We then constructed a regression IGEN model to characterize the molecular mechanisms including miRNA regulation, PPIs, transcription regulation, and DNA methylation in cells. To prune the false positive connections in the candidate IGEN and identify the model parameters of the IGEN in the real human bladder cells, we used methylation profiles, NGS expression data, and miRNA profiles in normal bladder cells and stage 1 and stage 4 bladder cancer cells. We then applied the constrained least squares method and the Akaike information criterion (AIC) [23], a system order detection method, to prune the false positive connections for obtaining the real IGENs in the three stages of human bladder carcinogenesis. The three genome-wide real IGENs in normal bladder cells and stage 1 and stage 4 bladder cancer cells were then projected into the three core networks of the three stages of bladder carcinogenesis, respectively. Because the core networks contain the identified signal transduction pathways, that is, the receptors and TFs of the core network can be directly or indirectly connected by the core proteins/TFs, the proteins/TFs, and the corresponding genes that participate in the identified signaling pathways of the core networks are considered as the core network biomarkers for normal and cancerous cells, respectively. The miRNAs with very different connections in regulating the genes of the core network biomarkers between two cells are also involved in the core network biomarkers. By comparing the identified connections of the IGENs, we investigated how the connection changes of the core network biomarkers from normal bladder cells to stage 1 bladder cancer cells and from stage 1 bladder cancer cells to stage 4 bladder cancer cells contribute to bladder carcinogenesis.
We also investigated how the module network of the core network biomarkers, including the KEGG pathways and biological processes, participates in bladder carcinogenesis.  According to the information on the biological processes and signaling pathways in the GO database, the NCBI Entrez Gene database, and the KEGG pathway database, the roles of the TFs/proteins in the core network biomarkers are projected into three pathways: the SUMOylation, ubiquitination, and proteasome (SUP) pathway; the tumor necrosis factor (TNF) signaling pathway; and the endoplasmic reticulum (ER) signaling pathway. The roles of the downstream genes in the core network biomarkers are projected into three biological processes: cell proliferation, DNA repair, and metastasis. The module network, including the KEGG pathways, TFs, miRNAs, and biological processes, is connected according to the three identified IGENs in the three types of bladder cell. By comparing the connection changes of the module networks from normal bladder cells to stage 1 bladder cancer cells, and from stage 1 bladder cancer cells to stage 4 bladder 4 Disease Markers cancer cells, we ultimately unraveled the cellular mechanisms behind bladder carcinogenesis and proposed two multiple drug combinations for treating stage 1 and stage 4 bladder cancers, respectively. Additionally, to determine how the two major risk factors, aging and smoking, influence bladder carcinogenesis, we highlighted not only the significantly expressed genes between smokers and nonsmokers, but also the significantly expressed genes between young (≤45 years old) and old (>45 years old) people in the core network biomarkers of bladder carcinogenesis. Finally, we investigated the carcinogenic mechanism of human bladder cells by which the identified major factors, including downregulated miR-1-2, aging, and smoking, lead to the progression from normal bladder cells to stage 1 bladder cancer cells through the SUP and ER signaling pathways. The smoking-related protein HSP90AA1 and DNA methylation of ECT2 mediate the progression from stage 1 bladder cancer cells to metastasis in stage 4 bladder cancer. Activated DNA repair and accumulated epigenetic alterations lead to the phenotypic changes of bladder cells from normal to cancerous, and from cancerous to metastatic cells owing to the immortality of cancer cells. Based on the core network biomarkers in bladder carcinogenesis, a multiple drug combination comprising gefitinib, estradiol, yohimbine, and fulvestrant was designed for treating stage 1 bladder cancer with minimal side-effects, while a multiple drug combination comprising gefitinib, estradiol, chlorpromazine, and LY294002 was designed for treating stage 4 bladder cancer with minimal side-effects.

Materials and Methods
According to the flowchart in Figure 1, we constructed a candidate human IGEN by mining large databases, including BioGRID, TargetScan, HTRIdb, TRANSFAC, and ITFP. However, many false positive and insignificant connections existed in the candidate human IGEN for normal and cancerous bladder cells. Using the NGS expression data, miRNA profiles, and the methylation profiles of normal and cancerous bladder cells in TCGA, we identified the association parameters of the network connections. We also applied AIC to detect the systems order, that is, the number of connections, and to delete the insignificant connections that were out of system order to prune the false positive connections in the candidate IGEN and obtain the two real IGENs for normal and cancerous bladder cells, respectively. By applying PGNP to the two real IGENs in normal and cancerous cells, we first identified the core proteins/TFs that played a major role in the principal networks of the IGENs, constituting the core IGENs in normal and cancerous cells. To determine how the signaling cascades from the core receptor proteins to the core TFs participate in bladder carcinogenesis, the core proteins, which mediate the signal transductions from the core receptor proteins to the core TFs, and their corresponding genes were considered the core network biomarkers of the normal and cancerous cells.
The miRNAs with very different connections in regulating the genes of the core network biomarkers between normal and cancerous cells were also involved in the core network biomarkers. Finally, by comparing the connection changes of the core network biomarkers from normal cells to stage 1 cancer cells, and from stage 1 cancer cells to stage 2 cancer cells, we investigated the cellular mechanisms of bladder carcinogenesis.

Data Preprocessing of Omics Data.
We downloaded the genome-wide mRNA and miRNA NGS data and the methylation profiles from TCGA, including 17 samples for normal bladder cells, 348 samples for stage 1 bladder cancer cells, and 56 samples for stage 4, that is, metastatic stage, bladder cancer cells. The data also contained 6 samples for young (≤45 years old) people, 477 samples for old (>45 years old) people, 98 samples for nonsmokers, and 323 samples for smokers. We used one-way analysis of variance (ANOVA) to identify significant differences in gene expression between smokers and nonsmokers, and between young and old people ( value < 0.05). We used the gene symbols of the human gene information data downloaded from the NCBI FTP site as standard human gene names to integrate the omics data, including NGS data, methylation profiles, drug response genome-wide microarray data in CMAP, drug-gene interaction DGIdb data, miRNA-target gene association data in TargetScan, PPIs in BioGRID, transcription regulations in HTRIdb, and ITFP and TRANSFAC data. We also used the GO database, the NCBI Entrez Gene database, and the KEGG pathway database to find the biological processes and pathways of each gene. We used Matlab's text-file and string manipulation tools for text mining.

Construction of the Stochastic Regression Models for the IGEN System.
The goal of the stochastic regression model is to characterize molecular mechanisms, including PPIs, transcription regulations, miRNA regulations, and epigenetic regulations via DNA methylation, by NGS data through detecting false positives of candidate IGENs in human cells.
For the stochastic regression model of the gene regulatory subnetwork in the candidate human IGEN, including transcription regulations, miRNA regulations, and epigenetic regulations via DNA methylation, we identified the regulation capabilities of TFs and miRNAs in the GRN of the candidate IGEN. For the expression levels of the th gene, its DNA methylation and its th TF/protein and th miRNA in the th sample are denoted by ( ), ( ), ( ), and ( ), respectively. Then, the stochastic regression model of GRN is described by the following stochastic regression equation: ; Ω and denote the candidate regulations based on the databases of transcription regulation and miRNA-target association, respectively; indicates the regulatory ability from the th TF ( ) to the th gene; V ( ) represents the stochastic noise due to the modeling residue and fluctuation in the th gene; and , , and are the total number of TFs, miRNAs, and data samples in the omics data, respectively.
( ) denotes the effect of methylation ( ) on the binding affinity of TFs, miRNAs, or RNA polymerase on the th gene which also represents the impact of DNA methylation of the th gene on the binding affinities of miRNAs, RNA polymerase, and TFs in the gene expression process. The effect on binding affinities ( ), for = 1, . . . , , ranged between 0.2 and 1, while the expression range of the genome-wide DNA methylation ( ), for = 1, . . . , , is between 0 and 1. If DNA methylation of the th gene is close to 1, the effect on the binding affinity to the th gene is close to 0.2, which implicates the impact of DNA methylation on the binding affinities of miRNAs, RNA polymerase, and TFs to be like an inhibitor. The th mRNA expression results from transcription regulations ∑ ≡Ω , ̸ = ( ) ( ), miRNA repressions ∑ ≡ ( ) ( ) ( ), the mRNA basal expression ( ), and the stochastic noise due to measurement and random fluctuations ] ( ). In model (1), the TF regulations, miRNA regulations, and basal levels are all influenced by the DNA methylation ( ) on the th gene.
For the stochastic regression model of the miRNA regulatory subnetwork in the candidate IGEN, the expression levels of the th miRNA and its th target gene in the th sample, denoted by ( ) and ( ), respectively, could be described by the stochastic regression model of miRNA regulatory network (MRN) as the following stochastic regression equation: where the repression ability of the th miRNA to the th gene ≤ 0; the basal level of the th miRNA expression ≥ 0; ⊂ ≡ {1, . . . , }; denotes the candidate regulations of the th miRNA based on the database of miRNA-target gene association; ( ) represents the stochastic noise owing to the modeling residue and fluctuation in the th miRNA. The th miRNA expression in (2) results from miRNA-gene interactions ∑ ≡ ( ) ( ) ( ), the miRNA basal expression , and the stochastic noise ( ).
For the stochastic regression model of the PPI subnetwork in the candidate IGEN, the expression level of the th protein and its th connecting protein in th sample, denoted by ( ) and ( ), respectively, could be described by the stochastic regression model of PPIN as the following stochastic regression equation: where the basal level of the th protein expression ℎ ≥ 0; Ω ⊂ Ω ≡ {1, . . . , }; Ω denotes the candidate interactions of the th protein based on the PPI database; indicates the interaction ability of the th protein to the th protein; and ( ) represents the stochastic noise owing to the modeling residue and fluctuation in the th protein. The th protein expression in (3) results from the rate of formation of the protein complex ( ) ( ) proportional to the product of the concentration of each protein ∑ ≡Ω , ̸ = ( ) ( ) [24,25], the protein basal expression ℎ , and the stochastic noise ( ).
We proposed general stochastic regression models to characterize cellular mechanisms, including genetic and epigenetic regulations, in human cells. A number of parameters, including the TF regulatory ability , the miRNA repression ability , and the protein interaction ability , needed to be estimated and were determined using the databases of PPI, miRNA-target gene association, and transcription regulation.

Identification of the TF Regulatory Ability , the miRNA
Repression Ability , and the Protein Interaction Ability and Their Statistical Significance Testing. We used the mRNA and miRNA expression data from the NGS as the expression levels for ( ) and ( ), respectively, and used DNA methylation profiles as the expression level of ( ) to identify the model parameters , , , , , and ℎ in (1)-(3). Because large-scale measurement of protein activities has yet to be realized and 73% of the variance in protein abundance can be explained by mRNA abundance [26], mRNA expression profiles were always used to substitute for the protein expression profiles. Therefore, we also applied mRNA expression levels in the NGS data as the expression levels of ( ) to identify the parameters in (1)-(3). If the simultaneously measured genome-wide protein expression data and the mRNA expression data in each bladder cancer stage are available, the general models in (1)-(3) can also be applied to identify the real IGEN of the cancer more precisely. Because the parameters in (1) have certain constraints, such as the nonpositive miRNA repressions and nonnegative basal levels, the regulatory parameters were identified by solving the constrained least square parameter estimation problem in the following.
In order to identify the parameters in (1), the stochastic regression model of GRN was rewritten as the following linear regression form: ] where ( ) denotes the regression vector and 1 is the parameter vector of target gene to be estimated. ( ) and ( ) are available in the omics data. The regression model (4) at different data samples can be rearranged as follows: . . .
where denotes the number of data samples in the NGS data of a bladder cancer stage. For simplicity, we define the notations , Φ , and to represent (5) as follows: The constrained least square parameter estimation problem of 1 is formulated as follows: ] .
This gives the constraints to force the miRNA repression to be always nonpositive and the basal level to be always nonnegative in (1); that is, ≤ 0 and ≥ 0. The constrained least square problem was solved using the active set method for quadratic programming [27,28]. Similarly, the stochastic regression model of the miRNA regulatory subnetwork in (2) was rewritten in the following regression form: where ( ) indicates the regression vector and 2 is the parameter vector to be estimated. For simplicity, we define the notations , Ψ , and to represent (8) as follows: The parameter identification problem is then formulated as follows: ] .
This gives the constraint to force the miRNA repression to be always nonpositive and the basal level to be always nonnegative in (2); that is, ≤ 0, and ≥ 0. Finally, the protein model (3) uses the same way like above to make sure ℎ ≥ 0.

Disease Markers 7
Furthermore, in order to extract the core network biomarkers from normal and cancerous cells, we first used NGS data and methylation profiles in the normal and stage 1 and 4 bladder cancer cells to identify an IGEN for normal bladder cells and a general IGEN for bladder cancer cells. The two identified IGENs were used to extract the core network biomarkers in bladder carcinogenesis. We then used the association parameters in the general IGEN of bladder cancer cells as the initial condition of the constrained least square parameter estimation and applied the data on stage 1 and 4 bladder cancer cells to identify the IGENs for stages 1 and stage 4, respectively. According to the three identified IGENs in normal bladder cells, and the stage 1 and 4 bladder cancer cells, we determined the cellular mechanisms of the core network biomarkers in bladder carcinogenesis. The proposed methodology to identify the IGENs for normal bladder cells and stage 1 and 4 bladder cancer cells was summarized in the flowchart in Figure 2.
By applying Student's -test to the parameter estimation method [29], the values for the estimated parameters, including the TF regulatory ability , the miRNA repression ability , and the protein interaction ability , were calculated to determine the significance of the parameters. Additionally, to determine the significance of expression level and DNA methylation profile of a gene/miRNA between normal bladder cells and cancerous bladder cells, we applied one-way ANOVA to calculate the value.
After the parameter identification problem had been solved, we identified the IGEN for each bladder cell type. For example, we identified the regulatory parameter RPS20,JUN = 0.26 from the TF JUN to the target gene RPS20 ( value < 0.02) in stage 4 bladder cancer cells, the interaction parameter HUWE1,ADRM1 = 1.2 between the two proteins ADRM1 and HUWE1 ( value < 0.005) in stage 1 bladder cancer cells, and the coupling rate RPS20,MIR155 = −1.2 between the miRNA miR155 and the mRNA RPS20 in stage 4 bladder cancer cells ( value < 0.07).

Principal Genome-Wide Network Projection (PGNP).
After the identification of the IGENs in normal and cancer cells, we extracted the core network biomarkers of the IGENs based on the perspectives of the functional modules and pathways to reveal the cellular mechanisms behind bladder carcinogenesis. To extract the core network biomarkers, including the core proteins, their corresponding genes, and their upstream miRNAs, from an IGEN on a genome-wide scale, we first decomposed the combined network matrix of the IGEN to left-and rightsingular vectors and singular values based on singular value decomposition (SVD). The top left-and right-singular vectors with the top singular values constitute the principal network of the IGEN. The projection distance of each gene/protein/miRNA to these top singular vectors represents the significance of this gene/protein/miRNA in the IGEN. The genes/proteins/miRNAs with the top projection distance ultimately constitute the core network biomarkers of the IGEN. Let the combined network matrix of the TF regulatory ability , the miRNA repression ability , and the protein interaction ability of the IGEN in (1) ] .

-(3) be represented by
By applying PGP, the matrix is then be decomposed as follows: where , V ∈ R are the th left-and right-singular vectors of , respectively. The diagonal entries of are the singular values of in descending order, 1 ≥ ⋅ ⋅ ⋅ ≥ . The eigenexpression fraction ( ) is defined as We choose the top singular vectors of such that ∑ =1 ≥ 0.85, with the minimal , so that the top principal components contain 85% of the IGEN from an energy point of view. The principal projections of to the top singular vectors of , or similarities, are defined as follows: ( , ) = ⋅ V , for = 1, . . . , (2 + ) , = 1, . . . , , where and V denote the th row vector of and the th singular vector of , respectively. Furthermore, we defined the 2-norm distance from the target genes, miRNAs, and proteins/TFs to the top singular vectors, respectively, as follows: where ( ) for = 1, . . . , , for = + 1, . . . , + , and for = + + 1, . . . , 2 + are the 2-norm distances from the target genes, miRNAs, and proteins/TFs to the top singular vectors, respectively. According to ( ) for = + + 1, . . . , 2 + , we can identify the core proteins/TFs that play a major role in the principal networks of the IGENs, constituting the core IGENs in normal and cancer cells. The identified core proteins/TFs contain receptors that mediate the signaling cascades connected to core TFs. The core proteins, which participate in signal transduction from core receptors to core TFs, and their corresponding genes, were considered the core network biomarkers for normal and cancerous cells. The miRNAs with very different connections in regulating the genes of the core network biomarkers between two cells were also involved in the core network biomarkers.

Design of a Multiple Drug Combination with Minimal Side-Effects for the Treatment of Bladder Cancer.
To design a multiple drug combination with minimal sideeffects for the treatment of bladder cancer based on the core network biomarkers of the IGEN, we considered two databases, CMAP and DGIdb. CMAP contains the genomewide microarray data in response to 1327 drugs in five cell lines, while DGIdb comprises a drug-gene interaction database. Multiple drug therapy induces a genome-wide Disease Markers 9 response. The strategy of multiple drug screening is that the multiple drugs should inhibit the highly expressed genes, activate the reduced expression of the genes, and not influence the nondifferentially expressed genes in the core network biomarkers of bladder cancer cells compared with normal bladder cells. The binding protein of the designed multiple drug combination can also be obtained using the DGIdb. The strategy leads to improved drug safety and efficacy in the treatment of bladder cancer.

Construction of IGEN.
We first used NGS expression data and methylation profiles in normal bladder cells and stage 1 and 4 bladder cancer cells to identify a real IGEN for normal bladder cells and a general real IGEN for bladder cancer cells (see Section 2). By applying PGNP to the real IGEN of the normal bladder cells and the general real IGEN of the bladder cancer cells, we then obtained 115 core proteins/TFs for the core IGEN of the normal bladder cells and 138 core proteins/TFs for the core IGEN of the bladder cancer cells. To determine how the signaling cascades from the core receptor proteins to the core TFs participate in bladder carcinogenesis, the core proteins, which mediate the signal transductions from core receptor proteins to core TFs, and their corresponding genes were considered the core network biomarkers. The miRNAs with a high number of different connections regulating the genes of the core network biomarkers between normal and cancerous cells were also involved in the core network biomarkers. Moreover, to identify the mechanism of carcinogenesis from stage 1 to stage 4 bladder cancer, we used the identified parameters of models (1)-(3) in the general IGEN of bladder cancer cells as the initial condition of the constrained least square parameter estimation. We then applied the data for stage 1 and 4 bladder cancer cells to obtain the two real IGENs for stage 1 and 4 bladder cancer, respectively. Furthermore, we analyzed the connection changes of the core network biomarkers between normal bladder cells and stage 1 bladder cancer cells ( Figure 3) and between stage 1 and 4 bladder cancer cells ( Figure 4) to determine the mechanisms of bladder carcinogenesis and accordingly design multiple drug combinations for treating bladder cancer with minimal side-effects.
To investigate the impact of the major risk factors, aging and smoking, on the core network biomarkers of bladder carcinogenesis, we highlighted the significantly expressed genes between young and old people and between nonsmokers and smokers in the core network biomarkers ( value < 0.05). Additionally, the genes with changes in the basal level of (1) between normal bladder cells and stage 1 bladder cancer cells and between stage 1 and 4 bladder cancer cells were also highlighted in the core network biomarkers of Figures 3 and  4, respectively. The basal level change of a gene between two cell types has been implicated in the epigenetic regulation of gene expression. The expression of a gene that exhibits a basal level change and a significant change ( value < 0.05) of its methylation profile between the two bladder cell types is probably regulated by DNA methylation in bladder carcinogenesis.

Projection of the Core Network Biomarkers into Biological Processes and Signaling Pathways to Investigate Carcinogenic
Mechanisms of Bladder Cancer. According to the information of the biological processes and signaling pathways in the GO and KEGG pathway databases, the roles of the genes in the core network biomarkers (Figures 3 and 4) are projected into three pathways: the SUP, TNF signaling, and ER signaling pathways and three biological processes: cell proliferation, DNA repair, and metastasis.
It has been reported that the SUP pathway is associated with increased proliferation in urinary bladder carcinogenesis [30]. HuaChanSu (HCS), a class of toxic steroids, has been used to show that the TNF pathway mediates the inhibition of cell proliferation in bladder cancer [31]. Moreover, the viability of human bladder cancer cells is reduced by using cantharidin, a natural toxin, through the ER pathway [32]. Therefore, the proteins of the core network biomarkers participating in the SUP, TNF, and ER signaling pathways play an important role in bladder carcinogenesis. We then determined how the core network biomarkers mediate bladder carcinogenesis through the influences of aging, smoking, epigenetic regulation, and miRNA regulation.
The role of the SUP pathway is to degrade misfolded proteins, influence PPIs, translocate proteins, and stabilize protein structure. Owing to the accumulation of genetic mutations and epigenetic alterations in cancer cells, the SUP pathway plays a crucial role in the maintenance of many important cellular processes in cancer cells. The repressed activity of ubiquitin C (UBC), which encodes the polyubiquitin precursor, influences degradation and translation of several proteins in stage 1 and stage 4 bladder cancer cells. For example, the repression of UBC affects the signal transduction of RARRES3, a tumor suppressor, in bladder carcinogenesis. To maintain the cellular functions of cancer cells, the regulation of the SUP pathway adapts to the accumulated genetic mutations and epigenetic alterations.
In normal cells, the TNF pathway is critical for inducing inflammation, which can cause cell death. Accumulated DNA damage, epigenetic alterations, or stresses can induce the TNF pathway, and the pathway then triggers cell death. JUN, one of the TFs in the TNF pathway, plays an important role in promoting the invasion and migration of bladder cancer cells [33]. We determined that the repressed expression of JUN in stage 1 bladder cancer cells leads to cancer cell immortality and causes accumulated genetic mutations and epigenetic alterations. Additionally, the results revealed that JUN was activated in stage 4 bladder cancer cells to mediate metastasis. The role of JUN in the metastasis of bladder cancer cells can also be supported [33]. It has also been reported that the TNF pathway acts as a switch between inflammation and cancer [34]. Moreover, downregulated BCL3, which participates in the TNF pathway in the adipose tissue of the bladder wall, leads to reduced inflammation in bladder carcinogenesis [35].
The ER pathway participates in the regulation of protein folding, protein synthesis, and posttranslational modifications [36]. Misfolded proteins, arising from genetic mutations, epigenetic alterations, or stresses, induce the ER pathway to restore cellular homeostasis in normal cells. Owing to Bold lines indicate the high regulatory or interaction parameters, that is, , , and , identified in the stochastic regression models (1)-(3) of the IGEN. The bold proteins, including RARRES3, TUBA1C, PSMD8, HSPA1B, RPS20, CALR, PAAF1, and KPNA2, were the identified core network biomarkers. The major factors, including downregulated miR1-2, the aging-related proteins, HSP90B1, CALR, HSPA5, PDIA3, RPN1, and ECT2, the smoking-related proteins, HUWE1, HSPA5, and ECT2, and the epigenetic regulation of ENO1, HSP90B1, CALR, and PDIA3, lead to the progression from normal bladder cells to stage 1 bladder cancer cells through the SUP and ER signaling pathways.
the immortal nature of cancer cells, the accumulated genetic mutations and epigenetic alterations in bladder cancer cells can activate most of the genes that contribute to the ER pathway in bladder carcinogenesis (Figures 3 and 4). In the ER pathway, only RARRES3, a tumor suppressor gene, was downregulated in stage 1 bladder cancer cells.

The Impact of Aging, Smoking, and miRNA and Epigenetic
Regulation on Bladder Carcinogenesis through the Core Network Biomarkers. Major factors, including downregulated miR1-2 and aging-and smoking-related proteins, may lead to the progression from normal bladder cells to stage 1 bladder cancer cells through the SUP and ER signaling pathways.  It has been reported that aging and smoking are the major factors that accumulate genetic and epigenetic alterations and ultimately induce bladder carcinogenesis. In Figure 3, our results reveal that ADRM1 regulates KPNA2, which promotes proliferation, and is mediated by the aging-related proteins, HSP90B1, CALR, HSPA5, PDIA3, RPN1, and ECT2, the smoking-related proteins, HUWE1, HSPA5, and ECT2, and the epigenetic regulation of ENO1, HSP90B1, CALR, and PDIA3, through the SUP and ER signaling pathways. ADRM1 knockdown leads to a reduction of cancer cell proliferation and has been found in gastric [37], ovarian [38], liver [39], and colorectal cancers [40] and acute leukemia [41]. Therefore, the results support the hypothesis that aging is the most important factor in inducing bladder carcinogenesis through the SUP pathway. Additionally, our results ( Figure 3) show that the inhibited aging-related miRNA miR1-2 in stage 1 bladder cancer cells leads to miR1-2 dysregulation of genes including KPNA2, TUBA1C, HN1, PSMD11, PSMD12, and TK1, which influence cell proliferation, DNA repair, and metastasis. miR1-2 has also been identified as a tumor suppressor in bladder cancer cells [42].

miR1-2 and miR200b Mediate the Reduction of Cell
Proliferation and Metastasis through KPNA2 and ECT2, Respectively. The receptor ADRM1 signal triggers the signaling cascade from the smoking-related protein HUWE1 to the aging-related proteins HSP90B1 and RPS20 and the smoking-related TF COPS5. The TF COPS5 upregulates the metastasis-associated gene ECT2, which is suppressed by miR200b in stage 1 bladder cancer cells. The results show the cross-regulation between the transcription of the smokingrelated protein COPS5 and the aging-related protein ECT2. The aging-related miRNA miR1-2 and the smoking-related miRNA miR200b act as a switch to depress the proliferationassociated protein KPNA2 in stage 1 and stage 4 bladder cancer cells (Figures 3 and 4) and the metastasis-associated gene ECT2 in stage 4 bladder cancer cells (Figure 4), respectively.

The Smoking-Related Protein HSP90AA1 and DNA Methylation of ECT2 Mediate the Metastasis of Bladder Cancer.
Our results reveal that receptor RARRES3 signaling triggers the activated TF JUN mediated by the smoking-related protein HSP90AA1, and JUN then activates the metastasis-associated gene PSMD12 in stage 4 bladder cancer cells (Figure 4). Receptor ADRM1 signaling also triggers the metastasisassociated protein PSMD12 through the proteins PSMD8 and PAAF1 and epigenetic regulation in stage 4 bladder cancer cells. This shows that metastasis-associated ECT2 is activated by epigenetic regulation in stage 4 bladder cancer cells. Receptor RARRES3 signaling also triggers the agingrelated and proliferation-associated TF PSMD11 through the smoking-related protein HSP90AA1 in stage 4 bladder cancer cells. The activated TF JUN also regulates the proliferationassociated gene PSMD11 and the DNA repair-associated gene RPS20. There is evidence that curcumin (diferuloylmethane) can suppress tumor initiation, promotion, and metastasis. Curcumin can also inhibit the expression of JUN [43]. Additionally, the RNAi-induced induction of ECT2 suppresses cell migration, invasion, and metastasis [44]. Our results indicate that the upregulation of ECT2 in stage 4 bladder cancer cells is regulated by epigenetic regulation of ECT2 expression. This is also supported by the significant change in the DNA methylation profiles in ECT2 between normal bladder cells and stage 4 bladder cancer cells ( value < 0.007).

Functional Module Network Analysis in Bladder
Carcinogenesis. The activated DNA repair of bladder cancer cells leads to metastasis owing to the immortality of cancer cells.
According to the modular information in the GO database and the KEGG pathway database, the genes/proteins in the core network biomarkers (Figures 5 and 6) are projected into three pathways, the SUP pathway, the TNF signaling pathway, and the ER signaling pathway, and three biological processes: cell proliferation, DNA repair, and metastasis. The module networks in Figures 5 and 6 show that the activated TFs KPNA2, COPS5, PSMD12, and ECT2 play an important role in mediating the signal transduction of the SUP and ER pathways to activate cell proliferation and metastasis in stage 1 bladder cancer. The metastasis of the stage 1 bladder cancer is repressed by the activated miRNAs miR200a and miR200b, as shown in Figure 5. The activated signal transduction from SUP and ER pathways also triggers DNA repair through the epigenetically regulated TFs PSMD11 and RNF126.
Additionally, the activated TFs PSMD11, RNF126, and JUN mediate the signal transduction from SUP, TNF, and ER pathways to trigger cell proliferation, DNA repair, and metastasis in stage 4 bladder cancer, as shown in Figure 6. Although miR155 is activated in stage 1/4 bladder cancer, miR155 suppresses FOS and RPS20 in stage 1 bladder cancer, and miR155 only suppresses the DNA repair-associated gene RPS20 in stage 4 bladder cancer. Furthermore, we suggest that DNA repair may play a critical role in repairing DNA damage, which results from genetic and epigenetic alterations, leading to phenotypic change of the bladder cells from normal cells to stage 1 cancer cells, and from stage 1 cancer cells to metastatic cancer cells.
In summary, aging and epigenetic regulation dominate bladder carcinogenesis through CALR, PDIA3, DNAJB11, HSPA5, RPN1, HSP90B1, KPNA2, ECT2, and PSMD11 and through COPS5, PSMD8, RNF126, CALR, PDIA3, HSP90B1, PSMD12, PSMD11, JUN, HN1, and ENO1, respectively. Smoking promotes bladder carcinogenesis especially in metastasis. Finally, the cellular mechanisms from normal to stage 1 bladder cancer cells and from stage 1 to stage 4 bladder cancer cells are summarized in Figures 7(a) and 7(b), respectively. When the accumulated genetic mutations and epigenetic alterations lead to the dysregulation of the TNF pathway in inflammation, the accumulated misfolded proteins in the ER pathway induce cell proliferation in stage 1 bladder cancer (Figure 7(a)). Regulation of the ER and TNF pathways adapts to the accumulated genetic mutations and epigenetic alterations through the SUP pathway. The progression of DNA repair and cell proliferation in stage 1 bladder cancer ultimately results not only in the repression of miR200a and miR200b during metastasis, but also in the regulation of the TNF pathway to metastasis, cell proliferation, and DNA repair in stage 4 bladder cancer (Figure 7(b)).

Two Separate Drug Combinations for Treating Stage 1 and Stage 4 Bladder Cancer Cells with Minimal Side-Effects.
The design of a multiple drug combination for treating stage 1 bladder cancer depends on a strategy of inhibiting the highly expressed genes ADRM1, COPS5, PSMD8, SUMO2, CALR, PDIA3, DNAJB11, HSPA5, RPN1, CUL1, HSP90B1, KPNA2, PSMD12, ECT2, TK1, TUBA1C, HN1, and ENO1; activating the suppressed genes UBC, JUN, RARRES3, and FOS; and suppressing the drug's effect on the nondifferentially expressed genes BAG6, HUWE1, PAAF1, PSMD10, FAF2, PCYT1A, and PSMD10. According to the drug design strategy (see Section 2), a multiple drug combination comprising gefitinib, estradiol, yohimbine, and fulvestrant was obtained for treating stage 1 bladder cancer.   and between stage 1 and stage 4 cancer cells, we investigated how the genetic and epigenetic regulations, miRNA regulations, and aging-related and smoking-related genes affect the biological functions that lead to bladder carcinogenesis. According to gene expression changes in the core network biomarkers between normal bladder cells and stage 1 bladder cancer cells and between stage 1 and stage 4 bladder cancer cells, we then identified two separate drug combinations for treating stage 1 and 4 bladder cancer cells. Therefore, the proposed IGEN construction method and PGP provide potential network biomarkers for bladder cancer diagnosis and treatment.