Introduction

The critical role of cerebrovascular structure (CvS) in healthy evolution and maintenance of brain function has been increasingly highlighted by recent findings. These studies have underscored the importance of CvS in not only supplying oxygen and nutrients to the brain1 but also in removing metabolic byproducts2. Abnormalities in CvS, such as those leading to aneurysms3, stenosis4, and strokes5, have been recognized as significant indicators of cerebrovascular diseases with substantial implications for morbidity and mortality. Furthermore, the connection between vascular health and neurodegenerative diseases has been strengthened by new researches6. For instance, pericytes, cells essential for blood–brain barrier function and cerebral blood flow, have been identified as crucial in linking cerebrovascular health to neurodegenerative diseases, including Alzheimer’s disease (AD) and other forms of vascular dementia7,8,9. The increasing prevalence of cerebrovascular diseases has heightened societal concern and recognition of the imperative to understand CvS in both health and disease across different imaging scales. This understanding is crucial for hierarchically identifying life-threatening vascular structural abnormalities, as well as vascular factors causing severe cognitive impairment and dementia10,11,12.

Analogous to growth charts13, normative modeling framework is a rising alternative to mapping individual differences14. To evaluate brain growth with respect to aging at population-14,15 and individual-levels16, models were built for different hierarchies of spatial scales15,17,18. It is widely acknowledged that normal brain evolution is supported by healthy CvS. A demand of unraveling the complex CvS by growth charts or normative modeling arises, which remains largely unexplored, not even for the span of adulthood. This can be due to the lack of tailored pipeline for large-scale cerebrovascular data analysis and the difficulties in processing multi-center cerebrovascular data typically with inconsistent contrast and signal-to-noise ratio. To fill this crucial gap, in this work, we base our investigation of CvS on a large dataset consisting of 2841 paired T1-weighted and time-of-flight magnetic resonance angiography (TOF-MRA) volumes from multiple research centers. To realize all-round inspection, efficient reconstruction and analysis of CvS is a fundamental step19,20,21. However, accurately building and analyzing of 3D cerebral arteries is challenging. Yet the complete reconstruction of cerebral arterial structures from multi-center TOF-MRA volumes has not been achieved. Traditional cerebrovascular segmentation methods generally exploit shallow intensity features, such as intensity similarity among vascular regions22,23,24 and intensity difference between vessels and background25,26,27. These methods lack robustness against intensity inhomogeneity in vascular regions and interference from background tissues. Hence, they are incapable of segmenting arterial structures without interruption of vessel branches or confusion from noisy background. Past decades have witnessed the prosperity of deep-learning-based biomedical image segmentation approaches owing to their automatic feature engineering, powerful self-learning and generalization capabilities28,29,30. But a barrier remains, that is the degraded performance of deep-learning-based models when applied to multi-center TOF-MRA volumes acquired with different scanner vendors and centers in clinical practice31. Current available datasets for cerebrovascular segmentation often include only a few subjects or are scanned at a single clinical center using the same imaging protocol32,33, which compromised the performance of deep-learning methods for multi-center cerebrovascular segmentation. Therefore, a large-scale multi-center dataset, as well as a robust system, for automatic and accurate segmentation of cerebral arteries are required.

Directly addressing the above issues, we first present the multi-center cerebrovascular dataset consisting of 2841 MRA volumes. Then a deep-learning model CereVessSeg incorporating self-supervised learning (SSL) strategy named CereVessPro to take full advantage of unlabeled data is proposed, segmentation performance of which outperforms current state-of-the-art models28,32,34,35,36,37,38,39,40 by 0.5%–4.4% in MRA volumes from healthy subjects and 1.7%–10.22% in MRA volumes from subjects with brain pathology in terms of DICE similarity coefficient (DSC)41. Besides, an automatic deep-learning model-centric pipeline for hierarchical quantitative analysis of CvS is proposed to extract vascular (and also cortical) features hierarchically from (1) whole brain, (2) four typical vascular regions including anterior cortical area (ACA), posterior cortical area (PCA), middle cortical area (MCA), and circle of Willis (CoW), and (3) conventional Brodmann areas after registration of each MRA volume to atlas defined in Montreal Neurological Institute (MNI) space. Utilizing this pipeline, we investigate hierarchically gender-stratified normative trajectories, including CVs versus aging and AVs versus aging for both CN and UK subjects. The evolutionary trends of CVs, showing a similar near-linear decreasing trend to previous findings15, were highly consistent across CN and UK subjects in each spatial hierarchy. However, AVs showed different trends across different hierarchies. To be specific, starting from the age of 20, normative trajectories of AVs for the whole brain and three typical brain regions, including ACA, MCA, and PCA, exhibited near-linear decreasing trends, in contrast to a rising trend of AVs in CoW. In general, AVs in most Brodmann areas showed near-linear decreasing trends. But parabolic or rising trends were also observed, suggesting desynchronized evolutionary patterns among different Brodmann areas. Furthermore, while CVs exhibited a clear gender gap, AVs displayed less obvious gender-related differences, thus distinguishing evolutionary patterns of regional AV from regional CV. Further investigations contrasting CN healthy subjects with those diagnosed with large territorial stroke (LTS), lacunar stroke (LS), AD revealed different levels in reduction of CVs and AVs wherein AD subjects suffer the most. These findings suggest that CVs and AVs could serve as complementary biomarkers in understanding the brain pathological consequences.

Results

Framework

The proposed pipeline for hierarchical quantitative analysis of cerebrovasculature is illustrated in Fig. 1. Paired T1 and TOF-MRA images were acquired from each subject. TOF-MRA volumes were utilized for cerebrovascular reconstruction using CereVessSeg, a segmentation model refined from nnUNet28 by incorporating attention module and loss function emphasizing learning from hard examples. To further improve the model performance, a contrastive learning strategy, CereVessPro, was adopted during model pretraining. TOF-MRA volumes were inputted into this model for individualized cerebrovascular reconstruction. T1 volumes were used to generate an average group-level structural template in MNI space and individualized deformation field in subject space. The deformation fields characterize the biprojection between MNI and individual spaces. CVs and AVs were quantified and contrasted across three spatial hierarchies, including a large spatial hierarchy defined by the whole brain, a medium spatial hierarchy defined by an arterial atlas with four typical areas, and a small spatial hierarchy defined by the conventional Brodmann atlas consisting of 82 Brodmann area. Normative models were further constructed, using locally estimated scatterplot smoothing15,42, for investigating arterial and cortical evolutionary trajectories across these spatial hierarchies.

Fig. 1: The proposed pipeline for hierarchical quantitative analysis of cerebrovasculature.
figure 1

To start with, subjects were admitted into different hospitals to create records. Basic demographics were measured and recorded, followed by guiding the subjects to different MR scanners within each hospital for image acquisition. Typically, paired T1 and MRA volumes were acquired. Each T1 image was used to generate a deformation field for individualized biprojection between the subject and MNI space. Each subject’s cerebrovascular segmentation was obtained by applying the deep-learning segmentation model to the corresponding TOF-MRA volume. For feature quantification, the atlas defined in MNI space was transformed to subject space using the previously created deformation field. Normative models were then built using group-wise statistical vascular and cortical features.

Full size image

Dataset

