Characterizing the spatial patterns and determinants of cerebrospinal fluid pseudorandom flow in the human brain with low b-value diffusion MRI Open Access

The cerebrospinal fluid (CSF) circulates through a complex network of inter-connected basilar cisterns, ventricular system, and extra-axial CSF spaces, which are semi-compartmentalized by brain parenchyma and meningeal membranes (Lu et al., 2022; Møllgård et al., 2023; Rasmussen et al., 2022). Cardiovascular pulsatility drives CSF flow and propels it along the perivascular spaces (Lloyd et al., 2020), where CSF is in continuous exchange with the interstitial fluid and plays a vital role in maintaining brain homeostasis and clearance system (Klostranec et al., 2021; Rasmussen et al., 2022). Disruptions in CSF circulation can lead to hydrocephalus and intracranial pressure disorders (Bothwell et al., 2019; Reeves et al., 2020). Moreover, abnormalities in CSF flow dynamics have been observed in the context of aging, neurodegenerative disorders, and small vessel disease (Rasmussen et al., 2022).

Prior work has used invasive and non-invasive techniques to study CSF flow in the human brain. However, studies using imaging approaches that involve intrathecal injection of contrast agents have often been limited by the number of participants or inclusion of healthy individuals due to the invasive nature of the procedure (Eide et al., 2021). Additionally, while CSF flow imaging techniques, such as phase-contrast imaging (Edelman et al., 1986; Enzmann & Pelc, 1991) and spin-labeling MRI (Yamada et al., 2008, 2015), are non-invasive, they cannot provide whole-brain maps of CSF flow in relatively short acquisition times. Moreover, phase-contrast imaging techniques that effectively capture CSF flow in regions with high velocities, such as the basilar cistern, their sensitivity to slower CSF flow—particularly within the lateral ventricles and extra-axial spaces—are limited. This reduced sensitivity arises from the use of gradient echo imaging and the need for prolonged echo times to measure slow-velocity CSF flow, which consequently results in signal loss due to T₂* decay (Fig. 1A). These limitations have resulted in studies that have primarily focused on flow patterns across specific CSF passageways (e.g., cerebral aqueduct, foramina of Monro, and craniocervical junction) (Yamada et al., 2015).

In contrast to phase-contrast imaging, diffusion-weighted MRI at low b-values (low-b dMRI) is sensitive to complex and slow CSF motion (Bito et al., 2021; Harrison et al., 2018; Taoka et al., 2019, 2021; Wen et al., 2022). Low-b dMRI employs a pulsed gradient spin-echo sequence, allowing it to capitalize on the long T₂ relaxation times of the CSF and enhance its sensitivity to slow-velocity CSF motions. Importantly, low-b dMRI is able to provide whole-brain maps of CSF motion utilizing rapid echoplanar imaging (EPI) readouts (Fig. 1A). Rather than changes in signal phase, low-b dMRI relies on signal magnitude that decreases with pseudorandom fluid flow (more than expected by simple diffusion; Fig. 1). At low b-values, CSF mean pseudo-diffusivity (MΨ, analogous to mean diffusivity at high b-values) reflects the variance in flow velocities, capturing both laminar and/or mixing flow, thus serving as a biomarker of CSF effective motility within the voxel (Bito et al., 2021, 2023).

Understanding CSF flow dynamics is further hindered by its complex nature. CSF flow is driven in tandem by oscillatory changes in brain volume across the cardiac and respiratory cycles and the pumping effects of cerebral vasculature (Kedarasetti et al., 2020; Mestre et al., 2018; Wagshul et al., 2011). Importantly, these effects are moderated by intracranial anatomy and compliance of the CSF spaces (Preuss et al., 2013; Rasmussen et al., 2022). Aside from cerebral aqueduct CSF flow studies using phase-contrast MRI (Rohilla et al., 2024; Sartoretti et al., 2019), most prior studies in healthy populations involved small cohorts, leaving demographic and intracranial anatomical effects on CSF flow in other brain regions largely unexamined. Therefore, comprehensive multimodal imaging studies are required to elucidate the effects of demographics, intracranial anatomy, and physiologic drivers on CSF flow patterns.

In this study, we sought to gain insight into the spatial organization of intracranial CSF circulation and characterize its underlying structural and physiological underpinnings. We hypothesized that CSF flow dynamics in the cranial cavity exhibit region-specific covariance patterns that are shaped by similar physiological drivers and structural constraints, enabling the derivation of a data-driven atlas for parts-based analysis of CSF flow. To this end, we capitalized on recent advances in low-b dMRI to evaluate CSF flow in the human brain (Bito et al., 2021; Harrison et al., 2018; Taoka et al., 2019, 2021; Wen et al., 2022). Accordingly, we derived maps of CSF MΨ in a large multi-b-value cohort, part of the Open Access Series of Imaging Studies (OASIS). We examined CSF MΨ throughout the brain using a whole-brain voxel-based analysis framework we developed, termed CSF pseudo-diffusion spatial statistics (CΨSS). Employing a hypothesis-driven approach, we applied seed-based correlation analysis to determine if the covariance patterns of CSF MΨ are consistent with the expected intracranial CSF flow patterns. Next, we utilized non-negative matrix factorization (NMF) to examine the covariance structure of CSF MΨ, enabling an unbiased, data-driven identification of the spatial patterns of CSF circulation. In multiple independent low-b dMRI datasets, we showed that these CSF circulation spatial patterns represented distinct pseudo-diffusion magnitude and flow directionality. Finally, by employing multiple neuroimaging modalities, we identified the major determinants of CSF MΨ and showed that individuals with aberrant CSF flow patterns tended to have structural incidental findings on their neuroradiological examinations.

