MVP-VSASL: measuring MicroVascular Pulsatility using Velocity-Selective Arterial Spin Labeling
==============================================================================================

* Conan Chen
* Ryan A. Barnes
* Katherine J. Bangen
* Fei Han
* Josef Pfeuffer
* Eric C. Wong
* Thomas T. Liu
* Divya S. Bolar

## 1 Abstract

**Purpose** By leveraging the small-vessel specificity of velocity-selective arterial spin labeling (VSASL), we present a novel technique for measuring cerebral MicroVascular Pulsatility named MVP-VSASL.

**Theory and Methods** We present a theoretical model relating the pulsatile, cerebral blood flow-driven VSASL signal to the microvascular pulsatility index (PI), a widely used metric for quantifying cardiac-dependent fluctuations. The model describes the dependence of PI on bolus duration *τ* (an adjustable VSASL sequence parameter) and provides guidance for selecting a value of *τ* that maximizes the SNR of the PI measurement. The model predictions were assessed in humans using data acquired with retrospectively cardiac-gated VSASL sequences over a broad range of *τ* values. In vivo measurements were also used to demonstrate the feasibility of whole-brain voxel-wise PI mapping, assess intrasession repeatability of the PI measurement, and illustrate the potential of this method to explore an association with age.

**Results** The theoretical model showed excellent agreement to the empirical data in a gray matter region of interest (average R2 value of 0.898 *±* 0.107 across six subjects). We further showed excellent intrasession repeatability of the pulsatility measurement (ICC = 0.960, *p<* 0.001) and the potential to characterize associations with age (*r* = 0.554, *p* = 0.021).

**Conclusion** We have introduced a novel, VSASL-based cerebral microvascular pulsatility technique, which may facilitate investigation of cognitive disorders where damage to the microvasculature has been implicated.

Keywords
*   velocity selective arterial spin labeling
*   VSASL
*   ASL
*   pulsatility
*   microvascular

## 2 Introduction

Pulsatile blood flow driven by the cardiac cycle has been increasingly linked to structural damage in the cerebral microvasculature1–4 and cognitive disorders such as mild cognitive impairment, Alzheimer’s disease, and other dementias1–7. In young, healthy subjects, this flow pulsatility is dampened by compliant arteries as the pulse wave travels toward the distal microvasculature and brain parenchyma. However, if this dampening is insufficient, for example, due to pathologic changes in the vasculature, then excess pulsatile energy will reach the smaller vessels, where the resulting exposure and subsequent damage is thought to be a precursor to the aforementioned cognitive disorders3,8. There has been much work assessing the pulsatility4,8–12 and vessel wall compliance10,13–17 at the large cerebral arteries. However, techniques measuring pulsatility in the microvasculature itself, the site where damage is thought to primarily occur, have received less attention. Such techniques will be crucial to clarifying the exact mechanistic links between microvascular pulsatility, tissue damage, and eventual neurodegeneration, which are not yet fully understood. Furthermore, since neurovascular factors such as flow pulsatility are viewed as early and modifiable risk factors contributing to these disorders, the ability to measure these biomarkers may be important for developing strategies for early detection and intervention5–7,18.

Historically, microvascular pulsatility has been challenging to measure due to the very small size of microvascular vessels and their relatively slow flow. Recent approaches using phase contrast MRI have leveraged an ultra-high 7T field strength to measure pulsatility in small perforating arteries19–22. Additional advances to data acquisition (higher temporal resolution20, dual velocity encoding22) and post-processing (automated vessel detection20) have improved the usability and robustness of the technique even further. However, this approach remains challenging at the lower field strengths common on clinical scanners, which limits its potential for clinical translation. As an alternative approach, we utilize a technique called velocity-selective arterial spin labeling (VSASL)23 to measure microvascular pulsatility. VSASL can be performed on clinical 3T scanners with a simple whole-brain scan prescription and has the potential to generate cerebral microvascular pulsatility maps on a voxel-wise basis.

VSASL is a variant of arterial spin labeling (ASL), a family of methods used to measure perfusion by magnetically labeling a bolus of arterial blood, allowing the bolus to flow into the microvasculature or tissue of interest, and then acquiring images sensitized to the blood flow24. While the traditional brain ASL variants generate the magnetic label within the feeding extracranial carotid and vertebral arteries, VSASL is unique in its ability to generate the label in smaller, more distal arteries23. In VSASL, the user-specified sequence parameter of cutoff velocity (*v*cut) determines the location along the vasculature where the leading and trailing edges of the bolus are defined, as illustrated in Figure 1. Using a typical setting of *v*cut = 2 cm/s23, the sequence will define the boundary of the labeled bolus in small arterioles with vessel diameters of about 50 µm25, which lies within the range of vessel diameters in the microvasculature. In the standard VSASL sequence, a velocity-selective (VS) label/control module (LCM) will first label blood flowing faster than *v*cut, thus defining the leading edge of the labeled bolus (Figure 1A). This is followed by a delay corresponding to the bolus duration *τ* (also a user-specified sequence parameter)23, during which the bolus flows distally toward its target tissue while decelerating below *v*cut in the process (Figure 1B). A second VS module called the vascular crushing module (VCM) is then applied to saturate remaining labeled blood still flowing faster than *v*cut, thus defining the trailing edge of the labeled bolus (Figure 1C). The labeled bolus signal is proportional to the volume of blood that flows across the *v*cut boundary in the time between LCM and VCM, thus making the labeled bolus signal (i.e. VSASL signal) sensitive to flow rate and its variation across the cardiac cycle.

![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F1.medium.gif)

[Figure 1:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F1)

Figure 1: 
Illustration of the VSASL labeling process. The top most diagram shows the VSASL pulse sequence, with events (A)-(C) corresponding to the sub-figures below. Below the sequence diagram is an example CBF(*t*) waveform, with the area shaded in red representing the VSASL signal *S*. Events (A)-(C) are described as follows: (A) Blood signal flowing faster than the cutoff velocity *v*cut is labeled by the LCM, which defines the leading edge of the bolus. Since *v*cut is specified to the velocity in small arterioles (using a typical *v*cut = 2 cm/s for our study), this leading edge is defined in the microvascular regime. (B) After labeling, there is a waiting period of *τ* seconds, where *τ* is the user-specified bolus duration. During this period, blood continues to flow downstream, past the *v*cut boundary and toward the capillaries and tissue. (C) Blood signal flowing faster than *v*cut is crushed by the VCM. This action defines the trailing edge of the bolus, and thus fully defines the temporal duration of the bolus as *τ*, with the image then acquired shortly thereafter.

Previous work by Franklin et al. using a single VS labeling module alone (the LCM) measured fluctuations of up to 36% in the amount of arterial label generated26. This measurement was made in an arterial ROI and was weighted by macrovascular blood volume and flow, since the LCM was applied without an accompanying VCM. Nevertheless, the results suggested the potential of VSASL to measure effects of cardiac pulsations. In this study, we use a standard VSASL sequence design that includes both the LCM and VCM23 to specifically achieve microvascular blood flow weighting. By then using retrospective cardiac gating and leveraging the microvascular specificity of VSASL, we realize the potential of measuring blood flow pulsatility in the microvasculature. We dub this technique MicroVascular Pulsatility using VSASL (MVP-VSASL).

We first present theory relating the pulsatile VSASL signal to the underlying pulsatile blood flow, and then examine the pulsatility of the VSASL signal using the pulsatility index (PI) metric. We derive a theoretical model describing the dependence of PI on bolus duration *τ*, and use the model to determine a theoretically optimal *τ* value for maximizing the SNR of the PI measurement. Using experimental VSASL data acquired in human subjects with a broad range of ages and heart rates, we then validate the predicted *τ*-dependence of PI and assess the intrasession test-retest repeatability of the PI measurement. In addition, we examine the association between pulsatility and age and demonstrate the feasibility of a voxel-wise pulsatility mapping approach that may facilitate regional pulsatility measurements in future applications.

## 3 Theory

### 3.1 Pulsatility of VSASL Signal

The control-label subtraction signal (denoted *S*) from a VSASL scan represents the signal of the labeled blood being delivered to the microvasculature. This signal *S* is typically modeled as being proportional to CBF · *τ* · exp (−*τ*/*T*1*b*), where CBF reflects uniform (non-pulsatile) cerebral blood flow, *τ* is bolus duration, and exp (− *τ/T*1*b*) is the *T*1 decay weighting factor (where *T*1*b* is the *T*1 of blood)23.

In the case of time-varying, pulsatile blood flow, the product CBF · *τ* becomes an integration of CBF(*t*) over the duration of the bolus, and the continuous-time VSASL signal *S*(*t*) can be described as: ![Formula][1]</img>  where *** denotes convolution and rect(*t/τ*) is a rectangular function of width *τ* representing the bolus. To examine cardiac-driven pulsatility, we model CBF(*t*) with a 2nd-order Fourier model, which has previously been used for cardiac-gated measurements in ASL13,17,27. This model of CBF(*t*) is: ![Formula][2]</img>  where *b* (= CBF), *b*cos,1, *b*sin,1, *b*cos,2 and *b*sin,2 are the Fourier coefficients, and r represents the mean cardiac period (RR interval, or inverse heart rate). Evaluating the convolution in Equation 1 yields: ![Formula][3]</img>  a smoothed version of the CBF(*t*) driving function with the coefficients modified by an order-dependent sinc(*kτ/*Γ) term, where ![Graphic][4]</img>.

The pulsatility of a given waveform *S* can be quantified via the pulsatility index (PI)9,10,19,28, which is given by: ![Formula][5]</img>  where *S*max is the maximum of *S*(*t*), *S*min is the minimum and *S*mean is the mean.

By applying Equation 4 to the VSASL signal in Equation 3, a model for PI as a function of *τ* can be derived, with the full expression of PI(*τ*) given in Equation S.10 of Section S1. Although useful to examine, Equation S.10 is not an ideal platform for further analysis as it involves evaluating extrema of a second-order Fourier function (which do not simplify further) and is parametrized by CBF coefficients (whose values are generally not known *a priori*). However, in physiological flow waveforms, the fundamental (1st-order) frequency is typically the largest harmonic in the power spectrum29. By assuming ![Graphic][6]</img> and neglecting the 2nd-order terms of *S*(*t*) (see Section S1 for mathematical justification of this approximation), Equation 4 can then be applied to the VSASL signal in Equation 3 to yield: ![Formula][7]</img>  where ![Graphic][8]</img> is a lumped fitting parameter. This simple expression describes a sinc-shaped dependence of PI on the ratio of the bolus duration *τ* (an adjustable VSASL scan parameter) to cardiac period Γ. Note that this form predicts that PI will vanish when *τ* = Γ, a feature that will be examined later with in vivo data.

### 3.2 Optimal Choice of Tau

We also determine the value of *τ* that maximizes the SNR of PI. In Section S3, we show that the SNR of the PI measurement is approximately proportional to the product of the SNR of the VSASL signal (∝*τ* · exp (−*τ*/*T*1*b*)) and the magnitude of PI (Equation 5): ![Formula][9]</img>  

