Ultrasound Localization Microscopy Reveals Microvascular Heterogeneity for Accurate Stratification of Glioblastoma Malignancy

1. Introduction

Glioblastoma (GBM) is a highly aggressive primary malignant brain tumor, with an incidence of 3.27 per 100,000 and a 5-year survival rate of 9.8%, reflecting its poor prognosis. Its incidence has increased by 1–2% annually over recent decades [1,2,3,4]. A key feature of GBM is microvascular proliferation, which, together with diffuse infiltration, promotes abnormal vessel formation [5,6]. The World Health Organization emphasizes the importance of microvascular and cellular proliferation in the diagnosis and classification of GBM [7]. Accurate assessment of these features is essential for predicting disease progression, guiding treatment, and enabling personalized care [4]. Treatment failure in GBM is often linked to significant intratumoral heterogeneity [8,9], which encompasses variations in cellular metabolism, the immune microenvironment, and angiogenic phenotypes, extending beyond molecular and genetic alterations. This multidimensional heterogeneity significantly impacts patient prognosis [9,10]. Notably, GBM subtypes (classical, neural, mesenchymal/telomerase reverse transcriptase-related, and proneural) exhibit marked differences in microvascular morphology and function [11,12]. Therefore, quantifying microvascular heterogeneity should extend beyond morphology to include functional assessments of vascular performance [13,14]. A comprehensive assessment is essential for accurate grading, predicting disease progression, and guiding anti-angiogenic therapy [15,16].

In clinical imaging, noninvasive in vivo acquisition of microvascular information across multiple spatial scales remains challenging [17]. Traditional methods such as CT, MRA, and fMRI can assess large-vessel anatomy and hemodynamics but lack the spatial resolution required to visualize small-caliber microvessels and capillary networks, often relying on indirect indicators such as vascular density [18,19]. Although TOF-MRA can delineate vascular structures, its sensitivity to slow-flowing blood remains limited [20]. Techniques such as PET-CT using tracers like ¹⁸F-FDG provide insights into angiogenesis and tissue metabolism but fall short in high-resolution microvascular imaging [21,22]. Microscopic techniques such as multiphoton microscopy and optical coherence tomography provide high-resolution vascular imaging but are constrained by shallow imaging depth (typically < 2 mm) and limited fields of view, thereby restricting clinical applicability [23,24]. Therefore, new methods are urgently needed to enhance image quality, temporal resolution, and data analysis, ultimately improving personalized care.

Over the past decade, ultrasound localization microscopy (ULM) has significantly advanced in vivo microvascular imaging [25]. Inspired by single-molecule localization techniques in fluorescence microscopy, ULM combines intravenous microbubble contrast agents with ultrafast ultrasound imaging to track individual microbubbles with subwavelength precision in vivo [26,27]. This approach surpasses the diffraction limit, achieving lateral and axial resolutions finer than λ/10 (≈10 μm) while preserving penetration depth and signal-to-noise ratio, thereby enabling direct visualization of capillary-level vasculature [28,29]. In data analysis, decision trees classify data or predict outcomes using if-then rules, offering efficient pattern discovery, high interpretability, and transparent visualization [30,31]. They are widely used in disease diagnosis, prognostic modeling, and identification of key factors influencing clinical outcomes [32,33]. Accordingly, we applied ULM to characterize and quantify microvascular heterogeneity, differentiate GBM malignancy strata, and compare imaging-derived features with histopathological measures, thereby evaluating their predictive performance and diagnostic value.

2. Materials and Methods

2.1. Ethics Statement

All procedures were approved by the Animal Care and Use Committee of Fudan University (Approval No. 202408008S) and conducted in accordance with ARRIVE guidelines. Each rat was treated as an experimental unit, and all procedures complied with applicable ethical regulations.

2.2. Animal Model

