Cardiac-induced cerebral pulsatility, brain structure, and cognition in middle and older-aged adults.

Changes of cardiac-induced regional pulsatility can be associated with specific regions of brain volumetric changes, and these are related with cognitive alterations. Thus, mapping of cardiac pulsatility over the entire brain can be helpful to assess these relationships. A total of 108 subjects (age: 66.5 ± 8.4 years, 68 females, 52 healthy controls, 11 subjective cognitive decline, 17 impaired without complaints, 19 MCI and 9 AD) participated. The pulsatility map was obtained directly from resting-state functional MRI time-series data at 3T. Regional brain volumes were segmented from anatomical MRI. Multidomain neuropsychological battery was performed to test memory, language, attention and visuospatial construction. The Montreal Cognitive Assessment (MoCA) was also administered. The sparse partial least square (SPLS) method, which is desirable for better interpreting high-dimensional variables, was applied for the relationship between the entire brain voxels of pulsatility and 45 segmented brain volumes. A multiple holdout SPLS framework was used to optimize sparsity for assessing the pulsatility-volume relationship model and to test the reliability by fitting the models to 9 different splits of the data. We found statistically significant associations between subsets of pulsatility voxels and subsets of segmented brain volumes by rejecting the omnibus null hypothesis (any of 9 splits has p < 0.0056 (=0.05/9) with the Bonferroni correction). The pulsatility was positively associated with the lateral ventricle, choroid plexus, inferior lateral ventricle, and 3rd ventricle and negatively associated with hippocampus, ventral DC, and thalamus volumes for the first pulsatility-volume relationship. The pulsatility had an additional negative relationship with the amygdala and brain stem volumes for the second pulsatility-volume relationship. The spatial distribution of correlated pulsatility was observed in major feeding arteries to the brain regions, ventricles, and sagittal sinus. The indirect mediating pathways through the volumetric changes were statistically significant between the pulsatility and multiple cognitive measures (p < 0.01). Thus, the cerebral pulsatility, along with volumetric measurements, could be a potential marker for better understanding of pathophysiology and monitoring disease progression in age-related neurodegenerative disorders.

Introduction

Aging and aging-related diseases such as Alzheimer’s Disease (AD) have been previously shown to be associated with reduced the tissue volume in specific brain regions and enlargement of the brain ventricles (Apostolova et al., 2012; Coupe et al., 2019). These volume losses, as measured by magnetic resonance imaging (MRI), can serve as important early imaging biomarkers to identify risk for AD, especially with hippocampal atrophy being a key supportive feature for the early stage of the disease (Dubois et al., 2007; Jack et al., 1992; Laakso et al., 1995; Schuff et al., 2009). However, potential mediators for the alteration in brain volume remain poorly understood.

The brain is perfused by high blood flow throughout the cardiac pulsation to maintain its metabolism. With aging, an increase in arterial stiffness and a decrease in arterial compliance are associated with increases in cerebral pulsatility (Mitchell et al., 2011; Townsend et al., 2015). The persistent exposure of small vessels to such high pulsation results in cerebral microvascular damage that can lead to structural and functional changes in the brain (Mitchell et al., 2011). Thus, it is important to measure the cerebral pulsatility and its relation to brain volume. The pulsatility has been assessed by flow-driven approaches, e.g. phase contrast MRI or transcranial Doppler ultrasound (Naqvi et al., 2013; Wahlin et al., 2014), which, however, can only measure limited locations of pulsatility.

The cardiac pulse generated by the heartbeat is delivered to large vessels in the brain, as well as the vascular trees, and ventricles. Each brain region has locally different arterial supplies and venous drainages. Thus, transmission of the cardiac pulsation to the brain may change while passing through different brain regions, and its effect on the brain tissue may be locally different. Thus, some brain regions may be more affected by the pulsatility changes. In other words, all voxels of pulsatility in the brain may not be associated with entire brain structural changes, i.e., specific regions of cerebral pulsatility, may be associated with specific regional brain volumes. Therefore, imaging pulsatilty for the whole brain is necessary to evaluate these regional relationships.

The cardiac pulsation and its variability contribute to the signal intensity fluctuations in the resting-state blood oxygen level-dependent (BOLD) signal (Tsvetanov et al., 2020). Thus, the amplitude of the cardiac-induced temporal fluctuation in resting-state functional MRI (rs-fMRI) time-series can be expressed as the cerebral pulsatility. Identification of cardiac-induced signal changes, i.e., cerebral pulsatility, has been successfully obtained from fMRI time-series by using the phases of cardiac cycle at fMRI time-series with the aid of an external device (e.g. pulse-oximeter) (Dagli et al., 1999; Glover et al., 2000). However, in our study, older adults often were not cooperative with the pulse oximeter and many failures arose to get the cardiac cycle. Thus, we developed a method to directly estimate the phase of cardiac cycle from the rs-fMRI data itself, rather than relying on an external device, for the calculation of the cerebral pulsatility. rs-fMRI is commonly included in most neuroimaging studies as it is easy to implement. Therefore, this method can be readily used to easily assess the pulsatility for the entire brain, especially when pulse-oximeter acquisition is missing or has failed.