This expression is maximized at *τ*opt: ![Formula][10]</img>  yielding an optimal value slightly below Γ*/*2. Figure 2A is a graphical demonstration of optimizing Equation 6 to obtain Equation 7. For an example cardiac period of Γ= 1, the optimal tau is computed as *τ*opt = 0.437 s and indicated on the graph, and Figure 2B then plots Equation 7 to show *τ*opt over a range of cardiac periods. In Figure 2A, we also note that the SNR(*τ*) curve is relatively flat around *τ* = *τ*opt. For example, the nearby point of SNR(Γ*/*2) = 0.219 (indicated by the black circle) is within 2% of SNR(*τ*opt) = 0.223. Figure 2C shows that the discrepancy between SNR(*τ*opt) and SNR(Γ*/*2) is small for a range of Γ (*<* 3% for typical Γ in [0.6, 1.2] s), suggesting that nearby *τ* values such as Γ*/*2 can also roughly optimize SNR with relatively little tradeoff.

![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F2.medium.gif)

[Figure 2:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F2)

Figure 2: 
(A) Graphical representation of determining optimal tau *τ*opt following Equations 6 and 7, for an example cardiac period Γ = 1 s. The red dashed line indicates the VSASL signal SNR model *τ* · exp (− *τ /T*1*b*). The blue dashed line indicates the PI sinc model of Equation 5. The purple solid line indicates the combined PI SNR model (Equation 6), which is optimized at *τ*opt (computed via Equation 7), as indicated by the green circle. (B) A plot of *τ*opt vs Γ per Equation 7. Several points are highlighted to indicate *τ*opt for values of Γ in a typical physiological range. (C) A plot of SNR(*τ*opt) and SNR(Γ*/*2). The discrepancy between SNR(*τ*opt) and SNR(Γ*/*2) is shown to be fairly small (*<* 3% for r in [0.6, 1.2] s), suggesting that values near *τ*opt, such as Γ*/*2, could also yield near-optimal SNR.

A supplementary analysis on optimal *τ* is provided in Section S3, which presents the theory behind Equation 6 and the approximations used in its derivation. As additional support, the theoretical SNR model of Equation 6 is also compared with a reference SNR curve derived from Monte Carlo simulations. For those simulations, measurement noise was first added to simulated VSASL signals, distributions of PI measurements were computed, and then SNR was calculated from those PI distributions, with the process repeated over a range of *τ* values. As shown in the Section S3.3 results, Equation 6 shows excellent agreement with the simulated SNR curve, and *τ*opt indeed yields close to the maximum SNR of the PI measurement.

## 4 Methods

### 4.1 Overview

A primary cohort of 7 healthy subjects (3 females and 4 males, aged 38.1 *±* 14.9 years), plus a secondary cohort of 10 relatively older healthy subjects (9 females and 1 male, aged 65.9 *±* 10.0 years), were enrolled in this study. The study was approved by the UCSD Institutional Review Board (IRB), and informed consent was obtained from all participants. The experiment consisted of a series of MRI scans, all acquired on a 3T MAGNETOM Prisma (Siemens, Erlangen, Germany) using a 64ch head/neck receiver coil. During the MRI scans, photoplethysmography (PPG) data were collected from the subject’s index finger using the built-in system sensor. These PPG data were used for retrospective cardiac gating.

### 4.2 MRI Acquisition

The scanning protocol for the primary cohort consisted of a localizer, a *T*1-weighted structural scan, and then a series of six VSASL scans: two at *τ* = 500 ms (to evaluate testretest repeatability near *τ*opt per Equation 7), and then one each at progressively longer values (750/1000/1250/1500 ms) for a *τ* -stepping experiment to test for the predicted sinc-dependence described in Equation 5. All primary cohort subjects completed the full protocol, with a few exceptions noted in Table 2. The secondary cohort of 10 subjects underwent a short scanning protocol consisting of a *T*1-weighted structural scan and a single VSASL scan at *τ* = 500 ms, with all subjects completing this short protocol. *T*1-weighted structural scan parameters are provided in Section S4.

View this table:
[Table 1:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/T1)

Table 1: 
Table of VSASL scan parameters, indicating bolus duration *τ*, background suppression inversion times (BGS TIs), number of label/control pairs, repetition time *TR*, and total scan duration. All scans used *T**sat* = 1500 ms and *PLD* = 100 ms. Also, all VSASL scans contained one *M* (proton density-weighted) repetition at the beginning of the sequence with *TR* = 5000 ms. The BGS timings were configured via a MATLAB optimization routine to target complete suppression (0% signal) of GM and CSF tissue signals at 50 ms before the start of readout. This routine also considered *T*2-decay during the label-control module (LCM) and vascular crushing module (VCM) assuming an effective echo time (eTE) of 22 ms, and *T*1,*GM* = 1300 ms, *T*2,*GM* = 100 ms, *T*1,*CSF* = 4000 ms and *T*2,*CSF* = 2000 ms.

View this table:
[Table 2:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/T2)

Table 2: 
Table of subject demographics and acquired data. Test-Retest refers to subjects who acquired two repeats of the *τ* = 500 ms scan. *τ* -stepping refers to subjects who completed the *τ* -stepping protocol consisting of scans at *τ* = 500*/*750*/*1000*/*1250*/*1500 ms. Only one *τ* = 500 ms scan was acquired on Subject 1 due to time constraints. This was rectified with longer scheduled scan sessions for all subsequent subjects. Both *τ* = 500 ms scans were acquired on Subject 3, but the first was excluded from analysis due to head motion. Subject 7 faced time constraints from late arrival, so only the two scans at *τ* = 500 ms were collected.

Note: The genders were redacted and the ages were replaced by a range for the medRxiv submission.

The VSASL pulse sequence consisted of a velocity-selective preparation module and an accelerated 3D gradient- and spin-echo (GRASE) readout. VSASL scans were configured for the five different bolus durations of *τ* = 500/750/1000/1250/1500 ms, with repetition time (TR) adjusted to accommodate the *τ* setting of a given scan. Other sequence timing parameters *T**sat* and post-labeling delay (PLD) were kept fixed across all scans at *T**sat* = 1500 ms and PLD = 100 ms, respectively. The PLD was kept to the minimum necessary to accommodate two spectrally-selective fat-saturation modules and two inferior saturation modules (to minimize CSF inflow effects) immediately prior to readout. Two background suppression (BGS) pulses were inserted for each scan between the LCM and the VCM. Additional sequence timing details are provided in Table 1.

An eight-segment ![Graphic][11]</img> insensitive rotation (BIR-8) train23,30 was used for velocity-selective saturation for both the LCM and the VCM, with velocity weighting in the inferior-superior direction. The BIR-8 train was configured with a cutoff velocity of *v*cut = 2 cm/s, with cutoff velocity defined as the first zero-crossing of the laminar flow response of the label condition23.

The 3D GRASE readout was configured as: matrix size = 56 *×* 56 *×* 24; voxel size = 4 *×* 4 *×* 6 mm3; FOV = 224 *×* 224 *×* 144 mm3; 10% Slice Oversampling; 6/8 Phase Partial Fourier; 2 *×* 2 GRAPPA Acceleration; Turbo Factor = 13; EPI Factor = 21; 1 segment (single-shot acquisition); flip angle = 120 degrees; TE = 12.62 ms.

### 4.3 Data Pre-Processing

For each subject, all VSASL scans were co-registered and motion-corrected using AFNI’s 3dvolreg command31. The *T*1-weighted structural scans were processed using FSL’s fsl_anat command32,33, producing a whole-brain mask and a partial volume map of gray matter (GM). The GM partial volume map was thresholded at 0.8 (80%) to produce a GM mask. The brain and GM masks were nearest neighbor-resampled to the VSASL scan resolution using AFNI’s 3d_resample31.

To minimize CSF contamination23, the VSASL labels and controls were pair-wise subtracted in MATLAB R2021A (Natick, MA) to form a perfusion-weighted time series, and then the median absolute deviation (MAD) was computed voxel-wise across the time dimension. To identify high-variance voxels contaminated by CSF, a whole-brain MAD threshold was defined at the 80th-percentile, along with a set of slice-specific MAD thresholds defined at the 80th-percentile within every brain slice. Voxels with MAD values above either the whole-brain or slice-specific thresholds were excluded from the GM mask, and the resulting mask served as the final GM ROI for the subsequent analysis.

Data were denoised using the iterative denoising method described in Power et al.34, with a few modifications: incorporating 2nd-order Fourier regressors (to preserve desired cardiac-driven fluctuations) and censoring volumes with outlier Γ values (e.g. due to finger motion). Supplementary details on denoising are available in Section S5.

### 4.4 Data Analysis

For each VSASL scan, the denoised data were spatially averaged within the GM ROI to produce time series for control and label measurements, respectively. A cardiac phase ϕ was assigned to every data point via PPG-based retrospective gating15,27. The control and label time series were separately fit to 2nd-order Fourier models based on prior work27. The VSASL signal was then obtained by subtracting the control and label Fourier coefficients to yield: ![Formula][12]</img>  where ![Graphic][13]</img> and ![Graphic][14]</img> are the resultant Fourier coefficients of the signal *S*(*ϕ*). After obtaining *S*(*ϕ*), PI was quantified using Equation 4. The stability of the PI measurements was then assessed using a residuals permutation approach35. This procedure was repeated over 1000 iterations, and the resulting distribution of PI values was used to compute 95% confidence intervals. Additional details on the *S*(*ϕ*) computation and residuals permutation approach are provided in Section S2.

The measured PI values from the *τ* -stepping VSASL scans were fit to a version of Equation 5 modified to account for heart rate variability: ![Formula][15]</img>  where *К* (*τ*) is the average of |sinc ![Graphic][16]</img>| functions evaluated over all the Γ*i*’s observed during the scanning session: ![Formula][17]</img>  and *N*Γ is the number of cardiac periods Γ observed. This effectively weights the model based on the distribution of Γ*i*’s observed during the scanning session. The model fits were evaluated using R2.

Using subjects with two runs of the *τ* = 500 ms scan, the test-retest repeatability of the GM ROI PI measurement at *τ* = 500 ms was assessed by using the Pearson correlation coefficient and intraclass correlation coefficient (ICC).

The association of PI with age across subjects was also assessed. Because PI depends on the ratio of *τ /*Γ, the comparison across subjects and conditions is facilitated by evaluating PI at a common reference ratio of *τ /*Γ = 0.5 (or equivalently, *τ* = Γ*/*2), which was chosen to minimize the extrapolation from the *τ* = 500 ms measurements (given that most observed cardiac periods were around Γ= 1000 ms). Based on Equation 9, we computed: ![Formula][18]</img>  to yield a PI metric that was adjusted for individual cardiac periods and therefore comparable across subjects. For subjects with two *τ* = 500 ms scans, we performed the computation twice and averaged the PI(*τ* = Γ*/*2) values. Then, the relationship between PI(*τ* = Γ*/*2) vs age was assessed using the Pearson correlation coefficient.

## 5 Results

The fits of the subject data to the model in Equation 9 are shown in Figure 3, with R2 values ranging from 0.735 to 0.991 and minima occurring near *τ* = Γ as predicted by Equation 5. Notably, these subjects represent a diverse sample of cardiac periods, with values of Γ ranging from about 0.7 s to 1.1 s. Figure 4 serves as a complementary figure demonstrating the dependence of PI on *τ* for a representative subject (Subject 2) with a cardiac period around Γ = 1081 ms. In both the images and the curves, the pulsatility observed in the *τ* = 500 ms scans largely vanishes at *τ* = 1000 ms, where the image intensities and signal curves flatten out across the cardiac cycle. This was not unexpected given that the value of *τ* = 1000 ms is close to r, where the sinc ![Graphic][19]</img> term of Equation 5 evaluates to 0.

![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F3.medium.gif)

[Figure 3:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F3)

Figure 3: 
Fits of the modified sinc model (Equation 9) to the PI values measured from the *τ* -stepping protocol on each subject. The blue data points indicate the PI measured at each *τ* scan, with the black curve representing the best fit modified sinc model of Equation 9. The model fits the data quite well across a range of heart rates represented by these subjects. We consistently see PI reach a minimum around *τ* = Γ as predicted by Equation 5.

![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F4.medium.gif)

[Figure 4:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F4)

