MRI‐visible dilated perivascular spaces in healthy young adults: A twin heritability study

Abstract We investigated the narrow‐sense heritability of MRI‐visible dilated perivascular spaces (dPVS) in healthy young adult twins and nontwin siblings (138 monozygotic, 79 dizygotic twin pairs, and 133 nontwin sibling pairs; 28.7 ± 3.6 years) from the Human Connectome Project. dPVS volumes within basal ganglia (BGdPVS) and white matter (WMdPVS) were automatically calculated on three‐dimensional T2‐weighted MRI. In univariate analysis, heritability estimates of BGdPVS and WMdPVS after age and sex adjustment were 65.8% and 90.2%. In bivariate analysis, both BGdPVS and WMdPVS showed low to moderate genetic correlations (.30–.43) but high shared heritabilities (71.8–99.9%) with corresponding regional volumes, intracranial volumes, and other regional dPVS volumes. Older age was significantly associated with larger dPVS volume in both regions even after adjusting for clinical and volumetric variables, while blood pressure was not associated with dPVS volume although there was weak genetic correlation. dPVS volume, particularly WMdPVS, was highly heritable in healthy young adults, adding evidence of a substantial genetic contribution in dPVS development and differential effect by location. Age affects dPVS volume even in young adults, while blood pressure might have limited role in dPVS development in its normal range.

To understand the pathophysiology of dPVS and to define degree of influence by modifiable environmental determinants for PVS dilatation, knowledge of the possible genetic contribution for PVS dilatation would be helpful by assessing its heritability, which is defined as the proportion of variation observed within a phenotypic trait that is due to genetic variation. To date, only one study has investigated the heritability of dPVS with genome-wide genotyping in a large population-based elderly cohort (Duperron et al., 2018).
The authors of this previous study found differential heritability patterns between BGdPVS and WMdPVS, in which higher heritability was seen in WMdPVS than that in BGdPVS (79 vs. 49%) and the overlap of genetic effects between the two phenotypes, in other words, genetic correlation, was larger with WM hyperintensities in BGdPVS than in WMdPVS.
However, this past study was performed with an elderly population (mean age 72.8 years) and even though the study population had no history of stroke and dementia, it had a high prevalence of hypertension (76.6%) which is a well-known risk factor for dPVS (Duperron et al., 2018). Furthermore, dPVS was assessed with a subjective 4-grade scale, and this type of assessment is time-consuming and susceptible to inevitable intra-/inter-rater variability.
As heritability has a dynamic nature that is influenced by changes in gene expression and environmental exposure as well as various genetic-environmental interactions over time (Bergen, Gardner, & Kendler, 2007;Eaves, Long, & Heath, 1986), it is also important to understand that heritability in elderly people might be more influenced by environmental factors accumulated over a longer period of time than in young healthy adults. Moreover, a high prevalence of vascular risk factors such as high blood pressure (BP) and hypercholesterolemia could affect heritability estimates (Eaves et al., 1986). Therefore, our main goal was to investigate the heritability of PVS dilatation in healthy young adults objectively using a fully automated dPVS quantification method. In doing so, we used the Human Connectome Project (HCP) dataset which consisted of twins and nontwin siblings with various clinical information and high-resolution three-dimensional (3D) structural MRI available.

| METHODS
This retrospective study was approved by our institutional review board. The requirement for informed consent was waived as a publicly available dataset was used.

| Subjects
Heritability analysis was performed using structural images and demographic data included in the WU-Minn Human Connectome Project dataset which recruited 1,200 healthy adult twins and nontwin siblings 22-35 years old to characterize the relationships between brain circuits, behavior, and genetics (Van Essen et al., 2013).
Among 1,206 individuals available in the March 2017 (S1200) release, 1,113 individuals who had 0.7 mm isotropic 3D T1-and T2-weighted images obtained on 3 T MRI were initially screened. Individuals with the following criteria were excluded: (a) severe neurodevelopmental, documented neuropsychiatric, or neurological disorders; (b) high BP, diabetes, or significant cardiovascular disease; (c) born preterm-defined as born before 34th weeks of gestation for twins and before 37 weeks of gestation for nontwins; (d) zygosity not verified by genotyping; and (e) nontwin siblings without their respective pairs or with different biological fathers or mothers. Detailed information about recruitment-as well as the inclusion and exclusion criteria-of the Human Connectome Project is described in a previous study (Van Essen et al., 2013). For this study, a total of 700 subjects (350 twin or nontwin sibling pairs), consisting of 138 monozygotic (MZ) twin pairs, 79 dizygotic (DZ) twin pairs, and 133 nontwin sibling pairs, were finally included in the heritability analysis.
To assess the reliability of the automated dPVS segmentation method, 45 subjects who underwent two MRI acquisitions on different dates (mean interval of 134.8 days [range, 18-328 days]) were also included from the WU-Minn Human Connectome Project Retest dataset.

