Introduction

Critical brain hypothesis posits that the brain operates near a critical regime, i.e., boundary between different phases showing qualitatively different behaviors1,2,3,4,5,6. This hypothesis has been investigated for more than two decades including criticisms such as the presence of alternative mechanisms explaining power law scaling in the relevant observables7,8,9,10. Experimental evidence such as the recovery of critical behavior after interventions, which is difficult to explain by alternative mechanisms, lends supports to the hypothesis9.

Theoretical and experimental work has shown that neural systems operating near criticality are advantageous in information transmission, information storage, classification, and nonlinear input filtering1,3,5,11,12,13,14. These findings align with the idea of edge-of-chaos computation, with which computational ability of a system is maximized at a phase transition between a chaotic phase and a nonchaotic phase15,16,17. These findings are also in line with a general contention that cognitive computations occur as neural dynamical processes18,19.

A prediction from the critical brain hypothesis is that neural dynamics in individuals with higher cognitive abilities should be closer to criticality than in those with lower cognitive abilities. However, whether high cognitive skills are associated with criticality has not been empirically proven. In fact, recent emerging evidence suggests that human cognitive performance is associated with appropriate transitions between relatively discrete brain states during rest20,21,22, working memory tasks23, and visual perception tasks24. Furthermore, these and other studies18,19,25 support that state-transition dynamics in the brain involve large-scale brain networks. These arguments are consistent with the proposal that many cognitive functions seem to depend on network connectivity among various regions scattered over the whole brain26. On these grounds, in the present study we hypothesize that complex and transitory neural dynamics of the brain network (i.e., dynamic transitions among discrete brain states) that are close to criticality are associated with high cognitive performance of humans.

Two major conventional methods for examining criticality and edge-of-chaos computation in empirical neural data are not capable of testing this hypothesis for their own reasons. First, many of the experimental studies testing the critical brain hypothesis have examined neuronal avalanches11,12, including the case of humans5,27,28. Neuronal avalanches are bursts of cascading activity of neurons, whose power-law properties have been related to criticality. However, studies of neuronal avalanches have focused on their scale-free dynamics in space and time, with which statistics of avalanches obey power laws. Scale-free dynamics of neuronal avalanches is a question orthogonal to patterns of transitions between discrete states. Second, nonlinear time series analysis has found that electroencephalography signals recorded from the brains of healthy controls are chaotic and that the degree of chaoticity is stronger for healthy controls than individuals with, for example, epilepsy, Alzheimer’s disease, and schizophrenia29. However, this method is not usually for interacting time series. Therefore, it does not directly reveal how different brain regions interact or whether possible critical or chaotic dynamics are an outcome of the dynamics at a single region or interaction among different regions.

In the present study, we develop a data-driven method to measure the extent to which neural dynamics obtained from large-scale brain networks are close to criticality and complex state-transition dynamics. The method stands on two established findings. First, statistical mechanical theory of the Ising spin-system model posits that the so-called spin-glass phase corresponds to rugged energy landscapes (and therefore, complex state-transition dynamics)30 and chaotic dynamics31,32,33. Therefore, we are interested in how close the given data are to dynamics in the spin-glass phase. Second, the Ising model has been fitted to various electrophysiological data6,34,35,36 and fMRI data recorded from a collection of regions of interest (ROIs)20,21,24,37,38 during rest or tasks with a high accuracy. Therefore, we start by fitting the Ising model to the multivariate fMRI data. Then, we draw phase diagrams of functional brain networks at a whole-brain level. By construction, the dynamical behavior of the system is qualitatively distinct in different phases. The method determines the location of a brain in the phase diagram and thus tells us whether the large-scale brain dynamics of individual participants are ordered, disordered, or chaotic (i.e., spin-glass) dynamics as well as how close the dynamics are to a phase transition curve, on which the system shows critical behavior.

We deploy this method to resting-state fMRI data recorded from human adults with different intelligence quotient (IQ) scores. As a cognitive ability of interest, we focus on fluid intelligence, which refers to the ability to think logically and solve problems with a limited amount of task-related information39. Fluid intelligence is strongly related to the general intelligence factor, g39 and predictive of real-world outcomes such as job performance40. We examine our hypothesis that large-scale brain dynamics of individuals higher in the intelligence score that measures fluid intelligence are closer to critical.

Results

Brain dynamics are close to the spin-glass phase transition

We first fitted the pairwise maximum entropy model (PMEM), which assumes pairwise interaction between ROIs and otherwise produces a maximally random distribution, which is a Boltzmann distribution. The PMEM is equivalent to the inverse Ising model, where the parameters of the Ising model are inferred from data. Because the model assumes binary data, we binarized the resting-state fMRI signals obtained from 138 healthy adults. The binarized activity pattern at N( = 264) ROIs41 at time t (\(t=1,\ldots ,{t}_{\max }\); \({t}_{\max }=258\)) is denoted by S(t) = (S1(t), …, SN(t)) ∈ {−1, +1}N, where Si(t) = 1 and Si(t) = − 1 (i = 1, …, N) indicate that ROI i is active (i.e., the fMRI signal is larger than a threshold) and inactive (i.e., smaller than the threshold), respectively. We fitted the following probability distribution to the population of the 138 participants by maximizing a pseudo likelihood (see “Methods”)24,34:

$$P({\bf{S}}| {\bf{h}},{\bf{J}})=\frac{\exp \left[-E({\bf{S}}| {\bf{h}},{\bf{J}})\right]}{\sum _{{\bf{S}}\in {[-1,1]}^{N}}\exp \left[-E({\bf{S}}| {\bf{h}},{\bf{J}})\right]}.$$
(1)

