Functional MRI brain state occupancy in the presence of cerebral small vessel disease—a pre-registered replication analysis of the Hamburg City Health Study

Cerebral small vessel disease (cSVD) is an arteriolopathy of the brain associated with age and common cardiovascular risk factors (Wardlaw, C. Smith, & Dichgans, 2013). cSVD predisposes patients to ischemic stroke (in particular lacunar stroke) and may lead to cognitive impairment and dementia (Cannistraro et al., 2019). Neuroimaging findings in cSVD reflect its underlying pathology (Wardlaw et al., 2015) and include white matter hyperintensities (WMH), lacunes of presumed vascular origin, small subcortical infarcts and microbleeds, enlarged perivascular spaces, as well as brain atrophy (Wardlaw, E. E. Smith, Biessels, et al., 2013). However, the extent of visible cSVD features on magnetic resonance imaging (MRI) is an imperfect predictor of the severity of clinical sequelae (Das et al., 2019) and our understanding of the causal mechanisms linking cSVD-associated brain damage to clinical deficits remains limited (Bos et al., 2018).

Recent efforts have focused on exploiting network aspects of the structural (Tuladhar, et al., 2016, 2020; Lawrence, Zeestraten, et al., 2018) and functional (Dey et al., 2016; Schulz et al., 2021) organization of the brain to understand the relationship between cSVD and clinical deficits in cognition and other domains that rely on distributed processing. Reduced structural network efficiency has repeatedly been described as a causal factor in the development of cognitive impairment, particularly executive dysfunction and reduced processing speed in cSVD (Lawrence, Chung, et al., 2014; Shen et al., 2020; Reijmer et al., 2016; Prins et al., 2005). Findings with respect to functional connectivity (FC), however, are more heterogeneous than their SC counterparts, perhaps because FC measurements are prone to be affected by hemodynamic factors and noise, resulting in relatively low reliability, especially with resting-state scans of short duration (Laumann et al., 2015). This problem is exacerbated in the presence of cSVD and worsened by arbitrary processing choices (Gesierich et al., 2020; Lawrence, Tozer, et al., 2018).

As a promising new avenue, time-varying, or dynamic, functional connectivity approaches have recently been explored in patients with subcortical ischemic vascular disease (Xu et al., 2021; Yin et al., 2022). Although the study of dynamic FC measures may not solve the problem of limited reliability, especially in small populations or participants with extensive structural brain changes, it adds another—temporal—dimension to the study of functional brain organization, which is otherwise overlooked. Importantly, FC dynamics not only reflect moment-to-moment fluctuations in cognitive processes, but are also related to brain plasticity and homeostasis (Laumann & Snyder, 2021; Laumann, et al., 2017), which may be impaired in cSVD.

In the present paper, we aimed to replicate and extend the main results of Schlemm et al. (2022). In this recent study, the authors analyzed MR imaging and clinical data from the prospective Hamburg City Health Study (HCHS, Jagodzinski et al., 2020) using a coactivation pattern approach to define discrete brain states, and found associations between the WMH load, time spent in high-occupancy brain states characterized by activation or suppression of the default mode network (DMN), and cognitive impairment. Specifically, every 4.7-fold increase in WMH volume was associated with a 0.95-fold reduction in the odds of occupying a DMN-related brain state; every 2.5 seconds (i.e., one repetition time) not spent in one of those states was associated with a 1.06-fold increase in TMT-B completion times.

The fractional occupancy of a functional MRI-derived discrete brain state is a participant-specific measure of brain dynamics and is defined as the proportion of BOLD volumes assigned to that state relative to all BOLD volumes acquired during a resting-state scan.

Our primary hypothesis for the present work was that the volume of supratentorial white matter hyperintensities is associated with fractional occupancy of DMN-related brain states in a middle-aged to elderly population mildly affected by cSVD. Our secondary hypothesis was that fractional occupancy is associated with executive dysfunction and reduced processing speed, measured as the time to complete part B of the Trail Making Test (TMT).

