Participants diagnosed with IS at the First Affiliated Hospital of Guangxi Medical University from January 2019 to March 2025 were included in the study. IS was diagnosed according to the American Heart Association and American Stroke Association criteria (13) and confirmed by MRI or brain computed tomography. Acute/sub-acute (AIS) and chronic IS (CIS) were identified according to characteristic imaging features observed across various MRI sequences (14). Exclusion criteria were as follows: (a) severe aphasia and/or impaired consciousness that compromised the ability to understand and follow instructions; (b) a history of major systemic diseases, including severe cardiac, hepatic, or renal dysfunction, autoimmune disorders, or other psychiatric illnesses; (c) coexisting severe neurological disorders, such as brain tumors or traumatic brain injury, which could potentially confound the study outcomes; (d) contraindications to MRI scanning, or poor image quality that did not meet the requirements for analysis. The healthy control (HC) group was composed of 51 individuals with comparable demographic characteristics to the patient group and had received a thorough medical examination. The study procedures received approval from the Ethics Committee of the First Affiliated Hospital of Guangxi Medical University (2025-E0302), and informed consent was obtained from each participant.

Imaging data were acquired using a 3 T MRI scanner (SIEMENS, Germany). High-resolution T1-weighted anatomical images were acquired using a magnetization-prepared rapid gradient-echo (MPRAGE) sequence with the following parameters: matrix size = 256 × 256, field of view = 256 × 256 mm², voxel size = 1 × 1 × 1 mm³, slice thickness = 1.0 mm with no inter-slice gap, echo time = 2.98 ms, repetition time = 2,300 ms, flip angle = 9°, and a total of 176 axial slices.

The three-dimensional T1-weighted images were preprocessed in surface-based space using the Computational Anatomy Toolbox (CAT12)¹ implemented in the Statistical Parametric Mapping software (SPM12)² framework. Briefly, preprocessing included skull stripping, tissue segmentation into gray matter, white matter, and cerebrospinal fluid, as well as estimation and reconstruction of the cortical surface. The cortical thickness was computed using the projection-based thickness method, which involves correction for topological defects, spherical mapping, and registration to a common surface template (15). In addition, CAT 12 was also used to estimate other morphological metrics, including fractal dimension, gyrification index, central surface, mean curvature, square root of sulcal depth, and cortical thickness (8, 16). We calculated these indices for each participant using default parameter settings. The morphometric maps were resampled to the fasverage template and smoothed using a Gaussian kernel with full width at half maximum (FWHM) after they were obtained. To be more precise, we took the cortical thickness maps for each person and made them smoother using something called a Gaussian kernel. The width of this kernel was 15 mm at full width. Other maps (such as fractal dimension maps) were made smooth using a Gaussian kernel with a FWHM of 25 mm, as explained by Ruan et al. (17).

We utilized the methods of Seidlitz et al. (5) to generate MSNs with the T1-weighted restricted parameters by using CAT 12. To be more specific, the outer surfaces of the brain were divided into 308 sections, of roughly equivalent dimensions (~500 mm²), by using a backtracking algorithm (18). These sections were based on 68 sections of the brain as shown in the Desikan–Killiany atlas (19). This Desikan–Killiany parcellated atlas was transformed onto each participant’s surface to obtain an individual surface parcellation. Six morphometric features (as detailed in the Imaging processing section) extracted previously from the T1-weighted images were used for each region. To take into account differences in the values of the features, each MRI feature in each region of each image was z normalize. Pearson correlation coefficients were computed between the morphometric feature vectors of each pair of regions, thereby generating a 308 × 308 node morphological similarity matrix for each subject. This matrix was thresholded in order to construct a weighted bipartite graph with an arbitrary connection density. This process is also known as the construction of MSN. Regional MSN strength was defined as the average weighted correlation coefficient between a given region and all other regions, serving as an index of region-specific morphometric connectivity.

