Image processing approaches to enhance perivascular space visibility and quantification using MRI

Imaging the perivascular spaces (PVS), also known as Virchow-Robin space, has significant clinical value, but there remains a need for neuroimaging techniques to improve mapping and quantification of the PVS. Current technique for PVS evaluation is a scoring system based on visual reading of visible PVS in regions of interest, and often limited to large caliber PVS. Enhancing the visibility of the PVS could support medical diagnosis and enable novel neuroscientific investigations. Increasing the MRI resolution is one approach to enhance the visibility of PVS but is limited by acquisition time and physical constraints. Alternatively, image processing approaches can be utilized to improve the contrast ratio between PVS and surrounding tissue. Here we combine T1- and T2-weighted images to enhance PVS contrast, intensifying the visibility of PVS. The Enhanced PVS Contrast (EPC) was achieved by combining T1- and T2-weighted images that were adaptively filtered to remove non-structured high-frequency spatial noise. EPC was evaluated on healthy young adults by presenting them to two expert readers and also through automated quantification. We found that EPC improves the conspicuity of the PVS and aid resolving a larger number of PVS. We also present a highly reliable automated PVS quantification approach, which was optimized using expert readings.

Imaging the perivascular spaces (PVS), also known as Virchow-Robin space, has significant clinical value. Many recent studies have shown pathological alteration of PVS in a range of neurological disorders [1][2][3][4][5][6][7][8][9] . It is also believed that PVS plays a major role in the clearance system, accommodating the influx of CSF to brain parenchyma through peri-arterial space, and the efflux of interstitial fluid to the lymphatic system through peri-venous space 6,[10][11][12][13] . MRI is a powerful tool that enables in vivo, non-invasive imaging of this less-known glia-lymphatic pathway.
In the clinical practice, PVS is quantified based on the number of visible PVS on the axial slice of a T2-weighted (T2w) image that has the highest number of PVS in the region of interest 14 . This process can be laborious and error-prone, so efforts to improve efficiency and accuracy have been made by using a wide range of automatic or semi-automatic segmentation techniques, from classical image processing approaches to deep neural network modelling [15][16][17][18][19][20][21][22][23][24][25] . Park et al. used auto-context orientational information of the PVS for automatic segmentation 25 . Recently, Ballerini et al. showed that Frangi filtering 26 could robustly segment PVS by extracting the vesselness map based on the PVS tubular morphology 20 . More recently, Dubost et al. used convolutional neural network with a 3D kernel to automate the quantification of enlarged PVS 27,28 .
While these methods have improved the automated segmentation of PVS, less effort has been made to enhance the visibility of PVS through postprocessing means. A typical MRI session includes a variety of different sequences and utilizing different intensity profiles of these sequences can potentially improve PVS detection rate for both visual reading and automated segmentation. Combining MRI signal intensities has been used for other applications to achieve tissue-specific sensitivity. Van de Moortele et al. combined T1w and proton density images by dividing the latter by the former to improve signal non-uniformity at ultra-high field and optimize vessels visualization 29 . MRI multi-modal ratio was also used to map cortical myelin content [30][31][32] (for a systematic evaluation of the contrast enhancement via the combination of T1w and T2w see 33 ). Additionally, Wiggermann et al. combined T2w and FLAIR to improve the detection of multiple sclerosis lesions 34 .
Here we describe a multi-modal approach for enhancing the PVS visibility, which was achieved by combining T1w and T2w images that were adaptively filtered to remove non-structural high frequency spatial noise. Furthermore, we present an automated PVS quantification technique, which can be applied to T1w, T2w or the enhanced contrast. The efficacy of the Enhanced PVS Contrast (EPC) on both visual and automated detection is assessed and the reliability of the automated technique is examined.

