Introduction

The intrinsic functional architecture of the brain and its changes due to cognitive engagement, ageing and diseases are nodal topics in neuroscience, attracting considerable attention from many disciplines of scientific investigation. Functional neuroimaging and electrophysiological monitoring provide multiscale and multilevel signatures of brain intrinsic functional connectomics, giving access to amount of data previously inconceivable and nowadays tricky to examine. In this framework complex network theory1,2,3, which represents the state-of-the-art of big data multivariate analytics, demonstrated to allow a description of both the structural and the functional organization of local cortical and subcortical mesoscopic interactions, returning a repertoire of unexpected properties of brain connectomics4. Accordingly the functional architecture of the brain has been described in terms of patterns of communication between brain regions, treated as evolving networks, whose dynamics was found to correlate with behavioral outcomes5. Brain networks are characterized by a balance of dense relationships among areas highly engaged in processing roles, as well as sparser relationships between regions belonging to systems with different processing assignments. This segregation facilitates communication among brain areas that may be distributed anatomically but are recruited for processing specific tasks6. Along with that the integrated functional organization of the brain involves each network component executing a discrete cognitive function mostly autonomously or with the support of other components where the computational load in one is not heavily influenced by processing in the others7. In this quest for a constantly improving quantitative description of the brain and of the cardinal features of its functioning complex networks play a cr ucial role8,9,10,11. Specifically, it proved to be able to elicit both the scaffold of the mutual interactions among different areas in healthy brains4, 12 and the local failure of the global functioning in diseased brains13,14,15. The use of complex network theory passed progressively from an initial assessment of basic topological properties16,17,18 to a more sophisticated description of global features of the brain, as small-worldness19, rich-club organization10 and topology20, 21. However, a clear understanding of the functional advantage of a network organization in the brain, the characterization of its substrate and a description of brain regions aggregation as a function of the level of their interaction are still missing. In this paper, we show, for the first time, the hierarchical participation of brain regions to the whole network as a function of the weights assigned to the links among them, as well as the existence of a linear backbone of the human resting state functional network. Here the weight of a link between two regions represents the strength of correlation between the time-series of spontaneous neural activity within the two regions, as measured by blood oxygen-level-dependent (BOLD) functional MRI (fMRI)22. Specifically we interpret a functional connection between two nodes as the magnitude of the synchrony of their low-frequency oscillations, which is associated with the modulation of activity in one node by another one23, 24 largely constrained by anatomical connectivity25, 26. We used a percolation analysis21, 27 to highlight the progressive engagement of brain regions in the whole network as a function of the connectivity strength. Subsequently, by means of the maximum spanning forest (MSF) and the maximum spanning tree (MST) representations1, 28 we obtained the basal scaffold of the brain network, whose structure is characterized by a linear backbone in which few nodes (cerebellar and occipital regions) play a central role, even progressively increasing the spatial resolution or reducing nodes size.

Results

Inter-subject variability

We explore the inter-subject variability with several approaches. We compute the standard deviation of pairwise correlation values across the 40 individuals and report it together with the average correlation matrix used for the analysis. Figure 1(a) shows that there are few negative correlation values, however small in absolute term. The standard deviation ranges in the interval [0.039, 0.55], with most of the highest values concerning correlation averages close to zero (Fig. 1(b)). Those are few links which show positive and negative values across the 40 individuals and average around zero. They mainly belong to the Cerebellum (see SI) and the inter-subjects variability in this case could be due to the specific anatomical location with consequences on the resulting fMRI. However, they do not affect our results as in our analysis only the greatest correlation values play a key-role (percolation, MSF, MST), while the smallest threshold considered in the hierarchical integration of brain areas is 0.33. In Fig. 1(c) we also report the coefficient of variation for the correlation matrix thresholded using the False Discovery Rate (FDR). Most of the values concentrate around zero, while the greatest ones correspond to correlation values close to zero.

Figure 1
figure 1

