Modular co-organization of functional connectivity and scale-free dynamics in the human brain

The human connectome is a comprehensive map of how brain regions are mutually connected and is fundamentally important for understanding neuronal communication and brain dynamics in both health and disease (Fornito, Zalesky, & Breakspear, 2015). Although functional connectivity (FC) in the human brain was initially characterized as a stationary network (Raichle, 2009), recent empirical (Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014) and theoretical (Deco, Jirsa, & McIntosh, 2013) studies have complemented this view by revealing considerable connectivity pattern fluctuations that are systematic beyond simply representing noise in weak stationary connectivity (Hansen, Battaglia, Spiegler, Deco, & Jirsa, 2015; Mitra, Snyder, Blazey, & Raichle, 2015). Such dynamic connectivity has been observed at a range of temporal scales from milliseconds to minutes (Larson-Prior et al., 2013). It has, nevertheless, remained unclear how these static and dynamic patterns of FC are related to the actual dynamical states of the brain and to the statistical patterns of activity propagation associated with this scale-free dynamics (Kopell, Gritton, Whittington, & Kramer, 2014).

The ubiquity of scale-free dynamics in neuronal activity suggests that neuronal systems could operate at a critical state (Chialvo, 2010). Critical neuronal dynamics are characterized, for instance, by the power-law scaling of neuronal activity avalanches at fast (10⁻³ to 10⁻¹ s) timescales (Beggs & Plenz, 2003), as well as by power-law long-range temporal correlations (LRTCs) in neuronal fluctuations at slow (10¹ to 10³ s) timescales (Linkenkaer-Hansen, Nikouline, Palva, & Ilmoniemi, 2001). However, the hierarchically modular rather than homogeneous structural connectivity of the brain may expand our concept of criticality from the system operating near a critical point into the system operating in an extended critical region or, more specifically, in a Griffiths phase (Hilgetag & Hutt, 2014; Moretti & Munoz, 2013). On the other hand, it is important to note that several kinds of noncritical processes are also associated with power-law LRTCs; hence, observations of scale-free dynamics do not, per se, indicate criticality (see, e.g., E. J. Friedman & Landsberg, 2013; Schwab, Nemenman, & Mehta, 2014).

Whether scale-free brain dynamics and/or neuronal criticality are related to FC has received attention only recently (Ciuciu, Abry, & He, 2014). Avalanche-like large-amplitude propagating patterns in BOLD signals have a connectivity structure similar to that of FC (Tagliazucchi, Balenzuela, Fraiman, & Chialvo, 2012). In the same vein, slow amplitude fluctuations of neuronal oscillations are characterized by both LRTCs (Linkenkaer-Hansen et al., 2001) and long-range interareal correlations with neuroanatomical structures similar to those of fMRI resting-state networks (Brookes et al., 2011). Hence, these converging, albeit indirect, relationships suggest that critical dynamics and the functional connectome could be linked, but in a manner that hitherto has been unclear. Recent computational-modeling studies have further elucidated this link by suggesting, on the one hand, that the modular network architecture of resting-state FC emerges specifically when the system is poised at a critical state (Haimovici, Tagliazucchi, Balenzuela, & Chialvo, 2013) and, on the other, that avalanches in a modular system propagate preferentially within connectivity modules (Russo, Herrmann, & de Arcangelis, 2014).

We advance here two hypotheses. First, we posit that neuronal avalanches propagate preferentially along the predominant pathways of FC and within the modules of the functional connectome. Second, we suggest that these modules are also endowed with partially independent dynamic states, which are reflected in stronger interareal correlations of LRTC exponents within the modules than between them. To test these hypotheses, we assessed the relationship between scale-free dynamics and FC by using intracranial stereotactic electroencephalography (SEEG) and magnetoencephalography (MEG). To ensure that the hypotheses were testable, we also performed several lines of control analyses to establish that the indices of dynamics and connectivity could be estimated independently. We characterized scale-free dynamics by connectomes of avalanche propagation patterns and interareal correlations in the LRTC scaling exponents, and compared the connection strengths and modular organizations in these connectomes with those in connectomes of FC measured by the phase synchronization (Vinck, Oostenveld, van Wingerden, Battaglia, & Pennartz, 2011) and amplitude correlations (Brookes, Woolrich, & Barnes, 2012) of ongoing neuronal oscillations. We show that both neural dynamics and FC have, in a frequency-specific and nontrivial manner, similar constellations of strongest connections and shared modular structures, which provides empirical evidence that neuronal connectivity and scale-free dynamics are intimately linked.