Method
Four folds of evaluations were performed. First, the visibility of PVS in EPC was qualitatively and quantitively examined. Second, the visibility of PVS to expert readers was evaluated by comparing the number of PVS counted in EPC and T2w images. Third, PVS automatic counting was introduced and evaluated. Fourth, the reliability of the PVS automated quantification was assessed using scan-rescan MRI data.
MRi data. T1w and T2w images of the human connectome project (HCP) 35 were used in the analysis.
Structural images were acquired at 0.7 mm 3 resolution images. We used data from "S900 release", which includes 900 healthy participants (age, 22-37 years). The preprocessed T1w and T2w images [36][37][38] were used. In brief: the structural images were corrected for gradient nonlinearity, readout, and bias field; aligned to AC-PC subject space and averaged when multiple runs were available; then registered to MNI 152 space using FSL 39 's FNIRT. Individual white and pial surfaces were then generated using the FreeSurfer software 40 and the HCP pipelines 36,37 . Among HCP subjects, 45 were scanned twice with scan-rescan interval of 139 ± 69 days; these scans were used to assess the reliability of the automated PVS quantification. enhanced pVS contrast (epc). Figure 1 summarizes the steps required to obtain EPC. After preprocessing the data, T1w and T2w images were filtered using adaptive non-local mean filtering technique 41 . Non-local mean technique measures the image intensity similarities by taking into account the neighboring voxels in a blockwise fashion, where filtered image is . For each voxel (x j ) the weight (ω) is measured using the Euclidean distance between 3D patches. The adaptive non-local mean filtering technique adds a regularization term to the above formulation to remove bias intensity of the Rician noise observed in MRI. Therefore the expected Euclidian distances between two noisy patches N i and N j is defined as: The Rician noise of the MRI images, calculated using robust noise estimation technique presented by Wiest-Daessle et al. 42 , was used as the noise level for non-local filtering 41 . To preserve PVS voxels while removing the noise, filtering was applied only on high frequency spatial noises. This was achieved by using a filtering patch with a radius of 1 voxel, which removes the noise at a single-voxel level and preserves signal intensities that are spatially repeated 41 . Finally, EPC was obtained by dividing filtered images (i.e. T1w/T2w).
PVS visibility was qualitatively compared across T1w, T2w, and EPC images in white matter and basal ganglia. PVS conspicuity was also assessed by comparing the PVS-to-white matter ratio in EPC images with that in T2, which was shown to provide a higher PVS contrast compared to T1w 24 . A number of PVS and non-PVS white matter voxels were randomly selected across 10 subjects (more than 50 voxels per subject) and the average PVS-to-white matter ratio was measured. expert reading and clinical scoring. PVS were independently rated in 100 randomly selected subjects of HCP by two expert readers (GB and NSB) on axial T2w and EPC images using a validated 5-point visual rating scale 43 in basal ganglia and centrum semi-ovale (0: no PVS, 1: 1-10, 2: 11-20, 3: 21-40, and 4: >40 PVS). GB is a medical doctor with more than 5 years of neuroimaging research experience, and NSB is a neuroradiologist with 11 years of experience in radiology. Readers were blind to the rating results of each other. For validating the automated technique, the readers reached a consensus in few cases with different scores. One subject in which the PVS number was scored as ">100" by one of the readers was excluded from the statistical analysis. The total PVS score for each subject was calculated as the sum of the basal ganglia and centrum semi-ovale scores. The correlation between the number of PVS counted in T2w and EPC was calculated using Pearson correlation coefficient. PVS total counts were also compared using paired t-test. Lin's concordance correlation 44 was used to determine the concordance between the two raters. In addition, inter-class correlation (ICC) estimates and their 95% confident intervals were calculated based on a mean-rating (k = 2), absolute-agreement, two-way mixed-effects model, as recommended in 45 . The statistical analysis was performed using SciPy library (version 1.2.0) on Python 3 and MATLAB statistics and machine learning toolbox. P-values smaller than 1e-25 were reported as p = 0.
Subsequently, we applied Frangi filter 26 to T1w, T2w, and EPC images using Quantitative Imaging Toolkit 61 , which was implemented similar to 20 . Frangi filter estimates a vesselness measure for each voxel  s ( ) from eigenvectors λ of the Hessian matrix  of the image: Default parameters of α = 0.5, β = 0.5 and c were used, as recommended in 26 . The parameter c was set to half the value of the maximum Hessian norm. Frangi filter estimated vesselness measures at different scales and provided the maximum likeliness. The scale was set to a large range of 0.1 to 5 voxels in order to maximize the vessel inclusion. The output of this step is a quantitative map of vesselness in regions of interest, is taken to be the maximum across scales, as suggested in the original paper 26 . The range corresponds the specific levels in scale space that are searched for tubular structure feature detector. Thus, the outputs across voxels comprise vesselness measured across a range of filter scales.
In order to obtain a binary mask of PVS regions, the vesselness map should be thresholded. The binary mask enables automated PVS counting, volumetric, and spatial distribution analysis. Given that the vesselness value could vary across modalities, the threshold was optimized for each input image separately. We used the number of PVS counted by the experts for threshold optimization because of the absence of a ground truth. Vesselness values were standardized using robust scaling, in which values were scaled according to the inter-quartile range (IQR) to avoid the influence of large outliers: min Then, the binary image of PVS mask was obtained by thresholding  s ( ).
Then, the automated estimate of the total number of PVS was obtained by counting the number of connected components of the masked image P(s). Optimum thresholds were found by maximizing the concordance with the expert visual reading (concordance was used instead of absolute difference to optimize for global threshold and avoid biasing the threshold toward raters counting): a e a e t max ( , ) ( , ) t 0 where a and e are one-dimensional arrays of PVS counts across all subjects (n = 100), obtained from the automated (a) and expert (e) readings, respectively. Kandall's tau (τ) and Spearman's Rho (ρ) were used to measure concordance and correlation, respectively. Average of expert readings from EPC images were used for optimization. Throughout the manuscript, the optimum thresholds of 2.3, 2.7 and 1.5 were used for T1w, T2w, and EPC, respectively. After visual inspection, we noted that the imperfection of the white matter parcellation in periventricular and superficial white matter areas led to incorrected or missed PVS segmentation in white matter boundaries (an example is shown in Supplementary Fig. 1). Therefore, the subtraction of a dilated mask of ventricles from the PVS mask was applied to exclude the periventricular voxels and remove the incorrectly segmented PVS at the lateral ventricles-white matter boundary. After obtaining the final PVS mask, the number of PVS was obtained by counting the number of connected components of the PVS mask. Small components (<5 voxels) were excluded from automated counting to minimize noise contribution. The automated technique was applied on all subjects. Finally, one-way ANOVA was conducted to compare the effect of input image (T1w, T2w, and EPC) on the estimated number of PVS.
test-retest reliability. We also evaluated the test-retest reliability of the PVS quantification with MRI. Forty MRI data from HCP 35 that included scan-rescan were used for reliability analysis. EPC was derived, and the automated quantification pipeline was applied on the scan-rescan images (T1w, T2w and EPC images): identical parameters and threshold were applied on scan-rescan data. Then scan-rescan reliability was assed using ICC, Lin's concordance 44 and Pearson correlation analysis. ICC estimates and their 95% confident intervals were calculated based on a mean-rating (k = 2), absolute-agreement, two-way mixed-effects model.

