This study systematically curated Vascular Dementia (VaD)-related transcriptomic datasets from the Gene Expression Omnibus (GEO) database. Inclusion criteria encompassed: (i) GSE122063 dataset (11 healthy controls and 4 VaD patient frontal and temporal cortex samples) based on the Illumina HiSeq 2500 platform (GPL16699);(ii) GSE282111 dataset (4 healthy controls and 4 VaD patient white matter samples) generated using the NovaSeq 6000 platform (GPL24676); (iii) GSE186798 dataset, comprising 30 healthy controls and 30 VaD patient frontal cortex samples, profiled using the Affymetrix Clariom S Human array (GPL23159). Raw CEL files underwent standardized parsing via the GEOquery package, with subsequent elimination of low-quality probes and batch effect-confounded data (17, 18).

A multi-stage differential analysis was implemented using the R/Bioconductor computational framework (v4.4.2). Initially, stringent quality control was applied to the raw expression matrix, excluding gene rows with >5% missing expression values and probes lacking annotation information, while retaining unique gene identifiers corresponding to maximum expression values (19). To mitigate technical variability across different microarray and sequencing platforms, batch effect correction was applied using the ComBat algorithm from the sva package—an empirical Bayes-based method extensively validated for removing non-biological variance in transcriptomic data (20). The expression matrices were then subjected to quantile normalization and log2 transformation to reduce technical noise and approximate a Gaussian distribution ( Supplementary Figure S1 ).

Differentially Expressed Genes (DEGs) were identified using the limma package (v3.58.1) to construct linear models, with significance thresholds set at adjusted p<0.05 and log2FoldChange>0.656 (corresponding to a 1.6-fold expression difference) (21). Hierarchical clustering heatmaps were generated using the ComplexHeatmap package, while volcano plots illustrating the significance-effect size distribution of differential genes were created using the ggplot2 package (22).

To systematically decipher the biological significance of DEGs, an integrated enrichment analysis strategy was adopted. The clusterProfiler package (v4.10.0) was utilized to perform three-tier Gene Ontology (GO) annotation, encompassing Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) dimensions (23) ( Supplementary Tables 1–3 ). Significant signal transduction networks (FDR<0.05) were identified by mining the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database based on a hypergeometric distribution model (24, 25). Furthermore, a global Gene Set Enrichment Analysis (GSEA) strategy was implemented using predefined gene sets from the Molecular Signatures Database (MSigDB v7.5.1) to quantify the synergistic regulatory intensity of functional modules through a weighted enrichment scoring algorithm (26).

High-confidence protein-protein interaction networks (combined score ≥0.700) were constructed based on the STRING database (v11.5), with network visualization and topological parameter calculations performed using the Cytoscape platform (v3.8.2) (27–29). Key hub proteins were identified using the Molecular Complex Detection (MCODE) algorithm, with parameters set as follows: degree cutoff = 2, score cutoff = 0.2, k-core value = 2, and max depth = 100, to screen for topologically dense functional modules (30–32). Concurrently, core regulatory factors were cross-validated by integrating Degree Centrality and Betweenness Centrality metrics (33–35).

This study employed three widely-used machine learning algorithms to identify potential biomarkers for VaD. During model development, an imbalance in sample sizes between VaD and healthy controls was noted. To mitigate bias from this class imbalance, model performance was evaluated using stratified cross-validation and out-of-bag (OOB) error estimation, ensuring balanced representation across training and validation subsets. Firstly, the Least Absolute Shrinkage and Selection Operator (LASSO) regression model was applied. A five-fold cross-validation procedure was used to optimize the regularization parameter λ by minimizing the mean prediction error. The final model, retrained on the entire dataset with the optimal λ, yielded a subset of non-zero coefficient features as initial candidate genes. Leveraging L1 regularization, LASSO enables embedded feature selection while constraining model complexity. When combined with cross-validation, this approach effectively reduces overfitting and enhances generalizability. It offers strong stability and robustness in the context of high-dimensional, small-sample omics data, thereby reinforcing the reliability of the selected biomarkers (36, 37). The optimization objective function for LASSO regression is defined as:

Where ∑ i = 1 n ( y i − β 0 − ∑ j = 1 p x i j β j ) 2 represents the residual sum of squares, λ ∑ j = 1 p | β j | is the L1 regularization term, λ denotes the regularization parameter, facilitating the selection of genes with the highest predictive value through the introduction of L1 regularization.