OASIS-3 (Open Access Series of Imaging Studies: https://sites.wustl.edu/oasisbrains/) is a retrospective compilation of clinical and imaging data from >1,300 participants that were collected across several studies through the Charles F. and Joanne Knight Alzheimer Disease Research Center (Knight ADRC: https://knightadrc.wustl.edu/) at Washington University in St. Louis (LaMontagne et al., 2019). Participants include cognitively unimpaired adults and individuals at various stages of cognitive decline. This study utilized low-b dMRI data from two imaging subsets within OASIS-3 that employed distinct imaging protocols and MRI scanners (Table 1): (i) multi-low b-value subset: n = 187 (after exclusion of 10 participants due to motion artifacts), 103 female, average age: 71.5 years [range: 46.2–92.2 years], 163 cognitively normal (Clinical Dementia Rating [CDR]: 0) and 18 individuals with very mild cognitive impairment (CDR: 0.5); (ii) b100 subset: n=103 (after exclusion of a single participant due to motion artifact), 54 female, average age: 71.4 years [range: 46–91.6 years], 92 cognitively normal (CDR: 0) and 8 individuals with very mild cognitive impairment (CDR: 0.5). There were 26 overlapping participants between the two OASIS-3 subsets.

The Microstructure-Informed Connectomics (MICA) dataset (Royer et al., 2022) provides raw neuroimaging data collected from 50 healthy participants (Table 1). The imaging dataset is available on the Canadian Open Neuroscience Platform’s data portal (https://portal.conp.ca).

The Comprehensive Diffusion MRI dataset (CDMD) included imaging data from 26 healthy participants (Table 1). The imaging dataset is available as a figshare collection (https://doi.org/10.6084/m9.figshare.c.5315474.v1).

All study procedures were approved by the Institutional Review Board at our institution (IRB# 202404163), and written informed consent was obtained from all participants. High-resolution structural imaging data were collected from 11 healthy young adults (age: 25.5 ± 3.6 [range: 22–36] years; 7 males).

The imaging cohorts in the OASIS-3 dataset were defined based on their dMRI imaging schemes. For the subset of OASIS-3 participants in the multi-low b-value group, brain MRI scans were conducted using a 3 T BioGraph mMR scanner (Siemens, Erlangen, Germany). Axial diffusion-weighed pulsed-gradient spin-echo echoplanar images (EPI) were acquired using a multi-low b-value protocol with the following imaging parameters: TE, 86 ms; TR, 10,300 ms; voxel size, 2 × 2 × 2 mm³; field-of-view, 224 × 224 mm²; slice number, 80. Each diffusion gradient had a unique b-value ranging from 0 to 1400 s/mm², with a total of 26 volumes. For the subset of OASIS-3 participants in the b100 group, whole-brain diffusion-weighted EPI images were obtained using a 3 T Magnetom Vida scanner (Siemens, Erlangen, Germany) equipped with a 64-channel head coil with the following imaging parameters: TE, 79 ms; TR, 5,800 ms; voxel size, 2 × 2 × 2 mm³; field-of-view, 220 × 220 mm²; slice number, 80; multi-band factor, 2. Diffusion-weighted images were acquired in the axial plane along 66 gradient directions, with b-values of 100 s/mm² (16 directions), 250 s/mm² (10 directions), 500 s/mm² (12 directions), 1,000 s/mm² (12 directions), 1,500 s/mm² (10 directions), 2,000 s/mm² (6 directions), and 3 b0 volumes. Additionally, pulsed arterial spin labeling (PASL) imaging data were available for a subset of participants (n = 164) (LaMontagne et al., 2019). Axial 2D PASL images were acquired with proximal inversion with control of off-resonance effects (PICORE) tagging method and the following imaging parameters: TR, 3400 ms; TE, 13 ms; TI₁,700; TI₂, 1900 ms; in-plane resolution, 4 × 4 mm²; slice thickness, 5 mm; 9 slices; and 52 label/control pairs with an M0 reference image.

In the MICA study (Royer et al., 2022), brain MRI images were acquired with a 3 T Magnetom Prisma-Fit scanner (Siemens, Erlangen, Germany) equipped with a 64-channel head coil. Multi-shell diffusion MRI data were acquired using a diffusion-weighed pulsed-gradient spin-echo EPI sequence with the following parameters: 10 diffusion weighting directions at b = 300 s/mm² (only b = 300 s/mm² and b0 images were used for this study), TE, 64.40 ms; TR, 3500 ms; voxel size, 1.6 × 1.6 × 1.6 mm³, field-of-view, 224 × 224 mm²; and multi-band factor, 3. To correct for susceptibility induced image-distortions of, b0 images were also acquired in reverse-phase encoding direction.

In the CDMD study, all data were acquired on the 3 T Magnetom Connectome MRI scanner (Siemens, Erlangen, Germany) equipped with a maximum gradient strength of 300 mT/m and a custom-built 64-channel phased array head coil (Tian et al., 2022). Multi-shell diffusion MRI data were acquired using a sagittal diffusion-weighed pulsed-gradient spin-echo EPI sequence with the following parameters: TE, 77 ms; TR, 3800 ms; voxel size = 2 × 2 × 2 mm³, field-of-view, 216 × 216 mm, multi-band factor, 2. Diffusion-weighted images acquired at three different b-values were used in this study (b = 50, 350, 800 s/mm²; acquired along 32 diffusion encoding directions uniformly distributed on a sphere). Five b0 image volumes with reversed-phase encoding direction were acquired to correct for susceptibility-induced image distortions.

Brain MRI images were acquired with a 3 T Magnetom Prisma-Fit scanner (Siemens, Erlangen, Germany) equipped with a 64-channel head coil. Mid-sagittal, high-resolution balanced steady state free precession (b-SSFP) images were obtained using a 3D-CISS (constructive interference in steady state) sequence with the following parameters: TE, 1.98 ms; TR, 4.51 ms; flip angle: 33; voxel size, 0.5 × 0.5 × 0.5 mm3, field-of-view, 237 × 200 mm². 3D-CISS is a fully balanced and inherently flow-compensated gradient-echo sequence, providing fine anatomic details within the CSF spaces (Dinçer et al., 2009; Lu et al., 2022). This method achieves high signal-to-noise ratio in CSF due to its balanced steady-state acquisition and the high T₂/T₁ ratio characteristic of CSF.

At low b-values, the observed pseudo-diffusion tensor Ψ can be approximated as the sum of two components (Fig. 1): the covariance matrix of the 3D velocity probability density function, which represents pseudorandom flow (V), and the diffusion tensor (D, Fig. 1) (Bito et al., 2021, 2023). This relationship is expressed as:

Here, $τd$ represents the diffusion time. Hence, when the physical properties of CSF, such as temperature and viscosity, remain unchanged, the MΨ derived from the Ψ tensor is indicative of the flow velocity variance within a voxel caused by laminar and/or mixing movements of CSF.

After corrections for Gibbs ring artifact (mrdegibbs tool, part of MRtrix3) (Kellner et al., 2016; Tournier et al., 2019) and eddy current-induced distortions (eddy_correct tool with spline interpolation, part of FSL) (Jenkinson et al., 2012), skull stripping was performed (Brain Extraction Tool [BET], part of FSL) (Smith, 2002). The diffusion tensor model was fit to the low-b-value diffusion-weighted images (multi-low b-value subset: b: 0-350 s/mm² [7 non-b0 volumes]; b100 subset: b = 100 s/mm², 16 directions, 3 b0 volumes) to generate low b-value MΨ maps as a measure of pseudorandom flow magnitude. The diffusion tensor model was also fit into the intermediate/high b-value images in both datasets (multi-low b-value subset, b: 500–1,400 s/mm²; b100 subset, b: 500–2,000 s/mm²). The resulting FA and mean diffusivity images were fed into a 2-tissue segmentation algorithm (Atropos, part of ANTs; http://stnava.github.io/ANTs/) to generate white matter and CSF partial volume maps (Avants, Tustison, Wu, et al., 2011). Similar to gray matter-based spatial statistics (Nazeri et al., 2015, 2017), partial volume maps were used to generate pseudo-T₁-weighted images. MΨ and pseudo-T₁ images were used to create study-specific templates using the antsMultivariateTemplateConstruction2.sh workflow (Fig. 1). MΨ and CSF partial volume maps were subsequently registered to the template space. The final mask was generated by keeping voxels with a CSF fraction greater than 0.7 in more than 65% of the subjects. Voxels with CSF fraction < 0.7 were then filled with the average of the surrounding satisfactory voxels. Finally, spatial smoothing was applied within the CSF mask using the 3dBlurInMask function from the AFNI software package (Cox, 2012), with a full width at half maximum value of 3 mm. Building on tract-based spatial statistics (Smith et al., 2006) and gray matter-based spatial statistics (Nazeri et al., 2015) for gray and white matter analysis, CΨSS uses only diffusion-weighted MRI data for image segmentation and registration. By leveraging the intrinsic contrast between brain tissues in maps derived from diffusion-weighted MRI—specifically the distinct diffusion properties of CSF and parenchyma—this method enables accurate tissue segmentation in the native diffusion space, eliminating the need for registration to structural images and minimizing errors caused by EPI distortions.

Seed-based correlation analysis was performed separately in both OASIS cohorts. Regions of interest (ROIs) were defined in the template space by a neuroradiologist (A.N.) using ITK SNAP v3.8 (Yushkevich et al., 2016). The selected ROIs included the premedullary cistern, cerebral aqueduct-4^th ventricle interface, foramina of Monro, left Sylvian fissure, and right Sylvian fissure. Mean MΨ values were extracted from these ROIs for each subject using the fslmeants function in FSL. Voxel-wise Pearson’s correlation maps were generated for each ROI using the 3dTcorr1D function in AFNI. Voxelwise correlation values with false-discovery rate (FDR)-adjusted p-value < 0.01 were considered significant (corrected for five ROIs).

We utilized NMF to identify the CSF flow patterns in which the CSF MΨ covaried consistently among participants. NMF is an unsupervised machine-learning method, which is used for dimensionality reduction and matrix decomposition (Lee & Seung, 1999). In neuroimaging, NMF has been utilized to uncover interpretable covariance structures and to detect patterns within high-dimensional imaging data (Earnest et al., 2024; Nazeri et al., 2022; Sotiras et al., 2015, 2017). NMF decomposes a given non-negative matrix $X$ into two non-negative matrices $W$ (basis matrix) and $H$ (coefficient matrix) such that $X≈WH$⁠. For the purpose of this study, the two-dimensional non-negative data matrix $X$ was constructed by joining vectorized CSF MΨ maps in the template space (⁠$X=[x1,…,xN],xi∈ℝ+D$⁠, where $D$ is the number of CSF voxels in the template space and $N$ is the number of subjects). The resulting basis matrix $W(W=[w1,…,wK],wi∈ℝ+D)$ is composed of $K$ columns, where $K$ is the user-specified number of patterns. These columns represent the estimated CSF flow patterns. Each row in the coefficient matrix $W$ represents a voxel, and the weights in the row indicate the relative contributions of that voxel to the CSF flow patterns. The matrix $H (HT=[h1,…,hK],hi∈ℝ+N)$ contains subject-specific coefficients for each CSF flow pattern, which indicate the contribution of each pattern in reconstructing the original CSF MΨ map. The implementation of NMF used in this study (https://github.com/asotiras/brainparts) enforces orthonormality constraints for the estimated covariation patterns (⁠$WTW=I$⁠, where $I$ is the identity matrix) and projective constraints for their respective participant-specific coefficients $(H=WTX)$ (Sotiras et al., 2015; Yang & Oja, 2010). The implementation of NMF has been discussed elsewhere in more detail (Sotiras et al., 2015, 2017).

Consistent with prior studies using this technique (Nazeri et al., 2022; Sotiras et al., 2015, 2017), we ran multiple NMF solutions requesting 2–25 patterns to obtain a range of possible solutions for comparison. The optimal number of patterns (⁠$K$⁠) was selected based on the reconstruction error and stability of the solutions. The reconstruction error was calculated for each solution as the Frobenius norm between the data matrix and the NMF approximation. To determine the stability of our results (Nazeri et al., 2022; Sotiras et al., 2015, 2017), we performed a split-half reproducibility analysis with bootstrapping (with replacement) to examine the stability of NMF solutions across the range of possible resolutions (Ben-Hur et al., 2002; Lange et al., 2004; Thomas Yeo et al., 2011). To ensure comparable age and sex distributions, the split-halves were generated using the anticlust package in R (Papenberg & Klau, 2021). We utilized two criteria to assess the stability of the NMF solutions: reproducibility, defined as the solution with the highest mean adjusted Rand index (ARI) across 20 split-half bootstraps, and reliability, defined as the solution with the lowest standard deviation of ARI across the bootstraps. ARI is a measure of set similarity that accounts for chance, with a value of 1 indicating a perfect match and a value close to 0 indicating a random partitioning (Hubert & Arabie, 1985). An ARI value greater than 0.75 is considered excellent (Fleiss et al., 2013).

High-resolution 3D-CISS images from the BraFTI Healthy Neurofluid Imaging study were utilized to assess the fine anatomical structure of CSF spaces. Slice aliasing was removed using fslroi. A template was generated from the resulting images using antsMultivariateTemplateConstruction2.sh.

Low-b dMRI data from the MICA dataset (b = 300 s/mm²) were corrected for Gibbs ring artifact, susceptibility-induced distortion using the FSL’s topup function, and reverse-phase encoding b0 images. Data were additionally corrected for eddy-current-induced distortion and motion using FSL’s eddy function. Following skull-stripping with FSL’s Brain Extraction Tool (BET), the diffusion tensor model was applied to the preprocessed data to create MΨ maps.

For the CDMD dataset (Tian et al., 2022), the available preprocessed dMRI data were used for further analysis. MΨ and FA maps were generated by fitting the diffusion tensor model to b = 50 s/mm² dMRI data. Mean square displacement (MSD) maps were created by fitting the Laplacian-regularized mean apparent propagator (MAPL) MRI model (Fick et al., 2016) to the multi-shell dMRI data (b-values: 50, 350, 800 s/mm²) using dipy v1.5.0 (Garyfallidis et al., 2014). Given that the image acquisition was not gated, residual error in the diffusion tensor fit may reflect a combination of temporal variations from cardiac and respiratory cycles, in addition to noise. To quantify this, the error between the predicted and measured dMRI signal at b = 50 s/mm² was quantified by calculating the voxel-wise root mean square error (RMSE). This was normalized by the mean signal intensity of the corresponding b-value shell to generate the normalized RMSE (nRMSE), providing a relative measure of error across the image.

In both MICA and CDMD datasets, the CWW masks from the OASIS-3 template space were nonlinearly registered to the individual MΨ images using the antsRegistrationSyN.sh script (Avants, Tustison, Song, et al., 2011). Mean diffusion metric values from each CWW ROI were extracted in the native diffusion space.

T1-weighted images from the OASIS-3 dataset were segmented using FreeSurfer v5.3. (https://surfer.nmr.mgh.harvard.edu/) (Fischl, 2012). The image segmentations were reviewed by a trained lab member of the OASIS project (LaMontagne et al., 2019). TkMedit (part of FreeSurfer; http://freesurfer.net/fswiki/TkMedit) was used to revise images that failed the quality control. The revised images were rerun through the FreeSurfer pipeline. Total intracranial volume, brain parenchymal fraction (brain parenchymal volume divided by total intracranial volume), and ventricular volumes (normalized by total intracranial volume) were extracted for further analyses.

Cerebral blood flow (CBF) quantification was performed using the oxford_asl package v4.0.27, part of the BASIL toolbox within FSL (Chappell et al., 2009; Jenkinson et al., 2012). This procedure involved motion correction, automated spatial regularization, label-control subtraction, relative CBF quantification, and conversion of relative CBF to absolute physiological units (ml/100 g/min) using the M0 image (Chappell et al., 2009). Model parameters were selected based on the ASL Consensus Paper for a magnetic field strength of 3 T (Alsop et al., 2015): arterial blood longitudinal relaxation time (T_1b), 1650 ms; brain tissue longitudinal relaxation time (T_1t), 1300 ms; inversion efficiency (α), 0.98. The TI₂ value was adjusted slice by slice with an inter-slice acquisition time difference of 46 ms. A gray matter mask (gray matter partial volume effect > 0.5; derived from the fsl_anat script) was created in the ASL native space and used to extract mean gray matter CBF values.

In the OASIS study, aberrant CSF flow patterns were defined as either (i) one or more CWW with an MΨ deviating more than 3-standard deviations from the mean, or (ii) two or more CWWs with MΨ values deviating more than 2-standard deviations from the mean. In order to evaluate any related incidental findings and/or potential underlying etiologies of these aberrant CSF flow patterns, available MRI images of all OASIS study participants were evaluated by a neuroradiologist (A.N.). Common or minor incidental abnormalities such as white matter T₂/FLAIR hyperintensities, chronic lacunar infarctions, and/or small aneurysms (less than 2 mm) were not included in the search pattern.

All region-of-interest statistical analyses were performed with R v4.1.2. (http://www.r-project.org/). To determine the study-specific effect sizes (standardized β) of age, sex, and brain parenchymal fractions (BPF) on CWW MΨ values, general linear models were fitted separately for each of the two OASIS cohorts (Model 1: CWW_MΨ ~ age + sex + BPF). Effects of ventricular volumes and CBF were estimated while accounting for the effects of age, sex, and BPF (Model 2: CWW_MΨ ~ age + sex + BPF + [ventricular volume/CBF]). Results of the two OASIS cohorts were aggregated using linear mixed-effects models (lmer function, part of the lme4 package) with a random intercept for the cohort (Mixed Model: CWW_MΨ ~ age + sex + BPF + (1|Cohort)). For the aggregate analysis, individuals duplicated across subsets were removed from the larger, multi-low b-value OASIS-3 dataset. To mitigate the influence of outliers, we excluded CWW MΨ values falling beyond 1.5 times the interquartile range (IQR) from the analysis. The p-values of the linear mixed models and CBF general linear models were adjusted for multiple testing using the false discovery rate (FDR). FDR-corrected p-values lower than 0.05 were considered significant.

The OASIS-3 dataset included cognitively unimpaired participants as well as individuals with varying degrees of cognitive decline (Table 1). The multi-b-value imaging subset (n = 187) of the OASIS-3 dataset was used as the primary dataset for CSF seed-based correlation analysis and NMF. Single-shell low-b dMRI data from another OASIS-3 imaging subset (OASIS – b100 cohort, n = 103, b-value = 100 s/mm²), the Microstructure-Informed Connectomics (MICA – b300, 49 healthy participants, b-value = 300 s/mm²), and the Comprehensive Diffusion MRI Dataset (CDMD – b50, 26 healthy participants, b-value = 50 s/mm²) were used for replication analyses (Table 1).

A seed-based correlation analysis was conducted in the multi-b-value OASIS-3 subset to examine whether CSF MΨ covariance aligns with the expected similarity in flow characteristics, as driven by shared extracranial and intracranial mechanisms. To prepare for voxel-wise analysis, CSF MΨ maps were aligned to a common space using CSF pseudo-diffusion spatial statistics (see Fig. 1 for an overview of the CΨSS procedure). Five seed regions of interest (ROIs) were selected to explore CSF flow patterns in a hypothesis-driven way. These regions included the left and right Sylvian fissures to assess symmetry, the premedullary cistern to examine craniocervical flow, the cerebral aqueduct-fourth ventricle complex, and the foramina of Monro to analyze intraventricular flow (Fig. 2). CSF MΨ in the left and right Sylvian fissures displayed similar correlation patterns that demonstrated high bilateral symmetry with higher correlation values observed in regions ipsilateral to the seed than contralateral ones. The correlation patterns included the perivascular CSF spaces surrounding the main arteries of the anterior circulation in the ipsilateral and contralateral Sylvian fissures and carotid cisterns, as well as the anterior interhemispheric subarachnoid space. The correlation maps of the cerebral aqueduct and foramina of Monro MΨ were similar, reflecting the common pathway for trans-ventricular CSF flow, including the third and fourth ventricles. Premedullary CSF MΨ was positively correlated with CSF flow in the dorsal and ventral craniocervical junction columns extending to the cerebellopontine angles. Similar seed-based correlation maps were also observed in the b100 imaging cohort (Fig. 2) with Pearson r values between correlation maps generated from the multi-b-value and b100 cohorts ranging from 0.66–0.78.

Having demonstrated that the covariance of CSF MΨ is in accordance with the expected patterns of CSF flow, our next objective was to investigate the organization of intracranial CSF flow using an unbiased, data-driven approach. To achieve this, we applied NMF to the CSF MΨ maps from the OASIS-3 multi-b-value cohort. NMF is a multivariate pattern analysis technique, which has recently been adapted for use in neuroimaging (Nazeri et al., 2022; Sotiras et al., 2015, 2017) and allows for a parts-based representation of complex high-dimensional datasets (Brunet et al., 2004; Lee & Seung, 1999; Wu et al., 2016). We evaluated multiple NMF solutions with varying numbers of patterns (K = 2 to 25) to determine the most reproducible and reliable solution with the lowest reconstruction error. Split-half reproducibility analysis, as quantified using the ARI (Hubert & Arabie, 1985), demonstrated nonuniform reproducibility and reliability across NMF solutions (Fig. 3A). Reproducibility (i.e., high mean ARI across bootstraps) peaked at the 9- and 10-pattern solutions (mean ARI: 0.81 [9-pattern] and 0.82 [10-pattern]). The highest reliability (i.e., high stability with low ARI standard deviation across bootstraps) was also observed for the 9-pattern (SD = 0.017) and 10-pattern (SD = 0.018) solutions (Fig. 3A). The 10-pattern solution, which had a lower reconstruction error, was selected for the subsequent analyses (Fig. 3B).

The CSF waterways (CWWs; the NMF patterns) were symmetrically organized within the ventricular system and basilar cisterns (Fig. 3C). The ventricular system was subdivided into the lateral ventricle CWW1 and the trans-ventricular CWW10. CWW1 also included the extra-axial subarachnoid spaces, whereas the trans-ventricular CWW10 encompassed a small portion of the lateral ventricles adjacent to the foramina of Monro, the third and fourth ventricles, and the cerebral aqueduct. The Sylvian fissures were symmetrically divided into the deep cisternal (CWW8) and superficial opercular compartments (CWW7) (Tayebi Meybodi et al., 2019). Six other CWWs corresponded to the basilar cisterns: (i) cisterna magna, cervicomedullary cistern, and cerebellopontine angles (CWW2); (ii) superior premedullary cistern (CWW3); (iii) chiasmatic cistern (CWW4); (iv) prepontine cistern (CWW5); (v) perimesencephalic cisterns and choroid fissures (CWW6); and (vi) supracerebellar and velum interpositum cisterns (CWW9).

Figure 4 illustrates the structural and functional correlates of CWWs. Trans-ventricular CWW10 is bounded by the foramina of Monro and Magendie, and the interfaces of basilar cistern CWWs align with arachnoid membranes such as the anterior pontine membrane and the Liliequist membrane. The template generated from high-resolution 3D-CISS imaging revealed prominent arachnoid trabeculations within CWW9 and delineated the Liliequist membrane.

To investigate whether the identified CWWs also corresponded to CSF spaces with distinct pseudo-diffusion characteristics (across the population), we extracted the mean CSF MΨ values from each CWW for every participant. The CWWs displayed noticeable differences in their CSF MΨ distributions in the OASIS multi-b-value cohort (Fig. 5A). Specifically, the superior premedullary cistern (CWW3) exhibited the highest MΨ followed by the craniocervical junction (CWW2), the chiasmatic cistern (CWW4), and the prepontine cistern (CWW5). Conversely, the extra-axial/lateral ventricular CSF (CWW1) displayed the lowest MΨ. Among the basilar cisterns, the supracerebellar (CWW9) and perimesencephalic cisterns (CWW6) exhibited the lowest MΨ. Furthermore, the deep cisternal compartment of the Sylvian fissure (CWW8) consistently exhibited higher MΨ than the superficial opercular Sylvian fissure (CWW7). Similar mean MΨ distributions were observed across CWWs in three separate single-shell diffusion-weighted MRI datasets (Fig. 5A) with different low b-values (OASIS – b100, MICA – b300, and CDMD – b50). Notably, the MΨ values showed the lowest dynamic range in MICA – b300 cohort.

We leveraged the densely sampled diffusion MRI data from the CDMD to explore additional pseudo-diffusion characteristics of CSF (Fig. 5B). First, we calculated mean square displacement (MSD) values derived from the multi-shell low-b dMRI (b-values: 50 s/mm², 350 s/mm², 800 s/mm², 32 directions for each shell). MSD is a measure of the distance the water molecules travel during the diffusion time and is closely related to mean diffusivity. As expected, multi-shell MSD values were highly correlated with single-shell MΨ values across the CWWs (Pearson r: 0.83–0.99; p < 0.001), indicating that the use of single-shell low-b dMRI is sufficient for evaluation of pseudo-diffusion magnitude. Next, we assessed the directionality of CSF effective motility in different CWWs by calculating fractional anisotropy (FA) at b = 50 s/mm². Apart from the extra-axial/lateral ventricular CSF CWW1 that showed relatively isotropic pseudo-diffusion (FA: 0.17–0.35), other CWWs showed either high (CWW8) or intermediate anisotropy (Fig. 5B). Unlike MSD, the correlation between CSF FA and MΨ was highly variable across different CWWs (Fig. S1). Finally, to assess temporal variability of CSF effective motility, normalized mean square error (nRMSE) was calculated as the error between the predicted and measured dMRI signal at b = 50 s/mm². The basilar cistern CWWs exhibited the highest nRMSE values, likely reflecting the influence of cardiac pulsatility and/or respiratory cycle on CSF flow dynamics in these regions. The lowest nRMSE was seen in CWW1, suggesting low temporal variability and pulsatile motion in the lateral ventricles and extra-axial subarachnoid spaces. Higher nRMSE values were associated with higher MΨ across the CWWs (r > 0.65; except for CWW4 and CWW5), suggesting that greater average CSF effective motility is linked to greater temporal variability.

The data-driven NMF approach identified patterns of CSF MΨ covariance, but did not include information regarding participant demographics, intracranial anatomy, or cerebral perfusion. Accordingly, in both OASIS-3 cohorts, we examined how age, sex, brain atrophy, ventricular anatomy, and brain perfusion affect CSF MΨ (Fig. 6). Advancing age was associated with higher MΨ values in the Sylvian fissure CWWs (CWW7 and CWW8) and the chiasmatic cistern (CWW4), which are adjacent to the middle cerebral artery branches and circle of Willis, respectively. Women showed higher CSF MΨ values in the caudal basilar cisterns (CWW2 and CWW3).

Brain atrophy (i.e., lower brain parenchymal fraction) and larger ventricle volumes were associated with lower supracerebellar cistern CWW9 MΨ. Larger lateral ventricle volume was also associated with lower MΨ in the extra-axial/lateral ventricular CSF (CWW1) and superficial opercular Sylvian fissure (CWW7). Larger fourth ventricle volume was associated with higher transventricular CWW10 MΨ and lower craniocervical junction CWW2 MΨ.

In the multi-b-value OASIS cohort, we used the available arterial spin labeling (ASL) imaging data to investigate the association between cerebral perfusion and CSF effective motility within CWWs. Higher cerebral blood flow was associated with higher CSF MΨ values in the Sylvian fissure CWWs (CWW7 and CWW8), craniocervical junction (CWW2), and perivenous supracerebellar cistern (CWW9, Fig. 6B). Overall, these findings highlight the regional specificity of the effects of demographics, intracranial anatomy, and cerebral perfusion on CSF flow.

In the OASIS study, individuals were considered outliers if they exhibited (i) a single CWW with an MΨ above or below 3 standard deviations from the mean or (ii) two or more CWWs with MΨ values above or below 2 standard deviations from the mean. Out of thirty-eight participants with outlier CWW MΨ values (multi-b-value cohort: 23; b100 cohort: 15), five individuals had clinically consequential findings (Fig. 7) including: bilateral subdural effusions (IF001); fourth ventricle mass (IF002); multiple intracranial aneurysms (IF003); retinal detachment and new cerebellar and cortical cystic lesions (IF004); and cerebellar arteriovenous malformation (IF005). None of these clinically consequential findings were identified in other participants within the OASIS study subsets, indicating a significant enrichment of clinically consequential imaging findings in participants with aberrant CSF MΨ in the CWWs (Fisher Exact test: p = 3.04 × 10^-5). Three participants with aberrant CSF MΨ in the CWWs exhibited major neurovascular abnormalities, including the previously mentioned two cases (IF003 and IF005) and an additional case of dolichoectasia (IF006); no such abnormalities were observed in other participants (Fisher Exact test: p = 0.0021). Additional incidental findings included mega cisterna magna and cerebellar tonsillar ectopia, both of which were also observed in other participants.

In this work, we demonstrated the intricate interconnectivity of intracranial CSF circulation using a novel voxelwise CSF flow mapping framework and large whole-brain low-b dMRI datasets. We identified distinct spatial patterns of CSF effective motility using NMF, which is an unsupervised multivariate pattern analysis approach. We showed that the identified patterns correspond to distinct CSF pseudo-diffusion characteristics and are differentially impacted by age, sex, ventricular anatomy, and brain perfusion. Our results were replicated in independent datasets with different low-b dMRI protocols. Finally, we showed that individuals with aberrant CSF flow patterns tend to exhibit incidental neuroradiological findings.

Consistent with prior reports (Bito et al., 2021; Taoka et al., 2019, 2021; Wen et al., 2022), we demonstrate that low-b dMRI can rapidly probe CSF effective motility throughout the entire brain within a clinically feasible timeframe (1–3 minutes). Additionally, our study revealed that variations of the MΨ values across different CSF spaces are consistent across a range of imaging protocols, including those with different b-values (ranging from 50 s/mm²–350 s/mm²), single-shell and multi-b-value imaging schemes, various spatial resolutions, and MRI hardware. Further, we showed that other characteristics of CSF flow such as anisotropy and pseudo-diffusion variability can also be assessed using low-b dMRI with higher angular resolution. These findings suggest that low-b dMRI offers a practical, reliable, and versatile method for investigating CSF effective motility and may facilitate diagnosis and monitoring of brain disorders associated with altered CSF flow dynamics.

By providing a sparse and interpretable representation, NMF-derived CWWs enable the study of intracranial CSF flow characteristics while reducing multiple comparisons. Previous literature provides strong evidence that NMF-derived patterns reflect biologically meaningful processes, offering novel insights into brain organization in both health and disease (Earnest et al., 2024; Nazeri et al., 2022). Similarly, our findings suggest that CWWs reflect anatomical and physiological bases of intracranial CSF flow, rather than representing purely statistical constructs. The symmetry of CWWs with non-communicated compartments such as the lateral ventricles, extra-axial subarachnoid spaces, and Sylvian fissures suggests that CSF flow dynamics within these CSF spaces are driven and constrained by similar underlying mechanisms with left-right symmetry. In addition, the CWWs were segregated by functional and structural boundaries (Fig. 4). The trans-ventricular CWW10 was bordered by the foramina of Monro and Magendie, and was further demarcated by fine parenchymal structures such as the lamina terminalis (separating the third ventricle from the lamina terminalis cistern [CWW5]). In some other instances, the interfaces between basilar cistern CWWs aligned with the arachnoid membranes (Lü, 2015; Lu et al., 2022). The interface separating the chiasmatic cistern (CWW4) and prepontine cisterns (CWW5) approximated the sellar segment of the Liliequist membrane. The interface separating the interpeduncular cistern (CWW6) from the prepontine (CWW5) and chiasmatic (CWW4) cisterns followed the mesencephalic and diencephalic leaflets of the Liliequist membrane, respectively. The correspondence between the interfaces of basilar cistern CWWs and arachnoid membranes, such as the Liliequist membrane (Lu et al., 2022), underscores the influence of arachnoid matter’s delicate anatomy on shaping intracranial CSF flow patterns. Furthermore, supracerebellar and velum interpositum cisterns (CWW9) represent a distinct fine internal anatomy with pronounced arachnoid trabeculations. These findings further reinforce the notion that CWWs are not mere statistical artifacts but rather shaped by anatomical structures partially compartmentalizing the CSF spaces and reflect critical physiological processes governing CSF circulation.

Previous large-scale CSF flow studies have predominantly utilized phase-contrast MRI techniques to study CSF flow in the cerebral aqueduct (Rohilla et al., 2024; Sartoretti et al., 2019). Therefore, effects of demographics and intracranial anatomy on CSF flow in other brain regions have not yet been established. Our study addresses this gap by comprehensively examining the effects of demographics and intracranial anatomy on CSF flow across several brain regions in a relatively large sample size (n = 264 unique participants in the OASIS-3 cohorts). Our findings demonstrate an association between increased fourth ventricle volumes and higher trans-ventricular (CWW-10) MΨ values, corroborating previous phase-contrast MRI studies linking larger fourth ventricular volumes to higher aqueductal CSF flow velocities (Chiang et al., 2009; Magnano et al., 2012). However, our results diverge from these studies by not showing a similar increase in CSF flow with the enlargement of third and lateral ventricles or age (Rohilla et al., 2024; Sartoretti et al., 2019). The apparent discrepancies in our results could be attributed to the inclusion of brain parenchymal fraction—a measure of brain atrophy—in our regression models, suggesting that previously reported age effects likely reflect the influence of age-related brain atrophy rather than aging itself. Excluding this adjustment, we found associations between CWW-10 MΨ, age, and lateral and third ventricle volumes similar to those reported in prior studies (Fig. S2). These results indicate that, in contrast to the lateral and third ventricles, the association between fourth ventricle volume and trans-ventricular CSF flow is independent of brain atrophy.

We found that the ventricular system could be functionally divided into low isotropic CSF flow in the lateral ventricles (CWW1) and transventricular anisotropic and variable CSF flow (CWW10) encompassing the frontal horns of the lateral ventricles adjacent to the foramina of Monro, third ventricle, cerebral aqueduct, and the fourth ventricle. Moreover, these intraventricular CSF flow patterns were influenced by distinct intracranial anatomy. In the lateral ventricles, larger volumes were associated with lower CSF MΨ. This may indicate greater dissipation of the CSF flow within larger lateral ventricle volumes. In contrast, our results indicate that larger fourth ventricular volumes are associated with higher transventricular CSF MΨ, which is in concordance with the diffuser/nozzle pump model, which posits that fourth ventricular morphology affects CSF flow in the third ventricle–aqueduct–fourth ventricle complex (Longatti et al., 2019). Of note, these relationships, primarily reflecting normal variations in ventricle volumes, could potentially undergo changes in the setting of disease conditions such as hydrocephalus.

The association between cerebral blood flow and CSF MΨ in the caudal basilar cisterns likely reflects the coupling between the fluctuations in the cerebral blood volume and intracranial CSF volume throughout the cardiac cycle, as asserted by the Monro-Kellie doctrine (Mokri, 2001). Our results also show that higher cerebral blood flow is associated with increased MΨ in the peri-arterial CSF spaces. The peri-arterial CSF flow is propelled by arterial pulsations (Kedarasetti et al., 2020; Mestre et al., 2018), which could be amplified by higher cerebral blood flow. The observed age-related increase in periarterial CSF MΨ aligns with recent findings of elevated peri-arterial CSF flow and pulsatility in aging subjects (Wen et al., 2022). This rise in peri-arterial CSF flow is likely driven by increased intracranial arterial pulsatility with aging (Zarrinkoob et al., 2016).

By reviewing multimodal images of the OASIS study participants with aberrant CSF flow dynamics, we were able to identify various incidental findings ranging from vascular anomalies to mass lesions and intracranial fluid collections. In some instances, a direct causal link could be established between the observed CSF flow abnormality and the incidental findings. For instance, high CWW1 MΨ associated with subdural effusions suggests higher flow within the subdural collections compared to the extra-axial subarachnoid space, particularly adjacent to the bridging veins and arachnoid membrane. Taken together, our results highlight the importance of considering CSF dynamics in the context of a broad range of brain disorders and suggest that our approach may have implications for the diagnosis and treatment of a wide range of brain disorders.

There are some limitations that need to be considered. While in this study we capitalized on a large number of low-b dMRI images to evaluate the patterns of CSF flow and the factors influencing them, the low-b dMRI used in these datasets were acquired as a part of multi b-value datasets. Future studies can exploit the inherent high signal-to-noise ratio of low-b dMRI images permit CSF flow imaging with high spatial resolution with minimal image distortion. Such high-fidelity, high-spatial-resolution low-b dMRI is particularly crucial for evaluating relatively narrow sulcal subarachnoid spaces and brain areas adjacent to air-tissue interfaces. Moreover, in this study, the low-b dMRI images were not cardiac gated. Recent research has shown that gated low-b dMRI allows for characterizing the dynamics of CSF flow over the cardiac cycle (Wen et al., 2022). Of note, intrinsic physical properties of CSF such as temperature and viscosity could also affect the observed MΨ values (Hasan et al., 2015). This effect is likely negligible in a healthy population. However, caution is warranted in the interpretation of MΨ variations in pathological states that may alter physical properties of CSF (e.g., meningitis or hemorrhage) (Vyas et al., 2024; Yetkin et al., 2010). In our study, we examined the influence of demographics, intracranial anatomy, and cerebral perfusion on CSF dynamics. To further elucidate the underpinnings of intracranial CSF flow, future investigations need to consider additional factors such as brain compliance and viscoelastic properties (Gholampour et al., 2022). Finally, while our results provide valuable insights into CSF flow under normal conditions, it is underpowered to delineate the relationship between CSF flow and various brain diseases such as age-related neurodegenerative disorders and dementia. This is in part due to the small number of individuals with cognitive impairment in the OASIS study subsets, the inherent heterogeneity of dementias, and the significant variance in CSF flow that requires adjustment for demographic, structural, and physiological factors. Large-scale, integrative studies are necessary to explore potential CSF flow abnormalities in neurodegenerative disorders.

In conclusion, our work sets forth a new paradigm to study CSF flow by capitalizing on large sample size low-b dMRI datasets and unsupervised machine learning. Our study provides insights in the CSF flow organization, revealing distinct CSF waterways, which exhibit distinct CSF pseudo-diffusion characteristics and are differentially impacted by age, sex, ventricular anatomy, and brain perfusion. The data-driven CWW-based atlas of the intracranial CSF flow offers a systematic, parts-based approach to examining intracranial CSF flow patterns under physiological conditions and their abnormalities in the setting of various brain CSF disorders such as hydrocephalus, idiopathic intracranial hypertension, or Chiari 1 malformation.

The OASIS-3 data that were used in this work are available online (https://www.nitrc.org/projects/oasis). The MICA dataset is available on the Canadian Open Neuroscience Platform’s data portal (https://portal.conp.ca). CDMD imaging dataset is available as a figshare collection (https://doi.org/10.6084/m9.figshare.c.5315474.v1). The implementation of the non-negative matrix factorization used in this study is available at https://github.com/asotiras/brainparts.

A.N., T.L.S.B., and A.S. conceived and designed the study and analysis; A.N., H.H., P.L., and T.L.S.B., acquired the data; A.N., H.H., and T.D. analyzed the data; and A.N. and A.S. wrote the first draft of the paper. A.N., H.H., T.D., K.E.L., P.L., J.S.S., T.L.S.B., and A.S. provided suggestions and revisions for the final draft.

The authors declare no competing interests.

A.N. was partially supported by the National Institutes of Health (NIH) award P01NS131131. A.S. was partially supported by the NIH award R01AG067103. The facilities of the Washington University Center for High Performance Computing (CHPC) partially funded by NIH grants 1S10RR022984-01A1 and 1S10OD018091-01 were used for computational analyses.

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