Critical brain dynamics and human intelligence

According to the critical brain hypothesis, the brain is considered to operate near criticality and realize efficient neural computations. Despite the prior theoretical and empirical evidence in favor of the hypothesis, no direct link has been provided between human cognitive performance and the neural criticality. Here we provide such a key link by analyzing resting-state dynamics of functional magnetic resonance imaging (fMRI) networks at a whole-brain level. We develop a novel data-driven analysis method, inspired from statistical physics theory of spin systems, to map out the whole-brain neural dynamics onto a phase diagram. Using this tool, we show evidence that dynamics of more intelligent human participants are closer to a critical state, i.e., the boundary between the paramagnetic phase and the spin-glass (SG) phase. This result was specific to fluid intelligence as opposed to crystalized intelligence. The present results are also consistent with the notion of “edge-of-chaos” neural computation. Author summary According to the critical brain hypothesis, the brain should be operating near criticality, i.e., a boundary between different states showing qualitatively different dynamical behaviors. Such a critical brain dynamics has been considered to realize efficient neural computations. Here we provide direct neural evidence in favor of this hypothesis by showing that the brain dynamics in more intelligent individuals are closer to the criticality. We reached this conclusion by deploying a novel data-analysis method based on statistical-physics theory of Ising spin systems to functional magnetic resonance imaging (fMRI) data obtained from human participants. Specifically, our method maps multivariate fMRI data obtained from each participant to a point in a phase diagram akin to that of the Sherrington-Kirkpatrick model of spin-glass systems. The brain dynamics for the participants having high fluid intelligence scores tended to be close to the phase boundary, which marks criticality, between the paramagnetic phase and the spin-glass phase, but not to the boundary between the paramagnetic and ferromagnetic phases.


Introduction
Critical brain hypothesis posits that the brain operates near a critical regime, i.e., boundary between different phases showing qualitatively different behaviors [1][2][3][4][5] . This hypothesis has been investigated for more than two decades including criticisms such as the lack of power laws in the relevant observables [6][7][8][9] . Theoretical and experimental work has shown that neural systems operating near criticality are advantageous in information transmission, information storage, classification, and nonlinear input filtering 1, 4, 10-12 . These findings align with the idea of edge-of-chaos computation, with which computational ability of a system is maximized at criticality separating a chaotic phase and a non-chaotic phase [13][14][15] . These findings are also in line with a general contention that cognitive computations occur as neural dynamical processes 16,17 .
According to the critical brain hypothesis, 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 rest [18][19][20] , working memory tasks 21 and visual perception tasks 22 . On these grounds, in the present study we hypothesize that complex and transitory neural dynamics 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 avalanches 10,23 , including the case of humans 4,24,25 . However, neural avalanches do not imply transitory dynamics or their absence. Second, nonlinear time series analysis has found signatures of chaotic dynamics in electroencephalography (EEG) signals recorded from the brains of healthy controls more strongly than from those of patients with, for example, epilepsy, Alzheimer's disease, and schizophrenia 26 . However, this method 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 fact, many cognitive functions seem to depend on network connectivity among various regions scattered over the whole brain 27 . State-transition dynamics in the brain are also considered to involve large-scale brain networks 16,17,28 .
In the present study, we develop a novel method to measure the extent to which neural dynamics obtained from large-scale brain networks are close to criticality and complex transitory 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 transitory dynamics) 29 and chaotic dynamics [30][31][32] . 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 data 5,[33][34][35] and fMRI data recorded from a collection of regions of interest (ROIs) 18,19,22,36,37 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 novel method to resting-state fMRI data recorded from human adults with different degrees of intelligence. 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 information 38 . Fluid intelligence is strongly related to the general intelligence factor, g 38 and predictive of real-world outcomes such as job performance 39 . We examine our hypothesis that large-scale brain dynamics of more intelligent human individuals in the sense of fluid intelligence are closer to critical.

