The coupling of BOLD signal variability and degree centrality underlies cognitive functions and psychiatric diseases

Brain signal variability has been consistently linked to functional integration; however, whether this coupling is associated with cognitive functions and/or psychiatric diseases has not been clarified. Using multiple multimodality datasets, including resting-state functional magnetic resonance imaging (rsfMRI) data from the Human Connectome Project (HCP: N = 927) and a Beijing sample (N = 416) and cerebral blood flow (CBF) and rsfMRI data from a Hangzhou sample (N = 29), we found that, compared with the existing variability measure (i.e., SDBOLD), the mean-scaled (standardized) fractional standard deviation of the BOLD signal (mfSDBOLD) maintained very high test-retest reliability, showed greater cross-site reliability and was less affected by head motion. We also found strong reproducible couplings between the mfSDBOLD and functional integration measured by the degree centrality (DC), both cross-voxel and cross-subject, which were robust to scanning and preprocessing parameters. Moreover, both mfSDBOLD and DC were correlated with CBF, suggesting a common physiological basis for both measures. Critically, the degree of coupling between mfSDBOLD and long-range DC was positively correlated with individuals' cognitive total composite scores. Brain regions with greater mismatches between mfSDBOLD and long-range DC were more vulnerable to brain diseases. Our results suggest that BOLD signal variability could serve as a meaningful index of local function that underlies functional integration in the human brain and that a strong coupling between BOLD signal variability and functional integration may serve as a hallmark of balanced brain networks that are associated with optimal brain functions.


Introduction
Cognitive functions depend on an advanced functional system that requires a dynamic, context-sensitive balance between functional integration and segregation in the brain ( Chan et al., 2017 ;Fornito et al., 2012 ;Shine et al., 2018 ). Consistently, it has been revealed that anatomically separated brain regions fluctuate synchronously and show strong functional connectivity (FC), thus forming a complex functional network ( Biswal et al., 1995 ). Among this network, some regions, such as the posterior cingulate, lateral temporal, lateral parietal, and medial/lateral prefrontal cortices, have disproportionately higher numbers of connections and act as brain hubs to integrate diverse informational sources and balance functional integration and segregation in the brain ( Buckner et al., 2009 ;Sporns et al., 2007 ). Given the impor-dynamic changes (or variability) of local neuronal ensemble activities, which reflect the neural flexibility ( Nomi et al., 2017 ), or the dynamic range that affects the adaptability and efficiency of neural systems ( Garrett et al., 2013b ). First, developmental studies have revealed that both FC ( Cao et al., 2017 ;Toulmin et al., 2015 ) and brain signal variability ( McIntosh et al., 2008 ;Vakorin et al., 2011 ) increase with age, although elderly subjects show reduced FC ( Damoiseaux, 2017 ;Ferreira and Busatto, 2013 ) and signal variability ( Armbruster-Genc et al., 2016 ;Garrett et al., 2010Garrett et al., , 2011Garrett et al., , 2013a. Second, individuals with stroke ( Kielar et al., 2016 ), multiple sclerosis ( Petracca et al., 2017 ), Alzheimer's disease ( Scarapicchia et al., 2018 ) and other neurological disorders ( Zoller et al., 2017 ) have regional neural variability abnormalities and accompanying brain network dysfunction ( Colasanti et al., 2016 ;Dennis and Thompson, 2014 ). Third, dopamine and other neurotransmitters could boost signal variability ( Alavash et al., 2018 ;Garrett et al., 2015 ) and increase FC ( Alavash et al., 2018 ;Kelly et al., 2009 ;Nagano-Saito et al., 2008 ).
Several studies have directly examined the coupling of signal variability and FC using different brain image data and measures of signal variability. For example, one electroencephalography (EEG) study used multiscale entropy (MSE) to measure signal variability and revealed positive correlations between MSE and FC ( Misic et al., 2011 ). Using the blood oxygen level-dependent (BOLD) signal and amplitude of low-frequency fluctuation (ALFF) or fractional ALFF (fALFF) to measure brain signal variability, several studies revealed that fALFF/ALFF (f/ALFF) was positively correlated with FC ( Di et al., 2013 ;Sato et al., 2019 ;Tomasi et al., 2016 ;Yan et al., 2017 ). The variability of BOLD signals could be measured by the resting state fluctuation amplitude (RSFA) ( Kannurpatti and Biswal, 2008 ) or temporal standard deviation (SD BOLD ) ( Garrett et al., 2010 ;Nomi et al., 2017 ), and a recent study observed that greater SD BOLD values were associated with lower functional network dimensionality within functional networks, suggesting that local variability in BOLD signals is related to information communication ( Garrett et al., 2018 ).
Previous studies have led to several important questions that deserve further investigation. First, most of the above studies either examined coupling across subjects or at the group level, and very few studies have examined coupling across voxels (but see Yan et al., 2017 ). Second, considering the replication crisis in neuroimaging literature, it is important to use a large dataset and to examine the reliability of this coupling across both scanning sessions and different sites. Finally, and most importantly, previous studies have not examined the individual differences in the coupling between BOLD signal variability and FC; thus, the functional relevance of this coupling is unknown.
The current study aimed to address these questions using large resting-state functional magnetic resonance imaging (rsfMRI) samples from different scanners and acquisition parameters. It should be noted that most variability measures, such as SD BOLD , ALFF or RSFA, are affected by factors that include the echo time (TE) and the initial BOLD signal strength (S 0 , when the TE = 0) , which increases the difficulty of comparing results from different scanners and acquisition parameters. To improve the cross-scanner reliability, we used a new brain signal variability measure, i.e., the mean-scaled fractional standard deviation of BOLD signal (mfSD BOLD ), which was based on a previous study ( Zou et al., 2008 ). As a result, the first part of the current study was to compare the stability and reliability of measures of BOLD signal variability and explored the relationships among the BOLD signal variability, degree centrality (DC, a measure of FC strength) and CBF.
In the second part, we further examined the functional relevance of the coupling between BOLD signal variability and DC, which were both positively correlated with CBF. First, we systematically examined the reliability of cross-voxel and cross-subject couplings between mfSD BOLD and DC across different scanning sessions, datasets, and processing parameters. Second, we investigated the functional relevance of this coupling by exploring the association between coupling strength and individuals' cognitive function. Third, we further examined whether the degree of decoupling between mfSD BOLD and DC was associated with disease vulnerability. We hypothesized that the variability of BOLD signals and degree centrality were tightly coupled, and that the degree of coupling was related to cognitive function and susceptibility to brain diseases.

Materials and Methods
Three datasets were used in this study. The primary data were obtained from the publicly available Human Connectome Project (HCP) ( N = 927). The other two validation datasets were acquired from Beijing Normal University ( N = 416) and Hangzhou Normal University ( N = 29). All participants were healthy adults. These three datasets are described in detail below.

HCP dataset: participants and data preprocessing
We used rsfMRI data from the 1200 Subjects Release (S1200) dataset of the HCP . A total of 927 subjects (434 males; age = 22-35 years old, mean age = 28.60 years old) were used for analyzes after excluding participants who did not finish both sessions (REST1 and REST2) ( N = 176), those with incomplete fMRI data (the number of volumes was smaller than 1200) ( N = 21), those with excessive head motion (mean framewise displacement, mFD > 0.3, N = 51) or those without cognitive total composite scores ( N = 11). We also excluded participants aged 36 years old and above ( N = 13) to increase the age homogeneity. To reduce the potential influence of encoding directions, we only used the two rsfMRI scans with a left-to-right (LR) phase-encoding direction. To rule out the influence of twins and siblings ( McDonough and Siegel, 2018 ), a subgroup of 356 unrelated participants (170 males; age = 22-35 years old, mean age = 28.42 years old) were randomly selected from the 927 participants to verify the results. Each rsfMRI run was acquired at 2 mm isotropic resolution with a repetition time of 720 ms and lasted for 14 min and 33 s. For details of the data acquisition parameters, see Smith et al. (2013) .
The extensive preprocessing pipeline was performed by the HCP team , which includes echo-planar imaging gradient distortion correction, motion correction, field bias correction, spatial transformation into common Montreal Neurological Institute (MNI) space , and artefact removal using independent component analysis (ICA) with the MELODIC and FIX tools implemented in FSL Salimi-Khorshidi et al., 2014 ). We additionally constructed a single regression model to simultaneously perform nuisance regression (including the mean white matter (WM), cerebrospinal fluid (CSF), and global signal time series) and bandpass (0.01-0.1 Hz)/high-pass ( > 0.01 Hz) filtering ( Lindquist et al., 2019 ) using the 3dTproject function ( Caballero-Gaudes and Reynolds, 2017 ) provided by Analysis of Functional NeuroImage (AFNI) ( Cox, 1996 ). To examine the effect of global signal regression on the results, we performed another analysis using the model described above but without global signal regression.