Figure 4: 
The left-hand montage shows the demeaned and normalized perfusion-weighted images (PWI) for a select slice of Subject 2, with cardiac phase *ϕ/*2*π* indicated along the x-axis and *τ* values for each scan indicated along the y-axis. To produce the maps on the left-hand side, the data are fit on a voxel-wise basis to obtain *S*(*ϕ*) curves for every voxel, then voxel-wise demeaned, and finally normalized by the scalar GM ROI mean value. The voxel intensities thus represent the relative fluctuation amplitude around each voxel’s mean compared to the baseline GM value. The voxels in the middle of the slices were masked due to ventricular location and CSF contamination (see Section 4.3). The accompanying line plots on the right-hand side show the GM ROI-averaged signal *S*(*ϕ*). The red shaded regions represents the 95% confidence interval at each value of *ϕ/*2*π*, derived from the set of *S*(*ϕ*) curves resulting from the residuals permutation approach. The text indicates the PI value, along with its 95% confidence interval. Note that for better visual comparison, the curves (and images) were phase shifted so that the maximum amplitude roughly aligns across the scans (rows) in this figure.

Figure 5 demonstrates the computation of a voxel-wise PI map using data from Subject 2 (corresponding to the same underlying data used in the first row of Figure 4). Following the same analysis procedures described above in Methods for the GM ROI signal, the VSASL signal *S*(*ϕ*) was computed on a voxel-wise basis, followed by a voxel-wise PI calculation via Equation 4.

![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F5.medium.gif)

[Figure 5:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F5)

Figure 5: 
Demonstration of a voxel-wise PI map using Subject 2 data, *τ* = 500 ms, run 1. (A) Perfusion signal as a function of cardiac phase *S*(*ϕ*), computed as in Section 4.4 on a voxel-wise basis. Note that as compared to the demeaned maps shown in the top row of 4, the voxel-wise mean is preserved in these maps. (B) Individual components of the Equation 4 formula corresponding to *S*max, *S*min, *S*max − *S*min (the numerator of the PI formula), *S*mean (the denominator of the PI formula), and the computed PI map. (C) Select slices of the same PI map. Note that in (B) and (C), the PI slices are masked in order to focus on gray matter voxels.

Figure 6 shows the repeatability of PI at *τ* = 500 ms across subjects. All data points lie close to the line of unity, with a high ICC (ICC = 0.960, *p <* 0.001) and correlation coefficient (*r* = 0.986, *p* = 0.002). A representative visual example of this repeatability is demonstrated in Figure 4 for Subject 2, where the top two rows show the repeats at *τ* = 500 ms. The appearances of the perfusion maps are nearly identical across the cardiac cycle. The morphologies of the GM ROI-averaged signal curves (Figure 4, right) are also nearly identical, with consistent shape and amplitude in both runs.

![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F6.medium.gif)

[Figure 6:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F6)

Figure 6: 
Repeatability of PI(*τ* = 500 ms) across subjects. The x-axis denotes the PI of run 1, and the y-axis denotes the PI of run 2. The data points are the measured values from the 5 subjects with repeats of the *τ* = 500 ms scan. The error bars indicate the 95% confidence intervals of the values. The values follow the line of unity closely, indicating excellent test-repeat repeatability. This is further supported by the high Pearson correlation coefficient *r* = 0.986 (*p* = 0.002) and ICC = 0.960 (*p<* 0.001).

Figure 7 shows a significant positive correlation between PI(*τ* =−*/*2) and age (*r* = 0.554, *p* = 0.021), showing the potential of this method to explore an association with age.

![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F7.medium.gif)

[Figure 7:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F7)

Figure 7: 
Plot of PI(*τ* = Γ*/*2) vs age, showing a positive correlation between PI(*τ* = Γ*/*2) and age (*r* = 0.554, *p* = 0.021).

## 6 Discussion

In this study, we have presented MVP-VSASL, a novel approach for measuring cerebral microvascular pulsatility by leveraging the microvascular specificity of VSASL and retro-spectively gating the VSASL signal. A key innovation of our study is the simple theoretical model for VSASL signal pulsatility PI given in Equation 5, which provides an explicit form for PI as a function of bolus duration *τ* and cardiac period Γ. This sinc model is validated by its excellent fits to in vivo data from 6 subjects representing a broad range of cardiac periods. Furthermore, the model predicts a minimum in PI near *τ* = Γ, which is empirically observed in all subjects. This feature can be explained intuitively: A VSASL bolus of duration of *τ* = Γ integrates over one whole period of the periodic CBF signal, regardless of the position of the bolus in the cardiac cycle. This results in minimal fluctuation of the gated VSASL signal and thus a minimal value of PI.

Using this sinc model of Equation 5, we also derived a theoretically optimal value of *τ* in Equation 7 to maximize the SNR of the pulsatility measurement. While the optimal value *τ*opt is slightly below Γ*/*2 based on theory, choosing *τ* = Γ*/*2 represents an SNR tradeoff of less than 3% for typical cardiac periods in the range of [0.6, 1.2] s. In this study, we used a pre-defined set of *τ* values in order to keep the scan configurations consistent across subjects. Using a fixed value of *τ* = 500 ms for the repeatability and age-related scans in this study corresponded to a value of *τ* that was reasonably close to the sample (both cohorts combined) mean Γ*/*2 value of 460 ms (range 380 ms to 570 ms) and yielded predicted SNR values that were on average within 7% of the optimal SNR. Future work will be useful for better determining the tradeoffs between using a fixed *τ* (e.g. *τ* = 500 ms) or selecting a scan-specific *τ* = Γ*/*2 (which requires adjusting the scan timing parameters on a per-scan basis).

We also show the potential of the metric PI(*τ* = Γ*/*2) for detecting physiological differences across subjects, as the positive association of PI(*τ* = Γ*/*2) with age (Figure 7) agrees well with previous studies finding higher flow pulsatility (as measured by PC-MRI) in older vs younger subjects8,10,22. While future work is needed to validate this association in larger populations, the agreement with prior literature supports the validity of the proposed VSASL approach.

### 6.1 Comparison with other techniques

The strength of MVP-VSASL is the ability to measure pulsatility in the microvasculature with vessel diameters around 50 µm. This is notably different from the majority of current pulsatility approaches, which focus on larger arteries feeding the brain. This focus on large vessels is a common feature across different modalities and techniques (e.g. 2D PC-MRI8,9, 4D flow11,36, doppler ultrasound4,12, or ASL13–15,17) and different measures of pulsatility (e.g. pulsatility of blood volume13–15,17, vessel distension16, flow velocity12, or flow volume rate4,8,36). For blood volume pulsatility, the exceptions to our knowledge are the dynamic ASL approach by Yan et al.14, where significant differences in blood volume between systole and diastole were found in the small arteries and arterioles, and the 7T VASO approach by Guo et al.37, which can measure relative blood volume over the cardiac cycle to estimate a vascular compliance index.

To our knowledge, the only other method that measures flow pulsatility in the microvas-culature is PC-MRI at 7T19–22, which targets cerebral perforating arteries with diameters below 300 µm and flow velocities as low as 0.5-1.0 cm/s. This approach involves resolving and measuring the flow through individual vessels, averaging the flow curves across vessels, and computing pulsatility on the averaged flow signal. For an approximate comparison, the 7T PC-MRI microvascular pulsatility values in the basal ganglia (a deep gray matter region) were reported to be around 0.55-0.65 for young subjects (mean age around 25 years)22, and around 0.7-1.0 for older subjects (mean age around 60-65 years)21,22. Using an age cutoff of 45 years, our PI(*τ* = Γ*/*2) values in cortical GM are 0.591 *±* 0.165 for younger subjects (n = 5, age 30.0 *±* 6.48 years) and 0.835 *±* 0.217 for older subjects (n = 12, age 64.7 *±* 9.56 years), which agree well with the 7T PC-MRI range of values. While the pulsatility values from 7T PC-MRI serve as a useful reference point, a direct comparison of their pulsatility values to those of the current study should be interpreted with caution due to (1) the different ROIs used between studies, and (2) the effect of bolus duration in MVP-VSASL. Specifically, the 7T PC-MRI studies have focused on the centrum semiovale (a white matter region) and basal ganglia, whereas the current VSASL study focuses on cortical GM. Future work assessing the same ROIs and exploring an appropriate adjustment for the *τ* -dependent sinc factor would be useful for a more direct comparison between MVP-VSASL and the 7T PC-MRI approach.

Our VSASL approach builds upon the prior body of work examining cardiac-driven fluctuations26,38–41 and pulsatility13–15,17 of the ASL signal. While the majority of this literature focused on spatially-selective ASL techniques24, Franklin et al.26 examined VSASL signal fluctuations weighted by macrovascular blood flow and volume (due to the absence of a VCM), and did not explicitly compute a pulsatility index. Other ASL-based pulsatility measurements have used spatially-selective ASL approaches24 to assess cerebral blood volume (CBV) pulsatility and vascular compliance in mainly large arteries13–15,17, whereas MVP-VSASL focused on a flow pulsatility measurement in the microvasculature. For flow-weighted ASL signals, some studies recognized that selecting a bolus duration equal to the cardiac period (or a multiple) could mitigate cardiac-driven fluctuations27,38, which we noted (both theoretically and empirically) as well. However, we are the first to our knowledge to derive a model (Equation 5) explicitly describing the dependence of cardiac fluctuations on *τ* and Γ, providing a useful theoretical extension to predict the magnitude of such signal fluctuations.

### 6.2 Clinical Considerations

An accurate and reliable measurement of microvascular pulsatility can be valuable for clarifying the mechanisms connecting pulsatility, microvascular damage, and cognitive decline in various diseases. Previous studies have associated decreased performance over a range of cognitive domains with pulsatility in large cerebral arteries4,9,12, but these measurements remain distant from the microvasculature, where the damage linked to various cognitive disorders is likely occurring4. In this regard, microvascular pulsatility measurements may be more suitable for assessing the relevant environment. Furthermore, since vascular risk factors such as microvascular pulsatility may reflect early and potentially modifiable variables in disease progression5–7,18,42, MVP-VSASL could also help assess cognitive risk before structural changes and clinical symptoms of cognitive impairment arise.

Because large arteries each supply pulsatile blood to extensive vascular territories that include many regions of the brain, the large-artery pulsatility indices are inherently limited in their spatial specificity when interpreting their impact on specific brain regions. In comparison, pulsatility measurements in the microvasculature, being much more distal along the arterial tree and embedded in the tissue they are supplying, can potentially be used to probe specific areas implicated in disease pathogenesis (e.g. the parietal lobe and hippocampus for Alzheimer’s disease). MVP-VSASL, with its whole-brain coverage and the ability to simultaneously generate microvascular signals across the entire brain, offers a unique potential to facilitate such regional assessments. By computing PI on a voxel-wise basis, we have demonstrated the feasibility of generating a voxel-wise pulsatility map as shown in Figure 5. To our knowledge, this is the first demonstration of a microvascular pulsatility map, and the robustness and applicability of the approach for regional pulsatility assessments will be explored in future work.

Importantly, our technique has several practical features that may ease translation, including the ability to perform at 3T, a simple prescription with whole-brain coverage, and a clinical scan time of around 6 minutes. These features make MVP-VSASL an attractive option for integration into standard clinical and research brain MRI protocols.

### 6.3 Technical Considerations

In this study, we used the standard VSASL cutoff velocity of *v*cut = 2 cm/s23, which defines the VSASL bolus in the microvascular regime. However, the value of *v*cut can be adjusted to target other segments of the arterial network. Increasing *v*cut shifts the VSASL bolus-defining location more upstream into larger vessels, whereas decreasing *v*cut shifts the labeling region further distally toward the capillaries. Varying *v*cut could thus be useful for evaluating pulsatility along the arterial tree using a single technique (by adjusting *v*cut) and assessing metrics such as the vascular dampening factor8. Of note, decreasing *v*cut requires higher gradient strengths, which can exacerbate known technical issues like eddy currents, diffusion effects, and CSF contamination, but further improvements to the VSASL pulse train could make this exploration more feasible in the future23.

Prior studies on VSASL have explored its potential to measure venous flow43,44. The current analyses could also be applied to measure venous flow pulsatility, which has already received some attention using other techniques11.

A potential challenge of our approach is CSF contamination of the VSASL signal. This is an existing issue with VSASL in general, as CSF is difficult to suppress through standard background suppression techniques23. When using VSASL for quantifying mean cerebral blood flow, this manifests as an overestimation of mean perfusion in CSF-contaminated regions. For our study, CSF contamination may also affect signal fluctuations as CSF can also pulsate with the cardiac cycle45. Our pre-processing attempted to mitigate CSF effects as much as possible by aggressively excluding contaminated areas from the GM masks, but we acknowledge that it is not possible to completely eliminate CSF partial volume effects. Furthermore, these effects are challenging to remove through modeling due to the complex and heterogeneous nature of CSF flow. However, as VSASL sequences continue to improve, our approach can benefit from future advances related to CSF suppression46,47.

Finally, we used a standard BIR-8 train for VSASL labeling in this proof-of-principle study. The BIR-8 train is a velocity-selective saturation (VSS) technique noted for its robustness to ![Graphic][20]</img> inhomogeneity and eddy currents, making it well suited for this initial study. Velocity-selective inversion (VSI) approaches that provide increased SNR48 were attempted, but resulted in too many artifacts particularly at lower bolus duration values of *τ*. While the dynamic phase cycling method mitigates many of these artifacts49, the resulting four-fold decrease in temporal resolution renders the retrospective gating method impractical within a reasonable imaging time. More robust VSI techniques better suited for our pulsatility method are currently under investigation50.

## 7 Conclusion

We have introduced MVP-VSASL, a novel, practical and non-invasive approach of using VSASL to measure pulsatility in the microvasculature, including a theoretical framework relating the pulsatile VSASL signal to the pulsatility index. This technique may be used to probe the mechanisms underlying cognitive disease, monitor disease progression, and evaluate patient responses to therapy.

## Data Availability

All data produced in the present study are available upon reasonable request to the authors.

## Supporting Information

### S1 Evaluating Sinc Model Approximation

In this section, we evaluate the error associated with the sinc model of pulsatility (Equation 5), whose derivation involves neglecting the 2nd-order Fourier terms. To evaluate the error, an exact expression for PI(*τ*) of a 2nd-order Fourier signal is derived. Then the exact PI(*τ*) curves (the reference) are computed over a range of *τ* values, and Equation 5 (the model) is fit to the reference to evaluate the model error.

#### S1.1 Theory

To provide an exact expression of PI, we start with the 2nd-order Fourier model for *S*(*t*) given in Equation 3, with the summation expanded out for clarity in this section. Note the *T*1 exponential decay factor is excluded for simplicity, as it cancels out in the final PI calculation. ![Formula][21]</img>  

By combining the cos(*·*) and sin(*·*) terms of the same order, the equation is rewritten as: ![Formula][22]</img>  where ![Formula][23]</img>  represent the magnitudes of the 1st and 2nd-order Fourier components, respectively. To compute PI, we can apply Equation 4 (reproduced below): ![Formula][24]</img>  by evaluating the individual terms of the formula. The mean term *S*mean is simply the constant term: ![Formula][25]</img>  where *C* is the proportionality constant (which accounts for the switch from the ∝ symbol in Equation S.2 to the = symbol in Equation S.5). The max term *S*max is computed by taking the maximum of *S*(*t*): ![Formula][26]</img>  

Because the max[*·*] operation is insensitive to phase shifts, *S*max can be further rewritten as: ![Formula][27]</img>  where *ϕ*el is defined as: ![Formula][28]</img>  and represents the relative phase shift between the 1st- and 2nd-order sinusoids. Following the same logic, the min term *S*min can be represented in a similar way: ![Formula][29]</img>  

By substituting *S*mean, *S*max and *S*min into Equation 4, we obtain the following exact expression for PI(*τ*): ![Formula][30]</img>  

This is a formula for PI that depends on the scaling terms *b*, *W*1, *W*2 and a relative phase shift between 1st- and 2nd-order cosines given by *ϕ*rel. Note that if the 2nd-order component is neglected by setting *W*2 = 0, then the formula simplifies down to the sinc model stated in Equation 5 of the main text and reproduced below for convenience: ![Formula][31]</img>  where ![Graphic][32]</img> is the lumped fitting parameter described in the main text.

#### S1.2 Methods

To evaluate the errors associated with the sinc model approximation, we: (1) computed exact PI(*τ*) curves (Equation S.10), (2) fit the sinc model (Equation 5) to those exact curves, and (3) assessed the error of those fits.

To produce the exact PI(*τ*) curves, the constants (*b*, *W*1, *W*2) were set at physiologically reasonable values based on our empirical measurements and on prior literature. First, *b* = 1 and *W*1 = 0.5 were chosen to ensure the PI(*τ*) curves have values comparable to in vivo data measurements at *τ* = 500*/*750*/*1000*/*1250*/*1500 ms. For example, this choice yields a value of PI(*τ* = Γ*/*2) = 0.636, which is in the middle of the range of values shown in Figure 7. While the exact *b* and *W*1 values are arbitrary, their ratio controls the vertical scaling factor ![Graphic][33]</img> in the Equation 5 sinc model, which roughly scales the exact PI(*τ*) curves as well. To choose *W*2, prior estimates of blood flow power spectra in the internal carotid artery (ICA) were used to set the power ratio of the 1st-vs 2nd-order components ![Graphic][34]</img> as 5:11, which results in ![Graphic][35]</img>. (Note that the power ratio in the microvasculature may even be higher than 5:1, since arterial compliance generally behaves like a low-pass filter on flow waveforms. This would comparatively attenuate the 2nd-order power as blood travels from the ICA to the microvasculature and result in an even smaller relative value of *W*2.) To summarize, the constants are set at *b* = 1, *W*1 = 0.5 and *W*2 = 0.2236.