Inter-subjects variability. (a) Mean, (b) standard deviation and (c) coefficient of variation of pairwise correlation values across the 40 subjects involved in the study.

As next step, we compute the pairwise cosine similarity between the correlation matrices of subjects. We find high similarity between subjects (Fig. S1(a)), with few exceptions (mainly, subjects 35 and 38). However, their presence does not affect the average correlation matrix as explicitly shown by applying the Jackknife resampling technique. One at a time, we remove the correlation matrix corresponding to a subject from the sample and average the remaining 39 correlation matrices. Then, we compute the cosine similarity between the 40 resulting correlation matrices (Fig. S1(b)). The similarity values belong to [0.9996, 0.9999], revealing that the emerged diversity for some subjects is negligible when we compute the average correlation matrix. For this reason, in the rest of the paper we will focus on the average correlation matrix computed over the 40 subjects in the analysis.

Percolation Analysis

The representative human functional brain network is often analyzed introducing specific thresholds to map the fully-connected correlation matrix in a sparse binary matrix4. Here, to avoid any arbitrary assumption we perform a percolation analysis21, 27 on the whole network.

We rank correlations in increasing order. One at a time we therefore remove from the network the link corresponding to the actual correlation value in the list, then we delete it from the rank. Step by step, we explore the global organization of the network after link removals. Specifically, in Fig. 2(a) we show the number of connected components updated each time a new link is removed.

Figure 2
figure 2

Percolation analysis. (a) Distribution of the number of network components versus the threshold applied to the correlation matrix. (b) Distribution of plateaux length observed in Fig. 1(a) and computed as the increment of correlation value between two newborn components. All figures show the percolation curves for the real case (blue) and the ensemble of its randomizations (green). Moreover figure (b) reports upper confidence intervals: mean plus respectively, one (red), two (orange), three (light blue) and four (bleak) standard deviations.

When by going down in the list of correlation values, the removal of several links does not lead to the increase of the number of components, we have a flat horizontal pattern defined plateau (stable state). The emergence and the number of plateaux shed light on the hierarchical structure of the network, their length unveils the intrinsic stability of certain network configurations after the removal of links. For the human functional brain network a remarkable hierarchy in the disaggregation process emerges from the comparison with a proper null model. The same percolation analysis is performed on the ensemble of 100 randomizations (Methods) of the observed correlation matrix showing a faster disaggregation in disconnected components with plateaux absent or of small length.

We compute the distribution of plateaux length looking at the increment of correlation values when the network passes from n to n + 1 components and show it in Fig. 2(b). In the same figure we also report the average plateaux length computed on the ensemble of randomizations. A great variability characterizes the percolation curve of the real network with significant deviations of the plateaux length from the random case.

Chain-like modules

The percolation analysis highlights the existence of a not trivial functional organization of the human brain network. Here, we investigate it considering a filtered network where for each area all weighted links but the strongest are discarded. This approach gives rise to 36 disconnected components forming a Maximum Spanning Forest (Methods) with a new information on link directionality. It simply indicates that each source points toward its maximally correlated brain area and it is not necessarily reciprocated. Figure 3 shows the abundance of modules with size 2 and reciprocated links: those are meanly mirror-areas of the right and left brain hemispheres or adjacent/very close ROIs belonging to the same hemisphere. Nodes are coulored according to the anatomical regions reported in Fig. 3(a) (detailed names of ROIs in Table 1 in SI). A noteworthy result concerns groups of size greater than 3 exhibiting a chain-like structure, sometimes very long as for the Cerebellum. This implies that most of the nodes in the MSF have in-degree equal to one, few equal to 2 and very few greater than 2. Furthermore, nodes tend to connect with nodes belonging to the same anatomical region. The only exception is represented by the Temporal Lobe: ROIs in this region are linked with all the other anatomical areas except for the Cerebellum and the Deep Grey Matter ones.

Figure 3
figure 3

