The biological role of local and global fMRI BOLD signal variability in multiscale human brain organization

Domain: neuroscience
Published: 2026-01-30
Editorial outcome: minor_revision

Abstract
Variability drives the organization and behavior of complex systems, including the human brain. Understanding the variability of brain signals is thus necessary to broaden our window into brain function and behavior. Few empirical investigations of macroscale brain signal variability have been undertaken, given the difficulty in separating biological sources of variance from artefactual noise. Here, we characterize the temporal variability of the most predominant macroscale brain signal, the fMRI BOLD signal, and systematically investigate its statistical, topographical, and neurobiological properties. We contrast fMRI acquisition protocols, and integrate across histology, microstructure, transcriptomics, neurotransmitter receptor and metabolic data, fMRI static connectivity, and empirical and simulated magnetoencephalography data. We show that BOLD signal variability represents a spatially heterogeneous, central property of multi-scale multi-modal brain organization, distinct from noise. Our work establishes the biological relevance of BOLD signal variability and provides a lens on brain stochasticity across spatial and temporal scales.

Introduction
Variability is ubiquitous in our environment. It is a crucial feature of complex ecological and biological systems1: from day-to-day climate variability driving climate change and climate mitigation efforts2,3, to heart rate variability serving as a clinical tool to predict overall cardiovascular health and mortality4. Acting as a catalyst for system adaptability, variability determines the organization of complex systems, defines their spatiotemporal properties and guides their behavior5. Gaining insights into the variability of a system may thus unlock a deeper understanding of the system from which variability emerges.

The human brain is a complex system. From the processing of incoming information to the generation of motor outputs, variability is present at all levels of the central nervous system6,7. Brain activity varies across multiple timescales, from milliseconds to years, and these temporal changes can be observed across multiple spatial scales, from regions to networks8,9. Despite its pervasiveness, there continues to be resistance to exploring variability in human cognitive neuroscience. Human behavior has been shown to be stochastic7,10,11, and computational models have been operationalized the brain as a complex dynamical system12–16, yet corresponding empirical neuroimaging research still lags behind. Most functional MRI (fMRI) investigations are centered on static methodological approaches to brain function. Given the richness of information present in the fMRI blood oxygen level dependent (BOLD) signal, disentangling intrinsic biological sources of signal variance from extrinsic artefactual sources arising from the scanner has historically been a challenge17. fMRI is one of the most widely used and clinically tractable tools for the exploration of macroscale brain function. These realities position a systematic evaluation of BOLD signal variability as a research imperative, necessary to broaden our window into human brain function.

Thus far, BOLD signal variability has been primarily studied in relation to behavior, cognition, development, and clinical status18–24. Demonstrating such associations does not, however, directly establish BOLD signal variability’s biological role in whole-brain organization. As such, it remains unclear whether and to what extent BOLD signal variability investigations are capturing meaningful neural signatures or system noise. To this end, a few recent studies have begun to characterize BOLD variability’s neurobiological underpinnings, yet these investigations are typically restricted to comparisons across age and focus on single, mostly macroscale, neurobiological features25–29. A broader and more integrated assessment of BOLD signal variability in younger samples is therefore needed, one that examines its statistical, topographical and neuronal properties and systematically relates BOLD signal variability to a wider array of neurobiological features.

Evidence in favor of taking such an integrative multi-scale approach to BOLD signal variability comes from multiple lines of empirical and computational work collectively demonstrating the role of cellular, molecular, genetic, and metabolic factors in shaping local and global hemodynamic signals30–34.

First, BOLD signal variability is a statistical approximation of macroscale brain signal dynamics thus, both the estimation method and the properties of the fMRI data from which it is estimated will impact outcomes. Second, human brain function is organized along hierarchical modules, ranging from local functional units to global functional networks35,36, with heterogeneous topographies. Such organization enables information to be both segregated (local) and integrated (global) across units, subserving different neural and cognitive functions—as observed in various measures, such as mean functional activation, functional connectivity, computational, and BOLD signal complexity studies37–42. Along these lines, BOLD signal variability must therefore also be understood within a local-global framework43,44 that considers both the methodology and analysis level of interest. While previous fMRI reports have investigated the local-global nature of BOLD signal variability, they have generally focused on static global measures of brain organization (i.e., the correlation between local mean BOLD signals over time, functional connectivity45) and more local characterizations of variability on a single scale (e.g., voxel or region). The temporal variability of individual regions (local BOLD signal variability) must thus be evaluated alongside the dynamic coordination of functional neural units (dynamic functional connectivity, here global BOLD signal variability). If biologically relevant, local and global BOLD signal variability should present spatially heterogeneous topographies that recapitulate known neurobiological processes that unfold across spatial scales. Third, BOLD signal variability should reflect aspects of macroscale neuronal signals that occur at finer temporal scales and are captured by neuroscientific modalities with greater temporal resolution than fMRI, such as electrophysiology.

In this study, we defined measures of local and global BOLD signal variability as the regional moment-to-moment change in BOLD signal intensity between successive timepoints and as the similarity in inter-regional functional connectivity over time. To determine robust results dissociable from noise, we leveraged multiple state-of-the-art openly available fMRI datasets and assessed the validity and reliability of local and global BOLD variability across samples. To understand the spatial organization and biological properties of local and global BOLD variability, inspired by recent multiscale fMRI investigations32,46–49, we examined their topography within each fMRI dataset and interrogated associations with open-source data, including ex vivo histology50 and in vivo microstructure51, transcriptomics52,53, PET-derived neurotransmitter receptor and metabolic information47, and fMRI static connectivity data54,55. Finally, we leveraged magnetoencephalography (MEG) data and naturalistic electrophysiological simulations to mechanistically understand the temporal and neuronal properties of local multimodal signal variability. We found that measures of BOLD signal variability exhibited spatially heterogeneous topographies, were embedded within multi-scale brain organization, and were rooted in electrophysiological processes. Together, our work establishes local and global BOLD signal variability as biologically relevant, central features of multi-scale, multi-modal brain organization, distinct from noise. Our study represents a step forward towards understanding brain stochasticity across spatial and temporal scales.

Results
Quantification of local and global BOLD signal variability
We first sought to identify robust regional metrics of local and global BOLD signal variability. By local and global BOLD signal variability, we refer to changes in the dynamic properties of single time series (local) and of pairs of time series (global). The former case is what is traditionally referred to as “BOLD signal variability” in the literature, and the latter as “dynamic functional connectivity”. In using the general term “variability” for both measures, we aimed at conceptually juxtaposing two lines of research that have traditionally been treated independently but recently shown to be inherently intertwined18,44. In doing so, our work established the foundation for future studies to quantitatively compare them. For all analyses, brain regions were defined using the Schaefer 200 regions-17 networks parcellation solution56.

Local BOLD signal variability was calculated by taking the root mean squared successive difference (rMSSD) of each normalized regional time series. rMSSD quantifies moment-to-moment changes in the BOLD signal by measuring the mean of the squared differences in signal intensity between successive time points57. Greater local BOLD variability, therefore, is present in regions with greater rMSSD (Fig. 1A).

Global BOLD signal variability was estimated as dynamic functional connectivity using covSTATIS. As a three-way extension of multidimensional scaling and principal component analysis, covSTATIS is a multi-table method that linearly combines multiple similarity matrices (i.e., correlation/covariance tables) and uses eigenvalue decomposition to identify predominant structured patterns both at the group and individual level58–60. In our case, we used a combination of sliding window approaches and covSTATIS to examine, for each individual, the similarity of their windowed functional connectivity matrices over time and derive individual-level estimates of regional functional connectivity dynamics. Unlike most conventional fMRI dynamic connectivity methods61 that operate on coarser brain dynamics estimates, covSTATIS defines finer regional measures of global dynamics and does so with minimal user-dependent input, as it side-steps commonly adopted clustering approaches that have traditionally led to fractionated, study-specific definitions of dynamic functional connectivity62. For more information about covSTATIS and its implementation, we refer the reader to our dedicated recent publication60.

After partitioning each regional time series into equally sized windows via a sliding window approach63,64 (see “Methods” for details), we derived, for each window, functional connectivity measures for each region pair, as their product-to-moment correlation across timepoints (Fig. 1B step 1). This procedure resulted in NxNxT functional connectivity data tables for each individual, where N is the number of regions and T the number of windows. We next assessed, via the Rv similarity coefficient65, the similarity across all data tables across all individuals (Fig. 1B step 2). We then calculated their weighted average (a NxN data table), to obtain a group compromise space, where regional connections more similar across time and individuals were given a higher weight, since they were most represented in the sample (Fig. 1B step 3). We submitted the group compromise space to eigenvalue decomposition (Fig. 1B step 4) and obtained a multivariate connectivity space, wherein regions that showed similar connectivity values over time were closer together than regions with less similar connectivity values across windows (Fig. 1B step 5). covSTATIS next allowed us to back-project into this abstract multivariate Cartesian space, for every individual, each region’s mean connectivity value over time across all windows (Fig. 1B step 6, blue dot) and around it, the region’s connectivity value for each window (Fig. 1B step 6, green dots). Our last step involved calculating, for each individual, the area of the hull around each regional mean connectivity over time (Fig. 1B step 7). A greater area of the hull indicates greater distance/spread in connectivity across windows and is therefore characteristic of regions with greater dynamic functional connectivity, that is, greater global BOLD signal variability.

We next validated our method by assessing its consistency, test-retest reliability, and its ability to capture age effects typically reported in studies using traditional dynamic connectivity methods. We used resting-state fMRI data from one sample of 145 healthy young adults with two successive runs, and two independent cross-sectional healthy lifespan datasets (see “Methods” for details). For each sample and individual, we partitioned their regional time series data into windows of different length and derived covSTATIS’ area of the hull scores (i.e., global BOLD signal variability) for each iteration. We found global BOLD variability to be highly consistent across window lengths, both at the individual level (two-way mixed effects ICC model > 0.82 for all samples; Supplementary Fig. S1) and at the group level (0.92 < Spearman rho < 0.99, p < 0.001) in all samples. Global BOLD variability additionally showed high reliability at the group level (r = 0.98; p < .001; Fig. 1C) and moderate reliability at the individual level for all window lengths (average r ~ 0.30, p < .001, large Cohen’s d > 1.6). In line with current literature demonstrating an age-induced dampening of the brain’s dynamic range66, here we used partial least squares (PLS)67,68—a multivariate method that assesses the covariance between two or more sets of variables—and showed that covSTATIS-derived area of the hull scores also decreased with age, particularly in regions that preferentially show age effects in the literature69 (Fig. 1C; one significant latent variable (LV1) at p = 0.003 explaining 72% brain-age variance; lifespan sample 1 brain-age r = −0.19, lifespan sample 2 brain-age r = −0.30). Furthermore, age effects were stable across windows of different length, both in terms of their spatial location and directionality (LV1 at p < 0.001 explaining 58% brain-age variance; Supplementary Fig. S2). Together, these results highlight how covSTATIS is a valid and robust tool to estimate global BOLD signal variability.

Bridging across fMRI datasets: reliability and topography of local and global BOLD signal variability
We estimated local and global BOLD signal variability on two openly available resting-state fMRI datasets with diverse acquisition protocols. Given the heterogeneity of fMRI data used in the literature, it is imperative to quantify the dependency of local and global BOLD variability on the type of fMRI data used to extract them, to ensure generalizability. Here, we chose two datasets that differ in echo time and band acquisition. The number of echo times and bands influences the spatial and temporal resolution of fMRI, which in turn may impact the spatiotemporal properties of the underlying BOLD signal. Multi-echo acquisition allows for greater spatial coverage yet comes at the cost of slower acquisition time; multi-band acquisition allows instead for faster acquisition but is more susceptible to motion and scanner artifacts70–73. Throughout this paper, analyses were conducted on both samples in parallel. Young sample 1 refers to our multi-echo single-band fMRI dataset74, and young sample 2 refers to our single-echo multi-band fMRI dataset75 (see Fig. 2A and “Methods” for details about the samples).

Despite differences in the distribution of regional variability values across fMRI data type (Fig. 2A; rMSSD: young sample 1 mean(SD) = 23.73 (9.85), IQR = 12.75; young sample 2 mean(SD) = 13.93 (4.18), IQR = 5.16; covSTATIS: young sample 1 mean(SD) = 0.01 (0.002), IQR = 0.002; young sample 2 mean(SD) = 0.02 (0.003), IQR = 0.004; arbitrary units), local and global BOLD variability overall converged across samples (Fig. 2B, top). Greater regional and network reliability was found for global than local BOLD variability, as indicated by the greater number of regions and networks showing consistent mean values across fMRI samples (Fig. 2B, middle and bottom). For both local and global BOLD variability, reliability was lowest in sensorimotor (i.e., lower-order) areas and highest in heteromodal (i.e., higher-order), particularly default network, regions.

Building on these findings, we next described the spatial topography of local and global BOLD variability by fMRI data type. Understanding the spatiotemporal complexity of BOLD signal variability is necessary to build accurate computational and empirical models of brain function. Computational models typically treat macroscale signal variability as spatially homogenous to maximize mathematical tractability. Empirical fMRI studies mostly focus on the behavioral and clinical applicability of these measures. Furthermore, the known heterogeneity in spatial coverage across fMRI acquisition protocols is oftentimes overlooked.

For both fMRI data types, local and global BOLD variability showed a heterogeneous topography across brain regions (Figs. 3 and 4). Local BOLD variability presented greater topographical divergence across regions and networks than global BOLD variability, and its topography varied the most by fMRI data type (Fig. S3). Regional rank orders of local and global BOLD signal variability highlight the greater correspondence in higher-order cortices across fMRI data types (Fig. S4).

To contextualize these observed topographies, we mapped group-level regional profiles of local and global BOLD variability onto regional multi-scale maps of neocortical brain organization. Such maps were derived from open-source ex vivo cytoarchitectural50, in vivo microstructural51, ex vivo transcriptional (molecular)52,53 and static functional connectivity MRI data54. Both local and global BOLD variability were found to be organized along multiscale gradients: both measures significantly mapped onto more than one neurobiological system (Fig. 5). This multiscale mapping was consistent across fMRI data type. Local BOLD variability was situated along an anterior-posterior gradient that was maximally associated with laminar and cellular spatial organization (r = 0.52, p < 0.001 and r = 0.62, p < 0.001). Global BOLD variability instead evolved along a unimodal-transmodal gradient that preferentially related to underlying regional microstructure and static functional connectivity (r = −0.32, p = 0.03 and r = 0.47, p < 0.001). Remarkably, these two distinct axes of local and global BOLD variability provide direct replication of previous work showing how intrinsic fMRI dynamics, estimated as the optimum combination of more than 6000 temporal BOLD timeseries features, collectively evolve along these two axes of neocortical brain organization32.

Bridging across spatial scales: local and global BOLD signal variability sit at the intersection of multiscale neocortical organization
We next sought to comprehensively interrogate the multiscale properties of local and global BOLD signal variability. If variability is an emergent, crucial feature of complex systems, empirical fMRI measures of brain variability should recapitulate core aspects of multiscale brain organization. Building on this inference, we hypothesized that local and global BOLD variability would sit at the intersection of multiscale neocortical organization (Fig. 6, top left). To test this hypothesis, we investigated whole-brain relationships between our measures of local and global BOLD variability and open-source micro-, meso-, and macro-scale neurobiological variables, for each fMRI data type separately. Microscale neurobiological metrics included the first gradient of regional cytoarchitectural differentiation from ex vivo BigBrain histological data50, and the first gradient of microstructural differentiation from in vivo quantitative T1 imaging from the MICA-MICs dataset51. Mesoscale neurobiological metrics were calculated as composite scores on PET-derived whole-brain neuroreceptor density maps available through Neuromaps47. We derived composite scores, as opposed to individual neuroreceptor maps, given the heterogeneity of the molecular and chemical composition of each brain region76. Macroscale neurobiological metrics comprised: PET-derived whole-brain maps of oxygen metabolism, glucose metabolism, cerebral blood flow and cerebral blood volume; and large-scale gradients of brain organization, including fMRI-derived sensory-association axis, the principal gradient of fMRI static functional connectivity and MEG-derived intrinsic timescale, all downloaded from Neuromaps47. Based on recent work on the role of BOLD temporal autocorrelation properties in recapitulating macroscale brain organization77, for each of our fMRI samples, we additionally extracted group-level spatial maps of regional temporal autocorrelation scores from the BOLD signal. Such scores were derived as the product-to-moment correlation between successive (lag-1) and alternate (lag-2) timepoints of each regional time series tailored to each sample’s TR. Lastly, to expand on recent evidence on the link between arrhythmic oscillatory brain activity and brain stochasticity78, we leveraged open-source MEG data79,80 and parametrized neurophysiological spectra81, to spatially characterize whole-brain cortical arrhythmic brain activity and obtain a group-level spatial map. Each one of these group-level multiscale maps was separately correlated with group-level local and global BOLD variability maps from each of our fMRI samples, across regions. Statistical significance was assessed via 10,000 Hungarian spins on the regional parcellation of our local and global BOLD variability maps.