In Eq. (1),

$$E({\bf{S}}| {\bf{h}},{\bf{J}})=-\sum _{i=1}^{N}{h}_{i}{S}_{i}-\frac{1}{2}\sum _{i=1}^{N}\sum _{j=1,\, j\ne i}^{N}{J}_{ij}{S}_{i}{S}_{j}$$
(2)

is the energy of activity pattern S, \({\bf{h}}=\left\{{h}_{i}:1\le i\le N\right\}\), and \({\bf{J}}=\left\{{J}_{ij}:1\le i\,\ne\, j\le N\right\}\), where Ji = Jji. Although we refer to E as the energy, E does not represent the physical energy of a neural system but is a mathematical construct representing the frequency with which activity pattern S appears in the given data. Activity pattern S appears rarely in the data if E corresponding to S is large and vice versa. Parameter hi represents the tendency that Si = 1 is taken because a positive large value of hi implies that Si = 1 as opposed to Si = −1 lowers the energy and hence raises the probability that S with Si = 1 appears. Parameter Jij represents a functional connectivity between ROIs i and j because, if Jij is away from 0, Si, and Sj would be correlated in general. We denote the estimated parameter values by \(\hat{{\bf{h}}}\) and \(\hat{{\bf{J}}}\).

Then, to evaluate how close the current data are to criticality, we drew phase diagrams by sweeping values of J. In the phase diagrams, we fixed h at \(\hat{{\bf{h}}}\) following the theoretical convention30, including when the PMEM is applied to data analysis6. We set \({\bf{h}}=\hat{{\bf{h}}}\) also because changing the h values did not qualitatively change the phase diagrams (Supplementary Fig. 1). Then, we varied the mean μ and standard deviation σ of J by linearly transforming J, i.e.,

$${J}_{ij}=({\hat{J}}_{ij}-\hat{\mu })\frac{\sigma }{\hat{\sigma }}+\mu .$$
(3)

In (3), \(\hat{\mu }=1.57\times 1{0}^{-3}\) and \(\hat{\sigma }=3.57\times 1{0}^{-2}\) are the mean and standard deviation of the off-diagonal elements of \(\hat{{\bf{J}}}\) estimated for the empirical data. We chose the parametrization given in Eq. (3) motivated by the past investigation of the archetypical Sherrington-Kirkpatrick (SK) model of spin systems30. The SK model, a type of Ising model, is defined with parameters Jij (1 ≤ i ≠ j ≤ N) that are independently drawn from the Gaussian distribution with the tunable mean and standard deviation and has extensively been studied for investigating the so-called spin-glass phase transition owing to its theoretical tractability. In the spin-glass phase, the system shows a disorderly frozen pattern of spins rather than uniform or periodic ones. For each set of Jij values (1 ≤ i ≠ j ≤ N) specified by a (μ, σ) pair, we performed Monte Carlo simulations and calculated observables (see “Methods”). In this manner, we generated a phase diagram for each observable in terms of μ and σ.

Two primary observables (called order parameters in physics literature) employed in studies of spin systems are the magnetization, denoted by m, and the spin-glass order parameter, denoted by q. The magnetization is defined by m = \(\sum\)1≤i≤N〈Si〉∕N, where 〈 ⋅ 〉 represents the ensemble average, and quantifies the mean tendency that Si = 1 as opposed to Si = −1 is taken across the ROIs. The spin-glass order parameter is defined by \(q={\sum }_{1\le i\le N}{\langle {S}_{i}\rangle }^{2}/N\) and represents the degree of local magnetization at individual ROIs. We show m and q as functions of μ and σ in Fig. 1a, b, respectively. The obtained phase diagrams were qualitatively the same as those for the SK model of the same system size, which was given by Eqs. (1) and (2) with each Jij (=Jji, i ≠ j) being independently drawn from a Gaussian distribution with mean μ and standard deviation σ (Fig. 1e, f). The parameter space is composed of three qualitatively different phases30. The paramagnetic phase, characterized by m = 0 and q = 0 in the limit of N → ∞, represents the situation in which each Si randomly flips between 1 and −1, yielding no magnetization. The ferromagnetic phase, characterized by m ≠ 0 and q > 0, represents the situation in which (almost) all Si’s align in one direction (i.e., Si = 1 or Si = −1). The spin-glass (SG) phase, characterized by m = 0 and q > 0, represents the situation in which each Si is locally magnetized but not globally aligned to a specific direction30. Note that the finite size effect of our system blurred the boundaries between the different phases. The current data pooled across the participants lie in the paramagnetic phase and are close to the boundary to the SG phase (crosses in Fig. 1a, b). In theory, the spin-glass susceptibility, \({\chi }_{{\rm{SG}}}={N}^{-1}{\beta }^{2}{\sum }_{1\le i,j\le N}{c}_{ij}^{2}\), where cij = 〈SiSj〉 − mimj, diverges on the boundary between the paramagnetic and SG phases30. The empirical data yielded a relatively large χSG value in the phase diagram (Fig. 1c). In contrast, we did not find a signature of phase transition in terms of the uniform susceptibility defined by χuni = N−1β\(\sum\)1≤i,j≤N cij, which characterizes the transition between the paramagnetic and ferromagnetic phases30 (Fig. 1d). Note that the phase diagrams for χSG and χuni resemble those obtained from the SK model (Fig. 1g, h).

Fig. 1
figure 1