A cardiac period of Γ = 1 s was set based on the range observed in subjects (approximately [0.7, 1.1] s). Bolus duration *τ* was varied over *τ* ∈ [0.400, 1.600] s in steps of 0.001 s. The remaining unconstrained parameter of *ϕ*rel (the relative phase shift between 1st- and 2nd-order components) was varied between 0 to 2*π* in steps of 2*π τ* 0.001.

For each combination of *ϕ*rel and *τ*, the exact value of PI(*τ*) (Equation S.10) was computed by numerically evaluating its terms over *t* ∈ [0, Γ] in steps of 0.0001 s. This resulted in curves of PI(*τ*) for every value of *ϕ*rel. Finally, for each value of *ϕ* rel, the sinc model (Equation 5) was fit to the exact PI values at *τ* = 500*/*750*/*1000*/*1250*/*1500 ms, which correspond to the *τ* values used in the *τ* -stepping experiment.

Since the derivation of the sinc model (Equation 5) involves neglecting 2nd-order terms of the signal curves, a natural question is whether we could simply start with a 1st-order Fourier model of the signal *S*, which would obviate the need for any approximation step for deriving the PI(*τ*) expression. As a complementary empirical analysis, we also re-analyzed the *τ* - stepping data (Figure 3 of the main text) by using a 1st-order Fourier model to describe *S*(*ϕ*), and compared the results to those of the 2nd-order approach. The same procedures outlined in Section 4.4 were followed, with the exception of using a 1st-order Fourier model to fit the control and label data of each individual scan. After fitting the *τ* -stepping data to the sinc model (Equation 9), the resulting values of *A* (i.e. the scaling factor of the sinc model) and R2 of the fits were also compared between the 2nd-vs 1st-order approaches using a paired Wilcoxon signed rank test.

#### S1.3 Results

The results of fitting the sinc model (Equation 5) to the exact PI values at *τ* = 500*/*750*/*1000*/* 1250*/*1500 ms are shown in Figure S1, with Figure S1A showing the best case fit based on R2 (over values of *ϕ*rel), Figure S1B showing the worst case fit, and Figure S1C the error across the entire space of *ϕ*rel. There is very little error across the entire space of *τ* and *ϕ*rel, supporting the validity of the sinc model (Equation 5) for describing pulsatility as a function of *τ*.

Figure S2 shows the results of the empirical *τ* -stepping analysis comparison. As shown in the insets of Figure S2A and Figure S2B, the main qualitative difference between the 2nd-vs 1st-order approaches is that the 1st-order model of *S*(*ϕ*) has difficulty describing relatively sharper peaks (which are captured by the 2nd-order approach), an effect most clearly shown for the *τ* = 500 ms scans. As a result, the 1st-order approach tends to underestimate PI for each scan, which can be seen in the individual subject examples (i.e. the orange data points being consistently lower than the blue ones). This effect is also summarized in Figure S2C, which shows a scatter plot of the PI values. The x-axis indicates the PI values resulting from the 1st-order approach, and the y-axis indicating the PI values resulting from the 2nd-order approach. Almost all points lie above the line of unity, indicating consistently higher PI values with the 2nd-order approach. However, the ability of the sinc model to describe the PI(*τ*) data seems to remain generally strong regardless of the 2nd-vs 1st-order approach taken to analyze the data for each individual *τ* scan. As shown in Figure S2D, the R2 values of the fits are all around 0.75 or higher (with one outlier exception in Subject 3 for the 1st-order approach), and the Wilcoxon signed rank test shows no statistically significant difference between the R2 values (*p* = 0.094). In summary, while using the 1st-order approach of *S*(*ϕ*) seems to underestimate the PI values, the sinc model still appears to describe the PI(*τ*) shape well with either approach in these comparisons.

### S2 Computing VSASL Signal and Pulsatility

This section provides supporting details for Section 4.4 of the main text, and describes the computation of the VSASL signal *S*(*ϕ*) as a function of cardiac phase and the subsequent computations of PI and confidence intervals.

The computation of *S*(*ϕ*) involves (1) retrospective gating to map control and label data (acquired in time) to cardiac phase space, (2) fitting 2nd-order Fourier models to control and label data points, and (3) subtracting the fits to obtain *S*(*ϕ*). We start with (*N ×* 1) control and label measurement vectors: ![Formula][36]</img>  where *y*C,*i* and *y*L,*i* represent the *i*th control and label measurements, respectively, and *N* is the number of control/label pairs acquired from the scan. These measurement vectors can represent an ROI-averaged signal (as in the GM ROI-averaged signal used in Section 4.3) or the values from single voxels (as was used to produce the voxelwise PI map in Figure 5).

