Similar to previous studies (Argollo de Menezes and Barabasi, 2004), to model the signal S(t) we simulated the random diffusion of W walkers (messages) on a network of N nodes described by its adjacency matrix, A_ij. Each walker was placed at a randomly chosen network node from which it departed randomly along one of the edges of that node in the next time step. This diffusion process was independently repeated 1,200 times and we recorded the number of incoming visits by various walkers at each network node to compute the time-varying signal at each node, S_i(t). Temporal fluctuations in W were used to simulate externally induced modulations in S_i(t), which for random networks and scale-free networks results in b = 1 exponent in (1) (Argollo de Menezes and Barabasi, 2004). Thus we varied the number of walkers as a function of time as: W(t) = W + ξ(t), where ξ(t) was a uniformly distributed random variable in the interval [−ΔW, ΔW], with ΔW = k^*10³ and k = 0,1,2,…, 9, and W = 10⁴.

The FreeSurfer gray matter parcellations (wmparc.2.nii) for 7 randomly selected MRI datasets were used to determine imaging voxels in the occipital, cingulate and insular networks (Figure 1A). The occipital network comprised bilateral cuneus, lateral occipital, lingual, and pericalcarine cortices (number of nodes/voxels, N = 8,200 ± 600). The cingulate network comprised bilateral rostral anterior, caudal anterior, isthmus, and posterior cingulate (N = 3,100 ± 400). The insular network comprised the bilateral insula (N = 2,300 ± 100). The Pearson correlation was used to compute correlation matrices reflecting the functional connectivity between voxels within each network for each subject. A correlation threshold R = 0.2 (p < 0.05) was used to compute the corresponding binary adjacency matrices. We implemented the diffusion model described above (Argollo de Menezes and Barabasi, 2004). We assumed the signal is proportional to the rate of incoming messages at each node as a function of time, which was simulated using 1,200 steps.

To simulate the dynamics in the amplitude of the signal fluctuations, δ_i, at each node we segmented the S_i(t) data (1,200 time points) into 23 epochs (window length: 100 time points; window shift: 50 time points) using a popular rectangular sliding window approach (Chang and Glover, 2010). The temporal standard deviation of S_i(t) during each epoch was used to estimate δ_i. Degree, D_i, the number of links connected to a network node (Rubinov and Sporns, 2010), was computed for each of the 23 epochs of the synthetic S_i(t) data using a correlation threshold, R > 0.5 (p < 10⁻⁷). The linear model log(σ_X) = log(a) + blog〈X〉 with 2 freely adjustable parameters: log(a) and b, was used to fit the power law (1) to the temporal mean and dispersion values of the dynamic δ and D metrics (X).