a–d Phase diagrams for the empirical data. e–h Phase diagrams for the SK model. a, e: ∣m∣. b, f: q. c, g: χSG. d, h: χuni. In a, d the crosses represent the mean and standard deviation of the Jij estimated for the entire population of the participants, i.e., (\(\hat{\mu }\), \(\hat{\sigma }\)). In c a circle represents a participant. In a and e we plot ∣m∣ instead of m. This is because averaging over simulations and over realizations of J would lead to m ≈ 0 due to symmetry breaking, even if m ≠ 0 in theory such as in the ferromagnetic phase. i χSG as a function of σ, with \(\mu =\hat{\mu }\) being fixed. j χuni as a function of μ, with \(\sigma =\hat{\sigma }\) being fixed. In i and j, the curves are the cross-sectional view of c and d, respectively, along the dashed line in c or d. The circles in i and j represent the individual participants and are the projection of the circles in c and d onto the dashed line. k Scaling behavior of χSG when the system size \({N}^{\prime}\) is varied. The value of σ = σpeak that maximizes χSG is plotted against \(1/{N}^{\prime}\) in the inset. The dashed line is the linear regression based on the six data points, \({N}^{\prime}=40,60,80,120,160\), and 264. The coefficient of determination is denoted by R2.

Full size image

Next, we examined where brain activity patterns of each participant were located in the phase diagrams. We did so by finding the μ and σ values corresponding to the χSG and χuni values of each participant (see “Methods”). It should be noted that χSG and χuni can be calculated for each individual only from the covariance matrix of the data, without estimating the PMEM. The location of each participant in the phase diagram of χSG is shown by the circles in Fig. 1c. The cross section of this phase diagram for \(\mu =\hat{\mu }\) (along the dashed line shown in Fig. 1c) is shown in Fig. 1i. We also projected the χSG values for the individual participants (circles in Fig. 1i) based on the value of σ estimated for each individual (circles in Fig. 1c). Figure 1i suggests that the empirical data are located in a range of σ that constitutes a peak, further confirming that the brain dynamics of different participants are close to the paramagnetic–SG phase transition and to different extents. In contrast, the participants were far from the paramagnetic–ferromagnetic phase boundary. This is confirmed in Fig. 1j, which is a cross section of the phase diagram for χuni (along the dashed line shown in Fig. 1d) together with the χuni values for the single participants.

The χSG value for the individual participants was off the largest possible values in the phase diagram (Fig. 1i). To examine this point, we carried out a finite size scaling on χSG (Fig. 1k). To emulate systems of smaller sizes than N = 264, we selected \({N}^{\prime}\) out of the N ROIs uniformly at random and fitted the PMEM. The estimated parameter values are denoted by \(\hat{{\bf{h}}}\) and \(\hat{{\bf{J}}}\) without confusion. Then, we simulated the equilibrium state of the system by scanning J according to Eq. (3), where we varied σ while fixing \(\mu =\hat{\mu }\). In this manner, we sought to investigate how close the data were to the SG phase transition at each \(N^{\prime}\) value. As shown in Fig. 1k, the peak value of χSG increased as \(N^{\prime}\) increased, suggesting that the paramagnetic–SG phase transition is approached as the system size increases. In addition, the position of the peak, denoted by σpeak, shifted toward the value for the empirical data, \(\hat{\sigma }\), as N increased. By regressing \({\sigma }_{{\rm{peak}}}/\hat{\sigma }\) linearly on \(1/N^{\prime}\) (inset of Fig. 1k), we estimated \({\sigma }_{{\rm{peak}}}/\hat{\sigma }=1.45\pm 0.04\) in the limit \(N^{\prime} \to \infty\).

The performance IQ is associated with the criticality

To test our hypothesis that criticality of brain dynamics is associated with human fluid intelligence, we examined the correlation between χSG, which encodes the proximity of each participant’s neural dynamics to the paramagnetic–SG phase transition (Fig. 1c, i), and the performance IQ score. The performance IQ score is defined based on tasks that are reflective of fluid intelligence42,43. An enlargement of Fig. 1c is shown in Fig. 2a, where the participants are shown in different colors depending on whether they have a higher performance IQ score (defined by the score value larger than or equal to the median, 109, n = 68) and a lower score (n = 63). We found that higher-IQ participants tended to be closer to the paramagnetic–SG phase transition than lower-IQ participants, as measured by σ (t129 = 3.17, P < 0.002, Cohen’s d = 0.55 in a two-sample t test). The results were qualitatively the same when the outliers were excluded (t127 = 3.52, P < 10−3, d = 0.62). In contrast, the two groups were not different in terms of the distance to the paramagnetic–ferromagnetic phase transition as measured by μ (t129 = 0.77, P = 0.44, d = 0.13 with the outlier included; t127 = 0.85, P = 0.40, d = 0.15 with the outlier excluded).

Fig. 2: Association between the spin-glass susceptibility and the IQ scores.
figure 2

a Magnification of Fig. 1c. The blue and red circles represent participants with a high performance IQ score (≥109) and a low performance IQ score (<109), respectively. The two overlapping histograms on the horizontal axis are the distributions of \(\tilde{\mu }\) for each participant group. The histograms on the vertical axis are the distributions of \(\tilde{\sigma }\). b Relationship between χSG and the performance IQ. A solid circle represents a participant. The participants enclosed by the dashed circle represent outliers determined by Tukey's 1.5 quartile criteria45. The Pearson correlation value (i.e., r) and the P value shown in the figure are those calculated in the presence of the outliers. The solid line is the linear regression. c Relationship between χSG and the verbal IQ. d Relationship between χSG and the full IQ. The χSG and IQ values shown in b, c, and d are those after the effects of the age and the gender have been partialed out.

Full size image

