We utilized brain MRI, clinical, and demographic data from a randomly selected sample of participants of the Lothian Birth Cohort 1936 Study (https://en.wikipedia.org/wiki/Lothian_birth-cohort_studies), a longitudinal study of cognitive aging comprising community-dwelling individuals from in and around Edinburgh born in 1936. All participants provided written consent to take part in the study under protocols approved by Lothian (REC 07/MRE00/58) and Scottish Multicentre (MREC/01/0/56) Research Ethics Committees. The methods used for MRI acquisition and clinical data in this cohort have been reported previously (28–30). Participants had their first brain MRI scan at the second wave of testing at a mean age of 72.6 years and subsequent MRI examinations spaced at 3-year intervals (30). For this study, we randomly (i.e., using a random number generator function in MATLAB R2019b) selected a sample of 32 participants with brain MRI available for the first three consecutive scanning waves (i.e., a total of 96 individual MRI scans) with high-quality image data. Vascular risk factors, which included presence or absence of hypertension, hypercholesterolaemia, diabetes mellitus, stroke, and history of cardiovascular disease, were self-reported at each wave (30).

All brain images were acquired at the Western General Hospital of Edinburgh in a GE Signa Horizon HDx 1.5-T clinical scanner (General Electric, Milwaukee, WI) following the research acquisition protocol described in Wardlaw et al. (29) and ensuring full-brain coverage. In brief, the 3D T1-weighted sequence (160 slices) had an acquisition matrix of 192 × 192 and a voxel size of 1 × 1 × 1.3 mm³. The axial T2-weighted fast spin echo sequence (TR/TE = 11,320/102 ms, 80 slices) was acquired with a 256 × 256 matrix and a voxel size of 1 × 1 × 2 mm³. The axial fluid-attenuated inversion recovery (FLAIR) fast spin echo sequence (TR/TE/TI = 9,000/140/2,200 ms, 40 slices) was acquired with a 256 × 192 matrix and a voxel size of 1 × 1 × 4 mm³. The axial T2^*-weighted gradient echo sequence (TR/TE = 940/15 ms, 80 slices) was also acquired with a 256 × 192 matrix, but with a voxel size of 1 × 1 × 2 mm³. None of the sequences had an inter-slice gap.

Each study participant's image data were processed and checked individually. For each participant, all image sequences from all waves were co-registered to the T2-weighted image acquired in the first scanning wave using rigid-body registration (linear, six degrees of freedom) through FSL-FLIRT (31). We used existent binary masks of WMH from the three waves and PVS in the centrum semiovale from the first imaging wave generated as previously published (23, 32, 33) (Supplementary Figure 1). In brief, for the first scanning wave (i.e., referred hereby as wave 1) WMH masks were generated using a multispectral method that combines T2^*-weighted and FLAIR images mapped in the red–green color space and quantised to facilitate a robust thresholding using minimum variance quantisation (32). For the other two waves, WMH binary masks were generated from a statistical analysis of the FLAIR-normalized intensities, seeking full compatibility with the semi-automatic approach applied in the first wave, and a higher level of automation: while in the multispectral approach the Gaussians that describe the intensity distribution of the tissue classes in the FLAIR image were adjusted by the observer aided by a combination with another T2-weighted-based image, in the other this adjustment was performed automatically. In brief, hyperintense voxels on FLAIR were identified by thresholding intensity values equal to 1.69^*standard deviation above the mean, using an adjusted implementation from Zhan et al. (34). The resulting hyperintense areas unlikely to reflect pathology (i.e., artifacts and cortex) were removed automatically using a lesion distribution template generated from the segmentation results of the first wave. Further refinement was achieved by applying Gaussian smoothing, followed by the removal of voxels with intensity values below 0.1 and z-scores below 0.95. All WMH masks were checked and manually edited to ensure the segmentation was as accurate as possible. Bland–Altman analyses yielded a mean WMH volume difference (SD) of 0.38 (1.29) ml and an intra-class correlation coefficient of 0.938, in 15 randomly selected individual datasets that were segmented with both implementation approaches. Voxel-wise reliability analyses comparing the WMH masks manually edited after applying both approaches show only scattered differences in the periventricular boundaries and Dice similarity coefficient of 0.6 (SD = 0.128), similar to the published inter-observer differences using the same method (35). No clusters (i.e., of three or more voxels) of differences were identified.

Segmentation masks for PVS in the deep normal-appearing white matter (i.e., excluding the internal and external capsules that are part of the region clinically defined as basal ganglia for the purposes of PVS identification) were generated using a computational method described previously (23, 36), on T2-weighted images. Succinctly, images were resampled from 256 × 256 × 80 to 256 × 256 × 160, using spline interpolation, to perform the PVS segmentation in images with 1-mm³ isotropic voxels. According to Wardlaw et al. (1) and Valdés Hernández et al. (37), PVS were identified as tubular structures with lengths between 3 and 50 mm. Tubular structures in the “isotropic” T2-weighted images were enhanced using the Frangi filter in its 3D version, optimized for this purpose with the parameters described by Ballerini et al. (36), which, according to Ballerini et al. (23), enhances structures with widths above 0.5 and below 2.5 voxels. The output from the filter was thresholded, binarised, and quantified in the region of interest. We used 3D connected component analysis to computationally identify and assess the PVS. Segmentation masks for PVS inside WMH were generated using the same method (i.e., thresholding the output of the 3D Frangi filter) but in a fused image obtained by subtracting the FLAIR image from the T2-weighted image after both being corrected for bias-field inhomogeneities using FSL-FAST (38) and their intensities being normalized. This is to discern whether the low-intensity voxels within the WMH in FLAIR correspond to the PVS-like hyperintense structures identified using the T2-weighted image. Combining sequences to increase distinction of PVS-like structures is a well-established procedure (17, 21). The images were individually checked for noise and other artifacts that would affect accuracy of PVS segmentation. To avoid errors in misclassifying small lacunes resembling PVS as PVS, neuroradiological reports were checked.

Adjacency/closeness (or not) of PVS with WMH clusters [minimum size 3 mm diameter according to Wardlaw et al. (1)] was recorded while/by visually inspecting all slices where PVS were segmented. Initially, the inspection was done in axial slices, and they were double-checked in all radiological orientations throughout the sample in all waves. All visual assessments, performed with MRIcron v1.0.20190902 (https://www.nitrc.org/projects/mricron/), were blinded to visual rating scores for PVS and WMH and participants' clinical and demographic information. These were repeated by the same observer to ensure a perfect intra-observer agreement [mean differences in PVS count ± 95% confidence interval are equal to −0.207 ± 1.637, intra-class correlation coefficient (ICC) 0.99917 for wave 1 and −0.034 ± 1.361, ICC = 0.99961 for wave 2, see Bland–Altman plot in Supplementary Figure 2]. Deep WMH clusters were defined as regions or voxel clusters of WMH not contiguous with the WMH located in the periventricular lining. Periventricular WMH caps surrounding the horns of the lateral ventricles were counted as deep WMH clusters if they extend more than 13 mm from the ventricular surface into the deep white matter (39). Cross-sectional topological relationships were assessed on the first scanning wave. Each deep WMH cluster was recorded and classified as either “close” or “not close” to a perivascular spaces (PVS) in the first scan in co-registered T2-weighted and FLAIR images after superimposing the PVS and WMH binary masks in different contrasting colors (Figure 1). Deep WMH were defined as “close” to a PVS if their segmentation was overlapping or contiguous with a PVS and “not close” if it was not overlapping or contiguous. Adjacent slices were carefully inspected to ensure not to miss or double-count any PVS. Each deep WMH cluster was also labeled, and the region of the brain it was found in was recorded—either frontal, parietal, temporal, or occipital, determined by examining lobar segmentations from a digital anatomical atlas (40).

To identify topological relationships longitudinally, the total number of deep WMH clusters was counted again. This time, each was classified according to how their morphology changed relative to baseline PVS between consecutive waves (e.g., waves 1 and 2 and waves 2 and 3) in each participant. A further mask that represented the change in WMH between waves was also overlaid in the FLAIR image to facilitate this assessment (Figure 2). These relationships were defined as follows: “increase around”—deep WMH cluster increased so that more voxels of the WMH mask became contiguous with a PVS; “increase close”—deep WMH cluster previously considered “close” to a PVS, increased in size, but no more of it became contiguous with a PVS (i.e., the adjacency between both—PVS and WMH—masks observed was owed to the same number of voxels as in the previous wave); “increase not close”—deep WMH cluster increased but it was still not contiguous with a PVS; and “no increase”—no visible change in the deep WMH cluster (Figure 2). To avoid WMH cluster change in size being influenced by differences in WMH segmentation algorithms, this was only considered if the change was perceived as nearly having doubled (or halved) or more the original cluster size. Again, all deep WMH clusters counted were recorded according to the lobar brain region.

A statistical analysis was performed in R (version 3.6.2). First, the counted data per overall individual brain scan were summarized to show the median count and distribution within the sample, per category. Cross-sectional and longitudinal topological relationships were analyzed by calculating percentages per category, both overall and per brain region. To allow a longitudinal analysis of the cross-sectional relationships, a change in count of deep WMH clusters “close” and “not close” to baseline PVS between waves was simplified into a binary variable where “1” indicates an increase in the number of deep WMH clusters and “0” indicates otherwise (i.e., no change or a decrease). A binary logistic regression with random effects was used to model the association between the change in the number of deep WMH clusters and their location relative to a PVS. Random effects accounted for variation between participants and for correlation within participants, i.e., an individual could contribute more than one row of data to the analysis. Data per overall brain were used, and separate models were performed for waves 1–2 and waves 2–3. In addition to unadjusted models, we evaluated models that corrected for sex at birth and vascular risk factors, including hypertension, hypercholesterolaemia, diabetes mellitus, and history of cardiovascular disease, with a maximum of three predictors per model. We estimated our sample size guided by the only study that visually assessed the topological relationships between WMH and PVS (26), which randomly selected 600 deep WMH clusters from a cross-sectional sample for analysis. In our sample, at baseline (i.e., only in wave 1), the total number of deep WMH clusters counted surpasses 600. Resources available for this study allowed us to process data from 29 individuals. Nevertheless, given the exploratory nature of this study, the complexity, and nature of the assessments, to ensure reproducibility, comparability, trustworthiness, and objectivity in the analyses, this estimate was considered appropriate (41). We consider the primary value of this study to be an indication of feasibility. Moreover, although the LBC 1936 Study, which provided data for the current analysis, enrolls 728 individuals in its first MRI scanning wave, it is inappropriate to use own data for a sample size calculation for an analysis with the same data (41).

To further evaluate whether the spatial relationship between deep WMH clusters and PVS could have been (or not) influenced by the PVS burden, we analyzed the R² value of the univariate linear regression models that had the number of PVS close to WMH clusters as independent variables and the WMH volume at each wave as a dependent variable, and we further compared these with the R² values of the models that have, instead, the PVS volume and the average PVS size as dependent variables. As R² measures the strength of the relationship between the set of independent variables and the dependent variable, if the spatial relationship between PVS and WMH is determined by the PVS burden we would expect higher R² values (i.e., less variance around the fitted line) in the models that relate the PVS burden (i.e., total volume, mean PVS size) to the WMH volume. In addition, we calculated the density of PVS “close” and PVS “not close” to WMH at the three scanning waves and conducted a paired t-test (two-tailed) at each wave comparing the densities of PVS “close” and PVS “not close”. If our hypothesis on the existence of a spatial relationship between deep WMH clusters and PVS is true, the PVS density (i.e., number per ml) close to deep WMH clusters will be greater than (and significantly different from) the PVS density not close to WMH.

Data from three participants were excluded due to inaccuracies in their PVS segmentations mainly due to the presence of motion artifacts in the T2-weighted images. One participant was not included in the models performed between waves 1 and 2 as no deep WMH clusters were counted in either of these waves. The mean age of the sample at baseline was 71.87 years (SD = 0.38). The follow-up MRI scans from waves 2 and 3 were obtained three and six years later, respectively (i.e., at a mean age ~75 and 78 years). Table 1 shows the baseline characteristics of the sample including sex at birth, vascular risk factors, and Fazekas visual rating scores of WMH.

Figure 3 illustrates the variations in the cross-sectional topologies observed. For deep WMH clusters classed as “close” to a PVS, most bordered a PVS with a small area of overlap (Figure 3A). Within deep WMH clusters that were classed as “not close” to a PVS, there was considerable variability in the distance separating them (Figure 3D). Many of these different topologies existed simultaneously in the same participant (Figure 4).

The median count of deep WMH clusters “close” to baseline PVS (i.e., PVS identified at the first scanning wave) was greater than the median count of those “not close” at all time points (Figure 5). The counts were not normally distributed in this sample, with several participants having zero or very low counts of deep WMH clusters (see individual data in Supplementary files). A higher percentage of WMH clusters in wave 1 were found close to PVS (Table 2). In the rest of the waves, this pattern was also observed (wave 1: 77% “close,” 23% “not close”; wave 2: 76% “close,” 24% “not close”; wave 3: n = 69% “close,” 31% “not close”; Figure 6, Table 2). Despite far fewer deep WMH clusters being counted in the parietal and temporal regions than in the frontal region (Table 2), the percentages of these that were found “close” to baseline PVS (i.e., PVS identified at the first wave) were higher than those found “not close” in these regions (Figure 7; e.g., wave 1—frontal: 75% “close,” 25% “not close”; parietal: 81% “close,” 19% “not close”; temporal: 86% “close,” 14% “not close”). The occipital region was not included in Figure 7 as only one deep WMH cluster was counted in wave 1.

Between waves 1 and 2, the number of deep WMH clusters changed in 28 participants. There was a median increase of 2 in deep WMH clusters spatially close to baseline PVS (i.e., PVS detected at wave 1) and 1 in those not close. Between waves 2 and 3, the number of deep WMH clusters changed in all 29 participants. The median change in number spatially close to a PVS increased in 10, while the median change for those “not close” increased in 6. Logistic regression models with random effects found no significant associations between progression (defined by an increase in number) of deep WMH clusters and proximity of these clusters to PVS (i.e., either “close” or “not close”) at wave 1. Table 3 shows the odds ratios, 95% confidence intervals, and p-values for all models performed.

Figure 8 shows the counts of changes in deep WMH clusters' morphology in relation to baseline PVS. The percentages of the different changes between waves 1 and 2 and waves 2 and 3 are displayed in Figure 9 and Table 4. A greater percentage of deep WMH clusters increased in size (waves 1–2: 70% increased and 30% did not increase; waves 2–3: 77% increased and 23% did not increase). Of those that increased, the most frequent change was an increase in size around the location of the already existent (i.e., “baseline”) PVS (waves 1–2: 73% “increased around,” 10% “increased close,” 17% “increased not close”; waves 2–3: 72% “increased around,” 10% “increased close,” 18% “increased not close”). When lobar regions were analyzed separately, these trends were also observed (Figure 10, Table 4). The occipital region was not included in Figure 10 because, as previously referred, only one WMH cluster was identified in this region at wave 1.

Although far fewer deep WMH clusters were counted in the parietal and temporal regions compared with the frontal region, we more frequently observed larger WMH clusters forming around multiple PVS in the parietal and posterior temporal regions (Figure 11). Deep WMH clusters in the frontal region were typically smaller in volume with a small area of continuity with a PVS.

The PVS volume and mean size did not seem to have influenced the spatial relationships found. As shown in the Supplementary data spreadsheet “PVSandWMHburden_analysis.xlsx,” the R² value for the linear regression between number of CSO PVS “close” to WMH and WMH volume for the scanning wave 1 was 0.484, for wave 2 it was 0.421, and for wave 3 it was 0.542. On the contrary, the R² value for the linear regression between CSO PVS volume and WMH volume for wave 1 was 0.0539, and between mean CSO PVS size and WMH volume, it was 0.049, thus indicating that approximately half of the variation in WMH volume at each wave can be rather explained by the number of PVS “close” to WMH instead of by the PVS burden. Similar results were obtained while evaluating the univariate association between the number of PVS considered “close” or “around” WMH that increased from waves 1 to 2 (R² = 0.581) and from waves 2 to 3 (R² = 0.512), and the WMH volume at waves 2 and 3, respectively. Paired t-test analysis of the density of PVS “close” vs. “not close” showed significant differences between both variables at all waves [wave 1: mean (SD) PVS “close” = 0.114 (0.103) ml⁻¹, mean (SD) PVS “not close” = 0.0349 (0.0406) ml⁻¹, p = 1.259 × 10⁻⁶; wave 2: mean (SD) PVS “close” = 0.140 (0.117) ml⁻¹, mean (SD) PVS “not close” = 0.0450 (0.0560) ml⁻¹, p = 2.443 × 10⁻⁷; wave 3: mean (SD) PVS “close” = 0.216 (0.163) ml⁻¹, mean (SD) PVS “not close” = 0.100 (0.113) ml⁻¹, p = 1.133 × 10⁻⁸).

To classify the WMH clusters into “big” or “small” while avoiding inter-/intra-observer bias, we looked at the deep Fazekas scores (provided as part of the Supplementary data) given by the neuroradiologist. As defined in Fazekas et al. (49), the deep WMH scores are equal to 1 if small punctate and very few, 2 if bigger and in higher number, and 3 if large and confluent with the periventricular WMH. From the whole sample, only one subject (subject no. 16) has deep Fazekas scores equal to 3 and two subjects (subject nos. 7 and 11) have deep Fazekas scores equal to 2. Therefore, it is unlikely that the spatial relationship between WMH and PVS would have been driven by the size of the WMH clusters.