Scaling Properties of Human Brain Functional Networks

TL;DR

This study analyzes high-resolution human brain functional networks to determine their degree distribution, finding they generally do not follow a power-law but tend towards a generalized Pareto distribution, impacting understanding of brain hubs.

Contribution

It provides a comprehensive analysis using high-resolution data and multiple distribution models, clarifying the nature of brain network degree distributions.

Findings

Degree distributions tend towards the generalized Pareto model.

Results do not support power-law distribution in brain networks.

Implications for the number of hubs in brain networks.

Abstract

We investigate scaling properties of human brain functional networks in the resting-state. Analyzing network degree distributions, we statistically test whether their tails scale as power-law or not. Initial studies, based on least-squares fitting, were shown to be inadequate for precise estimation of power-law distributions. Subsequently, methods based on maximum-likelihood estimators have been proposed and applied to address this question. Nevertheless, no clear consensus has emerged, mainly because results have shown substantial variability depending on the data-set used or its resolution. In this study, we work with high-resolution data (10K nodes) from the Human Connectome Project and take into account network weights. We test for the power-law, exponential, log-normal and generalized Pareto distributions. Our results show that the statistics generally do not support a power-law, but instead these degree distributions tend towards the thin-tail limit of the generalized Pareto model. This may have implications for the number of hubs in human brain functional networks.

Tables

Table 1. Table 1 : Fit results of the exponential ( p ​ ( x ) = C ​ e − λ ​ x 𝑝 𝑥 𝐶 superscript 𝑒 𝜆 𝑥 p(x)=Ce^{-\lambda x} ) and the power law distributions ( p ​ ( x ) = C ​ a ​ x − α 𝑝 𝑥 𝐶 𝑎 superscript 𝑥 𝛼 p(x)=Cax^{-\alpha} )

Thr	$\lambda$	$x_{min}$	TL	KS	p   -value	TLr
$+0.8$	0.000	0.000	0.0	0.000	0.0000	0.714
$+0.7$	0.162	1.186	83.5	0.056	0.3125	0.500
$+0.6$	0.089	3.887	130.5	0.058	0.4840	0.306
$+0.5$	0.038	15.544	175.5	0.050	0.6170	0.187
$+0.4$	0.020	23.051	440.5	0.037	0.7305	0.160
$+0.3$	0.012	44.067	805.0	0.030	0.3305	0.167
$+0.2$	0.008	44.049	1552.5	0.024	0.0895	0.205
$+0.1$	0.006	25.520	4350.5	0.031	0.0015	0.433
$+0.0$	0.007	545.110	1597.0	0.037	0.0290	0.159
$-0.0$	0.021	123.255	1678.0	0.021	0.0485	0.168
$-0.1$	0.023	53.151	141.0	0.072	0.1835	0.014
$-0.2$	0.095	7.104	123.0	0.170	0.0170	0.105
$-0.3$	0.348	0.330	95.5	0.224	0.0000	0.534
$-0.4$	0.000	0.000	0.0	0.000	0.0000	0.876
$-0.5$	0.000	0.000	0.0	0.000	0.0000	0.881
$-0.6$	0.000	0.000	0.0	0.000	0.0000	0.510
$-0.7$	0.000	0.000	0.0	0.000	0.0000	-
$-0.8$	0.000	0.000	0.0	0.000	0.0000	-

Table 2. Table 2 : Fit results of generalized Pareto distribution ( p ​ ( x ) = 1 σ ​ ( 1 + k ​ x − x m ​ i ​ n σ ) − 1 − 1 k 𝑝 𝑥 1 𝜎 superscript 1 𝑘 𝑥 subscript 𝑥 𝑚 𝑖 𝑛 𝜎 1 1 𝑘 p(x)=\frac{1}{\sigma}(1+k\frac{x-x_{min}}{\sigma})^{-1-\frac{1}{k}} ) and of log-normal distribution ( C ​ 1 x ​ e ​ x ​ p ​ [ − ( l ​ n ​ ( x ) − μ ) 2 2 ​ σ 2 ] 𝐶 1 𝑥 𝑒 𝑥 𝑝 delimited-[] superscript 𝑙 𝑛 𝑥 𝜇 2 2 superscript 𝜎 2 C\frac{1}{x}exp[\frac{-(ln(x)-\mu)^{2}}{2\sigma^{2}}] ).

