Brain Network Dynamics Adhere to a Power Law

The temporal dynamics of complex networks such as the Internet are characterized by a power scaling between the temporal mean and dispersion of signals at each network node. Here we tested the hypothesis that the temporal dynamics of the brain networks are characterized by a similar power law. This realization could be useful to assess the effects of randomness and external modulators on the brain network dynamics. Simulated data using a well-stablished random diffusion model allowed us to predict that the temporal dispersion of the amplitude of low frequency fluctuations (ALFF) and that of the local functional connectivity density (lFCD) scale with their temporal means. We tested this hypothesis in open-access resting-state functional magnetic resonance imaging datasets from 66 healthy subjects. A robust power law emerged from the temporal dynamics of ALFF and lFCD metrics, which was insensitive to the methods used for the computation of the metrics. The scaling exponents (ALFF: 0.8 ± 0.1; lFCD: 1.1 ± 0.1; mean ± SD) decreased with age and varied significantly across brain regions; multimodal cortical areas exhibited lower scaling exponents, consistent with a stronger influence of external inputs, than limbic and subcortical regions, which exhibited higher scaling exponents, consistent with a stronger influence of internal randomness. Findings are consistent with the notion that external inputs govern neuronal communication in the brain and that their relative influence differs between brain regions. Further studies will assess the potential of this metric as biomarker to characterize neuropathology.


INTRODUCTION
During resting-state functional magnetic resonance imaging (rfMRI) (Biswal et al., 1995) the human brain sequentially engages in a series of diverse free-streaming subject-driven mental states supported by different brain networks (Mason et al., 2007;Doucet et al., 2012;Shirer et al., 2012;Liu and Duyn, 2013;Yang et al., 2014). These complex and time-varying functional operations require a dynamic brain network topology to support the context-dependent coordination of neuronal populations (Allen et al., 2014;Zalesky et al., 2014) and its characterization and measurement could facilitate development of clinical biomarkers in neurology and psychiatry (Hutchison et al., 2013). Thus, the temporal dynamics of the human brain connectome (Chang and Glover, 2010;Sakoglu et al., 2010) provides a new metric of brain function to assess healthy and disease conditions . However, our lack of understanding of the principles governing network dynamics may preclude the interpretation of the observed dynamics, which increases the within-subjects variability of the functional connectivity metrics (Tomasi et al., 2016a,b). A better understanding of how the collective behavior of neuronal communities contributes to the observable dynamics is crucial for the interpretation of the dynamics of functional connectivity.
Previous studies have shown that temporal mean S i , and dispersion, σ i , of the activity at a given node are related through a power law across network nodes (Argollo de Menezes and Barabasi, 2004) where the scaling exponent, b, is a property of the network. Based on theoretical grounds and independent from the topology of the network, b equals either ½ or 1, which reflect a competition between the system's internal collective dynamics and temporal changes in the external environment (Argollo de Menezes and Barabasi, 2004). Specifically, in the absence of external modulation, b = ½, but when external driving forces become dominant, b = 1. For instance, whereas the network of internet routers is characterized by b = ½, the network of highways and the World Wide Web are characterized by b = 1. However, empirical evidence from ecology, where (1) describes the spatiotemporal variability of natural populations, supports the existence of intermediate b-values (Taylor, 1961) suggesting that meaningful temporal dynamic require ½ < b < 1.
Inasmuch as brain networks have scale-free (Barabasi and Albert, 1999;Eguíluz et al., 2005) and small-world (Watts and Strogatz, 1998) properties exhibited by complex networks we hypothesized that the mean and σ of FC properties such as ALFF, the amplitude of the low frequency fluctuations (Yang et al., 2007) or lFCD, the local degree of connectivity (Tomasi and Volkow, 2010) would reveal the characteristic power scaling properties exhibited by other complex networks. Specifically, we hypothesized that the mean and σ would be related by the power law (1) and that different brain networks would exhibit different scaling exponents reflecting differential balance between internal randomness (random firing) and external inputs (non-random firing). We selected functional connectivity (ALFF and lFCD) metrics rather than raw signals because the mean and dispersion values of the BOLD-fMRI signals are not expected to be in agreement with Equation (1).