To assess group differences in MSN strength, we employed a linear regression model (LRM), with MSN strength as the dependent variable and group (AIS, CIS, or HC) as the main independent variable. Age, sex, and their interaction (age × sex) were included as covariates. This model was applied at both the global level (average MSN strength across all regions) and the regional level (MSNᵢ for each cortical region). For regional analyses, the full model was specified as: MSNᵢ = intercept + β₁ × age + β₂ × sex + β₃ × (age × sex) + β4 × group. Here, the group variable was dummy-coded, and β₄ represents the group effect (i.e., AIS vs. HC or CIS vs. HC). This model was independently fitted for each of the 308 cortical regions. Two-sided t-tests were used to assess statistical significance for group comparisons, and multiple comparisons were corrected using the Benjamini-Hochberg false discovery rate (FDR) method (p < 0.05). Regions showing significant differences in MSN strength were further visualized using brain surface maps. Spatial correlations between case–control t-values and healthy control MSN strength were assessed using spin permutation tests to account for spatial autocorrelation. The spin test is a spatial permutation method. It uses random rotations of spherical projections at the cortical surface. This method is more conservative than randomly shuffling locations. It preserves the data’s spatial symmetry and contiguity. Specifically, the initial procedure involved the generation of a null distribution through the implementation of 10,000 random spatial rotations of the cortical parcellation. Subsequently, the p-spin values were derived from the comparison to the null models.

To demonstrate the comparison of MSN strength between the HC group and the stroke subgroups (AIS and CIS), We used the Pearson correlation coefficient to compare MSN strength between the HC group and the stroke subgroups (AIS and CIS). The following terms were defined, similar to the descriptions by Dong D et al. (20). Decoupling: Positive correlation in the HC group and negative correlation in the stroke groups indicates that structurally similar regions become less similar. Dedifferentiation: Negative correlation in the HC group and positive correlation in the stroke groups indicates that structurally dissimilar regions become similar. Hypercoupling: Positive correlation in both the HC group and the stroke groups but stronger correlations in the stroke groups indicates that similar regions become more similar. Hyperdedifferentiation: Negative correlation in both the HC group and the stroke groups but weaker correlations in the stroke groups indicates that dissimilar regions become more dissimilar.

To investigate the associations between MSN features and clinical severity, we conducted correlation analyses between MSN strength and clinical scores, including the National Institutes of Health Stroke Scale (NIHSS) and modified Rankin Scale (mRS), for both AIS and CIS groups. Specifically, we examined: Global-level associations, by computing Spearman’s rank correlation coefficients between mean MSN strength across all 308 cortical regions and clinical scores; Regional-level associations, by calculating Spearman correlations between MSN strength in each cortical region (defined by the D-K308 atlas) and clinical scores. For the two subgroups, AIS and CIS, all relevant analyses were conducted respectively, and the Benjamini-Hochberg FDR correction method was adopted to control multiple comparisons. The correlation results were considered statistically significant when the p value was less than 0.05 (the p value after FDR correction).

To explore the molecular basis underlying cortical morphometric similarity, we utilized transcriptomic data from the AHBA,³ which provides genome-wide microarray-based gene expression profiles sampled from 3,702 anatomically localized sites across six adult postmortem donors (mean age = 42.5 ± 13.4 years; 5 males, 1 female) (10). In line with established preprocessing pipelines, we employed the Abagen toolbox⁴ to map and preprocess gene expression data with respect to the D-K308 region parcellation scheme (12, 21). The preprocessing procedure included the following steps: (a) updating probe annotations to gene-level identifiers; (b) excluding probes with expression levels below the background threshold in more than 50% of the samples; (c) retaining probes with the greatest spatial consistency for each gene; (d) assigning each sample to a cortical parcel based on a 2 mm Euclidean distance criterion; and (e) applying robust sigmoid normalization to expression values within donors. As the AHBA data comprises right hemisphere samples from merely two donors, our analysis was restricted to the left hemisphere regions (n = 152). Although output files retained the “DK308” label for consistency with the parcellation scheme, all subsequent analyses were strictly limited to 152 left-hemisphere regions. This preprocessing resulted in a gene expression matrix comprising 152 brain regions and 15,632 transcripts. To maintain whole-genome coverage, we did not further filter low-expression or low-variance genes; instead, robustness was ensured through subsequent bootstrapping and false discovery rate (FDR) correction.

