We retrieved Interstitial Cystitis/Bladder Pain Syndrome (IC/BPS)-related datasets from the Gene Expression Omnibus (GEO) database and downloaded three independent public datasets (GSE11783, GSE28242, and GSE57560). To address batch effects and ensure data consistency, we first performed batch effect correction using the “SVA” R package. The datasets from each group were merged, and principal component analysis (PCA) was applied to reduce dimensionality and eliminate potential batch effects. The integrated dataset includes 31 IC/BPS samples (10 from GSE11783, 8 from GSE28242, and 13 from GSE57560) and 14 normal bladder samples (6 from GSE11783, 5 from GSE28242, and 3 from GSE57560). Subsequently, differential expression analysis was conducted on the integrated dataset using the “Limma” R package, identifying differentially expressed genes (DEGs) that met the criteria of |logFC| > 1 and p.adjust < 0.05. These DEGs were visualized as a heatmap using the “Pheatmap” package.

Immune-related genes were obtained from the GeneCards database (https://www.genecards.org/) using the screening criteria of “Protein Coding” and a relevance score > 2, yielding 8,314 related genes. A Venn diagram was employed to identify the intersection between immune-related genes and IC/BPS-specific DEGs, thereby pinpointing immune-related DEGs associated with IC/BPS.

WGCNA constructs scale-free networks by correlating gene expression levels with clinical traits. Using the “WGCNA” R package, we first removed apparent outliers from the dataset. An optimal soft threshold was selected to generate an adjacency matrix based on the topological overlap matrix (TOM). Subsequently, genes with similar expression patterns were clustered into gene modules through average linkage hierarchical clustering, using TOM-based dissimilarity measures. Finally, the modules with the most significant positive and negative correlations to clinical traits were identified. Module membership (MM) and gene significance (GS) were calculated to assess the relationships between co-expressed genes and clinical traits. Genes with higher MM and GS values exhibited stronger correlations with their respective modules and clinical traits. Applying the screening criteria of MM > 0.5 and GS > 0.5, we identified key genes within the most relevant modules.

To gain deeper insights into the functional roles of immune-related IC/BPS-DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “ClusterProfiler” R package. GO enrichment analysis, a widely used bioinformatics tool categorizes gene functional annotations into molecular functions, biological processes, and cellular components (CC). Similarly, KEGG pathway enrichment analysis has been extensively employed to elucidate biological mechanisms and functions (13). A threshold of p < 0.05 was considered statistically significant for enrichment. The results were visualized using circle plots and the “GOplot” package.

The biological significance of IC core genes was evaluated using the “GSVA” R package based on gene sets from the MSigdb. This comprehensive approach provided a detailed understanding of the functional pathways and biological processes associated with immune-related IC/BPS-DEGs, offering valuable insights into the underlying mechanisms of IC/BPS.

Based on the intersection feature genes identified through WGCNA, we employed an integrated machine learning approach to screen for core genes associated with IC/BPS. A total of 11 classical algorithms were utilized, including Lasso, Stepglm, glmBoost, SVM, Ridge, Enet, plsRglm, Random Forest, LDA, XGBoost, and NaiveBayes. The process for generating immune-related IC/BPS signatures (IRICs) was as follows: Within the test set, prediction models were fitted using a leave-one-out cross-validation framework with 113 algorithm combinations; All models were validated across three GEO cohorts (3); For each model, the Harrell’s concordance index (C-index) was calculated across all test sets and GEO datasets, and the model with the highest average C-index was selected as the optimal one. Additionally, receiver operating characteristic (ROC) curves were constructed for each diagnostic feature in the training and validation sets. The area under the curve (AUC) values were computed to evaluate the diagnostic performance of IRICs in distinguishing normal samples. To further quantify model performance, a confusion matrix was introduced, comparing predicted values with actual values to optimize the gene selection process. This rigorous approach ensured that the selected genes were both biologically meaningful and statistically reliable.

To investigate the immunological characteristics between IC/BPS samples and controls, we employed the CIBERSORT algorithm to quantify the infiltration levels of 28 immune cell types and to assess the correlations among these cells. Furthermore, Spearman’s rank correlation test was utilized to evaluate the relationships between each core gene and immune cell types in IC/BPS samples. This comprehensive analysis provided insights into the immune microenvironment and its potential role in IC/BPS pathogenesis.

In this study, single-cell RNA sequencing (scRNA-seq) data from the GSE175526 dataset were processed using the “Seurat” R package. Cells were filtered based on gene counts and the percentage of mitochondrial genes, excluding those with too few genes or excessively high mitochondrial gene content. The “NormalizeData” function was employed for data normalization, and the “FindVariableFeatures” function was used to identify the top 1,500 highly variable genes. Principal component analysis (PCA) was performed on the normalized data using the “RunPCA” function. Batch effects across different samples were corrected using the “Harmony” function, followed by dimensionality reduction and cluster identification via uniform manifold approximation and projection (UMAP). The “SingleR” package was utilized to annotate different cell clusters into specific cell subpopulations. Subsequently, the expression levels of core biomarkers across various cell groups were analyzed and visualized. Finally, cell-cell interaction analysis was conducted using the “CellChat” package, providing insights into the intercellular communication networks within the IC/BPS microenvironment.

Total RNA was extracted from mouse bladder tissues using the SPARKeasyImprovedTissue/Cell RNA Extraction Kit II (AC0202, Shandong Sparkjade Biotechnology Co., Ltd.). Reverse transcription of mRNA was performed using the Evo M-MLV Plus 1st Strand cDNA Synthesis Kit (AG11615, ACCURATE BIOTECHNOLOGY (HUNAN) CO., LTD, Changsha, China). qRT-PCR was conducted using the GeniuScript III One Step Greener RT-qPCR Kit (Q1009A, U&G Biotech). Relative expression levels were calculated using the 2-ΔΔCT method. Specifically designed primers are listed in Supplementary Table 1 .

To further explore the clinical significance of IRICs, we utilized the DSigdb to predict small-molecule drugs with potential therapeutic effects. Protein crystal structures for docking were obtained from the AlphaFold database, while the 3D structures of small molecules were downloaded from the PubChem database and energy-minimized under the MMFF94 force field. Molecular docking was performed using AutoDock Vina 1.2.3. Prior to docking, all receptor proteins were processed with PyMol 2.5.5, including the removal of water molecules, salt ions, and small molecules. The docking box was then defined using the PyMol plugin center_of_mass.py, with the box center based on the active site location and the box edge length set to 22.5 Å. Additionally, ADFRsuite 1.0 was employed to convert the processed small molecules and receptor proteins into the PDBQT format required for AutoDock Vina 1.2.5. During docking, the exhaustiveness of the global search was set to 32, while other parameters remained at their default settings. The highest-scoring docking conformation was selected as the binding conformation for subsequent molecular dynamics simulations.

Based on the small molecule-protein complexes obtained from the aforementioned docking results, all-atom molecular dynamics (MD) simulations were performed using the AMBER 22 software. Prior to the simulation, the charges of the small molecules were calculated using the antechamber module and Gaussian 09 software with the Hartree–Fock (HF) SCF/6-31G* method. The small molecules and proteins were described using the GAFF2 force field and the ff14SB protein force field, respectively. Hydrogen atoms were added to the systems using the LEaP module, and a truncated octahedral TIP3P water box with a 10 Å cutoff was added to solvate the systems. Na+/Cl- ions were incorporated to neutralize the system charge, and the topology and parameter files for the simulation were generated.

The MD simulations were conducted using AMBER 22. The systems were first energy-minimized using 2,500 steps of steepest descent followed by 2,500 steps of conjugate gradient minimization. Subsequently, the systems gradually heated from 0 K to 298.15 K over 200 ps under constant volume and controlled temperature. After reaching the target temperature, the systems were equilibrated for 500 ps under the NVT (canonical) ensemble to ensure uniform solvent distribution. Finally, the systems were equilibrated for 500 ps under the NPT (isothermal-isobaric) ensemble. The production run was performed for 100 ns under the NPT ensemble with periodic boundary conditions. During the simulation, a non-bonded cutoff of 10 Å was applied, long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method, hydrogen bond lengths were constrained using the SHAKE algorithm, and temperature control was achieved using the Langevin thermostat with a collision frequency (γ) of 2 ps⁻¹. The system pressure was maintained at 1 atm, and the integration time step was set to 2 fs. Trajectories were saved every 10 ps for subsequent analysis.

The binding free energies between the proteins and ligands in all systems were calculated using the MM/GBSA method. To ensure accuracy in MM/GBSA calculations, MD trajectories from 90 to 100 ns were used for the analysis. The specific formula for the calculation is as follows:

In Equation 1, ΔE_internal represents the internal energy, ΔE_VDW denotes the van der Waals interactions, and ΔE_elec corresponds to the electrostatic interactions. The internal energy includes contributions from bond energy (E_bond), angle energy (E_angle), and torsional energy (E_torsion). ΔG_GB and ΔG_SA collectively represent the solvation free energy, where ΔG_GB is the polar solvation free energy and ΔG_SA is the nonpolar solvation free energy. For ΔG_GB, the GB model developed by Nguyen et al. (igb = 2) was employed. The nonpolar solvation free energy (ΔG_SA) was calculated based on the product of the surface tension (γ) and the solvent-accessible surface area (SASA), as ΔG_SA = 0.0072 × ΔSASA. Due to the high computational cost and low accuracy associated with entropy calculations, entropy changes were not considered in this study.

An interstitial cystitis/bladder pain syndrome (IC/BPS) model was established in 24 female C57BL/6 mice (10-week-old) through cyclic intravesical instillation of cyclophosphamide (CYP, 35 mg/kg) and lipopolysaccharide (LPS, 0.6 mg/kg) via transurethral catheterization under 1.5% isoflurane anesthesia. Animals were randomized into four groups (n=6/group): IC/BPS model, saline-instilled control, resiniferatoxin-treated (60 μg/kg), and acetohexamide-treated (30 mg/kg) (14). Following urethral catheterization using a 24G indwelling needle with stylet guidance, chemical instillation was performed every 3 days for five cycles, with saline flushing 30 min post-administration. Two weeks post-modeling, the same regimen was used for intravesical perfusion of the therapeutic agent. Urine samples were collected during 48-h metabolic cage housing post-treatment prior to CO₂ euthanasia. Bladder tissues and blood were harvested for pathological evaluation, with inflammatory cytokines (IL-6, TNF-α) quantified using QuantiCyto^® ELISA kits (Neobioscience: EMC004(H).96.2, EMC102a.96.2). Fixed paraffin-embedded bladder sections underwent hematoxylin-eosin and immunohistochemical staining, with microscopic imaging conducted using an Olympus BX53 system. All procedures complied with institutional animal ethics guidelines (Hospital Animal Ethics Committee Approval No. 20250304829577).

All statistical analyses were performed using Perl version 5.32.1 and R software version 4.4.1. Additionally, GraphPad Prism 8.0 (GraphPad Software Inc., USA) was employed for statistical evaluation and the creation of visual data displays. The ELISA results were analyzed by one-way ANOVA. A P-value of <0.05 was considered statistically significant for all analyses.

The three microarray datasets, GSE11783, GSE28242, and GSE57560, were normalized and merged into a large training/internal validation cohort. Batch effects were corrected using the “RemoveBatchEffect” function from the “Limma” package ( Supplementary Figures 1A, B ). After normalization, the data distribution of each dataset fell within a similar range, and batch effects between datasets were effectively mitigated ( Supplementary Figures 1C, D ).

DEGs were identified between 31 IC/BPS cases and 14 normal controls using the “limma” R package. To explore immune-related genes, we integrated genes obtained from the GeneCards database and relevant literature with the identified DEGs. As a result, 73 immune-related differentially expressed genes were identified. These findings were visually represented in a heatmap ( Supplementary Figure 2 ).

To comprehensively identify key genes associated with IC/BPS phenotypes, we performed WGCNA. Hierarchical clustering of all samples was conducted, and outliers were excluded, as shown in the sample dendrogram ( Figure 1A ). A soft threshold power (β) of 6 was selected to construct a scale-free topological network ( Figure 1B ). In this study, the minimum module gene count was set to 60, and the MEDissThres was set to 0.25, ultimately identifying 11 co-expression modules ( Figure 1C ). Our results revealed that the MEbrown module exhibited the strongest positive correlation with IC/BPS in the cohort (cor = 0.43, p = 0.003, Figure 1D ). Additionally, the gene significance (GS) of the MEbrown module was significantly correlated with module membership (MM) (cor = 0.73, p < 7.6e-139, Figure 1E ). These findings suggest that genes within the MEbrown module may have functional relevance to IC/BPS. Based on the criteria of GS > 0.5 and MM > 0.5, we screened 829 key genes from the MEbrown module. By intersecting these with previously identified DEGs, we filtered out 27 Immune-related IC/BPS-DEGs ( Figure 1F ).

To elucidate the immune response mechanisms in IC/BPS, functional enrichment analysis was performed on the initially screened genes ( Supplementary Figure 3 ). GO analysis revealed that immune-related IC/BPS-DEGs were significantly enriched in biological processes such as defense response to viruses, regulation of phagocytosis, and humoral immune response, and played important roles in immunoglobulin binding and cytokine receptor activity. Additionally, KEGG pathway analysis identified significant associations with ABC transporters and antifolate resistance. These enrichment results suggest the presence of a high level of stress and abnormal immune response activation during the progression of IC/BPS.

To construct an IRIC-based diagnostic signature, we employed 113 combinations of 11 machine learning algorithms for variable selection and model development. The integrated cohort was divided into a training set and an internal validation set in a 7:3 ratio to ensure balanced distribution of clinical features. Within the internal training set, 10-fold cross-validation was performed to evaluate each algorithm combination and calculate AUC values. The ranking of AUC values for all algorithms is shown in Figure 2A . Notably, the Stepglm[both]+Enet[alpha=0.5] combination demonstrated optimal performance across both internal and external datasets, with AUC values of 0.968 (95% CI: 0.912–1.000), 0.983 (95% CI: 0.900–1.000), 1.000 (95% CI: 1.000–1.000), and 0.949 (95% CI: 0.795–1.000) ( Figure 2B ). The confusion matrix can also be seen in Figure 2C , where the overall accuracy of the model on the training set is as high as (12 + 29)/(12 + 2 + 2+29) = 91.1%. The model shows strong classification ability on all the sets with high overall prediction accuracy. Eight core genes were identified: IFI27, CSF2RB, CDC25B, IGSF3, DUSP5, WNK3, DLG, and DCDC2. The AUC of IRICs surpassed that of individual gene diagnostic features, suggesting its potential benefit for IC/BPS patients ( Figure 2D ).

The volcano plot further revealed the differential expression of the 8 core genes, including 4 up-regulated and 4 down-regulated genes. The boxplot demonstrated that IFI27, CSF2RB, CDC25B, and DUSP5 were significantly up-regulated, while IGSF3, WNK3, DLG2, and DCDC2 were down-regulated ( Figures 3A, B ). The mRNA expression levels of the 8 core genes were further validated using qRT-PCR. Compared to the control group, the expression levels of IFI27 and DUSP5 were significantly up-regulated in the IC/BPS group ( Figures 3C, F ). Conversely, the expression of IGSF3 and WNK3 was significantly down-regulated in the IC/BPS group ( Figures 3G, J ). The mRNA expression of the remaining genes showed some differences but did not reach statistical significance ( Figures 3D, E, H, I ). These results indicate that the expression levels of the 8 core genes are largely consistent with the bioinformatics analysis.

To elucidate the potential signaling pathways and mechanisms underlying IC/BPS, pathway enrichment analysis was performed on both up- and down-regulated IRICs. GSVA of IRICs indicated that upregulated pathways were primarily involved in cytokine-cytokine receptor interaction, RIG-I-like receptor signaling, and NOD-like receptor signaling, implying their potential roles in metabolic-immune crosstalk or stress responses ( Figure 4 ).

We employed the CIBERSORT algorithm to assess the correlation between the expression profiles of 22 immune cell types in IC/BPS and normal groups. Based on immune infiltration analysis and intergroup comparison boxplots, we analyzed and visualized the infiltration abundance of 22 immune cell types in both groups ( Figure 5A ). Significant differences were observed in the expression of T cells CD4 memory activated, T cells follicular helper, and Macrophages M0 between the two groups. Subsequently, the correlation between each core gene and immune cells was displayed in a correlation network heatmap ( Figure 5B ). Notably, the IFI27 gene showed significant positive correlations with Macrophages M0, Macrophages M1, T cells CD4 memory activated, T cells follicular helper, and T cells regulatory (Tregs), while exhibiting a negative correlation with Macrophages M2 ( Figures 5C, D ).

To comprehensively characterize the feature genes, we performed single-cell analysis using the GSE175526 dataset ( Figure 6A ). The data were clustered using the “Seurat” R package. Visualization via UMAP identified 10 distinct cell types, including endothelial cells, epithelial cells, fibroblasts, monocytes, neuronal cells, hematopoietic stem cells, adipocytes, neutrophils, B cells, and CD8+ T cells ( Figure 6B ). The expression patterns of each core gene in different cell types were then investigated, and t-SNE visualization showed that IFI27 exhibited dense expression in most cell types ( Figure 6C ). Dot plot analysis further demonstrated that IFI27 was highly expressed in epithelial cells, endothelial cells, and fibroblasts, with partial expression in monocytes and CD8+ T cells ( Figure 6D ). To explore biological feature differences, we investigated complex communication networks among the 10 annotated cell types. The results indicated increased communication strength between CD8+ T cells (as ligand cells) and monocytes as well as neutrophils. Monocytes also exhibited strong communication with endothelial cells, neutrophils, and B cells ( Figure 6E ). Additionally, the number of interactions between fibroblasts and other cell types was notably increased ( Figure 6F ).

Molecular docking simulation is a convenient and effective approach for exploring the interactions between small molecules and target proteins. Here, we used Vina 1.2.3 software to perform docking studies on the compounds Acetohexamide, Resiniferatoxin, and the IFI27 protein. As shown in Figure 7A , hydrogen bond interactions were observed between the small molecules and THR-103 on the protein, while hydrophobic interactions occurred with ILE-106, LEU-99, VAL-41, ILE-33, and VAL-36. The formation of hydrogen bonds and hydrophobic interactions enhanced the binding affinity between the protein and small molecules. The interactions of the IFI27_Acetohexamide complex are illustrated in Figure 7B , revealing hydrogen bond formation with THR-103 and hydrophobic interactions with ILE-106, LEU-99, VAL-41, MET-44, and VAL-18 on the protein.

The root mean square deviation (RMSD) in molecular dynamics simulations reflects the motion of complexes, where higher RMSD values and more pronounced fluctuations indicate greater motion, and vice versa. As shown in Figure 7C , the RMSD changes of the IFI27/Acetohexamide and IFI27/Resiniferatoxin complexes during the simulation are presented. During the simulation, the RMSD of the IFI27/Acetohexamide complex increased rapidly in the 0–20 ns range and then stabilized, maintaining fluctuations within 7–10 Å throughout the simulation. In contrast, the RMSD of the IFI27/Resiniferatoxin complex started from a lower initial value, continued to rise within the first 40 ns, and eventually stabilized at approximately 8 Å. The two complexes exhibited different RMSD trends during the dynamics process, indicating differences in their dynamic stability. Although both complexes showed relatively stable RMSD values in the later stages, the IFI27/Resiniferatoxin complex exhibited significantly smaller fluctuations, suggesting higher structural stability.

The root mean square fluctuation (RMSF) reflects the flexibility of proteins during molecular dynamics simulations. Typically, drug binding reduces protein flexibility, thereby stabilizing the protein and enhancing enzymatic activity. As shown in Figure 7D , the RMSF values of all proteins bound to different small molecules were low (<5 Å), except at the terminal regions, indicating good rigidity of the protein core structure. In summary, the low flexibility of the protein is the foundation for stable binding between small molecules and the protein.

Based on the molecular dynamics simulation trajectories, we calculated the binding energy using the MM-GBSA method, which more accurately reflects the binding modes of small molecules to the target protein ( Figure 7E ). As shown in Table 1 , the binding energies of the IFI27/Acetohexamide and IFI27/Resiniferatoxin complexes were −19.91 ± 0.98 and -32.98 ± 1.74 kcal/mol, respectively. Negative values indicate binding affinity between the molecules and the target protein, with lower values indicating stronger binding. Our calculations clearly demonstrate that these molecules have a binding affinity with the corresponding proteins, with IFI27/Resiniferatoxin exhibiting higher binding energy. The binding energy of these complexes is primarily contributed by van der Waals and electrostatic interactions, while solvation energy is unfavorable for their binding.

Hydrogen bonds are one of the strongest non-covalent interactions, and a higher number of hydrogen bonds indicates better binding. As shown in Figure 7F , the IFI27/Acetohexamide complex had fewer hydrogen bonds (around 1) during the simulation, indicating a weak contribution of hydrogen bonds to their binding. In contrast, the IFI27/Resiniferatoxin complex had around 2 hydrogen bonds during the simulation, suggesting that hydrogen bonds play a significant role in their binding.

To validate the therapeutic efficacy of IFI27 and candidate drugs, histopathological evaluation of IC/BPS models revealed significant inflammatory cell infiltration within the lamina propria, accompanied by mild mucosal edema and partial epithelial denudation characterized by intercellular space widening ( Figure 8A ). Resiniferatoxin treatment marginally attenuated inflammatory infiltration, whereas acetohexamide markedly reduced mucosal edema and inflammatory infiltration. Mechanistically, upregulated IFI27 expression in IC/BPS models confirmed the pivotal role of innate immune activation, with immunohistochemical analysis demonstrating acetohexamide’s superior efficacy in downregulating IFI27 compared to resiniferatoxin ( Figure 8B ). Given the characteristic immune dysregulation in IC/BPS, cytokine profiling via ELISA revealed significantly elevated IL-6 and TNF-α levels in serum and bladder tissues of model groups. Both therapeutics exhibited anti-inflammatory effects, yet acetohexamide demonstrated significantly stronger suppression of proinflammatory cytokines than resiniferatoxin ( Figures 9A, B ).