A total of 103 male Sprague–Dawley rats (outbred SD; stock No. D000017; 6 weeks old; 200–250 g; Jicui Pharmakon Biotechnology Co., Ltd., China) were used to establish a GBM model. Animals were housed in a specific pathogen-free facility with ad libitum access to food and water under controlled temperature (23-25 °C) and humidity conditions (55 ± 5%), with a 12-h light-dark cycle.

C6 glioma cells (National Institutes of Health, USA) were prepared at 5 × 10^5 cells/10 µL. Rats were anesthetized with ketamine and xylazine (60 mg/kg) (Figure 1A). The head was secured in a stereotaxic frame. A midline scalp incision was made, and a burr hole was drilled 3 mm to the right of the midline. A 25-gauge flat-tip Hamilton syringe (Hamilton Company, Reno, NV, USA) was advanced vertically to a depth of 5.5 mm, and 20 µL of cell suspension was injected into the caudate nucleus at 4 µL/min (Figure 1B). The needle was left in place for 2 min before withdrawal. The injection site was disinfected and sealed, and the scalp was sutured. The rats were monitored for 14 days after implantation.

2.3. Pre-Imaging Preparation

Rats were anesthetized with 3% (v/v) isoflurane, and anesthetic depth was continuously monitored. Animals were placed in the supine position. The neck was shaved and disinfected with iodine, followed by 70% (v/v) ethanol. A polyethylene catheter (PE10; inner diameter, 0.28 mm; outer diameter, 0.61 mm) was inserted into the right jugular vein (Figure 1C). Analgesia was provided via subcutaneous administration of carprofen (Cat. No. MB1412, Meilun Biotechnology Co., Ltd., Dalian, China; 5 mg/kg once daily for 5 days).

The head was secured in a stereotaxic frame. Body temperature was maintained at 36–38 °C using a heating pad (RT-0501, Kent Scientific Corporation, Torrington, CT, USA). Ophthalmic ointment was applied to prevent corneal desiccation. After skin incision, the periosteum was removed. The skull was thinned with a cranial drill (78001, RWD) until pial vessels became clearly visible. Sterile saline was intermittently applied to prevent thermal injury and minimize bleeding and swelling. Under a surgical microscope (DOM-1001, RWD), the remaining bone flap was carefully lifted and removed to preserve vascular integrity and protect the superior sagittal sinus. Special care was taken to avoid dural tearing (Figure 1D). Minor bleeding caused by superficial tissue injury was controlled using a hemostatic sponge.

2.4. Multi-Plane RF Data Acquisition

An L22-14vX LF transducer (128 elements; transmit frequency, 15.625 MHz; MS200, VisualSonics, Toronto, ON, Canada) was mounted on a stereotaxic imaging frame equipped with a three-axis motorized positioning system (VT-80 linear stage, Physik Instrumente, Auburn, MA, USA) (Figure 1E). Anesthesia was maintained with 1.5% (v/v) isoflurane. The maximum tumor cross-section was identified using a portable ultrasound system (Mindray M10; Shenzhen Mindray Bio-Medical Electronics Co., Ltd., China). Co-registered data were acquired using a Verasonics Vantage ultrasound system (Vantage 256, Kirkland, WA, USA) connected to the L22-14vX LF transducer.

SonoVue (Bracco Imaging, Massy, France) was diluted to a final volume of 5 mL (11.8 mg/mL) and infused at 80 µL/min using a microinjection pump (R462, RWD). The solution was stirred every 2 min to maintain a uniform concentration (Figure 1F). The total injected contrast agent volume was 0.4 mL. RF data were acquired for 300 s using compounded plane waves at five angles (−5°, −2.5°, 0°, 2.5°, and 5°) with a 1,000 frames/s rate. The pulse repetition frequency was 5,000 Hz, and the effective post-compounding frame rate was 1,000 frames/s. The single-pulse duration was 1/5,000 s, and the transmit voltage was 25 V [34]. Hydrophone calibration indicated a mechanical index of 0.067, within safety limits [35]. All acquisition procedures were implemented in MATLAB (The MathWorks, Natick, MA, USA; R2020a).

2.5. Image Reconstruction