Maximum Spanning Forest of the real human functional brain network. (a) Anatomical grouping of AAL parcellation of the human brain. Frontal Lobe; Insula; Cingulate; Temporal Lobe; Occipital Lobe; Parietal Lobe; Deep Grey Matter; Cerebellum. (b) Maximum spanning forest components. Arrows indicate that the ROI-source is maximally correlated with the ROI-target. Link thickness represents the strength of correlation. (c) Node in-degree distribution of the MSF. Comparison between the real brain network and the ensemble of its randomizations (100).

We build the MSF of the 100 randomizations of the real network obtaining a number of components varying in the range12, 24. All randomized correlation matrices exhibit a star-like organization of components in the MSF, while the number of modules of size 2 is dramatically reduced. Moreover, colors of linked nodes are completely random. In Figure S1 we show the MSF performed on one of such randomizations. We compute the distribution of node in-degree of the observed MSF and of the ensemble of its randomizations (Fig. 2(c)). In the real case the in-degree is small with almost the 70% of values equal to 1. The ensemble of randomizations shows more variability and a left-skewed distribution of in-degree, with almost the 50% of values equal to 0, the 30% equal to 1 and a maximum of 15. The outcome confirms that the MSF of all the randomized correlation matrices shows similar features: star-like organization of modules with few in-degree hubs representing the strongest connected partner for several brain areas. This result highlights the fundamental hierarchical organization of brain functional areas.

From forest to tree

As next step, we link together all groups of nodes in Fig. 3(b) such that (i) we add only one link between two modules previously disconnected and with the greatest weight; (ii) we end up with a unique connected component (Methods). The resulting network is a Maximum Spanning Tree as all nodes are connected through a maximal weighted path without forming loops.

The MST in Fig. 4(a) preserves the chain-like organization of nodes and mainly reproduces the anatomical division in regions. Some star-like aggregations emerge revealing the centrality of certain nodes as ROI 99, belonging to the Cerebellum and 82, part of the Temporal Lobe. Remarkable is the separation of the Deep Grey Matter, which represents the ancient part of the human brain, connecting to the periphery of the tree with the Temporal Lobe, the Insula and Cingulate areas. Nor chain-like organization neither relevant agglomeration of nodes are observable in the tree obtained from the MSF of the randomized correlation matrices (one examples in Fig. S2). Also in this case the distribution of node degree characterizes the two kinds of network organization: star-like (null model) and chain-like (real brain network). The comparison between the real network and the ensemble of its randomization (Fig. 3(b)) confirms the star-like organization of the random MST, with half of nodes having degree equal to 1 and maximum degree equal to 17.

Figure 4
figure 4

Maximum Spanning Tree of the real human functional brain network. (a) Colors represent anatomical regions according to the grouping of AAL parcellation in Fig. 3(a): Frontal Lobe; Insula; Cingulate; Temporal Lobe; Occipital Lobe; Parietal Lobe; Deep Grey Matter; Cerebellum. (b) Node degree distribution in the MST. Comparison between the real brain network and the ensemble of its randomizations (100).

Hierarchical integration of brain areas

The maximum spanning tree approach doesn’t provide any evidence about the internal connectivity of components observed in the MSF, as loops are not allowed. Here we show a number of snapshots of the network associated with specific percolation phases.

We proceed in the following way: (i) we consider the correlation values associated with such percolation states as thresholds for the real network; (ii) we discard correlation values below them; (ii) we inspect the network formation process according to the strength of edges.

The rationale behind this approach comes from the possibility to explore the level of segregation of brain areas and how they functionally and hierarchically integrate.

The appropriate percolation states and related thresholds are chosen among the stable states of the percolation curve (Fig. 2(a)) considering only plateaux whose length significantly deviates from the random case (null-model) in Fig. 2(b).