More systematically, we found a mild positive correlation between χSG and the performance IQ score (r129 = 0.24, PBonferroni = 0.011; also see Fig. 2b). However, the verbal IQ score, which is based on individuals’ verbal knowledge42,43, was not correlated with χSG (r126 = 0.06, Puncorrected = 0.50, Fig. 2c). The correlation between χSG and the performance IQ score was also significantly larger than the correlation between χSG and the verbal IQ score (t121 = 2.33, P = 0.021, in the Williams t test for comparing two nonindependent correlations with a variable in common44). These results suggest that the criticality of brain dynamics plays more roles in fluid intelligence than when simply retrieving verbal knowledge. Note that we partialed out the effects of the age and gender in this and the following analysis unless we state otherwise.

The correlation between the full IQ score42,43 and χSG was intermediate between the results for the performance and verbal IQ scores (r130 = 0.19, P = 0.026; also see Fig. 2d), which is natural because the performance and verbal IQ scores are components of the full IQ score.

The association between the spin-glass susceptibility, χSG, and the different types of IQ scores were robust in the following four ways. First, the exclusion of the two outliers determined by Tukey’s 1.5 criteria45 did not affect the significance of the results (χSG vs performance IQ: r127 = 0.27, PBonferroni = 0.005; χSG vs verbal IQ: r124 = 0.13, PBonferroni = 0.27; χSG vs full IQ: r128 = 0.25, P = 0.005). Second, the results were robust against variation on the threshold value for binarizing the fMRI signal (Supplementary Fig. 2). Furthermore, changes in the threshold value did not substantially alter the phase diagrams (Supplementary Fig. 3). Third, the results were preserved even when the global signal (see Methods) was not subtracted from the fMRI signals (χSG vs performance IQ: r129 = 0.22, PBonferroni = 0.02; χSG vs verbal IQ: r126 = 0.046, Puncorrected = 0.61; χSG vs full IQ: r130 = 0.18, P = 0.043; the outliers were not removed). Fourth, we did not find a gender difference in the correlation coefficient between χSG and the IQ scores (performance IQ: Z = 0.33, P = 0.74 in a Z-test for a pair of correlation coefficients46; verbal IQ: Z = 0.43, P = 0.67; full IQ: Z = 0.17, P = 0.86). In this gender-difference analysis, we partialed out the effect of the age but not the gender.

Irrelevance of the paramagnetic–ferromagnetic transition

The IQ was not correlated with χuni (performance IQ: r129 = 0.10, Puncorrected = 0.27; verbal IQ: r126 = 0.093, Puncorrected = 0.30; full IQ: r130 = 0.10, P = 0.24, each test including the outliers; performance IQ: r124 = 0.013, Puncorrected = 0.89; verbal IQ: r121 = 0.039, Puncorrected = 0.67; full IQ: r125 = 0.020, P = 0.82, each test excluding the outliers). The specific heat (denoted by C; see “Methods” for definition) was only mildly correlated with the performance IQ score (performance IQ: r129 = 0.21, PBonferroni = 0.034; verbal IQ: r126 = −0.0056, Puncorrected = 0.95; full IQ: r130 = 0.13, P = 0.14, each test including the outliers; performance IQ: r125 = 0.16, Puncorrected = 0.08; verbal IQ: r122 = −0.016, Puncorrected = 0.86; full IQ: r126 = 0.10, P = 0.26, each test excluding the outliers). Because χuni and C diverge in the paramagnetic–ferromagnetic phase transition but not in the paramagnetic–SG phase transition30, these negative results lend another support to the relevance of the SG phase rather than the ferromagnetic phase to intelligence.

Consistency with the critical slowing down analysis

The previous literature used various measures of criticality. We measured for each participant such a measure, i.e., the scaling exponent of autocorrelation47,48. This measure quantifies the critical slowing down phenomenon, which has been observed in critical states of the brain48. Note that this index quantifies temporal correlation and is orthogonal to what we have measured. We computed the scaling exponent for the autocorrelation function of the fMRI signal at each ROI, using the detrended fluctuation analysis47,48. Then, we took the average of the scaling exponent over the N = 264 ROIs for each participant, which is denoted by α. The association between α and the IQ scores was consistent with the results for χSG (α vs performance IQ: r129 = 0.29, PBonferroni = 0.002; α vs verbal IQ: r126 = 0.19, PBonferroni = 0.068; α vs full IQ: r130 = 0.25, P = 0.003). These results were robust against the removal of outliers (α vs performance IQ: r128 = 0.28, PBonferroni = 0.002; α vs verbal IQ: r125 = 0.17, PBonferroni = 0.10; α vs full IQ: r130 = 0.25, P = 0.003).

We then performed a multivariate linear regression of the performance IQ with χSG and α being the independent variable. We found a significant regression equation (F2,128 = 8.0, P < 0.001, adjusted R2 = 0.11). Both χSG and α were significantly correlated with the performance IQ (χSG: β = 0.18, P = 0.039; α: β = 0.24, P = 0.0067). This result implies that the association between χSG and the performance IQ that we have found is not a byproduct of that between α and the performance IQ. The variance inflation factor for both independent variables was equal to 1.07; this value is small enough for justifying the use of the multivariate regression.

Effects of data length and individual variability

We examined if the limited data length and between-participant variability in our data influenced our results. First, we investigated how the estimation of the individual participant’s χSG and χuni depended on the length of her/his fMRI data (Supplementary Fig. 4a). The results were qualitatively the same as those obtained with all the data if we used approximately more than two thirds of the data (i.e., number of volumes per participant larger than ≈ 150). The correlation between χSG and the IQ scores and that between χuni and the IQ scores were also preserved with the aforementioned data length (Supplementary Fig. 4b, c). Therefore, our main results based on the χSG and χuni are considered to be reliable in terms of the data length.