Image reconstruction was performed on a desktop computer equipped with an Intel^® Core™ i7-10700 CPU (2.90 GHz) and 32 GB of random access memory. The system memory comprised KINGBANK KP426 DDR4 modules (2 × 16 GB; rated 3200 MHz, CL16, 1.35 V), operating at an effective speed of 2666 MHz (Figure 1G). Motion correction and microbubble trajectory analysis were performed using B-spline-based nonrigid registration. The reference frame for each B-mode acquisition was defined as the frame with peak contrast. Registration was performed using a control-point grid spacing of one-tenth of the image width and normalized cross-correlation as the similarity metric, with optimization performed by gradient descent (step size, 0.1; 100 iterations). During motion compensation, a displacement threshold of 1.5 times the point spread function (PSF) was applied. For most microbubbles, the PSF size was 3 × 3 pixels. Frames exceeding this threshold were excluded as outliers. The PSF was determined by analyzing the two-dimensional intensity distribution of a manually identified microbubble. A radial symmetry algorithm refined microbubble positions after initial localization. Microbubble tracks were linked using a Kalman filter (process noise covariance, 0.1 µm²; observation noise covariance, 1 µm²), followed by the Kuhn–Munkres algorithm, to estimate trajectories with a minimum length of 15 frames.

Singular value decomposition was used to separate microbubble and tissue signals in each dataset, with parameters selected empirically. Microbubble centers were localized with subpixel accuracy using a centroid-based method, and their trajectories were tracked frame by frame using the Hungarian algorithm to estimate flow velocities and directions. A temporal band-pass filter was applied to suppress both low-frequency tissue clutter and high-frequency noise. Tracking parameters were set as follows: gap = 0, maximum linking distance = 19.712 µm, and maximum detectable velocity = 197.12 mm/s. Tenfold bilinear interpolation was applied during trajectory accumulation to reconstruct vascular maps with a pixel size of 3.45 µm and an effective temporal resolution of approximately 8 min. The grid resolution was 9.856 µm. ULM reconstruction required approximately 100,000 frames and about 20 min for completion. Reconstruction quality was reduced in superficial cortical regions due to strong echoes from ear bars and acoustic shadowing caused by the skull [36].

2.6. Histopathological Analysis

After the ULM experiments, rats were euthanized, and brain tissues were harvested. Using the Mindray M10 ultrasound system, a 1-mm-thick tissue slice was obtained from the largest tumor cross-section with a brain slicing matrix (World Precision Instruments, RBMA-600S, Sarasota, FL, USA) (Figure 1H). Samples were decalcified, embedded in paraffin, and sectioned into 5-µm slices. Endothelial cells were immunostained using anti-CD31 (platelet endothelial cell adhesion molecule-1; Abcam, Cat. No. AB222783, Cambridge, UK). Biotinylated anti-α-smooth muscle actin (α-SMA; Abcam, Cat. No. AB124964) was used to label smooth muscle cells and pericytes. Whole-slide images were captured at 40× magnification and saved in a virtual slide format. Further examination at 200× magnification was performed using a MAGSCAN-NER scanner. Quantitative analysis was conducted using Tissue Studio^® (Definiens, Munich, Germany).

The α-SMA⁺/CD31⁺ ratio, defined as the vascular maturity index (VMI), was used to assess the coverage of endothelial cells by pericytes and smooth muscle cells within vessels and to represent the proportion of mature vessels among all vessels. Vessels were categorized as follows: (i) mature, characterized by an intact basement membrane and continuous pericyte and smooth muscle coverage, with an α-SMA⁺ signal intensity ≥ endothelial signal intensity (ratio ≥ 1); (ii) immature, characterized by an incomplete or absent basement membrane, exposed vessel walls lacking continuous support structures, and low pericyte and smooth muscle coverage (ratio ≤ 0.45); and (iii) transitional, representing an intermediate state between mature and immature vessels, characterized by a discontinuous basement membrane and partial, uneven coverage (ratio 0.45–1.00) [37]. The cell proliferation index (PI) was evaluated using a rabbit polyclonal anti-Ki67 antibody (Abcam, Cat. No. AB16667), which is commonly used to assess tumor growth and serves as a key indicator of cellular proliferative activity. Samples were classified as high-proliferation tumors (HPT; ≥10%) or low-proliferation tumors (LPT; <10%) [38].