Secondly, the Support Vector Machine - Recursive Feature Elimination (SVM-RFE) method was utilized. A five-fold cross-validation framework was applied for both model training and feature selection, using four folds for training and one for validation in each iteration to minimize information leakage and reduce overfitting risk. In each iteration, the feature with the lowest ranking coefficient is eliminated, ultimately yielding a descending order of all feature attributes (38). The feature rankings from each fold were then aggregated to derive an averaged consensus order. To further enhance model robustness and prevent redundancy-related overfitting, classification accuracy and error rates were evaluated across varying numbers of top-ranked features, enabling the selection of an optimal feature subset.

Lastly, the Random Forest (RF) algorithm, an ensemble learning method based on decision trees, was employed. Model performance was evaluated using the Out-of-Bag (OOB) estimation strategy, wherein each tree is trained on a bootstrapped subset (63% of the data), and validated on the remaining (37%) unseen samples. This ensemble voting mechanism effectively reduces variance, enhances generalizability, and mitigates overfitting, while simultaneously enabling the assessment of feature importance (39–41). The feature importance is calculated as follows:

Where T is the total number of decision trees, N j t represents the node set split by feature in the t tree. Genes consistently identified by both LASSO and RF algorithms were prioritized as candidate biomarkers, ensuring methodological robustness and biological relevance.

By integrating the results from LASSO, SVM-RFE, and RF, feature genes commonly identified by all three methods were selected as candidate biomarkers, thereby enhancing the reliability and robustness of the screening process (42).

This study established a multi-dimensional diagnostic efficacy evaluation system, implementing Receiver Operating Characteristic (ROC) curve analysis using the pROC package (v1.18.5) in R. Through a continuous threshold scanning strategy, the classification performance of the feature gene set was systematically quantified: constructing a parametric space trajectory with sensitivity [TP/(TP+FN)] as the ordinate and 1-specificity [FP/(FP+TN)] as the abscissa (43, 44). The Area Under Curve (AUC), serving as a non-parametric statistic, was calculated based on the Mann-Whitney U test principle. The DeLong algorithm was employed to compute 95% confidence intervals (45). An AUC value approaching 1 indicates perfect discrimination capability of the feature gene set; AUC>0.9 denotes excellent diagnostic efficacy, 0.8-0.9 good, 0.7-0.8 moderate, and <0.7 suggests limited clinical applicability.

Concurrently, the study utilized the caret package to construct confusion matrices, calculating metrics such as accuracy, recall, and F1 score. An F1 score approaching 1 represents higher accuracy and recall of the model, with these two values converging numerically (46). In sample classification, this implies a more perfect overall classification effect, with all predicted results closely approximating true values, indicating enhanced reliability and effectiveness of the model in identifying and distinguishing sample categories.

To further assess gene performance in classification tasks, visualization was conducted using the ggplot2 and ggpubr packages, facilitating a more intuitive and comprehensive analysis of the genes’ diagnostic efficacy (47–49).

The CIBERSORT algorithm was employed to quantitatively assess the immune cell composition within samples. This algorithm, based on support vector regression principles, utilizes deconvolution analysis to parse mixed gene expression data into relative proportions of specific immune cell subsets ( Supplementary Table 4 ). The CIBERSORT analysis was implemented using R, with boxplots generated to visualize differences in immune cell infiltration among different sample groups (50).

To explore potential associations between feature genes and the immune microenvironment, Pearson correlation analysis was conducted. The correlation coefficient r between feature gene expression levels and the abundance of various immune cell subsets was calculated using the following formula:

Where X and Y represent gene expression values and immune cell abundances, respectively, and X ¯ and Y ¯ denote their respective means. A threshold of P < 0.05 was set for statistical significance, elucidating potential regulatory relationships between feature genes and immune infiltration atterns (51).