METHODS
To interpret the observed power scaling law (1), we study a simple dynamical model based on random diffusion. Using this model and functional connectivity information extracted from rfMRI datasets, we assessed the validity of Equation (1) in the context of brain functional connectivity. However, since direct application of Equation (1) to the mean and dispersion values of the raw fMRI time series is meaningless (the MRI signal mainly reflects tissue properties such as water density and T1 and T2 relaxation rates, which do not change as a function of time; the BOLD signal is zero-mean by definition), we simulated the temporal dynamics of ALFF and lFCD.

Model
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 3 and k = 0,1,2,..., 9, and W = 10 4 .

Simulations
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 −7 ). The linear model log(σ X ) = log(a) + b log 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).

Preprocessing
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 ttest; 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 bandpass filtering was used to minimize physiologic noise of high frequency components.

Dynamic ALFF and lFCD
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.

Region-of-Interest (ROI) Analysis
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 regionsof-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).

Statistical Methods
The linear model log(σ ) = log(a) + b log 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 .

Simulations
The power law (1) fitted well (R 2 > 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.

Amplitude of Fluctuations
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 2 > 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).
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 2 > 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 .

Local Degree
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 2 > 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 −3 , 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.
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).

Effect of Bandpass Filtering
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 2 = 0.82 (without) and 0.96 (with bandpass filtering)].