We constructed a large-scale, multi-center dataset comprising 2841 pairs of T1 and TOF-MRA volumes, totaling 475,754 MRA slices and 322,944 T1-weigted slices. The dataset includes 1126 pairs acquired from healthy UK subjects (termed as UK-Hset), 950 pairs from healthy Chinese subjects (termed as CN-Hset), 642 pairs from Chinese subjects with brain pathology (termed as CN-Pset), and 123 additional pairs from Chinese subjects with Alzheimer’s Disease (termed as CN-ADset). Demographic information for each dataset is illustrated in Table 1. For UK-Hset, data of 570 subjects were collected from IXI dataset (http://brain-development.org/ixi-dataset/), wherein MRI scans were acquired at Hammersmith Hospital (IXI-HH), Guy’s Hospital (IXI-Guys), and the Institute of Psychiatry (IXI-IOP) in London, UK. The remaining data of 556 UK healthy subjects were sourced from OASIS-3 dataset43. For CN-Hset, data were obtained from the Affiliated Hospital of Guizhou Medical University (AHGMU) in Guizhou Province, China. For CN-Pset, data were collected from Datian General Hospital (DGH) in Fujian Province, China. For CN-ADset, data were collected from Zhongxiang Hospital of Traditional Chinese Medicine (ZHTCM) in Hubei Province, China. Based on the diagnosis information (see Supplementary Table 1), the subjects in CN-Pset can be subdivided into 165 subjects with large territorial stroke (termed as CN-LTSset) and 334 with lacunar stroke (termed as CN-LSset).

Table 1 Demographic information for UK-Hset, CN-Hset, CN-Pset, CN-ADset, respectively
Full size table

Moreover, we manually annotated 271 (healthy:pathological = 150:121) TOF-MRA volumes, including 28,128 slices. This annotated dataset, named CereVessMRA, supports the development of deep-learning methods for cerebrovascular segmentation. The 150 annotated healthy TOF-MRA volumes from UK subjects were sourced from three different institutions, as part of the UK-Hset. The 121 annotated pathological TOF-MRA volumes from Chinese subjects were provided by DGH in Fujian Province, China. This large-scale and multi-center cerebral arterial annotation dataset facilitated the development of generalizable deep-learning models for automatic and precise cerebral arterial segmentation.

Evolutionary patterns of AVs and CVs were consistent between UK and CN healthy subjects across the whole brain and the four typical brain regions

Both groups showed similar general trends in CVs and AVs evolution, as observed in the left panel of Fig. 2A, B (please refer to Supplementary Figs. S1, S2 for age and gender distributions of UK-Hset and CN-Hset. No statistical differences in gender or age were identified between these two groups, as shown in Supplementary Tables S1, S2). To be specific, at the whole-brain level, decreasing trends of CVs and AVs were observed. This declining evolutionary trend was preserved in other three typical brain regions, including ACA, MCA, and PCA. In contrast, while both groups showed a similar rising trend of AVs in CoW region, CVs in CoW region exhibited a decreasing trend. Comparatively, CN subjects displayed a slower trend than their UK counterparts in the evolution of both AVs and CVs. Additionally, larger variations in the evolutionary curves were observed in CN subjects compared to UK subjects for both AVs and CVs.

Fig. 2: Consistent evolutionary patterns of AVs and CVs for UK and Chinese healthy subjects.
figure 2

Data are presented as smooth curves obtained by locally estimated scatterplot smoothing (LOESS) together with 95% confidence interval (CI) curves. A Evolutionary patterns of AVs across the whole brain and four typical brain regions for CN (top panel) and UK (bottom panel) subjects. B Evolutionary patterns of CVs across the whole brain and four typical brain regions for CN (top panel) and UK (bottom panel) subjects. AV arterial volume, CV cortical volume, w/i within, ACA anterior cerebral area, MCA middle cerebral area, PCA posterior cerebral area, CoW circle of Willis. Source data are provided as a Source Data file.

Full size image

Within each group, statistical significance existed in the distributions of AVs and CVs across genders, as shown in Supplementary Tables S3–S6. To be specific, for UK healthy subjects, AVs distributions showed significant statistical differences across genders in spatial hierarchies including the whole brain (pcorrected < 0.01), CoW (pcorrected < 0.001), ACA (pcorrected < 0.01), and MCA (pcorrected < 0.05) across genders, with PCA (puncorrected = 0.0796) as an exception. For CN healthy subjects, significant statistical differences in AVs distributions were found in all these spatial hierarchies (pcorrected < 0.001, 0.05, 0.05, 0.001 for the whole brain, CoW, ACA, MCA, and PCA, respectively) across genders. For CVs distributions, gender differences were clearly recognized. Both groups exhibited the same level of statistical differences across genders (pcorrected < 0.001 for each spatial hierarchy).

The degree of variation in CVs and AVs differs among large territorial stroke, lacunar stroke, and AD

As displayed in Fig. 3A and Supplementary Fig. S3A, subjects with AD were accompanied by significant reduction of AVs and CVs (refer to Supplementary Fig. S4 for the demographic information of CN-ADset, Supplementary Figs. S5, S6 for histograms of age-matched subjects from CN-Hset and CN-ADset for further statistical analysis using two-sample t-test). In average, AD males were found to have AVs that were 23.5%, 36.4%, 12.0%, 19.9%, 30.2% smaller than those of their healthy control (HC) counterparts at the whole-brain level, CoW, ACA, MCA, and PCA, respectively, as detailed in Supplementary Table S7. For CVs, the reductions were 24.7%, 48.1%, 24.4%, 22.7%, 13.9%, as shown in Supplementary Table S8. For females, AVs were reduced by 19.8%, 40.0%, 17.2%, 16.4%, and 26.3% as shown in Supplementary Table S9, while CVs were reduce by 20.4%, 40.7%, 20.0%, 18.0%, 10.3% as shown in Supplementary Table S10.

Fig. 3: Comparisons of CN healthy subjects and subjects with brain pathology in terms of AVs.
figure 3

Data are presented as smooth curves obtained by locally estimated scatterplot smoothing (LOESS) together with 95% confidence interval (CI) curves. A Comparison of CN healthy subjects and AD subjects in terms of AVs. B Comparison of CN healthy subjects and subjects with large territorial stroke in terms of AVs. C Comparison of CN healthy subjects and subjects with lacunar stroke in terms of AVs. Corresponding comparisons in CVs are provided in Supplementary Fig. S3. AV arterial volume, w/i within, ACA anterior cerebral area, MCA middle cerebral area, PCA posterior cerebral area, CoW circle of Willis. Source data are provided as a Source Data file.

Full size image

Subjects with LTS also exhibited reduction in CVs and AVs as shown in Fig. 3B and Supplementary Fig. S3B, but at a relatively milder extent compared to those observed in AD subjects (refer to Supplementary Fig. S7 for the demographic information of CN-LTSset, Supplementary Figs. S8, S9 for histograms of age-matched subjects from CN-Hset and CN-LTSset for further statistical analysis using two-sample t-test). Significant differences in CVs and AVs were not consistently observed across hierarchies, as demonstrated in Supplementary Tables S11–S14. For both genders, AVs and CVs provided complementary information. To be specific, for males, significant differences were found in AVs between HC and LTS subjects within ACA, but not in CVs. For females with LTS, while no significant differences existed in CVs within ACA and PCA, AVs showed significant reductions, being 25.5% and 16.1% smaller in ACA and PCA, respectively, compared to HC.

Reduction of CVs and AVs were also observed in subjects with LS as displayed in Fig. 3C and Supplementary Fig. S3C (see Supplementary Fig. S10 for the demographic information of CN-LSset, Supplementary Figs. S11, S12 for histograms of age-matched subjects from CN-Hset and CN-LSset for further statistical analysis using two-sample t-test). As demonstrated in Supplementary Tables S15–S18, LS deteriorated AVs in both genders, but males experienced greater reductions in CVs compared to females, particularly within ACA and PCA.

Similar region-specific evolutionary trends in AVs were observed in healthy subjects from both CN and the UK, while some regions exhibited distinct patterns

While decreasing trends of AVs and CVs at large spatial scales, such as the whole brain and four typical brain regions were evident, these trends did not consistently apply to more localized regions, especially when examined at a higher resolution spatial hierarchy. As displayed in Fig. 4A, B, at the small spatial scale wherein Brodmann areas are defined, generally decreasing trends of AVs similar to those observed at the whole-brain level were identified in regions including PMC, SMC, APC, OA, PVC, SVC, AVC, MTG, STG, PiriC, PCC, DEC, FG, TA, POper, PT, DPC2, POrb, and RA (see Supplementary Table S19 for the full names of these abbreviations) for both CN and UK healthy subjects (refer to Supplementary Fig. S13 for normative models of AVs in each Brodmann area). However, this decreasing trend was not consistent across all Brodmann areas. Different evolutionary patterns of AVs, including parabolic or even increasing trends, were observed in other Brodmann areas including PSC, ITG, VACC, SA, VEC, RCC, DACC, EA, TA, AG, SG, AC1, and AC2. In contrast to CVs, where evolutionary trends generally showed consistent decreasing patterns across different Brodmann areas (refer to Supplementary Fig. S14 for normative models of CVs in each Brodmann area), AVs exhibited region-specific evolutionary trends. Evolutionary trends of AVs were not uniformly preserved across spatial hierarchies. Supplementary Videos 1–16 provide dynamic visualizations of the evolutionary trends of CVs and AVs. While CVs exhibited mostly consistent decreasing trends across Brodmann areas, AVs in some parietal lobe regions were found to peak during middle or even old ages.

Fig. 4: Mappings of regional arterial and cortical evolutionary trends for both genders at the small spatial scale defined by Brodmann areas.
figure 4

The horizontal axis represents increasing age. Value from each Brodmann area is color-mapped by the deviation from the regional average across the age range. A Lateral and medial views of arterial and cortical evolutionary trends for CN and UK healthy females. B Lateral and medial views of arterial and cortical evolutionary trends for CN and UK healthy males. Source data are provided as a Source Data file.

Full size image

Segmentation model equipped with channel-wise attention, hard example mining strategy and contrastive learning of voxel propagation

A supervised-learning model CereVessSeg with a contrastive learning method CereVessPro was developed for cerebral vascular segmentation. With nnUNet28 as the backbone, the model CereVessSeg additionally introduced attention mechanism-based module (AM) to enhance discriminative capacity of feature presentation for cerebral vessels. Further, a hard example mining strategy was incorporated in the conventional cross-entropy loss, named LmCE, to mitigate segmentation errors caused by complex background tissues as well as low-contrast vascular branches. The contrastive learning CereVessPro was proposed to pretrain CereVessSeg model, which exploits position consistency of cerebral vessels among different subjects, thus providing a robust initialization for CereVessSeg in the cerebrovascular segmentation task. Subsequently, CereVessSeg was finetuned using our annotated CereVessMRA dataset. We compared our method with several convolution-based medical segmentation methods, including 3D U-Net, V-Net, DenseVoxNet, VoxResNet, U-ception, DeepvesselNet, and nnU-Net28,32,34,35,36,37,38, and also transformer-based methods, including UNETR and SwinUNETR39,40. All these compared methods were trained from scratch using our CereVessMRA dataset.

Four metrics, including intersection over union (IoU), DSC, average symmetric surface distance (ASD) and 95% Hausdorff distance (HD95)41, were used to evaluate the segmentation performances of different methods as well as for quantitative comparisons among different methods. The IoU and DSC, defined as follows, quantify the degree of similarities between ground truth and model predictions:

$$\,{\mbox{IoU}}\,({{\bf{y}}},\hat{{{\bf{y}}}})=\frac{| {{\bf{y}}}\cap \hat{{{\bf{y}}}}| }{| {{\bf{y}}}\cup \hat{{{\bf{y}}}}| },$$
(1)
$$\,{\mbox{DSC}}\,({{\bf{y}}},\hat{{{\bf{y}}}})=\frac{2| {{\bf{y}}}\cap \hat{{{\bf{y}}}}| }{| {{\bf{y}}}|+| \hat{{{\bf{y}}}}| },$$
(2)

where y and \(\hat{{{\bf{y}}}}\) represent ground truth and the predicted binary segmentation results, respectively. ∣.∣ represents the number of voxels. ASD measures the average distance between the closest points on the surface of ground truth ys and the surface of predicted binary segmentation \({\hat{{{\bf{y}}}}}_{s}\), formulated as following:

$$\,{{\mbox{ASD}}}\,({{{\bf{y}}}}_{s},{\hat{{{\bf{y}}}}}_{s})=\frac{1}{| {{{\bf{y}}}}_{s}|+| {\hat{{{\bf{y}}}}}_{s}| }\left(\sum_{p\in {{{\bf{y}}}}_{s}}\mathop{\min }_{q\in {\hat{{{\bf{y}}}}}_{s}}d(p,q)+\sum_{q\in {\hat{{{\bf{y}}}}}_{s}}\mathop{\min }_{p\in {{{\bf{y}}}}_{s}}d(q,p)\right),$$
(3)

HD95 is similar to conventional Hausdorff distance, which describes the maximum minimal distance between points on ys and \({\hat{\bf{y}}}_{s}\) as follows:

$${\rm{HD}} \, ({\bf{y}}_{s},{\hat{\bf{y}}}_{s})= \max \left(\mathop{\max}\limits_{p \in {\bf{y}}_{s}} \mathop{\min}\limits_{q \in {\hat{\bf{y}}}_{s}}d(p , q) , \mathop{\max}\limits_q \in {\hat{{\bf{y}}}_{s}} \min_{p \in {\bf{y}}_{s}}d(q,p)\right),$$
(4)

Notably, HD95 eliminates the impact of a small subset of outliers by using the 95th percentile of distances.

The performances of different methods on TOF-MRA volumes from healthy subjects and brain pathological subjects are displayed separately in Table 2. Mean squared distance (MSD) accesses the boundary consistency on average, with smaller MSD values suggesting fewer boundary gaps between ground truth and predictions. HD95 measures the network’s capability to suppress the influence of non-vascular regions during segmentation, while allowing for some outliers. Thus, a combination of a low DSC and a high MSD suggests compromised segmentation performance due to low capability in accurately recognizing vascular boundaries. And a combination of a low DSC and a high HD95 implies that deteriorated segmentation performance due to limited capability in filtering out some non-vascular regions. For MRA volumes from healthy subjects, while nnUNet28 was able to outperform other existing popular methods including DenseVoxNet, VoxResNet, 3D U-Net, DeepVesselNet, U-ception, V-Net, UNETR, and SwinUNETR, the proposed method can further advance its capability in boundary detection and outlier suppression. Interestingly, high accuracy in MRA volumes from healthy subjects did not guarantee equally good performance on MRA volumes from subjects with brain pathology. Besides our methods, methods including DenseVoxNet, 3D U-Net, DeepVesselNet, V-Net, and nnUNet showed higher accuracies for MRA volumes from pathological subjects compared to MRA volumes from healthy subjects, whereas VoxResNet, U-Ception, UNERT, and SwinUNETR exhibited decreased performance. Brain pathology always includes a decrement in the size of brain cells (please refer to Supplementary Fig. S15 for an illustration of brain pathology). Pathology can be either generalized, which means that all of the brain has shrunk or it can be focal, affecting only a limited area of the brain. Thus, changed or collapsed cerebral vasculature, especially in distal part of vessels, may be observed in pathology brains. In return, this morphological change of vasculature may lower the difficulties of segmentation due to losses of distal small vessels, which may explain the reason why improved performance were observed as mentioned above. Overall, among all other methods, our proposed method also achieves best segmentation performance for MRA volumes from pathological subjects.

Table 2 Quantitative evaluation of different segmentation methods
Full size table

Supplementary Tables 20, 21 display the statistical significance results from t-tests comparing our proposed CereVessSeg model, with and without CereVessPro pretraining, against other methods on MRA volumes from healthy and pathological subjects. CereVessSeg, both with and without pretraining, performed significantly better than all compared methods for MRA volumes from subjects with brain pathology. But for MRA volumes from healthy subjects, CereVessSeg without pretraining was not significantly better than SwinUNETR and nnUNet in terms of IoU and MSD. With pretraining, CereVessSeg became significantly better than all other compared methods, indicating the boost by CereVessPro in generally improving model performance for multi-center data.

Visualization comparisons are shown in Fig. 5 (refer to Supplementary Fig. S16 for more visualized comparisons). Figure 5 reveals that DenseVoxNet, VoxResNet, 3D U-Net, DeepVesselNet, U-ception, UNETR, SwinUNETR, and V-Net exhibited many false positives (blue) and false negatives (green) in their segmentation results, compared to ground truth, particularly in the areas highlighted by arrows. In comparison, nnUNet produced fewer segmentation errors. Moreover, CereVessSeg, enhanced with the proposed AM and LmCE, further reduced both false positives and false negatives. Finally, with CereVessPro pretraining, CereVessSeg achieved the most complete and accurate cerebrovascular segmentation, closest to the ground truth compared to all other methods.

Fig. 5: Visualized segmentation results from different methods.
figure 5

For every method, the segmentation result for a whole TOF-MRA volume is presented on the middle raw. Then two regions marked by the yellow and gray boxes are zoomed in and then displayed above and below to show details of segmentation. In every sub-figure, true positives are marked in red; false positives are marked in blue; false negatives are marked in green. White arrows point to areas with significant segmentation differences.

Full size image

We conducted ablation studies to measure the effectiveness of the proposed AM and LmCE in boosting the performance of CereVessSeg. The results are exhibited in Table 2. As seen from Table 2, each of the two proposed modifications independently improved segmentation performance, surpassing other state-of-the-art models.

Additionally, we ablated different SSL methods for CereVessSeg, including masked volume inpainting44, BYOL45, and SimSiam46. The segmentation performances of CereVessSeg using these methods for pretraining are listed in Table 2. As shown, inpainting, BYOL or SimSiam did not enhance CereVessSeg’s performance for MRA volumes from healthy or pathological subjects. In comparison, the proposed CereVessPro significantly improved performance on both types of MRA volumes, indicating its superiority for our segmentation task.

Supplementary Fig. S17 shows some failed cases by our proposed CereVessSeg, pretrained with CereVessPro. One kind of failed cases are vascular branches with very low intensity contrast against the background. Despite incorporating the proposed AM and LmCE, our methods still struggle to detect a few dim peripheral vessels and vascular ends, leading to under-segmentation. Another kind of failed cases are vessels in the head scalp and skull, which resemble cerebral vessels in appearance but are not part of the cerebrovascular system. Due to the patch-wise segmentation strategy used for whole-brain MRA inference, the segmentation model has difficulty perceiving the spatial locations of each patch in relation to the entire brain, making it challenging to tell apart cerebral vessels from vessels in the head scalp or skull. However, it is important to note that the first kind of failed cases are rare and involve only small regions, exerting minimal effects on cerebrovasculature quantification. In the second kind of cases, the over-segmented vessels in the head scalp and skull are excluded from subsequent analysis of cerebral vascular structures. Thus, these failed cases would not affect the overall quantification of cerebrovasculature.

External validation

We conducted an experiment to validate our proposed CereVessSeg with CereVessPro on an external dataset, TubeTK (https://public.kitware.com/Wiki/TubeTK/Data). The TubeTK dataset contains MRA volumes from healthy subjects, acquired using a Siemens Allegra head-only 3T MR system. Voxel spacing of these MRA volumes is 0.5 × 0.5 × 0.8 mm3. Supplementary Fig. S18 presents segmentation examples of our segmentation method on TubeTK dataset. Compared to the provided annotations, which miss many vascular branches, our proposed CereVessSeg produced more complete cerebrovasculature segmentation. Besides, upon the training process not using MRA volumes from AD subjects, CN-ADset was utilized for external validation as well. Supplementary Fig. S19 presents segmentation examples on CN-ADset, showing our segmentation method accurately segmented almost every vascular region and thereby achieved highly complete and precise cerebrovascular segmentation results. The segmentation performances on both the TubeTK and CN-ADset demonstrated the generalizability of our CereVessSeg with CereVessPro, thanks to multi-center training data and the efficient contrastive learning method.

Discussion

In this work, we aimed to map the hierarchical evolutionary patterns of cerebrovasculature with respect to aging, based on a quantitative investigation of 2841 MRA volumes. A multi-center cerebrovascular dataset consisting of 271 MRA volumes with accurate manual annotation of cerebral arteries was built. To the best of our knowledge, this is the largest publicly available dataset. Utilizing this dataset, a deep-learning-based model for cerebrovascular segmentation that surpasses current state-of-the-art methods was proposed. This segmentation model enabled accurate segmentation of cerebrovasculature when applied to the multi-center dataset consisting of 2841 MRA volumes. Hierarchical vascular features were extracted from regions representing different spatial scales after individualized registration of masks from MNI space to subject space. Normative models for cerebrovasculature were built to reveal hierarchical evolutionary patterns. Previous studies have revealed the global and regional evolutionary trajectories of cerebral cortex across the whole lifespan, from 17 weeks post-conception to nearly 100 years of age15,17. However, as the only source of blood supply to neocortex47, the cerebrovasculature plays a central role in healthy human brain function48,49,50 and cerebrovascular diseases51,52. A healthy normative model of cerebral arterial volume (AV), essential for identifying potential abnormal evolutionary cases, remains largely unexplored. Our work filled this crucial gap by delineating the milestones and critical periods in maturation of the cerebrovasculature, revealing asynchronized patterns of cerebrovascular maturation, especially at the spatial scale of Brodmann area.

Different evolutionary trends between CVs and AVs

In this work, hierarchical normative models for CVs were built. Compared to a previous normative model for CVs15, the normative model of whole-brain cortical volume (WBCV) we built exhibited similar near-linear decreasing trends with respect to aging during adulthood. This decreasing trend was preserved across different hierarchies. Specifically, the normative models for bilateral Brodmann areas demonstrated similar near-linear decreasing trends at the whole-brain level. Evolutionary trends of CVs were also investigated in 34 bilateral brain regions defined in the Desikan–Killiany parcellation (refer to Fig. S8.1.2 of ref. 15), wherein similar preserved decreasing trends were observed, thus demonstrating the effectiveness of our pipeline for cortical (and arterial) quantitative analysis.

Comparatively, whole-brain arterial volume (WBAV) we got from 1.5T or 3T TOF-MRA volumes in this work also showed a similar near-linear decreasing trend (Fig. 2A) to that of WBCV (Fig. 2B). The TOF-MRA volume is acquired contrasting flowing blood against stationary tissues53, making it qualified for quantitatively measuring vascular morphology54 within large arteries55. If blood flow significantly decreases, especially in the arteries covering the cortex, the signal intensity captured by the MR will also decrease, leading to signal loss from the corresponding arteries. Hence, it may be argued that rather than reflecting aging-induced morphological variations, the decreasing trend we observed in AV was attributed to aging-induced declined in blood flow56,57. However, as examined in previous studies, using high-field (7T) TOF-MRA to improve detection ability resulted in quantitative measurements of small targets (e.g., aneurysms) in high agreement with those obtained using 3T58. Hence, aging-induced declined blood flow, if exists, still remained at a detectable level in 3T TOF-MRA images. Additionally, this argument is not supported by our findings that AVs in CoW regions exhibited a rising trend with respect to aging. At the small spatial scale, some normative curves from different Brodmann areas showed parabolic evolutionary trends, further negating the blood flow decline hypothesis. In fact, morphological losses of cerebrovasculature have been observed using contrast-enhanced MRA in normal aging animal models59. Therefore, the observed arterial evolutionary trends likely map the intrinsic morphological changes of cerebrovasculature.

A rising trend of AV was identified in CoW as shown in Fig. 2A. Morphological variations within CoW have been linked to lethal brain diseases60,61. Hence, understanding these variations and associated ischemic risks in clinical practice is necessary62 for effective patient care. However, a consensus has not been reached for age-related (positive63 versus negative64) arterial morphological differences within CoW, probably due to different experimental settings including sample size and data dimension. Our work, leveraging a large dataset, could serve as a reliable reference for evaluating the evolutionary patterns of CVs and AVs within CoW. Besides, our dataset supports further research at finer spatial resolutions to investigate the morphological and topological properties of different CoW segments.

Evolutionary correspondence between CVs and AVs

Our finding that CVs declined gradually throughout adulthood aligns with previous studies15,65. Additionally, our results support previous observations that the onset and extent of cortical tissue loss are region-specific17,65. However, it should be noted that aging-related cortical brain atrophy (as observed in Fig. S3) did not correspond to loss of neuronal number. It has been demonstrated that cortical thinning during normal aging is accompanied by the preservation of total neuronal number, especially in frontal and temporal regions66, suggesting that the reduction in neocortical volume loss with aging is more related to the loss of neuronal and dendritic architecture rather than the loss of neurons. The neuronal and dendritic architecture, i.e., the neuron’s dendritic arbor, determines the synaptic density67 which plays a central role in energy consumption within the brain68. Collapse of the neuronal and dendritic architecture would lead to the reduction of synaptic density, decreasing energy demand and subsequently, a proportional reduction in blood vessel69. This may be a primary reason of the general decreasing trend of AVs with respect to aging. Visualizations of regional evolutionary trends of AVs and CVs are shown in Supplementary Figs. S13, S14. Notably, AVs in Brodmann areas near the parietal and temporal lobes reached their peak at older ages compared to other brain regions. We employed a leave-one-out strategy to assess the stability of the normative model for each Brodmann area, as shown in Supplementary Figs. S20, S21. Both qualitative and quantitative assessments confirm that the normative model remains consistent across all leave-one-out rounds.

Complementary insights from CVs and AVs in brain pathologies

Our study investigated the hierarchical evolutionary patterns of CVs and AVs in subjects with various brain pathological conditions, including AD, LTS, and LS. For AD subjects, we observed a general collapse in cerebral vasculature, suggesting a widespread hypoperfusion and diminished blood flow throughout the brain. Specifically, AVs within the CoW were notably low for both genders, indicating significant hypoperfusion, accompanies by degenerative trends in CVs as well. These findings underscore the importance of AVs in AD pathogenesis, potentially outweighing CVs alone.

Comparative analysis between HC and stroke subjects revealed distinct patterns. Subjects with LTS showed moderate declines in CVs and AVs, while those with LS exhibited relatively milder decreases (refer to Supplementary Tables S11–S18). AD subjects, however, demonstrated the most pronounced reductions in both AVs and CVs across various brain regions, including the CoW, ACA, MCA, and PCA. These differences highlight the complementary roles of AVs and CVs in understanding cerebrovascular health and disease. Additionally, gender-specific differences were observed in the evolutionary trends of CVs and AVs. For instance, males with AD had significantly smaller AVs and CVs compared to healthy male controls across multiple regions. Females with AD similarly showed substantial reductions in both AVs and CVs, though the degree of change varied by region. These gender-specific trends further emphasize the necessity for personalized approaches in clinical diagnostics and interventions. Our findings fill a crucial gap in cerebrovascular research by providing hierarchical normative models of AVs and CVs, facilitating statistical inferences at the individual level. This work supports the development of automated pipelines for early detection and timely intervention of cerebrovascular abnormalities, ultimately advancing clinical practices in diagnosing and treating neurodegenerative diseases and cerebrovascular conditions in future.

Limitations

While our study provides valuable insights into the joint evolutionary patterns of CVs and AVs, several limitations should be noted. First, the deep-learning model’s performance may be influenced by variability in image quality across the multi-center dataset, potentially affecting segmentation accuracy. Second, while we analyzed a large cohort, the cross-sectional nature of the study limits our ability to infer causal relationships between cerebrovascular changes and brain pathologies. Additionally, the focus on a specific age range and disease conditions may not fully capture the complexity of cerebrovascular changes across different populations.

Methods

Construction of datasets

This work collected 2841 human cerebrovascular TOF-MRA volumes from four groups: healthy UK subjects, healthy Chinese subjects, Chinese subjects with brain pathologies, including AD, LTS, and LS. Figure S31 illustrates the data collection procedure. Volumes of healthy UK subjects were from two public datasets, IXI dataset (http://brain-development.org/ixi-dataset/) and OASIS-3 dataset43. IXI dataset totally contains 570 TOF-MRA volumes from multiple centers. The OASIS-3 dataset originally contains 2094 MR sessions from both cognitively normal subjects and subjects with cognitive decline (version 2022.04), with 556 MR sessions including TOF-MRA volumes. Nine hundred and fifty TOF-MRA volumes of healthy Chinese subjects were from the AHGMU in Guizhou Province, China, with informed consent from each participant and approval from the Institutional Review Board (IRB) of AHGMU. Cerebrovascular TOF-MRA volumes from subjects with brain pathology were collected from DGH in Fujian Province, China, with informed consent from all participants and IRB approval from DGH. The hospital provided 760 MRI scans, 642 of which included TOF-MRA imaging. One hundred and twenty-three TOF-MRA volumes of AD subjects were acquired from ZHTCM in Hubei Province, China, with informed consent from the participants and IRB approval from ZHTCM. For detailed acquisition parameters from each site, refer to Supplementary Table S22.

Data exclusion

Data exclusion was performed per the criteria shown in Supplementary Fig. S23. Specifically, 20 TOF-MRA volumes from IXI dataset were excluded due to poor image quality or severe artifacts via manual inspection. Additionally, 1 TOF-MRA volume without a corresponding good-quality T1 scan and 16 TOF-MRA volumes lacking subject age information were excluded. For the 556 TOF-MRA volumes from the OASIS-3 dataset, 44 volumes of poor image quality or severe artifacts were excluded via manual inspection. Using clinical diagnosis information and past medical history, 419 volumes of cognitively normal subjects were finally retained for healthy subject analysis. Among 950 volumes from healthy Chinese subjects, 21 TOF-MRA volumes of low quality were excluded. Thirty volumes without good-quality T1 images and 13 volumes missing age information were further excluded. For TOF-MRA volumes from subjects with brain pathology, we excluded 13 volumes of poor image quality or severe artifacts, 12 volumes without corresponding good-quality T1 scans, and 1 volume missing age information. For TOF-MRA volumes from AD subjects, after removing volumes of poor image quality or severe artifacts, 119 volumes remained.

Data annotations

We constructed a cerebral vessel segmentation dataset, CereVessMRA, to develop a deep-learning method for cerebrovascular segmentation. CereVessMRA dataset contains 271 TOF-MRA volumes of good image quality, with 28,128 slices in total. The whole CereVessMRA dataset involves 150 TOF-MRA volumes from healthy subjects and 121 TOF-MRA volumes from subjects with pathology. All TOF-MRA volumes were manually annotated cerebral vascular regions using an in-house LabelTool software (refer to LabelTool.docx) in the Supplementary Material for an introduction of this software.

For purpose of confirming annotation credibility, 3 clinicians with 1 year of clinical experience independently labeled TOF-MRA volumes from scratch to get coarse annotations during the first round. Annotating cerebral vessels in each volume from scratch took around 1 day. After the first-round annotation, clinicians modified these coarse annotations again to get preliminary annotation results. Supplementary Fig. S24 displays examples of bad annotations after the first round and good annotations after revision. Bad annotations include missing vascular regions, incomplete annotation of vascular regions, and incorrect annotations of vessel-like background tissues and bright tissues. Average time for revising these bad annotations for one TOF-MRA volume was about 3 h. Then, two radiologists with 5 years of experience double checked these preliminary annotations and revised them if necessary. Finally, the annotation was overlaid on every TOF-MRA volume for quality inspection to ensure precise annotation. Another 14 days were spent on this inspection and revision process. DSC score between the first-round coarse annotation and final qualified annotation ranges from 0.71 to 0.82, indicating the significant improvement through revisions and underscoring the necessity of these revisions.

Attention module-boosted model CereVessSeg

We developed a 3D convolutional neural network model, CereVessSeg, for automatic cerebral arterial segmentation in 3D TOF-MRA volumes. Detailed structure of CereVessSeg is shown in Fig. 6. Based on nnUNet28, the overall architecture of CereVessSeg is an encoder–decoder with six resolution stages and skip connections. Due to the morphological and intensity variations of cerebral arteries as well as noisy backgrounds, deep-learning model might produce segmentation errors, particularly false negative errors when segmenting thin vessels or false positive errors when segmenting some background tissues with low-contrast to vascular regions. To address this issue, we introduced visual attention mechanism and built an attention module (AM) in CereVessSeg, as displayed in Fig. 6. With AM, CereVessSeg automatically screens out non-discriminative low-level features extracted by the encoder path. The filtered low-level features are then integrated with high-level features, thereby valuable local spatial information and global context information are combined to produce fine-grained and accurate cerebrovascular segmentation.

Fig. 6: Detailed architecture of CereVessPro and CereVessSeg for cerebrovascular segmentation.
figure 6

The encoder and the first two stages of the decoder from the segmentation model CereVessSeg are initially pretrained using the contrastive learning method CereVessPro with a large-scale unannotated dataset. Afterward, the entire CereVessSeg model is finetuned using the annotated dataset.

Full size image

We employed a modified voxelwise cross-entropy loss to optimize CereVessSeg. Segmentation of cerebrovascular voxels that have low contrast with the background is particularly challenging for deep-learning models attempting to build complete arterial structures. Besides, a considerable number of background voxels with similar intensity or texture to arteries can degrade segmentation accuracy. These challenging voxels complicate effective learning for CereVessSeg. To mitigate this issue, we adopted a hard example mining strategy and developed the modified cross-entropy loss as follows:

$${L}_{{{\rm{mCE}}}}={L}_{{{\rm{L}}}1-{{\rm{CE}}}}+{L}_{{{\rm{Err}}}-{{\rm{CE}}}},$$
(5)
$${L}_{{{\rm{L}}}1-{{\rm{CE}}}}=- \sum _{c=0}^{1} \sum _{h=0,w=0,d=0}^{H,W,D}| | {{{\bf{y}}}}^{(c,h,w,d)}-{\hat{{{\bf{y}}}}}^{(c,h,w,d)}| {| }_{1}{{{\bf{y}}}}^{(c,h,w,d)}\log {\hat{{{\bf{y}}}}}^{(c,h,w,d)},$$
(6)
$${L}_{{{\rm{Err}}}{{-}}{{\rm{CE}}}}=- \sum_{h,w,d\in \,{\rm{FN}}}\frac{1}{| {\rm{FN}}\,| }{{{\bf{y}}}}^{(c=1,h,w,d)}\log {\hat{{{\bf{y}}}}}^{(c=1,h,w,d)} \\ - \sum_{h,w,d\in \,{\rm{FP}}}\frac{1}{| {\rm{FP}}\,| }{{{\bf{y}}}}^{(c=1,h,w,d)}\log {\hat{{{\bf{y}}}}}^{(c=1,h,w,d)},$$
(7)

where \(\hat{{{\bf{y}}}}\) represents the predicted segmentation result by CereVessSeg, and y represents the corresponding ground truth. Intuitively, a significant gap exists between the prediction and ground truth for a hard example, and the gap grows with the unidentifiable level of the example. Therefore, in LmCE, losses of voxels are reweighed in conformity to the L1 distance between prediction and ground truth. Moreover, losses of specific hard examples that cause segmentation errors get extra emphasis, in order to enhance the learning efficiency of CereVessSeg on these hard examples to realize accurate cerebrovascular segmentation.

Deep supervision mechanism was leveraged at multiple resolution stages of CereVessSeg’s decoder. As shown in Fig. 6, the final training loss of CereVessSeg is as follows:

$$L={L}_{{{\rm{Seg}}}}+\sum _{i=1}^{3}{L}_{{{\rm{Seg}}}}^{i},{L}_{{{\rm{Seg}}}}={L}_{{{\rm{DICE}}}}+\lambda {L}_{{{\rm{mCE}}}},$$
(8)

where LDICE is DICE loss, and λ is the weight coefficient of LmCE. \({L}_{{{\rm{Seg}}}}^{i}\) is the auxiliary losses of segmentation predictions from intermediate stages. LDICE is defined as follows:

$${L}_{{{\rm{DICE}}}}=1-\frac{2| {{\bf{y}}}\cap \hat{{{\bf{y}}}}| }{| {{\bf{y}}}|+| \hat{{{\bf{y}}}}| },$$
(9)

By exploiting deep supervision, CereVessSeg explicitly captures multi-resolution features of cerebral arteries, and refines segmentation from coarse to fine.

Contrastive learning of CereVessPro encourages arterial structural homogeneity

Conventionally, supervised-learning models for medical image segmentation suffer from scarcity of labeled training data, since dense annotation for medical image segmentation requires expert knowledge and careful identification. SSL presents as a desirable approach to relieve the large-scale annotation burden70,71,72,73. SSL leverages unlabeled data to learn underlying feature presentation by defining a surrogate task, which is different from but related to the supervised task. Meanwhile, SSL is capable of alleviating generalization errors without compromising model performance. Contrastive learning is a flourishing branch of SSL, which aims at grouping similar samples (positive pairs) closer while separating dissimilar samples (negative pairs) away74,75,76,77. Generally, a positive pair is created from two randomly augmented views of the same image, and a negative pair is created from views of two different images.

In this work, we employed contrastive learning to exploit information from unlabeled cerebral TOF-MRA volumes to boost cerebrovascular segmentation performance of CereVessSeg model as well as promote generalization ability. PixPro, a novel contrastive learning method proposed by Xie et al.76, encourages pixel consistency between positive pairs to boost pixel-level representation learning of 2D dense-prediction models. As PixPro does not require negative pairs, it is inherently well suited to pretrain 3D model for medical volume segmentation. In this work, we integrated cerebrovascular characteristic with PixPro to develop a contrastive learning method CereVessPro for our 3D segmentation network CereVessSeg. CereVessPro exploits position consistency of CvS, that is the same cerebral arterial branch tends to appear in certain area across different cerebral TOF-MRA volumes. Therefore, besides encouraging features of local voxels from two augmented views to be consistent as in PixPro, CereVessPro further emphasizes features of voxels in identical spatial regions across different TOF-MRA volumes to be close. We employed CereVessPro to pretrain CereVessSeg model on unlabeled MRA volumes.

Figure 6 displays the architecture of proposed CereVessPro. Concretely, for a batch of 3D cerebral TOF-MRA volumes \({{{\mathcal{X}}}}^{B}={\left\{{{{\bf{x}}}}^{i}\right\}}_{i=1}^{B}\), two augmented views are generated from each volume by random cropping and other image augmentation techniques to get two group of samples \({{{\mathcal{X}}}}_{1}^{B}={\left\{{{{\bf{x}}}}_{1}^{i}\right\}}_{i=1}^{B},{{{\mathcal{X}}}}_{2}^{B}={\left\{{{{\bf{x}}}}_{2}^{i}\right\}}_{i=1}^{B}\). \({{{\bf{x}}}}_{1}^{i}\) and \({{{\bf{x}}}}_{2}^{i}\) share equal size. The two views of every volume are separately inputted into a regular encoder network θ and a momentum encoder network ϕ to extract feature maps \({f}_{\theta }({{{\bf{x}}}}_{1}^{i})\) and \({f}_{\phi }({{{\bf{x}}}}_{2}^{i})\). The encoder network consists of a backbone network and a projection head network. In traditional way, only the encoder path in the segmentation network is used as backbone for pretraining. However, the final low-resolution feature maps extracted by encoder are insufficient for local representations during the contrastive learning. To address this issue, in our CereVessPro, we further enrolled the first two stages of decoder from CereVessSeg in the backbone module. The projection network is realized by two successive convolution layers with kernel size 1 × 1 × 1. We proposed to formulate the contrastive learning in a multi-scale manner. Feature of an example produced by the penultimate stage in the decoder path of backbone is also passed through a projection network. As a consequence, paired auxiliary features for two augment views are obtained, denoted as \({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i})\) and \({f}_{{\phi }^{{\prime} }}({{{\bf{x}}}}_{2}^{i})\). Then passing features extracted by the regular encoder network, i.e., \({f}_{\theta }({{{\bf{x}}}}_{1}^{i})\) and \({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i})\), through a voxel propagation module (VPM), the propagated feature maps \(P({f}_{\theta }({{{\bf{x}}}}_{1}^{i}))\) and \(P({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i}))\) are obtained.