Figures 5, 6 and 7 show the human functional brain network related to the aforementioned thresholds. As the threshold value decreases (from Figs 5 to 7 and from panel (a) to (c)), the number of links increases showing how groups of nodes appear and their between and within connections intensify. The sequence of snapshots reveal the hierarchical functional organization of the brain network. Nodes do not connect forming several large and disconnected components, but: (i) nodes belonging to the Occipital Lobe start to connect each other and the anatomical region is almost completely connected after few steps; (ii) some nodes belonging to the same region form small chains; (iii) several nodes directly connect to the Occipital Lobe. Therefore, the most central area is represented by the Occipital Lobe, while ROIs belonging to the Deep Grey Matter lie in the periphery connecting to the network at the last moment in the percolation phases. The Cerebellum forms a quite long chain before joining the Occipital Lobe through the ROI 99 already emerged in the MST for its centrality.

Figure 5
figure 5

Snapshots of the human functional brain network. Each snapshot is realized thresholding the correlation matrix. Thresholds are chosen looking at the most significant deviation of plateaux length of the real percolation curve from its ensemble of randomization (Fig. 2(b)) Colors represent anatomical regions according to the grouping of AAL parcellation in Fig. 3(a): Frontal Lobe; Insula; Cingulate; Temporal Lobe; Occipital Lobe; Parietal Lobe; Deep Grey Matter; Cerebellum. Link thickness represents the strength of correlation.

Figure 6
figure 6

Snapshots of the human functional brain network. Each snapshot is realized thresholding the correlation matrix. Thresholds are chosen looking at the most significant deviations of plateaux length of the real percolation curve from its ensemble of randomizations (Fig. 2(b)). Colors represent anatomical regions according to the grouping of AAL parcellation in Fig. 3(a): Frontal Lobe; Insula; Cingulate; Temporal Lobe; Occipital Lobe; Parietal Lobe; Deep Grey Matter; Cerebellum. Link thickness represents the strength of correlation.

Figure 7
figure 7

Snapshots of the human functional brain network. Each snapshot is realized thresholding the correlation matrix. Thresholds are chosen looking at the most significant deviation of plateaux length of the real percolation curve from its ensemble of randomization (Fig. 2(b)). Colors represent anatomical regions according to the grouping of AAL parcellation in Fig. 3(a): Frontal Lobe; Insula; Cingulate; Temporal Lobe; Occipital Lobe; Parietal Lobe; Deep Grey Matter; Cerebellum. Link thickness represents the strength of correlation.

Discussion

In this paper we have reported the results of a network theory-based investigation of a large human fMRI dataset, aimed at assessing the basal architecture of the resting state functional connectivity network. We have computed the Maximal Spanning Forest (MSF) and the Maximal Spanning Tree (MST) for our population-wise brain. The former showed a peculiar composition of the basic modules composing the network, by revealing their structure as linear chains. The MSF reveals both intra- and inter-hemispherical modules, and the presence of small modules alongside with larger sub-networks. The MST, on the other hand, enables the classification of connectors and provincial areas. We also used a modified percolation analysis that retains the information on all the connected components of a given network. Such variation of the percolation analysis does not require the application of a threshold to obtain binary connectivity networks. By means of this technique, applicable to experimental correlation matrices, we showed hierarchically the progressive integration of brain regions in the whole organization of the connectivity network, which does not appear in a null model defined by constraining the spectrum of the observed correlation matrix. The percolation analysis here used represents a generalization of the classical one27 and is characterised by some improvements21 among which it must be recognized the ability to detect newly disconnected modules from secondary components.

Aggregation stages as a function of connectivity strength