Both hypotheses were tested in an independent subsample of the HCHS study population using the same imaging protocols, examination procedures, and analysis pipelines as those in Schlemm et al. (2022). The robustness of the associations was explored using a multiverse approach by varying key steps in the analysis pipeline.

This study analyzed data from the Hamburg City Health Study (HCHS), an ongoing prospective, population-based cohort study aiming to recruit a cross-sectional sample of $45000$ adult participants from the city of Hamburg, Germany (Jagodzinski et al., 2020). From the first $10000$ participants of the HCHS, we planned to include those who were documented to have received brain imaging (n = 2648) and exclude those who were analyzed in our previous report (Schlemm et al., 2022) (n = 970). The ethical review board of the Landesärztekammer Hamburg (State of Hamburg Chamber of Medical Practitioners) approved the HCHS (PV5131), and all participants provided written informed consent.

From the study database, we extracted the participants’ age at the time of inclusion in years, their sex, and the number of years spent in education. During the visit to the study center, participants underwent cognitive assessment using standardized tests. From the database, we extracted their performance scores on the Trail Making Test part B, measured in seconds, as an operationalization of executive function and psychomotor processing speed (Arbuthnott & Frank, 2000; Tombaugh, 2004). For descriptive purposes, we also extracted data on past medical history and reported the proportion of participants with a previous diagnosis of dementia.

The magnetic resonance imaging protocol for the HCHS includes structural and resting-state functional sequences. The acquisition parameters for a $3 T$ Siemens Skyra MRI scanner (Siemens, Erlangen, Germany) have been previously reported (Frey et al., 2021; Petersen et al., 2020) and are given as follows:

For $T1$-weighted anatomical images, a 3D rapid acquisition gradient-echo sequence (MPRAGE) was used with the following sequence parameters: repetition time $TR=2500ms$⁠, echo time $TE=2.12ms$⁠, 256 axial slices, slice thickness $ST=0.94mm$⁠, and in-plane resolution $IPR=(0.83×0.83)mm2$⁠.

$T2$-weighted fluid attenuated inversion recovery (FLAIR) images were acquired with the following sequence parameters: $TR=4700ms$⁠, $TE=392ms$⁠, $192$ axial slices, $ST=0.9mm$⁠, $IPR=(0.75×0.75)mm2$⁠.

$125$ resting state functional MRI volumes were acquired (⁠$TR=2500ms$⁠; $TE=25ms$⁠; $flip angle=90°$⁠; slices = $49$⁠; $ST=3mm$⁠; $slice gap=0mm$⁠; $IPR=(2.66×2.66) mm2$⁠). The participants were asked to keep their eyes open and to think of nothing.

We verified the presence and voxel dimensions of expected MRI data for each participant and excluded those for whom at least one of $T1$-weighted, FLAIR, and resting-state MRI was missing. We also excluded participants with neuroradiologically confirmed space-occupying intra-axial lesion. To ensure reproducibility, no visual quality assessment of raw images was performed.

