The Gene Expression Omnibus (GEO) database of NCBI (https://www.ncbi.nlm.nih.gov/geo/), a freely accessible database that offers a full set of microarrays, high-throughput hybridization, chips, and other resources for research purposes, was used to collect the gene expression profile (GSE11100) of glioblastoma multiforme, which has a total of 22 samples, including from the healthy vs malignant regions of the human brain [2]. GSM280297, GSM280298, GSM280360, GSM280361, GSM280362, GSM280363, GSM280367, GSM280368, GSM280375, and GSM280376 were used as human brain healthy samples, and GSM279179, GSM279182, GSM279184, GSM279187, GSM280369, GSM280370, GSM280371, GSM280372, GSM280373, GSM280374, GSM280377, and GSM280378 were used as the human brain malignant samples. This information is based on the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2). In this study, we used a variety of analytical techniques. We first investigated differential gene expression and then used principal component analysis (PCA) to organize data with related properties. To demonstrate the genes that showed differential expression visually, we also used heatmaps. Additionally, we created a system of PPIs that was especially targeted at the DEGs. We located sub-networks within this network and highlighted significant hub genes. We developed a network that shows how proteins and drugs interact. Finally, as the final step in our thorough research, we performed molecular docking and MD simulation tests.

Using the Network Analyst web server tool (https://www.networkanalyst.ca/), the dataset was evaluated for genes displaying differential expression. The focus of this assessment was to locate genes with statistically significant expression variations (Clough and Barrett, 2016). Rows representing specific gene entries and columns representing the various samples made up the dataset’s structure. We used experimental data to determine the inclusion of 10 healthy samples and 12 malignant samples to achieve an equitable distribution. To provide uniformity and clarity throughout the dataset, all gene-probe IDs were also changed into Entrez IDs (Stacklies et al., 2007). We carried out a meta-analysis of the microarray data using NetworkAnalyst, an integrated meta-analysis web application (Xia et al., 2015). Then, we normalized each dataset using a different technique, such as the log2 transformation, quantile normalization, and variance stabilization normalization (VSN). The normality of the data was then confirmed by inspecting principal component analysis (PCA) plots and box plots (Khan et al., 2020). Using a significance threshold of p < 0.05, we carried out differential expression analysis on each individual dataset. The Benjamini–Hochberg approach was used to determine an FDR cutoff of ≥2, which we also used. In addition, a t-test using the Limma algorithm (LAT) was included in this analysis (Davis, 2007).

The initial ontology analysis of the DEGs was performed using the online bioinformatics program DAVID v6.8 (https://david.ncifcrf.gov/). This was done by tagging them with KEGG pathway information and Gene Ontology (GO) terminology, and annotations with a significance level of p < 0.05 were considered (Jiao et al., 2012). Both DEGs and hub genes were included in the modules that we uploaded. Next, we set up several characteristics, including mode-function, species (Homo sapiens), molecular function, cellular component, and ontology/pathway–biological process. Data from the Gene Ontology and KEGG databases were used as support for this analysis of functional roles and pathway enrichment (Ashburner et al., 2000; Hulsegge et al., 2009). Utilizing the Benjamini–Hochberg approach, we performed an enrichment analysis utilizing the two-sided hypergeometric test. Using a kappa score of 0.96, the analysis was conducted, and the enrichment was calculated using a threshold value greater than 0.005.

We created a PPI network including each DEG using the STRING database (https://string-db.org/) (Szklarczyk et al., 2015). The database currently comprises 67,592,464 human proteins and 20,052,394,041 documented interactions. The interactions were generated using Cytoscape_v3.10.1 and evaluated using a variety of essential variables. The protein interactions were first imported into Cytoscape version 3.10.1, and they were subsequently evaluated using the integrated tools (Chen Y.-C. et al., 2013). For the combined score to be deemed significant, it must be less than 0.75 (medium confidence score) (Oany et al., 2021). We searched for the direct interactions between each DEG and its first-level interactors. The network design was changed to only include the primary DEGs after a complete dataset with data on both DEGs was imported. Zero-order interactions were specifically chosen to enable exact PPI visualization while removing the crowded and complex network appearance commonly known as the hairball effect (Xia et al., 2014). The network can be analyzed and compared using several different topological variables. Even though Cytoscape is open-source, it has a “NetworkAnalyzer” capability that can be used to examine protein/gene networks. NetworkAnalyzer was used in this study to evaluate significant characteristics such as clustering coefficient, power law-conforming node distribution, node degree distribution, network centralization, and density. These evaluations were carried out to separate the features of the three generated networks (Ma et al., 2011).

In order to show and explore the generated PPI networks, Cytoscape version 3.10.1 (https://cytoscape.org/) was used (Shannon et al., 2003). The degree of nodes in the PPI network is determined by the number of edges to which they are connected. Hub genes are identified as nodes with high degree values (Ghafouri-Fard et al., 2023). We mapped the hub genes to examine their PPI details. Essentially, CytoHubba is a well-known integrated Cytoscape tool that analyzes the attributes and ranks the nodes correspondingly (Pan et al., 2023). It uses 11 approaches to examine the functionality of the network, including locating hub genes or network nodes. Therefore, we used CytoHubba to identify hub genes that would make viable novel therapeutic targets for glioblastoma multiforme tumor therapy.

The relationship between transcription variables and their related target genes was established using the X2K eXpression2Kinases (X2K) online application, which can be found at http://www.maayanlab.net/X2K/. This server received as input the whole set of DEGs, along with their individual gene symbols (Chen E. Y. et al., 2013). Established Fisher test p-values were used to determine the 10 most notable transcription factors (TFs) and kinases, as well as their enrichment scores. These results came from the TF and kinase module, which was built using information from the ChEA69 database to build the ChIP-X. A regulatory network was created, and Cystoscope was then used to view the “graphml” file (Brandes et al., 2001). It is ensured that the regulatory network has enough interconnected nodes throughout the network development process. If there is a gap in the pathway between kinases and TFs, the system automatically widens it to allow TFs to bind with enough intermediary proteins.

The Gene Expression Profiling Interactive Analysis (GEPIA2) (http://gepia2.cancer-pku.cn/) (Tang et al., 2017) is an extensive online platform that quickly and adaptably offers a variety of features based on information from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx). For genes that display differential expression in a particular cancer sample, GEPIA2 evaluates both the effect on survival and the analysis of expression patterns. The overall survival impact of hub genes in GBM was assessed using the GEPIA2 single-gene analysis by calculating the log-rank p-value and the 95% confidence interval-based hazard ratio (HR) (Schröder et al., 2011). The hub genes, on the other hand, were chosen according to their relative levels of expression, with a log2FC cutoff value of less than 1 and a q-value cutoff of less than 0.01.

The analysis focused on the top 10 hub genes to investigate gene–drug interactions. DrugBank version 5.1.10 (https://go.drugbank.com/), which was integrated with the NetworkAnalyzer program, provided information about the drugs and their corresponding targets. It has 16,222 medication entries in total, of which 2,751 are approved small-molecule pharmaceuticals, 1,604 are approved biologics (including proteins, peptides, vaccines, and allergens), and 134 are approved as nutraceuticals. More than 6,722 investigational medications are also now in the discovery stage. The data entry of each drug is detailed, offering a substantial quantity of information in its 200 or more data fields. This dataset includes information on a number of different topics, such as the chemical makeup and characteristics of the medicine and specifics regarding the target or targets (Wishart et al., 2018).

Prior to docking, the 2D structures of all possible compounds as inhibitors against the target protein were sketched and minimized using MM2 Force-Field of ChemOffice 2012. Under the Tripos force field (TFF) of UCSF Chimera, a monomer of the target MMP9 protein (pdb id: 5th6) was assigned with Gasteiger charges and minimized for 1,500 steps, which may be divided into 750 of conjugate and 750 of the steepest descent. Molecular docking performed via PyRx 0.9 was used. Molecular docking is an in silico approach for determining the binding modes and affinities of two molecules. In the current work, Arg56-N, a charged amino acid, was employed as a binding site residue to dock a molecule with a protein-specific site. A total of four top ligands were docked in the MMP9 protein active site. PyRx 0.9 (Dallakyan and Olson, 2015) software was used to investigate the binding free energies of the docked complexes. For each molecule, 10 solutions were prepared and screened for the best docked pose with high binding affinity. Further investigation was performed to check the binding interactions among the inhibitor moieties and the atomic residues of the protein active site. This was inferred via Discovery Studio (Biovia, 2017), where hydrogen bonds and other active bonds were observed.

First, docked complicated system preparatory simulation was carried out. For protein synthesis, AMBER16 (Case et al., 2016) force fields ff03.rl, GAFF (Sprenger et al., 2015), and ff14SB (Case et al., 2014) were utilized, and the system was solvated within a three-point transferable intermolecular potential (TIP3P) (Brice and Dominy, 2013) water box of 8.0. Minimization was used to remove unfavorable conflicts. For minimization, a conjugate gradient with 1,000 steps and the steepest decent technique were used at a value of 8. Heating was carried out for 10-ps using the Langevin dynamics method for temperature control; at a constant temperature of 300 K, 100-ps equilibration is necessary before the manufacturing run begins. The total energy remains constant during equilibration, whereas kinetic and potential energies change. The docked complex was produced in 200-ns increments, followed by equilibration. Periodic boundary conditions were suggested in the simulation box using a canonical ensemble. To keep the temperature constant, the Berendsen coupling integration approach was applied (Berendsen et al., 1984). The AMBER16 PTRAJ (Process TRAJectory) (Roe and Cheatham, 2013) module was used to generate output files for result analysis. PTRAJ was used to calculate four attributes, and graphical representations were examined in xmgrace (Srivastava et al., 2021).

The fundamental components of a biological system are genes and the products they produce, which interact randomly to form a complex network (Alamri and ul Qamar, 2023). Understanding immunology, the mechanisms of defense, signaling tasks, modes of transportation, and the onset of diseases more thoroughly can be accomplished by analyzing both genes and protein expression. Specifically, in situations where the variation is subject to fluctuations, the VSN, followed by the quantile normalization method, was used to standardize the expression data from microarray analysis. Box plots, density plots, and the data means are shown in Figure 2 both before and after the normalizing process. The figure illustrates how the normalization procedure produces constant means across all samples and successfully removes any noise that may have been present in the data. The PCA of the data revealed clear clusters. PCA was used to separate the controlled and mutant samples of the dataset into discrete clusters based on the levels of gene expression.

Several statistical tests were run on the datasets for both the control and disease state, including Student’s t-test, Benjamini–Hochberg technique, and Pearson’s correlation test. A total of 132 DEGs were discovered Utilising these investigations (Table 1). A total of 13 of these DEGs showed upregulation, while only 29 showed downregulation individually, and the rest of genes were insignificant. Using the adjusted p-value, the final ranking of the DEGs was established (Table 2). The network of PPIs is very important for understanding the workings of cellular networking. To assess how serious the issue is, changes in the protein cellular network in both healthy and diseased states are useful. In the present investigation, proteins and genes are represented by nodes, while interactions between these proteins and genes are represented by edges. We mapped these genes and created an interaction network for the identified DEGs. The mapped DEGs were shown graphically using Cytoscape. All interactions involving the DEGs were retrieved from the STRING database using a moderate confidence level of 0.40. The network created utilizing every DEG found is shown in Figure 3.

A created network consisted of 66 nodes linked by 105 edges. The network topological metrics showed that there were significant interactions within it, with an average node degree of 1.75 and an average local clustering coefficient of 0.304. Both a volcano plot and a heatmap are shown in Figure 4 to represent the DEGs.

The top 10 results of the enrichment analysis for DEGs that are upregulated and downregulated in the context of GO analysis were determined using the DAVID database, including topics such as cellular components, molecular functions, biological processes, and KEGG pathways, as described by GO. A highly significant criterion for the terms used in the analysis was the FDR. Our analysis of the GO cellular components (CCs) of the DEGs revealed a highly significant and distinctive distribution pattern for our proteins, including extracellular space, extracellular region, extracellular matrix, collagen trimer, cell surface, basement membrane, voltage-gated calcium channel complex, microtubule bundle, synapse, tertiary granule lumen, perinuclear region of the cytoplasm, GABAergic synapse, glutamatergic synapse, sarcoplasm, and voltage-gated sodium channel complex. We found 15 key functional categories for which our proteins actively participate in terms of molecular functioning (MF), such as extracellular matrix structural constituent, growth factor activity, heparin binding, carbohydrate binding, collagen binding, calcium ion binding, sequence-specific double-stranded DNA binding, identical protein binding, dehydroascorbic acid transporter activity, oxidoreductase activity, integrin binding, D-glucose transmembrane transporter activity, co-receptor binding, metalloendopeptidase inhibitor activity, and protein dimerization activity.

Furthermore, the significant KEGG pathways are primarily connected to cancer-related processes, and the reported pathways, which include ECM–receptor interaction, Notch signaling pathway, breast cancer, and endocrine resistance, are highly enriched.

The enhanced biological processes (BPs) identified in DEGs include extracellular matrix organization, cell adhesion, negative regulation of neuron differentiation, Notch signaling pathway, positive regulation of osteoblast differentiation, skeletal system development, negative regulation of cardiac muscle cell differentiation, regulation of neurogenesis, positive regulation of osteoblast proliferation, central nervous system myelination, cell–matrix adhesion, astrocyte differentiation, cell migration, inner-ear receptor stereocilium organization, positive regulation of neural precursor cell proliferation, oligodendrocyte development, positive regulation of oligodendrocyte differentiation, positive regulation of sprouting angiogenesis, nervous system development, positive regulation of angiogenesis, positive regulation of the transforming growth factor beta receptor signaling pathway, and negative regulation of the canonical Wnt signaling pathway.

The interactions of hub genes with the red nodes determined and ranked the nodes according to their degree of significant connection. Nodes having a degree value > 10 were classified as hub nodes using CytoHubba, which calculated the node degrees (Table 1). MMP9 had the highest degree, followed by 11 POSTN and Hes family BHLH transcription factor 5 (HES5) with a 9-degree score. Herein, achaete-scute family BHLH transcription factor 1 (ASCL1), collagen type-V alpha 2 (COL5A2), fibroblast growth factor 1 (FGF1), collagen type-VI alpha 2 (COL6A2), syndecan 4 (SDC4), bone morphogenetic protein 2 (BMP2), and delta-like canonical notch ligand 3 (DLL3) were among the additional genes revealed to have lower degree values. The shortest path network of these hub gene interactions showing strong interactions based on the degree of score and color is given in Figure 5.

The most well-known TFs and protein kinases connected to DEGs were identified in this analysis based on their contributions to the development of regulatory network activities. TFs, kinases, and transient proteins involved in their development were connected to create a network that served as a regulatory complex. First, we used ChIP-seq experiments (ChEA) to identify integrated target genes for transcription factors, which were then used to predict the most important TFs, which were then mapped onto PPI networks. Figures 6A,B depict the TFs predicted here and the PPI network. Based on the hypergeometric p-value, the main transcription factors in this case are SUZ12 polycomb repressive complex 2 subunit (SUZ12), enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), (tripartite motif containing 28 (TRIM28), SRY-box transcription factor 2 (SOX2), (RE1 silencing transcription factor (REST), and (SMAD family member 4 (SMAD4). We also identified and subsequently added into the PPI network kinases that are most likely to behave as regulators for the extended PPI network, as shown in Figures 6C,D. The major kinases in these DEGs identified based on the hypergeometric p-value were cyclin-dependent kinase 1 (CDK1), casein kinase 2 alpha 1 (CSNK2A1), mitogen-activated protein kinase 14 (MAPK14), cyclin-dependent kinase 2 (CDK2), glycogen synthase kinase 3 beta (GSK3B), homeodomain interacting protein kinase 2 (HIPK2), mitogen-activated protein kinase 1 (MAPK1), mitogen-activated protein kinase 3 (MAPK3), AKT serine/threonine kinase 1 (AKT1), cyclin-dependent kinase 4 (CDK4), cyclin-dependent kinase 2 (ERK2), mitogen-activated protein kinase 3 (ERK1), casein kinase 2 alpha 1 (CK2ALPHA), glycogen synthase kinase 3 beta (GSK3BETA), and component of inhibitor of nuclear factor kappa B kinase complex (IKKALPHA). All the TFs and kinases identified with their scores are given in Supplementary Table S1.

To analyze the overall impact on the survival of 10 important hub genes chosen from both the upregulated and downregulated DEGs, GEPIA survival evaluation was used. Figure 7 shows the findings of the current study. COL6A2 had a hazard ratio (HR) of 1.4 and was the only one of the 10 genes under investigation to show decreased overall survival in the group with increased expression. On the other hand, from the expression level of the hub genes, we identified that the nine hub genes MMP9, periostin (POSTN), ASCL1, COL5A2, FGF1, COL6A2, SDC4, BMP2, and DLL3 were downregulated in the normal cells, and only one hub gene, HES5, was upregulated in the tumor cells (Figure 8).

Following the identification of the top hub genes, the DrugBank database was searched for possible therapeutic candidates. In this study, we investigated each target and, using manual searches, identified the drugs related to them. The search panel included medication from a range of categories, including FDA-approved, investigational, and nutraceuticals. We manually searched for medications that interacted with each hub gene-recognized protein. In total, eight medications were found. Four drugs for MMP9, two for POSTN, and one each for COL5A2 and FGF1 were found, while no potential drugs were found for the rest of the targets. The identified drugs for all these targets are given in Table 2.

We obtained the four best compounds from the database against the target hub gene. All the compounds were docked into the target protein active site, with the coordinates of the Arg56: N atom set as the center of binding. A radius of 10 Å was set around Arg56: N, and 10 iterations were generated for each docked compound. Docked scores were used to assess chemical binding affinity for protein activity. The docked score for the complexes carmustine, lomustine, marimastat, temozolomide, and solasodine (control drug) with the MMP9 protein is −6.3, −7.4, −7.7, −8.7, and −5.7 kcal/mol, respectively. Furthermore, major residues of the protein active site implicated in hydrogen bonding, such as Glu47, Arg56, Tyr179, Pro102, Phe107, Arg106, Pro180, Asp187, Phe59, Gl57, Pro54, and Val53, were involved in binding forces, as shown in Figure 9.

Some of the residues indicated possible hydrophobic interactions. The theoretical activity of a compound requires both hydrogen and hydrophobic contacts. Furthermore, numerous additional interactions were observed, which must be considered for the computational efficacy of the compound. Thus, it has been inferred that these chemical moieties are oriented in such a way that they bind along the channel of the active pocket. Thus, it prevents the total substrate access to the active site during the metabolic pathways and inhibits the final production of cancerous substances.

Protein conformational dynamics are the most essential aspect linked with their function. The functional information of a protein molecule is encoded in its structure. To comprehend the functional variability, an understanding of the structure is required (Karplus and Kuriyan, 2005). The current work used MD simulation to investigate the conformational aspect of protein–ligand interactions and assess the stability of the inhibitor complex system in a real-time environment. Data reduction analyses such as root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF) were used to assess the conformational changes and stability index of secondary structure components of the simulated complexes. RMSD describes the backbone analysis Cα atom dynamics of the docked protein during a time scale of 100 ns, and fluctuation was observed at different time intervals among carmustine complex and marimastat complex systems, but a persisting graph of simulation stability was detected among the marimastat and temozolomide complex system. Figure 10 shows that the average RMSD value for the carmustine complex and marimastat complex system is 4.1 Å and 5.03 Å, respectively, whereas the average RMSD for the lomustine and temozolomide complex systems is 1.2 Å and 1.6 Å, respectively. To check the stability of the two carmustine complex and marimastat complex systems, simulation time intervals were extended for a further 100 ns. It was found that the system gains stability after 100 ns and remains equilibrated till 200-ns time intervals. Herein, the average RMSD recorded for both systems was 4.2 Å and 4.9 Å, respectively. It is possible that the natural flexibility of the N-terminal area is the cause of the oscillations shown in the simulation intervals. An overall structure may fluctuate because of the conformational changes introduced by the loop sections inside the N-terminal segment. These variations could be caused by the dynamic interactions the N-terminal residues with the surrounding solvent or other protein components. Furthermore, the RMSD for the protein target in Apo and complex states with the control inhibitor solasodine results in an average mean square value of 3.61 Å and 3.68 Å, respectively. Fluctuation was recorded in both Apo and complex states with a time scale of 200 ns. This overall fluctuation was then observed in the residual system of the RMSF, where the N-terminal region of the protein target showed high conformational changes. The average RMSF recorded for the Apo state was 2.8 Å. Meanwhile, for the control complex, it was measured as 3.0 Å with a time scale of 200 ns, as shown in Figure 11. The binding poses during simulation time intervals were retrieved for all the complexes at the end of the simulation, where it was inferred that strong hydrogen bond linkages were observed with a time scale of 200 ns, as shown in Figure 12. Overall, the RMSD graph pattern supports any large domain alterations within the protein–ligand complex structural framework of carmustine complex and marimastat complex systems. Herein, the inhibitors remained inside the active site during the initial time frame but moved as the time scale extended till the end of the simulation. On other hand, the ligands lomustine and temozolomide in the complex system were properly supplemented inside the binding site and did not destabilize the protein, as shown in Figure 10. Figure 13 shows RMSF values that detect the structure flexibility and volatility of Cα residues over time. A high pitch of fluctuation was recorded among the carmustine complex and marimastat complex systems. No such fluctuation was observed till 200-ns simulated time intervals, with a minor peak due to the loop region of the protein among lomustine and temozolomide in the complex system. Strong hydrogen bonds were monitored in the active zone of the MMP9 protein target. These residues include Asp187, Phe59, Gl57, and Pro54. The average RMSF of the carmustine complex and marimastat estimated from 200 ns was 1.4 Å and 1.7 Å, with a maximum peak of 4.1 Å and 6.3 Å, respectively, while substantial variations were noted at several residual points, and finally, at the conclusion of the graph for the other stable complex systems, a minor residual point was observed with an average RMSF of 1.0 Å and 0.7 Å, with no such displacement of atomic residues among the helix and sheet regions, which makes them highly stable systems. The binding pocket residues involved in RMSF include His188, Pro54, Asp187, Val53, Leu186, Gly185, Phe59, Phe62, Leu66, Ala143, Arg56, Arg51, Tyr50, Pro102, and Arg106. These residues may change as a result of ligands attaching to the active site. This implies that the flexibility or dynamics of these regions is higher. Active site residues may be involved in greater RMSF for a number of reasons; some of these residues have flexibility, whereas some are compact due to their carbon alpha background.