Next, a cardiac phase *ϕ* was assigned to every label and control data point by retrospectively gating on the photoplethysmography (PPG) waveform trace. First the peaks of the PPG waveforms were detected. Then, a cardiac phase *ϕ*C,*i* or *ϕ*L,*i* was assigned to every control and label volume, respectively, using: ![Formula][37]</img>  where *t* is set to *t*C,*i* or *t*L,*i* (defined as the center of the labeling/control period with width *τ*) for each volume, *t*1 is the preceding cardiac peak, and *t*2 is the subsequent cardiac peak2. For a given data point, time *t*C,*i* or *t*L,*i* was computed by subtracting *PLD* + *τ /*2 from the start time of the readout. The phase values were concatenated into (*N ×* 1) vectors for the controls and labels: ![Formula][38]</img>  where *ϕ*C,*i* and *ϕ*L,*i* represent the cardiac phase of the *i*th control and label measurements, respectively. Then (*N ×* 5) 2nd-order Fourier design matrices were constructed as: ![Formula][39]</img>  

To fit **Y**C and **Y**L to separate 2nd-order Fourier models, the data are modeled as: ![Formula][40]</img>  where **c** and **l** are the 2nd-order Fourier coefficients for controls and labels, respectively. The terms ***ϵ***C and ***ϵ***L represent measurement noise, which we assume are represented by independent and identically distributed (i.i.d.) normal random variables with variances ![Graphic][41]</img> and ![Graphic][42]</img>, respectively. (Normality of noise in MRI magnitude images is a reasonable approximation when the SNR of a given voxel is sufficiently high. The noise begins to resemble a normal distribution at SNRs as low as ∼33. The control and label data of this study comfortably lie within the high SNR regime, with empirical voxel-wise SNRs on the order of ∼50 for cortical gray matter voxels.) The coefficients **c** and **l** are then estimated via least squares: ![Formula][43]</img>  with **ĉ**= [*ĉ* *ĉ*cos,1 *ĉ*sin,1 *ĉ*cos,2 *ĉ*sin,2 ]*T* and ![Graphic][44]</img> denoting the cor-responding estimates. The difference (signal *S*) coefficients are then computed by control-minus-label subtraction: ![Formula][45]</img>  where ![Graphic][46]</img>. The signal curve *S*(*ϕ*) can then be obtained by inserting the coefficients into a 2nd-order Fourier model (Equation 8 of the main text, reproduced below): ![Formula][47]</img>  

(Equivalently, the control and label coefficients could be inserted into separate 2nd-order Fourier models, followed by subtracting the models.) To compute pulsatility, *S*(*ϕ*) are evaluated over a fine grid of *ϕ* values (for example, *ϕ* ∈ [0, 2*π*] in steps of 2*π* · 0.0001), and then PI can be numerically calculated with Equation 4 of the main text as: ![Formula][48]</img>  

The stability of the PI measurements can be assessed using a residuals permutation approach. We first computed the data estimates ![Graphic][49]</img> and ![Graphic][50]</img> as: ![Formula][51]</img>  and then computed the residuals ![Graphic][52]</img> C and ![Graphic][53]</img> as: ![Formula][54]</img>  

Then, we randomly permuted the residuals in time (separately for labels and controls) and added them back to the fitted data to create simulated data vectors **Y**C,*Perm,i* and **Y**L,*Perm,i*: ![Formula][55]</img>  where *P*C,*i*(*·*) and *P*L,*i*(*·*) denote the permuting/shuffling operators for the control and label data, respectively, for the *i*th iteration. Using these simulated data vectors, PI was computed again following Equations S.17-S.20. This procedure was repeated over 1000 iterations to generate a distribution of PI values, which were then used to compute 95% confidence intervals.

### S3 Impact of Measurement Noise on Pulsatility

This section assesses the bias and SNR of the PI measurement under measurement noise (i.e. noisy control and label data points). We derive the theoretical relationship between measurement noise and the distributions of the estimated difference coefficients ![Graphic][56]</img> and ![Graphic][57]</img> (which were estimated in Equation S.18). The impact of noise on ![Graphic][58]</img> is then assessed via simulations and theoretical approximations of the SNR of PI.

#### S3.1 Theory

##### S3.1.1 SNR and Bias of Pulsatility Index

We first approach this problem by assuming ground truth 2nd-order Fourier coefficients for the control, label and difference signals as follows: ![Formula][59]</img>  respectively. Following the framework outlined in Section S2, we assume a VSASL scanning experiment that has acquired (*N ×* 1) control and label data vectors **Y**C and **Y**L which follow the model described in Equation S.16. By substituting the expressions of **Y**C and **Y**L (Equation S.16) into the least squares coefficient estimates **ĉ** and ![Graphic][60]</img> (Equation S.17), we obtain: ![Formula][61]</img>  

These expressions can then be substituted into Equation S.18 to obtain our estimate of difference coefficients ![Graphic][62]</img> as: ![Formula][63]</img>  

In a realistic experiment, the entries of ***ϕ***C (and ***ϕ***L) have values in [0, 2*∈*] that are unevenly spaced across the interval. However, when we have a sufficient number of measurements (e.g. *N* = 72 as in our experiments), approximating ***ϕ***C (and ***ϕ***L) as vectors with uniformly spaced entries has negligible impact on the **ĉ** (and ![Graphic][64]</img>) estimation (results not shown). This evenly-spaced approximation is: ![Formula][65]</img>  whose entries are given by: ![Formula][66]</img>  

Consequently, the design matrices are also approximated as: ![Formula][67]</img>  and Equation S.26 simplifies to: ![Formula][68]</img>  where ![Graphic][69]</img> and ![Graphic][70]</img>. By assessing the expected value and covariance, we see that ![Graphic][71]</img> can be described as a normally distributed random vector: ![Formula][72]</img>  where ![Formula][73]</img>  with off-diagonal terms equal to 0 due to the orthogonality of the evenly-spaced 2nd-order Fourier regressors.

This expression of ![Graphic][74]</img> can then be used to examine the behavior of the pulsatility index measurement. Inserting the coefficients into *S*(*ϕ*) yields the estimated VSASL signal curve (Equation 8 of the main text, reproduced below): ![Formula][75]</img>  

Using similar arguments as in the derivation of Equation 5 of the main text to neglect second-order terms, the pulsatility index estimate ![Graphic][76]</img> is then calculated via Equation 4 of the main text as: ![Formula][77]</img>  

The SNR of ![Graphic][78]</img> is then defined as: ![Formula][79]</img>  where ![Graphic][80]</img> is the expected value and ![Graphic][81]</img>. The bias of ![Graphic][82]</img> is defined as: ![Formula][83]</img>  

##### S3.1.2 Approximation of SNR

When examining ![Graphic][84]</img> (Equation S.34), its numerator can be identified as a Rician random variable which we denote as ![Graphic][85]</img>: ![Formula][86]</img>  where ![Graphic][87]</img> and ![Graphic][88]</img> denote the non-centrality and scale parameters of the distribution, respectively. On the other hand, its denominator can be identified as a normal random variable as described by Equation S.32 ![Formula][89]</img>  

Thus, the pulsatility estimate of ![Formula][90]</img>  is the ratio of a Rician random variable and a normal random variable. To our knowledge, this exact distribution has not been analytically described in prior literature. However, through several simplifications, we will derive an approximation of ![Graphic][91]</img> that facilitates the PI SNR model presented in Equation 6 of the main text. The subsequent equations will focus on the derivations, and the errors associated with each step will be examined further in the Methods and Results.

To illustrate these approximations, we start with a random variable *Z* defined as the ratio of *X* and *Y* ![Formula][92]</img>  where *X* ∼ *Rice*(*µ**x*, *σ**x*) has a non-centrality parameter *σ**x* and scale parameter *σ**x*, and ![Graphic][93]</img> has a mean *µ**y* with standard deviation *σ**y*. Provided a sufficiently high value of *µ**x**/σ**x* (e.g. *>* 5), the Rician random variable *X* can be approximated as a normal random variable ![Graphic][94]</img>, which yields a ratio of two normal random variables. Provided a sufficiently high value of *µ**y**/σ**y* (e.g. *>* 10 4,5), which makes the denominator *Y* unlikely to observe negative values, the ratio of two independent normal random variables can be approximated as a single normal random variable5 of the form: ![Formula][95]</img>  

This approximation can be applied to ![Graphic][96]</img> by identifying the corresponding terms ![Graphic][97]</img> and ![Graphic][98]</img>. Before applying the ![Graphic][99]</img> approximation, we first re-express these corresponding terms to elucidate their dependence on *τ*. To do so, we compare the Fourier model of *S*(*t*): ![Formula][100]</img>  with the expanded form of Equation 3 of the main text: ![Formula][101]</img>  (Note that Equation S.42 is consistent with Equation 8 of the main text by using the sub-stitution *ϕ* = 2*πt/*Γ.) As indicated above by the underbraces, this comparison yields the following definitions for *d*, *d*cos,*k* and *d*sin,*k*: ![Formula][102]</img>  where *C* is the constant of proportionality. Using these definitions for *d*, *d*cos,1 and *d*sin,1 (with the latter two obtained by setting *k* = 1), the terms *µ**x*, *σ**x*, *µ**y* and *σ**y* can be reexpressed as: ![Formula][103]</img>  ![Formula][104]</img>  ![Formula][105]</img>  ![Formula][106]</img>  

Finally, these terms are applied to Equation S.41 to yield the ![Graphic][107]</img> approximation: ![Formula][108]</img>  where ![Graphic][109]</img> per Equation 5 of the main text (and the subscript ”approx,1” denotes the first of two approximations, with the second coming below). Using Equations S.45-S.48, Figure S3 serves to qualitatively show the regimes of *τ* where the approximations used in Equation S.41 (based on the ratios *µ**x**/σ**x* and *µ**y**/σ**y*) generally hold.

Computing the SNR of ![Graphic][110]</img> (via Equation S.35) yields: ![Formula][111]</img>  

By examining the ![Graphic][112]</img> factor (which is present in the standard deviation of ![Graphic][113]</img> and the denominator of SNRapprox,1) in a plausible physiological regime, we can make a further simplification. Assuming values of *b* = 1 and ![Graphic][114]</img> (as in SectionS1.2), the ![Graphic][115]</img> factor simplifies into ![Graphic][116]</img>. Note that the ![Graphic][117]</img> term has a maximum value of 1 (at *τ* = 0 s) that quickly approaches 0 with increasing *τ* due to its envelope of (Γ*/ πτ*)2, thus representing a negligible contribution to the overall value of the factor. For example, assuming Γ = 1 s, neglecting ![Graphic][118]</img> at values of *τ* = 0*/*0.5*/*1.5 s produces errors of 0.057*/*0.038*/*0.013 when expressed as a fraction of ![Graphic][119]</img> respectively. As empirical support, one can also note that the measured values of PI were all on the order of 1 or less (Figs. 3 and 7). Thus, for a reasonable physiological regime, we can further approximate ![Graphic][120]</img>, which yields another approximation of ![Graphic][121]</img>: ![Formula][122]</img>  and analytical model of the SNR: ![Formula][123]</img>  

Note that by expanding out the ![Graphic][124]</img> term, this final expression yields the SNR model presented in Equation 6 of the main text (reproduced below): ![Formula][125]</img>  

Thus, Equation S.52 is also optimized at *τ*opt, given by Equation 7 of the main text.

#### S3.2 Methods

Using the theory derived in Section S3.1.1, we use simulations to assess the impact of noise on the bias and SNR of PI measurements, particularly as a function of *τ*. To do so, we (1) assume a noiseless reference PI(*τ*) curve, (2) add measurement noise, (3) compute distributions of ![Graphic][126]</img> at every *τ* value, and (4) assess the bias and SNR of PI(*τ*) as a function of *τ*.

We begin by assuming reference CBF(*t*) coefficients *b*, *b*cos,1, *b*sin,1, *b*cos,2 and *b*sin,2. We assume similar values as in Section S1 (the sinc model error simulations) by using *b* = 1 and ![Graphic][127]</img>, with *b*cos,2 = *b*sin,2 = 0 also assumed to neglect the 2nd-order component for simplicity. The ![Graphic][128]</img> constraint was satisfied by specifying ![Graphic][129]</img> (any pair of values satisfying the constraint will yield equivalent results).

Then, over a grid of *τ œ* [0.001, 2.000] s in steps of 0.001 s, the corresponding *S*(*t*) coefficients **d** = [ *d* *d*cos,1 *d*sin,1 *d*cos,2 *d*sin,2 ]*T* are computed for every value of *τ* based on Equation S.44 using a cardiac period Γ = 1.0 s, and *T*1*b* = 1.6 s (with the proportionality constant *C* arbitrarily set to 1 for simplicity, since *C* cancels out in the eventual PI calculation). Once computed, *d*, *d*cos,1, *d*sin,1, *d*cos,2 and *d*sin,2 are inserted into the *S*(*ϕ*) 2nd-order Fourier model (Equation S.19) and then PI is computed (Equation S.20). For example, at *τ* = 500 ms, the coefficients are computed as *d* = 0.366, *d*cos,1 = *d*sin,1 = 0.0825 and *d*cos,2 = *d*sin,2 = 0, which produce a pulsatility value of PI = 0.636. By repeating this procedure for every *τ* value, the ground truth (reference) PI(*τ*) curve is obtained.