In this study, we tested the relationship between the cerebral pulsatility obtained from rs-fMRI and the brain volumetric measures from anatomical images, and identified the brain regions that are mostly associated with the pulsatility-volume relationship using the sparse partial least squares (SPLS) method. The sparsity automatically selects important features and reduces the risk of detecting spurious associations. Thus, it is desirable for better interpretation of high-dimensional variables, i.e., many voxels of pulsatility and multiple segmented brain volumes. A multiple holdout framework was used to optimize the regularization parameters for the sparsity and to improve the generalizability of the optimized SPLS models (Monteiro et al., 2016). Then, we tested significance in the pathway of observed volumetric measures between the pulsatility and multiple cognitive measures by mediation analysis.

Materials and methods

Participants

Data from 108 participants who had completed preprocessing in June 2019 from Connectomics in Brain Aging and Dementia, a Connectome Related to Human Disease study (https://www.humanconnectome.org/disease-studies), was used for this analysis. This represents almost a half of our targeted recruitment for the current ongoing study. The study was approved by the University of Pittsburgh Institutional Review Board and written informed consent was obtained from all participants. The individuals enrolled in the study did not have a self-reported history of major central nervous system pathology (e.g., epilepsy, tumor, Huntington’s disease, Parkinson’s disease, multiple sclerosis, stroke, head injury with loss of consciousness for more than 30 min), major psychiatric disease (e.g., severe major depression, bipolar disorder, substance abuse), or active cancer. We also excluded persons who had sensory or motor deficits sufficient to impair their ability to perform the core neuropsychological tasks. Individuals taking psychoactive medications at doses that could affect the neuropsychological testing or brain functional responses such as benzodiazepines (e.g., lorazepam) or narcotic analgesics (e.g., acetaminophen/propoxyphene) were likewise excluded from the study.

Participants completed a multidomain neuropsychological (NP) battery, including tests of memory (Consortium to Establish a Registry of Alzheimer’s Disease Word List Learning test; modified Rey Osterrieth figure recall); language (30-item Boston Naming Test; verbal fluency); attention (Trail Making Test A; digit span forward); executive functions (Trail Making Test B; digit span backward; clock drawing); and visuospatial construction (modified Block Design; Rey Osterrieth figure copy). The Montreal Cognitive Assessment (MoCA), a global screening measure, was also administered.

Cognitive classifications were made independently by two neuropsychologists (JTB and BES) and any differences were resolved by group discussion. Cognitive status categories included normal cognition (i.e., NP test performance within normal limits, no clinically significant complaints); subjective cognitive complaints (NP within normal limits, clinically significant complaints of cognitive decline); NP impaired without complaints, and mild cognitive impairment (MCI; NP impaired with clinically significant complaints of cognitive decline). Detailed diagnostic criteria have been previously reported (Snitz et al., 2015, 2018). Participants with AD were recruited from the University of Pittsburgh Alzheimer Disease Research Center (ADRC), having completed comprehensive diagnostic evaluation including a multidisciplinary diagnostic case conference (Lopez et al., 2005).

The characteristics of the study subjects included in this analysis are shown in Table 1. Additional information on race and medical history is listed in Supplementary Table S2.

Table 1

Study subject characteristics.

Category	# of subjects	Age	Sex (Fe/male)	Handedness (L/R)
Healthy control	52	67.6 ± 7.7	37/15	1/51
Subjective cognitive decline	11	68.0 ± 7.1	6/5	1/10
Impaired without complaints	17	60.4 ± 7.3	12/5	4/13
MCI	19	68.4 ± 11.6	11/8	1/17
AD	9	66.2 ± 9.2	2/7	3/6

MRI data acquisition

All subjects were studied with a 3T MRI system (Siemens Prisma, Erlangen, Germany) using a 64-channel head coil. MRI data was obtained with the HCP imaging protocols (http://protocols.humanconnectome.org/CCF) including T1- and T2-weighted structural scans and eight separate runs of rs-fMRI (four runs each on day1 and day2). Anatomical images were acquired using a 3D T1-weighted magnetization prepared rapid gradient echo (MPRAGE; TR = 2.4 s, TE= 2.22 ms, TI = 1000 ms, FA = 8, voxel size = 0.8 × 0.8 × 0.8 mm3), and a 3D T2-weighted sampling perfection with application optimized contrast using different flip angle evolutions image (SPACE; TR = 3.2 s, TE = 563 ms, voxel size = 0.8 × 0.8 × 0.8 mm3). For both sequences, integrated parallel acquisition (iPAT) reconstruction was performed by the generalized auto-calibrating partially parallel acquisition (GRAPPA) algorithm with an acceleration factor of 2. The rs-fMRI data was obtained with a 2D single-shot gradient-echo echo-planar imaging (GE EPI; TR = 800 ms, TE = 37 ms, FOV = 208 mm, voxel size = 2 × 2 × 2 mm3, FA = 52°, and 72 slices with 420 volumes). Two spin-echo EPI acquisitions were acquired with opposite phase-encoding polarities (anterior-posterior and posterior-anterior) to allow for correction of the EPI geometric distortion. The simultaneous multislice (SMS) imaging technique was applied with a multiband composite RF pulse using an acceleration factor of eight. The cardiac cycle and respirations were monitored using the scanner’s built-in pulse oximeter placed on the subject’s index finger and a respiratory belt placed around the chest. The scanner provided a transistor-transistor logic (TTL) signal output for each slice to synchronize physiology and MRI data.

Data processing

The MRI data was processed using AFNI (http://afni.nimh.nih.gov/afni/), Freesurfer (https://surfer.nmr.mgh.harvard.edu/), and FSL (http://fsl.fmrib.ox.ac.uk/fsl) programs combined with an in-house Matlab program (MathWorks, Natick, MA).

Brain volume segmentation

The volume segmentation for each brain region was performed with Freesurfer 6.0 after alignment of T1- and T2-weighted anatomical images, bias field correction, and removal of gradient nonlinearity and readout distortion (Glasser et al., 2013). A total of 45 brain regions, summing the left and right hemispheres, were segmented from T1-weighted anatomical MPRGAE image. All segmented regional volumes were corrected by total intracranial volume (ICV; corrected volume = volume × 100/ICV).

Technique for cardiac-induced cerebral pulsatility

The main purpose of this technique is to measure the cardiac phase (θ) from a rs-fMRI time series without using an external measurement of cardiac cycle (e.g. pulse oximeter). Next, the pre-existing method that fits fMRI time-series of each voxel to the Fourier series with the phase was used to identify cardiac-induced pulsation (Glover et al., 2000).

Motion correction was performed by rigid-body registration of every image volume of the time-series to the SMS reference volume. But, because cardiac phase for each MRI acquisition is affected by severe motion and is difficult to correct, data with more than 4 mm (more than two voxel size) of motion were excluded.

Each individual’s main cardiac coupling areas can be obtained by transformation from an MNI template of the pulsatility map to each individual’s original space. To create a template for the pulsatility map, each individual’s time-series of cardiac phase was obtained by RETROICOR with successful measurement of pulse oximeter data, randomly selected from 32 subjects. The cerebral pulsatility map of each individual was calculated using Eq. (1) (see below) and then transformed to the MNI space and averaged. In this case, the cardiac phase (θ) was estimated by assigning the timing of each acquired slice in a single volume to one systole and the subsequent systole cardiac cycle (Dagli et al., 1999; Glover et al., 2000). The cardiac phase is defined as θ = 2⋅π⋅(t − ta)/(tb−ta), where t is MRI acquisition point (blue bars in Fig 1A), ta is the time of the R-wave peak in the cardiac cycle preceding t, and tb is the time for the subsequent R-wave peak (green arrows in Fig. 1A).

Fig. 1.

The cardiac phase measurement. (a) With an external measurement of cardiac cycle. The cardiac pulses synchronized with the scanner triggers determine the phase of the cardiac cycle at the time that each slice is acquired. Green arrows indicate R–R interval and blue bars show the acquisition of each image. Filled red circles are the calculated phase of the first slice of each volume and open circles are the phase for other slices. (b) Without an external measurement of cardiac cycle. The signal intensity changes in the main vessels can be followed by cardiac cycles (black sinusoid). The Hilbert transform generates signal S(t) that is constructed from a real-valued cardiac-induced sinusoidal input signal; S(t) = s(t) + i* h(t), where, S(t) is the analytic signal constructed from s(t) and its Hilbert transform, s(t) is the input signal, and h(t) is the Hilbert Transform of the input signal. The cardiac phase was obtained from the angle of the complex (real and imaginary) value for each MRI acquisition (red circles).

In order to measure the cardiac phase (θ) without an external measurement of cardiac cycle, we used the main component of the signal intensity changes in major vessel areas where it can mainly express the cardiac pulsation. To derive the main cardiac-induced sinusoidal signal in the brain, voxels above a threshold > 3.3σ of the template (>99.9%) were selected from the MNI template of the pulsatility map and generated a binary mask. This mask was transformed to each individual’s native space to determine each individual’s high cardiac-induced pulsatile regions that include main large vessels (see Fig. 2A). The signal intensity changes in the main arteries follow the cardiac cycle. Thus, a principal component analysis can be used to derive the main component of the cardiac cycle at each data acquisition point (i.e. the slices in the same SMS acquisition) after detrending the signal and removing motion. Since our TR is slightly less than the cardiac cycle, the phases per TR can correspond to the same pulsation. Thus, the main component of pulsation was estimated from all slice acquisitions because more signals provide better estimation of the main component. The amplitude of this sinusoidal time-series was normalized to −1 to 1.

Fig. 2.

Multiple holdout framework. (a) Regularization parameter optimization framework. (b) Permutation for statistical significance. (1) Nine spilt datasets were generated. For each spilt data, (2) the regularization parameters were optimized with 100 further randomly split data of trainingsub and validation from the training dataset. The subscript sub was used to distinguish from the name “training” previously used in (1). The correlation between the projection of validation data onto the weighting vectors by SPLS from trainingsub was obtained for 40 × 40 grid search. (3) Grids of correlation values from 100 subsamples were averaged. The highest correlation was selected from the grid search as the optimal constraint parameters. (4) The weighting vectors were obtained by SPLS from the training dataset with optimal regularization parameters, and the correlation between the pulsatility and volume components was calculated by projecting the weighting vectors onto the holdout dataset. (5) The same correlation procedure was applied to 10,000 randomly permuted null datasets. (6) p-value was calculated by the fraction of permutation values obtained in (5) that are at least as extreme as the original value obtained in (4).

In order to calculate the cardiac phase (θ) from the sinusoidal time-series, we used the angle of complex at each data acquisition point of sinusoidal time-series (red circles in Fig. 1B). That is, the phase for each bin in the FFT is the same as the relative phase shift of the sinusoid that represents the bin in the time domain. The Hilbert transform is performed with a real-valued cardiac-induced sinusoidal input signal and returns a complex (real and imaginary) result of the same length. The phase was calculated by the angle of the complex data at each volume of the time-series. Note that the Hilbert transform method carries the assumption of a monocomponent signal in order to obtain instantaneous phase.

The model of cardiac pulsation, the 2nd-order Fourier series of the cardiac phase (θ) time courses, was fit to each voxel of the imaging data with 6 motion parameters as nuisance variables. Then, the cerebral pulsatility map was calculated from the coefficients (ai and bi) of the linear regression divided by its standard error (SE, i.e., temporal variance that is not explained by the model of cardiac pulsation at each voxel) (Beall, 2010), using

To determine a meaningful cardiac pulsatile map, a null statistic distribution was created by a Monte Carlo technique of regressing random phase data for an arbitrary noise model (Beall, 2010). The null statistic distribution was used to find the threshold for 3σ significance, where the meaningful cardiac pulsatile signal is larger than 99.73% of the null couplings distribution.

The test-retest reliability of the pulsatile maps between two separate days was measured by the image intra-class correlation coefficient, which expresses a global measure of reliability for high-dimensional multivariate imaging data (Shou et al., 2013). The voxels within the brain mask were evaluated. A 95% confidence interval was calculated by 5000 bootstrapping repetitions. This method was also used to compare the pulsatility maps obtained by pulse-oximeter and by our approach. All measurements obtained from two separate days were averaged for each subject.

Relationship between cerebral pulsatility and brain volume

Each subject’s pulsatility map was transformed to the MNI space and voxels within the brain mask were spatially normalized as Z-scores so that each individual’s map was on the same scale. The cardiac pulsatile data were then generated as a two-dimensional matrix (subjects × voxels in the brain = 108 × 258,370) and assigned to independent variables (X). The multiple segmented brain volumes from Freesurfer were used as dependent variables (Y). Age, sex and handedness were included as covariates.

SPLS was performed to determine relationships between voxel-wise cerebral pulsatility (X) and multiple segmented brain volume (Y) to overcome PLS’s difficulty of interpreting the high-dimensional model (Mihalik et al., 2020; Monteiro et al., 2016). SPLS provides a simpler model by automatically selecting important features (most relevant variables) in the model through sparse weighting vector pairs. The sparsity was adjusted by l1-norm constraints. Thus, the model does not include irrelevant features by forcing the coefficients of those variables to be zero. This approach can attain statistically significant associations between subsets of brain voxels and subsets of brain volumes. A multiple holdout framework was used to optimize the regularization parameters and to test the generalizability of the optimized SPLS models to reduce computational burden of a regular nested cross-validation (Monteiro et al., 2016). The data was divided as training, validation and holdout (test) sets. The procedures are summarized in Fig. 2. Specifically, (1) 9-fold cross-validation (i.e. 9 datasets) was performed, where data were randomly divided into 9 subsets, and each fold was set aside as a holdout dataset (=12 subjects). (2) For each split (cross-validation) dataset, while excluding the holdout dataset, the remaining data (referred to as training data) were further randomly split 100 times to two subsamples, training (75% of the data = 72 subjects, referred to as trainingsub) and validation (25% of the data = 24 subjects), without overlapping and replacement. The regulation parameters were estimated and validated using the training and validation subsamples in each of the 100 subsamples. (3) A grid search, 40 equidistance points in and (Witten et al., 2009), which makes a total of 1600 regulation parameter combinations, was used to find the best regulation parameter set for optimization of the sparsity, where c and c are the regulation parameters and p and q are the numbers of features for X (i.e. the number of brain voxels) and Y (i.e. the number of brain volumes and covariates), respectively. For each pair of (c, c), the correlation between the projection of validation data onto the weighting vectors was obtained as ρ = |corr(Xu(−, Yv(−)|, where X and Y are the validation datasets and u(− and v(− are the pair of weighting vectors that calculated from SPLS using the trainingsub datasets for each c and c,. Then, the average correlation from the 100 subsamples was calculated . The optimal constraints of (c, c) that have the highest correlation, were selected from the 40 × 40 grid of correlations. (4) The optimal regularization parameters were used to fit SPLS for the training data, and the weighting vectors for pulsatility (u) and brain volumes (v) were acquired. Then, the holdout data were projected onto these weighting vectors, and the absolute correlation between two scores was calculated as ρ = |corr(Xu, Yv)|, where X and Y are the holdout dataset and u and v are the weighting vectors obtained by SPLS from the training dataset. (5) The statistical significance of ρ, as expressed in a P-value, was assessed by the following permutation test (Nichols and Holmes, 2001). 10,000 permutation datasets were generated (excluding the holdout data) by randomly changing the order of rows (subjects) in pulsatility data while leaving the brain volumes data unchanged (Monteiro et al., 2016). This permutation process breaks the relationship between pulsatility (X) and brain volumes (Y). The weighting vectors (upm and vpm) were acquired by SPLS from each permuted dataset, and correlations were calculated by ρpm = |corr(Xupm, Yvpm)|. (6) The statistical significance of the dependence between the two weighting vectors is tested under the null hypothesis, i.e., no relationship between the two sets of variables. A P-value was calculated by counting the number of ρpm that were greater than the observed ρ (i.e. estimated correlation coefficient from the non-permuted dataset) and dividing it by the number of permutations. A pseudo-count (= 1) was added to the numerator and denominator to avoid P-value being zero (i.e. P > 0), which arises when the test statistic is never surpassed by any of the permutation correlations. Therefore, the resulting smallest achievable P-value is ~1/number of permutations (Monteiro et al., 2016; Nichols and Holmes, 2001). A low P-value indicates that the null hypothesis is rejected, and the finding is statistically significant. Since the P-value is obtained by a small sample of holdout, the results may be variable depending on how each spilt is divided. All these processes were repeated for 9 cross-validated datasets. Cross-validation is a resampling procedure used to evaluate the model on a limited data sample. It estimates how the model (from the training data) is expected to perform in general by predicting holdout data (considered as unseen data). The cross-validation provides an equal opportunity for all data to be the holdout data. It generally results in a less biased estimate of the model than a simple training/test split. The concept of omnibus hypothesis is used to evaluate the significance of the test. Omnibus hypothesis is that the null hypotheses (i.e. there is no significant correlation) of all spilt data are true. In other words, if any of split data is statistically significant, the omnibus hypothesis is rejected. The Bonferroni correction was performed for multiple comparisons. Thus, if any of the 9 splits has P-value < 0.0056 (= 0.05/9), the omnibus hypothesis is rejected, and the effect found by SPLS is considered statistically significant.

For each split data, the brain regions of pulsatility that associated with the brain volume were identified by the voxels that satisfied p-value < 0.01 for 10,000 permuted data of the pulsatility weighting vectors.

Moreover, in order to observe additional maximum variances (the second weighting vectors) from a dataset, the projection deflation method, which removes the effect explained by the first weighting vectors (u and v) from the dataset, was used (Monteiro et al., 2016). The permutation was repeated for this projection deflated dataset to obtain P-value.

In addition, a single permutation test was performed to evaluate the average association across the 9 splits. That is, the whole data (no split) was projected onto the averaged weighting vectors from the 9 splits, and then the absolute correlation between X and Y scores was calculated. 10,000 permutations were performed with the whole data in the same manner as described above, using the average weighting vectors. Then, a P-value was calculated.

Mediation analysis

To evaluate the indirect effect of cerebral pulsatility on cognitive score through the pathway of brain volume, mediation analysis was performed (Preacher and Hayes, 2008) using the BROVO toolbox (https://sites.google.com/site/bravotoolbox). The statistical significances of the direct and indirect pathways were determined by 10,000 permutation iterations. The latent variables (L and L, which are also called scores of X and Y) were computed by projecting the pulsatility (X) and brain volume (Y) dataset onto their weighting vectors by SPLS with averaged optimal constraints (c and c), as L = Xu and L = Yv. L was treated as independent variable and each cognitive measure was assigned to dependent variable, and L was used as a mediator. Age, sex and handedness were included as nuisance covariates. Since SPLS maximizes the covariance between two datasets, it could be biased for their relationship. Thus, PCA was also applied. Unlike SPLS, the scores of the PCA best characterizes each dataset individually. The score of the pulsatility data in PCA was obtained from a subset of pulsatility, which is any voxel satisfying p < 0.01. The score of brain volume was acquired from a subset of brain volume that was identified in the SPLS analysis and assigned as a mediator.

Results

Cerebral pulsatility map

Of the total 864 MRI runs analyzed, 88 had severe motion and 6 were stopped by the subjects during the scan (11%). In 254 of the runs the pulse oximeter was either incompletely or poorly measured (29% of data). Thus, our method was used to obtain the cardiac phase for the pulsatility calculation, without relying on external cardiac cycle data.

Fig. 3A shows the ROI for the very high pulsatility region that was transformed from the MNI template to each individual’s space and overlaid on the EPI. Most of the regions are occupied by major large blood vessels, indicating that the transformation of large vessel areas from the template to subject space was successful. Fig. 3B shows the first principal component of the cardiac-induced signal acquired from the ROI. The phases calculated by our approach were slightly shifted, but in good agreement with the phases calculated by data from the pulse oximeter (Fig. 3C). Note that the phase shift shown in Fig. 3C does not change the sum of square of the coefficients, although each coefficient value is different (Eq. (1)). The correlation coefficients between both phase values were r > 0.8. We compared the cerebral pulsatility maps across the whole brain acquired with our proposed technique and pulse-oximeter methods. The image intraclass coefficient between the pulsatility maps obtained by the pulse oximeter and our approach is 0.85, with a 95% confidence interval of (0.84, 0.86). A small confidence interval indicates high certainty in the estimation. This means that the method calculates the phases accurately from rs-fMRI data without the external pulse oximeter data. Fig. 3D shows the cerebral pulsatility map obtained by our method from one participant with MCI. All subjects showed similar results in that the areas of high cardiac-induced signal were localized to the regions of the major arterial blood vessels, the venous sinuses, and the ventricles. For the reliability of pulsatility measured on two separated days (10.3 ± 13.2 days apart), the image intraclass coefficient was 0.68, with a 95% confidence interval of (0.66, 0.70).

Fig. 3.

(a) ROI for large vessels transformed from the MNI template to individual’s space. (b) The normalized first component of cardiac pulsatile signal obtained from the ROI of a. (c) The comparison of cardiac phase calculated from our method vs. pulse-oximeter from 4 separated runs. Different colors indicate four different runs. (d) Pulsatility map calculated from one MCI subject. Pixels with larger than 3σ (99.73%) of the null distribution are displayed. High signal intensity indicates high pulsatility.

Relationship of the pulsatility and brain volume

The averaged optimal regularization parameters that control the sparsity of the SPLS weights were 10 ± 3.1% for the pulsatility variables and 8.3 ± 5.5% for the brain volumetric variables across the splits. The averaged grids of correlation values from 9 splits are shown in Supplementary Fig. S1.

Table 2 shows the absolute correlation values between the latent variables of the pulsatility and brain volume that were calculated by projections of weighting vectors onto the holdout dataset. P-values for 9 splits are also shown. The omnibus null hypothesis (Homni) was rejected for both the first and second associations. Thus, the associations between subsets of pulsatility and brain volumetric measures are statistically significant. In addition, the permutation tests with the averaged weighting vectors on the whole data (across 9 splits) were also found to be significant for both the first and second pulsatility-volume relationships (p < 0.001).

Table 2

SPLS p-values computed with 10,000 permutations (statistically significant results are shown in bold, p < 0.0056). Holdout correlation was calculated by |Corr(X)|. The second relationship was performed by projection deflation.

split		1	2	3	4	5	6	7	8	9	Reject Homni
holdout corr.		0.74	0.49	0.81	0.75	0.72	0.90	0.58	0.70	0.89
p-value	1st	0.0048	0.1449	0.0054	0.0050	0.0252	0.0001	0.0978	0.0321	0.0006	Yes
	2nd	0.1600	0.1061	0.0069	0.0338	0.0012	0.0068	0.1666	0.0106	0.0001	Yes

The brain volumetric weights of 9 splits identified by SPLS were averaged. For the first pulsatility-volume relationship, the pulsatility was positively associated with the volumes of lateral ventricle, followed by choroid plexus, inferior lateral ventricle, and 3rd ventricle and negatively associated with hippocampus, ventral DC, and thalamus (Fig. 4a). For the second pulsatility-volume relationship, the pulsatility has an additional negative relationship with the amygdala and brain stem volumes and small positive relationship with age (Fig. 4b). The statistically significant voxels (p < 0.01) of averaged pulsatility weights for the first and second pulsatility-volume relationship are displayed in Fig. 4c and d. Overall, positive relationship of pulsatility in the arteries and negative relationship of pulsatility in venous sinus and CSF were observed; the anterior cerebral artery, the middle cerebral artery, insular artery and posterior cerebral artery were positively associated and ventricle and sagittal sinus were negatively associated with pulsatility. The voxel-wise pulsatility weights of the second pulsatility-volume relationship were similar, but noisier and less localized than the first relationship because the projection deflation removed the effect explained by the first weighting vectors from the dataset (Fig. 4d).

Fig. 4.

(a) The averaged volume weight vectors for the first pulsatility-volume relationship are displayed. The lateral, interior lateral and 3rd ventricles and choroid plexus are found to be positively associated with the pulsatility, while hippocampus, ventral DC, and thalamus show negative relationship with the pulsatility (b) Additionally, the amygdala and brain stem volumes are negatively associated for the second pulsatility-volume relationship. Age shows small positive relationship. Blue labels: 1. Hippocampus, 2. Lateral ventricle, 3. Inferior lateral ventricle, 4. Ventral DC, 5. 3rd ventricle, 7. Brain stem 8. Amygdala, 9. Thalamus, 13. Choroid plexus, 46. age. Labels for other volume variables are listed in supplementary materials (Table S1). (c) Spatial distribution of the averaged cerebral pulsatility correlated with brain volume for the first pulsatility-volume relationship. The statistically significant voxels were displayed (p < 0.01). The pulsatility of arteries shows positive relationship, while that of ventricles and venous sinus shows negative correlation with the brain volumes. (d) Spatial distribution of the pulsatility for the second pulsatility-volume relationship. Color bar: the scale of weighting vectors. Error bars: SEM.

Mediation pathway between the pulsatility and cognitive measures

The results from the mediation analysis, the statistical significance (p-values) of each pathway, are summarized for the first component in Table 3. Since the volumes identified by SPLS are closely correlated with each other, instead of using each volume separately, we analyzed the indirect pathway using the first component of the brain volumes by PCA and the latent variable of volume by SPLS as a mediator. The first PCA-driven brain volume component has 97.4% of the shared variance between the identified brain volumes. It indicates very high correlations between the volumes. The indirect pathway between the pulsatility and MOCA, using the brain volumetric component as a mediator, was statistically significant for both PCA and SPLS approaches (p < 0.001). The indirect pathways between the pulsatility and each of NP batteries were also statistically significant (p < 0.001). The direct pathway, which is the relationship between pulsatility and each of cognitive measures after controlling for the brain volume mediator, was also statistically significant for both approaches (p < 0.01), except in visuospatial function. However, the indirect pathway for 1.53% of the shared variance of the second PCA-driven brain volume component was not statistically significant for all mediation analysis, while the second volume component of SPLS showed statistically significance for all cognitive measures (Supplementary materials, Table S3). Similarly, the direct pathway of PCA showed no significance, while SPLS showed statistical significance for all cognitive measures (p < 0.01), except for visuospatial function.

Table 3

Results of mediation analysis are summarized for the first component (p-value for each pathway).

	MOCA	Memory	Language	Attention	Executive	Visuospatial
Indirect pathway: relationship between pulsatility and each cognitive measure via segmented brain volume
PCA	5.2 × 10−5	7.5 × 10−5	7.7 × 10−5	6.8 × 10−5	6.4 × 10−5	0.0014
SPLS	1.5 × 10−4	1.3 × 10−3	1.1 × 10−4	1.2 × 10−4	1.3 × 10−4	8.3 × 10−5
Direct pathway: relationship between pulsatility and each cognitive measure while controlling for segmented brain volume
PCA	3.7 × 10−4	0.004	0.0012	0.0091	0.0029	0.3025
SPLS	1.1 × 10−4	1.1 × 10−4	9.5 × 10−5	1.2 × 10−4	1.3 × 10−4	0.3827

Discussion

We have successfully obtained a map of cerebral pulsatility directly from rs-fMRI data with good reliability. We applied the holdout framework (for unseen data) to avoid overfitting the data and repeated the process for nine randomly divided splits to test each holdout data by SPLS, which is more robust than the approach with a single split. This holdout framework approach increases generalizability and reproducibility. We found statistically significant association between the detected cerebral pulsatility regions and volumes of hippocampus, amygdala, thalamus, Ventral DC, brain stem, choroid plexus, and ventricles. In addition, these identified brain volumes were found to be significant mediators in the pathway of the pulsatility to cognitive measures.

Technical aspect of pulsatility map

Our main purpose of using fMRI data is to map cerebral pulsatility, unlike the original purpose of RETROICOR (Glover et al., 2000), which has been for removing the physiological noise from fMRI. Similar approaches of pulsatility measurement have been attempted using the signal intensity of rs-fMRI (Theyers et al., 2019; Tong et al., 2014). Theyers et al. sorted the signal intensities of image volumes with their corresponding positions in the cardiac cycle. The resampled time series were then fitted with a Fourier series to obtain the pulsatility map. Tong et al. used multiple time-shifted signals generated from a pulse-oximeter as cardiac signal regressors to acquire dynamic cardiac-related pulsation from BOLD data with a very short TR (=0.4 s) using SMS imaging. This sampling rate is sufficient to meet Nyquist sampling requirements for heart rate. However, this approach is limited by a relatively low spatial resolution (3.4 × 3.4 × 3.5 mm3) when used to cover the whole brain. At longer TR in fMRI experiments with higher spatial resolutions (e.g. 2 mm isotropic, as used this study), the cardiac pulsations can be undersampled, and aliasing of cardiac signal is difficult to avoid. With a heart rate in the range of 1–1.6 Hz (60–100 BPM), a TR of 0.312 s (=3.2 Hz) or shorter is required to meet Nyquist sampling requirements. Most of all, the RETROICOR method has demonstrated excellent identification of cardiac-induced signals with a relatively long TR (Dagli et al., 1999; Glover et al., 2000; Lund et al., 2006). However, one common prerequisite of all these methods is an external measurement of cardiac cycle, generally acquired by the scanner’s built-in pulse oximeter. Thus, if an external cardiac pulsation measurement is unavailable or inaccurate, then the pulsatility cannot be measured. In this work, the rs-fMRI derived pulsatility map shows excellent agreement with the pulsatilty map obtained by an external measurement of cardiac cycle. Some participants of our study, especially those that are elderly or impaired had trouble with measurements of pulse oximeter signal. The optical plethysmograms can be easily distorted and/or the device can fall off from the subject’s finger with small hand movements. The technique described here overcomes this issue by obtaining cardiac phase time-series directly from the rs-fMRI data itself. Note that our method can also work for long TR, as the RETROICOR method identified the aliased components (Glover et al., 2000; Restom et al., 2006).

A simple method using the amplitude of fluctuations in the rs-fMRI BOLD signal (a.k.a. RSFA) has also been suggested to evaluate vascular function and white matter pulsatility on aging (Makedonov et al., 2013; Tsvetanov et al., 2015). However, this is valid only when a fraction of cardiac-induced signals to the total fluctuation is consistent with aging. The BOLD fluctuation is affected by not only cardiac pulsation and its variability, but also other sources as well, such as neurovascular coupling, spontaneous oscillation of vessels, respiratory-induced changes, and head motion. Since the highest fractional variance of physiological noise to total BOLD signal variance was ~0.7 (Behzadi et al., 2007), except in large vessels and ventricle regions, the fraction of cardiac pulsation to the total fluctuation is expected to be low. Thus, it is better to extract the pulsatile signals from the complex mixture of the BOLD signal to avoid other component contaminations.

We did not compare the accuracy of our pulsatility map to other pulsatility methods that more directly measured flow velocity at a restricted location, such as phase contrast MRI or transcranial Doppler. A recent study demonstrated the pulse wave velocity (PWV) measured by the phase-contrast MRI was related to the variance of rs-fMRI (Hussein et al., 2020). Thus, we expect that our pulsatility is associated with PWV because our pulsatility is also derived from the rs-fMRI signal fluctuations. In addition, our findings agree well with multiple sources of pulsatility measurements in the previous study; arterial pulsatility index, CSF volume pulsatility, arterial flow volume pulsatility and CSF pulse pressure showed positive relationship with ventricular volume and negative relationship with hippocampus volume (Wahlin et al., 2014). Thus, our pulsatility map can well characterize the cardiac-induced pulsatile signal in the brain.

Relationship between the pulsatility and volumes

The PLS method is a computationally efficient tool to discover features of multivariate covariance between two datasets when the number of variables (many voxels of brain pulsatility and multiple segmented brain volumes in this study) are much larger than the number of samples (Krishnan et al., 2011). Although the latent sources of variability are not directly observable in the original data, the strength of this approach is the reduction of the high dimensional variable sets of the original data to the compact linear representations that best characterize relationship between two sets of variables. However, one drawback of PLS is that it is difficult to interpret the high-dimensionality of variables in the model because all features in both views (X and Y) are included. Thus, it is difficult to separate irrelevant features (noise) from the model. We applied SPLS to select subsets of pulsatility and brain volume based on the lasso constraint, which improves the interpretability. The sparsity of the variables is determined by the regularization parameters (c and c), i.e., the lower values of c and c determine the stronger l1-norm constraint. It induces less numbers of pulsatility voxels and volumes for the association. A pair of regularization parameters was automatically determined by the grid search. In the pulsatility-brain volume relationship, 3–8 brain volumes were selected by SPLS depending on the sparsity of each split. For example, higher values of c and c in a split resulted in a greater number of brain volumes associated with the pulsatility, compared to other splits. However, most of identified brain regions obtained by SPLS that were associated with pulsatility in both the first and second relationships overlapped across splits.

The volumes of choroid plexus and ventricles were mostly associated with the pulsatility. The choroid plexus, located in brain ventricles, is responsible for production of cerebrospinal fluid (CSF) and its volume is related with pathogenesis of neurodegenerative diseases, such as the clearance of CSF proteins (Tadayon et al., 2020). Larger choroid plexus volume is also related to smaller total gray matter and larger lateral ventricle volumes (Lizano et al., 2019). Recently, there is a report that an increased volume of choroid plexus is associated with increased vascular burden (Egorova et al., 2019), which could be related to an increase in cerebrovascular pulsatility. Additionally, the pulsatility was observed to be associated with hippocampus, amygdala and thalamus volumes. These regions are responsible for emotion and learning and formation of memory, and these volume changes have been related to cognitive and psychiatric dysfunctions as well as a neuroimaging evidence of characterizing AD progression (de Jong et al., 2008; Miller et al., 2013; Poulin et al., 2011; Wachinger et al., 2016).

The brain is vulnerable to increased arterial pulsatility (Webb et al., 2012; Wohlfahrt et al., 2014). With aging, increased arterial stiffness, that is inversely related to the vessel compliance, can be associated with increased arterial pulsatility (Wagshul et al., 2011). The compliance of the central arteries usually dampens hemodynamic pulsations in an effective way to protect the cerebral microcirculation from high pulse flow. However, the arterial stiffness can induce less effective dampening and the elevated pulsatility induces increasing stress on the capillary microvasculature system that may lead to brain tissue damage. These changes in capillary hemodynamics may be related to reduction in regional volumes that are observed in this study.

Age, one of the nuisance variables, was also found to be associated with the pulsatility in the second pulsatility-brain volume relationship, but the amplitude of relationship was much smaller than the identified brain volumes. It can be interpreted that the effect of pulsatility on brain volume goes beyond the effect on age.

The spatial distribution of pulsatility was found in adjacent feeding arteries of detected volumes by SPLS. The hippocampus region is supplied from many arteries, most of which branch off from either the posterior cerebral artery or the anterior choroidal artery (Erdem et al., 1993), and the thalamus derives its blood supply from all branches of the posterior cerebral artery. Supply to the choroid plexus of the lateral and third ventricles is provided by the anterior choroidal artery and by several choroidal branches of the posterior cerebral artery. On the other hand, decreased pulsatility of venous sinus and CSF areas were also observed. The pulsation of the venous sinus tree and ventricular areas is associated with continuous removal of metabolic waste. A reduction of CSF circulation (formation and resorption) is related with cognitive decline (Benveniste, 2018; Iliff et al., 2013). Reduced pulsatility in CSF also suggests that increased ventricular volume may be related with decrease in the gray matter volumes, without an elevated pressure in the brain (Wilson, 2016). Otherwise, the pulsatility in CSF would be increased by the enlargement of the ventricles that is seen in normal pressure hydrocephalus. This effect is different from that observed in AD (Bateman, 2003).

Although the relationship between the pulsatility and brain volume was robustly evaluated by the multiple holdout method, one limitation of this approach is the long computing time. The main bottleneck of this computational burden is the permutation testing. It took ~11 h for 10,000 permutations on Intel Core i7 CPU at 2.6 GHz. However, multiple runs can be simultaneously executable for multiple splits in multi-core process to reduce the total compute time.

Mediation analysis

The mediation analysis showed the identified brain volumes are statistically significant mediators between the pulsatility and each cognitive measure. This result agrees well with the previous report that excessive arterial pulsatility has been linked to poorer cognitive performance and brain atrophy among the elderly (Mitchell et al., 2011). We hypothesized a pathway from the pulsatility through brain volume to cognition. As a cross-sectional examination, this study cannot prove causality, and the direction of pathway can be alternative. The indirect pathway of pulsatility between brain volume and each cognitive measure was also found to be statistically significant (p < 0.001), except for visuospatial function. Therefore, longitudinal approaches are needed to determine the causality. Although the causality order of pulsatility and volume could be alternated, the brain volumetric measure in MRI reveals a relatively early indicator than the severity of cognitive outcomes in AD pathological cascade model (Balsis et al., 2018; Jack et al., 2013). Thus, the cerebral pulsatility could also be a sensitive cerebrovascular marker in the relatively early stage of the disease, and therefore helpful to better understanding of early pathophysiological mechanisms of AD progression.

Conclusion

The cardiac-induced pulsatility in the brain regions of its supply arteries, venous sinus, and CSF was positively associated with the ventricles and choroid plexus, and negatively related with most of the deep subcortical volumes. The indirect mediating pathway of these volumes are statistically significant between the pulsatility and each of the cognitive measures. Thus, the region-specific changes of pulsatility over the entire brain can provide unique information for cerebrovascular alteration and should be considered as a potentially important factor for cognitive decline in future studies.