To explore the molecular correlates of regional MSN differences across the course of ischemic stroke, we employed PLS regression to link spatial patterns of MSN change with transcriptomic data (22). Specifically, we extracted t-statistics representing group-level differences in MSN features between acute/subacute and chronic stroke patients across 152 cortical regions of the left hemisphere. These t-values served as the response variable, while gene expression levels derived from the AHBA were used as predictors. Expression data were matched to the D-K308 parcellation, restricted to the left hemisphere regions for consistency.

The number of PLS components was determined by examining the cumulative variance explained across the first 15 components. Models including 1–15 components were constructed using the plsregress function, and the percentage of variance explained in the dependent variable was recorded for each component. A null distribution was then established through 1,000 permutations, in which the dependent variable was randomly shuffled and PLS was recalculated 1,000 times. The variance explained by each component was obtained. In the original data, the first component significantly exceeded the random level (p = 0.001), while subsequent components contributed less than 2% incremental variance. Therefore, only the first two components were able to account for a substantial proportion of variance in the response variable and exceeded the permutation-based null distribution, which led us to focus on PLS1 and PLS2. Ultimately, the first partial least squares component (PLS1), representing a weighted combination of gene expression most strongly associated with regional MSN differences in ischemic stroke, was retained for further analysis.

In addition, we estimated the stability of each gene’s loading on PLS1 using bootstrap resampling. For each gene, a Z score was calculated as the ratio of its mean weight to the bootstrap-derived standard deviation. Although the calculation of Z scores assumes approximate normality, bootstrap resampling provides an empirical estimate of weight variability, thereby reducing over-reliance on distributional assumptions. To address the inherent sign indeterminacy of PLS, bootstrap weights were aligned to the direction of the original PLS component to ensure consistency across resampling iterations. Since the sign of PLS components is mathematically arbitrary, our interpretation focused on the relative ranking of gene weights rather than their absolute direction. Confidence intervals were not computed for ROI-level PLS scores. Instead, uncertainty was primarily quantified at the gene-weight level using bootstrap standard errors, Z scores, and FDR-corrected p-values. The associations between ROI-level PLS scores and MSN t-values were further validated using spin permutation tests to account for spatial autocorrelation.

Genes were then ranked according to their Z scores, and statistical significance was determined using FDR correction with a threshold of p < 0.05. Genes with PLS1 Z scores greater than 3.0 were designated as PLS1+, while those with Z scores less than −3.0 were classified as PLS1−. Only genes significantly associated with PLS1 were included in subsequent spatial expression mapping and enrichment analyses. PLS1 scores were also correlated with MSN t-values to validate the explanatory power of the PLS model.

To explore the potential biological functions of genes associated with structural network differences in ischemic stroke, we conducted GO and pathway enrichment analyses on the ranked gene lists derived from PLS regression. Specifically, genes exhibiting statistically significant differences in bootstrap-derived Z scores for each condition (acute/subacute and chronic stroke) were analyzed separately to identify enriched biological processes and signaling pathways. Furthermore, to pinpoint genes that may play a role in both early and late stages of stroke, we performed an intersection of PLS1 genes from the two time points, followed by functional enrichment analyses on the overlapping set. All results were plotted by the Bioinformatics online platform (last accessed April 1, 2025)⁵, which provides integrated tools for GO term classification and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (23). Statistically significant terms were determined with a threshold of p < 0.05.

A total of 170 patients with IS, including 60 in the AIS and 59 in the CIS, along with 51 HCs, were included in the final analysis. The demographic and vascular risk profiles of all participants are presented in Table 1. No statistically significant differences were observed in age (AIS: 55.98 ± 8.91 years; CIS: 56.41 ± 11.43 years; HC: 54.06 ± 8.91 years) or gender distribution across the groups (p > 0.05). Similarly, comparisons among groups revealed no significant differences in the prevalence of major vascular risk factors, including smoking, drinking, hypertension, diabetes mellitus, dyslipidemia, coronary heart disease, and atrial fibrillation (all p > 0.05). Clinical severity, as measured by the NIHSS and mRS, showed a trend toward higher scores in the AIS group compared to the CIS group, although these differences did not reach statistical significance (NIHSS: p = 0.053; mRS: p = 0.060). These findings suggest comparable baseline demographic and risk factor distributions among groups, minimizing potential confounding effects in subsequent morphometric and transcriptomic analyses.