To test the predictions of the random diffusion model we analyzed rfMRI datasets drawn from the Human Connectome Project (HCP; http://www.humanconnectome.org/). No experimental activity with any involvement of human subjects took place at the author's institutions. The 66 participants (age: 30 ± 3 years; 32 females; Subject IDs: 100408, 103515, 103818, 105115, 105216, 106319, 110411, 118730, 118932, 119833, 120212, 122317, 123117, 125525, 127933, 128632, 129028, 130013, 131924, 133625, 133827, 133928, 134324, 136833, 137128, 138231, 138534, 140824, 142828, 143325, 144226, 149337, 149539, 150423, 151526, 153429, 156637, 158540, 159239, 159340, 160123, 161731, 162329, 163129, 165840, 167743, 172332, 178950, 182739, 191437, 192439, 192540, 194140, 197550, 199150, 199251, 200614, 201111, 210617, 217429, 249947, 250427, 255639, 304020, 307127, 329440) of the WU-Minn HCP Q1 data release included in this study provided written informed consent according to procedures approved by the IRB at Washington University in St. Louis.

Resting-state (eyes open) functional images were acquired using a gradient-echo-planar (EPI) sequence with multiband factor 8, TR 720 ms, TE 33.1 ms, flip angle 52°, 104 × 90 matrix size, 72 slices, 2 mm isotropic voxels, and 1200 timepoints (Smith et al., 2013; Uğurbil et al., 2013). Scans were repeated twice using different phase encoding directions (LR and RL) on each of two imaging sessions (REST1 and REST2) collected on different days. The “minimal preprocessing” datasets, which include gradient distortion correction, rigid-body realignment, field-map processing, spatial normalization to the stereotactic space of the Montreal Neurological Institute (MNI), high pass filtering (1/2,000 Hz frequency cutoff) (Glasser et al., 2013), independent component analysis-based denoising (Salimi-Khorshidi et al., 2014), and brain masking were used in this study.

Framewise displacements, FD, computed for every time point from head translations and rotations using a radius of r = 50 mm (Power et al., 2012) did not differ between MRI sessions or phase encoding directions across subjects (p > 0.2, paired t-test; 〈FD〉 = 0.176 ± 0.05 mm). Scrubbing was not implemented to preserve the frequency spectra used for the computation of ALFF. Multilinear regression of head translations and rotations were used to minimize motion related fluctuations in the MRI signals (Tomasi and Volkow, 2010). Standard 0.01–0.08 Hz band-pass filtering was used to minimize physiologic noise of high frequency components.

The average of the power spectrum's square root in the 0.01–0.08 Hz low frequency bandwidth was used to compute the ALFF (Yang et al., 2007). Functional connectivity density mapping was used to compute the lFCD (Tomasi and Volkow, 2010) at three different thresholds R > 0.3, 0.4 and 0.5. A sliding window approach (Chang and Glover, 2010) with two different window lengths (72s and 144s) and two different window shapes (rectangular and Hamming) was used to compute dynamic ALFF and lFCD maps with 2-mm isotropic resolution at two different temporal resolutions (36s and 72s). The window shift was set as half of the window length.

To test the power law (1) we contrasted scaling factors for the simulated signal fluctuations (δ) and degree (D) against those for ALFF and lFCD. Since lFCD has high sensitivity and specificity for gray matter (Tomasi et al., 2016c), the FC metrics were averaged within the anatomical gray matter regions-of-interest for each individual to minimize confounds arising from the variability of the folding patterns of cortical gray matter. Specifically, the FreeSurfer gray matter parcellations (wmparc.2.nii) were used as ROIs to compute the averages of the temporal mean and dispersion values of ALFF(t) and lFCD(t) within 34 cortical and 9 subcortical gray matter regions in each brain hemisphere. A probabilistic atlas for each of the gray matter parcellations was developed by averaging each of the gray matter parcellations across subjects independently, and used to display ROI results (i.e., b_ALFF or b_lFCD) in the MNI stereotactic space (Figures 2A, 3A).

The linear model log(σ) = log(a)+blog〈X〉 with 2 free adjustable parameters: log(a) and b, was used to fit the power law (1) to the temporal mean and dispersion values of the dynamic ALFF and lFCD metrics (X). Paired t-test was used to assess within subjects differences in b_ALFF and b_lFCD as a function of session, phase encoding direction, correlation threshold, and window length and shape. Two samples t-test and Pearson correlation were used to assess gender and aging effects on b_ALFF and b_lFCD.

The power law (1) fitted well (R² > 0.8) the temporal mean and standard deviation values of S_i(t) across nodes. The b exponents increased monotonically with ΔW, which is consistent with the notion that internal randomness (diffusion) and external modulation (ΔW) proportionally alter S_i(t) in the network (Argollo de Menezes and Barabasi, 2004). Thus, ΔW contributed to the temporal variability of the signal at each network node, gradually increasing b from ½ to 1 in all three brain networks (Figure 1B) as it occurs in other complex networks. Thus, if its magnitude is significant (ΔW ~ ½ 〈W〉), the external modulation can dominate the dynamics of S_i(t).

The mean and dispersion values of δ_i computed across epochs were also in good agreement with the power law (1). Our simulations suggest that, b_δ ~ 1, when internal randomness dominates over the external modulations (Figure 1C). However, b_δ decreased with the amplitude of the external modulation and was constrained in the interval [0.5, 1]. Similarly, the mean and dispersion values of D_i computed across epochs were in good agreement with the power law (1). Our simulations suggest that b_D ~ 1 when the external modulation dominates over the internal randomness, but b_D increases significantly above 1 when the relative weight of the external modulation decreases (Figure 1D). The power law failed to fit the data when internal randomness dominated over the external modulation (ΔW/W > ½) suggesting lack of association between the mean and dispersion values of D in this regime.

A linear fit of whole-brain average and dispersion values of ALFF on a log-log plot computed across nodes demonstrated good agreement between Equation (1) and the dynamic amplitude of the signal fluctuations in each of the individual ROI (b_ALFF = 0.66 ± 0.16, mean ± standard deviation; 28 < t-score <294; P < 1E-37; Figure 2A). Consistent findings emerged from average and dispersion values within anatomical regions, independently for each of the 86 gray matter ROIs (R² > 0.8; Figures 2B,C). The average scaling exponent was not different for subcortical (cerebellum, thalamus, caudate, putamen, pallidum, hippocampus, amygdala, accumbens and ventral diencephalon; b_ALFF = 0.78 ± 0.04, mean ± standard error) than for cortical (b_ALFF = 0.81 ± 0.07) regions (p > 0.4), independent of session, phase encoding direction (LR vs. RL), sliding window length (72s vs. 144s) and shape (rectangular vs. Hamming).

There were no significant differences in b_ALFF across subcortical regions. However, b_ALFF varied significantly across cortical regions. Specifically, the scaling exponent was higher for limbic (cingulum, orbitofrontal, parahippocampal and entorhinal) and visual (lingual, fusiform and pericalcarine) areas, the temporal and insular cortices and pars orbitalis (b_ALFF = 0.89 ± 0.02) than for occipital (cuneus, lateral occipital), parietal (inferior, superior, precuneus, postcentral), language (opercularis, triangularis, supramarginal) and prefrontal (paracentral, precentral, rostral, middle and superior frontal) areas (b_ALFF = 0.72 ± 0.03; p < 10⁻⁹; Figure 2B). The scaling exponent had normal distribution (center b_ALFF = 0.80; width = 0.16) across the 86 gray matter ROIs (R² = 0.999, Gaussian fit; Figure 2C).

Significant between-subjects variability in the scaling exponent emerged from the data when we fitted Equation (1) to the mean and dispersion values of ALFF across the 43 ROIs, independently for each individual (b_ALFF = 0.66 ± 0.05; Figure 2D) and with similar robustness (R² > 0.96). The scaling exponent slightly decreased with age (slope = −0.03/decade; R = −0.234; p = 0.03, one-tailed). However, there were no significant gender differences (p > 0.77; two-tailed two-sample t-test) in b_ALFF.

Similar to ALFF, a linear fit of whole-brain average and dispersion values of lFCD on a log-log plot computed across nodes demonstrated good agreement between Equation (1) and the local degree of brain functional connectivity in each of the individual ROI (b_lFCD = 1.05 ± 0.17, mean ± standard deviation; 43 <t-score <179; P <3E-49; Figure 3A). Average and dispersion values within anatomical regions showed consistent findings with those from the whole-brain analysis (R² > 0.8; Figure 3B). The average scaling exponent was higher for subcortical (b_lFCD = 1.23 ± 0.09, mean ± standard error) than for cortical (b_lFCD = 1.06 ± 0.10) regions (p < 10⁻³, two-tailed two-sample t-test), independent of the correlation threshold used in the computation of the lFCD (R > 0.3, 0.4, or 0.5), session, phase encoding direction, window length (72s vs. 144s) and shape.

For lFCD, the scaling exponent was higher for limbic (cingulum, orbitofrontal, parahippocampal, and entorhinal), language (opercularis, orbitalis, triangularis), temporal (inferior, middle superior), and frontal (paracentral, superior and pole), insula and fusiform gyrus (b_lFCD = 1.13 ± 0.07) than for occipital (cuneus, lateral occipital, lingual and pericalcarine), parietal (inferior, superior, precuneus, supramarginal, paracentral, postcentral), prefrontal (precentral, rostral, middle, and superior) and temporal (entorhinal temporal pole, transverse) areas (b_lFCD = 0.98 ± 0.06; p < 10⁻⁶; Figure 3B). Across the 86 gray matter ROIs the scaling exponent had a right-skewed distribution with peak at b_lFCD = 1.03 and width = 0.17 (R² = 0.95, Gaussian fit; Figure 3C). Fitting mean and dispersion values of lFCD across the 43 gray matter ROIs, independently for each subject, revealed modest between-subjects variability in the scaling exponent (b_lFCD = 1.09 ± 0.06; Figure 3D), and b_lFCD did not show significant age or gender differences (p > 0.23).

In visual areas (pericalcarine, lateral occipital, and cuneus) the standardized scaling factors b_z were lower than average and were significantly lower for lFCD than for ALFF (p < 0.0005, t-test; Figure 4, left). In prefrontal regions (middle Frontal, superior frontal, precentral, paracentral, pars opercularis, and caudal anterior cingulate), the lower than average b_z was lower for ALFF than for lFCD (p < 0.001). In frontal and temporal poles, entorhinal and lingual cortex showed the higher than average b_z was higher for ALFF than for lFCD (p < 0.0001; Figure 4 right). In anterior (rostral) and posterior (isthmus) cingulate, fusiform gyrus and subcortical regions (hippocampus, thalamus and cerebellum) the higher than average b_z was higher for lFCD than for ALFF (p < 0.0006).

Given that frequency information may be of interest and that the ICA-FIX denoising procedure can remove a significant fraction of the physiological noise of respiratory origin (Salimi-Khorshidi et al., 2014), we also computed dynamic lFCD measures without 0.01–0.08 Hz bandpass filtering to assess the effect of higher frequencies on the power scaling law (Equation 1). Without bandpass filtering the scaling exponent b of the dynamic lFCD metrics was significantly larger than with bandpass filtering (p < 0.0001; Figure 5), and the agreement between the data and Equation (1) was significantly reduced [R² = 0.82 (without) and 0.96 (with bandpass filtering)].

Note that b = ½ emerges either from diffusion or from flow models, independently of the number of steps in the diffusion model, and from random networks as well as from scale-free networks. This indicates that b = ½ is not a particular property of the random diffusion model, but it is shared by several dynamic processes (Argollo de Menezes and Barabasi, 2004). Our computational resources did not allow demanding whole brain network simulations at 2-mm isotropic resolution (~10⁵ nodes/voxels). Thus, our simulations suggesting that when internal randomness dominates over the external modulations (ΔW/W ~ 0) b_δ ~ 1 and b_D > 1, but when external modulations dominate over internal randomness (ΔW/W ~ 1) b_δ ~ 0.5 and b_D ~ 1 are limited to the 3 exemplary networks in this work. However, it is likely that they apply also to the whole brain. Instrumental noise likely resulted in overestimations of intrinsic randomness in subcortical regions for which the 32 channel RF coil used by the HCP has low sensitivity. Since the theoretical model was developed across network nodes, the interpretation of the power law across ROIs and subjects could be considered controversial. Our empirical evidence, however, suggests that the temporal mean and standard deviation values of dynamic functional connectivity metrics also adhere to a power law computed across ROIs or subjects, which are consistent with the power law computed across nodes (i.e., across nodes of each individual ROI, b_lFCD = 1.05 ± 0.17 mean ± standard deviation; across the 86 gray matter, ROIs b_lFCD = 1.03 ± 0.17; across 66 subjects, b_lFCD = 0.98 ± 0.16). This suggests similar effects of randomness and external modulators on power scaling factors computed across network nodes, ROIs or subjects.

Dynamic lFCD is restricted to the local functional connectivity cluster. We did not assess the dynamics of global functional connectivity density (gFCD) because at high spatiotemporal resolution gFCD is extremely demanding and beyond our computational resources. However, this is not a strong limitation because previous studies have shown that the lFCD and gFCD metrics are proportional to one another (Tomasi and Volkow, 2011).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.