Second, as we did in our previous studies21,37, we divided the participants into two subgroups of the same size and ran some of the main analyses for the subgroups. We started by comparing the pairwise activity correlation, 〈SiSj〉, for each (i, j) pair between the two subgroups. The 〈SiSj〉 values were strongly correlated between the subgroups and also between the empirical data and estimated PMEMs for the two subgroups (Supplementary Fig. 5). We further confirmed that the phase diagrams were similar between the two subgroups (Supplementary Fig. 6). Moreover, we estimated \(\tilde{\mu }\) and \(\tilde{\sigma }\) for each participant only using the subgroup of participants to which the focal participant belongs. The results were similar to those estimated based on all the participants (Supplementary Fig. 7). Therefore, we conclude that the estimation of the phase diagrams (Fig. 1a–h) and their derivatives (i.e., \(\tilde{\mu }\) and \(\tilde{\sigma }\)), which are based on the estimated phase diagrams, are robust enough against fluctuations in data, such as those caused by a reduced number of participants.

Discussion

We provided empirical support that neural dynamics of humans with higher intellectual ability are closer to critical. The present results are consistent with the standing claim of the “critical brain hypothesis” and “edge-of-chaos computation”, which jointly dictate that the brain is maximizing its computational performance by poising its dynamics close to the criticality, particularly the criticality involving a chaotic regime.

Here we presented an explicit, albeit only moderate, correlation between the IQ scores and the distance from criticality at an individual’s level. Human intelligence has been shown to be associated with genetic factors, brain size, the volume of specific brain regions49, and the structure of brain networks26,49. The present results derived from dynamic fMRI signals provide an orthogonal account of human intelligence as compared with these previous studies and are consistent with the view that cognition is a dynamical process linked to neural dynamics18,19.

A previous study showed that sleep deprivation pulls the brain dynamics away from the criticality50. This result is consistent with ours because sleep deprivation generally compromises one’s cognitive and intellectual functions51.

Previous studies showed that the functional connectivity between particular pairs of ROIs or between subsystems of the brain in the resting state was correlated with intellectual ability49,52. These previous results are consistent with ours in the sense that the SG susceptibility can be regarded as the square sum of a type of functional connectivity over the pairs of ROIs and the intellectual score was positively correlated with the SG susceptibility in our analysis. In contrast to these previous studies, which looked at individual connectivity between two regions or subnetworks, we considered N = 264 ROIs scattered over the brain41 as a single functional network. We took this approach for two reasons. First, intelligence is considered to depend on large-scale brain networks26,52,53,54. Second, phase diagram analysis ideally requires a thermodynamic limit, i.e., infinitely many ROIs. One strategy to further approach the thermodynamic limit is to use a single voxel acquired by MRI as a node, significantly scaling up N. In this case, spatial correlation among ROIs, which we have ignored in the present study, would be prominent. Because the spatial dimensionality affects the phase diagrams even qualitatively30, this case may require two- or three-dimensional SG models. We leave this as a future problem. The literature also suggest that specific brain systems such as the fronto-parietal network55 and the default-mode network56 predict intelligence of humans. Running the same analysis for these and other brain systems to seek specificity of the results warrants future work. Because the present method requires hundreds of ROIs, we may benefit from considering voxel-wise networks of a specific brain system that allow many ROIs for particular brain systems.

In our previous paper, we posed the limited accuracy of fitting the PMEM to fMRI data when N is large38. The argument was based on the probability that each of the 2N possible activity patterns appears compared between the empirical data and the estimated PMEM. In the present manuscript, we have not used this accuracy measure, because it cannot be calculated when N is large. Instead, we validated the model by confirming that the difference between the empirical data and estimated PMEM in terms of the signal average, 〈Si〉, and the pairwise correlation, 〈SiSj〉, is small (Supplementary Fig. 8). This approach is based on the assumption that the average and second order correlation of signals explain most of the information contained in the given data, which has been confirmed for smaller N in previous studies using fMRI data21,24,37,38. Although only comparing 〈Si〉 and 〈SiSj〉 between the data and model is a weaker notion of accuracy of fit than using the accuracy measure38, the former approach has widely been accepted, explicitly or implicitly, in the literature57,58. However, we point out that how to justify the use of PMEMs when N is large remains an open issue.

There are various types of criticality, corresponding to different types of phase transitions. Within the framework of the Ising model, we showed that human fMRI data were in the paramagnetic phase and were close to the boundary with the SG phase but not to the boundary with the ferromagnetic phase. Furthermore, high fluid intelligence was associated with the proximity to the boundary between the paramagnetic and SG phases. In theory, the SG phase yields chaotic dynamics in spin systems including the SK model31,32,33, whereas the ferromagnetic phase is obviously non-chaotic. Therefore, although the definition of the chaos in the SG phase is different from that observed in cellular automata15 and recurrent neural networks16,17, our results are consistent with the idea of enhanced computational performance at the edge of chaos.