We first constructed individual MSNs by computing interregional Pearson correlations across six cortical features derived from T1-weighted MRI data. At the group level, no statistically significant differences in global MSN strength were observed between AIS or CIS patients and healthy controls (Figure 1A). Therefore, we further examined region-specific differences in MSN to identify localized changes associated with stroke. Regional case–control t-maps (AIS-HC and CIS-HC) were generated by conducting linear regression analyses at each cortical region, controlling for age and sex. For each comparison, two-sided t-statistics were computed to quantify differences in regional MSN strength between stroke patients and healthy controls. Positive t-values indicate regions with increased MSN in AIS or CIS relative to controls, while negative t-values reflect decreased MSN. The spatial patterns of these t-values revealed distinct region-specific disruptions in morphometric similarity associated with the acute/subacute and chronic phases of IS. In the AIS group, significant differences in MSN strength relative to HC were predominantly observed in the frontal and temporal cortices, with additional involvement of parietal, occipital, and limbic regions (Figure 1B; Supplementary Table S1). Notably, regions showing reduced MSN strength included bilateral superior and middle frontal areas, the precentral gyrus, middle temporal cortex, and lingual gyrus. Conversely, increased MSN strength was identified in regions such as the temporal pole, insula, entorhinal cortex, cingulate cortex, and inferior parietal lobule. These findings suggest widespread and regionally divergent MSN difference during the acute/subacute phase of IS. In the CIS group, widespread differences in MSN strength were observed across multiple cortical regions compared to HC (Figure 1C; Supplementary Table S1). Reduced MSN strength was mainly observed in the bilateral prefrontal cortex, including the caudal middle frontal, rostral anterior cingulate, superior frontal, orbitofrontal, and medial prefrontal areas, as well as in the precentral gyrus, posterior cingulate cortex, parietal regions (including precuneus, supramarginal, inferior parietal), and occipital regions (including lingual and lateral occipital gyri). Increased MSN strength was primarily located in right hemisphere regions, including the caudal middle frontal, fusiform, entorhinal, insula, inferior parietal, and lateral orbitofrontal cortices, as well as the frontal pole and supramarginal gyrus. These results indicate that distinct MSN patterns are observed in the chronic phase, characterized by both enduring disruptions and regionally elevated similarity in specific cortical hubs. In both AIS and CIS groups, altered MSN connections were predominantly characterized by decoupling (39%), followed by dedifferentiation (AIS: 25%, CIS: 26%), hypercoupling (AIS: 28%, CIS: 21%), and hyperdedifferentiation (AIS: 8%, CIS: 7%) (Figure 1D). A significant negative correlation was observed between MSN strength in healthy controls and group differences (AIS: Spearman’s r = −0.35, p-spin < 0.001; CIS: Spearman’s r = −0.58, p-spin < 0.001), indicating that in the comparison between the HC group and the stroke patient group (AIS and CIS), areas with higher MSN strength in the HC group often exhibit lower MSN strength in the stroke patient group, and conversely. This negative correlation suggests that regional differences in MSN strength between the HC group and the stroke patient group are negatively correlated with MSN strength itself in the HC group. Most significantly changed connections were located in the decoupling and dedifferentiation quadrants, reflecting disrupted integration of structurally coherent regions following stroke.

To explore the associations between MSN features and clinical severity, we performed Spearman correlation analyses between both global mean MSN strength and regional MSN values (based on D-K308 atlas) and clinical scores (NIHSS and mRS), separately for the AIS and CIS groups. After applying FDR correction (p > 0.05, Supplementary Table S2), no statistically significant correlations were observed between global mean MSN strength and either NIHSS or mRS scores in the AIS or CIS groups. At the regional level, prior to FDR correction, region-specific associations between MSN values and clinical scores were observed (Figure 2). However, after correcting for multiple comparisons, none of the associations between regional MSN values and clinical severity remained statistically significant (p > 0.05, Supplementary Table S3).