2.7. Imaging-Derived Parameter Analysis

Image features were extracted using Fiji/ImageJ v1.54f (Wayne Rasband, National Institutes of Health, USA) with scale-invariant feature transform and wavelet-based methods. Supervised segmentation was performed using the Trainable Weka Segmentation tool, including preprocessing, image modulation, feature extraction, and deep learning–based vessel segmentation. The three-dimensional vascular volume was converted into two-dimensional representations, enabling extraction of vascular skeletons, centerlines, and features from binarized deep-learning model outputs [39].

Segmentation performance was evaluated using the out-of-bag error, which measures the error rate for samples not used in training the random forest model. This metric, kept below 5%, supports unbiased estimates and enhances generalizability and stability.

Orientation variance (OV) quantifies the coherence of vessel orientations and was computed as the mean of unit direction vectors across all vessel segments (Eq. 1):

$$OV=\frac{1}{N}\sum _{i=1}^N{({\theta }_i-\overline{\theta })}^2$$

(1)

where N is the number of vessels, θᵢ denotes the orientation of the i-th segment, and $\overline{\theta }$ represents the mean orientation.

The fractal dimension (FD) characterizes vascular network complexity and branching patterns. A higher FD indicates greater structural complexity. FD was estimated using the box-counting method across multiple scales, followed by linear regression of the resulting log–log relationship. The goodness of fit was required to satisfy R² > 0.99, ensuring consistency between the model and measured data (Eq. 2):

$$FD=\underset{r\to \infty }{\lim }\sum _{i=1}^N{({\theta }_i-\overline{\theta })}^2$$

(2)

where r denotes the box size, and M(r) represents the minimum number of boxes required to cover the fractal set.

To assess morphological heterogeneity, curvature was calculated for each vessel segment (Eq. 3):

$$\mathit{Curvature}=\frac{L\mathrm{c}}{L}$$

(3)

where Lc represents the actual path length of a vessel segment, and L denotes the Euclidean distance between its endpoints.

Network connectivity reflects topological heterogeneity and structural compactness and was computed using the Connectivity plugin (Eq. 4):

$$\mathit{Connectivity}=1-\frac{N_{\mathrm{c}}}{N_{\mathrm{n}}}$$

(4)

where Nc denotes the number of connected components (isolated network blocks), and Nn represents the total number of nodes, including endpoints and bifurcations.

Local thickness (LT) was used to estimate the spatial distribution of vessels, pores, and tubular diameters, thereby capturing developmental heterogeneity [40]. LT was assessed using the Local Thickness plugin (Eq. 5):

$$LT\left(x\right)=2·\underset{r}{\max }\left\{r\right|B(x,r)\subseteq S\}$$

(5)

where LT(x) denotes the local thickness at position x, B(x,r) represents a sphere centered at x with radius r, and S denotes the structure of interest.

The spatial distribution index (SDI) was used to quantify vascular spatial distribution and heterogeneity and was computed using the Simple Analysis 2D/3D plugin (Eq. 6):

$$\mathit{SDI}=\alpha ·\mathit{AI}+\beta ·\mathit{SHI}+\gamma ·K\left(\gamma \right)$$

(6)

where AI denotes the aggregation index, SHI denotes the spatial heterogeneity index, and K(γ) represents Ripley’s K function, which describes multi-scale spatial clustering. α, β, and γ serve as weighting coefficients that regulate the relative contributions of each indicator to the overall index and are determined by experimental data and model fitting [41] (Figure 1I,J).

2.8. Statistical Analysis