Results
Evaluation 1: comparing the EPC with T1w and T2w. The PVS were more visible in EPC compared to T1w and T2w images (Figs 2 and 3, and Supplementary Fig. 2). The superiority of the EPC was evident in both white matter (Fig. 2) and basal ganglia (Fig. 3). EPC allowed the detection of PVS that were hardly identifiable in T1w and T2w (see yellow arrows in Fig. 2: PVS that could barely be spotted in T1w and T2w were evident in EPC). PVS were more conspicuous in EPC compared to T1w and T2w. Moreover, PVS-to-white matter contrast was significantly higher in EPC compared to that in T2w ( Supplementary Fig. 3: ~2 times higher).
There are a few cases, concordantly found by both readers, where PVS count is higher in T2w than EPC: in most of those cases, the difference was minimal (1-2 PVS) and did not change the PVS class. Only 4 cases had 6-10 more PVS in T2w than EPC. Upon revisit we noted that these are related to cases in which the detection of the PVS from noise by the visual reader was challenging (due to the small size of the PVS and given that those cases had particularly noisy T2 images). Readers confirmed that the EPC results for these cases are more reliable. ( 2) 2 7 (Fig. 6). Visual inspection showed that, as expected, the Frangi filter was able to detect the tubular structures of the PVS (Fig. 7). No statistical evidence was found to suggest the automated quantification of PVS using EPC is superior to those derived from T1w and T2w. PVS quantification (number of PVS) were significantly correlated across T1w, T2w, and EPC results (all at p < 0.0001) and all the automatic measurements reported a similar concordance level with the expert scores. Lin's concordance coefficient between automated PVS counts and expert scores was 0.81, the bias correction value (a measure of accuracy) was 0.88, and the Pearson correlation coefficient was 0.61   www.nature.com/scientificreports www.nature.com/scientificreports/ (p = 1.5e-11). The analysis of variance showed that the effect of input image (T1w, T2w, and EPC) on the number of PVS measured was significant F(2,294) = 48.56, p = 6.3e-19 which suggests that the automated technique for PVS quantification needs to be applied on the same image modality across study data.

Evaluation 4: test-retest reliability of automated quantification. Excellent test-retest reliability
was observed in PVS automated quantification regardless of the input image used (Fig. 8 and Table 1). Same thresholds, optimized on different subjects, were used for scan-rescan data. Lin' concordance coefficient between scan-rescan PVS measurements were 0.89, 0.94 and 0.83 for T1w, T2w and EPC, respectively. T2w images showed the highest concordance compared to other inputs. PVS measures were significantly correlated between scan-rescan images (r = 0.90, p = 2.8e-14, r = 0.95, p = 1.4e-20, and r = 0.85, p = 2.8e-11 for T1w, T2w and EPC, respectively). For all three inputs, no difference between the number of PVS in scan and rescan was observed.