As a complex system, the human brain is intrinsically organized into modules, creating a high degree of flexibility and adaptability of the system without fundamentally altering underlying structure29. Remarkably, our percolation analysis unearthed the existence of a hierarchical organization of areas, clusters of areas and modules within the human functional brain network. The process of network aggregation from isolated components appeared slow - compared to the null-model - with the emergence of significant stable configurations identified by several plateaux in the percolation curve. This outcome suggests a progressive engagement of brain regions according to the local amount of connectivity, namely a network organization in hierarchical modules whose participation to the entire structure is progressively included as the links strength decreases. The sequence of stable network configurations that, from the single areas, leads to the fully connected network, can be achieved by choosing values associated with specific plateaux of the percolation curve. Although this is technically topological modularity it is remarkable that it reflects functional and anatomical features. For instance bilateral correspondence of regions belonging to the same module was obtained without the introduction of symmetry constraints imposed in the analyses, supporting the conclusion that a hallmark of large-scale resting-state networks is their inter-hemispheric symmetry30. Interestingly, descending the hierarchical staircase of the percolation we found that occipital regions as the cuneus, the lingual gyrus, the occipital superior and middle gyrus and the calcarine sulcus tended to form a densely interconnected module. Along with that, frontal, cerebellar and temporal regions resulted involved in different modules characterised by linear chains, which joined the occipital cluster, in a progressive integration process, through a nodal bridge represented by the precuneus. In particular, temporal and occipital lobes resulted connected thanks to the bridging action of the fusiform and the inferior occipital gyrus as expected by considering the local cytoarchitecture31. Thalamus, putamen, pallidus, hippocampus, amygdala and caudate (here refereed to as deep grey matter) was found to form binary bilateral modules that joined the whole structure at the end of the integration process. If on one hand the remarkable connectivity within the occipital lobe can be ascribed to an efficient coupling via the interplay of tangential intracortical and callosal connections32, on the other hand such connectedness could be ascribed to fixational eye movements33. In support to this observation, it must be noticed that in Parkinson’s disease patients, which present with eye movement disturbances that accompany the cardinal motor symptoms, the worse the oculomotor performance is, the more the regional functional connectivity is decreased34.

The basal scaffold of the brain functional architecture