According to PixPro, voxels from the same local area in an original MRA volume are defined as positive pairs. We adapted the pixel-to-propagation consistency loss proposed in PixPro to fit 3D representation learning of volumetric images, so as to encourage features of positive-pair voxels in MRA computed from two augmented views to be consistent. The loss is as follows:

$${L}_{{{\rm{VoxPro}}}} = - \sum _{i=1}^{B} \sum_{(m,n)\in {\Omega }^{+}}\{\cos [P{({f}_{\theta }({{{\bf{x}}}}_{1}^{i}))}_{m},{f}_{\phi }{({{{\bf{x}}}}_{2}^{i})}_{n}]+\cos [P{({f}_{\theta }({{{\bf{x}}}}_{1}^{i}))}_{n},{f}_{\phi }{({{{\bf{x}}}}_{2}^{i})}_{m}]\}\\ - \sum _{i=1}^{B} \sum_{(m,n)\in {\Omega }^{+}}\{\cos [P{({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i}))}_{m},{f}_{{\phi }^{{\prime} }}{({{{\bf{x}}}}_{2}^{i})}_{n}]+\cos [P{({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i}))}_{n},{f}_{{\phi }^{{\prime} }}{({{{\bf{x}}}}_{2}^{i})}_{m}]\},$$
(10)