| Imaging acquisition
All MRI data were obtained at Washington University in St. Louis, MO, on a Siemens 3T MR scanner "Connectome Skyra," which was customized with a 100mT/m gradient coil and inner bore diameter of 56 cm, using a standard 32-channel head coil.
More detailed information regarding imaging protocols can be found in the WU-Minn Human Connectome Project S1200 Release Reference Manual (Van Essen et al., 2013).

| Automated dPVS segmentation method
When measuring dPVS volume, 3D T1-and T2-weighted images were used along with Freesurfer generated sub-region masks from the Human Connectome Project database. dPVS volumes were calculated fully automatically through the following two steps: 2.3.1 | Potential dPVS voxel extraction from T2-weighted images Potential dPVS voxels were extracted from T2-weighted images using 3D Frangi filtering (Frangi, Niessen, Vincken, & Viergever, 1998), which calculates the vesselness measure based on the eigenvalues of the Hessian. This method is known to be particularly useful for delineating enhancing vascular structures (Jimenez-Carretero, Santos, Kerkstra, Rudyanto, & Ledesma-Carbayo, 2013;Xiao et al., 2011) and also for detecting dPVS on medical images (Ballerini et al., 2016;Zong, Park, Shen, & Lin, 2016). To reduce the sensitivity of the signal intensity variations across subjects, the signal intensity of the T2-weighted image was normalized by subtracting the mean of the signal intensities, and then divided by the standard deviation of signal intensities before the Frangi filter was applied. After normalizing signal intensity, to calculate the vesselness measure, a 3D Frangi filter was applied to the T2-weighted image and thresholding was performed to extract potential dPVS voxels. The processed T2-weighted images had a large number of false positives outside the brain region. To reduce these false positives, the WM and BG masks from the Freesurfer segmentation results were used.

| Refinement by convolutional neural network-based classification
To reduce the remaining false positives around the ventricles or subarachnoid space, we trained a 3D deep convolutional neural network (CNN) to classify three classes (Class 1: dPVS, Class 2: false positives around the ventricles, and Class 3: false positives around the subarachnoid space) using a 24 × 24 × 24 3D patch as input data as described in Figure 1.
A kernel size for 3D convolution was 3 × 3 × 3 for all convolutional layers and the number of nodes for the fully connected F I G U R E 1 Workflow diagram of fully automated dilated perivascular space (dPVS) segmentation. (a) Extraction of potential dPVS voxels using threedimensional Frangi filtering and thresholding with subsequent falsepositive reduction via white matter (WM) masking, and (b) three-dimensional convolutional neural network (CNN) machine learning classifiers with a 24 × 24 × 24 three-dimensional input data patch. FP, false positives; ML, machine-learning; SAS, subarachnoid space layer was 128. A rectified linear unit activation was applied after each convolution layer and 3D max-pooling (2 × 2 × 2) was performed after two consecutive convolutional layers. The number of channels for the convolutional layers was 8, 16, 32, and 32 in order. A crossentropy loss between the network output and the label was used to train the network. For this purpose, three-class labeling on the processed images in Section 2.3.1. was manually performed for 20 randomly selected nontwin subjects. To increase the robustness of the input variations, several data augmentation steps such as flipping, adding randomly generated noises or biases, and multiplying randomly generated scaling factors were applied for each patch data. To reduce the effects of the unbalanced numbers between individual classes, each class was fed to the network at uniform frequency in the training process. The trained network was used to generate the final dPVS voxels from the potential dPVS voxels for the test dataset. Based on the output of the 3D CNN algorithm, regional dPVS voxels were counted from the BG and WM segmented masks. The volume of dPVS was calculated by the product of the number of dPVS voxels and volume of 1 voxel (i.e., 0.7 × 0.7 × 0.7 = 0.343 mm 3 ).