Statistical analyses were performed using IBM SPSS Statistics 22.0 (IBM, Armonk, NY, USA) and GraphPad Prism 9.0 (GraphPad Software Inc., La Jolla, CA, USA). Normality was assessed using the Shapiro–Wilk test for groups with sample sizes ≤ 50 and the Kolmogorov–Smirnov test for groups with sample sizes > 50. Normally distributed variables are presented as the mean and 95% confidence interval (CI). Non-normally distributed data are presented as the median (P25, P75). Spearman’s rank correlation was used to assess associations between ULM-derived metrics and histopathological measures. To improve robustness, correlation results were reported as Spearman’s r, p values, and the coefficient of determination (R²). Group differences were evaluated using the Kruskal–Wallis H test or the Mann–Whitney U test, with Bonferroni correction applied for post hoc multiple comparisons. Statistical significance was defined as p < 0.05.

For classification analysis, decision-tree models were constructed using the classification and regression tree (CART) algorithm, with Gini impurity as the splitting criterion. Cases were recursively partitioned according to the values of the imaging-derived variables to maximize between-group separation at each node. To improve model interpretability and reduce overfitting, tree growth was automatically terminated according to the internal CART stopping rules implemented in SPSS when further splitting no longer produced a sufficient reduction in impurity or when terminal-node requirements were met. The final tree was selected based on the lowest cross-validated risk estimate. Model performance was evaluated using overall classification accuracy, class-specific sensitivity, specificity, precision, recall, and F1 score, along with the estimated misclassification risk and its standard error. Feature importance was quantified based on the total reduction in Gini impurity attributable to each predictor in the final tree. In addition, 1,000 bootstrap resamples were used to assess the robustness of the classification results.

3. Results

3.1. Associations Between Microvascular Heterogeneity and Pathology

OV (r = −0.827, 95% CI [−0.733, −0.901], p < 0.001), FD (r = −0.815, 95% CI [−0.702, −0.901], p < 0.001), Cur (r = −0.883, 95% CI [−0.799, −0.937], p < 0.001), Con (r = −0.831, 95% CI [−0.728, −0.910], p < 0.001), LT (r = −0.809, 95% CI [−0.714, −0.889], p < 0.001), and SDI (r = −0.807, 95% CI [−0.706, −0.889], p < 0.001) were all negatively correlated with VMI. In contrast, OV (r = 0.889, 95% CI [0.822, 0.929], p < 0.001), FD (r = 0.867, 95% CI [0.799, 0.910], p < 0.001), Cur (r = 0.855, 95% CI [0.767, 0.911], p < 0.001), Con (r = 0.865, 95% CI [0.795, 0.908], p < 0.001), LT (r = 0.903, 95% CI [0.842, 0.941], p < 0.001), and SDI (r = 0.892, 95% CI [0.829, 0.934], p < 0.001) were all positively correlated with PI (Figure S1). Linear regression analysis further showed that the FD-based model explained the greatest proportion of variance in VMI (R² = 0.547). For PI, the LT-based model explained the greatest proportion of variance (R² = 0.915), followed by the SDI-based model (R² = 0.827) (Figure 2).

3.2. Microvascular Heterogeneity Across Different VMI Levels

Compared with immature microvessels, mature and transitional microvessels showed significantly lower values for all parameters (p ≤ 0.007). Immature microvessels exhibited disorganized trajectories and a highly complex branching architecture, characterized by excessive branching, pronounced tortuosity, and multiple abnormal anastomoses, thereby forming shunt-like perfusion pathways. The lumina exhibited irregular dilation, the spatial distribution was patchy and clustered, and the basement membrane was largely absent. Compared with mature microvessels, transitional microvessels showed significantly greater tortuosity (p = 0.002), whereas the other parameters showed nonsignificant upward trends (p = 0.074–0.529) (Table S1; Figure 3).

3.3. Microvascular Heterogeneity Across Different PI Levels