The previous accounts of the critical brain or critical neural circuits are mostly concerned with phase transitions different from the paramagnetic–SG phase transition or its analogs. Examples include phase transitions between quiescent (i.e., subcritical) and active (i.e., supercritical) phases as an excitability control parameter changes11,12,59,60,61, between ordered and chaotic phases as connectivity parameters change17, between a low-activity monostable state and a high-activity multistable state62, and the divergence of heat capacity5,6,35,36. Note that, in the theory of the Ising models, the heat capacity diverges on the boundary between the paramagnetic and ferromagnetic phases, whereas it increases without diverging on the boundary between the paramagnetic and SG phases30. Most of these previous results based on the Ising model related neural dynamics to the paramagnetic–ferromagnetic phase transition rather than the paramagnetic–SG transition. Roughly speaking, paramagnetic and ferromagnetic phases correspond to active and quiescent phases, respectively. Computational studies also support the ferromagnetism13,63,64. In contrast, we provided a signature of the paramagnetic–SG phase transition, not the paramagnetic-ferromagnetic transition. Fraiman et al. reported that the Ising model at the paramagnetic–ferromagnetic phase transition explains properties of functional networks based on fMRI data63. They used a two-dimensional Ising model with a uniform strength of interaction between pairs of nodes that are adjacent on a square lattice (and Jij = 0 for the rest of pairs). Another study that suggested the paramagnetic–ferromagnetic phase transition for fMRI signals also assumed a uniform Jij64. In contrast, we did not constrain the Jij values and instead inferred the Jij values (i.e., structure of functional network) using the PMEM. Because these previous studies63,64 did not assume heterogeneity in Jij as we did, their results do not contradict ours. In fact, the assumption of a uniform Jij corresponds to setting σ = 0 in our phase diagrams. If one varies μ under the condition σ = 0, the only possible phase transition is the paramagnetic–ferromagnetic transition (Fig. 1a–d). However, that phase transition point, which is derived under the condition σ = 0, is far from the location of the empirical data when σ is allowed to deviate from 0 (crosses in Fig. 1a–d). Therefore, allowing heterogeneity in Jij may be key to further clarifying the nature of critical neural dynamics.

We showed that neural dynamics for each participant were close to but substantially off the criticality separating the paramagnetic and SG phases. Other studies using the PMEM65 and other models66 also support off-critical as opposed to critical neural dynamics in the brain. A study applying the PMEM to local field potentials suggested that such off-critical dynamics may potentially have functional advantages because the off-critical situation would prevent the dynamics to get past the phase boundary to enter the other phase under the presence of noise66. The other phase may correspond to pathological neural dynamics such as epilepsy. The off-critical neural dynamics that we found for our participants, regardless of their IQ scores, may benefit from the same functional advantage.

Applying the current analysis pipeline to various neuroimaging and electrophysiological data in different contexts, from health to disease, and during rest and tasks, to evaluate the relevance of the different types of phase transitions warrants future work. For example, as a disease progresses, the brain dynamics may be gradually altered to transit from one phase to another, or to approach or repel from a phase transition curve. In fact, the method is applicable to general multivariate time series. Deployment of the present method to other biological and nonbiological data may also be productive.

One could classify the data from participants with high and low IQ scores using a simple multivariate Gaussian decoder67. Such a decoder would assume as input the mean and covariance of the fMRI data for each participant or its random samples having the same mean and covariance. In fact, multivariate Gaussian distributions having the same covariance structure as the empirical data yielded similar results (Supplementary Fig. 9). Because our PMEM also assumed the same input but was not optimized for classifying the participants, an optimized Gaussian decoder will probably be more efficient than our PMEM in explaining the IQ scores of the participants. This approach is conceptually much simpler than the present one, which employ the PMEM and its phase diagrams. However, the aim of the present study was to find empirical support of the critical brain hypothesis by relating the fMRI data to the phase diagrams of an archetypal spin system rather than to efficiently classify participants.

We found that the SG susceptibility was positively, although not strongly, correlated with individual differences in the performance IQ score but not in the verbal IQ score. The verbal IQ reflects individuals’ knowledge about verbal concepts and crystalized intelligence43; crystalized intelligence refers to one’s cognitive functioning associated with previously acquired knowledge and skills. In contrast, the performance IQ reflects fluid intelligence, which refers to active or effortful problem solving and maintenance of information39. Our results imply that the critical brain dynamics may be particularly useful for active and flexible cognitive functions.

Methods

Participants

One-hundred thirty eight (n = 138) healthy and right-handed adult participants (54 females and 84 males) in the Nathan Kline Institute’s (NKI) Rockland phase I Sample68 were analyzed. The data collection was approved by the institutional review board of the Nathan Kline Institute (no. 226781). Written informed consent was obtained for all the participants. Although the data set contains a wide range of the age (18–85 yo), the present results were not an age effect because the IQ values are standardized for age42 and because we have partialed out the effect of age (and the gender) in the present analysis. Participants’ IQ scores were derived from the Wechsler Abbreviated Scale of Intelligence42. We used the full scale IQ (full IQ for short), performance IQ, and verbal IQ.

Preprocessing

We used the same MRI data and the same preprocessing pipeline as our previous study’s69, except that we used resting-state fMRI signals from 264 ROIs, whose coordinates were derived in the previous literature41. In short, we submitted the resting-state fMRI data in the NKI Rockland phase I Sample with TR = 2500 ms and for 10 min 55 s for each participant to our preprocessing pipeline in FSL and applied band-pass temporal filtering (0.01–0.1 Hz).

The obtained fMRI signals xi(t) (\(i=1,\ldots ,N;t=1,\ldots ,{t}_{\max }\), where \({t}_{\max }=258\)) were transformed into their z-values using zi(t) = (xi(t) − μ(x(t)))∕σ(x(t)), where μ(x(t)) and σ(x(t)) represent the average and standard deviation of xi(t) over the N ROIs, respectively. Note that μ(x(t)) is the global signal. When we tested the robustness of the results by not removing the global signal, we set zi(t) = xi(t). We binarized the signal as follows:

