RNA-seq data from GSE4290 was downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/gds) database. The screening criteria were as follows: First, the dataset including GBM samples and non-tumor samples. Second, both non-tumor samples and GBM samples should be greater than 20 to ensure the quality of WGCNA. The GSE4290 dataset includes 81 samples from patients with GBM as well as 23 normal samples from healthy individuals (Sun et al., 2006). After preprocessing, we finally obtained 77 GBM samples and 23 normal samples with clinical data.

The detailed information of imipramine was obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/). The SwissTargetPrediction database (http://www.swisstargetprediction.ch/) (Daina et al., 2019) and the TargetNet database (http://targetnet.scbdd.com/home/index/) (Yao et al., 2016) as well as CTD (https://ctdbase.org/) (Davis et al., 2023), and BindingDB (https://bindingdb.org/bind) (Chen et al., 2002) were utilized for identifying potential target genes. A further search was conducted to access target proteins names using the UniProt database (https://www.uniprot.org/uploadlists/).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed by clusterProfiler R package with a threshold of adjusted p-value <0.05 (Wu et al., 2021).

The R package “limma” was used to identify the differentially expressed genes (DEGs) between the GBM samples and healthy controls (Ritchie et al., 2015). |log2 (foldchange)| > 2 and the adjusted p-value <0.05 were used as screening criteria.

GSEA is used to determine whether a defined set of genes shows a statistically significant difference in two different traits. Pathways that may play a role in the pathogenesis of GBM were explored using GSEA. The R package clusterProfiler was used for GSEA analysis. The following settings were used for the analysis: false discovery rate (FDR) < 0.25, adjusted p-value <0.05, |normalized enrichment score (NES)| > 1, which to be considered significant.

All genes were selected to construct the weighted gene co-expression network. First, a hierarchical clustering analysis was conducted to exclude the outlier samples. The Pearson correlation value was computed for each pair of genes. Then, an adjacency matrix was conducted. The WGCNA package’s “pickSoftThreshold” function was utilized to select the optimal β that satisfies the scale-free distribution (Langfelder and Horvath, 2008). The adjacency matrix was converted into two matrices: the topological overlap matrix (TOM) and 1-TOM. The TOM reflects the similarity between genes, while 1-TOM reflects the dissimilarity among them. Last, hierarchical clustering was used to classify genes into distinct modules. The module eigengene (ME) was computed to reflecting the gene expression profiles for each module. The settings of parameters were as follows. The soft threshold β = 6, minModuleSize = 50, mergeCutHeight = 0.25.

Screened genes were analyzed in STRING (version 11.5, https://string-db.org/) to investigate protein-protein interactions (PPI) (Szklarczyk et al., 2019). PPI was constructed and visualized using Cytoscape software (version 3.9.1) (Shannon et al., 2003). A Molecular Complex Detection algorithm (MCODE) was used to detect dense regions of tightly bound proteins or PPI. The MCODE algorithm is used to identify critical sub-networks that contribute to the development of GBM and to identify critical subpopulation genes. The parameters of MCODE were all set using default settings.

TISCH2 is a database that utilizes single-cell RNA sequencing (scRNA-seq) technology to concentrate on the tumor microenvironment (TME) and its immune response against tumors (Han et al., 2023).

For a single gene, TISCH2 supports exploring its expression level across different single-cell datasets or cancer types. We explored the expression level of EGFR gene in different single-cell datasets with TISCH2.

The scTIME Portal is a database and a portal for single cell transcriptomes of tumor immune microenvironment (Hong et al., 2021). We employed it to obtain the expression of pathway genes in all cell types in GSE84465. We also used scTIME Portal to perform the cell-to-cell communication analysis in GSE84465.

CAMOIP is a tool for analyzing the expression data and mutation data from the TCGA. All analyzes related to wild-type EGFR and mutant EGFR were performed on CAMOIP (http://www.camoip.net/) (Lin et al., 2022).

Firstly, the protein was optimized before docking using the Protein Preparation Wizard, and the ligand was prepared with LigPrep tool. We utilized the OPLS-2005 force field to provide partial atomic charge attribution, protonation states generation, and energy minimization. Grid-Based Ligand Docking with Energetics (Glide v11.5, Schrödinger) was used to dock all ligands into the catalytic pocket of the RPT protein in “extra precision mode” without applying any constraints. Finally, the complex with the lowest score was selected for further study. The molecular dynamics simulation protocol can be found in previous article (Wu et al., 2023).

R software (version 4.2.0) (http://www.r-project.org) was used for statistical analysis. Wilcoxon tests were used to analyze subgroup differences. p < 0.05 was used as the threshold for statistical significance.

The SwissTargetPrediction, CTD, BindingDB, and TargetNet databases were used to identify target genes associated with imipramine. Finally, 275 target genes associated with imipramine were identified and retrieved from these databases (Figure 1A; Supplementary Table S1). KEGG enrichment analysis shows that drug target genes were mainly enriched in neuroactive ligand−receptor interaction, serotonergic synapse, and calcium signaling pathway (Figure 1B). GO enrichment analysis shows that the enriched pathways were mainly involved in BP: response to xenobiotic stimulus, cellular response to catecholamine stimulus; CC: synaptic membrane, neuronal cell body; MF: neurotransmitter receptor activity, G protein−coupled amine receptor activity (Figure 1C). In total, the enrichment results showed that imipramine might be affecting neurotransmitter signaling circuits.

In GSE4290 datasets, first, we removed the batch effect between the samples (Supplementary Figure S1). According to the filtering criteria (|log2 (foldchange)| > 2 and the adjusted p-value <0.05), a total of 777 DEGs were screened, of which 207 DEGs were upregulated and 570 DEGs were downregulated (Supplementary Table S2). The result of the differential analysis was visualized in Figure 2A. Some genes participated in the pathways including Gabaergic Synapse, Synaptic Vesicle Cycle were significantly downregulated in GBM samples (Figures 2A–C). This phenomenon is consistent with previous reports (Jung et al., 2019).

GSEA analysis was performed to investigate the differences in pathway landscapes between GBM and normal samples. We noticed that Synaptic Vesicle Cycle, GABAergic synapse, Glutamatergic Synapse and cGMP−PKG Signaling Pathway were significantly downregulated in GBM group (Figures 2B–E). These results indicate synthesis and vesicle-mediated release of the neurotransmitter GABA as well as signal transduction activities may be closely related to the pathogenesis of GBM.

The GBM samples and healthy control subjects from the GSE4290 dataset were analyzed utilizing the WGCNA package in R software. The outlier samples were removed (Figure 3A). The sample dendrogram and trait heatmap was shown in Figure 3B. The cluster dendrogram was shown in Figure 3C. The matrix was transformed into an adjacency matrix using the soft-thresholding parameter, denoted as β. The R ² achieved more than 0.85 when β = 6, which became the power of our adjacency matrix (Figure 3D). After merging modules with similarity >0.75, nine co-expression modules were identified using the hierarchical clustering method of which the turquoise module exhibits the highest correlation with GBM (R = 0.78, p < 0.001), as shown in Figure 3E. Meanwhile, the scatter plot of MM relative to GS showed a good correlation between GS and MM within the turquoise module (R = 0.93, p < 0.001) (Figure 3F). Thus, the turquoise module could be a pivotal module closely associated with GBM. Additionally, the network heatmap was shown in Figure 3G.

The overlap between the DEGs and hub genes in the turquoise module was shown in Figure 4A, which was named GBM genes (Supplementary Table S3). The overlap between the GBM genes and Drug targets was shown in Figure 4B, which was named intersect genes. Then, we performed KEGG analysis on intersect genes. KEGG enrichment analysis showed that the enriched pathways were also mainly involved in Calcium signaling pathway, Gap junction, Glutamatergic synapse, Neuroactive ligand−receptor interaction, Serotonergic synapse, Glioma, GABAergic synapse (Figure 4C).

It is widely recognized that drugs elicit diverse effects on their target genes or target proteins. These effects may encompass the inhibition or upregulation of target gene expression, as well as the enhancement or suppression of target protein activity. Importantly, the precise regulatory effect of a drug on its target requires experimental validation. In this study, we solely analyze and predict the regulatory impact based on existing literature reports, with the aim of offering theoretical insights for subsequent researchers.

Traditional network pharmacology articles construct PPI with the intersect genes of drug targets and disease-related targets. However, this approach may not be entirely accurate, because there is a regulatory effect between proteins, and the drug targeting its own target protein may also affect the activity of other proteins. So, we do not only focus on the intersect genes and we selected all drug target genes and GBM genes to perform a large-scale PPI network and screen the key cluster with most interactions to explore the regulatory relationship between drug targets and disease-related targets.

A total of 883 drug target genes and GBM genes were imported into the STRING online database (version 11.5) to construct the PPI network. Then, the key cluster proteins were identified by Cytoscape’s plugin code “MCODE”. The key cluster proteins with the most interactions were visualized in Figure 5A. It is not difficult to infer that imipramine affects the progression of GBM by regulating the release and transmission of the neurotransmitter, such as GABA, thereby affecting MAPK signaling and PI3K-AKT signaling. In Figure 5A, the size of protein represents the degree of the protein, we noticed that EGFR is both a drug target and a GBM gene, while also exhibits a high degree. So, we infer EGFR as a key target for imipramine, which plays important role in MAPK signaling, PI3K-AKT signaling, Neuroactive ligand−receptor interaction, Synaptic Vesicle Cycle and GABAergic synapse.

To explore the distribution of EGFR genes in tumor microenvironment (TME), online datasets TISCH2 was utilized. The expression level of EGFR in different cells was shown in Figure 6A. Interestingly, we noticed that EGFR was expressed on neurons with a relative high level in GSE84465. In this way, the GBM single-cell RNA-seq dataset GSE84465 was chosen for further analysis. The Uniform Manifold Approximation and Projection (UMAP) plot of GSE84465 was shown in Figure 6B. The distribution of EGFR in different cell types was shown in Figure 6C. We can see that EGFR was expressed in neurons and astrocytes, implicating its potential role in the regulation of neurotransmitter.

Then, we visualize four pathways involving in neurotransmitter activity (Figures 7A–D). Genes involving in GABAergic synapse mainly distributed in fibroblast, neurons, astrocytes, cycling-T, cycling-Myeloid, macrophage, implicating these genes may play important role in immune cell fate (Figure 7A). Genes involving in Glioma pathway were also distributed in neurons, fibroblast, and macrophage, B-memory, which confirmed the importance of neurons, stromal cells and immune cell in the pathogenesis of glioma (Figure 7B). Consistently, the genes involving in Neuroactive ligand−receptor interaction were also distributed in astrocytes, neurons, fibroblast, macrophage and CD8 cells (Figure 7C). Genes involving in Synaptic Vesicle Cycle mainly distributed in neurons, astrocytes, cycling-T, cycling-Myeloid, macrophage (Figure 7D). So, the neurons might play a pivotal role in reshape TME by regulating crosstalk between tumor cell, immune cell and stromal cell.

To confirmed our conjecture, cell-to-cell communication analysis was conducted. The ligand and receptor pair strength was shown in Figures 8A–C. We noticed that EGFR gene plays important role in the communication between neurons and neoplastic cells (Figure 8A). Additionally, neurons and cycling-myeloid, neurons and macrophage also exhibit communications. Surprisingly, EGFR, COPA, GRN, TNFRSF1B were key ligands and receptors participant in communications between neurons, macrophage and neopaltic cells. So, we confirmed the key role of EGFR again, which may play an important role in regulating the crosstalk between neurons and TME. The communications between all cells were shown in Figure 8D.

The status of EGFR plays a dominant role in the therapeutic effect, which might attribute to the immune function difference. Here, we explore the landscape of immune genes in TCGA-GBM cohort (Figures 9A, B). We noticed that CD160, HVCN1, IL13, IL12A were upregulated in EGFR mutated samples. However, EGFR status might not affect the aspect of immunogenicity in GBM (Figures 10A–C). Additionally, the pathways including PD-L1 expression and PD-1 checkpoint pathway, Neuroactive ligand−receptor interaction, GABAergic Synapse and Synaptic vesicle cycle showed no difference between the EGFR-mutated patients and EGFR-wildtype patients (Figures 10D–G). Interestingly, we noticed that some immune-related pathways were downregulated in the EGFR-mutated group, including Cytokine-cytokine receptor interaction, Intestinal immune network for IgA production, Viral protein interaction with cytokine and cytokine receptor, Chemokine signaling pathway, Phagosome, IL-17 signaling pathway and Toll-like receptor signaling pathway (Figures 10H–N). So, EGFR-mutated patients might suffer a worse immune status. In addition to mutation status, EGFR expression also affects immune checkpoint-related pathways. The PD-L1 expression and PD-1 checkpoint pathway was upregulated in the patients with high EGFR expression level (Figure 10O). In conclusion, if imipramine can target EGFR and EGFR mutants, it may improve the immune microenvironment of patients.

The binding ability of the key targets was predicted using molecular docking technology. According to the general consensus, binding energies below −4.25 kcal/mol indicate specific binding between the ligand and the receptor. A binding energy below −5.0 kcal/mol indicates enhanced binding capability. Through molecular docking, we found that both EGFR (−4.291 kcal/mol) and EGFRvIII (−5.242 kcal/mol) have an excellent binding ability to the imipramine (Figures 11A, B).

The root-mean-square deviation (abbreviated as RMSD) is a measure of the deviation of atomic coordinates relative to a reference structure, and it is often used to evaluate whether a simulation system has reached stability. A stable RMSD indicates that the corresponding atoms have become stable, while a fluctuating RMSD implies fluctuations. The RMSD value of EGFR-imipramine is higher than that of EGFRvIII-imipramine, which indicates that the stability of the EGFRvIII-imipramine composite system is higher than that of the EGFR-imipramine composite system (Figure 12A).

The solvent accessible surface area (SASA) is calculated by considering the interactions between van der Waals forces and solvent molecules on the solute surface. The SASA of a protein decreases as the compactness of the protein increases, so changes in SASA can be used to predict alterations in protein structure. The fluctuations in SASA values for EGFR-imipramine are greater than those for EGFRvIII-imipramine, indicating that the stability of the EGFRvIII-imipramine composite system is higher than that of the EGFR-imipramine composite system (Figure 12B).

To investigate the interactions between proteins and ligands, we first performed a hydrogen bond analysis on the protein-ligand complexes. We found that the average number of hydrogen bonds for EGFR-imipramine and EGFRvIII-imipramine were 0.09 and 0.13, respectively. This indicates that the hydrogen bond interaction strength between EGFRvIII and imipramine is greater than that between EGFR and imipramine (Figure 12C).

To better explain the interaction energy between ligands and receptors, we used the gmx_mmpbsa script (https://jerkwin.github.io/gmxtool/) to determine the binding energies of all protein-ligand complexes during the equilibrium phase. In the application of the MMPBSA method, the total binding energy is decomposed into four independent parts (electrostatic interactions, van der Waals interactions, and polar and nonpolar solvation interactions). The results of the protein-ligand binding energies are shown in Table 1. In the protein-ligand complex systems, the main interaction energies are van der Waals and electrostatic energy. The binding free energies of EGFR and EGFRvIII proteins with imipramine are −55.882 and −64.786 kJ/mol, respectively, further indicating that the interaction strength between EGFRvIII and imipramine is greater than that between EGFR and imipramine.