| Statistical analysis
2.4.1 | Reliability assessment of the automated dPVS segmentation method The reliability of the automated dPVS segmentation method was assessed using the intraclass correlation coefficient (ICC) for the automatically measured dPVS volume in 45 subjects with test and retest datasets.

| Validity assessment of the automated dPVS segmentation method
For validity assessment, fractions of false-positive and false-negative volumes from the automated dPVS segmentation method were compared with those from manual segmentation according to the STandards for ReportIng Vascular changes on nEuroimaging (STRIVE) criteria (Wardlaw, Smith, Biessels, et al., 2013). Axial T2-weighted images of 10 randomly chosen subjects were visually assessed by two board-certified radiologists (BLINDED and BLINDED with 6 and 13 years of experience in neuroradiology, respectively) who were blinded to the volumetric results of the automated dPVS segmentation method. dPVS were defined as lesions showing high signal intensity on T2-weighted images similar to that of cerebrospinal fluid and following the course of penetrating vessels with a diameter of less than 3 mm. dPVS in BG and WM were manually segmented on the axial images and their volumes were calculated at the level of anterior commissure and at the upper end of corpus callosum, respectively, first independently and then in consensus by the two neuroradiologists. ICCs between the manually and automatically measured dPVS volumes were also calculated.

| Group comparisons of clinical and volumetric variables and their association with dPVS volumes
According to the results of the normality tests, the clinical and volumetric variables of each group were demonstrated as medians (interquartile ranges) for continuous variables and numbers (percentages) for categorical variables. The continuous and categorical variables were compared among groups using linear and generalized linear mixed models, respectively, thereby accounting for the within-pair phenotypic correlations. In these models, subject group was set as a fixed effect while twin pair (i.e., Family ID provided by Human Connectome Project dataset) was set as a random effect. P values for pair-wise group comparisons were corrected for multiple comparisons using the Bonferroni method. To explore possible genetic influences, within-pair differences of clinical and volumetric variables were compared between groups via the Kruskal-Wallis test with posthoc Tukey's test for continuous variables and the Chi-squared test with Bonferroni correction for categorical variables. To explore factors that may possibly affect dPVS volume and adjust for their effects on the following heritability analyses, linear mixed models were used for BGdPVS and WMdPVS with each variable set as a fixed effect and twin pair (i.e., Family IDs provided by Human Connectome Project dataset) set as a random effect. Age and sex were included as covariates in these models except when analyzing the variables of age and sex themselves. For a simplified statistical analysis, we dichotomized ethnic groups into White and nonWhite, with the latter group including Black or African American and Asian, natural Hawaiian, or other Pacific Islanders. Subjects who did not know or report their ethnic group or belonged to more than one ethnic group were excluded when the ethnic group was analyzed as a covariate.

| Heritability analysis using univariate and bivariate models
We implemented the classical twin study design, in which narrowsense heritability (h 2 ) was defined as the ratio of variance from a phenotypic measurement, which can be shown as additive genetic varia-