where (m, n) is a positive voxel pair from two augmentation views, and Ω+ is the set of all these positive voxel pairs. Following PixPro, paired voxels with normalized distance smaller than threshold \({{\mathcal{T}}}=0.5\) are treated as a positive pair.

Since cerebral arteries within the same anatomical area of different subjects share structure similarities, positive pairs are supposed to include voxels from the same local area of different TOF-MRA volume as well. Given this idea, we further proposed LVessPro in our CereVessPro contrastive learning method to encourage features of voxels in the same area of different MRA volumes to be consistent. The formulation of LVessPro is as follows:

$${L}_{{{\rm{VessPro}}}}= - \sum _{i=1,i\ne j}^{B} \sum _{j=1}^{B}{\sum}_{({m}^{{\prime} },{n}^{{\prime} })\in {\Omega }^{{\prime}+}}\{\cos [P{({f}_{\theta }({{{\bf{x}}}}_{1}^{i}))}_{{m}^{{\prime} }},{f}_{\phi }{({{{\bf{x}}}}_{2}^{j})}_{{n}^{{\prime} }}]+\cos [P{({f}_{\theta }({{{\bf{x}}}}_{1}^{i}))}_{{n}^{{\prime} }},{f}_{\phi }{({{{\bf{x}}}}_{2}^{j})}_{{m}^{{\prime} }}]\}\\ - \sum _{i=1,i\ne j}^{B} \sum _{j=1}^{B} \sum_{({m}^{{\prime} },{n}^{{\prime} })\in {\Omega }^{{\prime}+}}\{\cos [P{({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i}))}_{{m}^{{\prime} }},{f}_{{\phi }^{{\prime} }}{({{{\bf{x}}}}_{2}^{j})}_{{n}^{{\prime} }}]+\cos [P{({f}_{{\theta }^{{\prime} }}({{{\bf{x}}}}_{1}^{i}))}_{{n}^{{\prime} }},{f}_{{\phi }^{{\prime} }}{({{{\bf{x}}}}_{2}^{j})}_{{m}^{{\prime} }}]\},$$
(11)