Next, we add noise to the simulations and compute ![Graphic][130]</img>. In order to add a realistic amount of noise, we used the in vivo GM ROI *τ* = 500 ms scan analysis results to choose an appropriate value of ![Graphic][131]</img>. First, ![Graphic][132]</img> and ![Graphic][133]</img> were estimated as the variance of the in vivo control and label fit residuals (Equation S.22) at *τ* = 500 ms, and then ![Graphic][134]</img> was computed as ![Graphic][135]</img>. Then, the ratio of *σ* D*/d* was estimated to be *≈* 0.10 (on average across subjects). Since *d* = 0.366 at *τ* = 500 ms in these noise simulations, we chose *σ* D = 0.0366 to produce an equivalent ratio of *σ* D*/d*. This value of *‡*D = 0.0366 (i.e. noise level) was then kept constant across all values of *τ*, consistent with the empirical observation that the residuals at *τ* = 750*/*1000*/*1250*/*1500 ms were comparable in magnitude to the residuals at *τ* = 500 ms.

The number of measurements was set at *N* = 72. The covariance matrix **C** was computed using Equation S.32, and then 105 realizations of ![Graphic][136]</img> were simulated following ![Graphic][137]</img> (Equation S.31). We then computed ![Graphic][138]</img> for each realization (Equation S.20) to produce a distribution of ![Graphic][139]</img> values at every value of *τ*.

We then evaluated ![Graphic][140]</img> based on its SNR (Equation S.35) and bias (Equation S.36) at every value of *τ*, with the SNR curve then numerically assessed for an optimal value of *τ* : ![Formula][141]</img>  

The analytical expressions of SNRapprox,1 (Equation S.50) and SNRapprox,2 (Equation S.52) were also evaluated at every value of *τ*, and similarly assessed for their optimal values as well: ![Formula][142]</img>  ![Formula][143]</img>  

#### S3.3 Results

In Figure S4A, we see generally good agreement between the noiseless reference curve PI and the expected value ![Graphic][144]</img>, with little spread in ![Graphic][145]</img> as indicated by the shaded areas. However, ![Graphic][146]</img> becomes less accurate as *τ* approaches 0 s and the noise begins to dominate. While both the numerator and denominator of the calculation are affected by noise, a noisy denominator ![Graphic][147]</img> can be particularly problematic for the ![Graphic][148]</img> calculation, as individual realizations of *d* can approach 0 (or even cross into negative values) and make the ![Graphic][149]</img> quotient approach infinity. As a result, ![Graphic][150]</img> can become unstable in very low-*τ* regimes, which is not surprising due to the low mean signal.

In Figure S4B, the asymptotic behavior of the bias ![Graphic][151]</img> and Std ![Graphic][152]</img> curves for *τ ⟶* 0 reflects the aforementioned stability issues. Another defining feature is the spike in bias around *τ* = Γ (and *τ* = 2 Γ). Here, the reference ![Graphic][153]</img> is 0, but ![Graphic][154]</img> values remain strictly positive due to the *S*max − *S*min numerator of Equation 4, resulting in the positive bias. However, the bias is otherwise negligible for most values of *τ*, including regimes around *τ*opt.

Figure S4C shows the simulation-based SNR curve (Equation S.35) alongside SNRapprox,1 (Equation S.50) and SNRapprox,2 (Equation S.52), demonstrating excellent agreement in their overall shapes. Furthermore, their respective maxima occur at very similar values of *τ* (represented by the vertical lines). Of note, the difference between SNR(*τ*opt) and SNR(*τ* *ú*) is less than 0.1% (18.89 vs 18.90), showing that *τ*opt indeed approximately maximizes SNR. Altogether, the agreement in curve shapes and negligible difference between SNR(*τ*opt) and SNR(*τ* ***) supports the use of the simple theoretical model (Equation 6) to describe the *τ* - dependence of SNR and guide the selection of an optimal *τ* value.

Examining the curves more closely, SNRapprox,1 is nearly identical to SNR except around *τ* = Γ and *τ* = 2 Γ. The difference in these regimes is expected, as the assumption of normality of the ![Graphic][155]</img> numerator, built into Equation S.41 (and thus into SNRapprox,1), diverges from its Rician nature in these low-SNR regimes3. Figure S5 serves as a supporting figure by comparing the simulated ![Graphic][156]</img> distribution (blue histogram) and the analytical PDF of ![Graphic][157]</img> (purple solid line), which approximates the ![Graphic][158]</img> distribution poorly at *τ* = Γ and *τ* = 2 Γ, but otherwise shows excellent agreement for the *τ* values shown. Examining SNRapprox,2, its curve (green dotted line) is nearly identical to SNRapprox,1 but is observed to be slightly higher (noticeable in the regime around 0.2 s to 0.6 s) due to the PI2 term that was neglected in the derivation step from Equation S.50 to Equation S.52.

### S4 T1-weighted Structural Scan Parameters

Two *T*1-weighted structural scan configurations were used with similar parameters and both based on a Magnetization-Prepared Rapid Acquisition Gradient Echo (MPRAGE) sequence. The main difference was the acceleration factor, with one using 4*×*-acceleration for shorter scan time, and the other using 2*×*-acceleration to maintain higher SNR for a separate project. The parameters for the former configuration (with the latter configuration noted in paren-theses wherever different) were: resolution = 0.9375 mm isotropic (1 mm isotropic), matrix size = 176 *×* 256 *×* 256, 176 slices in the sagittal orientation, flip angle = 8 degrees, TR = 2300 ms (2500 ms), TE = 2.29 ms (2.88 ms), inversion time = 900 ms (1060 ms), GRAPPA = 4*×* (2*×*) acceleration in the A-P phase-encode direction, acquisition time = 3:08 (7:30).

### S5 Iterative Denoising

This section provides supporting details for the iterative denoising method referenced in the main text Section 4.3, which was adapted from Power et al.6 with a few modifications. Overall, the algorithm involves identifying outlier volumes over successive iterations, followed by a final censoring of outlier timepoints and subtraction of nuisance variance from the data.

An initial temporal censoring mask was constructed by identifying volumes with framewise displacement (FD) exceeding 0.75mm (with FD computed as in6), and volumes with a cardiac period Γ more than 3 scaled median absolute deviations away from the median of the Γ timecourse. The outliers in Γ often corresponded to index finger motion, which distorted the PPG trace and caused unreliable retrospective gating.

For each iteration, the algorithm (1) censors timepoints identified by the (current) temporal censoring mask from the original data and regressors, (2) regresses out nuisance variance and (3) updates the temporal censoring mask with additional volumes identified based on the (current) regression residuals. In our implementation, nuisance regressors were constructed (motion timecourses, motion first-derivative timecourses, and first- and second-order low-frequency drift). To preserve cardiac-driven fluctuations, we also included 2nd-order Fourier regressors of controls and labels, which were constructed by starting with **X**C and **X**L (derived in Equation S.15 in SI Section S2) and interleaving with 0’s to account for the alternating control/label ASL acquisition scheme. For example, 0’s were inserted into the **X**C regressors for every row corresponding to a label data point (and vice versa for the **X**L regressors at control data points). Figure S6A shows how they are interleaved with 0’s and then concatenated alongside the nuisance regressors. Figure S6B shows **X**C and **X**L separately and sorted in cardiac phase *ϕ* order to better illustrate the shape of these regressors.

The nuisance and Fourier regressors were both input into the iterative denoising algorithm, and additional timepoints were censored with each iteration based on DVARS and SD metrics as in the original Power et al. paper6. The algorithm was stopped when no new timepoints were censored for a given iteration. Then the nuisance regressors were orthogonalized with respect to the Fourier regressors, and the orthogonalized nuisance fits were subtracted out.

## S7 Figures

![Figure S1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F8.medium.gif)

[Figure S1:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F8)

Figure S1: 
Evaluating Sinc Approximation Error: (A) Best case fit of Equation 5 to the exact PI(*τ*) curve. The curves are nearly identical, with negligible absolute error (*<* 0.01) over all values of *τ* ∈ [0.4, 1.6] s and R2 = 0.9999. (B) Worst case fit of Equation 5 to the exact PI(*τ*) curve. Even in this worst case, the agreement between curves is excellent, with minimal absolute error (*<* 0.05) over all values of *τ* ∈ [0.4, 1.6] s and R2 = 0.9922. (C) Image of the absolute error across all values of *ϕ* rel ∈ [0, 2 *π*].

![Figure S2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F9.medium.gif)

[Figure S2:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F9)

Figure S2: 
Comparison of the *τ* -stepping results when using 2nd-vs 1st-order Fourier models to describe *S*(*ϕ*). Shown in (A) and (B) are results from two representative subjects (4 and 5). Plotted in blue are the *τ* -stepping results using a 2nd-order Fourier model (identical to those shown in Figure 3 of the main text). Plotted in orange are the PI values when using a 1st-order Fourier model. The insets show comparisons of the *S*(*ϕ*) curves for each approach. Shown in (C) is a scatter plot of PI values, with the 1st-order approach on the x-axis and the 2nd-order approach on the y-axis. Shown in (D) is a comparison of the R2 of the sinc model fits.

![Figure S3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F10.medium.gif)

[Figure S3:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F10)

Figure S3: 
A supporting figure for qualitatively understanding the regimes of *τ* where the approximations used in Equation S.41 generally hold. The expressions of *µ**x*, *σ* *x*, *µ**y* and *σ* *y* were evaluated using Equations S.45-S.48, with values of ![Graphic][159]</img>, Γ = 1 s, *σ* D = 0.0366 and *N* = 72 that were used in Section S3.2. The blue regions indicate where *µ**x**/ σ* *x* is less than a threshold of 5, where the normal approximation of the Rician numerator of ![Graphic][160]</img> becomes relatively poor. The orange regions indicate where *µ**y**/ σ* *y* is less than a threshold of 104,5, where the approximation of the ratio of two normal random variables as a single normal random variable becomes relatively poor. The green regions indicate the regimes where both thresholds are exceeded and the approximations hold relatively well. These green areas are consistent with the regimes in Figure S4C showing excellent agreement between SNR and SNRapprox,1.

![Figure S4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F11.medium.gif)

[Figure S4:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F11)

Figure S4: 
Noise simulations for ![Graphic][161]</img>. (A) Reference PI curve (red), expected value of ![Graphic][162]</img> (blue dashed line), and the spread of ![Graphic][163]</img> (shaded regions corresponding to 68% and 95% confidence intervals). (B) The bias (orange) and standard deviation (cyan) of ![Graphic][164]</img> as a function of *τ*. (C) The SNR of the ![Graphic][165]</img> measurement (blue solid line) and approximations SNRapprox,1 (Equation S.50, purple dashed line) and SNRapprox,2 (Equation S.52, green dash-dotted line), along with the *τ* values where their maxima occur (i.e. ![Graphic][166]</img> and *τ*opt, respectively) denoted by the vertical lines. The plotted points indicate the SNR (blue solid line) evaluated at those *τ* values.

![Figure S5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F12.medium.gif)

[Figure S5:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F12)

Figure S5: 
This is a supporting figure for understanding the SNR curves shown in Figure S4C. Shown are the simulated distributions of ![Graphic][167]</img> (blue histogram) compared to the analytical PDFs of ![Graphic][168]</img> (purple solid line) and ![Graphic][169]</img> (green dotted line), for a range of selected *τ* values in steps of 0.250 s (except the first step at 0.100 s). Also indicated are the reference PI value (red vertical line) and ![Graphic][170]</img> (blue vertical dashed line). At each value of *τ*, the texts show the values involved in computing the SNR. The histogram and its analytical approximations (![Graphic][171]</img> and ![Graphic][172]</img>) generally agree well, except at *τ* = Γ and 2 Γ, where the built-in normal approximation of the Rician numerator of ![Graphic][173]</img> no longer holds.

![Figure S6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/23/2024.06.21.24309261/F13.medium.gif)

[Figure S6:](http://medrxiv.org/content/early/2024/06/23/2024.06.21.24309261/F13)

Figure S6: 
Illustration of the 2nd-order Fourier models in the iterative denoising algorithm. (A) The overall design matrix used in the iterative denoising algorithm of Section 4.3 of the main text. The 2nd-order Fourier matrices, representing separate models for controls and labels are sorted in time according to the temporal acquisition order of the control and label volumes and interleaved with 0’s according to the label/control alternation. They are then concatenated alongside the nuisance regressors (drift, motion and motion derivatives) and input into the iterative denoising algorithm. (B) The 2nd-order Fourier design matrices (constructed as in Equation S.15) separated out and sorted in cardiac phase *ϕ* order for visualization.

## 8 Acknowledgments

The authors would like to thank Maria Bordyug, Lauren C. Edwards, Alin Alshaheri Durazo, Amanda I. Gonzalez, and Mary Ellen Garcia for their assistance in collecting the data.

*   Received June 21, 2024.
*   Revision received June 21, 2024.
*   Accepted June 23, 2024.


*   © 2024, Posted by Cold Spring Harbor Laboratory

This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at [http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/)

## 9 References

1.  1.Duprez DA, De buyzere ML, Van den noortgate N, et al. Relationship between periventricular or deep white matter lesions and arterial elasticity indices in very old people. Age Ageing. 2001;30:325–330. doi: 10.1093/ageing/30.4.325.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/ageing/30.4.325&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=11509311&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

2.  2.Henskens LH, Kroon AA, Van oostenbrugge RJ, et al. Increased aortic pulse wave velocity is associated with silent cerebral small-vessel disease in hypertensive patients. Hypertension. 2008;52:1120–1126. doi: 10.1161/HYPERTENSIONAHA.108.119024.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1161/HYPERTENSIONAHA.108.119024&link_type=DOI) 

3.  3.Mitchell GF. Effects of central arterial aging on the structure and function of the peripheral vasculature: implications for end-organ damage. J Appl Physiol. 2008;105:1652–1660. doi: 10.1152/japplphysiol.90549.2008.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1152/japplphysiol.90549.2008&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18772322&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

4.  4.Mitchell GF, Van buchem MA, Sigurdsson S, et al. Arterial stiffness, pressure and flow pulsatility and brain structure and function: The Age, Gene/Environment Susceptibility-Reykjavik Study. Brain. 2011;134:3398–3407. doi: 10.1093/brain/awr253.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/brain/awr253&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22075523&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000296976500031&link_type=ISI) 