Overall, across fMRI data types, local and global BOLD variability were both associated with several neurobiological measures belonging to each spatial scale (Fig. 6, middle and right). At the microscale, greater local BOLD variability was associated with greater laminar differentiation (r = 0.52, p = 0.007 young sample 1). Specifically, greater rMSSD was present in regions showing heightened differentiation in cell size and density and clearer cortical layer separation. Motivated by previous work on the cytoarchitectural properties of static BOLD signal measures33, we next derived cortical thickness for supragranular (mean across layers I–III), granular (layer IV) and infragranular layers (mean across layers V–VI) from BigBrain data50, and related their whole-brain spatial distribution to our local BOLD variability spatial maps. We predicted local BOLD variability to be associated with the expression of granular layer IV in particular. Layer IV receives feedforward thalamo-cortical inputs82,83. If regions with a prominent layer IV need to orchestrate incoming feedforward projections, then such regions may exhibit greater variability in their functional activity. Theories about the thalamo-cortical pathway have proposed high local BOLD variability along these connections84. Additionally, layer IV is thickest in visual areas and absent in motor regions, precisely where we observed the highest and lowest levels of local BOLD variability. In line with our predictions, we found that greater local BOLD variability was significantly associated with stronger layer IV expression (r = 0.39, p = 0.01 young sample 1; Fig. S5). At the mesoscale, greater local BOLD variability was related to greater ionotropic receptor density (r = 0.39, p = 0.002 young sample 1), reduced metabotropic receptor density (r = −0.23, p = 0.008 young sample 2), decreased receptor diversity (r = −0.25, p = 0.002 young sample 1), and decreased excitation/inhibition (E/I) ratio (r = −0.47, p < 0.001 young sample 1; r = −0.31, p = 0.04 young sample 2). At the macroscale, heightened local BOLD variability was associated with increased oxygen (r = 0.41, p = 0.01 young sample 1; r = 0.54, p < 0.001 young sample 2) and glucose metabolism (r = 0.62, p < 0.001 young sample 2), along with greater cerebral blood flow (r = 0.50, p < 0.001 young sample 2) and volume (r = 0.41, p = 0.009 young sample 1).

At the microscale, greater global BOLD variability covaried with thinner infragranular layers (r = −0.45, p = 0.004 young sample 1; Fig. S5). Greater global BOLD variability was also related to overall decreased microstructural differentiation (r = −0.32, p = 0.03 young sample 2). At the mesoscale, greater global BOLD variability was associated with increased metabotropic receptor density (r = 0.33, p < 0.001 young sample 1; r = 0.26, p = 0.003 young sample 2). At the macroscale, it was positively related to the principal gradient of fMRI static functional connectivity (r = 0.47, p < 0.001 young sample 2) and the sensory-association axis (r = 0.35, p = 0.01 young sample 2).

To quantify the robustness of these multiscale correlations across fMRI data type, we next quantified the degree of overlap in the neurobiological correlates of local and global BOLD variability across our two fMRI samples. We computed rank correlations on the Fisher-z transformed correlation vectors characterizing the relationships between local and global BOLD variability and multiscale variables for each fMRI sample. Both local and global BOLD variability exhibited strong overlap in their associations with multiscale neurobiological factors across fMRI data types (Fig. 6, bottom).

To quantify the central role of local and global BOLD variability in multiscale brain organization, we next ran cartographic analyses on our two sample-specific correlation matrices. Cartographic analyses are commonly used to derive graph theory metrics from brain networks35. We first assigned local and global BOLD variability, micro-, meso- and macro-scale measures to four different communities. We then took the absolute value of the reported correlations and calculated, for each fMRI sample, the participation coefficient of local and global BOLD variability39,85. Participation coefficient scores allowed us to determine how evenly distributed across spatial scales were the correlations of local and global BOLD variability, for each fMRI data type. Scores closer to 1 indicate greater multiscale participation86. We found high participation scores for both metrics and fMRI samples (young sample 1: local = 0.65, global = 0.71; young sample 2: local = 0.49, global = 0.59).

To further validate the multiscale nature of local and global BOLD variability, we used dominance analysis (DA) (see “Methods” for details) to build a predictive model for each fMRI dataset, where we estimated the unique contribution of each neurobiological measure in predicting local and global BOLD variability18,87. In line with our correlational results, we found that both local and global variability were predicted by a combination of neurobiological variables within each sample (Fig. 7).

As a final step, we repeated all analyses above, including the second gradient of functional connectivity organization54. Broad patterns of association were largely maintained across the two gradients, yet the principal gradient showed more reliable inter-sample results (Fig. S6).

Bridging across temporal scales: local BOLD signal variability is anchored in electrophysiological processes
To further understand the temporal and neuronal properties of BOLD signal variability, we turned to electrophysiological data and used a combination of open-source MEG source-modeled, broadband (1–150 Hz) resting-state data79,80 and simulations of naturalistic electrophysiological time series, to derive and characterize local measures of brain variability. We decided to compute only measures of local variability and not of global variability, in light of the known methodological variance of functional connectivity derivatives in the MEG literature88.

We derived two measures of local brain signal variability on MEG regional time series, based on how variability is independently characterized in fMRI and MEG: (1) moment-to-moment signal intensity changes via rMSSD (MEG signal variability), and (2) 1/f exponent, that is, the slope of the MEG power spectrum, shown to capture background arrhythmic activity89,90. We first looked at the spatial topography of both measures. For the former, we observed the highest variability levels in somato-motor areas, lowest values in visual and dorsal attention regions, and mid-levels in higher-order cortices (Fig. 8A; network maps in Fig. S3). For the latter, we replicated previous reports showing steeper slopes (i.e., greater 1/f exponent values) in posterior cortical regions89 (Fig. 8B, left).

Next, for each individual, we obtained whole-brain measures of local MEG variability for both rMSSD and 1/f exponent measures by averaging across brain regions. We then related the two variables across individuals and found a strong negative association between them (r = −0.6; p < .001, Fig. 8B middle): individuals with a flatter 1/f exponent showed heightened levels of local rMSSD-derived MEG variability. To expand on these findings, we built a hierarchical linear model where we tested for rMSSD-1/f exponent relationships while accounting for regional heterogeneities in both measures (Fig. S7). The model showed regional diversity in the association between rMSSD and the 1/f exponent (Fig. S7). To further get a mechanistic understanding of this relationship, we simulated naturalistic electrophysiological time series using the NeuroDSP toolbox91 (see “Methods” for details). We manipulated the steepness of the 1/f exponent at various parameters ranging from −0.7 to −1.5 in steps of −0.1. 10 simulations were run per step, and local rMSSD-derived MEG variability was calculated for each manipulation. We found that flatter 1/f exponents were strongly predictive of greater local MEG variability (R2 = 0.97; β(SE) = −9.86 (0.02); CI [−9.51; −9.10]; p < 0.001; Fig. 8B right panel). These results highlight the arrhythmic nature of electrophysiological signal variability.

As a final step, we tested cross-modal relationships between local fMRI and MEG signal variability. We report both cross-sample and within-sample associations (i.e., fMRI-MEG data obtained on the same individuals from the MEG sample with available fMRI data). At the group level, we found a consistent, negative association across and within sample between fMRI and MEG variability (Fig. 8C). A similar negative association was also observed at the individual level (average Fisher z-to-r = −0.14; t(102) = −11.91, p < 0.001, 95% C.I. = [−0.17, −0.12], Cohen’s d = −1.17). Together, these findings reinforce the biological multimodal nature of local brain signal variability.

Methods
Neuroimaging datasets
All main analyses were carried out on two open-source resting-state fMRI datasets: young samples 1 and 2. To validate covSTATIS as a global BOLD signal variability estimation method, we leveraged two adult lifespan datasets: lifespan samples 1 and 2. To evaluate the temporal and neuronal properties of BOLD variability and to bridge across neuroimaging modalities, we used an open-source resting-state MEG dataset. Below, we briefly describe each neuroimaging dataset.

fMRI young sample 1
A total of 150 healthy young individuals ages 18–34 (M age = 22 years, SD age = 3 years, 55% F) from the Neurocognitive Aging Dataset74 were included in this study. Gender information was not available, and all analyses were carried out on the whole sample (both female and male data together). All individuals in the study were right-handed, healthy, with no evidence of neurological, psychiatric, or other medical conditions impacting neural functioning69. Participants underwent two multi-echo resting-state fMRI scans of 10-min duration each within the same session. Analyses conducted in this paper were performed on the first run of data. The second run was used to assess test-retest reliability of global BOLD signal variability on 145 individuals, since 5 participants out of the total sample did not have a second run of data (M age = 22 years, SD age = 3 years, 55% F). Data were collected at the Cornell Magnetic Resonance Imaging Facility, at Cornell University (New York, US). All participants provided written informed consent and were compensated for their time. Research protocols were approved by the Cornell University Institutional Review Board.