In our MEG and SEEG data, a neuronal avalanche was defined as a set of suprathreshold peaks in waves of broadband (1–120 Hz) ongoing activity occurring in consecutive time bins (Figure 1A). The statistical properties of avalanches can be described by the distributions of their sizes—that is, the total numbers of peaks in each avalanche—and lifetimes—that is, the durations of the avalanches. In this study we used time bin widths of 8 and 16 ms for the MEG and SEEG data, respectively, and thresholds of two, three, and four standard deviations (SDs) to characterize the propagation of neuronal avalanches in different scaling regimes (Zhigalov, Arnulfo, Nobili, Palva, & Palva, 2015). To evaluate the overall statistical nature of these avalanches, we first assessed the goodness of fit between the size distribution data and a truncated power-law model by computing the Kolmogorov–Smirnov distance (D_KS). The largest D_KS for the different thresholds was 0.04 (p > 0.85) suggesting that the distributions were approximated well by the truncated power-law function (Figure 1B). The scaling exponents of the avalanche sizes were estimated using the maximum-likelihood approach (see Materials and Methods). We found similar mean power-law scaling exponents in the MEG and SEEG data, with values close to 1.5 at the three-SD threshold (Figure 1C). We also assessed the avalanche branching ratios at different thresholds and observed values close to 1 at three SDs in both the MEG and SEEG data. The scaling exponent of 1.5 and branching ratio of 1 (Figure 1C) correspond to a critical branching process and are well in line with a body of prior observations on neuronal avalanche dynamics (D. Plenz, 2012).

To corroborate the truncated power-law fitting of the distributions of avalanche sizes, we applied a likelihood ratio test to compare the truncated power-law model with an exponential model. The results showed that the truncated power-law model outperformed the exponential model (p < 0.001, log-likelihood ratio test) for all subjects and thresholds.

The distributions of the avalanche lifetimes were biased by the numbers of samples and were less robust than the size distributions in individual subjects. To assess the lifetime distributions, we thus concatenated avalanches across subjects. Consistent with the prior experimental and theoretical literature, the power-law exponents of the avalanche lifetimes were near 2 at the three-SD threshold, where a size distribution exponent of 1.5 was found (Zhigalov, Arnulfo, Nobili, Palva, & Palva, 2017, Figure S1).

Finally, we assessed the scaling function for the shapes of the neuronal avalanches (Figure 1D), by averaging avalanches within ranges of lifetimes (Figure 1D, upper panel) and then rescaling the functions by normalizing the lifetimes (Figure 1D, lower panel). At rescaling, a collapse of the avalanche shapes was observed, which is again well in line with previous studies (Beggs & Timme, 2012; N. Friedman et al., 2012) and further supports the hypothesis that the neuronal avalanches observed here arose from a critical process.

We characterized the spatial-propagation patterns of neuronal avalanches by compiling a connectome of the empirical probabilities with which the activity peaks of the cortical parcels in the first time bin were followed by activity peaks of the parcels in the second time bin (Figure 1E and Zhigalov et al., 2017, Figure S2; see Materials and Methods).

We hypothesized that the neuronal avalanches observed in MEG and SEEG (Figure 2A) would propagate preferentially between brain areas coupled by FC. To test this hypothesis, we measured electrophysiological FC by means of the pairwise phase synchronization and amplitude correlations (Figures 2B and 2C) of narrow-band neuronal oscillations. Furthermore, to assess the interareal relationships of scale-free dynamics per se—that is, to measure how local dynamical states were correlated among all cortical regions—we obtained LRTC connectomes (Figure 2D) by correlating across subjects the LRTC scaling exponents of orthogonalized amplitude time series for all pairs of parcels. The LRTC exponents are proportional to the system’s proximity to the critical point (Poil, Hardstone, Mansvelder, & Linkenkaer-Hansen, 2012), and hence, if brains behaved like a “single” near-critical system, this connectome should be uniform. However, already a visual inspection of the LRTC connectome (see Figure 2D) showed that, like the other connectomes, it was highly nonhomogeneous, and although many brain areas were significantly correlated in both MEG (r > 0.54, p < 0.05, Pearson’s correlation test) and SEEG (r > 0.39, p < 0.05, Pearson’s correlation test), also a considerable fraction of regions were below these nominal significance thresholds. Hence, at least among some subsets of cortical areas, the local dynamics were only weakly correlated.

To test the hypothesis of colocalized avalanche propagation pathways and FC, we first measured the connection-level similarity—that is, the edge similarity—of avalanche propagation with phase synchronization and amplitude correlations by means of the Pearson correlations of these connectomes (see Materials and Methods). For phase-correlation-based FC, the results showed that the edge similarity was highly significant (well above the 99.9% confidence intervals of the surrogate data, corresponding to uncorrected p < 0.001) between avalanche propagation at three SDs and phase synchronization. This relationship was found in a limited range of frequencies around the α and β bands (i.e., 8–14 Hz and 15–30 Hz, respectively; Figure 3A, green lines), with peaks at around 16 Hz in MEG and at around 32 Hz in SEEG. The similarity between phase correlations and avalanche propagation at both two and four SDs was much weaker (p < 0.001, Wilcoxon test) in both the α and β bands than at three SDs, and significant only at four SDs in the β band (Figure 3A, blue lines) in MEG. The edge similarity between avalanche propagation and the phase connectomes was also highly significant in SEEG (Figure 3A), but the differences between thresholds were not significant in SEEG (p > 0.26, Wilcoxon test).