For the remaining participants, structural and resting-state functional MRI data were preprocessed using FreeSurfer v6.0 (https://surfer.nmr.mgh.harvard.edu/), and fMRIPrep v20.2.6 (Esteban et al., 2019), using default parameters. Participants were excluded if automated processing using at least one of these packages failed.

For our primary analysis, the extent of ischemic white matter disease was operationalized as the total volume of supratentorial WMH obtained from automated segmentation using a combination of anatomical priors, BIANCA (Griffanti et al., 2016), and LOCATE (Sundaresan et al., 2019), post-processed with a minimum cluster size of $30$ voxels, as described in Schlemm et al. (2022). In an exploratory analysis, we partitioned voxels identified as WMH into deep and periventricular components according to their distance to the ventricular system (cut-off $10mm$⁠, Griffanti et al., 2018).

The output from fMRIprep was post-processed using xcpEngine v1.2.3 to obtain de-confounded spatially averaged BOLD time series (Ciric et al., 2017). For the primary analysis, we used the 36p regression strategy and the Schaefer-$400$ parcellation (Schaefer et al., 2018), as in Schlemm et al. (2022).

Different atlases and confound regression strategies, as implemented in xcpEngine, were included in an exploratory multiverse analysis.

Co-activation pattern (CAP) analysis was performed by first aggregating parcellated, de-confounded BOLD signals into a $(nparcels×∑intime points,i)$ feature matrix, where $ntime points,i$ denotes the number of retained volumes for participant $i$ after confound regression. Clustering was performed using the $k$-means algorithm (⁠$k=5$⁠) with a distance measure given by 1 minus the sample Pearson correlation between points, as implemented in Matlab R2021a. We estimated the participant- and state-specific fractional occupancies, which are defined as the proportion of BOLD volumes assigned to each brain state (Vidaurre et al., 2018). The two states with the highest average occupancies were identified as the basis for further analysis.

For demographic (age, sex, and years of education) and clinical (TMT-B) variables, the number of missing items is reported. For non-missing values, we provide descriptive summary statistics using median and interquartile range. The proportions of men and women in the sample are reported. Since we expected based on our pilot data (Schlemm et al., 2022) that the proportion of missing data would be small, primary regression modelling was carried out as a complete-case analysis.

As an outcome-neutral quality check of the implementation of the MRI processing pipeline, brain state estimation, and co-activation pattern analysis, we compared fractional occupancies between brain states. We expected that the average fractional occupancy in the two high-occupancy states would be higher than the average fractional occupancy in the other three states. Point estimates and 95% confidence intervals are presented for the difference in average fractional occupancy to verify this assertion.

For further analyses, non-zero WMH volumes were subjected to logarithmic transformation. Zero values retained their value of zero; to compensate, all models included a binary indicator for zero WMH volume if at least one non-zero WMH value was present.

To assess the primary hypothesis of a negative association between the extent of ischemic white matter disease and time spent in high-occupancy brain states, we performed a fixed-dispersion Beta regression to model the logit of the conditional expectation of the average fractional occupancy of two high-occupancy states as an affine function of the logarithmized WMH load. Age and sex were included as covariates. The strength of the association was quantified as the odds ratio per interquartile ratio of the WMH burden distribution, and is accompanied by a 95% confidence interval. Significance testing of the null hypothesis of no association was conducted at the conventional significance level of 0.05. Estimation and testing were carried out using the “betareg” package v3.1.4 in R v4.2.1.

To assess the secondary hypothesis of an association between time spent in high-occupancy brain states and executive dysfunction, we performed a generalized linear regression with a Gamma response distribution to model the logarithm of the conditional expected completion time in part B of the TMT as an affine function of the average fractional occupancy of two high-occupancy states. Age, sex, years of education, and logarithmized WMH load were included as covariates. The strength of the association was quantified as a multiplicative factor per percentage point and accompanied by a 95% confidence interval. Significance testing of the null hypothesis of no association was conducted at the conventional significance level of 0.05. Estimation and testing were performed using the glm function included in the “stats” package from R v4.2.1.

The analysis plan was pre-registered on June 27, 2023 at https://osf.io/fcqmb. The sample size calculation was based on an effect size on the odds ratio scale of 0.95, corresponding to an absolute difference in the probability of occupying a DMN-related brain state between the first and third WMH-load quartile of 1.3 percentage points, and between the 5% and 95% percentile of 3.1 percentage points. Approximating half the difference in fractional occupancy of DMN-related states between different task demands (rest vs n-back) in healthy participants, which was estimated to lie between 6 and 7 percentage points (Cornblath et al., 2020), this value represented a plausible choice for the smallest effect size of theoretical and practical interest. It also equals the estimated effect size based on the data presented in Schlemm et al. (2022).

Simple bootstrapping was used to create $10000$ hypothetical datasets of size $200$⁠, $400$⁠, $600$⁠, $800$⁠, $900$⁠, $910$⁠,… , $1090$⁠, $1100$⁠, $1200$⁠, $1400$⁠, $1500$⁠, and $1600$⁠. Each dataset was then subjected to the estimation procedure described above. For each sample size, the proportion of datasets in which the primary null hypothesis of no association between fractional occupancy and WMH load could be rejected at $α=0.05$ was computed and recorded as a power curve in Figure 1.

A sample size of $960$ would have allowed the replication of the reported effect with a power of $80.2%$⁠. We had anticipated a sample size of $1500$⁠, which would have yielded a power of $93.8%$⁠.

In both (Schlemm et al., 2022) and our primary replication analysis, we made certain analytical choices in the operationalization of brain states and ischemic white matter disease, namely the use of the 36p confound regression strategy, the Schaefer-$400$ parcellation, and a BIANCA/LOCATE-based WMH segmentation algorithm. The robustness of the association between WMH burden and time spent in high-occupancy states with regard to other choices was explored in a multiverse analysis (Steegen et al., 2016). Specifically, in an exploratory analysis, we estimated brain states from BOLD time series processed according to a variety of established confound regression strategies and aggregated over different cortical brain parcellations (Table 2, Ciric et al., 2017, 2018). The extent of cSVD was additionally quantified by the volume of deep and periventricular white matter hyperintensities.

For each combination of analytical choice of confound regression strategy, parcellation, and subdivision of white matter lesion load (⁠$9×9×3=243$ scenarios in total), we quantified the association between WMH load and average time spent in high-occupancy brain states using odds ratios and $95%$ confidence intervals as described above.

No hypothesis testing was performed for these multiverse analyses. Rather, they serve to inform about the robustness of the outcome of the test of the primary hypothesis. Any substantial conclusions about the association between the severity of cerebral small vessel pathology and the time spent in high-occupancy brain states were drawn from the primary analysis using pre-specified methodological choices, as stated in the Scientific Question in Table 1.

In previous work, two high-occupancy brain states have been related to the default mode network (Cornblath et al., 2020). We further explored this relationship by computing, for each individual brain state, the cosine similarity of the positive and negative activations of the cluster’s centroid with a set of a priori defined functional ‘communities’ or networks (Schaefer et al., 2018; Yeo et al., 2011). The results were visualized as spider plots for the Schaefer atlases.

In further exploratory analyses, we describe the associations between brain state dynamics and other measures of cognitive ability such as memory and language.

Summary data from the first $1000$ imaging data points of the HCHS have been published with Schlemm et al. (2022) and formed the basis for the hypotheses tested in this replication study. Before preregistration, we had implemented our prespecified analysis pipeline described above in R and Matlab, and applied it to this previous sample. Data, code, and results from this pilot analysis have been stored with the archived Stage 1 report on GitHub (https://github.com/csi-hamburg/HCHS-brain-states-RR, v1.5) and preserved on Zenodo.

At the time of planning of this study, all demographic, clinical, and imaging data used in this analysis had been collected by the HCHS and were held in the central trial database. Quality checks for non-imaging variables had been performed centrally. WMH segmentation based on structural MRI data of the first $10000$ participants of the HCHS had been performed previously using the BIANCA/LOCATE approach (Rimmele et al., 2022). Functional MRI data and clinical measures of executive dysfunction (TMT-B scores) had not previously been analyzed by the pre-registering author (ES).

For deconfounding and aggregating BOLD data at brain parcellation level, the software xcpEngine was used in version 1.2.3, not 1.2.1, to ensure that that the correct MNI reference template (MNI152NLin2009cAsym) is used for registration of brain atlases. This decision was made before analyzing the data.

For this replication study, a total of $2648$ datasets were available, of which $970$ were already included in our previous analysis and thus discarded. In $13$ of the resulting $1678$ datasets, one or more MRI sequences were missing. Of the complete datasets (n = $1665$⁠), we excluded $5$ participants due to intra-axial space-occupying lesions. An additional $9$ participants were excluded because of unsuccessful preprocessing, WMH segmentation, or xcpEngine failure, resulting in $1651$ datasets for analysis. A study-flowchart is provided in Figure 2.

Baseline demographic and cognitive values, including the number of missing items, are reported in Table 3.

WMH volumes (median $1.05$ $mL$⁠, IQR $0.47$ $mL$ to $2.37$ $mL$⁠), motion estimates, and fractional occupancies of brain states 1 through 5 are reported in Table 4.

In an outcome-neutral quality check of the implementation of (i) the MRI processing pipeline, (ii) brain state estimation, and (iii) co-activation pattern analysis, the mean difference in fractional occupancy between high- and low-occupancy states was consistently maintained, with a point-estimate of the separation between two high-occupancy and three low-occupancy states of $6.7%$ (⁠$95%$ confidence interval, $6.2%$ to $7.1%$⁠) in the 36p paradigm. This indicates that the implementation of the pipeline was correct and that the brain state estimation and co-activation pattern analysis worked as intended.

The results of the test of our primary preregistered hypothesis of an association between supratentorial WMH volume and the time spent in high-occupancy brain states are shown in Figure 3 and Table 5.

Adjusted for age and sex, there was a $0.94$-fold reduction in the odds of occupying a high-occupancy brain state for every $5.1$-fold increase in WMH load (P = $5.01×10−8$⁠).

The results of the test of our secondary preregistered hypothesis of an association between time spent in high-occupancy brain states and executive function as measured by the complete part B of the TMT are shown in Figure 4 and Table 6.

Adjusted for age, sex, WMH volume, and years of education, there was a $0.98$-fold reduction in the time to complete the TMT-B for every $5%$ increase in the time spent in high-occupancy brain states (P = $0.0116$⁠).

In a multiverse analysis, the main findings of associations between WMH load and FO and, to a lesser extent, between FO and TMT-B were robust with respect to the processing choices of brain parcellation and confound regression strategy.

A nominally statistically significant negative association between the total WMH load and time spent in high-occupancy states was observed in 48 out of 81 scenarios, with 8 out of 81 significant positive associations occurring with the Desikan–Killiany parcellation only (Fig. 5A). For periventricular (deep) WMH volume, the results were similarly robust with 49 out of 81 (39/81) negative and 8 out of 81 (0/81) positive associations of nominal statistical significance, respectively (Appendix Figures 1 and 2).

The secondary finding of an association between greater TMT-B times and lower fractional occupancy was less robust with only 16 out of 81 nominally statistically significant negative and no significant positive associations, irrespective of operationalization of cSVD (total vs. periventricular vs. deep WMH volume) (Fig. 5B, Appendix Figures 1 and 2).

Based on the cosine similarity between positive and negative activations of cluster centroids and indicator vectors of pre-defined large-scale brain networks, network activation profiles were computed for brain states estimated from Schaefer parcellations of varying spatial resolutions.

Figure 6 shows the corresponding spider plots, identifying states characterized by activation (DMN+) or suppression (DMN-) of the default mode network as states with the highest fractional occupancy.

Associations between the time spent in high-occupancy DMN-related brain states and cognitive measures beyond TMT-B are shown in Figure 7.

Adjusted for age, sex, WMH load, and years of education, FO in DMN-related states appeared to be associated with better word recall (adjusted OR 1.19, nominal P 0.013), but not with global cognitive functioning (MMSE, adjusted OR 1.09) or vocabulary (aOR 1.09), nor with verbal fluency (animal naming, adjusted $exp(β)$ 1.04), or pure processing speed (TMT-A, adjusted $exp(β)$ 0.97).

In this pre-registered cross-sectional study, we replicated the key findings of Schlemm et al. (2022) in an independent population-based sample of $1651$ middle-aged to elderly participants of the Hamburg City Health Study.

First, we confirmed that the severity of cerebral small vessel disease is associated with the time spent in high-occupancy brain states, defined by functional MRI. More precisely, we showed that every 5.1-fold increase in the volume of supratentorial white matter hyperintensities of presumed vascular origin (WMH) was associated with a 0.95-fold reduction in the odds of occupying a brain state characterized by activation or suppression of the default-mode network, at any given time during the resting-state scan.

Second, we confirmed that the time spent in high-occupancy brain states at rest is associated with cognitive performance. More precisely, a 5%-reduction in the fractional occupancy of DMN-related brain states was associated with a $1.02$-fold increase in the time to complete part B of the trail making test (TMT).

In a pre-planned multiverse analysis, findings relating to our primary and, to a lesser extent, secondary hypotheses were robust with respect to variations in brain parcellations and confound regression strategies. Inconsistent results were found with the Desikan–Killiany parcellation, likely reflecting the notion that the spatial resolution and functional specificity of this coarse, structurally defined atlas are inadequate for analyzing functionally defined brain states. Across brain parcellations, effect sizes were smaller with the ICA-AROMA confound regression strategy and failed to reach nominal statistical significance. This might be due to a relatively large residual motion component in measures of dynamical functional connectivity after de-noising with ICA-AROMA, as described previously (Lydon-Staley et al., 2019).

We also confirmed across several brain parcellation resolutions that high-occupancy states at rest are characterized by either activation or suppression of the default mode network, reflecting its role as the predominant task-negative brain network.

In unplanned, exploratory analyses, we described the association between brain state dynamics and cognitive measures other than executive function and processing speed and reported a strong, preliminary association between time spent in high-occupancy states and delayed word recall.

We further explored, and report in Appendix Tables 1 and 2, the effect of motion; results relating to our primary and, to a lesser extent, secondary, hypotheses were robust to additional, unplanned adjustments for DVARS, RMSD, and mean framewise displacement.

The presented results provide robust evidence for a behaviorally relevant association between cerebral small vessel disease and functional brain network dedifferentiation.

Further research is required to replicate our findings in different populations, such as those affected more severely by cSVD or cognitive impairment, or being studied using different imaging protocols, to determine the generalizability of our findings with respect to varying operationalizations of the notions of cSVD, brain state, and cognition, and to understand the mechanisms underlying the reported associations.

Preprocessed data and code are available at https://github.com/csi-hamburg/HCHS-brain-states-RR, v2.1.1. https://doi.org/10.5281/zenodo.7337577.

Conceptualization: T.I., B.C., G.T., and E.S.; Data curation: C.M., M.P., B.M.F., and E.S.; Formal analysis: T.I., M.P., B.M.F., and E.S.; Funding acquisition: B.C., G.T., and E.S.; Investigation: T.I., C.M., U.H., C.G., B.C., G.T., and E.S.; Methodology: T.I., B.C., G.T., and E.S.; Project administration: J.F., U.H., S.K., J.G., R.T., C.G., B.C., and G.T.; Resources: J.F., J.G., R.T., C.G., B.C., and G.T.; Software: T.I., M.P., B.M.F., and E.S.; Supervision: B.C., G.T., and E.S.; Validation: T.I., M.P., and E.S.; Visualization: T.I., M.P., and E.S.; Writing—original draft: T.I. and E.S.; Writing—review & editing: T.I., C.M., M.P., B.M.F., J.F., U.H., S.K., J.G., R.T., C.G., B.C., G.T., and E.S.

The authors of this article declare that they have no financial conflict of interest with the content of this article.

This preprint was created using the LaPreprint template (https://github.com/roaldarbol/lapreprint) by Mikkel Roald-Arbøl.

Here we present, in analogy to Figure 5, the results of the multiverse analyses of the association between cSVD burden, FO of DMN-related states, and executive function, when cSVD is operationalized as the volume of deep or periventricular white matter hyperintensities, respectively.

We also present, in analogy to Tables 5 and 6, regression tables for the association between time spent in DMN-related brain states (FO) and WMH volume, and between TMT-B and FO, adjusted for DVARS, RSMD, and framewise displacement, in addition to age, sex, and, in the latter case, years of education.