This study systematically analyzed the single-cell transcriptome landscape of peripheral blood mononuclear cells (PBMCs) in vascular dementia (VaD) using the GSE282111 dataset, encompassing 4 VaD patients and 4 healthy controls. A multi-stage data processing approach was implemented based on the Seurat framework (v5.0.1): Initially, the LogNormalize algorithm was applied to preprocess and normalize the raw UMI count matrix ( Supplementary Figure S2 ). Subsequently, the IntegrateLayers function was employed to execute a cross-sample integration strategy based on Canonical Correlation Analysis (CCA), mitigating technical variations between batches ( Supplementary Figure S3 ). The PercentageFeatureSet function was utilized to quantify the proportion of mitochondrial gene expression (percent.mt), with stringent selection of high-quality cells (percent.mt ≤ 5%) for downstream analysis (52) ( Supplementary Figure S4 ). Highly Variable Genes (HVGs) were identified using the FindVariableFeatures function, employing a variance stabilizing transformation (vst) strategy to select the top 2000 HVGs, capturing the primary sources of cellular heterogeneity ( Supplementary Figure S5 ). Principal Component Analysis (PCA) was performed on the HVG matrix, extracting the first 50 principal components, with dimensional significance validated using the JackStraw algorithm. A shared nearest neighbor graph was constructed using the FindNeighbors function (k=20), followed by unsupervised cell clustering using the Louvain algorithm (resolution=0.8). Finally, the FindAllMarkers function (min.pct=0.25, logfc.threshold=0.25) was employed to identify differentially expressed markers for each cell cluster, with cell type annotation performed in reference to the Cell Marker database (v2.0) (53).

Male C57BL/6 mice (8–10 weeks old, 22–25 g body weight) were obtained from SpeiPharm Biotechnology Co., Ltd. (Beijing, China) and maintained under specific pathogen-free (SPF) conditions in the Animal Science Laboratory of Inner Mongolia University for Nationalities with ad libitum access to food and water under a 12-h light/dark cycle. The bilateral common carotid artery stenosis (BCAS) mouse model was established as previously described (54), wherein microcoils (0.18 mm internal diameter) were surgically implanted to induce stenosis in both common carotid arteries (CCAs). Animals were randomly assigned using a random number table to either the BCAS group (vascular dementia model, VAD) or the control group, with six mice per group.

Cognitive function was assessed using the Morris water maze (MWM) paradigm, performed in accordance with established protocols (55, 56). The apparatus consisted of a circular pool (maintained at 24 ± 2°C) divided into four quadrants, with a 10-cm diameter escape platform submerged 1 cm below the water surface in the center of the southwest (target) quadrant. Following a 30-s habituation period, mice underwent a 6-day testing protocol comprising: (i) 4 days of directional navigation training, (ii) 1 day of navigation testing, and (iii) 1 day of spatial exploration test. All sessions were conducted at consistent diurnal timepoints. During acquisition, each mouse performed four consecutive trials daily from randomized starting quadrants (inter-trial interval: 60 min), with the platform positioned in quadrant IV. Trials were automatically terminated upon platform localization (dwell time >5 s) or after 60 s (guided to platform with 15-s rest). Navigation testing consisted of single 60-s trials initiated from quadrant II (maximal distance from target), with recording of: (i) escape latency, (ii) path length, and (iii) swimming velocity. The probe test (120-s duration) evaluated spatial memory retention by quantifying: (i) target quadrant occupancy time, and (ii) platform location crossings, following platform removal. All behavioral parameters were quantified using automated video-tracking software (EthoVision, Noldus, Netherlands).

Gene expression profiles of HSP90AA1, HSPA1B and DNAJB1 in whole blood specimens were quantitatively assessed through reverse transcription-quantitative polymerase chain reaction (RT-qPCR) methodology. Total RNA isolation was achieved using TRIzol reagent (Invitrogen), with subsequent cDNA synthesis performed employing the RevertAid First Strand cDNA Synthesis Kit (Yeasen Biotechnology) in strict adherence to the manufacturer’s protocol. Amplification reactions were conducted utilizing SYBR Green Real-Time PCR Master Mix (Yeasen Biotechnology) on a StepOne Plus Real-Time PCR System (Applied Biosystems), with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) serving as the endogenous control for data normalization. All qPCR assays were performed in technical triplicate across six biologically independent RNA preparations, with primer sequences detailed in Table 1 .