Thr	$k$	$\sigma$	$x_{min}$	TL	KS	p   -value	TLr
$+0.8$	0.000	0.000	0.000	0.0	0.000	0.0000	0.696
$+0.7$	0.000	3.663	1.539	58.5	0.049	0.1315	0.531
$+0.6$	0.118	11.094	2.943	172.5	0.048	0.3305	0.440
$+0.5$	0.031	28.656	20.072	338.5	0.034	0.7380	0.229
$+0.4$	0.001	62.035	34.191	372.5	0.027	0.6155	0.174
$+0.3$	-0.128	102.619	85.215	638.5	0.021	0.7975	0.176
$+0.2$	-0.178	152.365	136.945	938.5	0.014	0.8275	0.135
$+0.1$	-0.226	241.640	201.475	2387.0	0.010	0.7555	0.238
$+0.0$	-0.224	218.500	380.050	4681.5	0.009	0.6750	0.467
$-0.0$	0.182	32.619	127.480	1195.5	0.012	0.8425	0.119
$-0.1$	0.426	12.883	9.690	1469.0	0.020	0.1725	0.149
$-0.2$	0.439	2.549	0.763	820.5	0.030	0.1890	0.457
$-0.3$	0.433	0.996	0.649	101.5	0.062	0.0105	0.405
$-0.4$	0.000	0.000	0.000	0.0	0.000	0.0000	0.253
$-0.5$	0.000	0.000	0.000	0.0	0.000	0.0000	1.000
$-0.6$	0.000	0.000	0.000	0.0	0.000	0.0000	1.000
$-0.7$	0.000	0.000	0.000	0.0	0.000	0.0000	-
$-0.8$	0.000	0.000	0.000	0.0	0.000	0.0000	-

Full text

1 Introduction

Much interest in theoretical neuroscience has revolved around graph-theoretic scaling properties of the network of structural and functional correlations in the human brain. Some authors have described the degree distribution of nodes in brain functional networks as scale-free; that is, these networks follow a power-law degree distribution $P(k)\sim k^{-\alpha}$ with an exponent close to 2 [7, 11], indicating the presence of a small number of hub-nodes that connect widely across the network. Other studies have suggested that functional brain networks are not scale-free, but instead are characterized by an exponentially truncated distribution [1, 8, 10]. The scaling characteristics of these networks are associated with the number and organization of network hubs and consequently may have implications with our understanding of how the brain responds to disease or damage [2, 1]. This calls for a rigorous statistical methodology to infer underlying models that best describe the degree distribution of brain functional networks.

Initial studies were based on least-square fitting of log-log plots of frequency distributions to answer this question. This approach, although seemingly straightforward, is inadequate from a statistical point of view as elaborated in [6]. Least-square fitting may give systematically biased estimates of the scaling parameters and most of the inferential assumptions for regression are violated. Moreover, in all these studies no statistical testing was mentioned to measure the goodness-of-fit of each fitted degree distribution. As an alternative, Maximum Likelihood Estimation (MLE) of the scaling parameters should be used and alternative distributions should be also tested. In [6] an analytical framework for performing such tests for power-law models is provided, which has been subsequently extended for testing other distributions as well. However, it has been noted that results are still very much dependent on the way the data is preprocessed, how the network is extracted, its dimensions and whether one uses region or voxel-based networks. For instance, Hayasaka et al. [10] found that, although degree distributions of all analysed functional networks followed an exponentially truncated model, the higher the resolution, the closer the distribution was to a power-law.

In this work, as a first step to address this issue we analyzed the resting-state fMRI (rs-fMRI) data of 10 subjects obtained from the Human Connectome Project database. Using the MLE method, advocated in [6], we estimate the scaling parameters for the best possible fit for a model distribution and then check the goodness-of-fit for this distribution by comparing it to synthetic generated data. We do this for four model distributions: power-law, exponential, log-normal and generalized Pareto. The reason for choosing the generalized Pareto model is due to the fact that it interpolates between fat-tail and thin-tail distributions, including the power-law and exponential as special cases. In what follows, we find that at a resolution of 10K nodes, the statistics favor the generalized Pareto thin-tail distributions.

2 Materials and Methods

2.1 Subjects, imaging data and network extraction

High-quality, high-resolution resting state fMRI scans of 10 subjects from the Human Connectome Project (Q1 data released by the WU-Minn HCP consortium in March 2013 [13]) were analysed in this study (age range: 26–35, 16.7% male). Individual rs-fMRI data were acquired for $\sim$15 minutes providing a total of $\sim$98,300 grayordinates time–series of 1,200 time points each. A schematic illustration of the process used to build the networks is provided in figure 1A. Building and visualizing functional networks was done using the BrainX3 platform [3, 5, 4]. For all the subjects, the original data-set was downsampled to $\sim$10,000 nodes by averaging the time-series of neighbouring grayordinates within a 5 mm3 cube.

Pearson’s correlation coefficients were calculated between each possible pair of nodes to build a $N\times N$ functional connectivity matrix, which is symmetric by construction and with self-connections set to zero. The matrix was then thresholded to derive weighted undirected adjacency matrices. We examined a range of 18 different thresholds (R) between $-0.7$ and $0.8$, at $0.1$ steps. Outside this range, the functional matrices become too sparse for meaningful analysis. For a positive threshold, each entry in the correlation matrix is set to 0 if its value is less than the threshold value and maintains its value otherwise. For a negative threshold, absolute values of the entries less than the threshold are maintained, while others are set to 0. In a weighted network, the weighted degree of a node is defined as the sum of all weighted edges connected to that node. Figure 1B illustrates the degree distributions of extracted networks across three different thresholds for a representative subject and averaged over all 10 data-sets.