A second striking result has been obtained from MSF and MST analysis of the whole brain network. For the first time we report a non-trivial and ordered structure for the basal scaffold of the functional network of the brain. Keeping track only of the strongest connection for each region, the network appeared organised in components exhibiting a linear chain-like structure. When the same procedure has been applied to a proper randomisation of the real human brain functional network, the outcome strongly differed: hubs and star-like structure emerged. This result was tested for robustness by improving the number of regions included in the grey matter parcellation. We calculated correlation matrices from two further templates composed of 276 and 531 regions of interest, respectively, each obtained as a sub-parcellation of the AAL template. The chain-like organization has been proved to persist as the topological outcome of the MSF and the MST analyses of the human brain functional network, against the randomised counterpart that showed a star-like organization (Figures S3S14 in the Supplementary Information). Firstly, modules composition within the MSF deserves a comment. We found linear clusters of regions grouped according neuroanatomical criteria, except the temporal and insular regions. While it is expected to find close portions of the cortex strongly coupled, it is not obvious that the same amount of coupling encompasses regional bilaterality embedded in the modularity of the whole brain network. The exception represented by the temporal and insular regions can be ascribed to different reasons. The insular lobes showed to be more connected with frontal regions, through the Rolandic opercula, than with each other (each one respecting its own laterality). Concurrently, the fronto-insular cortex, wherein generation and experience of emotion are found35, is implicated in the elaborate circuitry associated with awareness36 and uniquely (together with the anterior cingulate cortex and the dorsolateral prefrontal cortex) characterized by the presence of Von Economo neurons37, 38. On the other hand the fragmentation of the temporal lobe can follow the mosaic structure of its functional specialization. In fact it includes auditory sensory and association cortices, part of the posterior language cortex, visual and higher order association cortices, primary and association olfactory cortices and enthorinal cortex39. The grouping of other heterogeneous clusters is due both to spatial closeness (precentral and postcentral gyrus) and to homofunctional specialization (hippocampus and parahippocampal gyrus). The MST representation of the functional brain network scaffold was found characterized by a long linear structure, whose skeleton follows the anatomical spatial distribution of lobes (from the frontal to the occipital lobe, passing through the temporal and parietal lobes), small chains linked alongside and a star-like structure (composed of cerebellar regions) at the end of one side. The main connectors along the main chain included: the cuneus/precuneus, the superior temporal gyrus, the superior occipital gyrus, the inferior parietal lobule, the medial orbitofrontal cortex/gyrus rectus and the VI lobule of the cerebellum. Notably, precuneus, inferior parietal cortex and medial orbitofrontal cortex constitute the Default Mode Network40, 41, which is a crucial module of the whole brain network both for its functioning and as the first target of metabolic alterations both in neurological and psychiatric diseases42, 43. The large connectivity with the Default Mode Network showed by the superior temporal gyrus, which confers to this region a central role in the MST, has been observed in fMRI data44 but never found in FDG-PET data. Since one of the differences between the two techniques is the level of noise produced during the recordings (relatively silent during FDG-PET vs intensively noisy during MRI), the absence of overlap of the patterns of connectivity reported in that region might be explained by the auditory engagement included in the MRI experiments45. The same argument works as far as the central role of the superior occipital cortex, being the visual activity engaged by the fixation during the resting-state fMRI experiment. A final comment must be reserved to the cerebellar module. The understanding of human cerebellar function has undergone a paradigm shift. No longer considered purely devoted to motor control, a wide role for the cerebellum in cognitive and affective functions is supported by anatomical, clinical and functional neuroimaging data46, 47. Recent evidence from functional connectivity studies in humans indicates that the cerebellum participates in functional networks with sensorimotor areas engaged in motor control and with association cortices that are involved in cognitive processes48,49,50. The evidence of the cerebellar engagement in higher cognitive functions as well as in sensorimotor processes might explain its central role in the basal scaffold of the whole functional brain network. Specifically, it has been reported that the lobule VI of the cerebellum, here found as a main connector of the network, has been found to be involved in different cognitive tasks, from mental rotation to working memory processing51. In addition, it has been showed that during the first stages of consciousness loss, induced by mild sedation, the pattern of thalamic functional disconnections included regions belonging to the cerebellum52, supporting the evidence of a strong cerebello-thalamic interaction showed by the MST analysis.

Choice of the similarity measure

Lastly it must be noticed that both in MSF and MST analysis, as weights of our graphs, the square of the correlation coefficient (r2) was used. Positive and negative correlations, characterized by the same absolute value, refer to the same amount of information in terms of connectivity strength. Since in this study we considered undirected weighted graphs, the phase-shift between time series, which gives rise to the sign of the correlation coefficient, was not significantly informative. Indeed53, showed that both positive and negative correlations are related to synchronized neural activity and that the sign of correlations can result from neural-mediated hemodynamic mechanisms, which lead to temporal and spatial heterogeneity. Moreover, r2 was chosen since it has a strict relationship with the coefficient of determination R2. Specifically, after a simple linear regression, R2 equals the square of the correlation coefficient between the observed and modelled data values54. In terms of a correlation analysis, r2 estimates how well a time series can predict another one and can be considered as a measure of the magnitude of the relationship between time series.

Conclusion

In conclusion we have reported that the whole architecture of the brain activity at rest is characterised by chain-like structures organized in a hierarchical modular arrangement. This kind of organization is thought to be stable under large-scale reconnection of substructure55, and efficient in terms of wiring costs56. In addition the chain-like feature of the basal modules facilitate the wiring and reconnection processes, not asking for a specific region as a hub, but functioning the whole module as a multiple access point for the functional plasticity.

Methods