Similar to the phase-correlation FC, amplitude-correlation-based FC was also significantly correlated with avalanche propagation at three SDs for frequencies up to ∼30 Hz in MEG and up to ∼60 Hz in SEEG (Figure 3B, green lines). The spectral profiles of these similarities were, however, distinct from those for phase FC in extending to the lowest frequencies (4−8 Hz), even clearly peaking at 8 Hz (in MEG). This suggests that the amplitudes but not the phase couplings of 4- to 8-Hz oscillations are associated with the propagation pathways of broadband avalanches. The similarity was again significantly greater for avalanche propagation at three SDs than at two or four SDs (Figure 3B) across the entire range of frequencies in MEG (p < 0.001, Wilcoxon test), whereas these differences in SEEG were not significant. These data thus suggest that neuronal avalanches indeed propagate preferentially between brain areas that are coupled by phase synchronization and/or amplitude correlations.

To corroborate that the avalanche propagation and phase or amplitude connectomes can be estimated independently, and hence to ensure that the hypotheses tested above were falsifiable, we performed two types of control analyses (see Materials and Methods). In the first approach, we deleted from the narrow-band data all data segments in which the broadband amplitude exceeded the avalanche detection threshold—that is, all data used in the avalanche analyses—and recomputed the phase and amplitude connectomes. Evaluation of their similarity with the original avalanche propagation connectome revealed high similarity that was not significantly different from what we had obtained with the uncut data (p > 0.99, KS test; Figure 4A, solid lines). Reversing this analysis, we estimated the phase-synchronization and amplitude-correlation connectomes by only using the suprathreshold data. Here, the similarity of the suprathreshold functional connectomes with the original avalanche propagation connectome (Figure 4A, dashed lines) was at chance level. These analyses thus show that the similarity of FC and avalanche propagation does not arise from temporally colocalized time windows in the data: avalanches neither bias the quantification of FC nor are trivially produced by concurrent FC.

To corroborate these insights, in the second approach, we swapped all the time windows of suprathreshold avalanche segments between random parcel pairs, and thereby systematically reordered the avalanche propagation connectome. We then refiltered the data and recomputing the phase-synchronization and amplitude-correlation connectomes. Their similarity with the original avalanche propagation connectome did not decrease significantly much below the original values (p > 0.99, KS test; Figure 4B, solid lines). On the other hand, the reordered avalanche propagation connectome was uncorrelated with the functional connectomes estimated from the same reordered data (Figure 4B, dashed lines). Although this result was not surprising, in the sense that the original correlation of FC and avalanches was broken by reorganizing the avalanches, this analysis confirmed the first one and showed conclusively that the bias that the avalanche segments impose on FC estimates from the same data is negligible, and hence that the avalanches and FC can be estimated independently.

These data, together with the prior findings of (i) avalanche–FC correlations being distinct for amplitude and for phase FC; (ii) avalanche–FC correlations being spectrally limited, whereas FC is not (see below and in Figure 6); and (iii) avalanches at two SDs not being more similar to FC than at three SDs, thus strongly suggest that the broadband avalanches are neuro physiologically and phenomenologically distinct from the multitude of phase-/amplitude-coupled narrow-band oscillations.

We then assessed how the interareal relationships in local LRTCs were correlated with the patterns of FC. To this end, we first estimated the edge similarity between the LRTC and phase/amplitude connectomes and found robustly significant, but again band-limited, similarities (Figures 5A and 5B) with spectral profiles akin to those found for avalanche propagation (see Figures 3A and3B). Finally, the comparison of avalanche propagation at three SDs with the LRTC connectomes also revealed significant similarity in both the MEG and the SEEG data (Figure 5C). The similarity was significantly larger for the avalanche propagation at three than at two and four SDs in the MEG data (p < 0.001, Wilcoxon test). Moreover, the phase-synchronization and amplitude-correlation connectomes were clearly correlated in both MEG and SEEG (Figure 5D), with 16 and 32 Hz as the modal frequencies, respectively. Throughout these similarity estimates, and also in the mutual correlation of phase and amplitude FC, the frequency band of significant correlations extended to ∼60 Hz in SEEG and only to ∼30 Hz in MEG, which is likely attributable to the approximately millimeter versus approximately centimeter spatial scales across which SEEG and MEG measure coherent neuronal population activity.