Discussion
In this article, we presented a combined T1w-T2w approach to enhance the visibility of the PVS. EPC was utilized for both expert and computer-aided readings and was evaluated qualitatively and quantitatively. The proposed map (EPC) enhances the contrast and improves the conspicuity of the PVS, resulting in detection of a significantly larger number of PVS identified by expert readers. EPC benefits from the inverse signal profile of fluid on T1w and T2w images: when these images are combined together, a magnified PVS-tissue contrast can be obtained. EPC also uses a spatial non-local mean filtering technique, which has shown to be effective for mapping PVS 19 . We noted that even in the high-resolution T2w images of the human connectome project (0.7 mm 3 resolution), it is difficult to detect small PVS ( Supplementary Fig. 2), while they could be identified with this new technique.
According to the current standard visual rating scale for PVS 43 , most of the subjects in this cohort belonged to the class with the highest amount of PVS (category 4: >40 PVS). We noted that >80% of subjects on T2 are categorized as having a PVS class 3 (>50%) or 4 (>30%) according to the rating scale. Moreover, the inter-subject variability of PVS was large (standard deviation of ~15), despite the fact that subjects are all young healthy adults. This highlights the fact that the current rating approach has inherent limitations, because (1) a counting approach is highly dependent on the image resolution and quality and is difficult to generalize, (2) dichotomizing the PVS count reduces the statistical power and could underestimate the extent of variation, because considerable variability may be absorbed within each group 62 , (3) it does not consider the morphometric and spatial information of the PVS. In order to detect the intra-and inter-subject PVS alteration and ultimately to determine the role of PVS in different pathologies, there is a need to improve the current imaging and rating techniques.
With the advancement of the MRI technology, structural imaging in submillimeter resolution is achievable in a plausible time frame; for example, T1w MRI with 0.7 mm3 resolution can be obtained in 8 minutes at 3T 63 . Such imaging resolution enables to visualize PVS that were otherwise not apparent due to partial volume effect. With www.nature.com/scientificreports www.nature.com/scientificreports/ the added visibility of PVS comes the challenge of counting and mapping, as the visual rating becomes extremely laborious. Here we presented a pipeline that can be used to automatically map PVS.
The scan-rescan experiment showed that EPC is highly reliable, with no observed statistical difference across scan-rescan results. A trivial difference between scan-rescan PVS maps was observed, which are most likely due to (1) segmentation imperfection and image intensity differences of scan-rescan signal (e.g. due to subject motion) 64,65 , (2) normal physiological changes of PVS in the same subject, such as potential effects of time-of-day, sleep, and hydration on morphometric estimates of PVS [66][67][68][69] . It should also be noted that the preprocessing could affect the presence of the PVS (e.g. transformation of the MRI to a common space). Here we analyzed the data in the subject space, in which the MRI were AC-PC aligned using spline interpolation during the preprocessing and artifact correction steps 36 . In addition, T1w and T2w were co-registered, which also involves interpolation. These interpolations could affect PVS quantification. We noted that the automated PVS count was slightly higher in raw T2w images compared to the AC-PC aligned T2w images if the same threshold is used (the number of PVS was on average 1.1 higher in raw images compare to the aligned images, but not statistically significant p > 0.05). The optimum threshold, which depends on the quality of the data, therefore should be optimized based on the study data.
Whether PVS quantification is done by a neuroradiologist or automatically, an image with high PVS-tissue contrast is ideal. Currently, the modality of choice for PVS analysis is T2w, because it offers a higher contrast of PVS-tissue compared to T1w 24 . Yet, small PVS (<image resolution) are difficult to detect or separate from noise. Hence, PVS quantification is often limited only to enlarged PVS, despite the fact that physiological or pathological changes are expected to initiate in submillimeter scale 70,71 . In order to improve the mapping of the PVS, the MRI contrast of the PVS and the neighboring parenchyma should be increased. Besides the image processing approaches, PVS contrast can be enhanced through MRI technological improvement such as optimizing imaging sequence 24 and employing ultra-high field technology 19,24,25,72 .
Our quantitative PVS mapping and previous works showed that PVS can be mapped from an individual MRI contrast [15][16][17][18][19][20][21][22][23][24][25] . However, it should be noted that in the presence of pathology, additional image sequences (e.g. FLAIR) can be useful to discern PVS from pathological changes, such as white matter hyperintensities (WMH). Another advantage of a multi-modal approach for PVS quantification is the improvement of misclassification. For instance, vessels not surrounded by PVS are not easily distinguishable from vessels with PVS if only T1w is used, because both appear hypointense in this modality. T2w and EPC are able to solve this issue: in fact, in the absence of PVS, vessels appear hypointense in T2w and hyperintense in EPC, unlike vessels with PVS, which appear hyperintense in T2w and hypointense in EPC (see Supplementary Fig. 4). In addition, the correct identification of PVS can be achieved via the analysis of its morphometric characteristics, including size, shape, and anatomical location. Ballerini et al., for example, argued that Frangi filter ensures specificity in PVS segmentation given the tubular structure of the PVS.
It should also be noted that many current projects (e.g. human connectome project 35 ) already acquire both T1w and T2w modalities. Therefore, a technique that could benefit from the existing data could be valuable. We showed that in the presence of T1w and T2w, the combined contrast proved to enhance PVS visibility, which is of high clinical reading value. We did not aim to convey the message that one has to acquire both modalities to be able to resolve PVS. In fact, we noted that the automated segmentation was slightly less stable when EPC was used compared to T1w or T2w automated segmentation.
Two recent studies have also used multi-modal techniques for PVS segmentation and showed that it outperformed segmentation derived from a single modality 15,20 . Ballerini et al. applied segmentation on each modality and combined segmentations outputs using an AND operation. Boespflug et al. applied multivariate clustering, followed by morphometric filtering, to extract PVS from T1w, T2w, FLAIR and proton density images. It should be noted that the aim of these techniques was to improve the accuracy of the automated segmentation, but our study primarily aimed to propose a map that improves the visibility and detectability of the PVS, which can also make the visual scoring more accurate. Enhancement of PVS was also proposed using Haar transformation 19 and more recently using conventional neural network 73 , both on T2w images acquired using 7T MRI. Our approach focuses on enhancing the contrast of the PVS by combining T1w and T2w images and was optimized for 3T MRI,  www.nature.com/scientificreports www.nature.com/scientificreports/ which is more accessible in comparison with 7T MRI. Regardless of the differences, we anticipate that EPC could be used as an input to aforementioned techniques. Here we also performed test-retest comparison to analyze the reliability of PVS automated quantification.
In addition to its potential clinical relevance, mapping PVS can be useful to improve the accuracy of other quantitative MRI techniques as well, because these could be affected by the partial presence of PVS in image voxels. Recently, ignoring PVS could systematically affect how quantitative MRI measures such as diffusion tensor imaging (DTI) and spin echo dynamic susceptibility contrast (DSC) measures can be interpreted 74,75 . Such contribution and its potential confounding effect may be ameliorated if PVS is mapped and/or included in the analysis.
A limitation of multi-modal combination techniques is that it requires additional scan time and therefore is more prone to subject motion, which could negatively affect the co-registration. An interleaved acquisition was shown to be highly effective to ameliorate the co-registration issue 29 . Another limitation of EPC is that it requires the same image resolution for T1w and T2w. These sequences are often acquired in different resolutions, particularly in clinical practices (T2w are often acquired with thicker axial slices). Future investigations could focus on determining the extent to which the resolution difference between T1w and T2w affects the EPC quality and whether an intra-subject co-registration could amend this limitation.
A limitation of the automated quantification approach is that its performance depends on the parcellation accuracy. Imperfection of the brain parcellation could affect the automated quantification (see an example in Supplementary Fig. 1). While false positive PVS in periventricular area could be removed by applying a dilated mask of the ventricles, the false negative PVS in the superficial area of the white matter are more challenging to remove. While superficial white matter missed PVS is not expected to affect the count, it could affect the volumetric estimates. Further efforts are required to explore the effect of brain parcellation on PVS mapping or to build computational tools that minimize the parcellation dependency.

conclusions
Our combined T1w-T2w approach (EPC) has demonstrated to enhance the visibility of the PVS, resulting in improvement of PVS mapping. EPC allowed both the expert readers and the computer-aided algorithm to achieve a more accurate and precise quantification of PVS. This novel method, which can be easily applied to a number of MRI datasets, aims to overcome the limitations of current MRI sequences in PVS detection and quantitative analysis. This is relevant not only to better characterize the role of PVS when they are enlarged in pathological conditions, but especially to perform quantitative research on PVS when they are small, such as in physiological and prodromal states.

Data Availability
We have used human connectome project (HCP) dataset, which is already available to researchers.