Avalanche criticality in individuals, fluid intelligence, and working memory.

The critical brain hypothesis suggests that efficient neural computation can be achieved through critical brain dynamics. However, the relationship between human cognitive performance and scale-free brain dynamics remains unclear. In this study, we investigated the whole-brain avalanche activity and its individual variability in the human resting-state functional magnetic resonance imaging (fMRI) data. We showed that though the group-level analysis was inaccurate because of individual variability, the subject wise scale-free avalanche activity was significantly associated with maximal synchronization entropy of their brain activity. Meanwhile, the complexity of functional connectivity, as well as structure-function coupling, is maximized in subjects with maximal synchronization entropy. We also observed order-disorder phase transitions in resting-state brain dynamics and found that there were longer times spent in the subcritical regime. These results imply that large-scale brain dynamics favor the slightly subcritical regime of phase transition. Finally, we showed evidence that the neural dynamics of human participants with higher fluid intelligence and working memory scores are closer to criticality. We identified brain regions whose critical dynamics showed significant positive correlations with fluid intelligence performance and found that these regions were located in the prefrontal cortex and inferior parietal cortex, which were believed to be important nodes of brain networks underlying human intelligence. Our results reveal the possible role that avalanche criticality plays in cognitive performance and provide a simple method to identify the critical point and map cortical states on a spectrum of neural dynamics, ranging from subcriticality to supercriticality.

INTRODUCTION

The critical brain hypothesis states that the brain operates in close vicinity to a critical point that lies between order and disorder. This is characterized by a power law form of the event size distribution (Cocchi et al., 2017; Hesse & Gross, 2014). This hypothesis is supported by a set of observations of power law scaling in many different neural systems using various approaches (J. M. Beggs & Plenz, 2003; Gal & Marom, 2013; Meisel et al., 2013; Plenz, 2012; Shriki et al., 2013; Solovey et al., 2012; Enzo Tagliazucchi et al., 2012). Arguments in favor of this hypothesis have been strengthened by advantages in information transmission, information storage, and dynamic range, in neural systems operating near criticality (Shew et al., 2009; Shew et al., 2011; Yang et al., 2012), with evidence arising in both theoretical and experimental work (Shew & Plenz, 2012). Meanwhile, this hypothesis still faces challenges from several perspectives (J. Beggs & Timme, 2012). For example, computational studies suggested that power laws may emerge from simple stochastic processes or noncritical neuronal systems (Touboul & Destexhe, 2010), so power laws alone are prerequisite but not sufficient evidence for criticality. Meanwhile, it has been asked: “If the brain is critical, what is the phase transition (Fontenele et al., 2019)?” Indeed, the observation of power law avalanche activity along with a phase transition between order and disorder would be more persuasive for criticality. Furthermore, though previous studies have associated supercriticality with reduced consciousness (Meisel et al., 2013; Scott et al., 2014), near‐critical dynamics with rest (Priesemann et al., 2014), and subcriticality with focused cognitive states (Fagerholm et al., 2015), there remains a gap between the specific brain state and efficient information processing endowed by criticality as predicted by theory (He, 2011). To fully understand the functional roles of critical and noncritical dynamics, more research is required to relate brain states and cognitive performance to neural dynamics that lie on a spectrum, ranging from subcriticality to supercriticality. To obtain a deeper understanding of this phenomenon, it is necessary to develop data analysis methods to represent this phase spectrum with high resolution and characterize the subsequent reorganization of brains with the transition in this spectrum (Fontenele et al., 2019).

With advances in brain imaging techniques such as functional magnetic resonance imaging (fMRI), the critical brain hypothesis has found roles in interpreting fundamental properties of large‐scale brain networks in the context of structure–dynamics–function relationships (Karahanoğlu & Van De Ville, 2017; Lee et al., 2019). For example, it has been shown that structural connections of brains are mostly reflected in functional networks, and this structure–function coupling is disrupted when brains move away from criticality during anesthesia (Enzo Tagliazucchi et al., 2016). Another application of criticality is to explain the dynamic basis of brain complexity (Popiel et al., 2020; E. Tagliazucchi & Chialvo, 2013; Timme et al., 2016). In particular, functional connectivity (FC) complexity, which is an umbrella term describing the variability, diversity, or flexibility of functional connections in brain networks, has been associated with cognitive performance from many perspectives, such as high‐order cognition, aging, and cognitive impairment in brain disorders (Ahmadlou et al., 2014; Anokhin et al., 1996; Omidvarnia et al., 2021; Smyser et al., 2016; B. Wang et al., 2017). Studies have suggested that the FC complexity may possibly be at its maximum at the critical point, while the FC capacity arises from special topological properties of the structural network, such as hierarchical modular organization (Song et al., 2019; R. Wang et al., 2019).

To validate these applications, both computer modeling methods and experimental data analysis methods were used. Computer modeling utilizes structural imaging data to model large‐scale brain dynamics and functional networks (Deco et al., 2011; Nakagawa et al., 2013). However, there is still disagreement on which type of phase transition should be adopted for large‐scale brain networks, for example, first‐order discontinuous versus second‐order continuous phase transitions and edge of chaos criticality versus avalanche criticality (Kanders et al., 2017; Scarpetta et al., 2018). Experimental studies usually take advantage of dynamic changes caused by interventions, such as deprived sleep, anesthesia, or brain diseases, to show deviations from criticality and subsequent reorganization of FC networks (Hobbs et al., 2010; Meisel et al., 2013; Meisel et al., 2012; Rolls, 2021; Enzo Tagliazucchi et al., 2016). However, deviations caused by these interventions are usually unidirectional, either in the subcritical or supercritical directions. Furthermore, it is not clear whether deviations from criticality caused by different intervention methods follow an identical phase transition trajectory. Recent studies have proposed the concept of a “critical line” instead of a “critical point” and suggested that multiple phase transition trajectories may exist (Kanders et al., 2020). Therefore, the successful retrieval of the phase transition trajectory from the large‐scale brain networks will not only help to answer key questions regarding what the phase transition is if the brain is critical but also have important implications in brain functional imaging and large‐scale brain modeling.

In this study, we applied a large number of criticality‐related metrics to investigate the spatiotemporal dynamics of large‐scale brain networks as well as their variation in individuals and explored associations with three different cognitive abilities in a sample of 295 healthy young adults from the Human Connectome Project (HCP) 1200‐subject release (Van Essen et al., 2013). Firstly, we performed group‐level analysis with the classic avalanche criticality analysis method (J. M. Beggs & Plenz, 2003; Enzo Tagliazucchi et al., 2012) to show that there was a mismatch between data analysis and theoretical prediction. We then retrieved an inverted‐U curve by mapping individuals' brains onto the phase plane between the mean synchronization (MS) and synchronization entropy (SE) of blood oxygenation level‐dependent (BOLD) signals. We found that for subjects who exhibited moderate mean and maximal variability in synchrony of BOLD signals (i.e., located around the tipping point of the inverted‐U curve), their avalanche distribution was better fitted by a power law, suggesting that they are more likely to have critical dynamics. This is consistent with previous findings that the neural systems operate neither at the synchronous nor at the asynchronous ends of the spectrum but rather near the critical point between them (Fontenele et al., 2019). Meanwhile, the large individual variation around the critical point gave us a chance to examine previous conjectures on criticality in large‐scale brain networks. And we indeed found that both FC complexity and structure–function coupling were maximized around the criticality. We also utilized a sliding window approach to observe “instantaneous” phase transition occurring in individual brains. We found that brains persisting in the subcritical regime exhibited longer dwell times than those in other regimes. Finally, we found that the critical dynamics were associated with high scores in fluid intelligence and working memory tests but not with crystallized intelligence scores. Additionally, the critical dynamics in the frontal cortex, superior parietal lobule, angular gyrus, supramarginal gyrus, and so forth, which were believed as vital regions in the networks of Parieto‐Frontal Integration Theory (P‐FIT) for intelligence (Jung & Haier, 2007), exhibited significant correlations with fluid intelligence performance.