tion [A]/total observed variation (genetics [A] + shared environment
[C] + unique environment [E]). Since the common or shared environment is considered to be the same within a family, [C] was set as 1 for all twins and nontwin siblings (Patel et al., 2017;Schneider et al., 2019). The genetic variation [A] was set as 1 for MZ twin pairs and 0.5 for DZ twins and nontwin siblings based on previous literature assuming that DZ twins and nontwin siblings share on average 50% of their genes (Patel et al., 2017(Patel et al., , 2018 In the univariate analysis, both BGdPVS and WMdPVS volumes were first adjusted for age and sex, hereafter named as unadjusted. In addition to the adjustment for age and sex, the volumes were sequentially adjusted for (a) regional volume (BG volume for BGdPVS and WM volume for WMdPVS) or (b) ICV-termed as partially adjusted.
For completely adjusted univariate models, age, sex, regional volumes, and ICV were adjusted. In addition to these covariates, pulse pressure and ethnic group were also included as additional covariates to assess their influence. We assessed variable inflation factors to ensure that there was no issue with multicollinearity in our models.
In the bivariate analysis, after adjusting for age and sex, genetic correlation and shared heritability of dPVS volumes with pulse pressure, regional volume, ICV, and counterpart dPVS (BGdPVS to WMdPVS and vice versa) were calculated, thereby accounting for the shared genetic variation between the two phenotypes. Genetic correlation represents the measurement of genetic overlap between two traits, which indicates the similarity of the genetic variation between two phenotypes. Genetic correlation ranges from −1 to 1, with a correlation of 1 implying that the two traits in question are influenced by the same genetic factors, while a correlation of 0 represents independent genetic effects. Shared heritability assesses the amount of genetic variation shared between two traits; that is, the contribution of the shared genetic effect to the observed phenotypic correlation.
The Cholesky decomposition ACE models enable the calculation of genetic correlation and shared heritability.
From the models, we could obtain three paths each for genetics, shared environment, and unique environment: a 11 denoted the path from the first latent genetic variable to the first trait, a 12 the path from the first latent genetic variable to the second trait, and a 22 the path from the second latent genetic variable to the second trait; c 11 , c 12 , and c 22 denoted the shared environment and e 11 , e 12 , and e 22 the unique environment. Genetic correlation and shared heritability were then calculated as: genetic correlation = a 12 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 12 2 + a 22 2 p shared heritability = a 11 a 12 j j a 11 a 12 j j+ c 11 c 12 j j+ e 11 e 12 j j ð Þ For shared heritability, absolute values of the path coefficients were used to calculate the proportion of genetic contribution against the total contribution, which was calculated as the sum of each factor's contribution (Coan, 2003).
For both univariate and bivariate heritability analysis, full ACE models were compared to CE models, which did not include genetic factors, to measure genetic influence. Furthermore, the full ACE models were also compared with AE and E models and corresponding Akaike Information Criteria were calculated.
All statistical analyses including heritability analyses were performed using R statistical software (version 3.6.2, Vienna, Austria). All twin heritability analyses were performed using OpenMx (version 2.12.2) and umx (version 2.9.9) packages within R, which allowed structural equation modeling to estimate heritability (Neale et al., 2016

| Characteristics of the study cohort
The baseline characteristics of all subjects are summarized in Table 1.
Group-wise comparisons of the clinical and volumetric variables among the three groups are summarized in

| Group comparisons of within-pair differences for subject characteristics
The nontwin siblings had greater differences in age and sex than MZ or DZ twin pairs (all, p < .001). The differences in volumes of BGdPVS, WMdPVS, BG, and WM, and ICV were greatest in nontwin siblings followed by DZ and MZ twin pairs in that order (all, p < .001). The within-pair differences in diastolic BP and HbA1c levels were also significantly different between groups (p < .001 and p = .004, respectively). No significant within-pair differences were found for the other characteristics (Table 2).

Representative images and correlation plots of BGdPVS and
WMdPVS volumes between twin and nontwin siblings are shown in

| Effect of clinical and volumetric variables on dPVS volumes
Older age and male sex were significantly associated with larger dPVS volume (both, p < .001). The effect of sex was no longer significant after adjusting for ICV (p = .503 for BGdPVS; p = .143 for WMdPVS), while the effect of age remained significant after adjusting for ICV or all other variables (Ps < .001). The White ethnic group adjusted for age and sex was associated with larger BGdPVS than the nonWhite group (p = .004). However, after adding ICV as a covariate, the significance disappeared (p = .114). All volumetric variables adjusted for age and sex also significantly affected dPVS volume (all, p < .01). However, none of the BP variables adjusted for age and sex were associated with dPVS volume (Tables S4 and S5). Therefore, in subsequent heritability analyses we adjusted for age, sex, regional volumes, ICV, and ethnic group, all of which were significantly associated with dPVS volume. Besides, we also included pulse pressure as a covariate in heritability analyses as hypertension and arterial pulsation are wellknown key factors in the development of dPVS.

| Univariate heritability analysis
The results of the univariate heritability analyses are shown in ( Tables S6 and S7). The heritability estimates of the White ethnic group alone also resulted in slightly lower estimates but overall were within an acceptable range (Table S8). The standardized path coeffi-

cients of genetic [A], shared environment [C], and unique environment
[E] factors for dPVS volumes in each model were plotted in Figure S1.
In all analyses, heritability measures were significant when the full ACE models were compared to the model without the genetic factor (CE) (p < 0.001; Table S9).

| Bivariate heritability analysis
The results of the bivariate heritability analyses are shown in Table 4.
As in univariate analysis, heritability measures were significant when the full ACE models were compared with the CE model without the genetic factor (p < .001;  Figure S2.

| DISCUSSION
The current study assessed the heritability of dPVS volume in healthy young twins and nontwin siblings in a well-controlled large dataset with high quality images. There were three main findings in this study.
First, our results support a previous study's finding by also showing highly heritable dPVS burden (Duperron et al., 2018). In the current study, the heritability estimates of both BGdPVS and WMdPVS volumes were high, and even higher than those reported in the elderly cohort (Duperron et al., 2018). Second, the heritability estimates were higher in WMdPVS volume than in BGdPVS volume, supporting differential genetic influences on dPVS by location. Lastly, among the well-known risk factors for dPVS, age was significantly associated with dPVS volume even in young subjects, but the effect of BP on dPVS volume was limited in its normal range.  in clinical settings, but is time-consuming, rater dependent, and prone to ceiling and floor effects, and these inherent limitations prevent the precise assessment of dPVS burden and its heritability. Therefore, we developed a fully automated segmentation algorithm for the quantitative and objective evaluation of dPVS burden. Although, the dPVS volumes measured by our automated method were slightly smaller than manually segmented dPVS volumes in some subjects, the automated method revealed results strongly consistent with those of manual assessment and without false-positive findings. Moreover, the automated method showed excellent reliability in 45 subjects who repeated MRI scanning within relatively short intervals. Therefore, this study demonstrated that the automated segmentation algorithm could accurately and reliably assess dPVS burden in the HCP dataset.
When comparing within-pair differences, differences in dPVS volume between the MZ twin pairs were the smallest, followed by those between the DZ twin pairs and nontwin siblings, which suggests genetic effects on dPVS volume. In line with this result, the heritability estimates were high in both BGdPVS and WMdPVS in our healthy young cohort. The heritability estimates were higher than those reported in a previous study on elderly subjects that had a high frequency of vascular risk factors [BGdPVS (current: 61.0-65.8% vs. previous 40%); WMdPVS (current: 87.4-90.3% vs. previous: 79%)] (Duperron et al., 2018). This difference might be due to the relatively brief exposure of healthy young adults to environmental impacts, allowing genetic contribution to remain significant in the younger population. However, as genome-wide genotyping which was used in a previous study (Duperron et al., 2018) reveals lower heritability estimates compared to twin studies (Visscher & Goddard, 2015), findings of these studies and ours need to be compared and interpreted with caution. Future studies with the same heritability analytic methods may be able to draw more convincing conclusions on this issue.
Similar to the previous study (Duperron et al., 2018), the heritability of the dPVS volume was higher in WM than in BG in our study, suggesting a distinct genetic influence on dPVS according to location.
Growing evidence indicates different risk factors and clinical manifestations for BGdPVS and WMdPVS (Andreas Charidimou et al., 2017;Duperron et al., 2018;Martinez-Ramirez et al., 2013;Roher et al., 2003). BGdPVS burden has been associated with hypertension, hypertensive angiopathy, and vascular cognitive impairment (Chen et al., 2019;Hansen, Cain, Thomas, & Jackson, 2015), while WMdPVS burden has been associated with cerebral amyloid angiopathy and Alzheimer's disease (Chen et al., 2019;van Veluw et al., 2016). In Alzheimer's disease and cerebral amyloid angiopathy, amyloid β deposition mainly occurs in the arteries of the cerebral cortex and leptomeninges (Vinters & Gilbert, 1983), potentially increasing arterial stiffness and compromising normal arterial pulsation in these areas, while brain regions commonly affected by hypertensive arteriopathy (i.e., BG, brain stem and cerebellum) (Andreas Charidimou et al., 2017) are usually free from amyloid β deposition (Mandybur, 1975). As the arterial pulsation is a key driver of fluid movement through PVS (Mestre et al., 2018), arterial amyloid β deposition could be a possible explanation for the fluid stagnation and PVS dilatation in WM near the cortex. Likewise, the arterial stiffness in BG caused by hypertension, might, at least in part, cause BGdPVS. Given that Alzheimer's disease shows higher heritability (60-80%) than hypertension (30-60%; Gatz et al., 2006;Kupper et al., 2005;Scherbakov, Von Haehling, Anker, Dirnagl, & Doehner, 2013;Shih & O'Connor, 2008), the distinct heritability of associated diseases might be related to the higher heritability of WMdPVS than BGdPVS. Furthermore, although there have been conflicting data, a previous study even found a direct association between WMdPVS and the apolipoprotein E ε4 allele (Roher et al., 2003). Meanwhile, the lower heritability of BGdPVS might also be due T A B L E 4 Bivariate shared heritability estimates (shared h 2 ) and genetic correlations (r g ) between dPVS and regional volume, intracranial volume, pulse pressure, and different regional dPVS to its higher susceptibility to environmental factors. Although the underlying process is still unknown, the anatomical differences of PVS according to location might be responsible for its distinct vulnerability to the environment. BGdPVS has two layers while WMdPVS has one for the pial membrane and surrounds more proximal arterioles while WMdPVS has smaller arterioles or capillaries (Pollock, Hutchings, Weller, & Zhang, 1997;Zhang, Inman, & Weller, 1990). As proximal arterioles have looser endothelial tight junctions than smaller arterioles or capillaries (Bechmann, Galea, & Perry, 2007), they may be affected more by endothelial damage induced by hypertension or other environmental insults .
Moreover, BGdPVS directly communicates with the adjacent basal subarachnoid space (Bouvy et al., 2014) and lenticulostriatal arteries are exposed to higher BP than cortical arteries (Yamaguchi, Fukuyama, Yamauchi, & Kimura, 1991), which might also make BGdPVS more vulnerable to hypertension or toxic materials circulating in the cere-  (Lau et al., 2017). In line with this result, our study also found larger BGdPVS volume in the White ethnic group after adjusting for age and sex. Ethnic differences in brain size (Hsu et al., 2018;Rushton & Ankney, 1996)  was short-as we demonstrated that age affected dPVS dilatation even in young adults-the dPVS volumes might differ between the two MRI scans. In our study, dPVS volumes on the second scan were larger than the first scan (mean volume increases were 0.013 cm 3 for BGdPVS and 0.708 cm 3 for WMdPVS, data not shown). We believe that this difference might have lowered the reliability of the automated segmentation method, and eliminating this confounding factor by having MRI scans performed on the same day would further confirm the reliability of our algorithm. Third, the automatic segmentation algorithm was built based on T2-weighted images only. Although our algorithm did not degrade its function in our young healthy population, a future study employing multiple sequences including FLAIR images for the automatic segmentation algorithm might help differentiate dPVS mimicking WM hyperintensities or old lacunes, and therefore, improve the dPVS segmentation accuracy, particularly in the elderly population. Fourth, heritability is a ratio, namely, a relative measure. It is affected by variations in cohort and does not directly reflect genetic determination. Lastly, we had an issue with negative path coefficients on bivariate heritability analysis. Although we adopted a previously reported method to calculate shared heritability using absolute values (Coan, 2003), there is not enough research on how to handle and interpret these negative values in heritability analysis (Steinsaltz, Dahl, & Wachter, 2017). As seen in our study, a few reports have also demonstrated negative path coefficients (Kendler, Myers, Dick, & Prescott, 2010;Tanguay-Garneau et al., 2020). Moreover, a recent study has suggested that while these negative values might be random noise or results of model misspecification, they might also reflect a real physical feature of the biological process (Steinsaltz et al., 2017). We hope our study can motivate future research to unveil the significance of these negative values. Despite these limitations, we believe that our results might support the substantial genetic contribution of PVS dilatation.
In summary, the present study provides evidence for a substantial genetic contribution to dPVS burden in young healthy adults. Older age was associated with higher dPVS burden even in young adults,