$${S}_{i}(t)=\left\{\begin{array}{ll}\! +1 & {\rm{if}}\ {z}_{i}(t)\ge 0,\\ \! -1 & {\rm{if}}\ {z}_{i} (t)\, <\, 0.\hfill \end{array}\right.$$
(4)

Estimation of h and J by pseudo-likelihood maximization

The probability of each of the 2N activity patterns is equal to its frequency of occurrence normalized by the \({t}_{\max }\) time points and 138 participants. We fitted the Ising model to this probability distribution on the 2N activity patterns.

We estimated the parameter values of the Ising model (i.e., h and J) by maximizing a pseudo-likelihood (PL)38,70. We approximate the likelihood function by

$${\mathcal{L}}({\bf{h}},{\bf{J}})\approx \prod _{t=1}^{{t}_{\max }}\prod _{i=1}^{N}\tilde{P}({S}_{i}| {\bf{h}},{\bf{J}},{{\bf{S}}}_{/i}(t)),$$
(5)

where \(\tilde{P}\) represents the conditional Boltzmann distribution for a single spin, Si ∈ { −1, 1}, when the Sj values (j ≠ i) are equal to S∕i(t) ≡ (S1(t), …, Si−1(t), Si+1(t), …, SN(t)), i.e.,

$$\tilde{P}({S}_{i}| {\bf{h}},{\bf{J}},{{\bf{S}}}_{/i}(t))=\frac{\exp \left[{h}_{i}{S}_{i}+\sum _{j=1,\,j\ne i}^{N}{J}_{ij}{S}_{i}{S}_{j}(t)\right]}{\sum _{{S}_{i}^{\prime}=-1,+1}\exp \left[{h}_{i}{S}_{i}^{\prime}+\sum _{j =1,\,j\ne i}^{N}{J}_{ij}{S}_{i}^{\prime}{S}_{j}(t)\right]}.$$
(6)

In Eq. (6), one determines the probability of each activity pattern under the assumption that Sj (j ≠ i) does not change when drawing the value of Si (i = 1, ⋯ , N). We ran a gradient ascent updating scheme given by

$${h}_{i}^{{\rm{new}}}-{h}_{i}^{{\rm{old}}}=\epsilon \left({\langle {S}_{i}\rangle }_{{\rm{empirical}}}-{\langle \overline{{S}_{i}}\rangle }_{\tilde{P}}\right),$$
(7)
$${J}_{ij}^{{\rm{new}}}-{J}_{ij}^{{\rm{old}}}=\epsilon \left({\langle {S}_{i}{S}_{j}\rangle }_{{\rm{empirical}}}-{\langle \overline{{S}_{i}{S}_{j}}\rangle }_{\tilde{P}}\right),$$
(8)

where \({\langle \overline{{S}_{i}}\rangle }_{\tilde{P}}\) and \({\langle \overline{{S}_{i}{S}_{j}}\rangle }_{\tilde{P}}\) are the mean and correlation with respect to distribution \(\tilde{P}\) (Eq. (6)) and given by

$${\langle \overline{{S}_{i}}\rangle }_{\tilde{P}}=\frac{1}{{t}_{\max }}\sum _{t=1}^{{t}_{\max }}\tanh \left[{h}_{i}+\sum _{j^{\prime} =1,\,j^{\prime} \ne i}^{N}{J}_{i{j}^{\prime}}{S}_{{j}^{\prime}}(t)\right]$$
(9)

and

$${\langle \overline{{S}_{i}{S}_{j}}\rangle }_{\tilde{P}}=\frac{1}{{t}_{\max }}\sum _{t=1}^{{t}_{\max }}{S}_{j}(t)\tanh \left[{h}_{i}+\sum _{j^{\prime} =1,\,j^{\prime} \ne i}^{N}{J}_{i{j}^{\prime}}{S}_{{j}^{\prime}}(t)\right],$$
(10)

respectively. It should be noted that this updating rule avoids the calculation of 〈Si〉 and 〈SiSj〉 with the original spin system, Eqs. (1) and (2), which is computationally formidable with N = 264. As \({t}_{\max }\to \infty\), the estimator obtained by the PL maximization approaches the exact maximum likelihood estimator70. In fact, the Ising model with the estimated parameter values \({\bf{h}}=\hat{{\bf{h}}}\) and \({\bf{J}}=\hat{{\bf{J}}}\) produced the mean and correlation of spins in the empirical data with a sufficiently high accuracy (Supplementary Fig. 8).

We previously provided MATLAB code for estimating the Ising model from data by PL maximization38. The code is publicly available on GitHub repository (https://github.com/tkEzaki/energy-landscape-analysis).

Monte Carlo simulation

We used the Metropolis method71 to calculate the observables of the Ising model estimated from the empirical data and the SK model. In each time step, a spin Si was chosen uniformly at random for being updated. The selected spin was flipped with probability \(\min \{{e}^{-\Delta E},1\}\), where ΔE = E(Sflipped) − E(S), S is the current spin configuration, and Sflipped is the spin configuration after Si is flipped. The initial condition was given by Si = 1 with probability 1/2 (and hence Si = −1 with probability 1/2), independently for different i’s. We recorded the spin configuration S every N time steps.

For the empirical data, we discarded the first 106 × N time steps as transient and then recorded 107 samples of S in total. Based on the 107 samples, we calculated the averages of the observables (i.e., ∣m∣, q, χSG, χuni, and C). For drawing the phase diagrams with the N = 264 ROIs, we further averaged each observable over ten independent simulations starting from different initial spin configurations. In Fig. 1k, we averaged the χSG value over 40 combinations of \(N^{\prime}\) ROIs out of the 264 ROIs as well as over 107 samples and ten initial conditions.

For the phase diagram for the SK model, we discarded the first 104 × N time steps as transient and then collected 5 × 104 samples of S from each of 103 realizations of J. We drew the phase diagrams on the basis of the 5 × 104 × 103 = 5 × 107 samples.

Estimation of μ and σ for single participants

The estimation of the empirical interaction matrix, \(\hat{{\bf{J}}}\), requires a large amount of data, or practically, concatenation of fMRI data across different participants. Therefore, one cannot directly compute the mean and standard deviation of \(\hat{{\bf{J}}}\) (i.e., μ and σ) for each participant. Given this constraint, we estimated μ and σ for each participant (denoted by \(\tilde{\mu }\) and \(\tilde{\sigma }\)) using the χSG and χuni values for the participant (denoted by \({\tilde{\chi }}_{{\rm{SG}}}\) and \({\tilde{\chi }}_{{\rm{uni}}}\)) as follows (Supplementary Fig. 10).

First, we examined the phase diagrams in terms of χSG and χuni generated for the collection of all participants (Fig. 1c, d). Specifically, we calculated χSG(μ, σ) and χuni(μ, σ) values at μ = μk (k = 1, …, 25), where μ1 = −0.002, μ2 = −0.0015, …, μ25 = 0.01, and σ = σℓ (ℓ = 1, …, 21), where σ1 = 0, σ2 = 0.0075, …, σ21 = 0.15.

Second, at each μk(k = 1, …, 25), we computed the value of \({\check{\sigma }}_{k}\) satisfying \({\chi }_{{\rm{SG}}}({\mu }_{k},{\check{\sigma }}_{k})={\tilde{\chi }}_{{\rm{SG}}}\) (Supplementary Fig. 10a, c) using a linear interpolation of χSG(μk, σℓ) (ℓ = 1, …, 21), i.e., \({\check{\sigma }}_{k}=\alpha {\sigma }_{\ell ^{\prime} }+(1-\alpha ){\sigma }_{\ell ^{\prime} +1}\), where \(\ell ^{\prime}\) (\(1\le \ell ^{\prime} \le 21\)) is the integer satisfying \({\chi }_{{\rm{SG}}}({\mu }_{k},{\sigma }_{\ell ^{\prime} })\le {\tilde{\chi }}_{{\rm{SG}}}\, <\, {\chi }_{{\rm{SG}}}({\mu }_{k},{\sigma }_{\ell ^{\prime} +1})\), and \(\alpha =[{\chi }_{{\rm{SG}}}({\mu }_{k},{\sigma }_{\ell ^{\prime} +1})-{\tilde{\chi }}_{{\rm{SG}}}]/[{\chi }_{{\rm{SG}}}({\mu }_{k},{\sigma }_{\ell ^{\prime} +1})-{\chi }_{{\rm{SG}}}({\mu }_{k},{\sigma }_{\ell ^{\prime} })]\). Because χSG(μk, σℓ) increases with ℓ in the paramagnetic phase, the \(\ell ^{\prime}\) value is uniquely determined for each k, if it exists. In this manner, we obtained a piecewise linear curve whose knots were (\({\mu }_{k},{\check{\sigma }}_{k}\)) (k = 1, …, 25). On this curve, χSG(μ, σ) is approximately equal to \({\tilde{\chi }}_{{\rm{SG}}}\) (Supplementary Fig. 10e, g). It should be noted that we have assumed that \((\tilde{\mu },\tilde{\sigma })\) to be estimated is near \((\hat{\mu },\hat{\sigma })\) computed for the entire population (represented by the cross in Fig. 1a–d). More precisely, we are searching \((\tilde{\mu },\tilde{\sigma })\) in the vicinity of the paramagnetic–SG phase boundary on the paramagnetic side. This assumption is supported by the empirical values of m and q for individual participants, i.e., m = −8.0 × 10−3 ± 7.8 × 10−3 (mean  ±  SD) and q = 3.4 × 10−3 ± 0.4 × 10−3.

Third, we calculated a piecewise linear curve on which χuni(μ, σ) was approximately equal to \({\tilde{\chi }}_{{\rm{uni}}}\) (Supplementary Fig. 10f, g). To this end, we applied the same algorithm as the one used in the previous step but by fixing σℓ (Supplementary Fig. 10b) and finding \({\check{\mu }}_{\ell }\) (Supplementary Fig. 10d), exploiting the fact that χuni(μk, σℓ) monotonically increases with μ in the paramagnetic phase.

Finally, we computed (\(\tilde{\mu }\), \(\tilde{\sigma }\)) for the individual as the intersection of the two piecewise linear curves (Supplementary Fig. 10g).

Specific heat

The specific heat is defined by

$$C=\frac{\langle {E}^{2}\rangle -{\langle E\rangle }^{2}}{N{T}^{2}},$$
(11)

where T is the temperature. We set T = 1 because we implicitly did so in Eqs. (1) and (2).

To compute C for each participant, we first drew a phase diagram for C in terms of μ and σ for the entire population (Supplementary Fig. 11a). The obtained phase diagram was similar to that for the SK model (Supplementary Fig. 11b). Then, we determined the C value for each participant as the point in the phase diagram corresponding to the (μ, σ) for the participant. Because the phase diagram for C is drawn for discrete values of μ and σ, we applied the standard bilinear interpolation to determine the C value corresponding to a given (μ, σ).

Statistics and reproducibility

Statistical tests were performed using SPSS 24.0. The details of each analysis are found in prior sections.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.