MATERIAL AND METHODS

Data acquisition and preprocessing

fMRI data acquisition and preprocessing

We used resting‐state fMRI (rfMRI) data from the HCP 1200‐subject release (Van Essen et al., 2013). Each subject underwent two sessions of rfMRI on separate days, each session with two separate 14 min 24 s acquisitions generating 1200 volumes on a customized Siemens 3T Skyra scanner (, , , , , , , ). The rfMRI data used for our analysis were processed according to the HCP minimal preprocessing pipeline (Glasser et al., 2016; Glasser et al., 2013) and denoising procedures. The denoising procedure pairs the independent component analysis with the FSL tool FIX to remove nonneural spatiotemporal components (Smith et al., 2015). And as a part of cleanup, HCP used 24 confound time series derived from the motion estimation (the 6 rigid‐body parameter time series, their backwards‐look temporal derivatives, plus all 12 resulting regressors squared). Note that the global component of the fMRI fluctuations measured during the resting state is tightly coupled with the underlying neural activity, and the use of global signal regression as a preprocessing step in resting‐state fMRI analyses remains controversial and is not universally recommended (Liu et al., 2017). Therefore, the global whole‐brain signal was not removed in this work. We used the left‐to‐right acquisitions from the first resting‐state dataset (i.e., resting‐state fMRI 1 FIX‐denoised package).

The first 324 subjects in the dataset entered into our study, and we excluded 29 subjects for missing data. This left us with 295 subjects for further analysis, and 162 of them were females. All the participants were between the ages of 22 and 36, 58 were between the ages of 22 and 25, 130 were between the ages of 26 and 30, 104 were between the ages of 31 and 35, and 3 were 36 years old.

For further analysis, the whole cortex was parcellated into 96 regions using the Harvard‐Oxford atlas (Makris et al., 2006), and the details are provided at https://identifiers.org/neurovault.image:1699, and from each region the time series averaged across the voxels were extracted and Z‐normalized to construct region of interest (ROI) signals (atlas96 signals). To test the results for different parcellation, we also used the Human Brainnetome atlas that comprised 246 regions (Fan et al., 2016) and Zalesky atlas that comprised 1024 regions (Zalesky et al., 2010); the resulting signals were termed as atlas246 and atlas1024 signals, respectively.

Diffusion tensor imaging (DTI) data acquisition and preprocessing

The diffusion MRI images used in this study were also from the HCP 1200‐subject release (Sotiropoulos et al., 2013). Briefly, the diffusion data were collected using a single‐shot, single refocusing spin‐echo, echo‐planar imaging sequence (, , , , , , , , ). The full diffusive session includes 6 runs, representing 3 different gradient tables, with each table acquired once with right‐to‐left and left‐to‐right phase encoding polarities, respectively. Each gradient table contains roughly 90 diffusion weighting directions plus 6 b=0 acquisitions spaced throughout the run. Within each run, diffusion weighting comprised of three shells with b=1000, 2000, and 3000 s/mm2. All diffusion data were preprocessed with the HCP diffusion pipeline updated with EDDY 5.0.10 (Sotiropoulos et al., 2013), and the details are provided at https://www.humanconnectome.org. In this study, from the 295 selected subjects, only 284 subjects were entered into our DTI data analysis because 11 of them were missing the corresponding DTI data.

Cognition measures

We examined associations between cognitive ability and critical dynamics conducted in our rfMRI analysis. Three relevant behavioral tasks were used as a measure of cognitive ability, including fluid intelligence, working memory, and crystallized intelligence. The same 295 subjects were also entered into our cognitive ability analysis, but for the fluid intelligence analysis, five subjects were excluded due to missing their intelligence scores or information about age and education.

The fluid intelligence scores in the HCP data release were measured using the number of correct responses on form A of the Penn Matrix Reasoning Test (PMAT, , , ), which had 24 items and 3 bonus items, using nonverbal visual geometric designs with pieces to assess reasoning abilities that can be administered in under 10 min (Barch et al., 2013; Hearne et al., 2016). The PMAT (Bilker et al., 2012) is an abbreviated version of Raven's Standard Progressive Matrices test (Wendelken et al., 2007), which comprises 60 items.

Crystallized intelligence was measured using the picture vocabulary test (picture vocabulary, mean = 116.8205, SD = 10.1977, range = 92.3914–153.0889) from the National Institutes of Health (NIH) toolbox (Barch et al., 2013; Hearne et al., 2016). This measure of receptive vocabulary was administered in a computer‐adaptive testing (CAT) format. The participant was presented with four pictures and heard an audio recording saying a word and was instructed to select the picture that most closely showed the meaning of the word. Because the test used a variable length CAT with a maximum of 25 items, some participants had fewer items, and the presented words depended on the participant's performance.

Working memory was assessed using the List Sorting Working Memory test (list sorting, mean = 111.2075, SD = 12.0946, range = 84.63–144.50) from the NIH Toolbox (Barch et al., 2013), in which the participants were required to sequence sets of visually and a small number of orally presented stimuli in size order from smallest to biggest. Pictures of different foods and animals were displayed with both a sound clip and a written test that names them and involved two different conditions. In the 1‐list condition, participants ordered a series of objects, either food or animals, but in the 2‐list condition, participants were presented with both animal and food lists and asked to order each list by increasing size. The number of list items increased in subsequent trials, and the task was discontinued after two consecutive incorrect trials.

Data analysis methods

Synchrony and variability in synchrony

We measured the mean and variability in synchronization with a previously described approach (Meisel et al., 2013; Yang et al., 2012). First, we obtained the phase trace from the signal using its Hilbert transform :

Next, we calculated the Kuramoto order parameter as follows: in which is the number of ROIs in global network analysis or the number of voxels in a particular region in regional analysis. The Kuramoto order parameter was used as a time‐dependent measure of phase synchrony of a system. The MS of a time period was calculated as follows: where is the length of the time period. In this study, we calculated static MS of the entire scan period with time points. We derived the entropy of as the measure of variability in synchronization (SE): where is the probability that falls into a bin between . In this study, we chose the number of bins , and the robustness of our results was also tested within an interval between 5 and 100.

Avalanche analysis

In our avalanche analysis, the ROI signals were reduced to a spatiotemporal point process by detecting the suprathreshold peak positions intermediate between two above‐threshold time points, as shown in the example in Figure 1a. By binning the binary sequences with appropriate time resolution (time bin), we obtained a spatial pattern of active ROIs within consecutive time bins. An avalanche was defined as a series of consecutively active bins, which were led and followed by blank bins without activation. The size S and duration T of the avalanches were then defined as the total number of activations and total number of time bins during this avalanche, respectively (J. M. Beggs & Plenz, 2003; Enzo Tagliazucchi et al., 2012).