All microvascular heterogeneity parameters were significantly higher in the high-PI group than in the low-PI group (p < 0.001) (Table S2; Figure 4). Morphologically, high PI was associated with disoriented and irregular vessel alignment, greater branching complexity, more tortuous vascular courses, and increased numbers of bifurcation points and abnormal perfusion shunts. Numerous large-caliber, developmentally immature vessels were observed, along with an uneven and distinctly patchy spatial distribution.

3.4. Classification of Microvascular Maturity

A decision-tree model based on microvascular heterogeneity features performed well in stratifying microvascular maturity, achieving an overall accuracy of 90.29% (93/103, p = 0.029). For mature microvessels, the model achieved high accuracy (92.16%), sensitivity (94.00%), and recall (94.00%), indicating robust identification of the mature phenotype. For transitional microvessels, the model achieved high specificity (95.95%), indicating good discriminative ability; however, sensitivity was lower (86.21%), suggesting that a subset of cases was misclassified. For immature microvessels, performance was strong, with an accuracy of 96.12%, sensitivity of 95.45%, and specificity of 96.30%, indicating effective recognition of the immature phenotype (Table 1). Notably, FD served as the root node of the decision tree, suggesting that vascular network geometric complexity was the dominant determinant of microvascular maturity (Figure 5). Cross-validation indicated a low estimated misclassification risk (stratified risk = 0.097) and revealed no evidence of overfitting at the current sample size.

3.5. Classification of Cellular Proliferative Status

Although OV differed significantly between proliferative states, it showed substantial multicollinearity with the other variables (Table S3). In contrast, the remaining metrics showed more linear and stable relationships with PI. Therefore, after excluding OV, the optimized model achieved a higher split gain and an improved overall accuracy of 92.23% for proliferative-status stratification (95/103, p = 0.026). These results indicate that the decision-tree model effectively distinguishes between proliferative states. Recall (93.65%) and precision (93.65%) were identical, further supporting the model’s stability and reliability (Table 1). SDI was the primary feature at the initial split and served as the key imaging-derived indicator for differentiating high- and low-proliferative states, followed by LT (Figure 6). Cross-validation showed a low estimated misclassification risk (stratified risk = 0.078) and revealed no evidence of overfitting at the current sample size.

4. Discussion

Biopsy-based histopathological features remain essential for personalized therapy; however, their utility is often limited by intratumoral heterogeneity. In GBM, pathological angiogenesis and microvascular heterogeneity are key indicators of disease progression. Park et al. [42] proposed a vascular heterogeneity model highlighting its role as a marker of vascular maturity and perfusion. Our findings demonstrate that ULM can effectively monitor and quantify microvascular heterogeneity in a rat model while providing sufficient imaging depth at a relatively low cost. This approach may facilitate GBM malignancy stratification and support the clinical translation of ULM-derived biomarkers.

Our study demonstrates that noninvasive ULM can effectively detect spatial variations in microvessel distribution. We observed an inverse association between microvascular heterogeneity and vascular maturity, with mature microvessels showing greater stability and a more uniform distribution. Lower vascular maturity was associated with greater heterogeneity, likely reflecting immature microvessels within dysfunctional stroma that generate an unstable and frequently remodeled environment [43,44]. This heterogeneity may promote immune evasion and accelerate tumor growth [45]. FD, which reflects vascular network complexity and remodeling, captures features of vascular stability and maturity that are central to our analysis. Rapid tumor growth increases oxygen and nutrient demands, thereby triggering disordered angiogenesis and stromal remodeling, which in turn increase geometric complexity and spatial heterogeneity. Consequently, microvascular heterogeneity increases with proliferative activity [46,47,48]. LT and SDI measure vascular size and distribution, respectively, correlating with perfusion and proliferation, thereby explaining a substantial proportion of the variance.