Figure 1 outlines the complete workflow. This study employed multidimensional bioinformatics analysis methods to systematically identify the characteristic transcriptome differential expression profile of vascular dementia (VaD). Following standardized data processing and statistical screening, 897 differentially expressed genes (DEGs) with significant biological relevance were ultimately identified, comprising 356 upregulated and 541 downregulated genes. A visualization analysis platform, constructed based on the R language ecosystem, generated a hierarchical clustering heatmap of the top 50 DEGs using the pheatmap package ( Figure 2B ), complemented by a volcano plot created with the ggplot2 package ( Figure 2A ).

This study systematically integrated multi-level functional enrichment analysis strategies to annotate the biological functions of the 897 identified DEGs. Gene Ontology (GO) three-dimensional classification system analysis ( Figure 2C ) revealed that in the Biological Process (BP) dimension, DEGs were significantly enriched in key processes such as T lymphocyte activation regulation (p.adjust=2.84×10^-3), positive regulation of cytokine generation (p.adjust=4.09×10^-3), leukocyte-mediated immune effector processes (p.adjust=2.71×10^-3), and hematopoietic lineage differentiation (p.adjust=2.71×10^-3). At the Cellular Component (CC) level, DEGs were specifically localized to collagen matrix (p.adjust=3.01×10^-2) and secretory granule membrane structures (p.adjust=3.01×10^-2). Molecular Function (MF) analysis unveiled their significant involvement in signal receptor activation (p.adjust=1.62×10^-2) and ligand-receptor interactions (p.adjust=1.62×10^-2). These data suggest that differential genes may participate in VaD pathological processes by regulating innate immune responses, transmembrane signal transduction, and hematopoietic stem cell-directed differentiation.

KEGG pathway enrichment analysis further elucidated the pivotal roles of these genes in molecular interaction networks, with significant enrichment in pathways such as neuroactive ligand-receptor interaction (EF≈2.58, p.adjust=6.01×10^-4), phagosome maturation (EF≈3.57, p.adjust=2.76×10^-5), and Staphylococcus aureus infection (EF≈4.55, p.adjust=1.74×10^-5). Notably, systematic exploration based on Gene Set Enrichment Analysis (GSEA) ( Figure 2D ) revealed significant clustering tendencies in GO functional sets related to immune pathways, including leukocyte activation (NES=3.84, p.adjust=1.34×10^-2) and positive regulation of immune system processes (NES=4.20, p.adjust=1.34×10^-2). The KEGG pathway set highlighted coordinated regulation of pathological processes such as osteoclast differentiation (NES=2.83, p.adjust=8.07×10^-3), alcoholic liver disease (NES=2.59, p.adjust=8.07×10^-3), and tuberculosis infection (NES=2.70, p.adjust=8.07×10^-3). This multi-level enrichment profile suggests that the onset and progression of VaD may be closely associated with cross-system mechanisms including neuro-immune axis imbalance, aberrant pathogen pattern recognition, and bone metabolism homeostasis dysregulation.

This study systematically constructed a genome-wide protein-protein interaction (PPI) network, employing multidimensional association predictions for the 897 common differentially expressed genes based on the STRING database (v11.5). Network topological feature analysis was performed using Cytoscape (v3.8.2) ( Figure 2E ). Through combined screening of topological parameters including Betweenness, Degree, and Closeness, nine core regulatory genes were successfully identified: HSP90AA1 (Degree=22), HSPA1A (Degree=7), BAG3 (Degree=9), HSPH1 (Degree=6), HSPB3 (Degree=6), DNAJB1 (Degree=10), HSPA1B (Degree=9), TLR2 (Degree=24), and HSPB1 (Degree=7). Their interaction module exhibited significant small-world network characteristics (clustering coefficient=0.608, average path length=4.23).

HSP90AA1 and HSPA1B emerged as network hub nodes (Betweenness centrality > 0.01), mediating the ubiquitination and degradation of misfolded proteins through dynamic regulation of the assembly-disassembly equilibrium of molecular chaperone complexes. Experimental evidence indicates that HSP90AA1 (Heat Shock Protein 90α subtype), an ATP-dependent molecular chaperone, participates in the cascade reactions of ischemia-reperfusion injury by modulating conformational changes of transcription factors such as NF-κB and HIF-1α. Conversely, transcriptional upregulation of HSPA1B (an HSP70 family member) can effectively eliminate abnormal protein aggregates induced by oxidative stress through enhancement of the HSP70-proteasome axis function, thereby maintaining neuronal mitochondrial homeostasis ( Figure 2F ). The topological advantage of this interaction module suggests that dysregulation of the molecular chaperone-mediated protein quality control system may be a core molecular mechanism underlying the neurodegenerative pathology in VaD.