Beijing dataset: participants, data acquisition and preprocessing
A total of 493 healthy Chinese college students were recruited as a part of the Cognitive Neurogenetic Study of Chinese Young Adults (CN-SCYA) Project ( Feng et al., 2020 ). All of them had normal or correctedto-normal vision and no history of psychiatric or neurological diseases. Written consent was obtained from each participant after he/she received a full explanation of the study procedure. Seventy-seven participants were excluded due to excessive head motion (mFD > 0.3) ( N = 43) or incomplete cerebellar coverage ( N = 34). The final sample consisted of 416 participants (186 males; age = 17-29 years; mean age = 21.5 years). The study was approved by the Institutional Review Board at Beijing Normal University.
MRI data acquisition: Imaging data were acquired using a 3.0 T Siemens MRI scanner (Siemens Medical Systems, Erlangen, Germany) in the Brain Imaging Center at Beijing Normal University. Restingstate fMRI images were acquired using the Gradient Echo Type Echo-Planar Imaging (GRE-EPI) sequence. The acquisition parameters were as follows: TR/TE = 2000 ms/30 ms, flip angle (FA) = 90°, resolution matrix = 64 × 64, FOV = 200 × 200 mm 2 , thickness = 3.5 mm, and acquisition voxel size = 3.1 × 3.1 × 3.5 mm 3 . A total of 33 transverse slices were used to cover the whole brain. A total of 200 volumes were collected in a 6 min and 40 s scan, during which all participants were instructed to relax and keep their eyes closed but not to sleep ( Damoiseaux et al., 2006 ). The anatomical scan was acquired using the T1-weighted magnetization prepared rapid acquisition gradient-echo (MPRAGE) sequence with the following parameters: inversion time (TI) = 1100 ms; TR/TE/FA = 2530 ms/3.39 ms/7°, FOV = 256 × 256 mm 2 , matrix = 192 × 256, slice thickness = 1.33 mm, 144 sagittal slices, and voxel size = 1.3 × 1.0 × 1.3 mm 3 .
The first 5 volumes of each participant were removed to allow for T1 equilibration effects. The remaining 195 volumes were then corrected for intervolume head motion (six-parameter rigid-body transformation), intravolume temporal offsets (sinc interpolation) and high-pass filtered ( > 0.01 Hz). Data were then smoothed using a 6-mm Gaussian kernel. Subsequently, we implemented ICA denoising using the same approach as that used in the HCP dataset. Specifically, the classifier was trained on a small training set (20 participants). We first manually classified the "bad " components according to several primary criteria proposed by Griffanti et al. (2017 ). The training set was then used in a quadratic classifier to automatically separate components into artefact and nonartefact categories and applied to the new dataset. The identified artefact components were then filtered out from the corresponding fMRI data. All resulting images were spatially normalized to standard MNI space using DARTEL implemented in SPM and resampled to 3 × 3 × 3 mm 3 resolution. Finally, we performed the same nuisance regression and bandpass filtering that was performed on the HCP dataset.

Hangzhou dataset: participants, data acquisition and preprocessing
To examine the relationship among BOLD signal variability, CBF and FC, we analyzed another Hangzhou dataset with resting-state BOLD fMRI images and arterial spin labeling (ASL) fMRI images ( Sheng et al., 2018 ). Twenty-nine healthy right-handed participants were included in the current study (5 males; age = 23-66 years; mean age = 42.31 years). This study was approved by the Ethics Committee of the Sir Run Run Shaw Hospital, School of Medicine, Zhejiang University, and the Affiliated Hospital of Hangzhou Normal University. All participants provided written informed consent.
MRI data processing. The rsfMRI data were processed with the same procedures as those used on the Beijing dataset. For the ASL MRI data, individual CBF images were first obtained using Functool (version 12.2.01), an automated image postprocessing tool embedded in the GE healthcare MR-750 system. To correct for the partial volume effects (PVEs) introduced by the large voxel size used in ASL MRI, participants' structural images were segmented into grey matter (GM), WM, and CSF using SPM12. The corresponding GM, WM, and CSF probability maps were then registered and resampled to the native ASL MRI space using an inter-modality registration algorithm that uses normalized mutual information as the objective function. The PET-based method was then used to correct PVEs ( Hu et al., 2010 ;Wang et al., 2008 ). Finally, the PVE-corrected CBF images were spatially normalized to standard MNI space (using the transformation fields derived from the tissue segmentation data of structural images), resampled to 3 mm isotropic voxels, and spatially smoothed (Gaussian kernel with full width at half maximum (FWHM) = 6 mm).

