Arterial spin labelling perfusion MRI analysis for the Human Connectome Project Lifespan Ageing and Development studies Open Access

The Human Connectome Project (HCP) Lifespan studies extend the original HCP (collected in young adults) to characterise brain organisation and connectivity during the development (HCD, 5–21 years) (Somerville et al., 2018) and ageing (HCA, 36–100+ years) (Bookheimer et al., 2019) phases of life. Taken together, the studies acquired comprehensive neuroimaging datasets for more than 3,200 individuals to facilitate characterisation of the organisation and connectivity of the brain at various life stages, including childhood, puberty, middle age, menopause, later life and the “oldest old”. The Lifespan studies included multiple delay arterial spin labelling (ASL) MRI as an imaging modality, which permits, for the first time, perfusion measurements to be included in the HCP data releases.

ASL is a non-invasive imaging technique that uses magnetically labelled blood as an endogenous tracer; when acquired at multiple delays (as was the case here), it permits measurement of both cerebral perfusion (CBF) and arterial transit time (ATT). Cerebral perfusion is known to vary substantially during early life (Avants et al., 2015; Satterthwaite et al., 2014), to decline with age (Kiely et al., 2022; Leidhin et al., 2021), and to be affected by many of the diseases common in older age (Alosco et al., 2013; Binnewijzend et al., 2016; Hu et al., 2010; Melzer et al., 2011; Ortapamuk & Naldoken, 2006; Patel & Markus, 2011; Prohovnik et al., 2007; Shi et al., 2016). Though ASL has been acquired in previous population studies such as the UK Biobank (Miller et al., 2016), TILDA (Donoghue et al., 2018), GENFI (Mutsaerts et al., 2019), PREVENT (Mak et al., 2021), and Whitehall II Imaging Sub-study (Suri et al., 2019), the HCP Lifespan ASL dataset is unique in its size and quality, in particular due to the use of an unusually high spatial resolution multiple delay acquisition.

When the HCP was launched in 2010, the approach adopted was novel: to acquire a large and high-quality dataset, to maintain that quality during minimal pre-processing, and then release the outputs for others to analyse further. Many aspects of this “HCP-Style” to acquisition and analysis have been adopted by other studies (Glasser, Smith, et al., 2016); importantly, it has led to improvements in the reproducibility and transparency of research because all parties can access the underlying data (Elam et al., 2021). The objective of the HCP-ASL pipeline is to provide similarly high-quality pre-processed perfusion measures that others may use for downstream analysis, and many of the design decisions have been made with equivalence to existing HCP pipelines in mind in the hope of exploiting the data quality to its full extent, particularly with regard to cross-subject alignment of cortical areas.

A notable feature of existing HCP pipelines is the use of surface-based analysis to improve the alignment of cortical areas across individuals and reduce partial volume effects in group analyses (Coalson et al., 2018; Glasser et al., 2013; Glasser, Smith, et al., 2016), though exploiting the advantages of this paradigm requires accurate segmentations of the cerebral cortex. For the HCP studies, this was achieved via the use of high-resolution T1- and T2-weighted anatomical scans at 0.8 mm isotropic resolution or better, which may not be feasible for all studies (Glasser, Smith, et al., 2016). Conversely, the volumetric approach to group analysis relies on volumetric registration which is known to be inaccurate at aligning cortical areas across most of the cortical sheet and substantially worsens partial volume effects across tissues in group analyses (Coalson et al., 2018). In keeping with the existing HCP pipelines, the HCP-ASL pipeline produces perfusion estimates in a surface-based representation so that surface techniques may be used in downstream analysis. As few studies have yet explored surface representation of ASL-derived measures (examples include (Taso et al., 2021; Verclytte et al., 2015)), the HCP-ASL pipeline will enable the community to examine the utility of such an approach, in combination with other functional data also on the surface. The pipeline also produces volumetric outputs for compatibility with conventional analysis workflows.

This work details the processing steps performed by the HCP-ASL minimal processing pipeline and presents summary perfusion statistics calculated across the HCP Lifespan cohorts. There are numerous strategies that may be used for processing ASL data, particularly multi-delay data, and there remains debate as to which approach is best (Alsop et al., 2015; Fan et al., 2024; Lindner et al., 2023; Woods et al., 2024). We have not tried to systematically resolve these differences in this work, and instead present the pipeline as a best effort for the unique nature of the dataset. For an example analysis of perfusion measurements derived from the data, the reader is referred to Juttukonda et al. (2021).