Previous studies have investigated microvascular development and tumor microenvironment heterogeneity in GBM, emphasizing that vascular dysfunction induces hypoxia and acidosis, thereby reducing therapeutic efficacy [49,50]. Using multiparametric MRI and histology, these studies have shown how vascular changes influence the tumor microenvironment [8,51,52]. We propose three functional archetypes: mature microvessels exhibiting regular structure and low heterogeneity, immature microvessels demonstrating disordered structure and instability, and transitional microvessels with intermediate features, such as discontinuous basement membranes, uneven pericyte coverage, and reduced endothelial support, which may increase susceptibility to bending and distortion [44]. This tortuosity suggests a partially stable yet remodeling-prone state. The lack of changes in other heterogeneity metrics suggests these vessels may transition toward either maturity or immaturity, depending on environmental factors. Vascular normalization strategies combined with immunotherapy aim to stabilize vessel walls and improve tumor perfusion and drug delivery. High proliferative activity may induce vascular heterogeneity through the hypoxia–vascular endothelial growth factor pathway [44,45]. This process begins with increased metabolic demand, which exacerbates hypoxia and activates vascular endothelial growth factor signaling, leading to the rapid formation of abnormal neovessels. These vessels disrupt blood flow, thereby worsening hypoxia and reinforcing a feedback loop between hypoxia and angiogenesis, which results in abnormal vascular distribution in highly proliferative areas [50,51,52]. GBM-induced vascular heterogeneity, driven by high proliferation, metabolic alterations, and hypoxia-induced angiogenesis, represents a critical yet underexplored therapeutic target [53,54,55].

At the capillary level, ULM provides super-resolution information on vascular structure and flow, revealing subtle changes during the transition between mature and immature vasculature [26]. Decision trees, which capture nonlinear relationships, can prioritize variables in the presence of multicollinearity. In immature and highly proliferative states, most metrics tend to increase globally, allowing decision trees to identify these phenotypes with relatively few splits [33]. Compared with more complex black-box algorithms, the CART framework provides explicit hierarchical decision rules, thereby facilitating biological interpretation and potential clinical application. In our dataset, low cross-validation risk estimates for vascular maturity and proliferative-state classification support the robustness of the final model. Furthermore, removing OV from the proliferative-state model improved classification efficiency, suggesting that even in tree-based methods relatively tolerant to collinearity, reducing redundant information can enhance model stability. These characteristics make the proposed workflow attractive for future multimodal integration with molecular markers and radiomic features. This approach may improve predictive accuracy for pathological stratification by leveraging decision trees to simplify complex patterns into well-defined thresholds, thereby supporting preoperative assessment, radiotherapy planning, and follow-up [45,55]. Mapping ULM-derived heterogeneity features to tumor maturity or proliferation provides a more comprehensive perspective beyond localized biopsy samples, potentially guiding therapy, dose adjustment, and early detection of high-risk areas to improve treatment precision [29,50]. Integrating these features with decision-tree models yields an interpretable imaging workflow. However, further research on data standardization, validation, and clinical application is still needed.

Several limitations should be acknowledged. First, although transcranial ULM enables microcirculatory imaging, its clinical application requires lower ultrasound frequencies: 15 MHz for rodents, 6 MHz for neonatal fontanelle imaging, and 2 MHz for adult transcranial imaging [56]. Second, only male rats were used to minimize the effects of hormonal fluctuations; however, future studies should include both sexes to improve robustness [57]. Third, standard histological workflows typically use three slides per sample, which may introduce bias because of uneven vascular distribution. Whole-slide scanning may provide a more comprehensive and less biased assessment of vascular distribution [58]. Additionally, studies in human subjects are needed to validate the feasibility and effectiveness of this technique in clinical settings. Advances in imaging and computational modeling are facilitating the integration of multimodal, disease-specific data for personalized treatment, with tumor subtype, molecular alterations, and angiogenic features becoming key components of precision medicine.

5. Conclusions

Using ULM to assess microvascular heterogeneity, we identified distinct patterns across vascular maturity levels and proliferative states in GBM. These findings underscore the value of microvascular heterogeneity in malignancy stratification and suggest that it may serve as a feasible adjunct to histopathological diagnosis in GBM. We anticipate that these findings will facilitate clinical translation and support earlier detection of GBM.