FIGURE 1

Avalanche statistics obtained from group‐level analysis. (a) Example of a point process (red triangles) extracted from one normalized region of interest (ROI) blood oxygenation level‐dependent (BOLD) signal. (b) The probability distributions of group‐aggregated avalanche sizes for the threshold 1.4 SD and the time bin width of 1 volume (vol.) in the functional magnetic resonance imaging (fMRI) data. The distributions are well approximated by power law with an exponent of with Clauset's test , corresponding to (histogram bins = 40) and branching parameter . (c) The distribution of avalanche durations can be fitted well by a power law with an exponent of under the condition described in b, with Clauset's test . (d) There is a relation between the sizes and duration of the avalanches with a positive index , which is close to . In b–d, the gray open triangles were calculated from the surrogate data. (e) The branching ratio and power law scaling exponents of avalanche sizes for different thresholds used to define the point process

If a system operates near a critical point, the size distribution (), duration distribution (), and average size for a given duration () should be fitted into power laws: where , , and are critical exponents of the system (Friedman et al., 2012; Sethna et al., 2001). Furthermore, the following scaling relation was proposed as an important evaluation of the criticality (Fontenele et al., 2019; Friedman et al., 2012), namely, In this study, we defined to measure the distances of the systems from the critical point, so the smaller the is, the closer the systems are to the critical point.

The scaling exponents governing the power law distribution were estimated using the maximum likelihood estimator (MLE) (Clauset et al., 2009; Marshall et al., 2016). Briefly, the MLE procedure supposes that the empirical data sample is from a power law function in the range , with probability density (Fontenele et al., 2019; Marshall et al., 2016). We estimated critical exponents and by maximizing the likelihood function and via a lattice search algorithm (Marshall et al., 2016). We then used Clauset's goodness‐of‐fit test to quantify the plausibility of fits (Clauset et al., 2009; Deluca & Corral, 2013; Marshall et al., 2016). We used a power law model to produce data sets over the fit range and compared the Kolmogorov–Smirnov (KS) statistics between (1) the real data and the fit against and (2) the model data and the fit. If the real data produced a KS‐statistic that was less than the KS‐statistic found for at least 10% of the power law models (i.e., ), we accepted the data as being fit by the truncated power law because the fluctuations of the real data from the power law were similar in the KS sense to random fluctuations in a perfect power law model.

Surrogate data

To assess the statistical significance of the avalanche analysis results and MS–SE relationship, we generated comparable surrogate data and applied the analyses above to these data. Phase shuffling is often used in hypothesis testing for avalanche size distribution (Gireesh & Plenz, 2008; Shriki et al., 2013). Phase shuffling disrupts temporal as well as spatial correlations in multichannel time series.

Herein phase shuffling was done on the atlas96 signals. The phase randomization procedures were as follows (Prichard & Theiler, 1994): (1) the discrete Fourier transformation was taken to of each subject; (2) rotating the phase at each frequency by an independent random variable that was uniformly chosen in the range . Crucially, the different time series were rotated by the different phases to randomize the phase information; (3) the inverse discrete Fourier transformation was applied to these time series to yield surrogate data.

Branching parameter

The branching parameter , which is defined as the average number of subsequent events that a single preceding event in an avalanche triggers, is a convenient measure to identify criticality (J. M. Beggs & Plenz, 2003). In theory, the system is critical for and subcritical (supercritical) for ). In this study, was calculated as follows: where is the number of ancestors, is the number of descendants in the next time bin, and N is the total number of time bins with activations.

Definition of kappa

A nonparametric measure, κ, for neuronal avalanches was introduced by Shew and his colleagues (Shew et al., 2009). It quantifies the difference between an experimental cumulative density function (CDF) of the avalanche size, , and the theoretical reference CDF, , which is a power law function with theoretical expected exponent : where are avalanche sizes logarithmically spaced between the minimum and maximum observed avalanche sizes and is the number of histogram bins. The unit value of is characteristic of the system in a critical state, whereas values below and above 1 suggest subcritical and supercritical states, respectively.

Functional and structure connectivity (SC) matrix

We constructed an FC matrix from atlas96 signals by computing the Pearson correlation between ROI and ROI , and the mean FC strength was obtained by where means the absolute value.