The HCP Lifespan ASL data were processed with the HCP-ASL pipeline version 0.1.2, available at https://github.com/physimals/hcp-asl/releases/tag/v0.1.2. The following sections detail the data acquisition and the processing steps performed by the pipeline.

ASL data for the HCP Lifespan studies was acquired in 5.5 min of scanning on five 3T Siemens Prisma scanners at four sites (Harms et al., 2018). Data were collected for 1,306 subjects for HCD (655 female, ages 5–21 years) and 1,199 for HCA (681 female, ages 36–102 years); the distribution across sites is given in the supplementary material. Ethics approvals for the data collection were granted by the Washington University in St. Louis Institutional Review Board (HCA IRB ID #: 201603117, HCD IRB ID#: 201603135) and informed consent was obtained from all participants for being included in the study.

Pseudo-continuous ASL (PCASL) was used, with a 1.5 s label duration and multiple post-labelling delays (PLDs) of 0.2, 0.7, 1.2, 1.7, and 2.2 s, repeated 6, 6, 6, 10, and 15 times, respectively. A simultaneous multi-slice (SMS) 2D EPI acquisition with a multi-band factor of 6 was used without background suppression to achieve 2.5 mm isotropic resolution (2.27 mm slice thickness plus 10% gap), with 60 slices total, 59 ms readout time per slice, and 19 ms echo time (further details are given in Li et al. (2015)). Although 2D readouts have inherently lower signal-to-noise ratio than 3D readouts, they do not introduce substantial spatial blurring due to T₂ decay during the readout (Alsop et al., 2015; Woods et al., 2024). Such blurring would reduce the accuracy of registration between the ASL data and T1w anatomical image, negatively affecting the surface parameter maps that are the desired endpoint of the pipeline. Blurring can be reduced via the use of segmented 3D acquisitions, but such acquisitions are highly sensitive to motion. For this reason, after comparing 3D and 2D (with SMS) PCASL, HCP Lifespan adopted the latter (Harms et al., 2018).

To calibrate perfusion measurements into units of ml/100 g/min, two 2.5 mm isotropic PD-weighted M0 calibration images (TR > 8 s) were acquired at the end of the PCASL scan. In addition, a strong pre-saturation covering the imaging region was played out just before the PCASL module using three selective RF pulses with the spoiling gradients for each applied along two different directions (i.e., the first pre-saturation pulse was followed with spoiling gradients along the readout and slice directions; the second along the phase and slice directions; and the third along the readout and phase directions). For susceptibility distortion correction, two phase-encoding-reversed spin-echo images were acquired that were geometrically and distortion matched to the ASL data. T1- and T2-weighted anatomical and functional MRI images were acquired and pre-processed via the existing HCP minimal processing pipelines (Glasser et al., 2013, 2018, 2019; Harms et al., 2018; Robinson et al., 2018).

The ASL images contained several geometric and intensity artefacts that required correction before perfusion estimation could be performed. The use of a pre-saturation pulse in a multi-slice acquisition scheme contributed to slice-wise intensity variations in the inferior-superior direction in the ASL timeseries and calibration images, most notably between the last and first slices of adjacent bands, illustrated in Figure 1.

The perfusion signal in ASL is obtained by subtracting successive pairs of label and control images in the ASL timeseries. This means that ASL is sensitive to within-timeseries motion which can lead to spurious signal when mis-aligned images are subtracted, and it is, therefore, necessary to motion correct the data. Conventional methods for rigid-body motion correction such as FSL MCFLIRT (Jenkinson et al., 2002) may not be accurate when applied to banded data because the registration cost function may reward alignment of the intensity banding rather than anatomy, leading to anatomical misalignment, which will produce spurious signal after label–control subtraction. Furthermore, within a single label–control pair that have relative motion, even if it were possible to perfectly re-align anatomy there would remain artefactual intensity variations due to the bands occurring at different anatomical locations across the pair, resulting in imperfect elimination of the banding effects after subtraction. It was, therefore, decided to “de-band” the data before performing motion correction.

The first banding correction applied to the data accounted for the different relaxation times of each slice in the 2D readout, which gives rise to a saturation recovery effect in the signal. While this effect exists in a conventional non-SMS 2D ASL acquisition (Golay et al., 2004), with SMS acquisitions there are adjacent slices within brain tissue acquired in the first versus last SMS excitation, making the discontinuity more evident, and importantly, contributing to the prominent banding artefact that was observed. To remove the saturation recovery effect, a simplified model of the form $S(t)=1−e−tT1$ was fitted to the control images of the ASL timeseries using FSL FABBER, where T₁ is the longitudinal relaxation time of tissue and t the slice-specific interval from the pre-saturation module (which is the sum of the duration of the labelling module, the PLD, and the slice-specific delay in the readout module, accounting properly for the SMS nature of the acquisition). An example T1 map generated by FABBER is given in the supplementary material. The saturation recovery differences were then corrected (in both the control and label images) by normalising the intensities to the estimated value that they would have had if all the slices had been acquired simultaneously at the PLD.

During pipeline development, it was found that saturation recovery correction alone was unable to fully remove banding in both the calibration image and later PLDs of the ASL timeseries (illustrated in the Results section). Given the long TR of the calibration image (8 s), or long recovery times of the later PLDs (up to 3.7 s), theory predicts that longitudinal magnetisation will have largely recovered following those intervals and, therefore, the saturation recovery effect will be minimal. The presence of banding in the calibration images and later PLDs implied that there was a second banding mechanism (possibly a magnetisation transfer effect, addressed in the Discussion section), for which a correction was empirically derived using the calibration images of 80 subjects drawn equally from the HCA and HCD cohorts. After performing saturation recovery correction and masking the images to include grey and white matter only, a linear model for slice-wise mean intensity was fitted to each band and the coefficients averaged across subjects (see Supplementary Material for details).

Correction for gradient non-linearity distortion was performed using gradunwarp operating on the calibration image with the gradient coefficients file for the Prisma system (as used in the HCP PreFreeSurfer pipeline (Glasser et al., 2013)). Susceptibility distortion correction was performed using FSL TOPUP operating on a pair of phase-encode reversed spin-echo images that had undergone gradient distortion correction. Bias field (B1) estimation was performed using the spin echo-based approach (“SEBASED”, Glasser et al., 2016) operating on the calibration image that had undergone gradient distortion, susceptibility distortion, and banding correction.

Figure 2 shows how the aforementioned preprocessing corrections were derived and combined in the HCP-ASL pipeline to produce the fully corrected ASL and the calibration images on which perfusion estimation can be performed. Though two calibration images were acquired, only the first was used in the pipeline, though both were fully corrected. Some corrections were refined iteratively: for example, given that banding is expected to have a strong influence on the accuracy of motion estimation, each was performed twice on the assumption that an improvement in one should lead to an improvement in the other. Motion estimation was performed using MCFLIRT with the calibration image as the reference to yield a registration between the calibration image and ASL timeseries. In total, four subject-specific inputs and two generic inputs to the pipeline are required: for the individual, the phase encode reversed field maps, the calibration image, the T1w structural image, and the ASL timeseries, and generically, the Siemens Prisma gradient coefficients and the empirical banding correction coefficients.

In keeping with HCP fMRI pipelines, the HCP-ASL pipeline transforms ASL timeseries data into alignment with the structural images before analysis (Glasser et al., 2013), though retains the original 2.5 mm resolution (hereafter referred to as “ASL-gridded T1 space”). Multiple steps were used to obtain an accurate ASL to T1w registration. An initial registration was made using FreeSurfer’s boundary-based registration (BBR, bbregister) between the calibration image and the T1w image (Greve & Fischl, 2009; Jenkinson et al., 2002); this was combined with the motion correction matrices to obtain a registration between ASL and T1w. The calibration image has better grey–white tissue contrast than the un-subtracted ASL timeseries which facilitates accurate BBR. Then, perfusion estimation on the native ASL data with all banding and distortion corrections was performed to obtain a CBF map with increased grey–white tissue contrast, on which a second BBR was performed to the T1w image. The CBF map from this step was then discarded.

The transform between the ASL and calibration images and T1w alignment was merged with the motion correction, susceptibility distortion, and gradient distortion transforms and applied in a single step to the raw ASL and calibration images using the regtricks library to produce the fully transformed ASL and calibration images in ASL-gridded T1 space (Kirk, 2022). Jacobian intensity correction for distortion-induced signal pile-up was performed and cubic splines with pre-filtering were used to perform the interpolation (Unser et al., 1993). Minimising the number of interpolation operations aimed to reduce the introduction of additional partial volume effects that cannot be entirely removed with subsequent correction. SEABASED bias field estimation was re-run on the fully transformed calibration image to obtain bias field correction coefficients in ASL-gridded T1 space.

To intensity correct (banding and bias) the fully transformed ASL and calibration images, the calibration to T1w registration was applied to the empirical banding coefficients, and the ASL to T1w registration to the saturation recovery correction coefficients, to transform both banding corrections into ASL-gridded T1 space. Given the ASL and calibration image data were distortion corrected, no distortion correction was applied to the banding correction coefficients (it is only necessary to correct one or the other). The empirical banding coefficients were combined with the newly derived bias field estimates to intensity correct the calibration image, and both sets of coefficients were combined with the bias field estimates to intensity correct the ASL image.

In the presence of head motion, voxels travelling between neighbouring slices during the acquisition will receive differing intensity scaling during the banding corrections, which could lead to spurious signal after label–control subtraction. The general linear model (GLM) framework for motion-aware subtraction of banded and background-suppressed ASL data developed by Suzuki et al. (2019) was used in this work, with the small modification that it was not necessary to account for background suppression because it was not used during the image acquisition. The subtracted timeseries data were then used for perfusion estimation.

CBF and ATT estimation were performed using a variational Bayesian method via the oxford_asl script (Chappell et al., 2009, 2023). The aslrest Buxton model with CBF, ATT, and macrovascular components was used (Buxton et al., 1998; Chappell et al., 2011). A normal distribution prior with mean 1.3 s was used on ATT and an automatic relevancy determination (ARD) prior was used on macrovascular blood volume to remove this component from non-arterial voxels. Slice-timing correction was performed by adjusting the PLDs in each voxel by their slice-timing offset, and perfusion was converted from arbitrary units into ml/100 g/min using the mean signal value of CSF in the lateral ventricles from the calibration image (Pinto et al., 2020). Reference region calibration was used instead of a voxel-wise strategy due to the availability of a high-resolution calibration image, which allowed for accurate ventricular segmentation with minimal PVE, as well as mitigating the possibility of introducing (potentially uncorrected) banding artefacts from the calibration image if voxel-wise division were used. The tissue T1 values previously estimated for saturation recovery correction were passed to oxford_asl.

ASL is typically acquired with ~4 mm voxel sizes which gives rise to partial volume effects (PVE) caused by the coarse spatial resolution of the data in relation to the cerebral cortex, which has a mean thickness of around 2.5 mm. There is ongoing debate as to whether PVE correction (PVEc) should be routinely used in ASL studies. Though it is theoretically justified, particularly for studies covering a wide range of participant ages, there remain questions around the robustness and efficacy of currently available methods (Chappell et al., 2021). The ASL data acquired by the Lifespan HCP had a relatively high spatial resolution (2.5 mm isotropic), which not only reduced the PVE but also enabled more accurate PVEc to be performed. PVEc is performed in the voxel space using a spatial variational Bayes method implemented in oxford_asl (Chappell et al., 2011); the required partial volume estimates are obtained using Toblerone operating on the FreeSurfer-derived cortical surfaces and subcortical segmentations (Fischl, 2012; Kirk et al., 2020). For PVEc, a normal distribution prior with mean of 1.3 s was used for GM ATT and 1.6 s for WM ATT.

Both PVEc and non-PVEc perfusion estimates are generated so that the end user may decide which is most appropriate for their application. For non-PVEc, spatial regularisation is not used (deviating from the recommended settings for oxford_asl, which were otherwise used) to minimise the mixing of signal across the grey/white cortical boundary, at the cost of somewhat greater image noise. By contrast, for PVEc, spatial regularisation is an intrinsic feature of the spatial variational Bayes method used by oxford_asl.

Volumetric CBF and ATT maps from oxford_asl are produced in both ASL-gridded T1w space and MNI152 2 mm template space, via the existing FNIRT registration produced by the HCP PreFreeSurfer structural processing pipeline. To produce output parameter maps in greyordinates space (a combined surface and volumetric space represented in CIFTI format), the volumetric outputs of oxford_asl are projected onto the individual’s native cortical surface using the HCP ribbon-constrained method, registered with MSMAll multi-modal areal-feature-based registration and resampled to a common surface mesh, and finally smoothed with a 2 mm full-width half-maximum kernel using a surface-constrained method (Glasser et al., 2013, 2016; Robinson et al., 2014, 2018). The surface and MNI outputs (masked to consider subcortical structures only) were combined to produce the final CIFTI greyordinates files. Image-derived phenotypes (IDPs) were produced for each participant by computing the mean and standard deviation within each parcel of the HCP multi-modal parcellation (version 1.0) for both PVEc and non-PVEc variants of CBF and ATT.

Though the HCP-ASL pipeline does not perform any automated QC, it does produce intermediate outputs that enable researchers to perform visual checks. These are produced as a single Workbench scene file (a pre-specified visualisation) containing seven scenes, each of which is saved as a standalone image to facilitate quick inspection (although full exploration requires opening the scene in Workbench). The scenes illustrate (a) the bounding box of the ASL timeseries after motion correction within the acquisition field of view (voxels missing data due to their movement outside of the acquisition field-of-view for a single timepoint are excluded from perfusion estimation) and the quality of registration to the T1w scan, illustrated in Figure 4; (b) the raw ASL timeseries before and after all corrections are applied; (c) non-PVEc CBF displayed on the inflated cortical surface and in subcortical structures; (d) non-PVEc ATT displayed in the same manner; (e) PV estimates for GM and WM, and the ventricular CSF mask used for calibration; (f) a volumetric montage of non-PVEc CBF; and (g) a volumetric montage of non-PVEc ATT.

Figure 5 shows the effect of the two individual banding corrections (empirically derived and saturation recovery) individually and combined together on a slice-wise basis. The two banding mechanisms worked in opposing directions: at early PLDs, the saturation recovery effect was almost equal and opposite to the empirical banding effect, leading to a subtle overall correction, whereas for late PLDs where the saturation recovery effect was weaker, the empirical banding effect dominated.

Figure 6 shows the effect of the two banding corrections applied in isolation and together to the ASL timeseries from participant HCD0378150. The saturation recovery correction in isolation had minimal effect on the volumes of the ASL series with the longest (2.2 s) PLDs. This was consistent with theoretical prediction: at longer PLDs, the relative differences in slice-timing decrease, which will reduce the prominence of the saturation recovery effect, as was observed. Conversely, the empirical correction in isolation was unable to fully remove banding on the shortest (0.2 s) PLDs, for which the saturation recovery effect is substantial. For all PLDs, the combination of the two separate corrections was better able to remove banding than in isolation.

Figure 7 shows the group average CBF maps for the HCA cohort, both with and without PVEc. In these and the following figures, the HCP multi-modal cortical parcellation (version 1.0) has been overlaid in black on the surface (Glasser, Coalson, et al., 2016). Areas of differential CBF (colour transitions) were observed to follow areal boundaries on the cortical surface. PVEc notably increased CBF, particularly in cortical tissue.

Figure 8 shows the group average ATT maps for the HCA cohort, both with and without PVEc. Along with the volumetric representation given in Figure 11, ATT showed correspondence with known vascular anatomy, namely that the cores of the vascular territories (middle, anterior, and posterior cerebral arteries; MCA, ACA, PCA) were perfused before watershed regions that lie along the margins of these vascular territories; and anterior circulation (supplied by the MCA and ACA) was perfused before posterior circulation (supplied by the PCA). PVEc had a negligible effect on ATT.

Figures 9 and 10 show group average CBF and ATT maps for the HCD cohort, both with and without PVEc. The ATT maps showed longer transit times for posterior and watershed regions.

Figures 11 and 12 show group average maps for non-PVEc CBF and ATT in each cohort, shown in a volumetric representation. This is a conventional representation of the same data shown in the preceding figures.

Figure 13 shows distributions of each subject’s mean GM CBF and ATT before and after PVEc, grouped by cohort (the corresponding distributions for WM are given in the Supplementary Material).

Numerical values for each distribution mean are given in Table 1. Between cohorts, statistically significant differences were observed (all comparisons with p < 0.05). On the non-PVEc data, the HCA cohort had lower CBF than HCD (45 vs. 68 ml/100 g/min) and longer ATT (1.46 vs. 1.28 s).

Table 2 shows that PVEc significantly increased GM CBF in both cohorts: from 45 to 61 ml/100 g/min in HCA and from 68 to 85 ml/100 g/min in HCD. For both cohorts, PVEc led to a small but statistically significant decrease in GM ATT on the order of 0.01 s.

Figure 14 shows the effect of PVEc on GM CBF and GM ATT measurements, separated into subcortical and cortical greyordinates, for a single subject. For cortical greyordinates, the average increase in CBF following PVEc was around 15 ml/100 g/min, whereas for subcortical greyordinates, it was around 5 ml/100 g/min. For ATT, both cortical and subcortical greyordinates showed negligible average increases.

Figure 15 shows CBF and ATT maps for a single participant of the HCD cohort. Although the individual subject maps were noisier than the group average, it was possible to observe similarities with the group average maps; for example, elongated ATT in posterior and superior regions of cortex.

Figure 16 shows the same HCD individual as Figure 15, where CBF and ATT have been averaged within the parcels of the MCP multi-modal parcellation. This reduced the dimensionality of the data dramatically, enabling the data in Figure 12 to be represented in 379 parcels (whereas there are 91,282 greyordinates in Fig. 11).

Figure 17 shows volumetric CBF and ATT maps for a single subject of the HCD cohort. Though these appear to be noisy compared with typical ASL-derived maps, it should be noted they are at unusually high resolution and have not had any form of smoothing applied. Partial volume corrected variants of these maps are given in the supplementary material.

The HCP Lifespan ASL dataset and accompanying pipeline have yielded a dataset that is larger than existing ASL studies. The use of a high spatial resolution multiple delay acquisition has enabled CBF and ATT maps to be produced in both volumetric and surface representations. The group average surface-based ASL measures, which were generated using areal-feature-based surface registration with minimal smoothing, represent the most high-resolution maps of brain perfusion currently available, and will enable the wider exploration of surface-based analysis techniques for ASL. The utility of these surface maps has yet to be fully explored for ASL imaging, as it represents a significant departure for a modality that has nearly always been analysed in a volumetric manner.

CBF on the cortical surface was often observed to follow cortical areal boundaries (see Figs. 7 and 9), which suggests that cerebral perfusion may be further regulated within vascular territories. Given that prior results from HCP fMRI data have established functional and/or connective relationships between cortical areas, it would be interesting to investigate whether perfusion is a further component of these relationships. Furthermore, given perfusion is a component of the BOLD effect, greater understanding of this parameter in the cortex will be of utility to fMRI studies. Comparing between cohorts, the HCD cohort had higher CBF than HCA (a difference of 24.7 ml/100 g/min), though the difference was slightly reduced by PVEc (23.4 ml/100 g/min), which is consistent with age-related atrophy leading to increased PVE presenting a confounding factor for perfusion measurement with ASL on elderly cohorts.

The ATT maps, both surface and volumetric, evidenced longer transit times in posterior portions than the anterior portions, which is expected because the vascular calibre of the common and internal carotid arteries that feed the ACA and MCA is larger than that of the vertebral arteries that feed PCA through the basilar artery, resulting in faster flow at a given blood pressure. Further, the HCA cohort showed higher GM ATT than HCD by ~0.2 s, which is consistent with previously reported findings that ATT increases with age (though the difference was smaller in WM, as detailed in supplementary material). The group mean volumetric ATT maps showed some small differences in GM/WM contrast between the two cohorts. One explanation for this could be reduced contrast in core/watershed ATT in the HCD cohort, making GM/WM contrast more prominent, whereas in the HCA cohort, core/watershed contrast dominates, possibly due to vascular disease.

Due to the high spatial resolution of the acquisition, the volumetric maps produced by the pipeline appeared noisier than ASL acquired at typical resolution (e.g., 3.4 x 3.4 x 5 mm). This was particularly the case in WM (shown, for example, in Fig. 17). Not only does WM have intrinsically lower perfusion (and thus SNR) than GM, which poses a fundamental challenge for ASL, but also the acquisition itself was optimised for the cortex (high spatial resolution was used to facilitate accurate registration to the cortical surface, and the sensitivity profile for the multi-channel coil used preferentially increased SNR in the cortex). Under the Bayesian framework used by the FSL BASIL toolbox for perfusion quantification, parameter values will revert to their prior mean values when the data are uninformative. For WM ATT voxels, this is a likely explanation for why values tend to be near the prior mean of 1.6 s, although this would be expected to only have a small effect on the estimated perfusion value. Ultimately the pipeline presented here has not been optimised for WM voxels and more precise results might be achieved by averaging over regions in WM.

Nevertheless, recent work investigating high-resolution ASL has found it to be practical and feasible for perfusion measurement in the cortex and subcortex, despite the appearance of noise, because the reduction of inherent partial volume effects improves the localisation of perfusion without impacting sensitivity (Kashyap et al., 2024). A conventional mitigation for noise would be to perform spatial smoothing (similar to reducing spatial resolution), at the loss of spatial precision. The alternative strategy adopted here is to average perfusion measurements within neuroanatomically well-defined subdivisions to make efficient use of available SNR. This should increase statistical sensitivity for effects, while reducing dimensionality (thereby reducing the penalty for correction of multiple comparisons) without the need for indiscriminate spatial smoothing. These imaging-derived phenotypes should facilitate novel statistical analyses such as relating CBF and ATT to behaviour, cognition, genetics, or clinical conditions.

In both cohorts, PVEc was observed to lead to statistically significant increases in CBF estimates. This is despite the relatively high spatial resolution of the ASL data (2.5 mm isotropic), which reduces but does not eliminate the presence of PVE in the data. Indeed, Figure 14 shows that PVEc led to greater increases in cortical than subcortical CBF, which is consistent with there being more PVE present in cortical voxels. Recent work has noted the fundamental trade-off inherent to PVEc, while acknowledging that both PVEc and non-PVEc data have value depending on the research question (Chappell et al., 2021). Without PVEc, differences in GM PVE between individuals due to differences in anatomy are encoded into the perfusion values, thus becoming a source of between-individual variability. Though PVEc mitigates this, it also increases the complexity of the model by introducing tissue-specific parameters which could result in over-fitting, introducing a different source of between-individual variability.

The brainstem was observed to have relatively high CBF after PVEc. This could be due to inaccuracy of PV estimates when segmenting this region of the brain (which would directly impact PVEc), the presence of a stronger bias field in the middle of the brain, or macrovascular contamination from nearby arteries. At any rate, it is important to note that perfusion measurement via ASL in WM is inherently difficult due to the low SNR associated with this tissue.

The novel nature of the multi-band acquisition and the high demands for precision in spatial localisation in the HCP-Style approach to data acquisition and analysis (Glasser, Smith, et al., 2016) necessitated careful pre-processing of the ASL and calibration data to bring it into T1w alignment space at the original voxel resolution with properly normalised image intensities. The pipeline that has been implemented represents a best attempt at preprocessing given the unique nature of the acquisition, though it has not been compared against existing ASL pipelines. Recognising that users may wish to implement their own alternative kinetic modelling to the variational Bayesian method used in this work, the HCP-ASL processing pipeline also provides fully corrected ASL and calibration volumetric data, which means that a bespoke perfusion estimation may be performed without needing to re-create the prior correction steps.

Addressing the SMS banding artefact was a particular pre-processing challenge. The banding remaining after correction for saturation recovery is suspected to be caused by magnetisation transfer effects, namely that the acquisition of one slice has spill-over effects on the adjacent slice acquired next. In the derivation of the empirical correction used in this work, it was found that the coefficients differed slightly between samples of the HCA and HCD cohorts (illustrated in the Supplementary Material). This could be due to age-related differences in grey matter and white matter volumes. Ultimately it was decided to use a single set of correction coefficients derived over an equal-sized sample of both cohorts for simplicity, but in the future, a personalised correction technique to remove the need for a representative population might be feasible and would allow the pipeline to translate to unseen datasets with ease. A personalised approach could also be beneficial for application to subjects where tissue T1 may be modified by pathology, as this mechanism at least partially explains the banding artefact.

Across both cohorts, lower CBF was observed in the most inferior cortical regions such as the inferior temporal cortex and orbito-frontal cortex. It is likely that this was caused by uncorrectable susceptibility-induced signal loss. Though distortion correction can relocate signal to the correct position within the brain, it cannot correct the underlying loss of signal, which could only be solved by the use of a spin-echo acquisition.

An important limitation of the pipeline in its current form is that it does not implement any form of automated QC (instead providing some visual outputs that can be used as part of manual QC). In the context of the HCP ASL dataset, this means researchers must make their own judgement on which subjects to exclude from analysis. As for application of the pipeline to more heterogeneous or possibly pathological clinical data, this cannot be assumed to work without modification (particularly as some of the corrective steps implemented are acquisition specific). QC of motion estimation and correction would be a particularly important consideration for application of the pipeline to very young or elderly clinical subjects for which substantial motion is to be expected.

An area in which the pipeline could be improved is image denoising, potentially using independent component analysis (Carone et al., 2019). Preliminary work with this approach was unable to reliably identify noise components, so it was not adopted for the final pipeline. This could be due to the relatively limited number of timepoints in ASL data relative to fMRI, which reduces the performance of ICA (i.e., ICA’s performance is markedly improved in high spatial and temporal resolution fMRI relative to legacy low spatial and temporal resolution fMRI).

A further area of future development would be to adopt a fully surface-based approach to perfusion estimation. Since the BASIL toolbox operates in volume space, currently the pipeline produces all perfusion estimates in voxel space and then projects these onto the cortical surface. A more direct, and potentially more advantageous, approach would directly estimate perfusion on the cortical surface after projecting the ASL timeseries from volumetric data (Kirk, 2021). This would allow spatial regularisation (which is beneficial given the inherently low SNR of ASL data) to be applied in a manner that respects the existing HCP principle of constrained smoothing within neuroanatomically well-defined tissues or brain areas. Indeed, one could in theory directly model the ASL parameters within brain areas, a technique that has previously been shown to be beneficial for task-based fMRI GLM modelling (Glasser, Coalson, et al., 2016).

The HCP Lifespan ASL dataset is a large and high-quality resource that will enable perfusion to be studied in unprecedented detail during the development and ageing phases of life. The HCP-ASL pipeline has been developed as a best attempt to implement an HCP-Style data analysis while accounting for the specific nature of the acquisition with state-of-the-art correction techniques. Perfusion estimation has been performed on the pre-processed data using an established variational Bayesian method, though investigators can run their own perfusion estimation on the pre-processed data if they wish. Given ongoing debate on the utility of partial volume correction, particularly where ageing and pathology are concerned, both corrected and non-corrected perfusion estimates have been produced to provide investigators maximum flexibility. A preliminary group analysis revealed expected differences in haemodynamics between the development and ageing cohorts, specifically, reduced CBF and elongated ATT with increasing age.

The pipeline code is found at https://github.com/physimals/hcp-asl/releases/tag/v0.1.2 and the pre-processed data will be available at the NIH data archive (NDA): https://nda.nih.gov

T.F.K., F.A.K.M., and J.T. developed the pipeline code and prepared the manuscript. M.S.C. contributed to pipeline code. D.C., X.L., and Y.S. contributed to methodology and reviewed the manuscript. T.S.C. contributed to pipeline code and reviewed the manuscript. M.P.H. and M.F.G. contributed to methodology, performed data analysis, and reviewed the manuscript. M.F.G. also secured funding. M.A.C. supervised pipeline development, contributed to methodology, secured funding, and reviewed the manuscript.

The authors declare no competing interests.

J.T. and M.A.C. have received support from the Engineering and Physical Sciences Research Council UK (EP/P012361/1, EP/P012361/2). F.A.K.M. and M.A.C. were supported by the Beacon of Excellence in Precision Imaging, University of Nottingham. Y.S. is supported by the Royal Academy of Engineering under the Research Fellowship scheme (RF/201920/19/236). The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z). M.P.H. is supported by grants: U01MH109589 and U01AG052564. M.F.G. receives support from grants: R24MH108315, R24MH122820, and U19AG073585. T.F.K. was funded by the Bellhouse scholarship at Magdalen College, Oxford.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00444.