These findings together thus show that scale-free brain dynamics, including both power-law-scaled, millisecond-range activity avalanches and LRTCs spanning hundreds of seconds, are intimately linked with the architecture of FC in human cortex.

The findings so far converged on the idea that two facets of critical dynamics, avalanche propagation and interareal coupling of local LRTCs, were colocalized with phase- and amplitude-correlation-based functional connectivity. To corroborate this finding with an independent index of critical dynamics, we estimated LRTCs in the momentary dynamics of phase synchrony with phase detrended fluctuation analysis (DFA; Botcharova, Farmer, & Berthouze, 2014) between all pairs of cortical parcels in the MEG data. The adjacency matrices of the phase DFA LRTC scaling exponents were thus directly comparable with those of FC. We found phase DFA and FC to be highly significantly similar for both phase and amplitude FC (well above the 99.9% confidence intervals) in the α and β frequency bands and to have profiles similar to those of the avalanche propagation and functional connectomes (Figure 3). These observations show that temporal correlations and spatiotemporal interactions are linked at different timescales, which consolidates the prior findings (see Figures 4 and 5) and overall suggests that the underlying system is poised at criticality.

The association of scale-free dynamics and FC was found in a relatively narrow frequency range around the α and β frequency bands. Does this reflect a “special” and nontrivial emergent phenomenon or a trivial “by-product” of the scale-free dynamics being automatically associated with connectivity?

To address this question, we asked whether either FC or LRTCs were correspondingly more prevalent in this α/β range than in the surrounding bands. We first assessed the connection strengths of the phase-, amplitude-, and LRTC-correlation connectomes across the frequency range. For the phase and amplitude connectomes, the mean correlation strengths exhibited subtle peaks, at 9 Hz in MEG and at 7 Hz in SEEG, but overall they decreased monotonically as a function of frequency (Figures 6A and 6B) and remained well above the 99.9% confidence limits of the corresponding surrogate values throughout the studied frequency range. Both phase- and amplitude-based forms of FC thus appear robust for frequencies between 4 and 128 Hz. Similar averaging of the local LRTC scaling exponents (Figure 6C) showed that the LRTCs also peaked at around 10 Hz in MEG and at 8 Hz in SEEG, but again were well above the surrogate confidence limits from 4 to 128 Hz. Significant scale-free dynamics are thus not unique to any single frequency band. In the same vein, the interareal correlations of LRTC exponents were essentially constant and highly significant across the studied frequency range (Figure 6D). Hence, just as FC was a salient characteristic of all frequency bands, also scale-free dynamics and the interareal structure of LRTCs characterized the 8- to 32-Hz range over which most of the salient dynamics–FC coupling was found, as well as the 4- to 8-Hz and 32- to 128-Hz ranges over which the dynamics–FC correlations were much weaker or insignificant. Together with the previous analyses (see Figure 5), these results thus show that significant FC and scale-free dynamics in the human brain can be mutually coupled or uncoupled in a frequency-dependent manner, which indicates that this coupling is nontrivial and neither automatic nor epiphenomenal.

The similarity between avalanche propagation pathways, interareal correlations of LRTCs, and functional connectomes indicated a significant connection-level colocalization of these networks. We addressed next whether scale-free dynamics and FC had similar modular structures in the MEG data. The SEEG data were excluded from this analysis because their sparse spatial coverage did not allow for robust detection of the functional modules.

We identified network modules, subgraphs, by using agglomerative hierarchical clustering of the connectomes of scale-free dynamics and FC (see Figure 2), and then measured the similarity of these modular structures by quantifying the overlap among the subgraph assignments (see Materials and Methods). The results showed that the modular structures of both avalanche propagation at three SDs (Figure 7A) and the interareal correlations of LRTCs (Figure 7B) had highly significant subgraph similarities with both phase- and amplitude-based FC. This colocalization of network modules was found in a limited frequency range up to 40 Hz, with the greatest values again being between 8 and 30 Hz (Figures 7A and 7B). For completeness, we also measured the subgraph similarity between avalanche propagation at three SDs and LRTCs (Figure 7C), as well as between the levels of phase and amplitude FC (Figure 7D), and found that both were robust in the 8- to 32-Hz range. In this frequency range, scale-free dynamics and FC thus had both internally and mutually shared modular structures.

To identify the cortical regions that composed the shared functional modules in these connectomes, we used consensus mapping of the parcel-subgraph neighborhoods obtained with the hierarchical clustering (see Materials and Methods). This analysis showed that the most consistently shared modules were localized to sensorimotor, visual, temporal, and prefrontal areas (ps < 0.05, permutation test; see the cyan, blue, green, and red modules, respectively, in Figures 7E–7H).