Resting-state fMRI data were acquired on a 3T GE discovery MR750 using a multi-echo EPI sequence with online reconstruction (TR = 3000 ms; TE 1 = 13.7 ms, TE 2 = 30 ms, TE 3 = 47 ms; 83° flip angle; matrix size = 72 × 72; FOV = 210 mm; 46 axial slices; 3 mm isotropic voxels; 204 volumes) with 2.5 × acceleration and sensitivity encoding. Participants were instructed to keep their eyes open during the two scans. Multi-echo resting-state fMRI data were preprocessed and denoised using multi-echo independent component analysis (ME-ICA v3.2; https://github.com/ME-ICA/me-ica)70,71, as in our previous publication69. Briefly, functional images were minimally preprocessed (i.e., removal of the first 4 volumes, deobliquing, motion correction, anatomical-functional co-registration, spatial alignment of volumes across echoes-TEs). T2* maps were then used for anatomical-functional co-registration. Volumes were optimally combined across TEs and denoised. Image quality was lastly carried out on the denoised time series to check for co-registration and residual noise (framewise displacement > 0.50 mm, DVARS > 1, tSNR < 50, or fewer than 10 total included BOLD components). Functional images were parcellated using the group prior individual parcellation (GPIP)69,149, a participant-specific parcellation approach initialized on the Schaefer 200-17 network solution56. We selected GPIP as our parcellation method of choice to ensure high topographical precision of fMRI BOLD, and for consistency with previous work from our group using the same data69. Unlike standard parcellations, GPIP accounts for within-subject variation in parcel boundaries, offering substantially improved sensitivity to individual differences in brain activity and behavior149–152. Specifically, GPIP consists of an iterative Bayesian process: it starts from the user-provided reference group atlas (Schafer 200-17 in our case) and iteratively refines each individual’s parcel boundaries relative to their resting-state connectivity data149. It then computes concentration matrices (inverse covariance/partial correlation) for the whole sample via a group sparsity constraint. To ensure stability and precision, the algorithm iterates between the former and latter step for a total of 20 iterations or until no more than one vertex is changing per parcel149. For more details about GPIP, the reader is referred to our previous publication where we used GPIP on this sample69. Regional high kappa time series from ME-ICA were used to derive measures of local and global BOLD signal variability.

fMRI young sample 2
A total of 112 healthy young individuals age matched to young sample 1, ages 18–34 (M age = 23 years, SD age = 3 years, 54% F), from the Enhanced Nathan Kline Institute Rockland Sample75 were included in this study. All analyses were carried out on the whole sample (combined sex data). All individuals in the study were right-handed, healthy, with no current or past neurological, psychiatric diagnosis153. All participants provided written informed consent and were compensated for their time. Research protocols were approved by the NKI Institutional Review Board.

We included one run of resting-state fMRI data acquired on a 3T Siemens Trio scanner using a multiband (factor of 4) EPI sequence (TR = 1400 ms; TE = 30 ms; 65° flip angle; FOV = 224 mm; 64 axial slices; 2 mm isotropic voxels; 404 volumes). Participants were instructed to keep their eyes open during the scan75. Functional images were preprocessed in the same fashion as our previous publication153. Briefly, fMRI data were preprocessed using the data preprocessing assistant for resting-state fMRI advanced edition toolbox (DPARSF-A v4.3_170105; http://rfmri.org/DPARSF)153 as follows: removal of first 5 volumes, despiking, realignment, normalization to MNI template, smoothing (6 mm FWHM), ICA denoising via ICA-FIX trained on n = 24 participants, and nuisance regression (linear detrend, Friston 24 motion parameters, bandpass filter of 0.01–0.1 Hz). Data were lastly parcellated using the standard Schaefer 200-17 network solution56. Regional time series were used to calculate local and global BOLD signal variability.

fMRI adult lifespan sample 1
A total of 154 healthy adults ages 20–86 (M age = 49 years, SD age = 19 years, 62% F) from the greater Toronto area were included in this study. All analyses were carried out on the whole sample (both female and male data together). All individuals in the study were screened to be healthy with no major neurological and/or psychiatric condition, they were right-handed, fluent English speakers, cognitively normal (mini-mental state examination > 26), and with normal or corrected-to-normal vision19,20. All participants provided written informed consent and were compensated for their time. Research protocols were approved by the Research Ethics Board at Baycrest Health Sciences Center.

Resting-state fMRI data were collected on a 3T Siemens Trio scanner using an EPI sequence (TR = 2000 ms; TE = 27 ms; 70° flip angle; FOV = 192 mm; 40 axial slices; 3 mm isotropic voxels; 297 volumes). Participants were instructed to keep their eyes open during the scan. Functional images were preprocessed in the same fashion as our previous publication, using the optimizing of preprocessing pipelines for neuroimaging software package (OPPNI) (available at https://github.com/strotherlab/oppni)19. Briefly, brain masks were created for each participant in FSL, functional volumes were co-registered, corrected for motion, detrended and temporally filtered within 0.01–0.1 Hz. ICA-FIX was then used to denoise the data, trained on n = 40 participants, and further nuisance regression was applied on the denoised time series (motion, white matter, vessels, CSF). Lastly, data were smoothed (7 mm FWHM), normalized to MNI space, and resampled to 4 mm3 isotropic voxels. Data were parcellated using the standard Schaefer 200-17 network solution56. Regional time series were used to calculate local and global BOLD signal variability.

fMRI adult lifespan sample 2
A total of 154 healthy adults, age and sex-matched to our lifespan sample 1, ages (M age = 49 years, SD age = 19 years, 62% F), from the Enhanced Nathan Kline Institute Rockland Sample75 were included in this study. All analyses were carried out on the whole sample (combined sex data). All participants provided written informed consent and were compensated for their time. Research protocols were approved by the NKI institutional review board. Refer to fMRI young sample 2 for details about the resting-state data.

fMRI-MEG sample
Within-individual fMRI-MEG data were obtained from the openly available stage II of the Cambridge Centre for Aging and Neuroscience (Cam-CAN) data repository79,80. We included resting-state MEG data from a total of 104 healthy young adults in the same age range 18–34, as our primary fMRI datasets (M age = 28 years, SD age = 4 years, 56% F). All individuals in the study were screened to be healthy with no major neurological and/or psychiatric condition, they were right-handed, and fluent English speakers79,80. All MEG analyses, except for fMRI-MEG local signal variability correlations, were run on this sample. All analyses were carried out on the whole sample (combined sex data). fMRI-MEG associations were conducted on 103 out of the 104 individuals, that is on those subjects with available resting-state MEG and movie-watching fMRI data. All individual provided written informed consent and were compensated for their time. Research was conducted in compliance with the Helsinki Declaration and was approved by the Cambridgeshire 2 Research Ethics Committee.

Resting-state MEG data acquisition and processing. All participants underwent an approximately 8-min resting-state eye-closed MEG and a structural T1 MRI. MEG data were collected from a 306-channel VectorView MEG system (Elekta Neuromag, Helsinki) with 102 magnetometers and 204 orthogonal planar gradiometers at 1000 Hz sampling rate. The head position of participants was continuously monitored using four head-position indicator coils, while ocular (EOG) and cardiac (ECG) external electrodes were used to monitor physiological artifacts. See 79,80 for details on the dataset and data acquisition.

MEG data were preprocessed using Brainstorm81. The preprocessing pipeline followed previous published work154. Line noise artifact (50 Hz; with 10 harmonics) were removed using a bank of notch filters, in addition to 88 Hz noise—a common artifact characteristic to the Cam-CAN dataset. Slow-wave and DC-offset artifacts were removed with a high-pass FIR filter with a 0.3-Hz cut-off. Signal-space projectors were used to remove cardiac artifacts and attenuate low-frequency (1–7 Hz) and high-frequency noisy components (40–400 Hz) due to saccades and muscle activity. The MRI volumes of each participant were automatically segmented using Freesurfer155 and co-registered to the MEG recording using approximately 100 digitized head points. We constrained brain source models for each participant to their individual T1-weighted MRI data. We computed head models for each participant using the Brainstorm overlapping-spheres approach, and cortical source models using the Brainstorm implementation of linearly constrained minimum-variance beamforming (2018 version for source estimation processes); both processes were run using default parameters. MEG source orientations at 15,000 vertices were constrained normal to the cortical surface. We projected the resulting MEG source model for each participant onto the default anatomy of Brainstorm (ICBM152). We down-sampled the cortical source maps by averaging the time series within each parcel of the standard Schaefer 200-17 atlas56. This resulted in a matrix of 200 regions by 150s of scouted MEG time series data per person. Regional time series data were used to calculate local MEG variability.

We computed power spectral density (PSD) estimates across all vertices of the source model using the Welch method with 2s windows of 50% overlap. This resulted in PSD with a frequency resolution of 1/2 Hz. We down-sampled the PSD estimates to the Schaefer atlas by taking the mean spectral density within each ROI. We parametrized the resulting neural power spectra from 1 to 40 Hz using the specparam algorithm as implemented in Brainstorm81 with the following parameters: peak width limits [0.5 12]; maximum number of peaks: 3; minimum peak amplitude: 0.3 a.u.; peak threshold: 2.0 SDs; proximity threshold: 2.0 SDs; aperiodic mode: fixed. Finally, we extracted the aperiodic component from the resulting spectral models at each parcel, for each participant (i.e., 1/f exponent)89.

Movie-watching fMRI data acquisition and processing. Participants were shown a shortened 8-min version of the “Bang! You’re Dead” episode from the TV show “Alfred Hitchcock Presents” (1961). None of them had seen the movie before. Participants were told to pay attention to the movie while lying still in the scanner. Data were acquired on a Siemens 3T TIM Trio scanner using a 32-channel head coil. Functional data were collected with a multi-echo T2* EPI sequence with the following details: 193 total volumes, 32 axial slices, 3.7 mm thick, 0.74 mm gap, TR = 2470 msec, TE = [9.4, 21.2, 33, 45, 57] msec, flip angle = 78°, FOV = 192 × 192 mm, voxel size = 3 × 3 × 4.44 mm. Structural T1-weighted images were also used in this study for co-registration purposes. The latter images were acquired with a 3D MPRAGE sequence with the following details: TR = 2250 msec, TE = 2.99 msec, TI = 900 msec, flip angle = 9°, FOV = 256 × 240 × 192 m, 1 mm isotropic voxels, GRAPPA acceleration factor = 2.

Functional data were preprocessed using fMRIPrep v21.0.1156 (minimal preprocessing) and tedana v2.5157 for denoising (optimal combination of echoes, PCA using the “stabilized” Kundu decision tree, ICA to decompose BOLD TE-dependent and non-BOLD TE-independent signals)158.

Neurobiological data
To contextualize local and global BOLD signal variability with multiscale descriptions of brain organization, we relied on multiple open-source repositories and toolboxes.

Microscale data
We included the principal axes of cytoarchitectural and microstructural differentiation obtained via a nonlinear manifold learning technique called diffusion map embedding159,160, on histological staining of a ex vivo human brain and in vivo quantitative T1 (qT1) data.

The first gradient of cytoarchitectural differentiation was obtained from ultra high-resolution histological information from the BigBrain dataset50. BigBrain is a three-dimensional model of an adult human brain (Caucasian male, age 65), reconstructed from 20-mm-thick slices of a coronally-sectioned, Merker-stained post-mortem specimen. The profile of staining intensity along the thickness of the cortical sheet is thought to capture the composition of local cellular assemblies. To assure coverage along the thickness of the cortical sheet, cellular staining intensity was sampled and averaged across equivolumetric intracortical surfaces at every voxel step for (100-mm voxels) across 163,842 vertices per hemisphere. The resulting staining intensity surface was then parcellated using an anatomical atlas-constrained clustering approach, which produced a 1012-nodes solution for the BigBrain specimen, constrained by the Desikan–Killiany and Destrieux atlas boundaries. Pairwise product-to-moment correlation of nodal staining intensity profiles resulted in the affinity matrix used for gradient embedding.

The first gradient of microstructural differentiation was obtained from qT1 intensities sampled from and averaged across 50 healthy young participants (M age = 29.54 years, SD age = 5.62 years, 46% F) in an openly available dataset (MICA-MICs)51. All individuals were right-handed and without a history of neurological and/or psychiatric conditions. All analyses were carried out on the whole sample (combined sex data). All participants provided written informed consent and were compensated for their time. Research protocols were approved by the Ethics Committee of the Montreal Neurological Institute and Hospital. The profile of qT1 intensity along the thickness of the cortical sheet is shown to capture the variation in myelination in cortical regions. To assure coverage along the thickness of the cortical sheet on qT1 volumes, 14 equivolumetric intracortical surfaces were generated for each individual and combined into vertex-wise averages. Intensity profile maps were then parcellated according to the Desikan–Killiany and Destrieux atlas and averaged across 50 healthy young adults. Cross-correlation of nodal qT1 intensity profiles computed with partial correlation, with whole-cortex intensity profile controlled for, followed by log-transformation, resulted in the affinity matrix used for gradient embedding51.

These two principal gradients of cytoarchitectural and microstructural differentiation were released with the BigBrainWarp toolbox in BigBrain native volume161. We resampled both gradients from BigBrain native volume to the standard volume defined by the 1 mm isotropic ICBM 2009c Nonlinear Asymmetric brain template162 using linear interpolation. We finally parcellated both gradients using the standard Schaefer 200-17 parcellation solution56.

Equally from BigBrain Warp, layer thickness data from the BigBrain specimen for all six cortical layers were additionally included in this study. Briefly, staining intensity profiles were obtained from curvature-adjusted equidistant sampling at 200 points along the thickness of the cortical sheet. Layer transitions were identified by a convolutional neural network guided by expert neuroanatomists163. We processed and parcellated thickness maps for each layer as above. Thickness values were averaged across layers I–III to obtain supragranular estimations and across layers V–VI for infragranular estimations. Layer IV thickness was referred to as granular in the text.

Mesoscale data
We used openly available neurotransmitter receptor data from the Neuromaps toolbox47. Specifically, we included PET-derived receptor density distributions for the following neurotransmitter systems: serotonin (5-HT1a, 5-HT1b, 5-HT2a, 5-HT4, and 5-HT6)164, dopamine (D1165, D2166), GABA (GABAa, GABAbz)167, glutamate (mGluR5)168, and acetylcholine (a4b2169, M1170). In addition to the receptor data from Neuromaps, we also include a map of NMDA receptor density, fetched from https://github.com/netneurolab/hansen_receptors171–173. All maps were parcellated with the standard Schaefer 200-17 parcellation solution56 using Neuromaps.

Given the heterogeneity in the molecular and chemical composition of each brain region76, the following composite scores were computed: receptor diversity, excitation/inhibition ratio, ionotropic density and metabotropic density. Receptor diversity was estimated via the normalized Shannon entropy H derived as shown in Eq. 1. d indicates a region’s normalized receptor composition, a vector wherein each value reflects each receptor’s density value normalized to its maximum value across brain regions. L is the total number of receptors, 13. Greater H values signify greater receptor diversity for a brain area. Regional excitation/inhibition ratio was computed as the ratio between the mean density of excitatory receptors and the mean density of inhibitory receptors for each region (see Table 1 for a breakdown of receptor effects and types). Regional ionotropic receptor density was calculated as the mean density of ionotropic receptors within each area. Similarly, regional metabotropic receptor density was obtained as the mean density of metabotropic receptors within each region.

H = −sum(d * log(d)) / log(L) (1)

We computed the first principal component of gene expression from the Allen human brain Atlas52, what we referred in the study as transcriptional/molecular gradient. Regional microarray expression data were obtained from 6 post-mortem brains (1 female, ages 24–57 years, M age = 42.5 years, SD age = 13.38 years) and were processed using the Abagen toolbox53,174,175. To ensure complete coverage across both hemispheres, we mirrored samples bilaterally and interpolated missing voxels using a nearest-neighbor approach. All other parameters in Abagen were set to default. Finally, we only retained genes with differential stability >= 0.1, as a means of filtering out genes with high variability across donors176. Altogether, 15,633 genes were used in the principal component analysis. Note that we followed the same processing pipeline as in ref. 177, except we used the Schaefer 200-17 parcellation solution56.

Macroscale data
We retrieved macroscale brain maps available through the Neuromaps toolbox47 and parcellated them with the standard Schaefer 200-17 network atlas56. These maps include: PET-derived maps of oxygen metabolism, glucose metabolism, cerebral blood flow, and cerebral blood volume178; large-scale gradients of brain organization, specifically the sensory-association axis55, the principal gradient of fMRI static functional connectivity54, and the MEG-derived intrinsic timescale179. Lastly, we calculated temporal autocorrelation scores on our two primary fMRI datasets by taking the product-to-moment correlation between successive timepoints (lag-1) and alternate (lag-2) timepoints of each regional BOLD time series tailored to each sample’s TR, that is, lag-1 autocorrelation for young sample 1 and lag-2 for young sample 2, and created group-level spatial maps.

Local signal variability estimation
Local fMRI BOLD signal variability. Local BOLD signal variability was estimated for every individual within each main fMRI sample (young samples 1 and 2). Denoised regional BOLD time series were first mean-centered to a whole-brain mean of 0. Moment-to-moment temporal variability of the BOLD signal was then calculated by taking the rMSSD57 of each normalized regional time series, as shown in Eq. 2, where x represents a region’s BOLD signal intensity at two successive timepoints i and i+1, and n is the total number of timepoints for that region. At each region, we took the square root of MSSD to preserve its original units (MSSD is equivalent to variance, rMSSD to SD). Greater rMSSD values are indicative of regions with greater local variability levels.

rMSSD = sqrt( (1 / (n-1)) * sum((x_{i+1} - x_i)^2) ) (2)

Local MEG signal variability. The same formula as above was used to derive the regional moment-to-moment temporal variability of MEG time series via rMSSD.

Global fMRI BOLD signal variability estimation and validation
Estimation. Global BOLD signal variability was estimated via dynamic functional connectivity for every individual within each fMRI sample (young samples 1 and 2, lifespan samples 1 and 2). First, each regional BOLD time series was partitioned in equally sized windows via a sliding window approach within each sample. Following recent guidelines180 to ensure reliable estimations, window width was chosen to comprise between 20 and 40 timepoints (TRs), window length was between 45 and 60 s, window shape was set to squared, Leonardi high pass filtering was applied181, and windows were shifted by 1 TR. Given the heterogeneity in sampling rate across fMRI samples, this procedure resulted in a different number of timepoints per window and in a different number of total windows across datasets (young sample 1: 60s window, 20 TRs per window, 181 total windows; young sample 2: 45s window, 32 TRs per window, 368 total windows; lifespan sample 1: 45s window, 23 TRs per window, 275 windows; lifespan sample 2: 45s window, 32 TRs per window, 368 total windows). For each window and dataset, we calculated functional connectivity measures per region pair, as their product-to-moment correlation. NxNxT functional connectivity data tables for each individual were thus derived, where N is the number of regions and T the number of windows. Unlike most dynamic connectivity approaches that apply clustering methods on the data tables, in this study, we wanted to minimize additional user input and maximize data fidelity. We therefore introduced covSTATIS to derive dynamic functional connectivity, an approach that does not rely on clustering analyses.

As a three-way extension of multidimensional scaling and principal component analysis, covSTATIS is a multi-table method that linearly combines multiple similarity matrices (i.e., correlation/covariance tables) and uses eigenvalue decomposition to identify predominant structured patterns both at the group and individual level58–60. In our case, we used a combination of sliding window approaches and covSTATIS to examine, for each individual, the similarity of their windowed functional connectivity matrices over time and derive individual-level estimates of regional functional connectivity dynamics. First, we assessed, via the Rv similarity coefficient (i.e., squared product-to-moment correlation)65, the similarity across all Y data tables across all individuals within each sample and stored these values in a YxY Rv matrix. We then took the first eigenvector of the Rv matrix to generate weights for each connectivity table and used these weights to build an NxN group compromise space—the linear combination of all matrices where regional connections more similar across time and individuals were given a higher weight since they were most represented in the sample.

We then submitted the group compromise space to eigenvalue decomposition and obtained a multivariate connectivity space, wherein regions that showed similar connectivity values over time were closer together than regions with less similar connectivity values across windows. covSTATIS next allowed us to back-project into this abstract multivariate Cartesian space, for every individual, each region’s mean connectivity value over time across all windows and around it, the region’s connectivity value for each window. Our last step involved calculating, for each individual, the area of the hull around each region’s mean connectivity over time. Each convex hull was peeled on 95% of the data to control for outliers, similarly to traditional dynamic functional connectivity approaches64. Global BOLD variability thus corresponds to the regional area of the hull values: a greater area of the hull indicates greater distance, hence spread, in connectivity across windows for a specific region, and is therefore characteristic of regions with greater global BOLD variability.

Importantly, covSTATIS is not unlike other dynamic functional connectivity approaches, where window size selection depends on TR and scan length. We therefore advise the reader to follow the guidelines already put forward by others in this space180, as they apply to all methods that rely on sliding window approaches. To reiterate, the advantage of covSTATIS over other dynamic functional connectivity methods lies in the minimal number of user-dependent decisions post sliding window estimation (e.g., covSTATIS does not require any clustering). We refer the reader to our covSTATIS-dedicated publication where we contextualize our method, comprehensively walk through all its details, and offer a tutorial to facilitate its implementation using Open Data60.

Validation. To validate covSTATIS as a robust dynamic functional connectivity method, we partitioned regional time series data from our young sample 1 (for both run 1 and 2) and our two lifespan samples into windows of 45, 60, 80 and 100 s, and used covSTATIS to derive global BOLD signal variability scores for each iteration and sample. Table 2 provides a detailed breakdown of our windowing procedure per sample and window size. After obtaining global BOLD signal variability scores, we used these measures to assess (1) covSTATIS’ consistency across window lengths (all samples); (2) its test-retest reliability across runs for each window length (young sample 1, runs 1–2); and (3) its ability to robustly capture age-related alterations in dynamic functional connectivity across window lengths (lifespan samples 1–2). All our validation analyses were carried out both at the individual and at the group level.

For each sample, covSTATIS’ consistency was evaluated at the individual level via intraclass correlation coefficients (ICC) and their 95% confidence intervals based on a mean-rating (k=4 window lengths), consistency, two-way mixed effects model with window size as fixed effect and individuals as random term182. At the group level, consistency was assessed by relating in a pairwise fashion window-specific covSTATIS results across brain regions via Spearman rho statistics. At the individual level, covSTATIS’ test-retest reliability was determined for each window length via Fisher-z transformed Spearman rho correlations derived within subject, then averaged across individuals and finally converted back to sample-average rho estimates. We additionally performed t-tests contrasting Fisher-z transformed correlations against zero and reported t-test values, along with their confidence intervals and Cohen’s d effect sizes. At the group level, covSTATIS’ test-retest reliability was instead assessed by correlating group-level covSTATIS scores across regions between runs.

As a final validation step, we implemented a multivariate analysis technique, PLS67,68,183, on our two lifespan samples, to test whether covSTATIS area of the hull measures were sensitive to known age-related alterations in dynamic functional connectivity across the adult lifespan. In other words, we interrogated the covariance between covSTATIS measures and age, and we validated our results across windows of different lengths. Specifically, we ran two PLS analyses: (1) contrasting covSTATIS’ area of the hull scores calculated on 45s windows for both lifespan samples (i.e., a 2-group 1-condition design, these are the main results reported in Fig. 1C); (2) contrasting covSTATIS values for all four window lengths across the two lifespan samples (i.e., a 2-group 4-condition design, these results are shown in Supplementary Fig. S2).

Briefly, PLS calculates a covariance matrix between two (or more) sets of variables. This covariance matrix undergoes singular value decomposition and, as a result, orthogonal latent variables are generated (LVs; similar to principal components in PCA), which explain the covariance between the sets of measures. Resulting LVs consist of a left singular vector (U) of age weights, a right singular vector (V) of brain regions, and a diagonal matrix of singular values (S). Each element of V, also called loading/salience, represent the contribution of each region to each LV. To identify significant LVs, singular values were permuted 1000 times. For each individual, we obtained an estimate of the degree to which they expressed a particular LV’s spatial pattern (brain score) by multiplying each regional loading by their original value and summing over all brain regions. We then correlated subject-level brain scores with age to establish the age contribution to the observed spatial pattern. Significant correlations were identified by bootstrapping with 1000 resamples, the correlation values, and generating 95% confidence intervals around the original correlation values. Lastly, to determine the significance of brain regions to each LV, we applied 1000 bootstrap resamples on the regional loadings. A t-like statistic was derived, the bootstrap ratio (BSR), which is the ratio of each regional weight to its bootstrapped standard error. We applied a threshold to BSRs at a value of ±2.

Inter-sample reliability of local and global fMRI BOLD variability
To evaluate the reliability of local and global BOLD variability, we contrasted local and global BOLD signal variability between our two main fMRI samples, young samples 1 and 2. We calculated reliability across brain regions at the group level, and within each brain region and network across samples. Group reliability across regions was derived via product-to-moment correlation between samples for each metric. Regional- and network-level reliability were instead estimated via independent t-tests, one test per region and one per network. This allowed us to assess mean differences in local and global BOLD variability levels across samples for each region and network. Reliability was determined at p > 0.05.

Local and global fMRI BOLD variability topographies
To examine the spatial organization of local and global BOLD variability for each of our main fMRI samples, we calculated group-level local and global BOLD variability for the whole brain and for each functional network separately. Despite using the Schaefer 200-17 parcellation solution56, we decided to reduce the number of networks from 17 to 7 to facilitate interpretation, by merging together regions from different subnetworks into their principal network (e.g., visual central and peripheral into visual).

As a final step, we correlated group-level local and global BOLD variability spatial maps for each main fMRI sample, with 4 different spatial maps obtained on neurobiological data: the first gradient of cytoarchitectural differentiation, microstructural differentiation, transcriptional/molecular differentiation and static functional connectivity differentiation (details in previous sections). To test for significance, we applied 10,000 Hungarian spins on the regional labels of local and global BOLD variability spatial maps184. This allowed us to obtain 10,000 null models, while preserving spatial autocorrelation, to compare our original correlations against.

Multiscale analyses
We assessed the multiscale neurobiological correlates of local and global BOLD variability, separately for our two primary fMRI datasets, by relating group-level local and global BOLD variability maps to micro-, meso-, and macro-scale variables described in previous sections. We computed the product-to-moment correlation, across brain regions, between local and global BOLD variability and each neurobiological variable. To evaluate the significance of the correlations, we applied 10,000 Hungarian spins on the regional labels of local and global BOLD variability. Results are visualized in text as annoted heatmaps. For an alternative visualization, we entered all correlation values, separately for each fMRI sample, into the Fruchterman Reingold layout spring embedding algorithm and set an absolute threshold of 0.3.

To quantify inter-sample convergence, we ran Pearson’s product-to-moment correlations on the Fisher-z transformed correlations vectors characterizing the relationships between local and global BOLD variability, and multiscale variables for each fMRI sample.

To test for the central role of local and global BOLD variability in brain organization, we computed cartographical analyses, similarly to 39,85, on our sample-specific multiscale correlation matrices. We first assigned (1) local and global BOLD variability, (2) microscale variables, (3) mesoscale and (4) macroscale measures to four different communities. We then took the absolute value of the reported correlations and calculated, for each fMRI sample, the participation coefficient of local and global BOLD variability. Briefly, participation coefficient scores allowed us to determine how evenly distributed across spatial scales the correlations of local and global BOLD variability were, for each fMRI data type: the closer the score to 1, the greater their multiscale participation86.

Finally, to assess the unique contribution of each multiscale variable in predicting local and global BOLD variability, we ran a DA per each metric and sample. DA was run, predicting local and global BOLD variability from all multiscale variables at once. DA allowed us to estimate the relative importance of each predictor in a single multiple regression model. DA contrasts two predictors at a time against all possible sub-models (2p-1 sub-models, p being the total number of predictors) and quantifies the incremental contribution of each predictor when added to each subset of the remaining predictors, as the increase in R2. In our study, we derived percentage scores for each predictor, indicating their unique contribution to the prediction model87.

Multimodal analyses
Similarly to how we characterized the topography of local and global fMRI BOLD variability, we delineated the spatial organization of group-level local MEG variability for the whole brain and for each functional network separately. We also obtained individual-level local MEG variability scores and 1/f exponent measures per network, akin to our fMRI analyses (see previous section). For each subject, we then averaged both estimates across brain regions to obtain whole-brain local MEG variability and 1/f exponent per person, that we could then correlate to each other. To expand on this whole-brain analysis, we built a hierarchical linear model where we accounted for regional heterogeneities in local MEG variability and 1/f exponent values (Fig. S7). Each region was entered in the model as a random effect to investigate region-level relationships between local MEG variability and the 1/f exponent.

Next, we simulated naturalistic neurophysiological time series at various arrhythmic motifs using the NeuroDSP toolbox91. Each simulation was 150s long, linearly combined rhythmic and arrhythmic (i.e., periodic and aperiodic) components at random initial phases, and was sampled at 500 Hz. The simulated periodic component of each time series consisted of a fixed alpha (peak frequency of 10 Hz, amplitude of 0.7 a.u., and band width of 2 Hz) and beta peak (peak frequency of 19 Hz, amplitude of 0.4 a.u., and bandwidth of 5 Hz). We simulated 10 time series at various aperiodic slope parameters ranging from −0.7 to −1.5 in steps of −0.1. The parameters for simulations were selected based on previous literature89. We used these simulations to test for the linear relationship between the steepness of the 1/f exponent and local MEG variability.

Lastly, to investigate multimodal correlates of local brain variability, we related group-level spatial maps of local fMRI BOLD variability (young sample 1 and 2) to group-level spatial maps of local MEG variability, via product-to-moment correlations.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study showed that BOLD signal variability represents a spatially heterogeneous, central property of multi-scale multi-modal brain organization, which is distinct from noise and \"rooted\" in electrophysiological processes. Local BOLD signal variability is defined as the regional moment-to-moment change in BOLD signal intensity between successive timepoints while the global BOLD signal variability is defined as the similarity in inter-regional functional connectivity over time. The introduction is well-written, which establishes the rationale for this study. However, several methodological and interpretational issues mitigate my enthusiasm for the current findings.\n1. The claim that “BOLD signal variability is rooted in electrophysiological processes” needs more evidence and support. First of all, we know that BOLD and MEG have different neural underpinnings. BOLD is related to local post-synaptic potentials. There is also neurovascular coupling to be taken into account. The current data in Figure 5C only demonstrate the negative association between group-level local BOLD signal variability and group-level local MEG signal variability via one spatial correlation (if I understood correctly). There is no discussion on this either. It would be more convincing to use concurrent fMRI-EEG data on the same group of subjects and show similar associations at the individual level.\n2. It is interesting to see that covSTATIS-derived global BOLD signal variability has a high test-retest reliability of 0.98 using two successive runs. The sliding window size was chosen to be 45-60s. I wonder how this global variability measure, its reliability, and age effect change with window size from 40s to 100s, which is usually used in dynamic functional connectivity literature. It would be good to provide some guidance on the optimal window size for such a measure.\n3. I also wonder why “global signal variability” rather than “connectivity variability” or something similar is used here. The measure is not exactly based on the BOLD signal itself but rather on the correlations between two time series, i.e., connectivity.\n4. Discussion: “This statistical explanation elucidates why, in this study, lower-order regions exhibited lower reliability of variability measures, than higher-order areas across fMRI samples.” It is quite confusing to come to this conclusion as all previous text on statistical explanation is related to local and global BOLD signal variability (e.g., which one corresponds closer to the fMRI data etc). No lower-order and higher-order brain regions were mentioned before this. Please clarify and discuss.\n5. What are the unique features of specific low-level and high-level brain regions in the context of local and global BOLD signal variability? The word “heterogenous” is used frequently but there is a lack of discussion on the patterns common or unique to local and global signal variability.\n6. It is nice to see a trend of how local and global BOLD variability relate to different gradients in young samples (Figure 3C). However, the level of evidence is a bit weak as some did not pass spin test in both young samples. Is there any reason behind this? For functional gradient, besides the principal unimodal-transmodal gradient, it might be informative to check the unimodal gradient (2nd one). Moreover, what is the pattern in the old sample? Do we see the same dissociation?\n7. What is the additional value of the local and global BOLD signal variability proposed here on top of or in comparison with some other BOLD-based measures like static FC (the authors have found high correlation between static FC and global variability), ALFF, other DFC measures. It would be useful to demonstrate its behavioral relevance.", "reviewer_2": "• Key results: The manuscript examines the role of signal variability in brain function. Using several large fMRI datasets, the authors look at the association between variability and dynamic functional connectivity with microstructure, transcriptomics, neurotransmitter receptor and metabolic data, static connectivity, and empirical and simulated MEG data. They suggest that variability is structured and spatially heterogeneous, and is important in multi-scale characterizations of brain function.\n• Validity: I do not see any fatal flaws, but the reasoning behind methodological choices and variables included needs to be better decribed. These choices are especially important as they can significantly impact outcomes in variability work (the authors themselves suggest this). Specifically, in the literature, variability and dynamic functional connectivity (dFC) tend to be treated somewhat differently. I see no issues with choosing a dFC-type measure to index “global variability”, but the reasoning behind the choice needs to be framed in the context previous literature in the paper. In addition, local and global complexity has been discussed quite a bit with EEG/MEG (see for example https://pubmed.ncbi.nlm.nih.gov/21525281/; https://pubmed.ncbi.nlm.nih.gov/23395850/). I am not sure why these papers/ideas are not discussed as they are relevant, and why, in the context of MEG, they don’t appear to be considered in analytic choices. Is this because the focus is more on variability than complexity in signal?\n• Originality and significance: Several papers have examined the importance/unique contribution of variability measures in understanding brain functon(seminal papers include https://www.pnas.org/doi/full/10.1073/pnas.0901831106; https://www.jneurosci.org/content/30/14/4914; https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3922711/). The value added in this manuscript is the impressive amount of data examined, the gamut of micro and macroscale variables that brain signal variability was linked to, and the examination of local and global variability in BOLD fMRI. To better highlight significance, the authors should focus on these unique contributions rather than the idea that variability is a useful metric. They should better motivate their choice of macro and microscale variables, their analysis methods (some which do not follow directly from previous literature), and include greater discussion of the implications of their specific findings (i.e., why and how is brain signal variability an important measure for brain characterisation? How can it fit with / add to other more commonly used measures? How should knowing that variability has “spatially heterogeneous topography, encapsulate micro-, meso- and macroscale neurobiological phenomena” influence the use of this measure in cognitive network neuroscience?).\n• Data & methodology: The data and methodology are sufficiently described in text (with references to more detailed previously published descriptions) to enable reproducing the results. As a minor point, added details regarding Group Prior Individual Parcellation (GPIP), specifically how it may influence results over more commonly used parcellation methods would be useful.\n• Appropriate use of statistics and treatment of uncertainties: The statistical tests appear to be performed appropriately.\n• Conclusions: Data interpretation and related conclusions are valid.\n• Suggested improvements: No additional data or experiments are necessary.\n• References: Some of the references mentioned above need to be added. In addition:\no The following sentence is missing references “The human brain is a complex stochastic system, hence integrating functional variability of brain signals is essential for the modeling of human brain function.” Deco et al’s (2009) paper from above speaks to this and could be added.\no The authors suggest that “computational approaches to brain dynamics may thus benefit from introducing spatial heterogeneity in local measures of variability, to more precisely approximate empirical observations.” Has this not been done? And if not, this is really interesting and would benefit from a little more discussion.\no The finding that global variability increases with age could be better contextualised within the literature, as I believe this is the area in which fMRI signal variability analyses are most commonly performed. There variability papers already cited in the manuscript show consensus with overall increases in variability with age, but an acknowledgement that there are some papers that found general decreases or a mix of increases and decreases (e.g., https://pubmed.ncbi.nlm.nih.gov/32915200/; https://pubmed.ncbi.nlm.nih.gov/35901557/) could be included.\n• Clarity and context: The abstract is clear and accessible but the introduction and discussion are overly concise, especially for an experimental design as complex as the one used in this paper."}, "author_response": "REVIEWER 1, SUMMARY: This study showed that BOLD signal variability represents a spatially\nheterogeneous, central property of multi-scale multi-modal brain organization, which is distinct from\nnoise and \"rooted\" in electrophysiological processes. Local BOLD signal variability is defined as the\nregional moment-to-moment change in BOLD signal intensity between successive timepoints while the\nglobal BOLD signal variability is defined as the similarity in inter-regional functional connectivity over\ntime. The introduction is well-written, which establishes the rationale for this study. However, several\nmethodological and interpretational issues mitigate my enthusiasm for the current findings.\nWe thank the reviewer for their support for our study and their constructive feedback.\nREVIEWER 1, COMMENT 1: The claim that “BOLD signal variability is rooted in electrophysiological\nprocesses” needs more evidence and support. First of all, we know that BOLD and MEG have different\nneural underpinnings. BOLD is related to local post-synaptic potentials. There is also neurovascular\ncoupling to be taken into account. The current data in Figure 5C only demonstrate the negative\nassociation between group-level local BOLD signal variability and group-level local MEG signal\nvariability via one spatial correlation (if I understood correctly). There is no discussion on this either. It\nwould be more convincing to use concurrent fMRI-EEG data on the same group of subjects and show\nsimilar associations at the individual level.\nRESPONSE TO R1, COMMENT 1: We appreciate the reviewer’s positive feedback on our research\nquestion and approach, and we thank them for raising this point. As correctly noted, we claimed that\nfMRI BOLD variability is rooted in electrophysiological processes based on a significant spatial\ncorrelation between group-level fMRI and MEG variability measures. The reviewer is correct in that\nsuch measures were obtained in two different young adult samples. As the reviewer mentioned, fMRI\nand MEG signals arise from distinct neurobiological sources, making their association non-trivial. We\nagree that the use of concurrent fMRI-EEG data could therefore provide more robust, and clearer,\nevidence in support of this multimodal association, ultimately strengthening our claim. Given that the\nMEG dataset used in this work also had an fMRI component (Taylor et al., 2017 NeuroImage), we\ndecided to use these openly available within-subject fMRI-MEG data in this revision to address the\nreviewer’s comment. In sum, we used these data to empirically determine whether such an association is\nobservable within-subject, rather than a single association between-group results.\nConsistent with the between-group findings, we now show a similar and reliable negative relationship\nbetween fMRI BOLD variability and MEG variability at the within-subject level. While the magnitude\nof the group-level association is similar to the one reported in-text using two different samples (r = -0.2),\nthe within-subject relationship is lower (average r = -0.14). Nonetheless, it still shows a large effect size\n(Cohen’s d = -1.17). We have now included these new results, their interpretation and details on the\ndataset in our revised manuscript and here below.\nRevised Results:\n“As a final step, we tested cross-modal relationships between local fMRI and MEG signal variability.\nWe report both cross-sample and within-sample associations (i.e., fMRI-MEG data obtained on those\nindividuals from the MEG Sample with available fMRI data – n=103). At the group level, we found a\nconsistent, negative association across and within sample between fMRI and MEG variability (Figure\n5C). A similar negative association was also observed at the individual level (average Fisher z-to-r = -\n1\n0.14; t(102) = -11.91, p<0.001, 95% C.I. = [-0.17, -0.12], Cohen’s d = -1.17). Together, these findings\nreinforce the biological multimodal nature of local brain signal variability.” (p. 22)\n2\nvariability was estimated as MEG signal variability via rMSSD, and as the 1/f exponent in the MEG power\nspectrum. (A) Functional topography of local rMSSD-derived MEG variability. Central, grey-scale map represents\nthe group average map. Around it are the 7 canonical network-level group maps. Note that rMSSD values were z-\nscored, as units are arbitrary. | (B) Left: Functional topography of the 1/f exponent per network, across individuals.\nMiddle: whole-brain relationship between local rMSSD-derived MEG variability and 1/f exponent across\nindividuals. Right: We simulated naturalistic electrophysiological timeseries and varied the steepness of their 1/f\nexponent at various parameters from -1.5 to -0.7, shown as positive values in the graph, in steps of -0.1. 10\nsimulations were run per parameter. We calculated local rMSSD-derived MEG variability on each simulated\ntimeseries, and related exponent values with rMSSD scores. | (C) Cross-modal spatial relationships between local\nfMRI and MEG brain variability, between two independent samples (left and middle figures) and within the same\nsample using fMRI-MEG data (right figure).\nRevised Discussion\n[…] Via empirical and simulated MEG data, we first found local signal variability to be driven by the\nslope of the signal power spectrum […].\nVia cross-sample and within-sample comparisons, we observed a consistent negative association\nbetween fMRI-derived and MEG-derived local brain signal variability. These results expand on previous\nliterature showing a spatial correspondence between fMRI and MEG signals using static functional\nconnectivity137-141. Our findings further demonstrate that fMRI and MEG signals share meaningful\ninformation about human brain dynamics. Here, the negative direction of fMRI-MEG correlations may\nbe explained by differences in the topographies of the two signals, as they originate from distinct\nphysiological sources and have different frequency content142-144. The relatively modest magnitude of\ntheir correlation may be due to differences in study design between the two modalities: MEG data were\nacquired during rest, while fMRI data during movie watching79. While further work is needed to shed\nmore light on the nature of these relationships, our work provides preliminary evidence in favor of the\nbiological multimodal nature of local brain signal variability. Specifically, these findings urge the general\nneuroscience community to re-consider what is empirically deemed “noise” (i.e., slope of the\nelectrophysiological power spectrum, local BOLD signal variability) as biologically meaningful\nproperties of multimodal brain signals. (pp. 28-29)\nRevised Methods:\nfMRI-MEG Sample\nWithin-individual fMRI-MEG data were obtained from the openly available Stage II of the Cambridge\nCentre for Aging and Neuroscience (Cam-CAN) data repository79,80. We included resting-state MEG\ndata from a total of 104 healthy young adults in the same age range 18-34 as our primary fMRI datasets\n(M = 28y, SD = 4y, 56% F). All MEG analyses, except for fMRI-MEG local signal variability\nage age\ncorrelations, were run on this sample. fMRI-MEG associations were conducted on 103 out of the 104\nindividuals, that is on those subjects with available resting-state MEG and movie-watching fMRI data.\nAll individuals provided written informed consent. Research was conducted in compliance with the\nHelsinki Declaration and was approved by the Cambridgeshire 2 Research Ethics Committee.\nMovie-watching fMRI data acquisition and processing\nParticipants were shown a shortened 8-minute version of the “Bang! You’re Dead” episode from the TV\nshow “Alfred Hitchcock Presents” (1961). None of them had seen the movie before. Participants were\ntold to pay attention to the movie while lying still in the scanner. Data were acquired on a Siemens 3T\nTIM Trio scanner using a 32-channel head coil. Functional data were collected with a multi-echo T2*\nEPI sequence with the following details: 193 total volumes, 32 axial slices, 3.7 mm thick, 0.74 mm gap,\n3 × 3 × 4.44 mm. Structural T1-weighted images were also used in this study for co-registration purposes.\nThese latter images were acquired with a 3D MPRAGE sequence with the following details: TR =\n2250 msec, TE = 2.99 msec, TI = 900 msec, flip angle = 9°, FOV = 256 × 240 × 192 m, 1 mm isotropic\nvoxels, GRAPPA acceleration factor = 2.\nFunctional data were preprocessed using fMRIPrep 21.0.1156 and tedana157, as reported in our\nprevious publication158. (pp. 33-35)\nREVIEWER 1, COMMENT 2: It is interesting to see that covSTATIS-derived global BOLD signal\nvariability has a high test-retest reliability of 0.98 using two successive runs. The sliding window size\nwas chosen to be 45-60s. I wonder how this global variability measure, its reliability, and age effect\nchange with window size from 40s to 100s, which is usually used in dynamic functional connectivity\nliterature. It would be good to provide some guidance on the optimal window size for such a measure.\nRESPONSE TO R1, COMMENT 2: We thank the reviewer for their comment, which has prompted us to\nwrite a separate paper—now published in Aperture Neuro (Baracchini et al., 2024 Aperture Neuro)—\nabout covSTATIS as a new tool for cognitive network neuroscience.\nTo address the reviewer’s specific comment on window size, we assessed the consistency and reliability\nof covSTATIS-derived global BOLD BOLD signal variability results using sliding windows of different lengths.\nSpecifically, we partitioned regional timeseries data from our Young Sample 1 (for both run 1 and 2) and\nour two Lifespan Samples into windows of 45sec, 60sec, 80sec and 100sec, and used covSTATIS to\nderive global BOLD signal variability scores for each iteration and sample. We used these global BOLD\nsignal variability scores to assess (1) their consistency across window lengths (all samples); (2) their test-\nretest reliability across runs for each window length (Young Sample 1, runs 1-2); and (3) the stability of\ntheir age effects across window lengths (Lifespan Samples 1-2).\nWe found covSTATIS-derived global BOLD variability scores to be (1) highly consistent across window\nlengths both at the individual level (intraclass coefficient estimates and their 95% C.I.s obtained via a\nmean-rating (k=4 window lengths), consistency, two-way mixed effects model with window size as fixed\neffect and individuals as random term: ICC > 0.82 for all samples; Supplemental Figure S1) and group\nlevel (Spearman rho between pairwise group-level window results across regions: 0.92 < r < 0.99, p <\n.001) in all four samples; (2) moderately reliable at the individual level for all window lengths (average\nr ~ 0.30, p<.001, large Cohen’s d > 1.6); and (3) highly stable in their age effects across window lengths\nfor both aging samples (PLS analysis with 4 conditions and 2 groups: LV1, at p<.001 with 58% variance\nexplained, largely recapitulating Figure 1C results both in terms of location of effects and their\ndirectionality; Supplemental Figure S2). We report below, and add to our main text, all details about\nthe results and analyses. To note, these new reliability results (2) are lower than the ones reported in\nFigure 1C because they were obtained within individuals across regions, and not on averaged group\nvalues across regions (i.e., what we show in Figure 1C). Since all our analyses are conducted at the group\nlevel, we kept our group-level findings in text and added our individual-level reliability results to\nSupplemental.\nTogether, these analyses show that covSTATIS is a robust tool for global BOLD signal variability (i.e.,\ndynamic functional connectivity) estimations. While we agree with the reviewer that it would be\nbeneficial to provide guidance on the optimal window size for covSTATIS, covSTATIS is not unlike\nother dynamic functional connectivity approaches, where optimal window size selection depends on TR\nand scan length. We therefore advise the reader to follow the guidelines already put forward by others in\nthis space (Liégeois et al., 2017 NeuroImage), as they apply to all methods that rely on sliding window\n5\napproaches. To reiterate, the advantage of covSTATIS over other dynamic functional connectivity\nmethods lies in the minimal number of user-dependent decisions post sliding window estimation (e.g.,\ncovSTATIS does not require any clustering). In this revision, we refer the reader to our new covSTATIS\npublication where we contextualize our method, comprehensively walk through all its details, and offer\na tutorial to facilitate its implementation using Open Data. Together, we appreciate the reviewer's\ncuriosity about the method, which prompted us to provide the details of our approach in an appropriate\nvenue for a fulsome presentation (which was outside of the scope of the present manuscript).\nRevised Results\n“Global BOLD signal variability was estimated as dynamic functional connectivity using covSTATIS.\nAs a three-way extension of Multidimensional Scaling and Principal Component Analysis, covSTATIS\nis a multi-table method that linearly combines multiple similarity matrices (i.e., correlation/covariance\ntables) and uses eigenvalue decomposition to identify predominant structured patterns both at the group\nand individual level58-60. In our case, we used a combination of sliding window approaches and\ncovSTATIS to examine, for each individual, the similarity of their windowed functional connectivity\nmatrices over time and derive individual-level estimates of regional functional connectivity dynamics.\nUnlike most conventional fMRI dynamic connectivity methods61 that operate on coarser brain dynamics\nestimates, covSTATIS defines finer regional measures of global dynamics and does so with minimal\nuser-dependent input, as it side-steps commonly adopted clustering approaches that have traditionally\nled to fractionated, study-specific definitions of dynamic functional connectivity62. For more information\nabout covSTATIS and its implementation, we refer the reader to our dedicated recent publication60.” (p.\n5-6)\n“We next validated our method by assessing its consistency, test-retest reliability and its ability to capture\nage effects typically reported in studies using traditional dynamic connectivity methods. We used resting-\nstate fMRI data from one sample of 145 healthy young adults with two successive runs, and two\nindependent cross-sectional healthy lifespan datasets (see Methods for details). For each sample and\nindividual, we partitioned their regional timeseries data into windows of different length and derived\ncovSTATIS’ area of the hull scores (i.e., global BOLD signal variability) for each iteration. We found\nglobal BOLD variability to be highly consistent across window lengths both at the individual level (two-\nway mixed effects ICC model > 0.82 for all samples; Supplemental Figure S1) and at the group level\n(0.92 < Spearman rho < 0.99, p < .001) in all samples. Global BOLD variability additionally showed\nhigh reliability at the group level (r=0.98; p<.001; Figure 1C) and moderate reliability at the individual\nlevel for all window lengths (average r ~ 0.30, p<.001, large Cohen’s d > 1.6). In line with current\nliterature demonstrating an age-induced dampening of the brain’s dynamic range66, here we used Partial\nLeast Squares67,68—a multivariate method that assesses the covariance between two or more sets of\nvariables—and showed that covSTATIS-derived area of the hull scores also decreased with age\nparticularly in regions that preferentially show age effects in the literature69 (Figure 1C; one significant\nlatent variable (LV1) at p=0.003 explaining 72% brain-age variance; Lifespan Sample 1 brain-age r=-\n0.19, Lifespan Sample 2 brain-age r=-0.30). Furthermore, age effects were stable across windows of\ndifferent length both in terms of their spatial location and directionality (LV1 at p<.001 explaining 58%\nbrain-age variance; Supplemental Figure S2). Together, these results highlight how covSTATIS is a\nvalid and robust tool to estimate global BOLD signal variability.” (p. 6-7)\n6\nBOLD signal variability (i.e., covSTATIS’ area of the hull values) for each fMRI sample across windows of\ndifferent length (45sec, 60sec, 80sec, 100sec; see Methods for full details). We assessed within-sample inter-\nwindow consistency via intraclass coefficient statistics (ICC) and their 95% confidence intervals. ICCs were\nderived via a mean-rating (k=4 window lengths), consistency, two-way mixed effects model with window size as\nfixed effect and individuals as random term. ICC values were greater than 0.82 for all samples, indicating high\nconsistency of covSTATIS estimates.\nFigure S2. CovSTATIS-derived age effects are stable across windows of different length. To assess the\nstability of covSTATIS age effects, we ran our PLS analyses contrasting covSTATIS values for four window\nlengths across the two Lifespan Samples (i.e., a 2-group 4-condition design; see Methods for full details). We\nfound age effects to be stable both in terms of their spatial location (brain plot on the left) and directionality (bar\nplot on the right; LV1 at p<.001 explaining 58% brain-age variance).\nRevised Methods\nGlobal fMRI BOLD signal variability estimation & validation\nEstimation\n“[…] As a three-way extension of Multidimensional Scaling and Principal Component Analysis,\ncovSTATIS is a multi-table method that linearly combines multiple similarity matrices (i.e.,\ncorrelation/covariance tables) and uses eigenvalue decomposition to identify predominant structured\npatterns both at the group and individual level58-60. In our case, we used a combination of sliding window\napproaches and covSTATIS to examine, for each individual, the similarity of their windowed functional\nconnectivity matrices over time and derive individual-level estimates of regional functional connectivity\ndynamics. First, we assessed, via the R similarity coefficient (i.e., squared product-to-moment\nv\ncorrelation)65, the similarity across all Y data tables across all individuals within each sample and stored\nthese values in a YxY R matrix. We then took the first eigenvector of the R matrix to generate weights\nv v\nfor each connectivity table and used these weights to build an NxN group compromise space—the linear\n8\ncombination of all matrices where regional connections more similar across time and individuals were\ngiven a higher weight since they were most represented in the sample. […]\nImportantly, covSTATIS is not unlike other dynamic functional connectivity approaches, where\nwindow size selection depends on TR and scan length. We therefore advise the reader to follow the\nguidelines already put forward by others in this space180, as they apply to all methods that rely on sliding\nwindow approaches. To reiterate, the advantage of covSTATIS over other dynamic functional\nconnectivity methods lies in the minimal number of user-dependent decisions post sliding window\nestimation (e.g., covSTATIS does not require any clustering). We refer the reader to our covSTATIS-\ndedicated publication where we contextualize our method, comprehensively walk through all its details,\nand offer a tutorial to facilitate its implementation using Open Data60.” (p. 39-40)\nValidation\n“To validate covSTATIS as a robust dynamic functional connectivity method, we partitioned regional\ntimeseries data from our Young Sample 1 (for both run 1 and 2) and our two Lifespan Samples into\nwindows of 45sec, 60sec, 80sec and 100sec, and used covSTATIS to derive global BOLD signal\nvariability scores for each iteration and sample. Table 2 provides a detailed breakdown of our windowing\nprocedure per sample and window size. After obtaining global BOLD signal variability scores, we used\nthese measures to assess (1) covSTATIS’ consistency across window lengths (all samples); (2) its test-\nretest reliability across runs for each window length (Young Sample 1, runs 1-2); and (3) its ability to\nrobustly capture age-related alterations in dynamic functional connectivity across window lengths\n(Lifespan Samples 1-2). All our validation analyses were carried out both at the individual and at the\ngroup level.\nTable 2. Breakdown of our windowing approach per window length and sample.\n45sec window 60sec window 80sec window 100sec window\nYoung Sample\nwindow, 186 window, 181 window, 174 window, 168\n1 run 1\nwindows windows windows windows\nYoung Sample\nwindow, 186 window, 181 window, 174 window, 168\n1 run 2\nwindows windows windows windows\nLifespan\nwindow, 275 window, 268 window, 258 window, 248\nSample 1\nwindows windows windows windows\nLifespan\nwindow, 368 window, 357 window, 341 window, 329\nSample 2\nwindows windows windows windows\nFor each sample, covSTATIS’ consistency was evaluated at the individual level via intraclass correlation\ncoefficients (ICC) and their 95% confidence intervals based on a mean-rating (k=4 window lengths),\nconsistency, two-way mixed effects model with window size as fixed effect and individuals as random\nterm182. At the group level, consistency was assessed by relating in a pairwise fashion window-specific\ncovSTATIS results across brain regions via Spearman rho statistics. At the individual level, covSTATIS’\ntest-retest reliability was determined for each window length via Fisher-z transformed Spearman rho\ncorrelations derived within subject, then averaged across individuals and finally converted back to\nsample-average rho estimates. We additionally performed t-tests contrasting Fisher-z transformed\n9\ncorrelations against zero and reported t-test values, along with their confidence intervals and Cohen’s d\neffect sizes. At the group level, covSTATIS’ test-retest reliability was instead assessed by correlating\ngroup-level covSTATIS scores across regions between runs.\nAs a final validation step, we implemented a multivariate analysis technique, Partial Least\nSquares (PLS)67,68,183, on our two Lifespan Samples, to test whether covSTATIS area of the hull measures\nwere sensitive to known age-related alterations in dynamic functional connectivity across the adult\nlifespan. In other words, we interrogated the covariance between covSTATIS measures and age, and we\nvalidated our results across windows of different lengths. Specifically, we ran two PLS analyses: (1)\ncontrasting covSTATIS’ area of the hull scores calculated on 45sec windows for both Lifespan Samples\n(i.e., a 2-group 1-condition design, these are the main results reported in Figure 1C); (2) contrasting\ncovSTATIS values for all four window lengths across the two Lifespan Samples (i.e., a 2-group 4-\ncondition design, these results are shown in Supplemental Figure S2).\nBriefly, PLS calculates a covariance matrix between two (or more) sets of variables. […] (pp. 40-\n42)\nREVIEWER 1, COMMENT 3: I also wonder why “global signal variability” rather than “connectivity\nvariability” or something similar is used here. The measure is not exactly based on the BOLD signal\nitself but rather on the correlations between two time series, i.e., connectivity.\nRESPONSE TO R1, COMMENT 3: We thank the reviewer for raising this point. While we agree that\n“connectivity variability” could have been a viable term for our covSTATIS-derived measures, it would\nhave added to the confusion in terminology already present in the field of brain dynamics. By using the\nterm “global signal variability” instead, we were able to unify two lines of work thus far treated\nindependently in the literature: BOLD signal variability and dynamic functional connectivity. Building\non the lay definition of the term “variability” (i.e.., any measure of change), we decided to create a\ncontrast in the spatial nature of our measures of interest by referring to BOLD signal variability as a local\ndynamic property and to dynamic functional connectivity as a global one, in line with previous work\n(Krohn et al., 2023 Science Advances). In doing so, we aimed at bringing more clarity to the field of\nBOLD dynamics, and at pushing researchers to jointly analyse local and global aspects of BOLD signal\ndynamics. We are now working on a separate manuscript that specifically aims to disambiguate the many\npossible ways of quantifying dynamical properties in BOLD time series. We expand on this in our\nresponse to the reviewer’s comment 7.\nIn this revision, we clarified our use of the word “global signal variability” and report changes below.\nRevised results\nQuantification of local and global BOLD signal variability\nWe first sought to identify robust regional metrics of local and global BOLD signal variability. By local\nand global BOLD signal variability, we refer to changes in the dynamic properties of single timeseries\n(local) and of pairs of timeseries (global). The former case is what is traditionally referred to as “BOLD\nsignal variability” in the literature, and the latter as “dynamic functional connectivity”. In using the\ngeneral term “variability” for both measures, we aimed at unifying two lines of work thus far treated\nindependently yet shown to be inherently intertwined18,44. (p. 5)\nREVIEWER 1, COMMENT 4: Discussion: “This statistical explanation elucidates why, in this study,\nlower-order regions exhibited lower reliability of variability measures, than higher-order areas across\nfMRI samples.” It is quite confusing to come to this conclusion as all previous text on statistical\n10\nexplanation is related to local and global BOLD signal variability (e.g., which one corresponds closer\nto the fMRI data etc). No lower-order and higher-order brain regions were mentioned before this. Please\nclarify and discuss.\nRESPONSE TO R1, COMMENT 4: We thank the reviewer for bringing to our attention this source of\nconfusion. We have now clarified the use of these terms in the results and discussion sections, specifically\nmaking a direct connection to the visionary work of Mesulam (Mesulam, 1998 Brain).\nRevised Results\n[…] For both local and global BOLD variability, reliability was lowest in sensori-motor (i.e., lower-\norder) areas and highest in heteromodal (i.e., higher-order), particularly default network, regions. (p. 9)\nRevised Discussion\n[…] Levels of statistical approximation are intuitively linked to levels of brain organization: lower levels\nof statistical approximation allow for more specific local neural representations taking place in unimodal\nlower-order regions, while higher levels of statistical approximation reflect lower-dimensional global\nneural representations residing in heteromodal higher-order areas. Building on Mesulam’s visionary\nwork114, this statistical explanation elucidates why, in this study, lower-order regions exhibited lower\nreliability of variability measures, than higher-order areas across fMRI samples. (p. 25)\nREVIEWER 1, COMMENT 5: What are the unique features of specific low-level and high-level brain\nregions in the context of local and global BOLD signal variability? The word “heterogenous” is used\nfrequently but there is a lack of discussion on the patterns common or unique to local and global signal\nvariability.\nRESPONSE TO R1, COMMENT 5: We thank the reviewer for raising this point. We purposefully chose\nnot to discuss common/unique patterns of local and global BOLD signal variability, as setting up such\ncontrast would require appropriate in-depth methodological testing that is outside of scope of the present\nwork. Tracing common and unique local-global variability patterns would require extensive\nmethodological analyses relating such patterns to multiple fMRI data types – to ensure reliability and\ngeneralisability, as our reliability analyses pointed out. It is for this reason that our Discussion section\nsolely focuses on interpreting findings that are reliable and directly derived from our statistical approach\n(i.e., assessing local and global BOLD signal variability independently). We appreciate the reviewer’s\ninterest, and we acknowledge the lack of clarity in our study’s goal. As such, we have now revised our\nDiscussion section accordingly. Future work among the co-authors is seeking to determine such relations,\nand we hope that the present manuscript inspires such investigations as well.\nRevised Discussion\n[…] Here, we showed how the statistical dependencies of local and global BOLD variability\nmeaningfully interact with the statistical principles that govern the brain’s spatial organization: cortical\nareas devoted to local processing presented variability patterns that were more dependent on fMRI data\ntype for both local and global variability measures, than higher-order regions. These results call for more\ntargeted empirical research on the complex interaction between local-global cortical processing, local-\nglobal BOLD signal variability and fMRI data type. While direct in-depth investigations on the\ncommonalities and differences between low- and high-level cortical areas in the context of local and\nglobal BOLD variability are needed, here we provide a first general demonstration of the multiscale\nfeatures of local and global BOLD signal variability, as independently assessed. Reliable local and global\nBOLD variability features are thus separately discussed here in turn.\nclosely map onto local and global information processes in the brain. […] (pp. 25-26)\nREVIEWER 1, COMMENT 6: It is nice to see a trend of how local and global BOLD variability relate\nto different gradients in young samples (Figure 3C). However, the level of evidence is a bit weak as some\ndid not pass spin test in both young samples. Is there any reason behind this? For functional gradient,\nbesides the principal unimodal-transmodal gradient, it might be informative to check the unimodal\ngradient (2nd one). Moreover, what is the pattern in the old sample? Do we see the same dissociation?\nRESPONSE TO R1, COMMENT 6: Part 1 (sample differences): We agree with the reviewer that evidence\non the multiscale mapping of local and global BOLD signal variability appears mixed, as spin test results\nare dependent on fMRI sample. We justify these mixed findings based on the coarse sample-dependent\ndifferences in the topography of local/global BOLD variability (highlighted in Figure 2 and Figure 3A-\nB). Throughout the paper, between-map correlations were run at the group level across brain regions:\nbetween-sample differences in correlation values are thus driven by differences in the respective maps’\nspatial distributions. Despite these inconsistencies, overall multiscale association patterns of local and\nglobal BOLD variability were largely maintained across samples, as shown by our Spearman’s rho and\ndominance analysis results (Figure 3A bottom panel and Figure 3B). Together, while our work serves as\na first-time comprehensive picture of BOLD variability’s neurobiological features, we acknowledge that\nmore research is needed to understand how such features relate to the brain’s spatiotemporal architecture\non an individual-by-individual basis. Resulting findings will be crucial to disentangle the sources of\nbetween-sample variance reported here. Following these considerations and comments, we have edited\nour manuscript accordingly and report changes below.\nPart 1 (Sample Differences): Revised Discussion\n[…] We found local and global BOLD variability to be spatially heterogeneous multiscale processes that\nclosely map onto local and global information processes in the brain. While some inconsistencies\nemerged across fMRI samples due to inter-sample differences in the topography of local and global\nBOLD signal variability, their multiscale properties were largely conserved across samples. (p. 26)\n[…] Despite the impact of fMRI acquisition sequences on local and global BOLD variability, the\ncentrality of these measures in multiscale brain organization was reliable across fMRI datasets, further\nestablishing their biological relevance. More research is needed to understand how these multiscale\nfeatures relate to the brain’s spatiotemporal architecture on an individual-by-individual basis, given that\nour analyses were conducted at the group level and therefore disregarded intra-sample sources of\nvariance. (p. 28)\nPart 2 (functional gradient 2): We thank the reviewer for suggesting considerations of more than the\nprincipal gradient. Following this suggestion, we probed the association between local and global BOLD\nsignal variability with the second functional gradient. While the second gradient presented noisier, less\nreliable results in predicting local and global BOLD variability across samples, broad patterns of\nassociation were largely maintained across the two gradients. We provide these new results in-text and\nextensively report them in our Supplemental Materials. Edits are presented here below.\nPart 2 (Functional Gradient 2): Revised Results\n[…] As a final step, we repeated all analyses above including the second gradient of functional\nconnectivity organization54. Broad patterns of association were largely maintained across the two\ngradients, yet the principal gradient showed more reliable inter-sample results (Figure S6). (p. 18)\n12\nsecond functional gradient. | (A) Correlation matrices for each fMRI sample (upper triangle: Young Sample 1;\nlower triangle: Young Sample 2). Asterisks indicate correlations that survived significance testing (10,000\nHungarian spins of Schaefer’s regional labels). We found positive associations between local BOLD signal\nvariability and the second functional connectivity gradient (r=0.54, p =0.03 Young Sample 1), and global\n10k spin\nBOLD signal variability and the second gradient (r=0.70, p <0.001 Young Sample 1). As expected,\n10k spin\nparticipation coefficient scores of local and global BOLD variability including the second gradient largely\noverlapped with the values obtained by including only the first gradient (Young Sample 1: local = 0.62, global =\n0.66; Young Sample 2: local = 0.49, global = 0.57). Additionally, in probing the inter-sample convergence of\nmultiscale variability patterns, we found a significant inter-sample correlation only for local BOLD variability\n(r=0.64, p=0.007, 95% C.I. [0.22, 0.86]) and not for global BOLD variability (r=0.49, p=0.05, 95% C.I. [-0.001,\n0.80]), differently from what we reported on the first gradient alone. | (B) Dominance analysis results per metric\nand fMRI sample including the second functional gradient.\nPart 3 (Aging): We appreciate the reviewer’s interest in understanding how multiscale patterns present\nthemselves in older samples. We agree it is a valuable question for the field, especially given the number\nof investigations on age-related alterations in BOLD signal variability (e.g., Nomi et al., 2017; Boylan et\nal., 2020; Garrett et al., 2010; Garrett, Kovacevic, et al., 2013; Garrett et al., 2015, 2017, 2021; Goodman\net al., 2024; Grady & Garrett, 2014; Millar, Ances, et al., 2020; Armbruster-Genç et al., 2016; Grady et\nal., 2023), and age-related differences in the gradient architecture of the brain (e.g., Bethlehem et al.,\n2022; Setton et al., 2023). However, a systematic and robust investigation on this topic would require\nindividual-level data and analyses that, as mentioned in previous comments, are outside the scope of the\npresent study.\nREVIEWER 1, COMMENT 7: What is the additional value of the local and global BOLD signal\nvariability proposed here on top of or in comparison with some other BOLD-based measures like static\nFC (the authors have found high correlation between static FC and global variability), ALFF, other\nDFC measures. It would be useful to demonstrate its behavioral relevance.\nRESPONSE TO R1, COMMENT 7: We thank the reviewer for their enthusiasm and interest. These are\nfantastic questions that we (and others) are currently exploring in separate lines of work. Ultimately,\nintegrating these projects in the current study would be outside of scope, as it would surpass the journal’s\nword limit and additionally overload the reader, given how dense the paper is already. As such, we touch\non these lines of work here below and give attention to this important issue in the discussion section.\nFirst, we would like to address the issue of the metrics provided in the current manuscript. Overall, the\nvalue provided by our measures of local and global variability is their relatively straightforward\nimplementation, their mathematical simplicity (i.e., rMSSD captures processes related to basic\ndistributional properties of a timeseries), and user-independent nature. We expect our rMSSD findings\nto hold across other basic timeseries measures (i.e., standard deviation), and to relate to more complex\nmeasures given the slow, linear nature of fMRI (i.e., ALFF). We also expect our covSTATIS findings to\nbe recapitulated by more standard sliding window approaches, given their relatively similar\nimplementation.\nUnderstanding inter-metric differences in BOLD variability is however a huge endeavour, as the field\ncurrently lacks a systematic framework that allows for inter-metric comparisons. To address this gap, in\nan ongoing project, we are putting forward a logically structured hierarchical grouping of BOLD\nvariability measures. Such grouping organises all different metrics used in the literature in separate\ncategories based on their mathematical implementation (distributional measures such as rMSSD used in\nthis paper, linear – such as ALFF/fALFF – and non-linear measures, measures of stationarity and non-\n14\nstationarity). In doing so, our framework will help researchers better understand how existing methods\nfit in the broader literature, encouraging appropriate usage and interpretation of BOLD variability metrics\n(see Draft Figure below). In this context, inter-category comparisons, such as the one suggested by the\nreviewer, lose meaning if thought as comparisons between similar concepts: different categories instead\ncapture different timeseries properties. This work was presented at OHBM 2025 in Brisbane (poster\n1575, abstract available through Aperture Neuro here: https://zenodo.org/records/15641972).\nAs a technical side note to directly address the reviewer’s comment, (f)ALFF and variance (and\ntherefore SD – the square root of variance – and rMSSD – r > 0.9 to SD on fMRI data) are highly similar\nin fMRI data. Mathematically speaking, ALFF measures are retrieved from the Autocorrelation Function\n(ACF). According to the Wiener-Khinchin theorem, the ACF contains all power spectral information of\nthe data (stationary time series). As such, the area under the curve of the data’s power spectral density\n(PSD) equals to the data’s variance. If the data are filtered, as it is standard practice in fMRI, then the\narea under the PSD is mostly composed of low frequency fluctuations, precisely what ALFF is\nmeasuring.\nTo add to this theoretical work, we are also empirically testing for inter-metric differences in\nBOLD variability using synthetic timeseries. This work was presented at OHBM 2024. The abstract\n(#2415) is available on Aperture Neuro (Organization for Human Brain Mapping. Abstract Book 6:\nOHBM 2024 Annual Meeting. Aperture Neuro. 2024;4(Suppl 1)).\nWe would like to acknowledge the work of Wehrheim and colleagues on this topic as well\n(Wehrheim et al., 2024 Human Brain Mapping).\nDraft Figure of the framework we are proposing in ongoing work, and we presented at OHBM 2025. This work is\nan in-prep manuscript, further highlighting the importance of the issue the reviewer raised.\nbehavior and appreciate the reviewer's interest in this relationship. There is a rich literature on the\nbehavioral relevance of BOLD signal variability. In fact, the field of BOLD signal variability started\nfrom asking precisely these types of brain-behavior questions. Collectively, the field has demonstrated\nhow BOLD variability provides unique behavioral information across a variety of trait (e.g., domains\nof cognition) and state (e.g., emotion) variables in healthy young adults, aging samples, and clinical\npopulations. The first author of the present paper extensively reviews this literature in her PhD thesis,\nwhich has just recently been deposited and is accessible online\n(https://www.proquest.com/openview/a8c1eeca7ba5c433c1caa187ff142c96/1?pqorigsite=gscholar&cbl\n=18750&diss=y; Baracchini, 2024 McGill University). In the present submission, our goal was to\nprovide a neurobiological framework to help contextualise, interpret, and guide existing and future\nbehavioral research on the topic. We argue in favor of more investigations into BOLD variability's\nrelationship with behavior.\nRevised Discussion\n[…] Finally, we acknowledge the wide array of brain signal variability metrics in the fMRI and MEG\nliteratures18,103,145, yet in the absence of a systematic framework for inter-metric comparisons, this study\nfocused on two of the most popular signal variability metrics (i.e., rMSSD and dynamic functional\nconnectivity). A logically structured grouping of variability measures is necessary to disambiguate the\nontological confusion in the field and enable robust and interpretable inter-metric comparisons. Different\nvariability metrics are derived from different mathematical formulations with inherently different\nassumptions: in other words, variability measures have different epistemologies and therefore they lead\nto distinct ontologies146. In the presence of such framework, future research should test how the reported\nmultiscale patterns of local and global brain signal variability may depend on the choice of variability\nmetric.\nOverall, the value provided by our measures of local and global variability is their relatively\nstraightforward implementation, their mathematical simplicity (i.e., rMSSD captures processes related\nto basic distributional properties of a timeseries), and user-independent nature. We expect our rMSSD\nfindings to hold across other basic timeseries measures (i.e., standard deviation), and to relate to more\ncomplex measures given the slow, linear nature of fMRI (i.e., (f)ALFF147). We also expect our\ncovSTATIS findings to be recapitulated by more standard sliding window approaches64, given their\nrelatively similar implementation. (p. 29)\nREVIEWER 2, SUMMARY: Key results: The manuscript examines the role of signal variability in brain\nfunction. Using several large fMRI datasets, the authors look at the association between variability and\ndynamic functional connectivity with microstructure, transcriptomics, neurotransmitter receptor and\nmetabolic data, static connectivity, and empirical and simulated MEG data. They suggest that variability\nis structured and spatially heterogeneous and is important in multi-scale characterizations of brain\nfunction.\nWe thank the reviewer for their support for our study and their constructive feedback.\nREVIEWER 2, COMMENT 1: Validity: I do not see any fatal flaws, but the reasoning behind\nmethodological choices and variables included needs to be better described. These choices are especially\nimportant as they can significantly impact outcomes in variability work (the authors themselves suggest\nthis). Specifically, in the literature, variability and dynamic functional connectivity (dFC) tend to be\ntreated somewhat differently. I see no issues with choosing a dFC-type measure to index “global\nvariability”, but the reasoning behind the choice needs to be framed in the context previous literature in\nthe paper. In addition, local and global complexity has been discussed quite a bit with EEG/MEG (see\nfor example https://pubmed.ncbi.nlm.nih.gov/21525281/; https://pubmed.ncbi.nlm.nih.gov/23395850/).\nI am not sure why these papers/ideas are not discussed as they are relevant, and why, in the context of\nMEG, they don’t appear to be considered in analytic choices. Is this because the focus is more on\nvariability than complexity in signal?\nRESPONSE TO R2, COMMENT 1: We thank the reviewer for raising this important point and for\nbringing to our attention key references in this space. In this revision, we have better contextualised our\nchoices of local and global BOLD signal variability metrics within the broader literature. We report\nchanges below.\nWe additionally address the issue of complexity in this revision: in line with our response to reviewer\n1’s comments (see comment #7), we have expanded on the distinction between variability and other\nmeasures of signal dynamics, such as signal complexity. As correctly pointed out by the reviewer here,\nour research question focused on the biological nature of BOLD signal variability, and not of BOLD\nsignal complexity. Measures of signal variability and signal complexity capture different aspects of a\ntimeseries – whether it is on fMRI or MEG data, as these measures are derived from distinct mathematical\nformulations with inherently different assumptions (e.g., linearity vs non-linearity). The field currently\nlacks a systematic framework to allow for comparisons between these measures: in its absence, there\nremains significant ontological confusion around this topic. As such, we purposefully decided to solely\nfocus on measures of signal variability (for both fMRI and MEG). However, in an ongoing project, we\nare putting forward a logically structured hierarchical grouping of time-series measures. Such grouping\norganizes all different metrics used in the BOLD signal variability literature in separate categories based\non their mathematical implementation (distributional measures such as rMSSD used in this paper, linear\n– such as ALFF/fALFF – and non-linear measures, such as complexity measures). In doing so, our\nframework will help researchers better understand how existing methods fit in the broader literature,\nencouraging appropriate usage and interpretation of BOLD variability metrics. Inter-category\ncomparisons, such as the one suggested by the reviewer, lose meaning if thought as comparisons between\nsimilar concepts: different categories instead capture different timeseries properties. This work was\npresented at OHBM 2025 in Brisbane (poster 1575, abstract available through Aperture Neuro here:\nhttps://zenodo.org/records/15641972). As such, we touch on these lines of work here below and give\nattention to this important issue in the discussion section.\n[…] Second, human brain function is organized along hierarchical modules, ranging from local\nfunctional units to global functional networks35,36, with heterogeneous topographies. Such organization\nenables information to be both segregated (local) and integrated (global) across units, subserving\ndifferent neural and cognitive functions – as observed in various measures such as mean functional\nactivation, functional connectivity, computational and BOLD signal complexity studies37-42. Along these\nlines, BOLD signal variability must therefore also be understood within a local-global framework43,44\nthat considers both the methodology and the analysis level of interest. While previous fMRI reports have\ninvestigated the local-global nature of BOLD signal variability, they have generally focused on static\nglobal measures of brain organisation (i.e., the correlation between local mean BOLD signals over time,\nfunctional connectivity45) and more local characterization of variability on a single scale (e.g., voxel or\nregion). The temporal variability of individual regions (local BOLD signal variability) must thus be\nevaluated alongside the dynamic coordination of functional neural units (dynamic functional\nconnectivity, here global BOLD signal variability). If biologically relevant, local and global BOLD\nsignal variability should present spatially heterogeneous topographies that recapitulate known\nneurobiological processes that unfold across spatial scales. (p. 4)\nRevised results\nQuantification of local and global BOLD signal variability\nWe first sought to identify robust regional metrics of local and global BOLD signal variability. By local\nand global BOLD signal variability, we refer to changes in the dynamic properties of single timeseries\n(local) and of pairs of timeseries (global). The former case is what is traditionally referred to as “BOLD\nsignal variability” in the literature, and the latter as “dynamic functional connectivity”. In using the\ngeneral term “variability” for both measures, we aimed at unifying two lines of work thus far treated\nindependently yet shown to be inherently intertwined18,44. (p. 5)\nRevised Discussion\n[…] Finally, we acknowledge the wide array of brain signal variability metrics in the fMRI and\nMEG literatures18,103,145, yet in the absence of a systematic framework for inter-metric comparisons, this\nstudy focused on two of the most popular signal variability metrics (i.e., rMSSD and dynamic functional\nconnectivity). A logically structured grouping of variability measures is necessary to disambiguate the\nontological confusion in the field and enable robust and interpretable inter-metric comparisons. Different\nvariability metrics are derived from different mathematical formulations with inherently different\nassumptions: in other words, variability measures have different epistemologies and therefore they lead\nto distinct ontologies146. In the presence of such framework, future research should test how the reported\nmultiscale patterns of local and global brain signal variability may depend on the choice of variability\nmetric.\nOverall, the value provided by our measures of local and global variability is their relatively\nstraightforward implementation, their mathematical simplicity (i.e., rMSSD captures processes related\nto basic distributional properties of a timeseries), and user-independent nature. We expect our rMSSD\nfindings to hold across other basic timeseries measures (i.e., standard deviation), and to relate to more\ncomplex measures given the slow, linear nature of fMRI (i.e., (f)ALFF147). We also expect our\n18\ncovSTATIS findings to be recapitulated by more standard sliding window approaches64, given their\nrelatively similar implementation. (p. 29)\nREVIEWER 2, COMMENT 2: Originality and significance: Several papers have examined the\nimportance/unique contribution of variability measures in understanding brain function (seminal\npapers\ninclude https://www.pnas.org/doi/full/10.1073/pnas.0901831106; https://www.jneurosci.org/content/30\n/14/4914;https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3922711/). The value added in this\nmanuscript is the impressive amount of data examined, the gamut of micro and macroscale variables\nthat brain signal variability was linked to, and the examination of local and global variability in BOLD\nfMRI. To better highlight significance, the authors should focus on these unique contributions rather\nthan the idea that variability is a useful metric. They should better motivate their choice of macro and\nmicroscale variables, their analysis methods (some which do not follow directly from previous\nliterature), and include greater discussion of the implications of their specific findings (i.e., why and\nhow is brain signal variability an important measure for brain characterisation? How can it fit with /\nadd to other more commonly used measures? How should knowing that variability has “spatially\nheterogeneous topography, encapsulate micro-, meso- and macroscale neurobiological phenomena”\ninfluence the use of this measure in cognitive network neuroscience?).\nRESPONSE TO R2, COMMENT 2: We thank the reviewer for their valuable feedback and detailed\nsuggestions. While we agree with the reviewer that other important studies have examined the functional\nrelevance of brain signal variability, its characterisation in fMRI remains limited. Existing fMRI reports\nhave shown relationships between local BOLD variability and local BOLD activity (Garrett et al., 2011),\na region’s local variability and its whole brain functional connectivity (Garrett et al., 2018), local\nvariability and non-grey matter tissue (Millar, Petersen, et al., 2020), and local variability and system-\nlevel dopaminergic and GABAergic availability (Garrett et al., 2015; Lalwani et al., 2021). Most of these\nstudies however are carried out in aging populations, underscoring the need for broader examinations in\nyounger samples. Furthermore, while the majority of BOLD signal variability studies focus on\nvariability’s behavioural, cognitive, and clinical correlates, these associations do not directly establish its\nrole in brain organisation.\nFollowing this reviewer's suggestion, in this revision and provided below, we revised the manuscript to\nacknowledge the key issues raised regarding our study design, and we have included greater discussion\non the significance and implications of our findings. We acknowledge that this issue was raised by both\nreviewers, and believe the current revision appropriately contextualizes the current findings.\nRevised Introduction\n[…] The human brain is a complex system. From the processing of incoming information to the\ngeneration of motor outputs, variability is present at all levels of the central nervous system6,7. Brain\nactivity varies across multiple timescales, from milliseconds to years, and these temporal changes can be\nobserved across multiple spatial scales, from regions to networks8,9. Despite its pervasiveness, there\ncontinues to be resistance to exploring variability in human cognitive neuroscience. Human behavior has\nbeen shown to be stochastic7,10,11 and computational models have operationalized the brain as a complex\ndynamical system12-16, yet corresponding empirical neuroimaging research still lags behind. (p. 3)\n[…] Thus far, BOLD signal variability has been primarily studied in relation to behavior, cognition,\ndevelopment, and clinical status18-24. Demonstrating such associations does not, however, directly\nestablish BOLD signal variability’s biological role in whole-brain organization. As such, it remains\n19\nunclear whether and to what extent BOLD signal variability investigations are capturing meaningful\nneural signatures or system noise. To this end, a few recent studies have begun to characterize BOLD\nvariability’s neurobiological underpinnings, yet these investigations are typically restricted to\ncomparisons across age and focus on single, mostly macroscale, neurobiological features25-29. A broader\nand more integrated assessment of BOLD signal variability in younger samples is therefore needed, one\nthat examines its statistical, topographical and neuronal properties and systematically relates BOLD\nsignal variability to a wider array of neurobiological features. Evidence in favour of taking such an\nintegrative multi-scale approach to BOLD signal variability comes from multiple lines of empirical and\ncomputational work collectively demonstrating the role of cellular, molecular, genetic, and metabolic\nfactors in shaping local and global hemodynamic signals30-34. (pp. 3-4)\n[…] To understand the spatial organization and biological properties of local and global BOLD\nvariability, inspired by recent multiscale fMRI investigations32,46-49, we examined their topography\nwithin each fMRI dataset and interrogated associations with open-source data, including ex_vivo\nhistology50 and in_vivo microstructure51, transcriptomics52-53, PET-derived neurotransmitter receptor and\nmetabolic information47, and fMRI static connectivity data54,55. (pp. 4-5)\nRevised Discussion\nResponse to [Why and how is brain signal variability an important measure for brain characterisation?]\n[…] Brain signal variability has been increasingly recognized as a fundamental feature of neural\nfunctioning. At the cellular and molecular level, studies have shown how neural variability serves as\npreparatory activity for future stimulus processing. It allows for the recapitulation of previous sensory\nexperience97-99, predicts ongoing behavioral and cognitive states100-104, facilitates synaptic\nformation100,105-107 and drives signal propagation108. At the macroscale level in fMRI, accumulating\nevidence has highlighted BOLD signal variability as a key feature of optimal behavioral functioning and\na marker of neurocognitive aging and clinical status29,21,25-27,29,92-95,109-113. Our work builds on these\nmultiple lines of evidence and further expands the relevance of brain signal variability from a multiscale\nneurobiological standpoint. (pp. 24-25)\nResponse to [How can it fit with / add to other more commonly used measures?]\n[…] Overall, the value provided by our measures of local and global variability is their relatively\nstraightforward implementation, their mathematical simplicity (i.e., rMSSD captures processes related\nto basic distributional properties of a timeseries), and user-independent nature. We expect our rMSSD\nfindings to hold across other basic timeseries measures (i.e., standard deviation), and to relate to more\ncomplex measures given the slow, linear nature of fMRI (i.e., (f)ALFF147). We also expect our\ncovSTATIS findings to be recapitulated by more standard sliding window approaches64, given their\nrelatively similar implementation. (p. 29)\nResponse to [How should knowing that variability has “spatially heterogeneous topography,\nencapsulate micro-, meso- and macroscale neurobiological phenomena” influence the use of this measure in\ncognitive network neuroscience?]\n[…] By establishing BOLD signal variability as a spatially heterogeneous, multifactorial, multimodal\nproperty of brain organization integral to healthy brain function, this study underscores the value of\nsystems-level approaches to brain function148. Our work provides fertile soil for new theoretical and\nmethodological work in cognitive network neuroscience, by opening the door to neurobiologically\n20\ngrounded hypotheses on the role of local and global BOLD signal variability in cognition and behavior.\nMoreover, our findings offer empirically grounded targets for generative models of brain signal dynamics\nin health and disease. (p. 30)\nREVIEWER 2, COMMENT 3: Data & methodology: The data and methodology are sufficiently\ndescribed in text (with references to more detailed previously published descriptions) to enable\nreproducing the results. As a minor point, added details regarding Group Prior Individual Parcellation\n(GPIP), specifically how it may influence results over more commonly used parcellation methods would\nbe useful.\nRESPONSE TO R2, COMMENT 3: We thank the reviewer for their positive feedback on our data and\nmethodology. We appreciate their suggestion to provide more detail on the GPIP approach, and we have\nnow clarified this in the revised manuscript. Edits are provided here below.\nBriefly, the goal of our work was to understand how neural signal dynamics contribute to the rich spatial\narchitecture of the human brain. As such, we wanted to ensure high topographical precision of fMRI\nBOLD. Individualised parcellation approaches, including GPIP, precisely enable this. While standard\ngroup-level parcellation methods are more commonly used in fMRI research, they often yield coarse and\nimprecise estimations of individual-level regional boundaries and therefore of the spatial unfolding of\nfMRI BOLD. Compared to standard group-level parcellation methods, individualised parcellation\napproaches offer substantially improved sensitivity to individual differences in brain activity and\nbehaviour, by adjusting each individual’s regional boundaries to their own functional connectivity\norganisation (Chong et al., 2017 NeuroImage; Kong et al., 2018 Cerebral Cortex; Levi et al., 2023\nNetwork Neuroscience; Sassenberg et al., 2023 NeuroImage). To note, we selected GPIP as our\nparcellation method of choice also for consistency with previous work from our group using the same\ndata (Setton et al., 2023 Cerebral Cortex).\nRevised Methods\nfMRI Young Sample 1\n[…] Functional images were parcellated using the Group Prior Individual Parcellation\n(GPIP)69,149, a participant-specific parcellation approach initialized on the Schaefer 200-17 network\nsolution56. We selected GPIP as our parcellation method of choice to ensure high topographical\nprecision of fMRI BOLD, and for consistency with previous work from our group using the same\ndata69. Unlike standard parcellations, GPIP accounts for within-subject variation in parcel boundaries,\noffering substantially improved sensitivity to individual differences in brain activity and behaviour149-\n152. Specifically, GPIP consists of an iterative Bayesian process: it starts from the user-provided\nreference group atlas (Schafer 200-17 in our case) and iteratively refines each individual’s parcel\nboundaries relative to their resting-state connectivity data149. It then computes concentration matrices\n(inverse covariance/partial correlation) for the whole sample via a group sparsity constraint. To ensure\nstability and precision, the algorithm iterates between the former and latter step for a total of 20\niterations or until no more than one vertex is changing per parcel149. For more details about GPIP, the\nreader is referred to our previous publication where we used GPIP on this sample69. (pp. 31-32)\nREVIEWER 2, COMMENT 4: Appropriate use of statistics and treatment of uncertainties: The\nstatistical tests appear to be performed appropriately.\napproach.\nREVIEWER 2, COMMENT 5: Conclusions: Data interpretation and related conclusions are valid.\nRESPONSE TO R2, COMMENT 5: We appreciate the reviewer’s positive feedback on our conclusions.\nREVIEWER 2, COMMENT 6: Suggested improvements: No additional data or experiments are\nnecessary.\nRESPONSE TO R2, COMMENT 6: We thank the reviewer for their positive assessment of our study.\nREVIEWER 2, COMMENT 7: References: Some of the references mentioned above need to be added. In\naddition:\no The following sentence is missing references “The human brain is a complex stochastic system, hence\nintegrating functional variability of brain signals is essential for the modeling of human brain function.”\nDeco et al’s (2009) paper from above speaks to this and could be added.\no The authors suggest that “computational approaches to brain dynamics may thus benefit from\nintroducing spatial heterogeneity in local measures of variability, to more precisely approximate\nempirical observations.” Has this not been done? And if not, this is really interesting and would benefit\nfrom a little more discussion.\no The finding that global variability increases with age could be better contextualised within the\nliterature, as I believe this is the area in which fMRI signal variability analyses are most commonly\nperformed. There variability papers already cited in the manuscript show consensus with overall\nincreases in variability with age, but an acknowledgement that there are some papers that found general\ndecreases or a mix of increases and decreases\n(e.g., https://pubmed.ncbi.nlm.nih.gov/32915200/; https://pubmed.ncbi.nlm.nih.gov/35901557/) could\nbe included.\nRESPONSE TO R2, COMMENT 7: We thank the reviewer for their suggestions and comments.\n[Point 1]: We agree with the reviewer and have now incorporated the suggested reference in our revised\nmanuscript.\nRevised Discussion\n[…] The human brain is a complex stochastic system16, hence integrating functional variability of brain\nsignals is essential for the modeling of human brain function. (p. 24)\n[Point 2]: We thank the reviewer for their insightful comment. We agree that our sentence on\ncomputational approaches to local dynamics would benefit from a bit more discussion.\nComputational modelers typically favor model parsimony, that is model simplicity and tractability, over\nhigh biological precision. For this reason, early population models assumed and enforced spatial\nhomogeneity in local dynamics (Wilson and Cowan, 1972 Biophysical Journal; Jansen and Rit, 1995\nBiological Cybernetics; Robinson et al., 1997 Physical Review E; Jirsa and Haken, 1997 Physica D:\nNonlinear Phenomena). Only in recent years, modelers have begun integrating spatial heterogeneity –\nfrom receptor density, gene expression, myelination, cortico-thalamic functional gradient maps – in their\nmodels, to more accurately capture key features of brain dynamics observed empirically, including\n22\ncovariance, dynamic functional connectivity, regional times-scales, and stimulus sensitivities (Deco et\nal., 2013 Journal of Neuroscience; Breakspear, 2017 Nature Neuroscience; Wang et al., 2020 Nature\nReview Neuroscience; Herzog et al., 2024 Network Neuroscience; Froudist-Walsh et al., 2021 Neuron;\nKlatzmann et al., 2025 Cell Reports; Müller et al., 2023 Cell Reports; Müller et al., 2020 Nature\nCommunications). It is in this context that we argue that our results will help shape existing and future\nneurobiologically-informed computational models, by further constraining them with a new spatially\nheterogenous axis of brain dynamics.\nIn response to the reviewer’s suggestion, we have now expanded our discussion to clarify our point and\ncontextualise our argument in relation to past and current modeling approaches.\nRevised Discussion\n[…] Our study adds to this body of evidence by revealing that local BOLD variability is a structured,\nspatially heterogeneous process. Computational approaches to brain dynamics may thus benefit from\nintroducing spatial heterogeneity in local measures of variability, to more precisely approximate\nempirical observations. In doing so, empirical fMRI data will not only serve as a model validation tool\nbut also as a model optimization tool. While early population models assumed and enforced spatial\nhomogeneity in local dynamics120-123, recent years have progressively seen the integration of spatial\nhomogeneity – from receptor density, gene expression, myelination, cortico-thalamic functional gradient\nmaps – in neural mass models, to more accurately capture key features of brain dynamics observed\nempirically, including covariance, dynamic functional connectivity, regional times-scales, and stimulus\nsensitivities14,15,34,124-128. It is in this context that we argue that our results will help shape existing and\nfuture neurobiologically-informed computational models, by further constraining them with a new\nspatially heterogenous axis of brain dynamics. (p. 26)\n[Point 3]: We thank the reviewer for their suggestion. In our manuscript, we found age-related decreases\n(and not increases) in covSTATIS-derived global signal variability, in line with most BOLD signal\nvariability studies in the literature. The reviewer is correct in saying that BOLD variability investigations\nare commonly carried out in aging populations, as aging researchers were among the first to investigate\nBOLD signal variability as a potential meaningful construct in fMRI (e.g., Randy McIntosh, Cheryl\nGrady and Douglas Garrett). In this revision, in line with the reviewer’s suggestion, we have briefly\ncontextualized our aging-related findings. As these results primarily served as a validation for the\ncovSTATIS method, their contextualization remains concise in our revised discussion.\nRevised Discussion\n[…] Local and global BOLD signal variability are not merely sources of biologically irrelevant noise.\nThey are reduced with age in line with most prior investigations19-21,25-27,29,92-94, though some studies have\nreported opposite effects95,96. We have shown how measures of BOLD signal variability have a spatially\nheterogeneous topography, they encapsulate micro-, meso- and macro-scale neurobiological phenomena,\nand are related to underlying electrophysiological neuronal activity. (p. 24)\nREVIEWER 2, COMMENT 8: Clarity and context: The abstract is clear and accessible but the\nintroduction and discussion are overly concise, especially for an experimental design as complex as the\none used in this paper.\nRESPONSE TO R2, COMMENT 8: We thank the reviewer for their feedback. In responding to all of\ntheir comments and to the comments of reviewer 1, we have expanded the introduction and discussion\n23\nsections of our paper to better reflect the complexity of our experimental design and to provide additional\ncontext for our results. We aimed to strike a balance between general audience accessibility and scientific\nfidelity. We appreciate both reviewers for their insight and comments.\nReviewer #1 (Remarks to the Author):\nI appreciate the hard work from the authors. The revision has addressed some concerns, but there\nare still three main concerns left:\nREVIEWER 1, COMMENT 1: What are the unique features of specific low-level and high-level brain\nregions in the context of local and global BOLD signal variability? The authors did not perform any\nadditional analyses to compare them quantitatively. I think this is important even though they are\nderived independently. The newly added discussion is also non-specific and not satisfying. E.g.\n“cortical areas devoted to local processing presented variability patterns that were more dependent\non fMRI data type for both local and global variability measures, than higher-order regions.” What\nare the cortical areas devoted to local processing? How is it different from higher-order regions?\nGiven the rich literature on functional gradients and organization involving sensorimotor and\nassociation regions, more comparison and interpretation for local and global variability is needed.\nTo me, this response contradicts the authors’ claim in their earlier response that these two processes\nare linked, and in this study, they purposely named them differently from the literature to study both\n(single time-series variability and paired time series (connectivity) variability. Without fully\naddressing this point, providing additional insights for the field is hard.\nWe thank the reviewer for their constructive feedback. We agree that quantitively comparing low-\nlevel and high-level brain regions in the context of local-global signal variability metrics is an\nimportant direction for future work. Yet, the present study did not seek to directly address this\nquestion. We apologize if our aims were not clear in our earlier response and have revised our\nmanuscript to clarify our goals accordingly. See below.\nImportantly, performing additional quantitative analyses to contrast local-global variability would\nsignificantly blow out an already extensive and comprehensive assessment of variability and fall\noutside the scope of this paper. We are actively pursuing this investigation in a separate project.\nLastly, we respectfully disagree with the reviewer that failure to quantitively address their request\nwould undermine the impact of our work. As outlined in our manuscript, the field of BOLD signal\nvariability remains significantly detached from its neurobiological and statistical underpinnings. Our\nwork directly addresses this gap and provides a conceptual foundation for future studies to build upon.\nRevised Results (note that the abstract, revised introduction and revised discussion already\nstate that local and global variability are treated independently):\n[…] We first sought to identify robust regional metrics of local and global BOLD signal variability.\nBy local and global BOLD signal variability, we refer to changes in the dynamic properties of\nsingle timeseries (local) and of pairs of timeseries (global). The former case is what is traditionally\nreferred to as “BOLD signal variability” in the literature, and the latter as “dynamic functional\nconnectivity”. In using the general term “variability” for both measures, we aimed at conceptually\njuxtaposing two lines of research that have traditionally been treated independently but recently\nshown to be inherently intertwined18,44. In doing so, our work establishes the foundation for future\nstudies to quantitatively compare them. (p. 5)\nREVIEWER 1, COMMENT 2: Patterns in the sample with older adults. This question has not been\naddressed. I do not think it is unreasonable to repeat the same analyses in old samples.\nWe agree with the reviewer that investigations into aging and lifespan development are of particular\ninterest. Yet, they are out of scope of the present manuscript, which focuses primarily on younger\nadults and already entails a significant number of analyses. A proper treatment of aging would\nrequire an extensive set of additional analyses, as aging affects both neural signal properties and the\nunderlying neurobiology. In sum, repeating our analyses in an aging cohort would result in a\nseparate line of inquiry rather than a direct extension of this study and ensuing findings would\nwarrant substantial interpretation and discussion that extend well beyond the word limit of this\nmanuscript.\nREVIEWER 1, COMMENT 3: Regarding my last comment on the relationship and comparisons\nwith other measures, I am aligned with Reviewer 2 R2 that it is important to demonstrate why this\nmetric is better than others and their implications in understanding our brain development, aging\nand behavior associations.\nWe thank the reviewer for rehearsing this point. Reviewer 2 was satisfied with how we addressed\nthis specific comment. As discussed in our previous response and in the manuscript, different\nmeasures of signal variability and complexity capture conceptually different timeseries properties.\nWe recognize the importance of benchmarking these approaches, yet, doing so would significantly\nshift the focus of the current work and compromise its clarity and coherence. As noted previously,\nwe are actively investigating this in a separate, dedicated project.\nReviewer #2 (Remarks to the Author):\nThe authors have addressed all my concerns.\nWe thank the reviewer for their constructive comments."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the hard work from the authors. The revision has addressed some concerns, but there are still three main concerns left:\n1. What are the unique features of specific low-level and high-level brain regions in the context of local and global BOLD signal variability? The authors did not perform any additional analyses to compare them quantitatively. I think this is important even though they are derived independently. The newly added discussion is also non-specific and not satisfying. E.g. “cortical areas devoted to local processing presented variability patterns that were more dependent on fMRI data type for both local and global variability measures, than higher-order regions.” What are the cortical areas devoted to local processing? How is it different from higher-order regions? Given the rich literature on functional gradients and organization involving sensorimotor and association regions, more comparison and interpretation for local and global variability is needed. To me, this response contradicts the authors’ claim in their earlier response that these two processes are linked, and in this study they purposely named them differently from the literature to study both (single time-series variability and paired time series (connectivity) variability. Without fully addressing this point, providing additional insights for the field is hard.\n2. Patterns in the sample with older adults. This question has not been addressed. I do not think it is unreasonable to repeat the same analyses in old samples.\n3. Regarding my last comment on the relationship and comparisons with other measures, I am aligned with Reviewer 2 R2 that it is important to demonstrate why this metric is better than others and their implications in understanding our brain development, aging and behavior associations.", "reviewer_2": "The authors have addressed all my concerns."}, "author_response": "We thank the reviewer for their constructive feedback. We agree that quantitively comparing low-level and high-level brain regions in the context of local-global signal variability metrics is an important direction for future work. Yet, the present study did not seek to directly address this question. We apologize if our aims were not clear in our earlier response and have revised our manuscript to clarify our goals accordingly. See below.\nImportantly, performing additional quantitative analyses to contrast local-global variability would significantly blow out an already extensive and comprehensive assessment of variability and fall outside the scope of this paper. We are actively pursuing this investigation in a separate project.\nLastly, we respectfully disagree with the reviewer that failure to quantitively address their request would undermine the impact of our work. As outlined in our manuscript, the field of BOLD signal variability remains significantly detached from its neurobiological and statistical underpinnings. Our work directly addresses this gap and provides a conceptual foundation for future studies to build upon.\nRevised Results (note that the abstract, revised introduction and revised discussion already state that local and global variability are treated independently):\n[…] We first sought to identify robust regional metrics of local and global BOLD signal variability. By local and global BOLD signal variability, we refer to changes in the dynamic properties of single timeseries (local) and of pairs of timeseries (global). The former case is what is traditionally referred to as “BOLD signal variability” in the literature, and the latter as “dynamic functional connectivity”. In using the general term “variability” for both measures, we aimed at conceptually juxtaposing two lines of research that have traditionally been treated independently but recently shown to be inherently intertwined18,44. In doing so, our work establishes the foundation for future studies to quantitatively compare them. (p. 5)\nREVIEWER 1, COMMENT 2: Patterns in the sample with older adults. This question has not been addressed. I do not think it is unreasonable to repeat the same analyses in old samples.\nWe agree with the reviewer that investigations into aging and lifespan development are of particular interest. Yet, they are out of scope of the present manuscript, which focuses primarily on younger adults and already entails a significant number of analyses. A proper treatment of aging would require an extensive set of additional analyses, as aging affects both neural signal properties and the underlying neurobiology. In sum, repeating our analyses in an aging cohort would result in a separate line of inquiry rather than a direct extension of this study and ensuing findings would warrant substantial interpretation and discussion that extend well beyond the word limit of this manuscript.\nREVIEWER 1, COMMENT 3: Regarding my last comment on the relationship and comparisons with other measures, I am aligned with Reviewer 2 R2 that it is important to demonstrate why this metric is better than others and their implications in understanding our brain development, aging and behavior associations.\nWe thank the reviewer for rehearsing this point. Reviewer 2 was satisfied with how we addressed this specific comment. As discussed in our previous response and in the manuscript, different measures of signal variability and complexity capture conceptually different timeseries properties. We recognize the importance of benchmarking these approaches, yet, doing so would significantly shift the focus of the current work and compromise its clarity and coherence. As noted previously, we are actively investigating this in a separate, dedicated project.\nReviewer #2 (Remarks to the Author):\nThe authors have addressed all my concerns.\nWe thank the reviewer for their constructive comments."}]

Decision letter