This study systematically integrated three machine learning paradigms: LASSO (Least Absolute Shrinkage and Selection Operator) regression, Support Vector Machine-Recursive Feature Elimination (SVM-RFE), and Random Forest (RF) to conduct multidimensional feature selection on the previously identified nine key genes ( Figures 3A–F ). In the LASSO regression model, the optimal regularization parameter (λ=0.0124, SE=0.0727) was determined through 5-fold cross-validation, ultimately retaining five high-weight features: HSP90AA1 (β=1.804), HSPB3 (β=-0.140), DNAJB1 (β=2.014), HSPB1 (β=0.667), and HSPA1B (β=0.162). The SVM-RFE algorithm, based on a radial basis function kernel and optimized through five-fold cross-validation iterations, identified a subset of six features with maximum classification information entropy (HSPA1B, DNAJB1, HSP90AA1, HSPB3, HSPH1, HSPA1A). The Random Forest model (n_estimators=500) evaluated feature importance using the Gini index, determining HSP90AA1 (Gini=5.55), HSPA1A (Gini=10.89), DNAJB1 (Gini=8.05), and HSPA1B (Gini=6.84) as core predictive factors.

A systems biology evaluation framework based on inter-model feature overlap (Jaccard similarity coefficient=0.43) revealed that HSP90AA1, HSPA1B, and DNAJB1 exhibited significant enrichment trends across all three machine learning models, ultimately confirming them as robust diagnostic markers for VaD ( Figure 3G ). Notably, HSP90AA1 and HSPA1B, as core components of the molecular chaperone complex, may participate in VaD-related endoplasmic reticulum stress and mitochondrial autophagy imbalance by regulating the HSP70-HSP90 axis-mediated protein quality control mechanism. DNAJB1, a member of the HSP40 family, plays a synergistic role in the clearance of misfolded proteins by activating the ATPase activity of HSP70. Besides, Gene expression level analysis revealed significant upregulation of HSP90AA1, HSPA1B, and DNAJB1 in VaD ( Figure 3H ).

To systematically evaluate the translational value of key genes identified through machine learning in the clinical diagnosis of vascular dementia (VaD), this study employed Receiver Operating Characteristic (ROC) curve analysis to construct a molecular diagnostic efficacy assessment system, complemented by confusion matrix analysis for comprehensive evaluation ( Figure 4 ). Quantitative analysis revealed that HSP90AA1 demonstrated exceptional discriminatory performance (AUC=0.946, 95%CI:0.9035-0.9892), with sensitivity and specificity reaching 91.7% and 84%, respectively. The corresponding confusion matrix showed 33 true positives (correctly predicted VaD), 37 true negatives (correctly predicted controls), 3 false negatives, and 7 false positives, further illustrating the gene’s diagnostic performance from a classification perspective ( Figures 4A, B ). HSPA1B exhibited even more remarkable discriminatory capacity (AUC=0.951, 95%CI:0.9103-0.9925), achieving 88.9% sensitivity and 88.6% specificity simultaneously. Its confusion matrix revealed only 4 false negatives and 5 false positives ( Figures 4C, D ). DNAJB1 presented balanced detection performance (AUC=0.963, 95%CI:0.9302-0.9965), with both sensitivity and specificity exceeding 88%. Its confusion matrix differed from HSPA1B only in the number of false negatives ( Figures 4E, F ). These findings not only validate the clinical application potential of the HSP90AA1-HSPA1B-DNAJB1 molecular combination from a translational medicine perspective but also provide a theoretical basis for constructing a non-invasive diagnostic system for VaD. The three genes - HSP90AA1, HSPA1B, and DNAJB1 - exhibited remarkably robust classification performance. Their accuracy and recall rates were exceptional, with notably high F1 scores of 0.88 for all three genes. This outstanding metric indicates their significant efficacy in balancing classification accuracy and positive case identification, enabling precise sample categorization. Given the relatively small sample size in this study (8 VaD patients and 15 healthy controls), we acknowledge the potential risk of model overfitting during training and its impact on feature generalizability. To enhance the robustness and reliability of our findings, we incorporated an independent validation cohort from the GSE186798 dataset. Processing followed the standardized pipeline detailed in Section 2.1. The key heat shock protein genes demonstrated consistent classification performance across datasets, with AUCs of 0.744 (HSP90AA1), 0.728 (DNAJB1), and 0.726 (HSPA1B), underscoring their cross-cohort predictive stability ( Supplementary Figure S6 ). Based on these results, HSP90AA1, HSPA1B, and DNAJB1 show considerable potential as diagnostic and therapeutic targets for vascular dementia, offering crucial molecular-level support for clinical diagnosis and subsequent treatment strategy formulation in VaD. Table 2 presents the individual classification performance of these key genes.