Criticality governs the spatiotemporal dynamics of numerous natural and artificial systems composed of large numbers of nonlinearly interacting elements. At the critical operating point, the system exhibits long-range correlations with power-law decay concurrently across all temporal and spatial scales (Bak, Tang, & Wiesenfeld, 1988; Chialvo, 2010). Brains have been proposed to operate near the critical point (Beggs & Plenz, 2003; Linkenkaer-Hansen et al., 2001) or in an extended critical region (Hilgetag & Hutt, 2014; Moretti & Munoz, 2013), which has fundamental implications for neurophysiological phenomena and neuronal information processing across the scales from neurons to brain systems. Neuronal criticality has primarily been operationalized with two hallmark phenomena, millisecond-range neuronal avalanches (D. Plenz, 2012) and long-range temporal correlations spanning hundreds of seconds (Linkenkaer-Hansen et al., 2001), although power-law characteristics have also been reported for several other forms of spatiotemporal brain dynamics (Kitzbichler, Smith, Christensen, & Bullmore, 2009). Conversely, although scale-free dynamics are an indisputable characteristic of brain activity, whether they arise from criticality or other mechanisms has remained a matter of debate (Touboul & Destexhe, 2010).

Much of our understanding of spatiotemporal brain dynamics has arisen from studies of hemodynamic and electrophysiological FC, which have shown that hierarchically modular architectures of neuronal interactions govern both resting-state brain activity (Mitra et al., 2015) and ongoing activity during task performance (de Pasquale et al., 2012). It is theoretically conceivable and predicted by computational models (Haimovici et al., 2013) that neuronal criticality and FC—concurrent, yet distinct and dissociable, network phenomena—would be related. However, very few prior experimental studies have aimed at elucidating this relationship (Pajevic & Plenz, 2009; Rubinov, Sporns, Thivierge, & Breakspear, 2011).

In this study, we sought to bridge this gap by quantifying the interareal relationships of scale-free dynamics and FC. Toward this end, we used both meso- and macro-scale electrophysiological recordings of human resting-state brain activity, with invasive SEEG (Arnulfo, Hirvonen, Nobili, Palva, & Palva, 2015) and noninvasive MEG (S. Palva & Palva, 2012), respectively. We first found that at the level of the strongest connections, avalanche propagation was well colocalized with both phase-synchronization-based and amplitude-correlation-based FC. To investigate how these interactions were coupled with scale-free dynamics, we estimated the interareal correlations of local LRTCs and found that these correlations were also colocalized with both avalanche propagation and FC. Scale-free neuronal dynamics and FC thus have a shared backbone of interareal relationships.

Both functional (Gallos, Makse, & Sigman, 2012; Meunier, Lambiotte, & Bullmore, 2010) and structural (Hagmann et al., 2008) connectivity have hierarchically modular network organizations. As the second step in this study, we assessed the modularity in our estimates of critical dynamics and found that both the propagation of neuronal avalanches and the interareal correlations of LRTCs exhibited modular structures, and that these structures were neuroanatomically similar to those of FC. This strongly suggests that neuronal communities characterized by strong internal phase synchronization and amplitude correlations are also characterized by preferentially internal avalanche propagation and correlated local LRTCs. Beyond linking scale-free dynamics and connectivity, this finding also indicates that the brain can be envisioned as a constellation of mutually coupled and hierarchically organized modules distributed to distinct neuroanatomical substrates and to different frequency bands, rather than as a homogeneous critical-state system.

Overall, these observations strongly suggest that FC and scale-free neuronal dynamics, and hence, putatively, neuronal criticality as well, are intimately related, which opens new avenues for studies aiming to understand the physiological interactions between these phenomena.

Neuronal criticality is often considered in the context of models or physical systems that have a relatively homogeneous spatial structure of their interacting elements (N. Friedman et al., 2012; D. Plenz & Thiagarajan, 2007). In contrast, human brains have a markedly nonhomogeneous and hierarchically modular large-scale structure (Gallos et al., 2012; Meunier et al., 2010), which theoretically has a major influence on the emergence of criticality in neuronal activity, by transforming the system from having a critical point into exhibiting a critical regime—that is, a Griffith’s phase (Hilgetag & Hutt, 2014; Moretti & Munoz, 2013). One key implication of modular structural connectivity is that it could allow for topologically distant modules of the network to have semi-independent dynamics. In this study, the finding that the interareal relationships of LRTC scaling exponents were clearly nonhomogeneous provides evidence for this theory and opens the possibility that different brain systems could have distinct operating points (or regimes).