Forty healthy subjects (mean age = 38, SD = 10; mean educational attainment = 15, SD = 3; male N = 19) participated in this study. All subjects were carefully screened for a current or past diagnosis of any DSM-5 Axis I or II disorder using the SCID-5 Research Version edition (SCID-5-RV57) and the SCID-5 Personality Disorders (SCID-5-PD58). Exclusion criteria included: (1) a history of psychoactive substance dependence or abuse during the last one year period evaluated by the structured interview SCID-5-RV, (2) a history of neurologic illness or brain injury, (3) major medical illnesses, that is, diabetes not stabilized, obstructive pulmonary disease or asthma, hematological/oncological disorders, B12 or folate deficiency as evidenced by blood concentrations below the lower normal limit, pernicious anemia, clinically significant and unstable active gastrointestinal, renal, hepatic, endocrine or cardiovascular system disease, newly treated hypothyroidism, (4) the presence of any brain abnormality and microvascular lesion apparent on conventional FLAIR-scans. The presence, severity and location of vascular lesions was computed using a semi-automated method59, (5) IQ below the normal range according to TIB (Test Intelligenza Breve, Italian analog of the National Adult Reading Test NART) (6) global cognitive deterioration according to a Mini-Mental State Examination (MMSE)60 score lower than 26, (7) dementia diagnosis according with DSM-5 criteria61 and (8) non-Italian language native speaker. All participants were right-handed, with normal or corrected-to-normal vision. They gave written informed consent to participate after the procedures had been fully explained. The study was approved and carried out in accordance with the guidelines of the IRCCS Santa Lucia Foundation Ethics Committee.

Data acquisition and preprocessing