This study employed the CIBERSORT algorithm to systematically characterize the immune cell heterogeneity in vascular dementia (VaD) (ns indicates no statistical difference, *p<0.05, **p<0.01, ***p<0.001). Deconvolution analysis identified six significantly differential immune cell subpopulations: naive B cells, memory B cells, plasma cells, regulatory T cells (Tregs), monocytes, and resting dendritic cells ( Figure 5A ). Immune cell interaction network analysis unveiled strong positive correlations between plasma cells and naive B cells (r=0.66), while significant negative correlations were observed between plasma cells and Tregs (r=-0.82), memory B cells (r=-0.50), and monocytes (r=-0.65). Naive B cells exhibited negative regulatory relationships with memory B cells (r=-0.57), Tregs (r=-0.73), and monocytes (r=-0.46), whereas memory B cells demonstrated positive synergistic effects with Tregs (r=0.43) and monocytes (r=0.37) ( Figure 5B ). Molecular-immune interaction analysis revealed that HSP90AA1 expression positively correlated with Tregs (r=0.51, p=1.38×10^-6), memory B cells (r=0.40, p=2.15×10^-4), and monocytes (r=0.25, p=0.0282), while negatively correlating with plasma cells (r=-0.39, p=3.32×10^-4) and naive B cells (r=-0.61, p=2.54×10^-9) ( Figures 5C and 6A–E ). HSPA1B showed positive correlations with Tregs (r=0.46, p=1.69×10^-5) and memory B cells (r=0.32, p=4.09×10^-3), and negative correlations with plasma cells (r=-0.29, p=9.97×10^-3), naive B cells (r=-0.36, p=1.11×10^-3), and resting dendritic cells (r=-0.41, p=1.59×10^-4) ( Figures 5C and 6F–J ). DNAJB1 expression patterns exhibited positive correlations with Tregs (r=0.33, p=2.83×10^-3) and memory B cells (r=0.34, p=2.25×10^-3), and negative correlations with naive B cells (r=-0.31, p=4.55×10^-3) and resting dendritic cells (r=-0.33, p=2.48×10^-3) ( Figures 5C , 6K–N ). This multi-omics data reveals that the HSP90AA1-HSPA1B-DNAJB1 axis influences the immune microenvironment through a dual regulatory mechanism: (1) inhibiting the differentiation of naive B cells into plasma cells within the B cell differentiation spectrum, and (2) promoting the immunosuppressive function of Tregs via the IL-10/TGF-β signaling pathway. These findings provide a theoretical framework for developing novel therapeutic strategies targeting the molecular chaperone-immune regulation network.