In an exploratory PLS regression, spatial patterns of MSN disruption were associated with cortical gene expression gradients. Spatial distribution of PLS1-weighted gene expression maps found distinct anterior–posterior gradients, broadly mirroring the topography of MSN disruptions (Figures 3A,B). Furthermore, regional PLS1 scores were positively correlated with case–control t-values in both AIS (Spearman’s r = 0.33, p-spin < 0.001) and CIS (Spearman’s r = 0.28, p-spin < 0.001), indicating that transcriptional patterns aligned with the degree of MSN disruption across cortical regions (Figures 3C,D). These exploratory findings suggest that stroke-related MSN differences are spatially coupled with intrinsic cortical gene expression architecture. We identified differentially expressed genes based on PLS1 scores with a threshold of p-FDR < 0.005. In the AIS group, 37 PLS1⁺ genes (positive scores) were up regulated, while 215 PLS1⁻ genes (negative scores) were down regulated (Figure 3E). In the CIS group, a larger number of genes exhibited significant changes, with 148 PLS1⁺ genes upregulated and 648 PLS1⁻ genes downregulated (Figure 3E, Supplementary Table S4).

To further identify key genes involved in the pathogenesis of IS, we analyzed the intersection of differentially expressed genes (DEGs) between AIS and CIS. We identified 234 overlapping genes shared by AIS and CIS DEGs. Subsequently, we performed spatial correlation analysis between the expression levels of these overlapping genes and case–control t-values. Genes were ranked based on the Spearman’s r-value and statistical significance (FDR-adjusted p < 0.05) (Figure 3F, Supplementary Table S5). Notably, genes such as KIF5B (r = 0.30, p < 0.001), C4orf3 (Spearman’s r = 0.31, p-spin < 0.001), APMAP (Spearman’s r = −0.32, p-spin < 0.001) and STOML1 (Spearman’s r = −0.36, p-spin < 0.001) exhibited strong correlations with case–control t-values. The expression levels of the most highly ranked positively (or negatively) associated genes showed consistency with (or in contrast to) the distribution of variant regional changes in MSN strength (Figures 3G–J).

To elucidate the biological signatures across cortical MSN differences, GO and KEGG enrichment analyses were performed on DEGs identified in the AIS and CIS groups. In the AIS group, biological process (BP) terms were mainly enriched in dendritic extension, cellular amino acid metabolism, and insulin secretion regulation. Cellular component (CC) terms involved neuron projection cytoplasm and synaptic membranes, while molecular function (MF) terms were related to glycine binding and protein kinase activity (Figure 4A). KEGG analysis indicated significant enrichment in motor proteins, Parkinson disease, and neurodegeneration pathways (Figure 4B). In the CIS group, BP terms highlighted neuronal projection regeneration and microtubule-based transport. CC terms involved synaptic vesicle membranes and proteasome complexes, and MF terms were enriched in GTP binding and nucleotide binding (Figure 4C). KEGG analysis showed overlap with AIS, including motor proteins, phagosome, and neurodegenerative disease pathways (Figure 4D). Analysis of overlapping DEGs revealed common enrichment in dendritic extension, axon cytoplasm, and organic acid binding (BP, CC, MF respectively; Figure 4E). Shared pathways were mainly related to motor proteins, Parkinson disease, and multiple neurodegenerative disorders (Figure 4F). Gene-pathway network analysis identified hub genes bridging key biological functions. FZD1 and PRKCG were linked to motor protein regulation and neuronal signaling. KIF27, KIF5B, and ACTR1A were associated with axonal transport and cytoskeletal organization. TUBB6, TUBB4B, and TUBB2A were central to microtubule stability. PSMC1, NDUFC2, UCHL1, TXN, UBA52, and MAP2K7 were related to protein degradation and oxidative stress regulation. These hub genes are associated with cortical MSN differences in stroke and may reflect coordinated changes related to neuronal structure, intracellular transport, and protein homeostasis.