The SC matrix was constructed using DSI Studio (http://dsi-studio.labsolver.org) from DTI data. The DTI data were reconstructed in the Montreal Neurological Institute (MNI) space using q‐space diffeomorphic reconstruction (F.‐C. Yeh & Tseng, 2011) to obtain the spin distribution function (F. Yeh et al., 2010). A diffusion sampling length ratio of 1.25 was used. The restricted diffusion was quantified using restricted diffusion imaging (F.‐C. Yeh et al., 2017), and a deterministic fiber tracking algorithm (F.‐C. Yeh et al., 2013) was used to obtain one million fibers with whole‐brain seeding. The angular threshold was randomly selected from 15° to 90°. The step size was randomly selected from 0.1 to 3 voxels. The anisotropy threshold was automatically determined by DSI Studio. The fiber trajectories were smoothed by averaging the propagation direction with a percentage of the previous direction. The percentage was randomly selected from 0% to 95%. Tracks with a length shorter than 5 or longer than 300 mm were discarded. The SC matrix was calculated by using the count of the connecting tracks using 96‐region Harvard–Oxford atlas.

FC entropy

The FC entropy is calculated by where is the probability distribution of , i.e., (Yao et al., 2013). In the calculation, the probability distribution was obtained by discretizing the interval (0, 1) into 30 bins.

FC diversity

The functional diversity () of the FC matrix is measured by the similarity of the distribution to the uniform distribution (R. Wang et al., 2019): where is a normalization factor, is in the range [0, 1], and is the probability distribution of , which was obtained by discretizing the interval (0, 1) into M bins ( in this work). For completely asynchronous or synchronized states, the correlation values fall into one bin at 0 or 1, where reflects the simple dynamic interaction pattern. In an extreme case where all types of FC equivalently exist, would ideally follow a uniform distribution (i.e., probability in each bin ) and .

FC flexibility

To obtain the flexibility of FC in the whole brain, we utilized the sliding window method to calculate connectivity number entropy () for each region (Lei et al., 2020; Song et al., 2019). A nonoverlapping sliding window method was applied to the atlas96 signals. The choice of window size must be sufficient to yield a stable Pearson's correlation coefficient within each window yet small enough to reveal the temporal‐dependent variation in FC (Lei et al., 2020; Sakoğlu et al., 2010). We chose a window size in the range of 20–30, corresponding to the number of windows () in the range of 40–60.

Within each time window, we first acquired the FC matrix via their time series in this window. Then, the binary network matrix was obtained by binarizing the FC matrix with a threshold . Subsequently, we calculated the number of regions connected to a particular region () in each time window. Therefore, we could obtain , the probability for a particular connection number occurring, where indicated the ith connection number among all possible connection numbers. Then, for the region , is a complexity measure (i.e., Shannon entropy) for the disorder in the connection numbers over time: where the summation index runs from 1 to the number of all possible connection numbers.

For each subject, the at the whole‐brain level was obtained by simply averaging the regional values over 96 regions:

Similarity between functional and structural networks

To measure the similarity between functional and structural connection networks, the FC matrices were thresholded by to yield binary adjacency matrices such that if and otherwise. The parameter was chosen to fix link density , which was defined as the ratio of the connections in the network () to the total possible number of connections. It is important to fix the link density when comparing networks, as otherwise, differences could arise because the average of the respective is different (and therefore the number of nonzero entries in ) but not because connections are topologically reorganized across conditions (Enzo Tagliazucchi et al., 2016).

The binary FC networks for each subject were compared with the group‐aggregated binary SC network but not with the individual's SC network to avoid fluctuations in individual SC networks. First, the binary adjacency matrices of SC matrices were obtained for each subject such that if there were tracked fiber links; otherwise, . Then, the binary adjacency structural connection matrices were summed up and again thresholded by a thresholding value to yield a group‐aggregated binary SC network. In this way, high values would exclude connections that were shared by fewer subjects but preserve connections that were common in most subjects.

To estimate the similarity between the binary FC network of each subject and the group‐aggregated binary SC network, we computed the Pearson correlation and Hamming distance between these two networks (Enzo Tagliazucchi et al., 2016). Specifically, the Hamming distance is defined as the number of symbol substitutions (in this case 0 or 1) needed to transform one sequence into another and vice versa, and in this case, it is equal to the number of connections that must be rewired to turn the functional network connection into the structural network connection.

Dynamic analysis of phase transitions

We used the sliding window approach to capture the time‐dependent changes in measures used in this study. In the calculation, for the atlas96 signals, the length of the sliding window was set to (volumes), and the sliding step was set to (volumes). In each window, we calculated the corresponding dynamic measures, including dynamic MS , dynamic SE , and dynamic FC matrix . From the dynamic FC matrix , we further obtained dynamic FC entropy , FC diversity , Pearson correlation between FC and SC , and Hamming distance .

Data and code availability

The rfMRI data, DTI data, and cognitive data are available from the HCP at humanconnectome.org, WU‐Minn Consortium. The WU‐Minn HCP Consortium obtained full informed consent from all participants, and research procedures and ethical guidelines were followed in accordance with the Washington University Institutional Review Boards (IRB #201204036; Title: “Mapping the Human Connectome: Structure, Function, and Heritability”). MATLAB (https://www.mathworks.com/) and SPSS (https://www.ibm.com/analytics/spss-statistics-software) were used to conduct the experiments reported in this study. The datasets supporting this article and the codes required to reproduce them can be found online (https://github.com/longzhou-xu/data_and_code_sort.git).

RESULTS

The signature of critical dynamics in the cortical network

For the 295 available subjects, we first investigated the power law distribution of avalanche size at the population level. Here, we defined the activation as the time point when the BOLD signals reached their peak value, while the signals one step before and after this time point were above the chosen threshold (Figure 1a). After preprocessing, the atlas96 signals were converted into point processes in which each time point represented an activation. We then calculated the avalanche size distribution (Figure 1b), as well as the avalanche duration distribution (Figure 1c). First, from the estimated and values, we tested whether the relationship between the scaling exponents holds for different thresholds of ROI signals (Fontenele et al., 2019; Friedman et al., 2012). We found that the closest matching occurred when the chosen threshold was around 1.4 SD (Figure 1d). Second, the power law distribution of avalanche sizes with a slope of could be predicted by theory for a critical branching process with branching parameter (Harris, 1963; Zapperi et al., 1995). However, we found that the threshold of 1.4 SD yielded and (Figure 1e), which did not match well with the theoretical prediction. We ran the same analysis on both atlas246 signals (Figure S1a) and atlas1024 signals (Figure S1b) to find the mismatch still exists (Figure S1c,d).

As the above close check of hallmarks of criticality did not agree with each other well, we moved forward to investigate whether this mismatch could be a result of inter‐subject variability. We calculated both the MS and SE (see synchrony and variability in synchrony) using the atlas96 signals for each of the 295 subjects and characterized the brain states of each subject with points in the MS versus SE phase plane, as seen in the top panel of Figure 2a. We found that the value of MS from these subjects extended from 0.2 to 0.7, and the distribution of subjects was not even but exhibited a greater tendency to the low MS range (Figure 2a, bottom panel). This result suggested that even in the resting state, there is significant variability among the subjects' brain states. It is clearly seen that these state points formed an inverted‐U trajectory in the phase plane. The SE exhibited a maximum at the moderate value of MS, which implied the existence of a state with dynamic richness between order and disorder. We found that the inverted‐U curves and the calculation of SE were robust against different parcellation (Figure S2). We also performed a phase randomization method on the fMRI data and found that this inverted‐U curve disappeared in the randomized surrogate datasets (size = 500, identified by visual inspection; examples can be seen in Figure S3). Therefore, we argued that this inverted‐U curve reflected a special spatiotemporal structure of brain dynamics that did not exist in randomized data.

FIGURE 2

Signatures of criticality as a function of mean synchronization (MS) in resting‐state brain networks. (a) Top panel: The inverted‐U trajectory of the MS versus synchronization entropy (SE) . The red dashed line represents the quadratic fit of the data (, , adjusted ). Bottom panel: The frequency count for the distribution of MS . (b) The branching parameters versus for each subject. The green dashed line indicates . The Pearson correlation value and the value are shown in the figure. The red dashed line represents the linear regression. For further analysis, we selected three representative groups of subjects according to their synchronization level: namely, low mean synchronization (LMS) group (, blue open circles in a and b), moderate mean synchronization (MMS) group (, green open circles in a and b), and high mean synchronization (HMS) group (, red open circles in a and b). (c) Avalanche size distributions for the LMS group, MMS group, and HMS group. To show the difference between these groups, we used gray lines with to guide the eyes. The corresponding group‐aggregated branching parameters are for the LMS group, for the MMS group, and for the HMS group. (d) Avalanche duration distributions for three groups in c. To show the difference between these groups, we used gray lines with to guide the eyes. (e) Scaling relations for the three groups. The blue line and purple line correspond to and , respectively. In the inset, indicates the distance to the critical point. (f) The dependence of on the numbers of histogram bins for the three groups

As expected, with the increasing of MS, the spatiotemporal activation pattern defined before exhibited transitions from random states to ordered states (Figure S4). We then calculated the branching parameter for each subject. We found that with increasing MS, the branching parameter increased from less than 1 to higher than 1, crossing 1 at a moderate value of MS (Figures 2b and S5 for different parcellation).

Furthermore, we selected three groups from the above subjects: the low mean synchronization group (LMS group; the 20 most left subjects in Figure 2a with an MS value of ), the moderate mean synchronization group (MMS group; the 20 subjects located near the peak of curve in Figure 2a with an MS value of ), and the high mean synchronization group (HMS group; the 20 most right subjects in Figure 2a with an MS value of ). For each group, we performed avalanche distribution analysis to identify which group was closest to the critical point (Figure 2c–f). After obtaining scaling exponents and for each group with a threshold of 1.4 SD (Figure 2c,d), the scaling relationship showed the best match for the MMS group (Figure 2e), and more detailed analysis results can be found in Figure S6.

Previous study showed that the truncations of power law fit have a dramatic impact on power law exponents, particularly on the ratio , while barely changes (Destexhe & Touboul, 2021). To test the robustness of our results, we performed analysis in Figure 2c–e with different fitting windows for avalanche size S ( and ) and avalanche duration T (, and ). We found that the ratio of the number of fittings that met the critical criterion ( and Clauset's goodness‐of‐fit test ) to all power law fit samples is highest for MMS group ( for LMS, for MMS, for HMS).

We also calculated , an often‐used parameter that could distinguish the difference between data and the theoretically suggested power law distribution (Fagerholm et al., 2015; Palva et al., 2013; Poil et al., 2012; Shew et al., 2009; Shew et al., 2011). As shown in Figure 2f, as the discrete bin number increases, the values become stable. The stabilized is smaller than 1 for the LMS group but larger than 1 for the HMS and MMS groups. The value for MMS group was closest to 1. In Figure 2f, the minimal and maximal avalanche size used for calculation of is 2 and 150, respectively. The robustness of the results for different choice of cutting off avalanche size is demonstrated in Figure S7. Therefore, the above results suggested that subjects' brains with moderate MS and maximal SE are poised closest to the critical point, supported by consistent hallmarks of criticality. On the other hand, the large dispersion of subjects among the phase space between asynchronous (subcritical) and synchronous (supercritical) states also provides an opportunity to investigate the phase transition in brains.

The complexity in the FC network is maximized by criticality

Since the disorder–order phase transition could be observed, we investigated how this phase transition could impact the organization of FC networks. For convenience, we used MS to indicate this transition. We assessed how the variousness in FC strength changes as the brain undergoes a phase transition from the subcritical to supercritical states. We used FC entropy and FC diversity as measures of variousness in FC strength in the brain networks. FC entropy is a direct measure of Shannon entropy from the probability distribution of FC strength obtained from the FC matrix, whereas FC diversity measures the similarity between the distribution of real FC matrix elements and uniform distribution. In previous studies, the former had been associated with healthy aging (Yao et al., 2013), and the latter is predicted to be maximized at the critical point by a computer model with Ginzburg–Landau equations (R. Wang et al., 2019). We found that both FC entropy (Figure 3a) and FC diversity (Figure 3b) peaked at the moderate value of MS; however, the peak position for these two measures was more rightward than that of SE.

FIGURE 3

Dependence of complexity in the functional connectivity (FC) network on the mean synchronization (MS) of blood oxygenation level‐dependent (BOLD) signals. (a) FC entropy as a function of MS . Red dashed line: quadratic fitting (, , adjusted ). (b) FC diversity as a function of MS . Red dashed line: quadratic fitting (, , adjusted ). (c) FC flexibility as a function of MS . Red dashed line: quadratic fitting (, , adjusted ). The red open triangles represent participants with outliers in the quadratic fittings in a and b

The flexibility in dynamic FC reflects the extent of abundant connection patterns among regions and how frequent switching may occur between different patterns. In this work, we adopted the as a measure of flexibility in FC networks. Our previous study showed that this measure was maximized at the critical point in a large‐scale brain network model that combined DTI structural data and excitable cellular automaton (Song et al., 2019), and this measure could be reduced in the brains of patients with moyamoya disease (Lei et al., 2020). In this study, we found that the flexibility in FC was maximized with a moderate value of MS (Figure 3c). The maximization was robust in a wide range of thresholds and sliding window lengths (Figure S8). This result supported our previous conclusion (Lei et al., 2020; Song et al., 2019). Compared with FC entropy and FC diversity, the peak position for FC flexibility was nearer to the critical point.

FC entropy, diversity, and flexibility are often used in rfMRI studies to measure the complexity in the structure and dynamic reconfiguration of FC networks. Here, the study suggested that the complexity in FC networks is maximized by criticality.

The maximized structure–function coupling around the critical point

From the obtained FC matrix and SC matrix for each subject, we constructed the FC networks for each subject and a group‐aggregated SC network (see similarity between functional and structural networks). We used and to control the link density in the FC networks and the group‐aggregated SC network, respectively. We measured the similarity between the FC network and group‐aggregated SC network with Pearson correlation and Hamming distance. Figure 4a,b demonstrates the dependence of similarity on the MS of each subject with a link density of 0.7 in the FC networks. The similarity between the FC and SC was maximal for subjects with moderate synchrony, as the Pearson correlation was maximized (Figure 4a), while the Hamming distance was minimized (Figure 4b) for these subjects. This maximization of similarity between the FC and SC could be observed in a wide range of link densities in the FC and SC networks. To further consolidate the above results, we measured the similarity between the FC and SC for the three groups (LMS, MMS, and HMS) defined above as a function of the FC network link density. Figure 4c,d shows that as the FC network link density increased, the correlation coefficient between the FC and SC matrices first increased and then decreased, and consistently, the Hamming distance exhibited the opposite tendency. When the FC link density was large, the MMS group showed a significantly higher correlation and a lower Hamming distance between the FC and SC networks than the other two groups. Similarly, by varying , we found that the maximized similarity in the FC and SC at the critical point was robust in the wide range of link densities in the SC network (Figures S9–S11).

FIGURE 4

The dependence of structure–function coupling on the mean synchronization (MS) of brain networks. (a) Pearson correlation between anatomical and functional networks as a function of MS . The link density in the functional connectivity (FC) network and threshold in the group‐aggregated structure connectivity (SC) network = 40 (corresponding to ) are shown in the figure. Red dashed line: quadratic fitting (, , adjusted ), which is better than linear fitting (, , adjusted ). (b) Hamming distance between anatomical and functional networks as a function of MS . Red dashed line: quadratic fitting (, , adjusted ), which is better than linear fitting (, , adjusted ). (c) The Pearson correlation between anatomical and functional networks as a function of FC density ( corresponding to ) for the high mean synchronization (HMS), moderate mean synchronization (MMS), and low mean synchronization (LMS) groups. (d) The Hamming distance between anatomical and functional networks as a function of FC density ( corresponding to ) for the HMS, MMS, and LMS groups. In c and d, the stars indicate significant differences between the HMS and MMS groups (two‐tailed two‐sample t‐test, , uncorrected); the open triangles indicate significant differences between the LMS and MMS groups (two‐tailed two‐sample t‐test, , uncorrected)

We noticed that for a large link density of the FC network, the dependence of similarity on link density monotonically decreased (Figures 4c,d and S9a,b). Since lower link density conserved only stronger links in FC networks, we deduced that structural connections were mostly reflected in the strong functional connections. Meanwhile, the similarity also decreased as the threshold in SC networks decreased (Figure S10a,b), suggesting that the structural connections that were mostly reflected in the functional connections were those shared by most subjects because structural connections specified to individuals would be excluded with high .

The dynamic phase transition in individual subjects' brains

The observed individual brain states dispersed around the critical point provided an opportunity to investigate the dynamic phase transition in individual brains. To this end, for the LMS, MMS, and HMS groups defined above, we randomly selected two subjects from each group. We calculated the dynamical MS and SE for these six subjects with the sliding window approach (Figures 5a and S1–S6) from their Kuramoto order parameters (Figure 5b). We observed a time‐dependent change in individuals' brain states in the state space following the inverted‐U trajectory, as shown in the top panel of Figure 2a. In the time period limited by scan duration, we observed that subjects who were farther away from the critical point tended to stay in the regime decided by MS, and events of crossing the critical point (black lines at ) to the other regime seldom occurred (s1, s2, s5, and s6 in Figure 5a or Figure 5b, top and bottom panel). Subjects who were nearer the critical point were more likely to cross the critical point (s3 and s4 in Figure 5a or Figure 5b, middle panel).

FIGURE 5

The dynamic phase transition in individual subjects' brains. (a) The dependence of dynamic synchronization entropy (SE) on dynamic mean synchronization (MS) from six subjects selected randomly from the low mean synchronization (LMS), moderate mean synchronization (MMS), and high mean synchronization (HMS) groups. The enlarged dark markers indicate the mean position for corresponding subjects (markers with the same shape). (b) The time‐dependent changes in the Kuramoto order parameter for six subjects as demonstrated in a (with the same color). (c) The normalized frequency count of for different levels of , indicated by lines with different colors. (d) The dwell time (the time interval between two successive critical point crossing events) distribution for different levels of . (e, f) The distribution of vertical and horizontal moving distances of phase points in one step of the sliding window. (g, h) The vertical and horizontal velocities of state points of each subject as a function of their MS . The vertical and horizontal velocities were calculated by and , where the symbol indicates the absolute value and was the average across all the windows. is the step used to slide the windows. Here, =10 time points (volumes). Red dashed lines in g and h: quadratic fitting (, , adjusted in (g); , adjusted in h). Both quadratic fittings were better than linear fitting (adjusted in g and adjusted in h)

To validate the above observation at the population level, we divided the 295 subjects at hand into eight groups with different levels of synchrony and calculated the corresponding probability distribution of the Kuramoto order parameter . It is seen clearly from Figure 5c that as the synchrony level decreases, the distribution of the Kuramoto order parameter becomes narrow and less tilted. Meanwhile, we found that the dwell time, which referred to the time interval between two successive critical point crossing events, exhibited heavier tails in its distribution for low synchrony groups (Figure 5d). These results implied the higher inertness in the subcritical regime than others, and brains were more likely to stay in this regime with longer dwell times.

Next, we calculated the distribution of vertical and horizontal moving distances in state space in a fixed time interval (the time points or volumes of one step of the sliding window) for all subjects. We found that the distributions of vertical and horizontal moving distances were both symmetrical with a mean of zero (Figure 5e,f), suggesting that the inverted‐U trajectory in the state space was stable and unlikely to change its shape as time progressed. Furthermore, the position ()‐dependent velocity distribution is maximal for horizontal velocity () and minimal for vertical velocity () near the critical point (Figure 5g,h). The maximal horizontal velocity around the critical point implied that at this point, the systems were most sensitive to the perturbations due to internal fluctuations or external modulations. Meanwhile, the lower vertical and horizontal velocities in the subcritical regime compared with the supercritical regime also reflected the high inertness in the subcritical regime.

It was of interest to determine whether the maximization of FC complexity, as well as function–structure coupling, around the critical point could be realized dynamically when the individual brains endured phase transition. To this end, we obtained the time‐dependent FC matrices with the sliding window method and calculated FC entropy (Figure 6a), FC diversity (Figure 6b), and two measures for similarity between FC and SC (Figure 6c,d) as a function of instantaneous MS in each time window. The time‐dependent complexity and similarity measures followed almost the exact trajectories as those in the static measurements shown in Figure 3a,b, as well as Figure 4a,b. This result implied that FC complexity and similarity between FC and SC were indeed modulated by phase transition in brains, and their maximization could be realized dynamically by positioning the system around the critical point.

FIGURE 6

Dynamic modulations of functional connectivity (FC) complexity and structure–function coupling during the phase transition of brains. (a) The dependence of dynamic FC entropy as a function of instantaneous mean synchronization (MS); thick dashed white line: quadratic fitting (, , adjusted ). (b) The dependence of dynamic FC diversity as a function of instantaneous MS. Thick dashed white line: quadratic fitting (, , adjusted ). (c) The dependence of dynamic FC–SC correlation as a function of instantaneous MS; thick dashed white line: quadratic fitting (, , adjusted ), which is better than linear fitting (adjusted ). (d) The dependence of dynamic FC–SC Hamming distance as a function of instantaneous MS; thick dashed white line: quadratic fitting (, , adjusted ), which is better than linear fitting (adjusted ). In a–d, each dot represents a calculation from one window. The dots with the same color represent the calculation for one subject. However, due to the limited number of colors used, different subjects may share the same color. In c and d, a link density of 0.7 was used to obtain the binary FC network, and a threshold of 40 was used to obtain the group‐aggregated structural network

High fluid intelligence and working memory capacity were associated with critical dynamics

The results above support the hypothesis that large‐scale brain networks lie in the vicinity of a critical point which is associated with moderate MS and maximal SE. Another key prediction from the critical brain hypothesis is that brains that are closer to criticality should be better in cognitive performance. Here, to address this prediction, we assessed linear relationships between SE and intelligence scores of the available subjects. We found that SE values were significantly correlated with fluid intelligence scores (PMAT; Figure 7a) but not with crystallized intelligence scores (picture vocabulary; Figure 7b). Meanwhile, we found that working memory scores, which were assessed using the Listing Sorting Working Memory test from the NIH Toolbox, were significantly correlated with SE (list sorting; Figure 7c). We also noted here that these scores were significantly correlated with many other measures that were found to be maximized at the criticality, namely, FC entropy, FC diversity, and FC flexibility (Figure S12). Meanwhile, there were significant quadratic relationships between MS and fluid intelligence, as well as working memory scores but not for crystallized intelligence scores (Figure 7d–f). We also found that these results still held when potential confounds such as age and education achievements were regressed out (Figures S13 and S14). Therefore, these results imply that brains that are closer to criticality are associated with higher fluid intelligence and working memory scores.

FIGURE 7

Correlations between cognitive performance scores and synchronization entropy (SE), as well as mean synchronization (MS). (a–c) Correlation between SE and Penn Matrix Reasoning Test (PMAT) scores, picture vocabulary test scores, as well as the list sorting working memory test scores. Red dashed lines in a–c: linear fitting. (d) Scatterplot of the PMAT scores against the MS. The red dashed line represents the significant quadratic fit of the data (, , adjusted ), which is better than the linear fitting (adjusted ). (e) Scatterplot of the picture vocabulary test scores against the MS. Both the linear and quadratic regressions are not significant (linear: ; quadratic: ). (f) Scatterplot of the list sorting working memory test scores against the MS. The red dashed line represents the significant quadratic fit of the data (, , adjusted ), which is better than linear fitting (adjusted )

Since a wide variety of experiments have demonstrated that fluid intelligence is associated with a distributed network of regions in the P‐FIT, including frontal areas (Brodmann areas [BAs] 6, 9, 10, 45–47), parietal areas (BA 7, 39, 40), visual cortex (BAs 18, 19), fusiform gyrus (BA 37), Wernicke's area (BA 22), and dorsal anterior cingulate cortex (BA 32) (Jung & Haier, 2007; Nikolaidis et al., 2017), we decided to find more relationships between these regions with critical dynamics indicated by maximized SE.

To obtain the relevant regions in a fine‐grained division of the brain, here we used the Human Brainnetome Atlas, which contains 210 cortical and 36 subcortical subregions (Fan et al., 2016). We extracted from each brain region the voxel‐level BOLD signals and calculated the regional SE for these 246 regions. We found that regions whose SE exhibited significant (, FDR corrected) positive correlations with PMAT scores were located in the frontal areas (i.e., bilateral superior frontal gyrus [SFG], middle frontal gyrus [MFG], precentral gyrus [PrG], right inferior frontal gyrus [IFG], and paracentral lobule [PCL]), parietal areas (i.e., bilateral AG, SMG, Pcun, right SPL), right inferior temporal gyrus (ITG), superior occipital gyrus (sOcG), and left cingulate gyrus (CG) (Figure 8 and Table 1).

FIGURE 8

The brain map for correlations between regional synchronization entropy (SE) and fluid intelligence. The color bar indicated the Pearson correlation value between regional SE and Penn Matrix Reasoning Test (PMAT). The cortical and subcortical regions were defined by the Human Brainnetome Atlas (http://atlas.brainnetome.org/bnatlas.html). Data were visualized using BrainNet Viewer (Xia et al., 2013)

TABLE 1

The brain regions exhibited significant correlation between SE and fluid intelligence

Lobe	Gyrus	Left/right	Brodmann area	MNI coordinate	R value	FDR
Frontal lobe	Superior frontal gyrus (SFG)	SFG_L_7_3	8, 9, 10	[−11, 49, 40]	0.2058	0.0120
		SFG_R_7_3	8, 9, 10	[13, 48, 40]	0.1795	0.0264
		SFG_R_7_7	10	[8, 58, 13]	0.1784	0.0268
	Middle frontal gyrus (MFG)	MFG_L_7_1	9, 46	[−27, 43, 31]	0.1827	0.0264
		MFG_L_7_3	10, 11, 46	[−28, 56, 12]	0.1947	0.0176
		MFG_R_7_3	9, 10, 46	[28, 55, 17]	0.2706	0.0007
		MFG_L_7_5	8, 9, 44, 46	[−33, 23, 45]	0.1871	0.0229
		MFG_R_7_5	9, 44, 45, 46	[42, 27, 39]	0.2522	0.0017
		MFG_L_7_7	10, 11, 47	[−26, 60, −6]	0.1662	0.0429
	Inferior frontal gyrus (IFG)	IFG_R_6_5	45, 47	[42, 22, 3]	0.1820	0.0264
	Precentral gyrus (PrG)	PrG_L_6_1	4, 6	[−49, −8, 39]	0.1797	0.0264
		PrG_R_6_1	4, 6	[55, −2, 33]	0.1718	0.0344
	Paracentral lobule (PCL)	PCL_R_2_2	4, 6	[5, −21, 61]	0.1728	0.0344
Temporal lobe	Inferior temporal gyrus (ITG)	ITG_R_7_5	20, 21, 37	[54, −57, −8]	0.2066	0.0120
Parietal lobe	Superior parietal lobule (SPL)	SPL_R_5_5	7, 40	[31, −54, 53]	0.1717	0.0344
	Angular gyrus (AG)	IPL_R_6_1	39	[45, −71, 20]	0.1687	0.0390
		IPL_R_6_2	7, 39	[39, −65, 44]	0.2373	0.0028
		IPL_L_6_5	39	[−47, −65, 26]	0.2011	0.0140
		IPL_R_6_5	39	[53, −54, 25]	0.2284	0.0043
	Supramarginal gyrus (SG)	IPL_L_6_3	2, 3, 40	[−51, −33, 42]	0.1997	0.0140
		IPL_R_6_3	2, 3, 40	[47, −35, 45]	0.2061	0.0120
		IPL_R_6_4	2, 40	[57, −44, 38]	0.2369	0.0028
	Precuneus (Pcun)	Pcun_L_4_4	7	[−6, −55, 34]	0.1810	0.0264
		Pcun_R_4_4	7	[6, −54, 35]	0.1899	0.0219
Limbic lobe	Cingulate gyrus (CG)	CG_L_7_6	23	[−7, −23, 41]	0.1868	0.0229
Occipital lobe	Superior occipital gyrus (sOcG)	sOcG_R_2_1	18, 19	[16, −85, 34]	0.2051	0.0120

Abbreviations: IPL, inferior parietal lobule; MNI, Montreal Neurological Institute; SE, synchronization entropy.

DISCUSSION

In this study, using the large‐scale WU‐Minn HCP dataset and a large number of criticality‐inspired metrics, we provided evidence that though healthy young brains at rest are on average close to the critical point, as critical brain hypothesis has suggested, there is a considerable individual variation around this point. This gave us a chance to validate previous theoretical predictions of criticality on large‐scale brain networks, and indeed, we observed that the complexity of brain FC, as well as the structure–function coupling, was maximized around the critical point. We proceeded to observe a dynamic phase transition in individual subjects and found that their brains tended to stay subcritical, as indicated by a longer dwell time in this parameter region. Finally, we found that high fluid intelligence and working memory capacity were associated with critical dynamics rather than noncritical dynamics, not only globally but also regionally, suggesting the functional advantages of critical dynamics in resting‐state brains.

Balance between functional segregation and functional integration is a central organizing principle of the cerebral cortex. It has been argued that FC complexity characterizes the interplay of functional segregation and functional integration (Rolls et al., 2021b; Sporns, 2013). A comparison between simulated and empirically obtained resting‐state FC indicates that the human brain at rest lies in a dynamic state that reflects the largest complexity its anatomical connectome can host (Rolls et al., 2021a; G. Tononi et al., 1994). Recently, many studies have tried to link complexity with cognitive performance, human intelligence, and even consciousness, either measured by Φ (big phi) in integrated information theory or discriminated between levels of sedation (Ahmadlou et al., 2014; Duncan et al., 2017; Saxe et al., 2018; Giulio Tononi et al., 1998). Meanwhile, there is a growing awareness that complexity is strongly related to criticality. A recent study showed that criticality maximized complexity in dissociated hippocampal cultures produced from rats (Timme et al., 2016). Here, in this study, we measured FC complexity from different perspectives, either on its strength diversity or on its dynamic flexibility (Figure 3a–c, Figure 6a,b). With the observation of the phase transition trajectory, we demonstrated that these measures of FC complexity were maximized around the critical point. Therefore, the formulation that criticality maximizes complexity was supported in our work empirically with fMRI data at the whole‐brain network level.

It has been shown that human brains possess a stable set of functionally coupled networks that echo many known features of anatomical organization (Krienen et al., 2014). Several computational modeling studies have demonstrated that critical dynamics could best explore the repertoire provided by the structural connectome (Deco & Jirsa, 2012; Enzo Tagliazucchi et al., 2016). Recent studies also suggested that the capacity of repertoire provided by the structural connectome could be extended by the hierarchical modular structural organization (R. Wang et al., 2019). Therefore, structure–function coupling was believed to be at its maximal when the system is at criticality (R. Wang et al., 2019), and it could be disrupted by losing criticality (Cocchi et al., 2014) or disruption of hierarchical organization of structural networks. Previous studies in anesthetized human brains have found structure–function decoupling accompanied by unidirectional departure from a critical point (Enzo Tagliazucchi et al., 2016). It is possible that FC flexibility could be used as a measure of the extent that FC explores the repertoire provided by structural connectome, and the highest FC flexibility occurs when the system is at criticality (Song et al., 2019) (Figure 3c). Our work demonstrated that the maximal exploration of structural connections at the critical point occurs in resting‐state brains (Figure 4). However, since we used a group‐aggregated structural connection networks, we did not investigate how organization of structural connections could impact on the capacity of network repertoire. This issue will be investigated in the future.

Interestingly, although the brain hovers around the critical point, the brain prefers to stay in the subcritical region, as the subject distribution was skewed toward a disordered state, and the dwell time in the subcritical state was longer (Figure 5). Previous analysis of in vivo data has argued that the mammalian brain self‐organizes to a slightly subcritical regime (Priesemann et al., 2014). It was suggested that operating in a slightly subcritical regime may prevent the brain from tipping over to supercriticality, which has been linked to epilepsy. Meanwhile, with a slightly subcritical regime deviates only little from criticality, the computational capabilities may still be close to optimal. However, our results showed that the resting‐state brains could actually stay in the supercritical regimes. So, the preference of brains for subcritical regime may not be because of prevention of too ordered states. In another study, by relating the Electroencephalogram (EEG)‐domain cascades to spatial BOLD patterns in simultaneously recorded fMRI data, the researchers found that while resting‐state cascades were associated with an approximate power law form, the task state was associated with subcritical dynamics (Fagerholm et al., 2015). They argued that while a high dynamic range and a large repertoire of brain states may be advantageous for the resting state with near‐critical dynamics, a lower dynamic range may reduce elements of interference affecting task performance in a focused cognitive task with subcritical dynamics (Fagerholm et al., 2015). Therefore, there remains a possibility that the resting state is not “pure resting state” but mixed with some occasional “task state” for some subjects. However, further delicately designed experimental studies are required to test this conjecture. It remains to uncover the relationship between cognitive states and neural dynamics that lies on a spectrum. The method proposed in this study may be useful in future studies of this topic.

Recently, Ezaki et al. used the Ising model to map BOLD signals on a two‐dimensional phase space and found that human fMRI data were in the paramagnetic phase and were close to the boundary with the spin‐glass phase but not to the boundary with the ferromagnetic phase (Ezaki et al., 2020). Since the spin‐glass phase usually yields chaotic dynamics whereas the ferromagnetic phase is nonchaotic, their results suggested that the brain is around the “edge of chaos criticality” instead of “avalanche criticality.” However, our findings support that avalanche criticality occur in large‐scale brain networks. Therefore, it is interesting to investigate whether both kinds of criticality could co‐occur in large‐scale brain networks (Kanders et al., 2017). Ezaki et al. also found that criticality of brain dynamics was associated with human fluid intelligence, though they used performance IQ to reflect fluid intelligence, which refers to active or effortful problem solving and maintenance of information. In our work, we assessed the correlation between fluid intelligence and the critical dynamics indicated by SE for brain regions and found that regions that showed significant positive correlations were located in parietal–frontal network (Figure 8 and Table 1). These regions were most frequently reported in studies of intelligence and its biological basis, including structural neuroimaging studies using voxel‐based morphometry, magnetic resonance spectroscopy, and DTI, as well as functional imaging studies using positron emission tomography (PET) or fMRI (Jung & Haier, 2007). Also, in the P‐FIT of intelligence, these regions are considered as the most crucial nodes of the brain network underlying human intelligence (Jung & Haier, 2007; Nikolaidis et al., 2017).

Our study suggested that not only fluid intelligence but also working memory capacity was associated with critical dynamics. This is possibly because working memory may share the same capacity constraint through similar neural networks with fluid intelligence (Halford et al., 2007; Jaeggi et al., 2008; Kane & Engle, 2002). In our study, the critical dynamics in the frontal and parietal network also exhibited significant correlation with working memory capacity (Figure S15 and Table S1). Furthermore, it has been well established that working memory is strongly modulated by dopamine, and too strong or too weak dopamine D1 activation is detrimental for working memory, with the optimal performance achieved at an intermediate level (Cools & D'Esposito, 2011; Vijayraghavan et al., 2007; Zahrt et al., 1997). This inverted‐U dose–response has been observed in mice (Lidow et al., 2003), rats (Zahrt et al., 1997), monkeys (Cai & Arnsten, 1997), and humans (Gibbs & D'Esposito, 2005). Recent studies on neural network models have shown that the optimal performance of working memory co‐occurs with critical dynamics at the network level and the excitation‐inhibition balance at the level of individual neurons and is modulated by dopamine at the synaptic level through a series of U or inverted‐U profiles (Hu et al., 2019). Here in this study, we demonstrated that the optimal performance of working memory and criticality co‐occurs at the system level.

However, our study had several limitations. Firstly, the surrogate data test used in this study ruled out the possibility that the results we obtained could be explained by autocorrelations in the data. However, the long‐range spatial correlation of criticality cannot allow one to test the results by ruling out the effects of correlation across the time series. Secondly, though we used the denoising fMRI data from HCP with standard data preprocessing procedure, it is still interesting to investigate how the preprocessing procedure affects the results. Thirdly, in the avalanche analysis, the activation events defined in this study were slightly different from definition used by others, such as threshold‐crossing events (Enzo Tagliazucchi et al., 2012) or above‐threshold events (Bocaccio et al., 2019; R. Wang et al., 2019). We compared these different methods and found that all these methods could generate scale‐free avalanche activities, but unlike our method, the other two methods failed to generate critical branching process (See Section II in Supporting information). Therefore, it is interesting to investigate the correlations between neural activities and events detected by different detection methods from BOLD signals. Finally, in this study, we only focused on the cognitive abilities that are associated with critical dynamics and found significant but not strong correlations between fluid intelligence, working memory, and critical dynamics. Recent works demonstrated that the functional network segregation, integration, and their balance could predict different cognitive abilities (R. Wang et al., 2021). Therefore, future investigation on the relationship between this functional balance and criticality across individuals may reveal various associations between diverse cognitive abilities with not only critical dynamics but also noncritical dynamics. Despite all these shortages, we hope that this study may inspire future research work to validate our findings, for example, through observing not only the association between the departure of criticality and the decline of cognitive performance, either in aging or brain disease, but also the restoration of criticality and the improvement of cognitive performance with pharmacological or noninvasive brain stimulation.

CONCLUSIONS

In conclusion, we mapped individuals' brain dynamics from resting‐state fMRI scans on the phase transition trajectory and identify subjects who are close to the critical point. With this approach, we validated two predictions of critical brain hypothesis on large‐scale brain networks, that is, maximized FC complexity and maximized structure–function coupling around the critical point. We also observed the tendency of brain to stay in subcritical regime. Finally, we found that the critical dynamics in large‐scale brain networks were associated with high scores in fluid intelligence and working memory, implying the vital role of large‐scale critical dynamics in cognitive performance. We also identified key brain regions whose critical dynamics was highly correlated with human intelligence. Our findings support the critical brain hypothesis that neural computation is optimized by critical brain dynamics, as characterized by scale‐free avalanche activity, and could provide a solution for improving the effects of future interventions targeting aspects of cognitive decline (Reinhart & Nguyen, 2019), possibly by controlling the criticality through noninvasive stimulation (Chialvo et al., 2020).

CONFLICT OF INTEREST

The authors declare no competing interests.

Appendix S1: Supporting Information

Click here for additional data file.