Following stringent quality control screening of eight samples from the GSE282111 dataset, we successfully annotated cell clusters using lineage-specific marker genes. High-resolution single-cell transcriptome analysis identified six biologically significant cell subpopulations ( Figures 7A, B ): S100A8+LYZ+ monocytes, SYT1+SNAP25+ neurons, CD74+CSF1R+ microglia, AQP4+GFAP+ astrocytes, FLT1+CLDN5+ endothelial cells, and MBP+MOBP+ oligodendrocytes. This cellular delineation provided a foundation for subsequent in-depth analyses. Spatial transcriptome analysis revealed cell-specific distribution patterns of heat shock protein family members. HSP90AA1 exhibited a widespread but heterogeneous distribution, predominantly enriched in oligodendrocytes (72.23%), astrocytes (7.27%), and microglia (5.15%) ( Figure 7C ). HSPA1B demonstrated a similar distribution pattern to HSP90AA1, albeit with slight variations in enrichment (oligodendrocytes 43.28%, astrocytes 5.18%, microglia 3.44%) ( Figure 7D ). DNAJB1 displayed a more specific glial preference, primarily localized in oligodendrocytes (28.04%), astrocytes (4.12%), and microglia (1.90%) ( Figure 7E ). These findings suggest potential cell-specific functions of heat shock protein family members across different neuroglial cell types (57).Differential expression analysis uncovered alterations in heat shock protein family member expression in the vascular dementia (VaD) group ( Figures 8 ). HSP90AA1 showed significant upregulation across all cell types except endothelial cells (Log2FC=0.25-0.83, q<0.001). HSPA1B exhibited a complex expression pattern: relatively stable in astrocytes and endothelial cells (Log2FC<0.5), significantly downregulated in microglia (0.78-fold, p=1.04×10^-8), and markedly upregulated in other cell types (Log2FC=0.23-2.19, q<0.001). DNAJB1’s expression pattern closely mirrored that of HSP90AA1, showing significant upregulation in all cell types except endothelial cells (Log2FC=0.52-2.22, q=0.007-0.024). These differential expression patterns reveal potential regulatory roles of heat shock protein family members in VaD pathogenesis.

Quantitative single-cell analysis ( Figure 9A ) unveiled characteristic glial subpopulation depletion in the disease group: oligodendrocytes significantly decreased by 4.56% (p<0.01), astrocytes by 0.99% (p<0.05), and microglia by 17% (p<0.001). In stark contrast, neurons exhibited pathological proliferation of 144.23% (p=0.0032), while monocytes and endothelial cells showed abnormal expansions of 7.36% and 109%, respectively (p<0.0001).

Further analysis revealed significant heterogeneity in adaptive immune cell subpopulations ( Figure 9B ). Naive B cells and plasma cells decreased by 67.36% and 60.08%, respectively (p<0.001), while memory B cells demonstrated a compensatory upregulation of 122.08% (p<0.001). Regulatory T cells (Tregs) showed a significant increase of 55.53% (p<0.001), and the monocyte subpopulation exhibited a proliferative trend of 41.98% (p<0.05). Notably, resting dendritic cells displayed a dramatic depletion of 78.42% (p<0.001), suggesting potential impairment of antigen presentation function.

To determine whether BCAS surgery successfully induced cognitive impairment in mice, behavioral assessments were conducted using the Morris water maze (MWM). The experimental paradigm comprised a 4-day directional navigation training phase (Days 1–4), followed by a navigation test on Day 5 to evaluate spatial learning ability and a spatial exploration test on Day 6 to assess spatial memory retention ( Figure 10A ). During the directional navigation training phase, the vascular dementia (VaD) group exhibited significantly prolonged escape latency compared to the control group ( Figure 10B ). On Day 5, trajectory analysis of spatial learning performance revealed that control mice navigated efficiently, adopting near-linear paths to the target, whereas VaD mice displayed circuitous and disoriented swimming patterns. On Day 6, spatial memory assessment demonstrated that control mice predominantly concentrated their search in the original platform quadrant, whereas VaD mice exhibited random exploratory behavior ( Figure 10C ). Quantitative analysis of Day 5 spatial learning performance indicated that, relative to controls, the VaD group displayed: (i) a significant increase in escape latency ( Figure 10D ) and (ii) a marked reduction in target quadrant crossings ( Figure 10E ). Similarly, Day 6 spatial memory assessment revealed that VaD mice exhibited: (i) diminished dwell time in the target quadrant ( Figure 10D ) and (ii) a significant decrease in platform crossings ( Figure 10E ). These results demonstrate that bilateral common carotid artery stenosis (BCAS) successfully induces cognitive dysfunction in murine models. Notably, no intergroup differences in swimming velocity were observed ( Figure 10F ), thereby excluding potential confounding effects of motor impairment or motivational deficits on MWM performance.

To validate our bioinformatic findings, qRT-PCR was performed to quantify mRNA expression levels of key heat shock proteins. Results demonstrated significant upregulation of HSP90AA1, HSPA1B, and DNAJB1 transcripts in the VaD group compared to controls ( Figures 11A–C ), corroborating the reliability and biological relevance of our prior data mining analyses.