We found also that both the interareal relationships and modular structures of neuronal avalanches and LRTCs had colocalized neuroanatomical bases. This finding complements prior studies that have shown that scale-free dynamics observed at the meso- and macro-scales with SEEG and MEG, respectively, are phenomenologically very similar (Zhigalov et al., 2015) and that the avalanche and LRTC scaling exponents are corrected at both spatial scales (J. M. Palva et al., 2013).

The connectivity patterns of oscillations in and around the α (8–14 Hz) and β (15–30 Hz) bands have been suggested to provide a “cortical core” (de Pasquale et al., 2012) or “hubs” (Hipp, Hawellek, Corbetta, Siegel, & Engel, 2012) for neuronal integration and resting-state network dynamics. We found here that the strongest similarities between scale-free dynamics and FC were observed specifically in this frequency range. Likewise, the similarity between neuronal avalanches and LRTC was also maximal in this frequency range, and an even stronger similarity was found for both interareal phase synchronization and amplitude correlations (see Figure 5D). The α and β oscillations in the resting state thus appear to be central for both neuronal integration and the connectivity–dynamics association. Importantly, the lack of such association at the lower frequencies (4–7 Hz) was not explained by a lack of FC or by evidence for power-law scale-free dynamics therein, because highly significant FC as well as LRTCs were observed in this frequency range (see Figure 6). This shows that the association of FC and scale-free dynamics is neither trivial nor automatic in brain dynamics, but rather is a unique dynamic property of a subset of frequency bands in ongoing brain activity.

Two sources of circularity could theoretically confound our analyses of the similarities between broadband avalanche propagation and narrow-band phase synchronization or amplitude correlations. As a conceptual confound, one may argue that the broadband avalanche events and the narrow-band oscillations actually reflect the dynamics of a single underlying neuronal process, which would lead to artificial or trivial colocalization of their connectomes. As a technical confound, the estimates of narrow-band FC could be biased by the putatively large-amplitude avalanche events, and conversely, the avalanche detection might pick up momentarily large-amplitude oscillations because both phenomena were estimated from the same data. Such “technical” cross-bias in the estimators of avalanche propagation and FC could render the estimators mutually and artificially correlated by design. However, control analyses (see Figure 4) showed that (i) the rare avalanche-like events in ongoing brain activity have a negligible impact on estimates of FC, and (ii) both avalanche propagation and FC can be estimated independently of each other. Thus, no technical circularity prevented our testing the hypothesized relationship between avalanches and connectivity. Moreover, because (iii) narrow-band FC during the avalanches was uncorrelated with the avalanche propagation estimated from the same data, the propagating broadband avalanche events appear to be distinct neurophysiological phenomena from the interareal phase or amplitude couplings of narrow-band oscillations. Finally, the finding that the avalanche propagation patterns detected at two SDs were less similar to FC than were those detected at three SDs constitutes strong evidence that the avalanche–connectivity relationship is not a trivial by-product of avalanches simply reflecting phase-lagged connectivity, because if this were the case, the two-SD avalanches should be closer in similarity to the FC estimates (which effectively correspond to the “0-SD” data) than the avalanches at three SDs.

Several lines of evidence have shown that abnormalities in neuronal scaling laws are biomarkers of a number of brain diseases. In major depressive disorder, the LRTC scaling exponents for central/frontal theta oscillations are lower than those among healthy controls, whereas posterior α oscillations exhibit no abnormalities (Linkenkaer-Hansen et al., 2005). In schizophrenia, the LRTC scaling exponents are suppressed in both the α and β bands over parietotemporal areas (Nikulin, Jonsson, & Brismar, 2012). Similarly, in autism spectrum disorders (Lai et al., 2010) and Alzheimer’s disease (Montez et al., 2009), lower-than-normal LRTC scaling exponents are a robust characteristic. Conversely, in epileptic patients, the epileptogenic zone exhibits greater β-band LRTC exponents than do the surrounding cortical regions (Monto, Vanhatalo, Holmes, & Palva, 2007). These findings are well in line with our discovery that different brain systems are capable of operating concurrently in distinct dynamic regimes. The intriguing, and possibly functionally essential, implication is that abnormal dynamics affecting only specific brain systems and/or frequency bands could be predictive of corresponding pathological states and mental symptoms. Neuroanatomical specificity in critical dynamics is also salient in the relationship between healthy neuronal and behavioral fluctuations, in which neuronal criticality only in well-delineated brain systems is predictive of behavioral scaling laws (J. M. Palva et al., 2013). Future studies should thus examine the possibility of modulating cortically well-localized dynamics noninvasively—for example, by using neurofeedback (Ros et al., 2016; Zhigalov, Kaplan, & Palva, 2016).

We analyzed resting-state invasive stereotactic-electroencephalography (SEEG) recordings from a cohort of 27 epileptic patients (18 males, nine females; ages 16–21 years) and noninvasive magnetoencephalography (MEG) data recorded from 14 healthy subjects (seven males, seven females; ages 18–27 years).