2.2 Fitting parametric models to weighted degree networks

For every network generated from subject data, the vector of degrees $\textbf{x}=[x_{1},x_{2},...,x_{n}]$ is sorted in ascending order for each threshold. For every $x_{i}$, following [6], we use the method of maximum likelihood to estimate the scaling parameter $\alpha$ providing the best possible fit for the hypothetical power-law distribution $P(x)\sim Cx^{-\alpha}$ for the tail of the observed data in the range $x_{i}$ to $x_{n}$. Next, we calculate the Kolmogorov-Smirnov (KS) statistic for this power-law distribution with respect to $x_{i}$. Out of all possible $x_{i}$ from the data, the one with the smallest KS statistic corresponds to the lower bound $x_{min}$ for power-law behavior in the data. The next step is to verify whether this is indeed a good fit for the data. For that, a large number of synthetic data-sets are sampled from a true power-law distribution with the same scaling parameter $\alpha$ and bound $x_{min}$ as the ones estimated for the best fit of the empirical data. We fit each synthetic data-set to its own power-law model and calculate the KS statistic for each one relative to its own model. An empirical p-value is then calculated as the fraction of the time the empirical distribution outperforms the synthetically generated ones (by having a smaller KS statistic value). If p-value $\leqslant 0.1$, the power-law hypothesis can be ruled out as a non plausible explanation of the data. Nevertheless, large p-values do not guarantee that the power-law is the best model and the power-law fit has to be compared to a class of competing distributions.

3 Results

Power-law testing was performed on Matlab (Mathworks Inc., USA) using [6]. Further, for testing exponential, log-normal and generalized Pareto models we adapted the framework provided in [6] to include these competing hypothesis. For each subject, we analyzed thresholds in the range $-0.7$ to $0.8$, with $0.1$ increments. The parametric goodness-of-fit test was conducted over 1,000 repetitions, ensuring precision of p-value up to two decimal digits.

Our results are summarized in tables 1 and 2, respectively. An hypothesis is considered plausible if the p-value is larger than 0.1. Averaging over subjects, the p-values indicate that the power law hypothesis is rejected in the 83.3$\%$ of the analyzed thresholds. Instead, 61.1$\%$ of the examined thresholds are consistent with a generalized Pareto hypothesis, 55.6$\%$ with a log-normal hypothesis and in 33.3$\%$ of the cases with the exponential hypothesis, with several of these thresholds passing multiple tests. Median p-values are consistently larger for the generalized Pareto hypothesis (figure 2).

For each threshold examined, we then perform log-likelihood ratio tests to check which one among the consistent models is the most plausible in describing the empirical data. For all the positive thresholds up to 0.7 the evidence strongly goes in favor of the generalized Pareto distribution. Overall, the generalized Pareto model is outperforming the other candidate models in 41$\%$ of the examined cases (all subjects and all thresholds). In a 13$\%$ of the comparisons the log-normal distribution resulted in a better fit, 3$\%$ were better fitted by an exponential model, 2$\%$ by a power law, whereas the remaining 41$\%$ could not be explained by any model (due to insufficient data points and extreme thresholds).

For several positive thresholds the $k$ parameter is equal or close to zero, thus approaching an exponential distribution, whereas for other thresholds in the positive range, the generalized Pareto model passes with negative $k$, meaning a suppressed tail (figure 3).

4 Discussion

In this study we sought to systematically analyze scaling properties of human brain functional networks in the resting state, obtained from high-resolution fMRI data. We constructed networks of 10,000 nodes. Our analysis took into account actual weighted degree distributions from the data and we scanned through the full range of positive as well as negative correlation thresholds. For model selection, we imposed a criterion of p-value $>0.1$ and we conducted a log-likelihood ratio test among the different hypothesis.

We have shown that the degree distribution of the nodes does not follow a scale-free topology, as reported in [7]. The power law hypothesis is strongly rejected in the majority of the thresholds we examined. Indeed, it is the generalized Pareto distribution that is consistently preferable to the competing models for most of the thresholds.

These results suggest that after taking into account continuously weighted rather than binary networks, the dynamics of brain functional networks might not be governed by as many ultra-high degree hubs as a typical scale-free network might suggest. This bodes well for real brain networks when considering resilience to attacks, compared to their scale-free counterparts. For future work, we intend to test whether these distributions hold for different network resolutions and parcellations. Moreover, it would be interesting to see how these results compare to the “core-periphery” organization of brain structural networks [9, 12], which shows a preference for a distributed core, rather than few ultra-high degree hubs.

Acknowledgments.

The work has been supported by the European Research Council under the EU’s 7th Framework Programme (FP7/2007-2013)/ERC grant agreement no. 341196 to P. Verschure. Data were provided [in part] 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.