5.  5.Gorelick PB, Scuteri A, Black SE, et al. Vascular Contributions to Cognitive Impairment and Dementia: A Statement for Healthcare Professionals From the American Heart Association/American Stroke Association. Stroke. 2011;42:2672–2713. doi: 10.1161/STR.0b013e3182299496.
    
    [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6OToic3Ryb2tlYWhhIjtzOjU6InJlc2lkIjtzOjk6IjQyLzkvMjY3MiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDI0LzA2LzIzLzIwMjQuMDYuMjEuMjQzMDkyNjEuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 

6.  6.Tublin JM, Adelstein JM, Del monte F, Combs CK, Wold LE. Getting to the Heart of Alzheimer Disease. Circ Res. 2019;124:142–149. doi: 10.1161/CIRCRESAHA.118.313563.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1161/CIRCRESAHA.118.313563&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

7.  7.Zlokovic BV. Neurovascular pathways to neurodegeneration in Alzheimer’s disease and other disorders. Nat Rev Neurosci. 2011;12:723–738. doi: 10.1038/nrn3114.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nrn3114&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22048062&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

8.  8.Zarrinkoob L, Ambarki K, Wahlin A, et al. Aging alters the dampening of pulsatile blood flow in cerebral arteries. J Cerebr Blood F Met. 2016;36:1519–1527. doi: 10.1177/0271678X16629486.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/0271678X16629486&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26823470&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

9.  9.Pahlavian SH, Wang X, Ma S, et al. Cerebroarterial pulsatility and resistivity indices are associated with cognitive impairment and white matter hyperintensity in elderly subjects: A phase-contrast MRI study. J Cerebr Blood F Met. 2021;41:670–683. doi: 10.1177/0271678X20927101.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/0271678X20927101&link_type=DOI) 

10. 10.Heidari pahlavian S, Cen SY, Bi X, Wang DJ, Chui HC, Yan L. Assessment of carotid stiffness by measuring carotid pulse wave velocity using a single-slice oblique-sagittal phase-contrast MRI. Magn Reson Med. 2021;86:442–455. doi: 10.1002/mrm.28677.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.28677&link_type=DOI) 

11. 11.Rivera-rivera LA, Schubert T, Turski P, et al. Changes in intracranial venous blood flow and pulsatility in Alzheimer’s disease: A 4D flow MRI study. Journal of Cerebral Blood Flow & Metabolism. 2017;37:2149–2158. doi: 10.1177/0271678X16661340.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/0271678X16661340&link_type=DOI) 

12. 12.Lim JS, Lee JY, Kwon HM, Lee YS. The correlation between cerebral arterial pulsatility and cognitive dysfunction in Alzheimer’s disease patients. J Neurol Sci. 2017;373:285–288. doi: 10.1016/j.jns.2017.01.001.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.jns.2017.01.001&link_type=DOI) 

13. 13.Pahlavian SH, Jog M, Ma S, Wang DJJ, Yan L. Quantification of intracranial vascular compliance using multi-PLD pseudo-continuous arterial spin labeling with retrospective cardiac gating. In Proceedings of the 27th Annual Meeting of ISMRM. 2019;4950.
    
    

14. 14.Yan L, Liu CY, Smith RX, et al. Assessing intracranial vascular compliance using dynamic arterial spin labeling. NeuroImage. 2016;124:433–441. doi: 10.1016/j.neuroimage.2015.09.008.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2015.09.008&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26364865&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

15. 15.Warnert EA, Murphy K, Hall JE, Wise RG. Noninvasive assessment of arterial compliance of human cerebral arteries with short inversion time arterial spin labeling. J Cerebr Blood F Met. 2015;35:461–468. doi: 10.1038/jcbfm.2014.219.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/jcbfm.2014.219&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25515216&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

16. 16.Warnert EA, Verbree J, Wise RG, Van osch MJ. Using high-field magnetic resonance imaging to estimate distensibility of the middle cerebral artery. Neurodegener Dis. 2016;16:407–410. doi: 10.1159/000446397.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1159/000446397&link_type=DOI) 

17. 17.Li Y, Lim C, Schär M, et al. Three-dimensional assessment of brain arterial compliance: Technical development, comparison with aortic pulse wave velocity, and age effect. Magn Reson Med. 2021;86:1917–1928. doi: 10.1002/mrm.28835.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.28835&link_type=DOI) 

18. 18.Montagne A, Nation DA, Pa J, Sweeney MD, Toga AW, Zlokovic BV. Brain imaging of neurovascular dysfunction in Alzheimer’s disease. Acta Neuropathol. 2016;131:687–707. doi: 10.1007/s00401-016-1570-0.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00401-016-1570-0&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27038189&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

19. 19.Bouvy WH, Geurts LJ, Kuijf HJ, et al. Assessment of blood flow velocity and pulsatility in cerebral perforating arteries with 7-T quantitative flow MRI. NMR Biomed. 2016;29:1295–1304. doi: 10.1002/nbm.3306.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/nbm.3306&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25916399&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

20. 20.Geurts L, Biessels GJ, Luijten P, Zwanenburg J. Better and faster velocity pulsatility assessment in cerebral white matter perforating arteries with 7T quantitative flow MRI through improved slice profile, acquisition scheme, and postprocessing. Magn Reson Med. 2018;79:1473–1482. doi: 10.1002/mrm.26821.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.26821&link_type=DOI) 

21. 21.Geurts LJ, Zwanenburg JJ, Klijn CJ, Luijten PR, Biessels GJ. Higher Pulsatility in Cerebral Perforating Arteries in Patients with Small Vessel Disease Related Stroke, a 7T MRI Study. Stroke. 2019;50:62–68. doi: 10.1161/STROKEAHA.118.022516.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1161/STROKEAHA.118.022516&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30580730&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

22. 22.Tang J, Heidari pahlavian S, Joe E, et al. Assessment of arterial pulsatility of cerebral perforating arteries using 7T high-resolution dual-VENC phase-contract MRI. Magn Reson Med. 2024;92:605–617. doi: 10.1002/mrm.30073.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.30073&link_type=DOI) 

23. 23.Qin Q, Alsop DC, Bolar DS, et al. Velocity-selective arterial spin labeling perfusion MRI: A review of the state of the art and recommendations for clinical implementation. Magn Reson Med. 2022;88:1528–1547. doi: 10.1002/mrm.29371.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.29371&link_type=DOI) 

24. 24.Alsop DC, Detre JA, Golay X, et al. Recommended implementation of arterial spinlabeled perfusion MRI for clinical applications: a consensus of the ISMRM Perfusion Study Group and the European Consortium for ASL in Dementia. Magn Reson Med. 2015;73:102–116. doi: 10.1002/mrm.25197.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.25197&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24715426&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

25. 25.Piechnik SK, Chiarelli PA, Jezzard P. Modelling vascular reactivity to investigate the basis of the relationship between cerebral blood volume and flow under CO2 manipulation. NeuroImage. 2008;39:107–118. doi: 10.1016/j.neuroimage.2007.08.022.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2007.08.022&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17920935&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

26. 26.Franklin SL, Schmid S, Bos C,  Osch MJ van. Influence of the cardiac cycle on velocity selective and acceleration selective arterial spin labeling. Magn Reson Med. 2020;83:872–882. doi: 10.1002/mrm.27973.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.27973&link_type=DOI) 

27. 27. Wen-chau Wu, Edlow B, Elliot M,  Jiongjiong Wang, Detre J. Physiological Modulations in Arterial Spin Labeling Perfusion Magnetic Resonance Imaging. IEEE T Med Imaging. 2009;28:703–709. doi: 10.1109/TMI.2008.2012020.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1109/TMI.2008.2012020&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19150788&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

28. 28.Shi Y, Thrippleton MJ, Blair GW, et al. Small vessel disease is associated with altered cerebrovascular pulsatility but not resting cerebral blood flow. J Cerebr Blood F Met. 2020;40:85–99. doi: 10.1177/0271678X18803956.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/0271678X18803956&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=30295558&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

29. 29.1.   C Sommer
    
    Rodríguez A, Tembl J, Mesa-gresa P, Muñoz MÁ, Montoya P, Rey B. Altered cerebral blood flow velocity features in fibromyalgia patients in resting-state conditions. PLOS ONE. 2017;12. Ed. by  C Sommer: e0180253. doi: 10.1371/journal.pone.0180253.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0180253&link_type=DOI) 

30. 30.Guo J, Meakin JA, Jezzard P, Wong EC. An optimized design to reduce eddy current sensitivity in velocity-selective arterial spin labeling using symmetric BIR-8 pulses. Magn Reson Med. 2015;73:1085–1094. doi: 10.1002/mrm.25227.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.25227&link_type=DOI) 

31. 31.Cox R. AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res. 1996;29:162–173. doi: 10.1006/cbmr.1996.0014.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1006/cbmr.1996.0014&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8812068&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1996UV56700002&link_type=ISI) 

32. 32.Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM. FSL. NeuroImage. 2012;62:782–790. doi: 10.1016/j.neuroimage.2011.09.015.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2011.09.015&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21979382&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000306390600032&link_type=ISI) 

33. 33.Smith SM, Jenkinson M, Woolrich MW, et al. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage. 2004;23:S208–S219. doi: 10.1016/j.neuroimage.2004.07.051.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2004.07.051&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15501092&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000225374100020&link_type=ISI) 

34. 34.Power JD, Mitra A, Laumann TO, Snyder AZ, Schlaggar BL, Petersen SE. Methods to detect, characterize, and remove motion artifact in resting state fMRI. NeuroImage. 2014;84:320–341. doi: 10.1016/j.neuroimage.2013.08.048.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2013.08.048&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23994314&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000328868600030&link_type=ISI) 

35. 35.Winkler AM, Ridgway GR, Webster MA, Smith SM, Nichols TE. Permutation inference for the general linear model. NeuroImage. 2014;92:381–397. doi: 10.1016/j.neuroimage.2014.01.060.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2014.01.060&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24530839&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000335713000035&link_type=ISI) 

36. 36.Vikner T, Eklund A, Karalija N, et al. Cerebral arterial pulsatility is linked to hippocampal microvascular function and episodic memory in healthy older adults. J Cerebr Blood F Met. 2021;41:1778–1790. doi: 10.1177/0271678X20980652.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1177/0271678X20980652&link_type=DOI) 

37. 37.Guo F, Zhao C, Shou Q, Shao X, Wang DJ. Assessing Cerebral Microvascular Compliance with High-Resolution VASO MRI at 7T. In Proceedings of the 33rd Annual Meeting of ISMRM. 2024;1382.
    
    

38. 38.Wu WC, Mazaheri Y, Wong EC. The effects of flow dispersion and cardiac pulsation in arterial spin labeling. IEEE T Med Imaging. 2007;26:84–92. doi: 10.1109/TMI.2006.886807.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1109/TMI.2006.886807&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17243587&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

39. 39.Fushimi Y, Okada T, Yamamoto A, Kanagaki M, Fujimoto K, Togashi K. Timing dependence of peripheral pulse-wave-triggered pulsed arterial spin labeling. NMR Biomed. 2013;26:1527–1533. doi: 10.1002/nbm.2986.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/nbm.2986&link_type=DOI) 

40. 40.Li Y, Mao D, Li Z, et al. Cardiac-triggered pseudo-continuous arterial-spin-labeling: A cost-effective scheme to further enhance the reliability of arterial-spin-labeling MRI. Magn Reson Med. 2018;80:969–975. doi: 10.1002/mrm.27090.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.27090&link_type=DOI) 

41. 41.Verbree J,  Osch MJ van. Influence of the cardiac cycle on pCASL: cardiac triggering of the end-of-labeling. Magn Reson Mater Phy. 2018;31:223–233. doi: 10.1007/s10334-017-0611-6.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s10334-017-0611-6&link_type=DOI) 

42. 42.1.   JC De la torre
    
    De toledo ferraz alves TC, Ferreira LK, Wajngarten M, Busatto GF. Cardiac Disorders as Risk Factors for Alzheimer’s Disease. J Alzheimers Dis. 2010;20. Ed. by  JC De la torre: 749–763. doi: 10.3233/JAD-2010-091561.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3233/JAD-2010-091561&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20413875&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000279538600007&link_type=ISI) 

43. 43.Bolar DS, Rosen BR, Sorensen AG, Adalsteinsson E. QUantitative Imaging of eXtraction of oxygen and TIssue consumption (QUIXOTIC) using venular-targeted velocity-selective spin labeling. Magn Reson Med. 2011;66:1550–1562. doi: 10.1002/mrm.22946.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.22946&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21674615&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

44. 44.Guo J, Wong EC. Venous oxygenation mapping using velocity-selective excitation and arterial nulling. Magn Reson Med. 2012;68:1458–1471. doi: 10.1002/mrm.24145.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.24145&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22294414&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

45. 45.Wen Q, Wright A, Tong Y, et al. Paravascular fluid dynamics reveal arterial stiffness assessed using dynamic diffusion-weighted imaging. NMR Biomed. 2024;37:e5048. doi: 10.1002/nbm.5048.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/nbm.5048&link_type=DOI) 

46. 46.Wong EC. Time Efficient CSF Suppressed Velocity Selective ASL using a T 2-FLAIR Preparation. In Proceedings of the 12th Annual Meeting of ISMRM. 2004;711.
    
    