Calculation of DC, SD BOLD and mfSD BOLD
After the above preprocessing steps, we calculated the voxelwise full/long-range DC, SD BOLD and fSD BOLD /mfSD BOLD (m/fSD BOLD ).
The DC quantifies the overall connectivity of a node to all other nodes in a network. Briefly, for a given voxel, its time series were extracted and correlated with the time series of all other voxels in the brain. To improve the normality of the correlations, we transformed the correlation coefficients to Z-scores using Fisher's r-to-Z transformation. The summation of the resulting Fisher's Z-score across all other voxels was then calculated and assigned to this given voxel. To exclude the confounding effects of spurious correlations, correlation coefficients smaller than a cut-off threshold r were excluded prior to the summation. To make the network sparsity level comparable across different datasets, we used different cut-off correlation thresholds to calculate DC for different datasets (Table S1). Negative correlations were excluded due to their ambiguous interpretation ( Fox et al., 2009 ;Murphy et al., 2009 ). For full-range DC, correlations with spatially close voxels (Euclidean distance < 20 mm) were excluded ( Power et al., 2011 ). These correlations were excluded to avoid the influence of possible shared signals between nearby voxels on DC and to minimize the effects of data processing (e.g., blurring, reslicing) and head motion on DC ( Power et al., 2012 ). Additionally, we calculated the long-range DC, which only included the correlations with voxels that had an anatomical distance of more than 75 mm ( Achard et al., 2006 ;He et al., 2007 ) with a given voxel, allowing us to further explore the effects of anatomical distance in this study. Furthermore, we examined the impact of the global signal (without global signal regression, NGSR) and the network type (binary network, defined by an adjacent matrix using one (the edges that survived the abovementioned cut-off threshold r ) and zero (edges below threshold)) ( Bullmore and Sporns, 2009 ) on the results using the HCP dataset.
The SD BOLD for each voxel was quantified as the standard deviation across the time series of low-frequency (0.01 < f < 0.1 Hz) signals ( Garrett et al., 2018 ), which reflected the extent of the BOLD signal changes over time. To minimize the effects of the initial BOLD signal strength (S 0 , which is principally altered by structural noise sources) ( Power et al., 2018 ) and TE , following Zou et al. (2008) , we also calculated the fSD BOLD , which was quantified as the low-frequency (0.01-0.1 Hz) SD BOLD divided by the high-pass filtering frequency ( > 0.01 Hz) SD BOLD . Furthermore, we calculated the mean-scaled (standardized) fSD BOLD (mfSD BOLD ) to minimize the effects of in-scanner head movements and global effects by dividing the fSD BOLD of each voxel by the mean fSD BOLD of the whole brain ( Kublbock et al., 2014 ).
The resultant DC, SD BOLD and m/fSD BOLD maps were further spatially smoothed with a 5 mm Gaussian kernel for the HCP dataset only because the Beijing and Hangzhou data were already smoothed during preprocessing. Unless otherwise specified, the results of the HCP data were obtained based on the average of the two sessions (REST1 and REST2). Notably, all the DC, SD BOLD and m/fSD BOLD calculations were restricted within dataset-specific brain masks derived according to the following two criteria: (a) nonzero variance of BOLD signals for all participants and (b) > 20% grey matter tissue probability based on the tissue probability map attained using the SPM12 package.

Stability and reliability of BOLD signal variability
To examine the stability and reliability of the variability of BOLD signals, we calculated the Pearson correlation coefficient and intraclass correlation coefficient (ICC) across sessions (REST1 and REST2 for the HCP dataset), different datasets (HCP, Beijing, and Hangzhou), different cut-off thresholds (0.05, 0.075, 0.10 and 0.15) for DC calculations (HCP dataset), and different data processing strategies, such as network binarization and without global signal regression (HCP dataset). Using psych and lme4 packages in R 4.0.2, the ICC was calculated by fitting linear mixed-effects models based on the following formula ( Shrout and Fleiss, 1979 ): where MSB is the between-target mean square and MSE is the mean square of errors. An ICC value greater than 0.4 is commonly expected in practice ( Zuo and Xing, 2014 ). Due to the different spatial resolutions and brain coverage across datasets, the spatial pattern similarity of the three datasets could not be examined at the voxel level. Thus, the Fan 246 template ( Fan et al., 2016 ) was used to extract the signals in each brain region, which were then used to calculate the cross-site reliability. We retained 222 regions of interest (ROIs) that were shared by all participants in all three datasets. In addition, the Yeo 2011 template ( Yeo et al., 2011 ) with 114 parcels (108 remaining parcels after combining all three datasets) was used to verify this result.
In a further analysis, we examined the cross-site reliability/ICC within each brain region. First, the group-averaged BOLD signal variability maps of the HCP dataset were first resampled to 3 mm resolution. Then, we calculated their cross-site reliability at each brain region based on the Fan 246 template ( Fan et al., 2016 ). Since the number of voxels was different in each brain region, we restricted our analysis to regions with at least 50 voxels (170 ROIs). For the regions with more than 50 voxels, 50 voxels were randomly selected from each region to calculate the ICC, and this process was repeated 1000 times. The crosssite reliability for each brain region was obtained by averaging the ICCs from the 1000 repetitions.

Correlations among mfSD BOLD , DC and CBF
Pearson correlation analyzes were performed to examine the crosssubject or cross-voxel relationship among mfSD BOLD , DC and CBF. Neighboring voxels were spatially dependent mainly because of spatial smoothing ( Xiong et al., 1995 ); thus, the effective degree of freedom ( df eff ) was corrected by the following formula ( Liang et al., 2013 ): where v is the nominal volume of the voxels, N is the number of voxels used in the analysis, and FWHM represents the spatial smoothness estimated by using AFNI's 3dFWHMx. The df eff result was then used to estimate the p values for the cross-voxel correlation analyzes. For the cross-subject correlation analyzes, we controlled for age, sex, age 2 , age × sex, age 2 × sex and mFD.