where \(({m}^{{\prime} },{n}^{{\prime} })\) is a positive voxel pair from two different subjects and \({\Omega }^{{\prime}+}\) is the set of all these positive voxel pairs. In this work, volumetric cerebral angiograms were already roughly aligned, and subvolume size of two augmented views covers a considerable area in the original MRA volume. Therefore, for LVessPro, there is no need to perform any registration step to get the exact positive pairs from two different TOF-MRA volumes. With LVessPro, the network learns effective cerebrovascular features, benefiting the downstream cerebrovascular segmentation task. Final training loss for our CereVessPro contrastive learning is:

$${L}_{{{\rm{CereVessPro}}}}={L}_{{{\rm{VoxPro}}}}+{L}_{{{\rm{VessPro}}}}.$$
(12)

Implementation details of CereVessSeg and CereVessPro

We used PyTorch framework to implement the cerebrovascular segmentation network CereVessSeg and the contrastive learning method CereVessPro. Network training and testing was conducted on a platform with Intel i9-10980XE CPU & GIGABYTE GeForce RTX 3090 with 23 Gb graphics memory and 64 Gb random access memory. Weights in CereVesSeg were first initialized by He normal initializer78. The CereVessPro contrastive learning to pretrain CereVessSeg was conducted utilizing all TOF-encoder volumes from healthy subjects and subjects with brain pathology, after filtering out low-quality volumes. Background border of zero value was cropped out from every volume. The whole encoder together with the first two stages of decoder from CereVessSeg were enrolled during pretraining. The contrastive learning was run for 400 epochs driven by a Layer-wise Adaptive Rate Scaling (LARS) optimizer with a cosine learning rate scheduler. The base learning rate was set to 1, and the minimum learning rate was set to 0.001. Batch size was set to 2. For the momentum encoder, the momentum value started at 0.99 and was increased to 1.