Results
Human brain dynamics are in the paramagnetic phase and close to the spin-glass phase 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 Hamiltonian of the Ising model are inferred from data. To fit the model, we binarized the resting-state fMRI signals obtained from 138 healthy adults (see Methods). The binarized activity pattern at N(= 264) ROIs 40 at time t (t = 1, . . . ,t max ; t max = 258) is denoted by S S S(t) = (S 1 (t), . . . , S N (t)) ∈ {−1, +1} N , where S i (t) = 1 and S i (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) 22,33 : In (1), is the energy of activity pattern We denote the estimated parameter values byĥ h h andĴ.
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 h h atĥ h h following the theoretical convention 29 , including when the PMEM is applied to data analysis 5 . We set h h h =ĥ h h also because changing the h h h values did not qualitatively change the phase diagrams (Fig. S1). Then, motivated by the investigation of the archetypical Sherrington-Kirkpatrick (SK) model of spin-glass systems 29 , we varied the mean µ and standard deviation σ of J by linearly transforming J, i.e., In (3),μ = 1.57 × 10 −3 andσ = 3.57 × 10 −2 are the mean and standard deviation of the off-diagonal elements ofĴ estimated for the empirical data. Given that the effect of changing h h h is negligible, a large σ corresponds to a low temperature of a spinglass model 29 . For each set of J i j 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 σ . The phase diagrams for the magnetization, m (= ∑ 1≤i≤N ⟨S i ⟩/N), and the spin-glass order parameter, q (= ∑ 1≤i≤N ⟨S i ⟩ 2 /N), where ⟨·⟩ represents the ensemble average, are shown in Figs. 1A and 1B, 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 J i j (= J ji , i ̸ = j) being independently drawn from a Gaussian distribution with mean µ and standard deviation σ (Figs. 1E and 1F). The parameter space is composed of the paramagnetic phase (m = 0, q = 0), ferromagnetic phase (m ̸ = 0, q > 0), and spin-glass (SG) phase (m = 0, q > 0) 29 , although 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 transition to the SG phase (crosses in Figs. 1A and 1B). In theory, the spin-glass susceptibility, diverges on the boundary between the paramagnetic and SG phases 29 . 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 β ∑ 1≤i, j≤N c i j , which characterizes the transition between the paramagnetic and ferromagnetic phases 29 (Fig. 1D). Note that the phase diagrams for χ SG and χ uni resemble those obtained from the SK model (Figs. 1G and 1H).
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 µ =μ (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 the 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 ′ (< N) out of the N ROIs uniformly at random and fitted the PMEM. The estimated parameter values are denoted byĥ h h and J without confusion. Then, we simulated the equilibrium state of the system by scanning J according to (3), where we varied σ while fixing µ =μ. In this manner, we sought to investigate how the data were close to the SG phase transition at each N ′ value. As shown in Fig. 1K, the peak value of χ SG increased as N ′ 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,σ , as N increased. By regressing σ peak /σ linearly on 1/N ′ (inset of Fig. 1K), we estimated σ peak /σ = 1.45 ± 0.04 in the limit N ′ → ∞.

The performance IQ but not verbal IQ is associated with the transition between the paramagnetic and spin-glass phases
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 (Figs. 1C and 1I), and the performance intelligence quotient (IQ) score. The performance IQ score is defined based on tasks that are reflective of fluid intelligence 41,42 . 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 σ (t 129 = 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 (t 127 = 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 µ (t 129 = 0.77, P = 0.44, d = 0.13 with the outlier included; t 127 = 0.85, P = 0.40, d = 0.15 with the outlier excluded).
More systematically, we found a significantly positive correlation between χ SG and the performance IQ score (r 129 = 0.24, P Bonferroni = 0.011; also see Fig. 2B). However, the verbal IQ score, which is based on individuals' verbal knowledge 41,42 , was not correlated with χ SG (r 126 = 0.06, P uncorrected = 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 (t 121 = 2.33, P = 0.021, in the Williams t-test for comparing two nonindependent correlations with a variable in common 43 ). 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 effect of the age and gender in this and the following analysis unless we state otherwise.
The correlation between the full IQ score 41,42 and χ SG was intermediate between the results for the performance and verbal IQ scores (r 130 = 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 was robust in the following four ways. First, the exclusion of the two outliers determined by Tukey's 1.5 criteria 44 did not affect the significance of the results (χ SG vs performance IQ: r 127 = 0.27, P Bonferroni = 0.005; χ SG vs verbal IQ: r 124 = 0.13, P Bonferroni = 0.27; χ SG vs full IQ: r 128 = 0.25, P = 0.005). Second, the results were robust against variation on the threshold value for binarizing the fMRI signal (Fig. S2). Third, the results were preserved even when the global signal (see Methods) was not subtracted from the fMRI signals (χ SG vs performance IQ: r 129 = 0.22, P Bonferroni = 0.02; χ SG vs verbal IQ: r 126 = 0.046, P uncorrected = 0.61; χ SG vs full IQ: r 130 = 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 coefficients 45 ; 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.

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 autocorrelation 46,47 . This measure quantifies the "critical slowing down" phenomenon, which has been observed in critical states of the brain 47 . 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 analysis 46,47 . 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: r 129 = 0.29, P Bonferroni = 0.002; α vs verbal IQ: r 126 = 0.19, P Bonferroni = 0.068; α vs full IQ: r 130 = 0.25, P = 0.003). These results were robust against the removal of outliers (α vs performance IQ: r 128 = 0.28, P Bonferroni = 0.002; α vs verbal IQ: r 125 = 0.17, P Bonferroni = 0.10; α vs full IQ: r 130 = 0.25, P = 0.003).
We then performed a multivariate linear regression of the performance IQ with χ SG and α being the independent variables. We found a significant regression equation (F 2,128 = 8.0, P < 0.001, adjusted R 2 = 0.11). Both χ SG and α were significantly correlated with the performance IQ (χ SG : β = 0.18, P = 0.039; α: β = 2.4, 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.

Discussion
We provided neural support that more intelligent human brains 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.
The criticality view of the brain itself is not new. 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 regions 48 , and the structure of brain networks 27,48 . The present 4/12 results derived from dynamic fMRI signals provide an orthogonal account of human intelligence as compared to these previous studies and are consistent with the view that cognition is a dynamical process linked to neural dynamics 16,17 .
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 intelligent performance 48,49 . 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 intelligence 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 brain 40 as a single functional network. We took this approach for two reasons. First, cognition is considered to depend on large-scale brain networks 27 . 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 qualitatively 29 , this case may be require twoor three-dimensional SG models. We leave this as a future problem.
Balanced excitation and inhibition is another key mechanism leveraging neural computation. It posits that the strength of excitatory and inhibitory connection among neurons is balanced at a network level and that this balance enhances neural computation including the case of critical-state neural dynamics 50,51 . Our results are consistent with balanced excitation and inhibition. Our parameter J i j is a measure of functional connectivity 36 , although it is not direct evidence of excitatory or inhibitory synaptic connectivity. The SG phase is characterized by a large variance of J i j across 1 ≤ i ̸ = j ≤ N as compared to its mean (Fig. 1). Therefore, the present results can be also interpreted as a positive correlation between the intelligence and the balanced excitation/inhibition among ROIs.
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 model [30][31][32] , 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 automata 13 and recurrent neural networks 14,15 , 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 analogues. Examples include phase transitions between quiescent (i.e., subcritical) and active (i.e., supercritical) phases as an excitability control parameter changes 10,23,50,52,53 , between ordered and chaotic phases as connectivity parameters change 15 , between a low-activity monostable state and a high-activity multistable state 54 , and the divergence of heat capacity 4,5,34,35 . 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 phases 29 . Most of these previous results are probably related to the paramagnetic-ferromagnetic phase transition, or generally the transition between a quiescent phase and an active phase, rather than the paramagnetic-SG transition. Computational studies also support the ferromagnetism of fMRI data 11,55,56 . In contrast, we provided a signature of the paramagnetic-SG phase transition, not the paramagnetic-ferromagnetic transition. This discrepancy may be partly caused by the data; the previous research mostly used data where each node is active with a small probability (e.g., spikes in electrophysiological data), whereas the present research used fMRI data.
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 non-biological data may also be productive.
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 intelligence 42 ; 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 information 38 . Our results imply that the critical brain dynamics may be particularly useful for active and flexible cognitive functions.

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 Sample 57 were analyzed. Although the data set contains a wide range of the age (18-85 yo),

5/12
the present results were not an age effect because the IQ values are standardized for age 41 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 Intelligence (WASI) 41 . 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's 58 , except that we used resting-state fMRI signals from 264 ROIs, whose coordinates were derived in the previous literature 40 . In short, we used the data from the Nathan Kline Institute's Rockland phase I Sample 57 . Then, we submitted the resting-state fMRI data with TR = 2500 ms and for 10 m 55 s to our preprocessing pipeline in FSL and applied band-pass temporal filtering (0.01-0.1 Hz).
The obtained fMRI signals x i (t) (i = 1, . . . , N;t = 1, . . . ,t max , where t max = 258) were transformed into their z-values using z i (t) = (x i (t) − µ(x(t)))/σ (x(t)), where µ(x(t)) and σ (x(t)) represent the average and standard deviation of x i (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 z i (t) = x i (t). We binarized the signal as follows:

Estimation of h h h and J by pseudo-likelihood maximization
The probability of each of the 2 N 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 2 N activity patterns. We estimated the parameter values of the Ising model (i.e., h h h and J) by maximizing a pseudo-likelihood (PL) 37,59 . We approximate the likelihood function by whereP represents the conditional Boltzmann distribution for a single spin, As t max → ∞, the estimator obtained by the PL maximization approaches the exact maximum likelihood estimator 59 . We ran a gradient ascent updating scheme given by where ⟨S i ⟩P and ⟨S i S j ⟩P are the mean and correlation with respect to distributionP (Eq. (6)) and given by and respectively. It should be noted that this updating rule avoids the calculation of ⟨S i ⟩ and ⟨S i S j ⟩ with the original spin system, Eqs. (1) and (2), which is computationally formidable with N = 264. The Ising model with the estimated parameter values h h h =ĥ h h and J =Ĵ produced the mean and correlation of spins in the empirical data with a sufficiently high accuracy (Fig. S3).

Monte Carlo simulation
We used the Metropolis method 60 to calculate the observables of the Ising model estimated from the empirical data and the SK model. In each time step, a spin S i was chosen uniformly at random for being updated. The selected spin was flipped with probability min{e −∆E , 1}, where ∆E = E(S S S flipped ) − E(S S S), S S S is the current spin configuration, and S S S flipped is the spin configuration after S i is flipped. The initial condition was given by S i = 1 with probability 1/2 (and hence S i = −1 with probability 1/2), independently for different i's. We recorded the spin configuration S S S every N time steps.
For the empirical data, we discarded the first 10 6 × N time steps as transient and then recorded 10 7 samples of S S S in total. Based on the 10 7 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 10 independent simulations starting from different initial spin configurations. In Fig. 1K, we averaged the χ SG value over 40 combinations of N ′ ROIs out of the 264 ROIs as well as over 10 7 samples and 10 initial conditions.
For the phase diagram for the SK model, we discarded the first 10 4 × N time steps as transient and then collected 10 5 samples of S S S from each of 10 3 realizations of J. We drew the phase diagrams on the basis of the 5 × 10 4 × 10 3 = 5 × 10 7 samples.

Estimation of µ and σ for single participants
The estimation of the empirical interaction matrix,Ĵ, 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Ĵ (i.e., µ and σ ) for each participant. Given this constraint, we estimated µ and σ for each participant (denoted byμ andσ ) using the χ SG and χ uni values for the participant (denoted byχ SG andχ uni ) as follows (Fig. S4).
Third, we calculated a piecewise linear curve on which χ uni (µ, σ ) was approximately equal toχ uni (Figs. S4F and S4G). To this end, we applied the same algorithm as the one used in the previous step but by fixing σ ℓ (Fig. S4B) and findingμ ℓ (Fig. S4D), exploiting the fact that χ uni (µ k , σ ℓ ) monotonically increases with µ in the paramagnetic phase.
Finally, we computed (μ,σ ) for the individual as the intersection of the two piecewise linear curves (Fig. S4G).

Specific heat
The specific heat is defined by 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 (Fig. S5A). The obtained phase diagram was similar to that for the SK model (Fig. S5B). 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 (µ, σ ). , the crosses represent the mean and standard deviation of the J i j estimated for the entire population of the participants, i.e., (μ,σ ). 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 µ =μ being fixed. (J) χ uni as a function of µ, with σ =σ 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 ′ is varied. The value of σ = σ peak that maximizes χ SG is plotted against 1/N ′ in the inset. The dashed line is the linear regression based on the six data points, N ′ = 40, 60, 80, 120, 160, and 264. The coefficient of determination is denoted by R 2 .