Functional relevance of the coupling between mfSD BOLD and DC
The coupling between DC and mfSD BOLD was defined based on their correlation coefficient, with a higher correlation indicating tighter coupling. Before examining the functional relevance of DC-mfSD BOLD coupling, the individual correlation coefficients were transformed to Fisher's Z-scores to improve the normality of the correlations and to facilitate further statistical tests.
Age and sex differences in DC-mfSD BOLD coupling . To examine the age (22-35 years) and sex (male and female) effects of the coupling, we performed linear least-squares regression in R 4.0.2, with age, sex, and age by sex interaction as independent variables and mFD as a covariant.
DC-mfSD BOLD coupling and cognitive performance . We used the cognitive total composite scores (CogTotalComp) in the HCP dataset to quantify individuals' cognitive performance . The scores were derived by averaging the normalized scores of each of the fluid and crystallized cognition measures, which were then converted to standardized scores using the following formula: − ̄ × 15 + 100 , where ̄ is the sample mean of the averaged normalized scores and is the sample standard deviation of the averaged normalized scores. Partial correlation analyzes were performed to examine the relationship between DC-mfSD BOLD coupling (correlation coefficient) and cognitive total composite scores, with age, sex, age 2 , age × sex, age 2 × sex and mFD as control variables. Permutation tests ( N = 10,000) were performed to test the significance of the R values.
DC-mfSD BOLD mismatch and cognitive performance . The DC and mfSD BOLD values were first separately Z-transformed, and a leastsquares linear regression model was then applied to fit the DC and mfSD BOLD results across voxels for each participant. Residuals that could not be explained by mfSD BOLD were used as a measurement of DC-mfSD BOLD mismatch, with a higher absolute value of the DC residual reflecting a greater mismatch between DC and mfSD BOLD . Subsequently, for a given voxel, a partial correlation analysis was conducted between the individuals' absolute DC residuals and the cognitive total composite scores, controlling age, sex, age 2 , age × sex, age 2 × sex and mFD, and this analysis was performed in DPABI ( Yan et al., 2016 ). The permutation test ( N = 1000) with threshold-free cluster enhancement (TFCE) was used to correct whole-brain multiple comparisons, which achieved the best balance between the familywise error rate ( < 5%) and test-retest (TRT) reliability Winkler et al., 2016 ).
DC-mfSD BOLD mismatch and disease susceptibility. To explore whether brain regions exhibiting greater mismatch were also more vulnerable to brain diseases, we adopted the disease-susceptible areas from a previous study ( Sha et al., 2018 ). Notably, there were two disordergeneral activation likelihood estimation (ALE) maps, with one indicating the probability of reduced activity (i.e., disorders < controls; Fig. 6 A, left) and the other indicating the probability of increased activity (i.e., disorders > controls; Fig. 6 A, middle). Since some brain regions can show both increased or decreased activity, we generated a combined disease-susceptible map by averaging the two maps and intersecting them with a dataset-specific mask. We performed two analyzes to examine disease susceptibility and DC-mfSD BOLD mismatches. First, we conducted a paired t-test to compare the strength of the DC-mfSD BOLD mismatch (i.e., higher absolute DC residuals) between the brain regions that are susceptible and insusceptible to brain diseases. We further examined whether brain regions that are more susceptible to brain diseases showed higher absolute DC residuals within the disease-susceptible areas. To test this hypothesis, we grouped the disease-susceptibility voxels into five bins according to the ALE score in increments of 0.2 (0:0.2:1) and then averaged the long-range and full-range DC residuals within each bin. We then used repeated-measures ANOVAs to examine whether there were significant main effects of disease susceptibility on absolute DC residuals. Notably, unless otherwise specified, p values were adjusted using the Holm-Bonferroni correction ( Holm, 1979 ) for multiple comparisons.  (upper panel) and SD BOLD (lower panel) for three datasets. (B) Spatial similarity across sites for SD BOLD (upper triangle) and mfSD BOLD (lower triangle) based on the Fan 246 template. (C) Cross-site reliability of mfSD BOLD (left) and SD BOLD (right) for each brain region. (D) Correlations between signal variability and mFD (each dot represents one participant). (E) Reliability of mfSD BOLD cross two sessions (HCP: REST1 and REST2; orange) and different preprocessing parameters (HCP: with and without global signal regression; green) for each individual participant. SD BOLD , standard deviation of the BOLD signal; fSD BOLD , fractional SD BOLD ; mfSD BOLD , meanscaled fSD BOLD ; ICC, intraclass correlation coefficient; mFD, mean framewise displacement (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

mfSD BOLD is a reliable measure of BOLD signal variability
SD BOLD has been used to quantify the variability of the BOLD signal for data from individual scanners. Nevertheless, the intra-site and intersite reliabilities of BOLD signal variability have not been strictly examined. In particular, the overall signal variability could be affected by the noise level (such as the physiological noise and head motion) and the signal-to-noise ratio of the scanner, making it less reliable across sites. We proposed that the mfSD BOLD , which is the mean-scaled proportion of the low-frequency SD BOLD (0.01 < f < 0.1 Hz) related to the overall SD BOLD ( > 0.01 Hz), could reduce the effect of noise and head motion. We performed several analyzes to validate this measure.
First, we found that mfSD BOLD was highly consistent across datasets ( Fig. 1 A, upper panel) while SD BOLD was not ( Fig. 1 A, lower panel). In particular, a high and consistent cross-site correlation was observed for mfSD BOLD ( rs = 0.62, 0.53 and 0.94, ps < 10 − 10 ; ICC = 0.792 ) while negative correlations were found between the HCP data and the other two sites ( rs = -0.42 and -0.46, ps < 10 − 10 ; ICC < 0.001 ) for SD BOLD ( Fig. 1 B) based on the Fan 246 template. Direct comparisons revealed higher cross-site pattern similarity for mfSD BOLD than for SD BOLD ( Zs = 12.27,11.38 and 5.41, ps < 10 − 7 ). Similar results were found based on the Yeo 2011 template (Fig. S1A). Notably, although the spatial pattern of mfSD BOLD showed greater cross-scanner consistency than SD BOLD , we found that the spatial pattern similarity of mfSD BOLD between the HCP and Beijing ( r = 0.62 ) or Hangzhou ( r = 0.53 ) datasets was still lower than that between the Beijing and Hangzhou ( r = 0.94 ) datasets ( Fig. 1 B), which was probably due to the differences in scanner and scanning parameters.
Across different brain regions, we found that mfSD BOLD has high cross-site reliability ( ICC > 0.4 ) in the cingulate gyrus, inferior frontal gyrus, inferior parietal lobules and fusiform gyrus ( Fig. 1 C, left), while SD BOLD has low cross-site reliability ( ICC ≤ 0.2 ) in many cortical regions except the cingulate gyrus, medial frontal gyrus, and part of the middle temporal gyrus ( Fig. 1 C, right). In general, the cross-site reliability at different brain regions of mfSD BOLD was higher than that of SD BOLD (paired t = 12.23, p < 0.001 ).
Furthermore, we found that the spatial pattern between mfSD BOLD and SD BOLD was quite similar for the Beijing ( r = 0.61 ± 0.05 ; Fig. S1B, left) and Hangzhou ( r = 0.66 ± 0.08 ; Fig. S1B, left) datasets, whereas that for the HCP dataset was much lower ( r = 0.27 ± 0.06 ; Fig. S1B, left). This finding is largely because the spatial pattern of SD BOLD is highly variable across datasets while that for mfSD BOLD is more consistent ( Fig. 1 B).
Second, a previous study suggested that head motion could affect overall BOLD signal variability ( Kublbock et al., 2014 ). Consistently, using the HCP dataset, we found that SD BOLD ( r = 0.139, p < 0.001 ) and fSD BOLD ( r = -0.230, p < 10 − 10 ) were correlated with individuals' head motion (measured by mFD) ( Fig. 1 D), whereas mfSD BOLD was not ( r = 0.054, p = 0.101 ) ( Fig. 1 D). This result was replicated on the Beijing sample ( Table 1 ), which SD BOLD ( r = 0.19, p < 0.001 ) was affected by head motion, while mfSD BOLD was not ( r = 0.07, p = 0.164 ). Nevertheless, similar results had not been observed on the Hangzhou sample ( Table 1 ), probably due to the small sample size. It should be noted that mean scaling did not affect the spatial pattern and thus did not affect the cross-site reliability. Third, both mfSD BOLD ( Fig. 1 E, orange) and SD BOLD ( Table 1 ) showed remarkably high spatial pattern similarities across sessions from the HCP for each individual (all ICCs > 0.90 ). Finally, the spatial distribution of neither mfSD BOLD (all ICCs > 0.90 ; Fig. 1  As summarized in Table 1 , compared to SD BOLD , the spatial pattern of mfSD BOLD was more reliable across sites, although both measures were stable across different sessions and robust to different preprocessing parameters. In addition, the mfSD BOLD was less affected by head movement as compared to SD BOLD . Accordingly, we used mfSD BOLD as the measure of BOLD signal variability in the following primary analyzes.

Relationships among DC, mfSD BOLD and CBF
To understand the physiological basis of DC and mfSD BOLD , we also measured CBF in the Hangzhou sample. Consistent with previous observations ( Liang et al., 2013 ;Sheng et al., 2018 ), we found a significant positive correlation between DC and CBF at the group level ( Fig. 2 A; full-range DC & CBF: rs = 0.480-0.496 for different thresholds; longrange DC & CBF: rs = 0.286-0.396 for different thresholds; ps adj < 10 − 10 ; df eff = 993 ) based on the mean DC and CBF maps (Fig. S2C), and this positive correlation also occurred for each individual ( Fig. 2  Intriguingly, mfSD BOLD values were also positively correlated with CBF both at the group level ( r = 0.462 , p < 10 − 10 , df eff = 993 ; Fig. 2 C)

Fig. 2. Relationships among DC, mfSD BOLD and CBF in the Hangzhou sample.
Correlation between DC and CBF at the (A) group level and (B) individual level under four correlation thresholds of DC calculations. The black dots represent mean correlation coefficients. Each colored dot represents one participant, and the grey areas represent + /-one standard deviation. The correlation between mfSD BOLD and CBF at the (C) group level (each dot represents one voxel in the brain) and (D) individual level. Group-averaged brain maps of DC and CBF; please see Fig. S2C. DC, degree centrality; mfSD BOLD , mean-scaled fractional standard deviation of the BOLD signal; CBF, cerebral blood flow; Corr., correlation. and for each individual (mean r = 0.296, ranging from 0.138 to 0.473 ; all ps < 10 − 5 , df eff = 993 ; Fig. 2 D). The above results suggest that DC and mfSD BOLD share a common physiological basis. Similarly, SD BOLD was also associated with CBF ( r = 0.39 ± 0.09 , all ps < 10 − 5 ) (Fig. S3), which is not surprising given that the spatial patterns of SD BOLD and mfSD BOLD for the Hangzhou dataset were remarkably similar ( Figs. 1 A,  S1B). Scattered distribution of the mfSD BOLD and fullrange (upper) or long-range DC (lower) at the group level. A correlation threshold of 0.10 was used for DC. Each dot represents one voxel in the brain. Cross-subject correlations between (D) full-range or (E) long-range DC (correlation threshold of 0.10) and mfSD BOLD . The cross-subject results using other DC correlation thresholds were similar and are shown in Fig.  S12. For the group-averaged DC maps, please see Fig. S2. DC, degree centrality; mfSD BOLD , mean-scaled fractional standard deviation of the BOLD signal; Corr., correlation.

Tight coupling between DC and mfSD BOLD
Having verified the reliability of mfSD BOLD and revealed a common physiological basis of DC and mfSD BOLD , we then examined the correlation between DC and mfSD BOLD . This correlation was examined across voxels for each participant and across participants for each voxel.
Cross-subject coupling between DC and mfSD BOLD . Having shown the significant cross-voxel coupling between DC and mfSD BOLD , we further examined whether there was significant cross-subject coupling, i.e., whether participants showing greater DC also exhibited greater mfSD BOLD . Indeed, we found strong and positive cross-subject couplings in most brain regions for both full-range ( Fig. 3 D) and long-range DC ( Fig. 3 E). The strongest positive full/long-range DC-mfSD BOLD couplings ( r > 0.65 ) were found in the bilateral frontal, occipital, and parietal lobe regions (Table S2). Interestingly, we also found some negative correlations between full/long-range DC and mfSD BOLD correlations, with the strongest negative correlations ( r < -0.55 ) located in the left an- Note: CogTotalComp, cognitive total composite scores; DC, degree centrality; mfSD BOLD , mean-scaled fractional standard deviation of the BOLD signal.
terior cerebellum lobe, bilateral hippocampus, and left putamen (Table S3). In an analysis of 356 unrelated participants, the results were also stable (Fig. S9). Again, the spatial pattern of cross-subject coupling was quite reliable across sessions (

Functional relevance of DC-mfSD BOLD coupling
The above analyzes revealed tight cross-voxel and cross-subject couplings between DC and mfSD BOLD . In addition to the overall strong coupling, we also found significant differences in coupling strength across participants and brain voxels. If DC-mfSD BOLD coupling reflects a fundamental organizing principle of the human brain, deviation from this coupling should have functional consequences. In particular, we predicted that participants with stronger coupling would show better cognitive performance and brain regions deviating from this coupling would be more susceptible to malfunctions. We performed several analyzes to test these intriguing hypotheses using the HCP dataset.

Coupling strength related to individual cognitive performance
Two analyzes were performed to examine the relationship between individual coupling strength and cognitive performance. In the first analysis, we examined the effects of age and sex on cognitive performance and DC-mfSD BOLD coupling strength through linear least-square regression models. We found that the cognitive total composite scores significantly decreased with age ( t = -2.81, p = 0.005 ; Table 2 ). The main effect of sex and the interaction between age and sex on the cognitive total composite scores were not significant ( Table 2 ; ps > 0.05 ). In parallel, we also found that the coupling between long-range DC and mfSD BOLD significantly decreased with age ( ts ≥ -2.94, ps adj < 0.05 ; Table 2 ), although the age effect was not consistent for full-range DC-mfSD BOLD coupling across different cut-off thresholds ( Table 2 ). Additionally, the main effects of sex and age by sex interaction were not significant ( Table 2 ; all ps > 0.05 ). Consistent results were found when only unrelated participants were used (Table S4).
In the second analysis, we further examined whether the strength of DC-mfSD BOLD coupling and cognitive performance was correlated across participants after controlling for age, sex, age 2 , age × sex, age 2 × sex and mFD. We found that the coupling strength between long-range DC and mfSD BOLD was positively correlated with the cognitive total composite scores across four correlation thresholds for DC ( rs = 0.09 ~0.11, ps adj < 0.05 ; Fig. 4 A; Table 3 ). No significant correlation was found for the coupling strength of full-range DC and mfSD BOLD ( rs = 0.02 ~0.04, ps unadj > 0.05 ; Fig. 4 A; Table 3 ). We also found comparable correlation coefficients when using unrelated participants (CogTotalComp & full-range DC-mfSD BOLD : rs = 0.07-0.08, ps unadj > 0.05 ; CogTotalComp & longrange DC-mfSD BOLD : rs = 0.09-0.11, ps unadj = 0.017-0.032 ), although the statistical significance was lower due to the reduced sample size. However, the coupling of DC (correlation threshold = 0.1) and SD BOLD was not associated with cognitive performance (CogTotalComp & fullrange DC-SD BOLD : r = -0.004, p = 0.909 ; CogTotalComp & long-range DC-mfSD BOLD : r = -0.006, p = 0.848 ; Table 3 ).
Previous studies have implicated short-range connections (ranging from 20 mm to 75 mm) in functional segregation in the human cortex Sporns, 2013 ). We also examined the cross-voxel coupling of short-range DC (correlation threshold = 0.1) and mfSD BOLD and their functional relevance based on the HCP dataset. We found that although short-range DC was also positively correlated with mfSD BOLD ( rs = 0.686 ± 0.017, ps < 0.001 ) for all participants, the degree of coupling did not correlate with cognition ( r = 0.001, p = 0.973 ), suggesting the specific functional role of long-range DC and mfSD BOLD coupling. Across participants, we did not find any brain region where the DC, mfSD BOLD or SD BOLD values were significantly correlated with cognitive performance.

Fig. 4. Coupling strength is associated with individual cognitive performance. (A)
Relationships between the coupling of DC and mfSD BOLD with cognitive total composite scores evaluated separately for each DC threshold. (B) Scatter plot of the strength of long-range DC (with a threshold of 0.10) and mfSD BOLD coupling and the cognitive total composite scores after controlling for confounding variables (each dot represents one participant). CogTotalComp, cognitive total composite scores; DC, degree centrality; mfSD BOLD , mean-scaled fractional standard deviation of the BOLD signal. * p adj < 0.05, * * p adj < 0.01, * * * p adj < 0.001.

Table 3
Relationships between cross-voxel coupling of DC and signal variability (SD BOLD vs mfSD BOLD ) and cognitive performance for the HCP data. Note: The correlation threshold of DC was 0.1. SD BOLD , standard deviation of BOLD signal; mfSD BOLD , meanscaled fractional standard deviation of BOLD signal; DC, degree centrality; rSTG, right superior temporal gyrus; rMCC, right middle cingulate cortex. * * * p adj < 0.001; n.s. non-significant; -no statistically significant brain region.

Residuals of DC-mfSD BOLD coupling related to cognitive performance
In addition to the overall coupling strength, we further examined the regional differences in DC-mfSD BOLD coupling and their relationships to brain functions. In the cross-voxel analysis, we constructed least-squares linear regression models for DC and mfSD BOLD . The regional mismatches between DC and mfSD BOLD were defined as the residual DC after regressing out mfSD BOLD ( Fig. 5 A). We found that DC residuals were very stable and reliable across sessions ( Overall, the reliability was higher for longrange DC residuals than for full-range DC residuals (Fig. S14).
We then examined whether the mismatch between DC and mfSD BOLD was related to individual cognitive performance. To address this question, we performed partial correlation analyzes for full-range/longrange DC residuals and cognitive function composite scores after controlling for age, sex, age 2 , age × sex interaction, age 2 × sex interac- Fig. 6. Brain regions showing mismatches between long-range DC and mfSD BOLD were susceptible to brain diseases. (A) Susceptible maps of brain diseases from Sha et al. (2018) . The ALE maps of higher (left) and lower (middle) disease-associated activities were averaged to generate the common disease-susceptible map (right). (B) Absolute long-range DC residuals were higher in regions that are susceptible to diseases than in regions that are not susceptible to the HCP (correlation threshold of DC = 0.1), Beijing and Hangzhou (DC correlation threshold of DC = 0.2) datasets. Error bars represent within-subject standard errors. (C) Absolute long-range DC residuals in brain regions show different degrees of disease susceptibility for the HCP dataset. The brain susceptibility regions were divided into 5 bins based on susceptibility (i.e., ALE score). Error bars represent 99% confidence interval within-subject standard errors. * * * p adj < 0.001. tion and mFD. We found that the absolute long-range DC residuals in the right superior temporal gyrus (STG) and right middle cingulate cortex (MCC) ( Fig. 5 B & C) were significantly and negatively correlated with cognitive performance under four different thresholds of DC. None of the brain regions showed a significant correlation between absolute full-range DC residuals and cognitive performance. When the same analyzes were performed on unrelated participants from the HCP dataset ( N = 356), we also found a negative correlation between the absolute long-range DC residuals and cognitive performance in the right MCC and right STG (Fig. S15) across four thresholds of DC, although the correlations were not statistically significant after multiple comparison corrections. In contrast, we did not find any brain region where the residuals of DC-SD BOLD coupling were significantly correlated with cognitive performance ( Table 3 ).
To further examine whether the direction of mismatches affects cognitive performance, we used the raw residual instead of the absolute residual. We found that the right STG (Fig. S16A, left panel) showed positive long-range DC residuals (Fig. S16A, middle panel), which were significantly negatively correlated with cognitive performance (Fig. S16A, right panel). In contrast, the right MCC (Fig. S16B, left panel) showed overall negative long-range DC residuals (Fig. S16B, middle panel), which were significantly positively correlated with cognitive performance (Fig. S16B, right panel). The results together suggest that participants with greater deviation (positive or negative) from the long-range DC-mfSD BOLD coupling showed worse cognitive function.

Mismatched brain regions are susceptible to brain diseases
The above analyzes suggest that a greater mismatch in mfSD BOLD and DC was related to worse cognitive permanence. A further question follows: were the regions showing a greater mismatch also more vulnerable to brain diseases? To answer this question, we defined the diseasesusceptible areas (for details, please see Materials and Methods ) from the results of a meta-analysis ( Sha et al., 2018 ) and conducted a paired ttest to explore whether the strength of the DC-mfSD BOLD mismatches of the disease-susceptible areas is significantly higher than that of insusceptible regions. The results revealed that for HCP sample, the absolute DC residuals within the disease-susceptible regions (number of voxels = 40323) were significantly larger than those outside of the diseasesusceptible regions (number of voxels = 65878) for both long-range DC residuals ( paired-ts = 137.278-207.098 for different thresholds, ps adj < 10 − 10 ) ( Fig. 6 B; Table 4 ) and full-range DC residuals ( paired-ts = 33.480-70.838 for different thresholds; ps adj < 10 − 10 ) ( Table 4 ).

Table 4
The differences in DC and signal variability (SD BOLD vs mfSD BOLD ) residuals between disease susceptible and insusceptible brain regions. Note: The DC correlation threshold of the HCP was 0.1, and that of the Beijing and Hangzhou samples was 0.2; SD BOLD , standard deviation of BOLD signal; mfSD BOLD , mean-scaled fractional standard deviation of BOLD signal; DC, degree centrality; HCP, Human Connectome Project. * * * p adj < 0.001; n.s. non-significant.
Using the Beijing and Hangzhou datasets, we also found the absolute long-range DC (correlation threshold = 0.2) residuals within the disease-susceptible regions (Beijing: number of voxels = 12279; Hangzhou: number of voxels = 13829) were significantly larger (Beijing: paired-t = 4.85, p < 0.001 ; Hangzhou: paired-t = 5.95, p < 0.001 ) than those outside of the disease-susceptible regions (Beijing: number of voxels = 24086; Hangzhou: number of voxels = 27219) ( Table 4 ; Fig. 6 B), but the absolute full-range DC residuals were smaller within the disease-susceptible regions than those outside of the disease-susceptible regions for the Beijing dataset ( paired-t = -11.66, p < 0.001), and no difference was found for the Hangzhou dataset ( paired-t = 1.55, p = 0.132 ) ( Table 4 ). These results suggest that the long-range DC and mfSD BOLD coupling is a more reliable indicator of brain function, consistent with the finding on the strength of coupling and cognitive performance.
We further examined whether DC-SD BOLD mismatches were also related to disease susceptibilities. The results showed that, for full-range DC, the residuals in the disease-susceptible area were smaller than those outside the disease-susceptible area for the HCP ( paired-t = -110.32, p < 0.001 ) and Beijing ( paired-t = -12.17, p < 0.001 ) datasets, but not for the Hangzhou dataset ( paired-t = -1.46, p = 0.156 ) ( Table 4 ; Fig. S17A). For long-range DC, compared to the residuals outside the disease-susceptible regions, those within the disease-susceptible regions were smaller for the HCP dataset ( paired-t = -141.06, p < 0.001 ), greater for the Beijing dataset ( paired-t = 5.86, p < 0.001 ), and not different for the Hangzhou dataset ( paired-t = 1.24, p = 0.225 ) ( Table 4 ; Fig. S17B).
In a further analysis, we examined whether the more susceptible brain regions showed higher DC-mfSD BOLD mismatches within the disease-susceptible areas. To test this hypothesis, we grouped the disease susceptibility voxels into five bins according to the ALE score in steps of 0.2 (0:0.2:1) and then averaged the long-range and full-range DC residuals within each bin. A repeated-measures ANOVA revealed a significant main effect of disease susceptibility on long-range DC-mfSD BOLD residuals at all correlation thresholds for the HCP dataset ( F = 4435.9 5 6439, ps < 10 − 10 ) ( Fig. 6 C and also see Table S5 for details of the post hoc tests). However, no such relationship was found for the Beijing and Hangzhou datasets or for the full-range DC-mfSD BOLD residuals. These results further suggest that brain regions deviating from long-range DC-mfSD BOLD coupling are more susceptible to mental diseases.

Discussion
In the current study, we examined how the coupling between BOLD signal variability and degree centrality is related to brain functions. We found strong and reproducible coupling between BOLD signal variability and degree centrality across sessions and scanning sites. Furthermore, the strength of coupling was related to cognitive functions, and the regions with decoupling were susceptible to brain diseases. Altogether, these findings reveal an important balance of brain signal variability and informative integration, which is essential for brain functions.

mfSD BOLD was stable and reliable
The first goal of the current study was to validate a BOLD signal variability measure that would be reliable across datasets and sessions and less affected by head motion artefacts. We found that the traditional BOLD signal variation measure, i.e., SD BOLD , showed quite different or even opposite spatial patterns across different datasets. For example, relatively stronger SD BOLD values were found in the thalamus, hippocampus, and anterior cingulate cortex in the HCP dataset, while weaker SD BOLD values were found in these structures than other regions in the Beijing and Hangzhou datasets. A change in BOLD contrast could be described as a change in the transverse relaxation rate or R 2 * (the inverse of T 2 * ) due to changes in blood oxygenation ( Bandettini et al., 1994 ;Menon et al., 1993 ;Ogawa et al., 1990 ), although it was also related to S 0 (when TE = 0) ( Kundu et al., 2012 ) and TE ( Bandettini et al., 1994 ;Evans et al., 2015 ). In the single-echo fMRI data, whether the change is due to a changed decay rate or a changed starting intensity cannot be determined. Therefore, it has been difficult to remove the influence of these effects ( Birn et al., 2006 ;Power et al., 2015 ), even when performing bandpass filtering. Importantly, by calculating the ratio of low-frequency signal fluctuations to full-frequency signal fluctuations, we could effectively suppress these effects on BOLD signal variability and generate a more reliable spatial pattern across sites, which is consistent with a previous study that measured BOLD variability with ALFF ( Zou et al., 2008 ). It should be noted that multi-echo fMRI could also control S 0 and TE to minimize the influence of these unrelated variables Kundu et al., 2012 ).
In addition, we demonstrated that the mean fSD BOLD of grey matter is correlated with head motion across subjects, which may be partially related to the introduction of a strong high-frequency signal with head movement. Since the fSD BOLD reflected the ratio of the low-frequency signal variation to the variation in the overall signal ( > 0.01 Hz), it was negatively correlated with head motion. A previous study based on fALFF found similar results ( Satterthwaite et al., 2012 ). Since the effect of head motion on neural signals is global, a normalization process, i.e., division by the mean fSD BOLD of the whole brain, could retain the spatial patterns but reduce the effect of head motion on overall signal variability ( Kublbock et al., 2014 ;Zou et al., 2008 ). Finally, we found that the mfSD BOLD showed remarkably high reliability and was unaffected by global signal regression. All these characteristics suggest that the mfSD BOLD could be an ideal index for future studies to examine BOLD signal variability.

Tight coupling between mfSD BOLD and DC
The current study revealed a consistent and robust coupling between DC and mfSD BOLD in all three datasets, both cross-voxels and cross-subjects. The reliability of this coupling was further proven by its robustness to processing parameters. The tight coupling was also consistent with previous studies that used different indicators of brain signal variability, such as f/ALFF ( Di et al., 2013 ;Sato et al., 2019 ;Tomasi et al., 2016 ;Yan et al., 2017 ), sample entropy ( Misic et al., 2011 ), and SD BOLD ( Garrett et al., 2018 ). Although f/SD BOLD and f/ALFF have been separately used in previous studies, our data suggest that they are strongly associated ( r > 0.9 ) (Fig. S18). Compared to ALFF ( Zang et al., 2007 ), SD BOLD may be easier to calculate and better to understand ( Garrett et al., 2010 ;Kannurpatti and Biswal, 2008 ). Unlike ALFF, SD BOLD can be easily applied to time-discontinuous fMRI data, such as data using blocked design, since frequencies estimated from discontinuous, concatenated data could be contaminated by the block boundaries.
Sample entropy, on the other hand, is generally used as a measure of signal complexity, although it is sometimes also used to describe signal variability. Given the different theorization between the sample entropy and the SD BOLD , these two indices might capture different aspects of signal variability. For example, sample entropy and SD BOLD measured by BOLD fMRI showed a significant anti-correlation ( Shafiei et al., 2019 ). Moreover, mixed results have been found regarding the relationship between entropy and FC ( McDonough and Nashiro, 2014a ;Liu et al., 2019 ;;McIntosh et al., 2014 ;Misic et al., 2011 ;Shafiei et al., 2019 ). For example, an EEG study separately examined short-range/high frequency (i.e., fine scale) temporal variability and long-range/low-frequency (i.e., coarse scale) temporal variability and revealed positive correlations between entropy and DC at both time scales ( Misic et al., 2011 ). A study using fMRI, however, revealed a negative association between entropy and FC at fine scales but a positive association at coarse scales ( McDonough and Nashiro, 2014a ). Another fMRI study observed a weak and negative relationship between increased entropy and decreased FC as a result of dopamine depletion via acute phenylalanine/tyrosine depletion ( Shafiei et al., 2019 ). Therefore, whether SD or sample entropy captures different aspects of brain information processing at different timescales and different brain regions remains to be examined. Future studies integrating and comparing these measures will yield more insights into brain dynamics and their functions.
Notably, the cross-voxel correlations between mfSD BOLD and DC were positive for all participants, but the cross-subject correlations were negative for some brain regions in the HCP dataset. The brain regions showing negative correlations were mainly located in subcortical regions (such as the hippocampus, amygdala, putamen, thalamus, etc.) and the cerebellum, which exhibited worse temporal signal-to-noise ratios than cortical regions ( Autio et al., 2020 ). Furthermore, we were unable to replicate the negative correlations in the Beijing dataset with a large sample ( N = 416). Future studies should examine the reliability of this negative correlation.
Finally, both mfSD BOLD and FC were correlated with regional CBF. This result replicates and extends previous studies showing a close relationship between regional CBF and FC ( Liang et al., 2013 ;Sheng et al., 2018 ;Tomasi et al., 2013 ). Our results further indicate that FC and signal variability might have a common physiological basis ( Tomasi and Volkow, 2018 ). Given the relatively consistent pattern cross-scanner for mfSD BOLD , we predict that this correlation between mfSD BOLD and CBF could hold across scanners, although empirical data are required to verify this prediction. Together, these findings provide clear and strong evidence to support the robust coupling between DC and mfSD BOLD .

Functional relevance of DC-mfSD BOLD coupling
The strong coupling between signal variability and functional integration within-subjects might reflect a synergy in nodal function and its network role in information processing. On the one hand, extant studies suggest that brain signal variability may represent a potentially meaningful and stable property of the human brain ( Zuo et al., 2010 ), which reflects neural dynamic range or flexibility ( Garrett et al., 2014 ;Garrett et al., 2013b ;Nomi et al., 2017 ) and the capability of information processing . Increased SD BOLD has been found for younger adults compared to elderly adults ( Garrett et al., 2010( Garrett et al., , 2013a( Garrett et al., , 2017. Furthermore, SD BOLD increases from rest to task and increases with task difficulty ( Garrett et al., 2014 ). Studies have also reported reduced regional variability in a variety of brain diseases, such as stroke ( Kielar et al., 2016 ), Alzheimer's disease ( Scarapicchia et al., 2018 ), multiple sclerosis ( Petracca et al., 2017 ), and other neurological disorders ( Zoller et al., 2017 ). Dopamine can boost SD BOLD , thereby improving individuals' cognitive performance ( Garrett et al., 2015 ;Guitart-Masip et al., 2016 ).
On the other hand, the DC of a node refers to the number of edges attached to the node ( Buckner et al., 2009 ), which quantifies the communication between brain regions required to exchange information. A larger DC of a brain area may indicate its greater role in information processing. Collectively, these findings suggest that the human brain as a self-organizing system could achieve a high level of matching capability and responsibility.
To test these ideas, we directly examined the functional relevance of the coupling for the first time. We found that the strength of DC-mfSD BOLD coupling was associated with cognitive function, while deviation from this coupling was related to worse cognitive functions and higher susceptibility to brain diseases. First, we found that the strength of mfSD BOLD and long-range DC coupling decreased with increasing age, which was accompanied by decreased cognitive performance. Similarly, Yan et al. (2017) found that the concordance among f/ALFF, DC and regional homogeneity was negatively correlated with age (range from 8 to 86 years old). Second, the degree of mfSD BOLD and long-range DC coupling was significantly related to individuals' comprehensive cognitive performance, even after controlling for age and sex. Third, the greater mismatches between long-range DC and mfSD BOLD in the right STG and MCC were associated with worse comprehensive cognitive performance. Finally, the brain regions with DC-mfSD BOLD mismatches were vulnerable to brain disorders. Together, these results suggest that the synergy of functional integration and brain signal variability plays an important role in the well-functioning brains of individuals.
In addition, we found that the coupling between the mfSD BOLD and full-range DC is stronger than that with long-range DC, but the coupling strength of the mfSD BOLD and long-range DC is more reliable and is strongly and reliably predictive of cognitive functions and susceptibility to brain diseases. It should be noted that the contribution of ultrashort connections (Euclidean distance < 20 mm) to full-range DC was excluded to reduce noise ( Power et al., 2012( Power et al., , 2011. These results suggest that long-range DC might be more relevant to higher-level brain functions while short-range connections (20 mm ≤ Euclidean distance < 75 mm) present limited relevance for comprehensive cognitive ability. In particular, long-range connections could provide quick links to other brain regions in the network ( Achard et al., 2006 ) and integrate information between segregated modules ( He et al., 2009 ), thus enabling efficient information processing. As such, the coupling between long-range DC and mfSD BOLD is more critical to achieve optimal brain functions than full-range DC and mfSD BOLD coupling.

Implications and future directions
The strong coupling between signal variability and FC and its correlation with brain functions provides a new perspective to improve the present understanding of how the brain works. A functioning brain re-quires good synergy between node capability and responsibility in network information integration. Several lines of future studies could be fruitful. First, the current study revealed a small yet significant age effect on the strength of long-range DC and mfSD BOLD coupling in our sample with a narrow age span. A previous study found weak and heterogeneous coupling between DC and SD BOLD in children and young adolescents, with correlation coefficients varying from -0.2 to 0.5 ( Sato et al., 2019 ), and these findings were most likely due to the age effect. Future studies can expand the age range to examine life span changes in this coupling, which could improve the present understanding of the developmental changes in brain functions. Second, our findings suggest that the coupling between BOLD signal variability and degree centrality is robust, reliable, and functionally significant, which raises the question regarding the biological meaning of this coupling. Although a systematic examination of this important question is certainly beyond the scope of this study, we did find that both mfSD BOLD and DC were related to CBF, indicating that the coupling of mfSD BOLD and DC might share a common metabolic basis. Future studies along this direction could lead to a more mechanistic understanding of this coupling. Third, future studies should examine DC-mfSD BOLD coupling between healthy individuals and patients with brain diseases and examine the link between the degree of mismatches and the severity of brain diseases, which could improve our understanding of brain functional abnormalities and identify targets for intervention. Fourth, although our results indicate that the mismatch between mfSD BOLD and DC has functional consequences, precisely what could contribute to the mismatch and the neural mechanisms underlying the association between the mismatch and brain function still require further study. Finally, future studies could use pharmacological manipulation or noninvasive brain stimulations to examine the causal relationship between functional integration and signal variability, which could contribute to the understanding of the origins of the coupling.

Conclusions
To summarize, our findings suggest that the human brain, as a complex functional system, exhibits strong coupling between regional signal variability and functional integration (as measured by both full-range and long-range DC). Furthermore, the greater coupling is predictive of higher cognitive ability, and the higher decoupling is associated with worse brain function and susceptibility to brain diseases. Our study not only provides a novel brain marker for brain functions but also provides a new perspective on how the brain is organized to achieve efficient information transformation and integration.

Declaration of Competing Interest
The authors declare no conflict of interest.

Acknowledgment
This work was sponsored by the National Natural Science Foundation of China ( 31730038 ), the NSFC and the Israel Science Foundation (ISF) joint project ( 31861143040 ), the Sino-German Collaborative Research Project "Crossmodal Learning" (NSFC 62061136001/DFG TRR169 ), and the Guangdong Pearl River Talents Plan Innovative and Entrepreneurial Team grant # 2016ZT06S220 . One of the datasets was provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. We thank Dr. Li-Xia Yuan (Assistant researcher of Hangzhou Normal University) for her help on data analysis.

Code and data availability
Codes supporting this study are freely available at the GitHub repository ( https://github.com/Cynthia1229/mfSD-DC ). Data from HCP is publicly available at https://db.humanconnectome.org/ . The Beijing and Hangzhou datasets are available from the corresponding author upon reasonable request. Sharing and re-use of the datasets require a formal data sharing agreement, as well as approval from the relevant Institutional Review Boards.

Ethics statement
The reanalysis of the Human Connectome Project (HCP) Open Access and Restricted Data and the reporting of results pertaining to this dataset had been reviewed and approved by the Beijing Normal University Institutional Review Board. The collection and analysis of the Beijing sample data was approved by the Institutional Review Board at Beijing Normal University. The collection and analysis of Hangzhou sample data was approved by the Ethics Committee of the Sir Run Run Shaw Hospital, School of Medicine, Zhejiang University, and the Affiliated Hospital of Hangzhou Normal University.