Using the pretrained results by contrastive learning, the whole CereVessSeg was later finetuned by supervised learning. CereVessSeg can also directly be trained from scratch, without the pretrain step. Volumes from healthy subjects and those from brain pathological subjects with their annotation from our CereVessMRA dataset were both randomly divided into 10 parts, 6 parts for training, 2 parts for validation, and the remaining 2 parts for testing. To reduce computational load as well as fit the GPU memory, sub-volumes of size 56 × 244 × 192 (axial × coronal × saggital) were randomly cropped for training. Other augmentation techniques, including mirroring, gamma, contrast, brightness, low-resolution stimulation, Gaussian noise addition, were further applied to sub-volumes before being input to network. Batch size was set to 2. Segmentation probability maps predicted by segmentation network were converted into binary segmentation masks by taking the argmax of channel-wise probability vector for every voxel. During testing phase, segmentation mask of every TOF-MRA volume was predicted using sliding-window inference, where the window size equaled to the sub-volumes size for training and the overlapping size was the half window size. The network training was driven by a stochastic gradient descent with a poly learning rate scheduler. The initial learning rate was set to 1 ×  10−2 and the power was 0.9.

Implementation details of compared methods

We employed several advanced deep-learning methods for 3D medical image segmentation, including 3D U-Net, V-Net, DenseVoxNet, VoxResNet, U-ception, DeepvesselNet, nnU-Net, UNETR, and SwinUNETR, to compare with our CereVessSeg. To ensure a fair comparison, all methods used the same patch size and batch size as CereVessSeg. Except for the initial learning rate and the number of epochs, all other hyperparameters followed the original papers. To tune the hyperparameters, all compared methods used the same training and validation sets as CereVessSeg. We tuned the initial learning rate by either increasing or decreasing it by a factor of 10 or 100. For 3D U-Net, V-Net, DenseVoxNet, VoxResNet, U-ception, and DeepvesselNet, we used an early stopping criterion, stopping training if the validation metric (DSC) does not improve for 20 consecutive epochs. For UNETR and SwinUNETR, we set the maximum number of epochs to 2000, which for nnUNet and CereVessSeg was set to 1000. The final model for each method was selected based on the highest validation metric achieved during training.