FMRI data were collected using gradient-echo echo-planar imaging at 3T (Philips Achieva) using a (T2*)-weighted imaging sequence sensitive to blood oxygen level-dependent (BOLD) (TR = 3 s, TE = 30 ms, matrix = 80 × 80, FOV = 224 × 224, slice thickness = 3 mm, flip angle = 90, 50 slices, 240 vol). A thirty-two channel receive-only head coil was used. A high-resolution T1-weighted whole-brain structural scan was also acquired (1 × 1 × 1 mm voxels). Subjects were instructed to lay in the scanner at rest with eyes open. For the purposes of accounting for physiological variance in the time-series data, cardiac and respiratory cycles were recorded using the scanner’s built-in photoplethysmograph and a pneumatic chest belt, respectively. The human brain was segmented into 116 macro-regions from the AAL template62. Resting state fMRI signals were averaged across each region to generate 116 time-series, which in turn were pairwise correlated calculating the Pearson’s coefficient and organized in a symmetric matrix. Two further templates were used to test the robustness of results at increasing spatial resolution by subdividing AAL regions as described in ref. 63. The number of regions included in the two templates was 276 and 531 respectively. Several sources of physiological variance were removed from each individual subject’s time-series fMRI data. For each subject, physiological noise correction consisted of removal of time-locked cardiac and respiratory artifacts (two cardiac harmonics and two respiratory harmonics plus four interaction terms) using linear regression64 and of low-frequency respiratory and heart rate effects65,66,67. FMRI data were then preprocessed as follows: correction for head motion and slice-timing and removal of non-brain voxels (performed using FSL: FMRIB’s Software Library, www.fmrib.ox.ac.uk/fsl). Head motion estimation parameters were used to derive the frame-wise displacement (FD): time points with high FD (FD > 0.2 mm) were replaced through a least-squares spectral decomposition as described in ref. 68. Data were then demeaned, detrended and band-pass filtered in the frequency range 0.01–0.1 Hz, using custom software written in Matlab (The Math Works). For group analysis, a two-step registration process was performed. fMRI data were transformed first from functional space to individual subjects’ structural space using FLIRT (FMRIB’s Linear Registration Tool) and then non-linearly to a standard space (Montreal Neurological Institute MNI152 standard map) using Advanced Normalization Tools (ANTs; Penn Image Computing & Science Lab, http://www.picsl.upenn.edu/ANTS/). Finally data were spatially smoothed (5 × 5 × 5 mm full-width half-maximum Gaussian kernel). In order to create an average adjacency matrix at the population level, subject-wise matrices were Fisher-transformed, averaged across subjects and back-transformed.

Percolation analysis and thresholding

We used a variation of the percolation analysis proposed by ref. 27, recently developed by ref. 21: (i) we ranked all experimentally determined correlation coefficients in increasing order; (ii) one at a time we removed from the network the link corresponding to the actual correlation value in the list; (iii) then we deleted it from the rank. Finally, step by step we explored the global organization of the network after each link removal by considering the number of its connected components. This enabled us to calculate connected components by systematically varying the threshold on the network. Among the possible thresholds we looked at those values associated with the plateaux showed in the curve of Fig. 2(a) obtained by the percolation process. Specifically, we considered all correlation values corresponding to points in the distribution of the real plateaux length overcoming the mean plus four standard deviations of the distribution of plateaux length of the ensemble of randomizations. This procedure led to the definition of 23 values, for which stable configurations of brain networks were found.

Maximum spanning forest and maximum spanning tree

Starting from the initial correlation matrix, corresponding to the fully connected network, we built a subnetwork keeping only links associated with the maximum value of each row. In other words, for each node only the strongest link survives. The relation is not symmetric, indeed if w ij represents the strongest link for node i, this does not imply that w ji represents the strongest link for node j. We ended up with a directed network of exactly N links connecting each source-ROI with its maximally correlated target-ROI.

This approach is equivalent to create a network from scratch in three steps: (i) to rank in decreasing order all the pairwise correlation values; (ii) to look at the list from the top and adding a new link if at least one of the two involved nodes has degree zero, otherwise discarding such correlation; (iii) to update the node degrees. The procedure ends when all nodes have at least degree equal to one. In both cases, the resulting network is made up by several components of nodes connected through their strongest links and without forming loops. In other terms, we have a Maximum Spanning Forest (MSF) of the initial correlation network. Furthermore, the first procedure introduced a directionality, simply revealing for each brain area which is its maximally correlated counterpart.

As second step, we connected all components of the MSF to form a unique tree-like component, following this procedure: (i) we took the correlation values discarded during the construction of the MSF and ranked them in increasing order; (ii) starting from the top of the list, the link associated to the considered correlation value was added to the network if the two nodes involved did not belong to the same group and the two sets were merged into one; (iii) otherwise the correlation value was discarded. The procedure ended when all the nodes belonged to the same connected component. By construction the resulting subnetwork of the original network is a Maximum Spanning Tree (MST) as it does not contain loops and maximizes the total weight with the minimum number of links (2(N − 1)).

A statistical benchmark for human brains

In order to define to what extent the stepwise structure highlighted by the percolation analysis, the MSF and MSF were significant, we compared the results with a proper statistical benchmark. Specifically, in order to understand whether the results were a proper representation of the real intrinsic organization of the brain functioning at rest we needed to define an appropriate null model. A matrix has to satisfy some requirements in order to be a correlation matrix: it must be symmetric, with diagonal elements equal to 1, with off-diagonal elements in the range [−1, 1] and it has to be positive semidefinite. Not all standard randomization techniques are suitable in this case. Here, we kept fix the spectrum distribution of the observed correlation matrix applying a a random orthogonal similarity transformation to the diagonal matrix of eigenvalues and then a sequence of Givens rotations69.

This choice avoids to introduce further procedures to adjust reshuffled matrices in order to be positive semidefinite. We generated 100 randomizations of the observed correlation matrix. All comparisons between the ensemble of randomizations and the real network were performed on the pointwise squared matrix in order to neglect signs but to keep the ranking of correlation values. Relevant quantities were computed on each randomization, than averaged over the ensemble.

Data Availability

The data collected for the study and analysed in this work are available at the following link: http://networks.imtlucca.it/downloadable-files. The reader can find the time-series extracted from the AAL parcellation (116 ROIs) after that all the pre-processing steps were perfomed.