DISCUSSION
Here we show for the first time that the mean and the dispersion values of dynamic FC metrics such as ALFF or lFCD are linked by a power law (1). This characteristic of complex networks such as rivers and highways networks, the Internet and the World Wide Web (Argollo de Menezes and Barabasi, 2004), and many biological systems (Taylor, 1961), reflects the competition between the system's internal collective dynamics and changes in the external environment. This strongly suggests that the dynamics of the FC metrics embeds important functional information, a possibility previously highlighted (Hutchison et al., 2013;Calhoun et al., 2014;Rashid et al., 2014;Hutchison and Morton, 2015), which could help in the development of biomarkers of brain function.
Our simulations were based on a random diffusion model previously proposed by Argollo de Menezes and Barabasi to explain the power scaling between the mean and the dispersion of the signals observed in natural and technological networks (Argollo de Menezes and Barabasi, 2004). Whereas the approach by Argollo de Menezes and Barabasi was based on random and scale-free networks (Barabasi and Albert, 1999;Barabási, 2009), the present approach was based on real FC networks directly extracted from in vivo resting fMRI data. The scaling exponents for the brain in the present work are consistent with those obtained previously in random and scale-free networks (Argollo de Menezes and Barabasi, 2004).
Here we extended the random diffusion model in order simulate the amplitude of spontaneous signal fluctuation and the degree of connectivity. Our simulations suggest that under FIGURE 4 | Brain regions showing statistically significant differences in standardized scaling exponents, b z , for ALFF and lFCD at p < 0.001. ACC, anterior cingulate cortex.
FIGURE 5 | Scatter plot demonstrating the effect of bandpass filtering on the power scaling (Equation 1) of lFCD across ROIs. The fitted slope, b, the scaling factor in Equation (1), is significantly steeper without bandpass filtering (standard errors in parenthesis) suggesting increased level of randomness (see Figure 1D). pure randomness (i.e., without external driving forces, W = 0) the mean and the dispersion values of the amplitude of signal fluctuations and degree are associated by power laws with scaling exponents b δ = 1 and b D > 1, respectively. However, under the influence of dynamic external modulations ( W/W ∼ 1), b δ < 1 and b D = 1 characterize the dynamic behavior of the signal fluctuations and degree. The analysis of variability of restingstate fMRI datasets from the HCP database shows a range of scaling exponents for ALFF (0.5 < b ALFF < 1) and for lFCD (1 < b lFCD ), which is consistent with the presence of dynamic external modulations of brain activity (0.5 < b δ < 1) and the corresponding degree (1 < b D ). Overall, our findings are also consistent with the existence of dynamic modulations of brain activity that may reflect orchestrated dynamic neural processing (Yu et al., 2012;Allen et al., 2014;Gonzalez-Castillo et al., 2015). This is the first study to document differences in scaling exponents between brain regions. Multimodal association areas (opercularis, triangularis, rostral, middle and superior frontal, precentral and paracentral, inferior and superior parietal and precuneus), somatosensory (supramarginal, postcentral) and visual (cuneus, lateral occipital) unimodal association areas showed low scaling exponent both for ALFF (b ALFF ∼ 0.7) and for lFCD (b lFCD ∼ 1). These findings suggest that the dynamics of the FC metrics was driven by external inputs ( W/W > ½) rather than by internal random processes ( W/W < 0.5; Figures 1C,D), which is also consistent with the existence of dynamic modulations of resting brain activity (Yu et al., 2012;Allen et al., 2014;Gonzalez-Castillo et al., 2015). The multimodal cortex is highly interconnected with higher-order association areas involved in cognition and motor planning (Goldman-Rakic, 1988). Thus dynamic engagement of functional connectivity hubs in multimodal and unimodal association cortices may explain the low scaling exponent in these regions. On the other hand, limbic and subcortical regions exhibited relatively higher scaling exponents (b ALFF ∼ 0.8 and b lFCD ∼ 1.2) suggesting a stronger influence of internal randomness in the resting dynamics of the FC metrics in these regions. We identify regional differences in the influence of internal randomness for different FC metrics. The direct comparison of standardized measures suggests a weaker influence of randomness in visual areas for lFCD than for ALFF and in prefrontal areas for ALFF than for lFCD, and a stronger influence of randomness in subcortical and limbic regions for lFCD than for ALFF. ALFF and lFCD reflect different network properties. Whereas ALFF is proportional to the BOLD signal fluctuations that reflect neuronal communication (Logothetis et al., 2001), the synchronous fluctuations of local communities measured by lFCD reflects the local degree of connectivity (Tomasi and Volkow, 2010).
The scaling exponent for ALFF, and to a lesser extent for lFCD, showed significant variability ( b ALFF = 12%; b lFCD = 9%) across subjects suggesting that the dynamics of the b has potential as a biomarker for psychiatry and neurology. To illustrate the potential of this metric here we show that even in a relatively small sample (66 subjects) with narrow age range (22-35 years), b ALFF is sensitive to aging effects, consistent with previous studies in large samples (∼1000 subjects) with wide age range (17-82 years) that documented age-related decreases in FC (Biswal et al., 2010;Tomasi and Volkow, 2012).
The scaling exponent for lFCD increased significantly above 1 when frequencies other than those in the 0.01-0.08 Hz band were not removed from the data. At the same time, the agreement with a power scaling was reduced when Equation (1) was fitted to the data without bandpass filtering. This likely reflects the introduction of additional randomness and is consistent with increased noise level and lack of additional information at higher frequencies than those in the 0.01-0.08 Hz band. The brain normally operates under certain level of randomness that is important for multiple operation including perception and decision-making. The relevance of internal neuronal noise has been most extensively studied for visual perception (Brascamp et al., 2006;Kim et al., 2006). Theoretical studies have also shown that randomness may influence behavioral responses when there are multiple routes to action and suggested that noise generated by random firing rates of neurons can be used to predict a decision (Rolls, 2012). Since limbic and subcortical regions support automatic, implicit decision making (Floresco et al., 2008;Mitchell, 2015) the higher scaling exponents in these regions suggests an important role of randomness in implicit decision making processes. The sensitivity to randomness of b could be useful for studying psychiatric disorders such as autism, which is associated with increased randomness of endogenous brain oscillations (Lai et al., 2010).

Study Limitations
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 scalefree 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 5 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).

AUTHOR CONTRIBUTIONS
DT designed the study, carried the analyses and wrote the manuscript. ES developed imaging tools and wrote the manuscript. NV wrote the manuscript.