Registration between subject and reference space

In this work, the widely used MNI space was adopted as the reference space. A unified framework79 combining image registration (to obtain the deformation field between subject and reference space), tissue classification (to obtain the probability maps of gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) in subject space) and bias correction was adopted to create the population-level anatomical average template using T1 images from all the subjects. Individualized deformation fields for biprojection between subject (the individual T1 image) and reference space (the average template) were also derived. In the registration process, the Diffeomorphic Anatomical Registration Through Exponential Lie Algebra (DARTEL)80 was applied to precisely accommodate the aging-induced brain atrophy81. To balance the different resolutions among multi-center MRA volumes, the average template was upsampled to 0.5 mm3. Please refer to Supplementary Fig. S25 for a detailed procedure of registration between subject and reference space.

Quantification of arterial volume

To obtain subject-level cerebrovasculature in the individual space, each subject’s TOF-MRA volume, along with the predicted cerebrovascular segmentation volume, was first coregistered to the paired T1 image before the segmentation of vessel is performed. If the resolution of the original T1 image is low (e.g., thickness of slice is 6 mm in our case), an operation of resample is first performed before coregisteration between TOF-MRA and corresponding T1 volumes is performed. The individualized deformation field was then adopted for warping the atlas defined in MNI space back to individual’s space where quantification of vascular/cortical features is computed and statistical analysis is performed.