Resting-state SEEG was collected for 10 min with eyes closed and without external disturbance using a 192-channel SEEG amplifier system (NIHON-KOHDEN NEUROFAX-110) at a sampling rate of 1 kHz. The number of electrode contacts along each penetrating shaft varied from eight to 15. These contacts were 2 mm long and 0.8 mm thick and had an intercontact distance of 1.5 mm (DIXI medical, Besancon, France). The anatomical positions and amounts of electrodes varied according to surgical requirements (Cardinale et al., 2013). We excluded from further analyses all contacts that were located within the epileptic focus or that exhibited epileptiform activity such as interictal spikes (Zhigalov et al., 2015). The MEG data were collected in a resting-state condition for 10 min while the subjects looked at a fixation point on the monitor screen. We recorded 306-channel (204 planar gradiometers and 102 magnetometers) MEG (Elekta Neuromag Ltd.) at a sampling rate of 600 Hz.

The SEEG study was approved by the Ethics Committee of the Niguarda “Ca’ Granda” Hospital, and the MEG study was approved by the Ethics Committee of Helsinki University Central Hospital. The subjects gave written informed consent for participation in research studies and for publication of the data. The studies were performed according to the Declaration of Helsinki.

The location of each SEEG electrode contact was extracted with submillimeter accuracy from postimplant cone-beam computerized tomography scans (Arnulfo, Narizzano, Cardinale, Fato, & Palva, 2015) and was subsequently coregistered to the Freesurfer (http://surfer.nmr.mgh.harvard.edu/) geometrical space. We used a “closest-white-matter” referencing scheme for the SEEG, in which electrode contacts in gray matter were referenced to the closest contacts in the underlying white matter (Arnulfo, Narizzano, et al., 2015). The MEG data were corrected for extracranial noise and for cardiac and eye blink artifacts by using the signal space separation method (Taulu, Simola, & Kajola, 2005) and independent-component analysis (Bell & Sejnowski, 1995), respectively. For cortical surface reconstructions, T1-weighted anatomical images were recorded with a 1.5-T magnetic resonance imaging (MRI) scanner (Siemens). The FreeSurfer software was used for automatic volumetric segmentation of the MRI data, surface reconstruction, flattening, and cortical parcellation. The MNE software (www.martinos.org/mne/) was used to create three-layer boundary element conductivity models and cortically constrained source models with fixed-orientation dipoles, and for computing the forward and inverse operators (Hamalainen & Ilmoniemi, 1994). The MEG sensor time series were filtered in multiple frequency bands using a broadband finite impulse response filter and Morlet’s wavelets. The filtered time series were inverse-transformed and collapsed into time series of 219 cortical parcels derived from the individual MRI scans (Daducci et al., 2012) by applying collapsing operators that maximized individual reconstruction accuracy (Korhonen, Palva, & Palva, 2014).

The scaling exponent η is defined as the coefficient of the linear regression of the function F(τ) plotted in double logarithmic coordinates. The optimal fitting range, in terms of the linearity of the slope, was limited to 3–300 s for the present dataset.

As in the estimation of amplitude–amplitude correlations, orthogonalization was used to eliminate the pairwise effects of linear signal mixing on the LRTC estimates, and hence to allow for assessment of the relationship between local cortical scaling exponents without artificial coupling of these regions (Blythe, Haufe, Muller, & Nikulin, 2014).

To assess the statistical significance of the results, we generated comparable surrogate data and applied the analyses above to these data. The surrogate time series were derived from the original ones so that the signal of each channel was circularly shifted by a random number of samples (Zhigalov et al., 2015). This approach provides more accurate statistical inference for neuronal avalanches than does the conventional phase-shuffling method (Prichard & Theiler, 1994). The circular shift breaks spatial correlations but essentially fully preserves the temporal autocorrelation structure. For SEEG, the surrogate data were obtained by circular-shifting the white-matter-referenced gray-matter recordings and used as is. For MEG, the surrogate data were obtained by first circular-shifting the original source-reconstructed MEG parcel time series and then forward and re-inverse modeling these data (Korhonen et al., 2014). This approach is essentially equivalent to the one proposed by Shahbazi and colleagues (Shahbazi, Ewald, Ziehe, & Nolte, 2010) and accurately reconstructed the net residual signal mixing present after the MEG measurement (forward modeling) and after the partial separation of the original cortical sources achieved by source modeling (inverse modeling).

To assess whether the avalanche propagation and phase-synchronization or amplitude-correlation connectomes can be estimated independently, we performed several lines of control simulations. We assessed the impacts of suprathreshold events—that is, neuronal avalanches—on FC in MEG using two approaches. In the first approach (see Figure 4A), we deleted all suprathreshold data segments above three SDs and recomputed the phase and amplitude connectomes (see Figure 4A, solid lines). This deletion was carried out separately for each parcel pair in these connectomes by using logical disjunction (i.e., the OR operation) of the avalanche segments in these parcels. We also recomputed the phase and amplitude connectomes for the complement of these data—that is, by using only the suprathreshold segments (see Figure 4A, dashed lines). The suprathreshold, subthreshold, and original functional phase/amplitude connectomes were compared with the avalanche propagation connectome using the edge similarity index (Eq. 14 below). Any statistical difference between the edge similarities at multiple frequencies was assessed with the KS test.

In the second approach (see Figure 4B), we swapped all time windows of the suprathreshold avalanche segments above three SDs between randomly assigned parcel pairs. This operation thus randomly reordered the avalanche propagation connectome. Because the segments were only partially overlapping between the parcel time series, we again used disjunction to combine avalanche segments within the pairs of parcels. After the segment swaps, the time series were refiltered to account for the filter-related temporal spreading of the large-amplitude avalanche events. Using these data, we then recomputed the phase and amplitude connectomes and assessed their edge similarity with the original avalanche connectome (see Figure 4B, solid lines). To assess how independently the avalanches and FC can be measured, we also reanalyzed the avalanche propagation patterns from these data and compared this swapped avalanche propagation connectome with the swapped phase and amplitude connectomes using the edge similarity index (see Figure 4B, dashed lines).

In addition to classical avalanche and amplitude LRTC measures, we characterized criticality in the LRTCs with moment-to-moment fluctuations of interareal phase synchrony (Botcharova et al., 2014). The phase time series between pairs of cortical parcels (MEG) or contacts (SEEG) were extracted from the complex-valued narrow-band-filtered data (4–120 Hz). DFA (see the LRTC Connectome section above) was applied to the first time derivatives of the phase differences. The resulting connectome represented the LRTC scaling exponents of phase synchrony between each pair of cortical parcels.

To assess whether the connection strengths (connectome values) have a specific frequency profile that might explain the profile of connectome similarity, we estimated the connection strengths for the phase, amplitude, and LRTC connectomes for the entire range of frequencies (4–120 Hz). The connection strengths were estimated as a mean for the connectome (adjacency matrix) for each frequency (see Figure 6).

Consensus mapping was applied to assess the shared subgraphs of the connectomes. The subgraph identities (see the previous section) of all four connectomes at a certain frequency (e.g., 10 Hz; see Figures 7E and 7G) were iteratively reordered to maximize overlap between the subgraphs. To do this, the vectors of subgraph identities of size [N × 1] were transformed into binary arrays of size [N × K], where N is the number of parcels (N = 219) and K is the number of subgraph identities (K = 7). Nonzero items of the array represented associations between a parcel (row) and its subgraph identity (column), so that each row contained only one nonzero item. The binary arrays of subgraph identities were assessed for all four connectomes. The columns of these arrays were iteratively reordered to maximize the sum of the element-wise product between the four arrays. The reordered binary arrays were then wrapped back to the subgraph identity vectors and visualized (e.g., Figures 7E and 7G).

Alexander Zhigalov: Data curation; Formal analysis; Methodology; Software; Visualization; Writing – original draft; Writing – review & editing Gabriele Arnulfo: Data curation; Methodology; Validation Lino Nobili: Data curation Satu Palva: Project administration; Resources J. Matias Palva: Conceptualization; Funding acquisition; Project administration; Resources; Supervision; Writing – review & editing

This study was supported by the Academy of Finland, Grant Nos. 253130 and 256472 (to J.M.P.) and 1126967 (to S.P.); by the Helsinki University Research Funds and European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement No. 29 604102 (Human Brain Project); and by the Italian Ministry of Health, Targeted Research Grant No. RF-2010-2319316 (to L.N.).

Functional connectivity:

Defined as the correlation between neurophysiological time series

Scale-free dynamics:

Any fractally self-similar—that is, power-law scaled—dynamics of a system

Neuronal avalanche:

A cascade of neuronal activity with power-law size and lifetime distributions

Long-range temporal correlations:

Functional correlations characterized by power-law decay of the autocorrelation function

Criticality:

The state of a dynamical system at a boundary between two qualitatively different types of dynamics

Neuronal oscillations:

Transiently rhythmic neuronal population activities

Resting-state:

An experimental paradigm in which the participant does not receive any stimulation or task

SEEG:

An invasive, presurgical technique in which local field potentials are recorded intracerebrally to identify the areas where epileptic seizures originate

Critical branching:

A stochastic process characterized by a nearly constant mean population size—that is, without dying out or explosive growth