47. 47.Guo J. Optimizing background suppression for dual-module velocity-selective arterial spin labeling: Without using additional background-suppression pulses. Magn Reson Med. 2024;91:2320–2331. doi: 10.1002/mrm.29995.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.29995&link_type=DOI) 

48. 48.Qin Q,  Zijl PC van. Velocity-selective-inversion prepared arterial spin labeling. Magn Reson Med. 2016;76:1136–1148. doi: 10.1002/mrm.26010.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.26010&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26507471&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

49. 49.Liu D, Li W, Xu F, Zhu D, Shin T, Qin Q. Ensuring both velocity and spatial responses robust to B0/B1+ field inhomogeneities for velocity-selective arterial spin labeling through dynamic phase-cycling. Magn Reson Med. 2021;85:2723–2734. doi: 10.1002/mrm.28622.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.28622&link_type=DOI) 

50. 50.Bolar DS, Barnes RA, Chen C, Liu TT, Wong EC. Reduced B0/B1+ Sensitivity in Velocity-Selective Inversion ASL Using Adiabatic Refocusing Pulses. In Proceedings of the 33rd Annual Meeting of ISMRM. 2024.
    
    

51. 51.Kuethe DO, Caprihan A, Gach HM, Lowe IJ, Fukushima E. Imaging obstructed ventilation with NMR using inert fluorinated gases. J Appl Physiol. 2000;88:2279–2286. doi: 10.1152/jappl.2000.88.6.2279.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1152/jappl.2000.88.6.2279&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10846046&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000087603200045&link_type=ISI) 

52. 52.Díaz-francés E, Rubio FJ. On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables. Stat Pap. 2013;54:309–323. doi: 10.1007/s00362-012-0429-2.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00362-012-0429-2&link_type=DOI) 

## S6 References

1.  1.1.   C Sommer
    
    Rodríguez A, Tembl J, Mesa-gresa P, Muñoz MÁ, Montoya P, Rey B. Altered cerebral blood flow velocity features in fibromyalgia patients in resting-state conditions. PLOS ONE. 2017;12. Ed. by  C Sommer: e0180253. doi: 10.1371/journal.pone.0180253.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0180253&link_type=DOI) 

2.  2.Wen-chau Wu, Edlow B, Elliot M,  Jiongjiong Wang, Detre J. Physiological Modulations in Arterial Spin Labeling Perfusion Magnetic Resonance Imaging. IEEE T Med Imaging. 2009;28:703–709. doi: 10.1109/TMI.2008.2012020.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1109/TMI.2008.2012020&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19150788&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 

3.  3.Gudbjartsson H, Patz S. The rician distribution of noisy mri data. Magn Reson Med. 1995;34:910–914. doi: 10.1002/mrm.1910340618.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/mrm.1910340618&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8598820&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1995TH39600017&link_type=ISI) 

4.  4.Kuethe DO, Caprihan A, Gach HM, Lowe IJ, Fukushima E. Imaging obstructed ventilation with NMR using inert fluorinated gases. J Appl Physiol. 2000;88:2279–2286. doi: 10.1152/jappl.2000.88.6.2279.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1152/jappl.2000.88.6.2279&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10846046&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000087603200045&link_type=ISI) 

5.  5.Díaz-francés E, Rubio FJ. On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables. Stat Pap. 2013;54:309–323. doi: 10.1007/s00362-012-0429-2.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00362-012-0429-2&link_type=DOI) 

6.  6.Power JD, Mitra A, Laumann TO, Snyder AZ, Schlaggar BL, Petersen SE. Methods to detect, characterize, and remove motion artifact in resting state fMRI. NeuroImage. 2014;84:320–341. doi: 10.1016/j.neuroimage.2013.08.048.
    
    [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.neuroimage.2013.08.048&link_type=DOI) 
    
    [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23994314&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F23%2F2024.06.21.24309261.atom) 
    
    [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000328868600030&link_type=ISI)

 [1]: /embed/graphic-2.gif
 [2]: /embed/graphic-3.gif
 [3]: /embed/graphic-4.gif
 [4]: /embed/inline-graphic-1.gif
 [5]: /embed/graphic-5.gif
 [6]: /embed/inline-graphic-2.gif
 [7]: /embed/graphic-6.gif
 [8]: /embed/inline-graphic-3.gif
 [9]: /embed/graphic-7.gif
 [10]: /embed/graphic-8.gif
 [11]: /embed/inline-graphic-4.gif
 [12]: /embed/graphic-12.gif
 [13]: /embed/inline-graphic-5.gif
 [14]: /embed/inline-graphic-6.gif
 [15]: /embed/graphic-13.gif
 [16]: /embed/inline-graphic-7.gif
 [17]: /embed/graphic-14.gif
 [18]: /embed/graphic-15.gif
 [19]: /embed/inline-graphic-8.gif
 [20]: /embed/inline-graphic-9.gif
 [21]: /embed/graphic-21.gif
 [22]: /embed/graphic-22.gif
 [23]: /embed/graphic-23.gif
 [24]: /embed/graphic-24.gif
 [25]: /embed/graphic-25.gif
 [26]: /embed/graphic-26.gif
 [27]: /embed/graphic-27.gif
 [28]: /embed/graphic-28.gif
 [29]: /embed/graphic-29.gif
 [30]: /embed/graphic-30.gif
 [31]: /embed/graphic-31.gif
 [32]: /embed/inline-graphic-10.gif
 [33]: /embed/inline-graphic-11.gif
 [34]: /embed/inline-graphic-12.gif
 [35]: /embed/inline-graphic-13.gif
 [36]: /embed/graphic-32.gif
 [37]: /embed/graphic-33.gif
 [38]: /embed/graphic-34.gif
 [39]: /embed/graphic-35.gif
 [40]: /embed/graphic-36.gif
 [41]: /embed/inline-graphic-14.gif
 [42]: /embed/inline-graphic-15.gif
 [43]: /embed/graphic-37.gif
 [44]: /embed/inline-graphic-16.gif
 [45]: /embed/graphic-38.gif
 [46]: /embed/inline-graphic-17.gif
 [47]: /embed/graphic-39.gif
 [48]: /embed/graphic-40.gif
 [49]: /embed/inline-graphic-18.gif
 [50]: /embed/inline-graphic-19.gif
 [51]: /embed/graphic-41.gif
 [52]: /embed/inline-graphic-20.gif
 [53]: /embed/inline-graphic-21.gif
 [54]: /embed/graphic-42.gif
 [55]: /embed/graphic-43.gif
 [56]: /embed/inline-graphic-22.gif
 [57]: /embed/inline-graphic-23.gif
 [58]: /embed/inline-graphic-24.gif
 [59]: /embed/graphic-44.gif
 [60]: /embed/inline-graphic-25.gif
 [61]: /embed/graphic-45.gif
 [62]: /embed/inline-graphic-26.gif
 [63]: /embed/graphic-46.gif
 [64]: /embed/inline-graphic-27.gif
 [65]: /embed/graphic-47.gif
 [66]: /embed/graphic-48.gif
 [67]: /embed/graphic-49.gif
 [68]: /embed/graphic-50.gif
 [69]: /embed/inline-graphic-28.gif
 [70]: /embed/inline-graphic-29.gif
 [71]: /embed/inline-graphic-30.gif
 [72]: /embed/graphic-51.gif
 [73]: /embed/graphic-52.gif
 [74]: /embed/inline-graphic-31.gif
 [75]: /embed/graphic-53.gif
 [76]: /embed/inline-graphic-32.gif
 [77]: /embed/graphic-54.gif
 [78]: /embed/inline-graphic-33.gif
 [79]: /embed/graphic-55.gif
 [80]: /embed/inline-graphic-34.gif
 [81]: /embed/inline-graphic-35.gif
 [82]: /embed/inline-graphic-36.gif
 [83]: /embed/graphic-56.gif
 [84]: /embed/inline-graphic-37.gif
 [85]: /embed/inline-graphic-38.gif
 [86]: /embed/graphic-57.gif
 [87]: /embed/inline-graphic-39.gif
 [88]: /embed/inline-graphic-40.gif
 [89]: /embed/graphic-58.gif
 [90]: /embed/graphic-59.gif
 [91]: /embed/inline-graphic-41.gif
 [92]: /embed/graphic-60.gif
 [93]: /embed/inline-graphic-42.gif
 [94]: /embed/inline-graphic-43.gif
 [95]: /embed/graphic-61.gif
 [96]: /embed/inline-graphic-44.gif
 [97]: /embed/inline-graphic-45.gif
 [98]: /embed/inline-graphic-46.gif
 [99]: /embed/inline-graphic-47.gif
 [100]: /embed/graphic-62.gif
 [101]: /embed/graphic-63.gif
 [102]: /embed/graphic-64.gif
 [103]: /embed/graphic-65.gif
 [104]: /embed/graphic-66.gif
 [105]: /embed/graphic-67.gif
 [106]: /embed/graphic-68.gif
 [107]: /embed/inline-graphic-48.gif
 [108]: /embed/graphic-69.gif
 [109]: /embed/inline-graphic-49.gif
 [110]: /embed/inline-graphic-50.gif
 [111]: /embed/graphic-70.gif
 [112]: /embed/inline-graphic-51.gif
 [113]: /embed/inline-graphic-52.gif
 [114]: /embed/inline-graphic-53.gif
 [115]: /embed/inline-graphic-54.gif
 [116]: /embed/inline-graphic-55.gif
 [117]: /embed/inline-graphic-56.gif
 [118]: /embed/inline-graphic-57.gif
 [119]: /embed/inline-graphic-58.gif
 [120]: /embed/inline-graphic-59.gif
 [121]: /embed/inline-graphic-60.gif
 [122]: /embed/graphic-71.gif
 [123]: /embed/graphic-72.gif
 [124]: /embed/inline-graphic-61.gif
 [125]: /embed/graphic-73.gif
 [126]: /embed/inline-graphic-62.gif
 [127]: /embed/inline-graphic-63.gif
 [128]: /embed/inline-graphic-64.gif
 [129]: /embed/inline-graphic-65.gif
 [130]: /embed/inline-graphic-66.gif
 [131]: /embed/inline-graphic-67.gif
 [132]: /embed/inline-graphic-68.gif
 [133]: /embed/inline-graphic-69.gif
 [134]: /embed/inline-graphic-70.gif
 [135]: /embed/inline-graphic-71.gif
 [136]: /embed/inline-graphic-72.gif
 [137]: /embed/inline-graphic-73.gif
 [138]: /embed/inline-graphic-74.gif
 [139]: /embed/inline-graphic-75.gif
 [140]: /embed/inline-graphic-76.gif
 [141]: /embed/graphic-74.gif
 [142]: /embed/graphic-75.gif
 [143]: /embed/graphic-76.gif
 [144]: /embed/inline-graphic-77.gif
 [145]: /embed/inline-graphic-78.gif
 [146]: /embed/inline-graphic-79.gif
 [147]: /embed/inline-graphic-80.gif
 [148]: /embed/inline-graphic-81.gif
 [149]: /embed/inline-graphic-82.gif
 [150]: /embed/inline-graphic-83.gif
 [151]: /embed/inline-graphic-84.gif
 [152]: /embed/inline-graphic-85.gif
 [153]: /embed/inline-graphic-86.gif
 [154]: /embed/inline-graphic-87.gif
 [155]: /embed/inline-graphic-88.gif
 [156]: /embed/inline-graphic-89.gif
 [157]: /embed/inline-graphic-90.gif
 [158]: /embed/inline-graphic-91.gif
 [159]: F10/embed/inline-graphic-92.gif
 [160]: F10/embed/inline-graphic-93.gif
 [161]: F11/embed/inline-graphic-94.gif
 [162]: F11/embed/inline-graphic-95.gif
 [163]: F11/embed/inline-graphic-96.gif
 [164]: F11/embed/inline-graphic-97.gif
 [165]: F11/embed/inline-graphic-98.gif
 [166]: F11/embed/inline-graphic-99.gif
 [167]: F12/embed/inline-graphic-100.gif
 [168]: F12/embed/inline-graphic-101.gif
 [169]: F12/embed/inline-graphic-102.gif
 [170]: F12/embed/inline-graphic-103.gif
 [171]: F12/embed/inline-graphic-104.gif
 [172]: F12/embed/inline-graphic-105.gif
 [173]: F12/embed/inline-graphic-106.gif