We computed cerebral arterial features hierarchically (1) from whole brain, (2) from four typical vascular regions, and (3) from 82 Brodmann areas. For every hierarchy, AVs were computed for investigating the normal aging effect on CvS.

Prior to quantifying vascular features for UK-Hset, the common scanning ranges of TOF-MRA volumes from multiple centers were computed, as each center had different axial scanning ranges for subjects. Specifically, TOF-MRA/T1 volumes were registered to MNI space, and voxels outside the scanning area were set to zero. We defined a binary mask for the scanning area in each TOF-MRA/T1 volume. A statistical probability map of these scanning areas was calculated, and voxels with a possibility >0.7 were included in the final binary mask representing common scanning range. Vascular features were then computed within the common scanning range to carry out further statistical comparisons.

Volume is a crucial attribute for assessing morphological changes in the brain across various neurological applications. Previous researches have shown that substantial changes in CV, WM volume or functional area volumes are associated with normal aging15, cognitive impairment and a range of diseases3,4,10,11,12. Intuitively, morphological alterations in cerebral arteries, such as stenosis, vasodilation and malformation, may be reflected in AV. Therefore, in this work, we adopted volume attribute to measures global and regional changes in cerebral arteries from an explicit aspect. In experiments, we calculated AVs of every subject by multiplying the number of voxels in segmented vascular regions by the voxel size. For a detailed description of the AV quantification process, please refer to Supplementary Fig. S25.

Quantification of cortical volume

In the registration step, the probability map of GM was obtained in individual space. We computed the CVs of each subject by the same manner to the calculation of AVs, wherein a threshold of 0.5 was used to binarize the probability map of GM. Correspondingly, CVs were also derived hierarchically (1) from whole brain, (2) from four typical vascular regions, and (3) from 82 Brodmann areas.

Statistical analysis

Two-sample Kolmogorov–Smirnov test82 was adopted to contrast the gender distributions between UK-Hset and CN-Hset. For each gender, two-sample t-test was applied to evaluate the statistical difference in age between these two groups. Besides, for the same gender from CN-Hset and UK-Hset, two-sample Kolmogorov–Smirnov test was also included to evaluate the difference of age distribution. To examine the varying impacts of brain pathologies, a two-sided two-sample t-test was applied to contrast age-matched hierarchical cortical (or arterial) volumes from different pathological group and HC group for each gender. Bonferroni test83 was applied for multiple-comparison correction in above experiments. For N independent comparisons, three significant levels, including pcorrected < 0.001 (for puncorrected < 0.001/N), pcorrected < 0.01 (for puncorrected < 0.01/N), and pcorrected < 0.05 (for puncorrected < 0.05/N), were reported in this work.

Annotation Tool

An in-house software was used for vessel annotation. In the in-house software, we implemented several convenient tools to improve the efficiency of vascular annotation. Readers can access this software in the Supplementary Software and find a detailed description in Supplementary Data 2. This software is also publicly available at https://github.com/YingChen7/CereVessSeg_with_CereVessPro.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.