Adrián Ponce-Alvarez

Lynn Uhrig

Nikolas Deco

Miguel Camilo

M L Signorelli

Béchir Kringelbach

Gustavo Jarraya

Deco

Camilo M Signorelli

Morten L Kringelbach

Béchir Jarraya

Gustavo Deco

Macroscopic quantities of collective brain activity during wakefulness and anesthesia

Keywords: 5 figures), Appendix, and Supplementary

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Figures S1 to S8.

Introduction

Interesting phenomena in biological systems are usually collective behaviors emerging from the interactions among many constituents. Large-scale brain activity is not an exception: the brain's network continuously generates coordinated spontaneous patterns of activity among brain regions at multiple spatiotemporal scales [START_REF] Biswal | Functional connectivity in the motor cortex of resting human brain using echo-planar MRI[END_REF][START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]. Changes in spontaneous brain activity are observed in different brain states, the study of which is essential to understand the organizing principles of brain activity.

For instance, anesthesia has been used to transiently induce loss of consciousness and to investigate the neural correlates of the awake state. Previous studies showed that different anesthetics, acting on different molecular targets (4), similarly impact the strength and the structure of functional correlations [START_REF] Ferrarelli | Breakdown in cortical effective connectivity during midazolam-induced loss of consciousness[END_REF][START_REF] Boly | Connectivity changes underlying spectral EEG changes during propofol-induced loss of consciousness[END_REF][START_REF] Guldenmund | Thalamus, brainstem and salience network connectivity changes during propofol-induced sedation and unconsciousness[END_REF][START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF][START_REF] Uhrig | Cerebral responses to local and global auditory novelty under general anesthesia[END_REF][START_REF] Greicius | Persistent default-mode network connectivity during light sedation[END_REF], and their dependence on interareal anatomical connections [START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF][START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. However, how changes in local regions and subnetworks combine to affect the collective brain dynamics and to lose consciousness remains largely unknown. To answer this question, it is essential to precisely characterize the collective properties of different brain states and their dependence on parameters at the system's level. This dependence is likely not straightforward since, as for many complex systems, the system's behavior could be differently affected by changes in its parameters. In such a case, while some parameters can largely vary without affecting the system's behavior (so-called "sloppy" parameters), even small changes in some others can significantly modify it [START_REF] Machta | Parameter space compression underlies emergent theories and predictive models[END_REF][START_REF] Panas | Sloppiness in spontaneously active neuronal networks[END_REF][START_REF] Ponce-Alvarez | Cortical state transitions and stimulus response evolve along stiff and sloppy parameter dimensions, respectively[END_REF].

In recent years, statistical mechanics has proven to be more and more useful to describe collective neural activity. Statistical mechanics shows that the behaviors of complex systems can be captured by macroscopic properties, which emerge from the collective activity of the units, in a way largely independent of the microscopic details of the system. These emergent (macroscopic) behaviors can be classified into qualitatively different ordered or disordered phases. Of particular interest are dynamics poised close to phase transitions, or critical points, where order and disorder coexist.

Theoretical reasoning shows that complex dynamics and optimal information processing are expected at critical points, making criticality a candidate unifying principle to account for the brain's inherent complexity necessary to process and represent its environment [START_REF] Chialvo | Emergent complex neural dynamics: the brain at the edge[END_REF][START_REF] Shew | Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches[END_REF][START_REF] Shew | The functional benefits of criticality in the cortex[END_REF][START_REF] Hidalgo | Information-based fitness and the emergence of criticality in living systems[END_REF]. Following this view, studies of whole-brain and local circuits dynamics have proposed that anesthesia shifts the dynamics from the critical point [START_REF] Tagliazucchi | Large-scale signatures of unconsciousness are consistent with a departure from critical dynamics[END_REF][START_REF] Fekete | Critical dynamics, anesthesia and information integration: Lessons from multi-scale criticality analysis of voltage imaging data[END_REF]. This is supported by the reduction of several measures of brain dynamics complexity under anesthesia [START_REF] Casali | A theoretically based index of consciousness independent of sensory processing and behavior[END_REF][START_REF] Hudetz | Dynamic repertoire of intrinsic brain states is reduced in propofol-induced unconsciousness[END_REF][START_REF] Schartner | Complexity of Multi-Dimensional Spontaneous EEG Decreases during Propofol Induced General Anaesthesia[END_REF][START_REF] Solovey | Loss of consciousness is associated with stabilization of cortical activity[END_REF]. The global mechanisms underlying different conscious states have been recently investigated using an anatomically-constrained dynamical model with a global coupling parameter in combination with EEG recordings [START_REF] Lee | Relationship of critical dynamics, functional connectivity, and states of consciousness in large-scale human brain networks[END_REF]. However, it remains unknown which are the macroscopic properties and the relevant local/global parameters describing the transition of collective activity from the awake to anesthetized states. Indeed, different local/global network parameters are likely to jointly determine the different brain states and to differently contribute to the state transitions.

In this study, we addressed these questions by analyzing the brain's collective activity in different levels of arousal, i.e., during wakefulness and under anesthesia. Specifically, we analyzed restingstate fMRI dynamics of awake and anesthetized macaque monkeys [START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. Five different anesthesia protocols, involving 3 different anesthetics (propofol, ketamine, and sevoflurane), were used to 3 induce moderate sedation or deep anesthesia. First, we derived efficient statistics that distinguished between awake and anesthetized brain states. Second, we used these statistics and the maximum entropy principle to model the brain's activity and to derive important emergent properties that described the different brain states. These emergent properties provided information about the system's physical state, and about its capability to produce work, to contain information and to transmit it. Finally, we investigated the dependence of collective activity on the different model parameters.

Results

We analyzed the resting-state fMRI dynamics of five rhesus macaques (Macaca mulatta) under different levels of arousal: wakefulness (n = 24 scans), two levels of propofol sedation (light, LPP, n = 21, and deep, DPP, n = 23), sedation through ketamine (KETA, n = 22), and two types of sevoflurane anesthesia (SEV2, n = 18, and SEV4, n = 11) (see Methods and Appendix). fMRI MION time-series were obtained for 𝑁 = 82 previously defined regions of interest (ROIs) (CoCoMac Regional Map parcellation). Each scan was 20 min long and was acquired in time frames of 2.4 s (i.e., 500 time frames).

Coupling to population reliably distinguished between awake and anesthetized brain states

We were interested on collective patterns displayed among the 𝑁 ROIs, for the six different experimental conditions. For this, we first binarized the z-scored time-series of each ROI, 𝑥 𝑖 (𝑡), by imposing a threshold 𝜃 = -1 (Fig. 1A,B, see Methods). Binarization of time-series has proven to effectively capture and compress fMRI large-scale dynamics [START_REF] Tagliazucchi | Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis[END_REF][START_REF] Watanabe | A pairwise maximum entropy model accurately describes resting-state human brain networks[END_REF]. We concentrated on different statistics that described the resulting binary data: the activation rate of each ROI, i.e., 〈𝜎 𝑖 〉, the correlation between ROIs, i.e., 𝐶 𝑖𝑗 = corr(𝜎 𝑖 , 𝜎 𝑗 ), and the population activity, i.e., 𝐾(𝑡) = ∑ 𝜎 𝑖 (𝑡) 𝑁 𝑖=1 (Fig. 1C). We were particularly interested on the coupling of each ROI to the population activity, defined as:

𝑧 𝑖 = corr[𝜎 𝑖 (𝑡), 𝐾 ≠𝑖 (𝑡)], (1)
where 𝐾 ≠𝑖 (𝑡) is the sum activity of all but ROI 𝑖: 𝐾 ≠𝑖 (𝑡) = ∑ 𝜎 𝑖 (𝑡)

𝑗≠𝑖

. Recent findings showed that propofol anesthesia affects the coupling to global signal in human and rats [START_REF] Tanabe | Altered global brain signal during physiologic, pharmacologic, and pathologic states of unconsciousness in humans and rats[END_REF]. In the following we showed that the statistics 𝒛 = [𝑧 1 , … , 𝑧 𝑁 ] provides, with only 𝑁 parameters, a compact description of the binary collective activity and can be used to classify the brain states. .

The couplings to the population were highly predictive of the functional correlations (Fig. 2A-C).

Indeed, the product 𝜂 𝑖𝑗 = 𝑧 𝑖 × 𝑧 𝑗 highly correlated with the functional correlation (FC) between the fMRI time-series of ROIs (𝑖, 𝑗) (corr.: 0.65-0.78, p < 0.001). Moreover, we found that the vector 𝒛 correlated across scans within the same experimental conditions, with the average correlation coefficient being equal to 0.3 ± 0.01 (Fig. 2D, blue distribution). This correlation was significantly higher (p < 0.001, F(2,3533) = 976.5, one-way ANOVA followed by Tukey's post hoc analysis) than those obtained using the vectors representing the average activities and the correlations, i.e., vectors

𝝁 = [〈𝜎 1 〉, … , 〈𝜎 𝑁 〉] and 𝝆 = [𝐶 1,2 , 𝐶 1,3 , … , 𝐶 𝑁-1,𝑁 ]
, respectively (Fig. 2D, green and red distributions, corr.: 0.06 ± 0.01 and 0.11 ± 0.01). Furthermore, we found that the coupling to the population could be used to classify the awake and anesthesia states with high accuracy (Fig. 3A). We tested this by using a classifier based on k-means clustering (see Methods). Based on the statistic 𝒛 we were able to classify the scans of two categories, awake vs. anesthesia (independently of the anesthetic), with 96.6% of correct classifications (chance level: 50%). This classification performance was higher than the one obtained using the statistics 𝝁 and 𝝆, yielding 74.0% and 74.8% of correct classifications, respectively (Fig. 3B,C). Classification among the six experimental conditions yielded lower performances but was higher for statistic 𝒛 than for 𝝁 and 𝝆: 38.5%, 33.1%, and 29.9%, respectively (chance level: 16.7%). Similar differences in classification performances for population couplings and functional correlations were obtained using continuous (not thresholded) signals (Fig. S1). Altogether, these results show that the coupling to population is a reliable marker to distinguish between awake and anesthetized brain states. To examine which ROIs contributed the most to distinguish between the awake state and anesthesia based on 𝒛, we performed PCA on the collection of z-scored vectors 𝒛. The first principal component was sufficient to separate the awake and anesthesia conditions (Fig. 3D). This component had strong coefficients for brain regions located in the cingulate, parietal, intraparietal, insular cortices, and the hippocampus (Fig. 3E). Overall, changes in average couplings to the population with respect to awake values were similar for all anesthetics (Fig. 3F). We next asked how these changes affect the collective properties of brain dynamics.

Modelling collective activity using maximum entropy models

Collective activity is ultimately described by the probability of each of the binary patterns 𝝈 = .

(

) 2
The model parameters 𝛼 𝑖 , 𝛽 𝐾 , and 𝛾 𝑖 are Lagrange multipliers associated to the constrained observables 〈𝜎 𝑖 〉, 𝑃(𝐾), and 〈𝜎 𝑖 𝐾〉, respectively. The normalizing constant 𝑍 is the partition function,

given by 𝑍 = ∑ 𝑒 -𝐸(𝝈) {𝝈} , which contains information about useful statistics predicted by the model (see below). This model can be extended to include the non-linear coupling between 𝜎 𝑖 and 𝐾. Indeed, the maximum entropy distribution that is consistent with the joint probability distributions of 𝜎 𝑖 and 𝐾, i.e., 𝑃(𝜎 𝑖 , 𝐾), yields the following energy function [START_REF] Gardella | A Tractable Method for Describing Complex Couplings between Neurons and Population Rate[END_REF]:

𝐸(𝝈) = -∑ 𝐻 𝑖,𝐾(𝝈) 𝜎 𝑖 𝑁 𝑖=1 , (3)
where 𝐾(𝝈) is the number of active ROIs in pattern 𝝈 and the parameters 𝐻 𝑖,𝐾(𝝈) represent the tendency of ROI 𝑖 to activate when 𝐾(𝝈) ROIs are active. For both linear and non-linear coupling-MEMs the model parameters were inferred from the data using maximum likelihood [START_REF] Gardella | A Tractable Method for Describing Complex Couplings between Neurons and Population Rate[END_REF]. Notably, for the coupling-MEMs the partition function can be calculated directlysomething that is generally not the case for most MEMs, since its calculation involves summing over all possible states.

We used these coupling-MEMs to fit the binary single-scan fMRI data for the different experimental conditions. The models accurately estimated the distribution of population activity 𝑃(𝐾) (average Jensen-Shannon divergence 𝐷 𝐽𝑆 between the model and data distributions: 𝐷 𝐽𝑆 < 10 -6 for both the non-linear and linear coupling-MEM; Fig. 4A and Fig. S2). Moreover, the models were able to moderately predict the covariances of the data (Fig. 4B,C), which were not used to constrain the models. Across the different datasets, the average correlation between the data and predicted covariances was 𝑟 = 0.28 ± 0.03 for the linear coupling model and reached 0.40 ± 0.02 for the nonlinear coupling model (see also Fig. S2). Furthermore, scan-classification based on parameters 𝛾 𝑖 yielded 86% and 45% correct classifications between awake and anesthetized conditions and among the six experimental conditions, respectively (Fig. S3A,B). Using parameters 𝛼 𝑖 the classifier performance decreased to 75% and 28%, respectively (Fig. S3B). Thus, the learned linear coupling-MEM showed consistent variations in parameters 𝛾 𝑖 (associated to 𝑧 𝑖 ) across the different arousal states.

Collective activity indicated reduced free energy, susceptibility, and heat capacity under anesthesia

We can learn interesting features of collective activity using the estimated models. One important quantity is the system's Helmholtz free energy, which is given by the difference between the average energy (〈𝐸〉) and the entropy (𝑆), i.e., 𝐹 = 〈𝐸〉 -𝑆. The free energy quantifies the useful energy that is obtainable from the system. Using the Boltzmann distribution, the free energy can be directly obtained from the partition function as 𝐹 = -ln(Z). Thus, since 𝑍 is tractable for the coupling-MEMs, we can directly estimate 𝐹. We found that the free energy was significantly higher for the awake state compared to all anesthetized conditions for both the non-linear (Fig. 4D) and the linear (Fig S4A)

coupling-MEM (p < 0.001, one-way ANOVA followed by Tukey's post hoc analysis). This result is both interesting and reasonable because it indicates that more useful energy can be extracted from the awake state than from the anesthetized state.

Two other important statistical quantities can be derived from the model, namely the susceptibility and the heat capacity. The susceptibility 𝜒 relates to the diversity of population states, while the heat capacity 𝐶 ℎ quantifies the diversity of accessible energy states. Specifically, the susceptibility and the heat capacity of the model are given by the variances of the population activity and the energy, respectively, i.e., 𝜒 = var(𝐾) and 𝐶 ℎ = var(𝐸). We found that 𝜒 and 𝐶 ℎ were significantly higher for the awake state compared to all anesthetized conditions for both linear and non-linear coupling-MEM (non-linear model: Fig. 4E,F, linear model: Fig S4B ,C; p < 0.001, one-way ANOVA followed by Tukey's post hoc analysis). This indicates that the awake system had larger population fluctuations and a larger repertoire of energy states than the system under anesthesia.

We next tested whether the same differences in these statistical quantities were found using MEM constraint by other statistics. To build the models we estimated the maximum entropy distribution 𝑃(𝝈) under the constraint that the activation rates (< 𝜎 𝑖 >) and the pairwise correlations (< 𝜎 𝑖 𝜎 𝑗 >) are preserved. The energy of the Boltzmann distribution that is consistent with these expectation values is given by [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF]:

𝐸(𝝈) = -∑ ℎ 𝑖 𝑁 𝑖=1 𝜎 𝑖 - 1 2 ∑ ∑ 𝐽 𝑖𝑗 𝜎 𝑖 𝜎 𝑗 𝑁 𝑗=1 𝑁 𝑖=1 . ( 4
)
In this pairwise-MEM, the parameter ℎ 𝑖 , called intrinsic bias, represents the intrinsic tendency of ROI i towards activation or silence and the parameter 𝐽 𝑖𝑗 represents the effective interaction between ROIs i and j. The estimation of the model parameters Ω = {𝒉, 𝑱} was achieved through a pseudo-likelihood maximization (32) (see Methods). Since this model requires the precise estimation of < 𝜎 𝑖 𝜎 𝑗 >, it cannot be fitted to single-scan data and, for this reason, we used concatenated data from each experimental condition. The pairwise-MEM accurately predicted the observed correlations and, to a lower extend, it predicted the distribution of population activity 𝑃(𝐾) (average correlation fit: 𝑟 = 0.985 ± 0.002; average 𝐷 𝐽𝑆 = 0.006 ± 0.002; Fig. 4G-I and Fig. S2)this is expected, since 𝑃(𝐾)

was not used to constrain the model. We found that biases and couplings parameters were changed for different states, with some parameters increasing or decreasing, and with a reduction of the variance of couplings in the anesthetized states (Fig. S5A-D). Moreover, coupling parameters showed a higher correlation with the anatomical connectivity (or brain connectome) in the anesthetized states than in the awake state (Fig. S5E).

10

Using this model, we calculated the collective statistical quantities for the different experimental conditions. Since in the pairwise-MEM the partition function is not tractable, we calculated 𝐹, 𝜒 and 𝐶 ℎ using Monte Carlo simulations (see Methods). Consistent with the above results, we found that the awake system had larger available energy (free energy, Fig. 4J, see also Fig. S6), larger population fluctuations (susceptibility, Fig. 4K) and larger repertoire of states (heat capacity, Fig.

4L) than the system under anesthesia. Thus, the different versions of the MEM used here indicate the same results concerning the statistical properties of awake and anesthetized states. Furthermore, as shown in the Appendix, the susceptibility can be viewed as a measure of the network response to a small stimulus. Consequently, we found that application of an external stimulus elicited larger and more diverse responses for the pairwise-MEM corresponding to the awake state than for the models corresponding to the anesthetized states (Fig. S7).

Awake collective activity displayed critical dynamics that were shifted to a super-critical regime under anesthesia

The pairwise-MEM can be used to assess the physical state of the system. Indeed, by introducing a scaling parameter 𝑇, analogous to the temperature in statistical physics, one can obtain relevant features of the collective dynamics. For this, we scaled all model parameters as 𝛀 → 𝛀/𝑇 and calculated the heat capacity as a function of 𝑇, given by 𝐶 ℎ (𝑇) = var[𝐸]/𝑇 2 . The "temperature" 𝑇 controls the level of disorder and its effects can be understood by examining the system's energy levels (Fig. S8). Briefly, at low temperatures, interactions dominate over fluctuations making the system predominantly silent and ordered. In contrast, at high temperatures, the system is disordered and relatively uncoupled because fluctuations dominate over interactions. Both low and high temperatures lead to a low 𝐶 ℎ . However, for a specific temperature 𝑇 max , order and disorder coexist in the system and 𝐶 ℎ is maximal as expected for critical dynamics [START_REF] Tkačik | Searching for collective behavior in a large network of sensory neurons[END_REF][START_REF] Hahn | Spontaneous cortical activity is transiently poised close to criticality[END_REF]. Thus, a maximal heat capacity at 𝑇 max = 1 (corresponding to the model learned from the data) suggests that the system operated close to a critical state (whereas 𝑇 max < 1 and 𝑇 max > 1 indicates super-critical and subcritical dynamics, respectively).

We found that the heat capacity curve was maximal for a temperature equal to 1 for the awake state, while it peaked at 𝑇 max < 1 for the anesthetized conditions (Fig. 4M). These results suggest that the awake state displayed critical dynamics, while dynamics under anesthesia were super-critical, which indicates that the anesthetics had a disconnection effect.

Couplings to population relate to the sensitive parameters of the system

We next evaluated how the different parameters affected the model's collective behavior. In general, changes in parameters can differently affect the system's behavior, with some parameters (called "stiff" parameters) effectively modifying it, while others have little effect on it ("sloppy" parameters) [START_REF] Machta | Parameter space compression underlies emergent theories and predictive models[END_REF]. We used an information-theoretical approach based on the Fisher Information Matrix (FIM, noted 𝐺) to detect the parameters that have a strong effect on the collective activity (see Methods). To evaluate the importance of each of the parameters, we defined the parameter's sensitivity as its absolute contribution to the first eigenvector of the FIM (see Methods). The regions with the largest associated sensitivity for parameter 𝛾 were located in the cingulate, parietal, and insular cortices (Fig. 5B). Those that contributed the least were visual and prefrontal cortices. Interestingly, the regions presenting larger reductions of 𝛾 between awake and anesthesia tended to be those with higher associated sensitivity (corr: 0.74, p < 0.001; Fig. 5C).

Finally, we further examined how changes in pairwise correlations between awake and anesthesia related to changes in parameters of different sensitivity. We analyzed the average difference of correlation (∆𝐶) between awake and anesthesia. Two groups of ROIs were clearly separated according to their positive or negative contribution to the first eigenvector of the matrix ∆𝐶, respectively (Fig. 5D). Those that contributed positively were prefrontal and visual cortices, and those that contributed negatively were the cingulate, parietal, and insular cortices (Fig. 5E). Both groups presented a reduction of correlations under anesthesia, but prefrontal and visual cortices were related to parameters of low sensitivity (Fig. 5F-G). Hence, although prefrontal and visual areas changed their correlations, these changes were related to parameters that had a low impact on collective dynamics.

Discussion

In this study, we analyzed the fMRI binary collective activity of monkeys during wakefulness and under anesthesia. We showed that the coupling between each brain region and the rest of the population provides an efficient statistic that discriminates between awake and anesthetized states.

We built MEMs based on this and other statistics to derive macroscopic properties that described the different brain states, such as the free energy 𝐹, the susceptibility 𝜒, and the heat capacity 𝐶 ℎ . All these quantities were maximized in the awake state. By studying the heat capacity curve 𝐶 ℎ (𝑇) as a function of scaling parameter, controlling the disorder of the system, we showed that awake critical dynamics were shifted to supercritical ones under anesthesia. Finally, using the FIM, we showed that changes in brain state were primarily dependent on changes in the couplings to population which were associated with the sensitive parameters of the MEMs and with specific brain regions.

Population couplings and network sensitivity

Previous research at the microcircuit level showed that neurons differ in their coupling to the population activity, with neurons that activate most often when many others are active and neurons that tend to activate more frequently when others are silent [START_REF] Okun | Diverse coupling of neurons to populations in sensory cortex[END_REF]. Using the FIM analysis to detect sloppy and stiff parameters, it has been shown that these different types of neurons have a different impact on the network activity, different stimulus response properties, and different involvement in cortical state transitions [START_REF] Ponce-Alvarez | Cortical state transitions and stimulus response evolve along stiff and sloppy parameter dimensions, respectively[END_REF]. We showed that brain regions coupled differently to the rest of the whole-brain network, consistent with previous findings [START_REF] Tanabe | Altered global brain signal during physiologic, pharmacologic, and pathologic states of unconsciousness in humans and rats[END_REF], and, furthermore, that these couplings primarily determined the collective activity (i.e., they were associated to the stiff parameters of the model) and varied across arousal brain states. Overall, we proposed that the distribution of sensitivity of brain regions and their functional role could be a general principle of neural networks at different scales.

Using principal components, we detected the combination of ROIs that contributed the most to distinguish between the awake and the anesthetized states based on their population couplings. The brain regions that changed their population coupling from the awake state to the anesthetized state were the cingulate, parietal, and insular cortices (Fig. 3E). Notably, the model parameters associated to the couplings of these regions were among those impacting the most the collective dynamics (Fig.

5B-C).

Our results suggest that anesthesia modified some important local/global parameters that effectively induced a change of brain state. Our results highlight the key role of the parieto-cingulate cortex in the mechanism of anesthesia-induced loss of consciousness. Previous studies have shown that the parietal cortex (4, 9, 36) and the cingulate cortex [START_REF] Luppi | Consciousness-specific dynamic interactions of brain integration and functional diversity[END_REF] are most strongly affected by anesthetics. These cortices also present alterations in brain injury-induced unconsciousness in humans [START_REF] Luppi | Consciousness-specific dynamic interactions of brain integration and functional diversity[END_REF][START_REF] Juengling | Separating functional and structural damage in persistent vegetative state using combined voxel-based analysis of 3-D MRI and FDG-PET[END_REF]. Moreover, consistent with our results, it has been shown that the insula plays an important role in awareness and is a potential neural correlate of consciousness [START_REF] Craig | How do you feelnow? The anterior insula and human awareness[END_REF].

Interestingly, although some brain regions, such as visual and prefrontal cortices, had different correlations between awake and anesthesia, they were associated to parameters with low impact on collective activity. This highlights the importance of studying not only the change in statistics between brain states but also their sensitivity on network dynamics. Consistent with our findings, a recent study of neuronal activity from several brain regions and in different arousal states [START_REF] Afrasiabi | More than just front or back: Parietal-striatal-thalamic circuits predict consciousness level[END_REF] shows that parieto-basal ganglia circuits predicted the state of consciousness, while prefrontal activity failed.

In addition, it has been proposed that the prefrontal cortex is mostly involved in the report of consciousness, rather than in the conscious experiences per se [START_REF] Storm | Consciousness regained: disentangling mechanisms, brain systems, and behavioral responses[END_REF].

Macroscopic thermodynamic quantities

Using the MEMs, we learned interesting collective properties describing the different brain states.

We measured the susceptibility that quantifies the diversity of spontaneous population fluctuations.

The susceptibility can be viewed as a measure of the network response to a vanishing stimulus (Fig

S7

, see also the Appendix). Thus, the higher susceptibility observed in the awake state, compared to the anesthetized states, is consistent with Transcranial Magnetic Stimulation (TMS) studies showing that stimulation elicits a more diverse and complex response in the awake state than in low-level states of consciousness, such as sleep, anesthesia, and coma [START_REF] Massimini | Cortical mechanisms of loss of consciousness: insight from TMS/EEG studies[END_REF][START_REF] Sarasso | Quantifying cortical EEG responses to TMS in (un)consciousness[END_REF][START_REF] Sarasso | Consciousness and complexity during unresponsiveness induced by propofol, xenon, and ketamine[END_REF][START_REF] Napolitani | Transcranial magnetic stimulation combined with high-density EEG in altered states of consciousness[END_REF]. Our study predicts that the network response to a localized stimulation would covary with the population couplings and the associated parameter sensitivities of brain region.

Moreover, the models also allowed the estimation of the system's heat capacity, a measure that quantifies the extent of the accessible dynamical repertoire. Indeed, a maximal heat capacity not only indicates that the system can display a large number of energy states, but also that these states are distinguishable (Fig S8). Thus, a large heat capacity indicates a large capacity to represent information in numerous separable states. The observed reduction of heat capacity in the anesthetized states is consistent with previous studies showing that the repertoire of correlation states is limited during anesthesia [START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. Furthermore, by varying a scaling parameter analogous to temperature, the resulting heat capacity curves suggest that awake dynamics were critical, while anesthetized dynamics were supercritical, consistent with previous predictions [START_REF] Tagliazucchi | Large-scale signatures of unconsciousness are consistent with a departure from critical dynamics[END_REF][START_REF] Fekete | Critical dynamics, anesthesia and information integration: Lessons from multi-scale criticality analysis of voltage imaging data[END_REF]. The model used here gives an intuitive interpretation of the transition between critical to supercritical dynamics. Indeed, in the pairwise-MEM supercritical dynamics are associated with a regime in which random fluctuations dominate over interactions, which is consistent with a disconnection of effective couplings. It is important to note that the scale parameter 𝑇 is only introduced to assess the state, i.e., subcritical, critical, or supercritical, of the observed system (the one given for 𝑇 = 1, for which the pairwise-MEM fits the data). This does not mean that differences between awake and anesthetized states are due to a global reduction of interactions and biases, instead different arousal states yielded different biases and couplings (Fig S5) which, in combination, resulted in a change of the system's state. This means that the anesthesia reconfigured the system and not only scaled its parameters. Consistently, we found that effective couplings correlated more with the anatomical connections for the anesthetized states than for the awake state, an effect that has been observed in empirical data [START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF][START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF] and cannot be explained by changes in global connectivity alone [START_REF] Lee | Relationship of critical dynamics, functional connectivity, and states of consciousness in large-scale human brain networks[END_REF].

Lastly, we measured the Helmholtz free energy of the estimated models. The free energy measures the useful energy that can be extracted from the system to the environment, i.e., its ability to produce work. Reasonably, the awake state led to higher free energy than the anesthetized states. Another important property of free energy is that its change with respect to the model parameters is equal to the Fisher information and, thus, it relates to the sensitivity of collective dynamics on these parameters. This result provides a direct link between the sensitivity of parameters and the change of a macroscopic quantity, the free energy, the behavior of which is known to characterize the phase transition [START_REF] Ginzburg | On the theory of superconductors[END_REF]. For the linear coupling-MEM, we showed that the couplings to population (𝑧) were associated to the parameters that have the strongest impact on collective activity. Consistently, we found that that 𝑧 was an efficient observable to classify the arousal states which collective dynamics were qualitatively different (in terms of criticality and supercriticality). Thus, these results give a coherent theoretical justification of the relevance of the statistic 𝑧 to characterize the brain states and to estimate their free energy. Altogether, our findings represent a significant step in the understanding of brain states, resulting in a coherent explanation of the transition from awake to anesthesia: the phase transition between brain states is driven by those parameters that change the free energy, which are the "stiff" parameters of the systems and which relate to population couplings.

Implications for studies on pathological low-level states of consciousness

An interesting extension of this work could be to study brain dynamics in coma using the present statistical mechanics framework. Loss of consciousness due to anesthesia or coma share common features: complexity of dynamics and neural communication are generally reduced in low-level states of consciousness [START_REF] Schartner | Complexity of Multi-Dimensional Spontaneous EEG Decreases during Propofol Induced General Anaesthesia[END_REF][START_REF] Sitt | Large scale screening of neural signatures of consciousness in patients in a vegetative or minimally conscious state[END_REF]. Consequently, estimates of complexity of human brain activity have been used to assess the depth of anesthesia [START_REF] Zhang | EEG complexity as a measure of depth of anesthesia for patients[END_REF][START_REF] Singh | Entropy as an indicator to measure depth of anaesthesia for laryngeal mask airway (LMA) insertion during sevoflurane and propofol anaesthesia[END_REF] and to predict the recovery of consciousness in vegetative patients [START_REF] Sarà | Functional isolation within the cerebral cortex in the vegetative state: a nonlinear method to predict clinical outcomes[END_REF]. Reduction of complexity is consistent with a deviation from critical dynamics when consciousness is lost. Since the coupling-MEMs can be fitted to data from single scans and, as shown here, their parameters change in different brain states, future investigation could use these models and the statistic 𝑧 in the case of disorders of consciousness.

Methods

Animals

This study included a total of five rhesus macaques (Macaca mulatta; 4 females, 1 male, 5-8 kg;

Experimental procedures

Monkeys received anesthesia either with propofol, ketamine, or sevoflurane [START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. The details of the anesthesia protocols are described in the Appendix. Monkeys were scanned on a 3-T horizontal scanner (Siemens Tim Trio; TR, 2,400 ms; TE, 20 ms; and 1.5-mm 3 voxel size; 500 brain volumes per scan session). Before each scanning session, a contrast agent monocrystalline iron oxide nanoparticle (MION) was injected into the monkey's saphenous vein. Acquisition and preprocessing of functional images followed the standard steps described in [START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF] and in the Appendix. Time-series were obtained for 𝑁 = 82 previously defined cortical regions of interest (ROIs) (CoCoMac Regional Map parcellation).

Data binarization

fMRI time series were binarized to study the data statistics and to learn two different families of maximum entropy models (MEMs). While binarization was required to construct the MEMs, transformation of continuous fMRI signals into discrete point processes has proven to effectively capture and compress fMRI large-scale dynamics [START_REF] Tagliazucchi | Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis[END_REF]. Indeed, it has been shown that point process resulting from signal thresholding largely overlaps with deconvoluted fMRI signals using the hemodynamic response function and preserve the topology of the resting state networks (RSNs) [START_REF] Solovey | Loss of consciousness is associated with stabilization of cortical activity[END_REF].

We here discretized the signals as follows. For each scan, the z-scored time-series of each ROI, 𝑥 𝑖 (𝑡)

(1 ≤ 𝑖 ≤ 𝑁), was binarized by imposing a threshold 𝜃 = -1. Two binarization procedures were used.

The first method detects the threshold crossings: the binarized activity is 𝜎 𝑖 (𝑡) = 1 if 𝑥 𝑖 (𝑡) < 𝜃 and 𝑥 𝑖 (𝑡 -1) > 𝜃, and 𝜎 𝑖 (𝑡) = 0 otherwise. The second method assigns the values 1 and -1 to all time points below or above the threshold, respectively: 𝜎 𝑖 (𝑡) = 1 if 𝑥 𝑖 (𝑡) < 𝜃, and 𝜎 𝑖 (𝑡) = -1 otherwise.

The first and the second procedure result in sparse and dense binary activity, respectively. We used the sparse and dense methods to construct coupling-MEMs and pairwise-MEM, respectively. This was to meet the assumptions of the model inference (see Appendix).

k-means classification

We used k-means clustering to classify the scans based on different statistics. Let 𝒗 (𝑖) be a vector calculated from scan 𝑖, e.g., the vector containing all pairwise correlations among the ROIs. We used k-means to partition the collection of 𝒗 (𝑖) into a pre-specified number (𝑘) of clusters. k-means minimizes the within-cluster variation, over all clusters. We used 𝑘 = 2 to evaluate how well the scans corresponding to the awake state and those corresponding to the anesthetized states (independent of the anesthetic protocol) could be classified based on vectors 𝒗 (𝑖) . To classify the six different experimental conditions, we used 𝑘 = 6. The classification performance was given by the proportion of correctly clustered scans. We used 100 random initial conditions of the k-means algorithm to obtain the average classification performance and its uncertainty.

Maximum entropy models (MEMs)

MEMs estimate the probability of all possible binary patterns, 𝑃(𝝈), that matches the expectation of a set of data observables. Let 𝒪 1 (𝝈), … , 𝒪 𝐿 (𝝈) be the set of 𝐿 data observables we seek to preserve.

For example, if we were interested only on activation rates, 〈𝜎 𝑖 〉, we would need to consider 𝑁 observables 𝜎 1 , … , 𝜎 𝑁 . Under the model distribution 𝑃(𝝈), the observables' expectations are given as:

〈𝒪 𝑖 〉 model = ∑ 𝑃(𝝈)𝒪 𝑖 (𝝈) {𝜎} , ( 5
)
and should fit those of the data,

〈𝒪 𝑖 〉 data = 1 𝑛 ∑ 𝒪 𝑖 (𝝈 𝑡 ) 𝑛 𝑡=1
, where 𝝈 𝑡 is the observed pattern at time 𝑡 (1 ≤ 𝑡 ≤ 𝑛). We search for the model distribution 𝑃(𝝈) that does less assumptions, i.e., the one that has maximal entropy 𝑆 = ∑ 𝑃(𝝈)ln𝑃(𝝈) {𝝈} . Thus, the problem is equivalent to maximizing a function (the entropy) given some constraints on the expectation values of the observables, a problem that can be generally solved using Lagrange multipliers. The maximum entropy distribution has the general form:

𝑃(𝝈) = 1 𝑍 exp (∑ Ω 𝑖 𝐿 𝑖=1 𝒪 𝑖 (𝝈)) = 1 𝑍 exp(-𝐸(𝝈)), (6)
Where Non-linear coupling-MEM. The above model can be extended to include the non-linear coupling between 𝜎 𝑖 and 𝐾. The complete coupling between 𝜎 𝑖 and 𝐾 is provided by the joint probability distributions of 𝜎 𝑖 and 𝐾, i.e., 𝑃(𝜎 𝑖 , 𝐾), which is the target observable of the non-linear coupling-MEM. In this case, the energy is given as Pairwise-MEM. The third model we considered is the one that targets the activation rates (< 𝜎 𝑖 >)

𝐸(𝝈) = -∑ 𝐻 𝑖,𝐾(𝝈) 𝜎 𝑖 𝑁 𝑖=1 (29)
and the pairwise correlations (< 𝜎 𝑖 𝜎 𝑗 >) of the data. The resulting energy function of the maximum entropy distribution is given as [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF]. The model parameter ℎ 𝑖 , called intrinsic bias, represents the intrinsic tendency of neuron i towards activation (𝜎 𝑖 = +1) or silence (𝜎 𝑖 = -1) and the parameter 𝐽 𝑖𝑗 represents the effective interaction between neurons i and j.

𝐸(𝝈) = -∑ ℎ 𝑖 𝑁 𝑖=1 𝜎 𝑖 - 1 2 ∑ ∑ 𝐽 𝑖𝑗 𝜎 𝑖 𝜎 𝑗 𝑁 𝑗=1 𝑁 𝑖=1
The parameters of the coupling-MEMs and pairwise-MEMs were estimated from the data using likelihood [START_REF] Gardella | A Tractable Method for Describing Complex Couplings between Neurons and Population Rate[END_REF] and pseudo-likelihood (32) maximization, respectively (see Appendix).

Macroscopic quantities

The analysis of the learned MEMs provides relevant properties of the collective activity. These quantities derive from the Boltzmann distribution and they are interpretable in the framework of statistical physics. The description and calculation of these quantities are presented in the Appendix in detail. Briefly, we studied the system's Helmholtz free energy, susceptibility, and heat capacity.

The free energy 𝐹 is given by the difference between the average energy and the entropy, i.e. 𝐹 = 〈𝐸〉 -𝑆 = -ln(Z); it quantifies the useful energy that is obtainable from the system. The susceptibility 𝜒 relates to the diversity of population states, i.e. 𝜒 = var(𝐾), but, importantly, it also relates to the system's response to intrinsic or external inputs (see Appendix and Fig. S7). The heat capacity 𝐶 ℎ quantifies the diversity of accessible energy states, i.e. 𝐶 ℎ = var(𝐸). The heat capacity measures the size of the dynamic repertoire of the system. Furthermore, a parameter 𝑇, that scales all model parameters (𝛀 → 𝛀/𝑇), can be introduced to study the effect of a change in the system's disorder ("temperature") on the repertoire of accessible energy states, i.e., the function 𝐶 ℎ (𝑇) = var(𝐸)/𝑇 2 . This function is informative of the state of the system in terms of criticality: a maximum of the heat capacity close to 𝑇 max = 1 suggests that the observed system is likely to be close to a critical state, whereas 𝑇 max < 1 and 𝑇 max > 1 indicate super-critical and sub-critical dynamics, respectively [START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF][START_REF] Tkačik | Searching for collective behavior in a large network of sensory neurons[END_REF][START_REF] Hahn | Spontaneous cortical activity is transiently poised close to criticality[END_REF] (see Appendix and Fig. S8).

Fisher information matrix

We were interested in detecting which parameters have the strongest effect on the collective activity.

For this, we studied the Fisher information matrix (FIM, noted 𝑮) of the learned MEMs. The FIM represents the curvature of the log-likelihood of the model, log 𝑃(𝝈|𝛀), with respect to the model parameters, i.e., it quantifies the sensitivity of the model to changes in parameters. It is given as:

𝐺 𝑘𝑙 = -∑ 𝑃(𝝈|𝛀) 𝜕 2 ln 𝑃(𝝈|𝛀) 𝜕Ω 𝑘 𝜕Ω 𝑙 {𝝈} . ( 7
)
where 1 ≤ 𝑘, 𝑙 ≤ 𝐿, where 𝐿 is the number of parameters. As shown in Appendix, the FIM is given by the second derivatives of the free energy:

𝐺 𝑘𝑙 = 𝜕 2 ln 𝑍 𝜕Ω 𝑘 𝜕Ω 𝑙 = - 𝜕 2 𝐹 𝜕Ω 𝑘 𝜕Ω 𝑙 . ( 8
)
To quantify the sensitivity of the different parameters we decomposed the FIM into eigenvectors, noted 𝝂 1 , … , 𝝂 𝐿 , and measured the sensitivity of a given parameter 𝑖 by its absolute contribution to the first eigenvector, i.e. |𝜈 1 (𝑖)|.

Statistical analysis

We used one-way ANOVA followed by Tukey's post hoc analysis to compare the means of three or more distributions and Wilcoxon rank sum test to compare the medians of two distributions. We measured the dissimilarity between two distributions (i.e., data vs. model distribution) through the Jensen-Shannon divergence. Correlation matrices were analyzed using standard PCA. Statistical models (i.e., maximum entropy models) were estimated using likelihood and pseudo-likelihood maximization.

MNI space [START_REF] Frey | An MRI based average macaque monkey stereotaxic atlas and space (MNI monkey space)[END_REF], and band-pass filtered in the frequency range of interest (0.0025-0.05 Hz). Timeseries were obtained for 𝑁 = 82 previously defined cortical regions of interest (ROIs) (CoCoMac Regional Map parcellation). Scans that presented signs of artifacts in time-series or power spectral density were discarded. The procedure was based on the visual inspection of the time series for all the nodes, the Fourier transform of each signal. A total of 119 scans were kept for subsequent analyses, corresponding to different levels of arousal: wakefulness (n = 24 scans), two levels of propofol sedation (light, LPP, n = 21, and deep, DPP, n = 23), ketamine anesthesia (KETA, n = 22), and two types of sevoflurane anesthesia (SEV2, n = 18, and SEV4, n = 11).

Anatomical connectivity

We used a fully weighted whole-cortex macaque structural connectivity matrix (connectome) derived by combining the information from fiber-tracing and tractography [START_REF] Shen | A macaque connectome for large-scale network simulations in TheVirtualBrain[END_REF]. The connectome is publicly available here: https://zenodo.org/record/1471588#.X44C6dAzY2x. Briefly, the tractography algorithm was optimized to best reproduce the weighted but partial-cortex tracer connectome from Markov et al. [START_REF] Markov | A weighted and directed interareal connectivity matrix for macaque cerebral cortex[END_REF], before estimating whole-cortex connectome weights. The directed connectome weights between ROIs of the CoCoMac parcellation were given as the number of streamlines detected between them, divided by the total number of streamlines that were sent from the seed. Tractographyderived structural connectivity matrices were averaged across nine macaque monkeys. For details see Shen et al. [START_REF] Shen | A macaque connectome for large-scale network simulations in TheVirtualBrain[END_REF].

Macroscopic quantities

The analysis of the learned MEMs provides relevant properties of the collective activity. These quantities derive from the Boltzmann distribution and they are interpretable in the framework of statistical physics. Using the distribution 𝑃(𝝈) the mean energy and the entropy are given as:

〈𝐸〉 = ∑ 𝑃(𝝈)𝐸(𝝈) {𝝈} , (S1) 𝑆 = -∑ 𝑃(𝝈)ln𝑃(𝝈) {𝝈} . ( S2
)
An important quantity describing the system is the Helmholtz free energy defined as: 𝐹 = 〈𝐸〉 -𝑆.

Using the above expressions for the mean energy and the entropy we get: Thus, the free energy can be obtained either by calculating the mean energy and the entropy or directly by means of the partition function. These two strategies can be more or less difficult depending on the model. For the coupling-MEMs, we measured 𝐹 for each scan through 𝑍, since 𝑍 is tractable and the models can be fitted to single-scan data. For the pairwise-MEM, for which 𝑍 is not tractable, we performed Metropolis Monte Carlo simulations of the model (10 6 steps) to estimate 〈𝐸〉 and 𝑆. The Monte Carlo simulations were repeated ten times to estimate uncertainties on these quantities. The free energy quantifies the useful energy that can be extracted from the system, called 'work' in physics [START_REF] Chatelain | Probability distributions of the work in the 2D-Ising model[END_REF].

Another important quantity is the susceptibility 𝜒 that measures the fluctuations of the population activity, i.e. 𝜒 = var(𝐾). The susceptibility can be also viewed as a measure of the population response to a perturbation applied to all units. Assuming that the perturbation 𝐵 adds a term to the energy, i.e., 𝐸(𝝈) + ∑ 𝐵𝜎 𝑖 𝑁 𝑖=1

, the susceptibility is given by the derivative of 〈𝐾〉 with respect to 𝐵, The zero-field susceptibility measures the spontaneous fluctuations of the population activity. This quantity can be measured directly in the coupling-MEMs through the estimated distribution 𝑃(𝐾). In the pairwise-MEM, the variance of 𝐾 can be estimated using Metropolis Monte Carlo simulations (ten simulations of 10 6 steps). In addition, note that 〈𝐾〉 = A last important quantity is the heat capacity. The heat capacity 𝐶 ℎ quantifies the diversity of accessible energy states of the system, i.e., 𝐶 ℎ = var(𝐸). A useful related measure is the variation of the heat capacity as a function of a scale parameter 𝑇, analogous to temperature in statistical physics.

This parameter acts as a scaling factor for all model parameters as 𝛀 → 𝛀/𝑇. The "temperature" 𝑇 controls the level of disorder and its effects can be understood by examining the system's energy levels (Fig. S8). This creates a family of scaled models in which 𝑇 = 1 corresponds to the MEM that was fitted to the data. The heat capacity as a function of 𝑇 is given by 𝐶 ℎ (𝑇) = var(𝐸)/𝑇 2 and provides useful features of the system. Indeed, it is known that a maximum of the heat capacity close to 𝑇 = 1 suggests that the observed system is likely to be close to a critical state, whereas 𝑇 max < 1 and 𝑇 max > 1 indicate super-critical and sub-critical dynamics, respectively [START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF][START_REF] Tkačik | Searching for collective behavior in a large network of sensory neurons[END_REF][START_REF] Hahn | Spontaneous cortical activity is transiently poised close to criticality[END_REF]. Hence, the curve 𝐶 ℎ (𝑇) can be used as a tool to assess criticality. The heat capacity measures the size of the dynamic repertoire of the system. It not only provides a measure of the system's complexity, but also assess whether the complexity is maximized and whether any reduction of complexity is due to a transition to subcritical or supercritical regimes (both regimes result in a decrease of complexity with respect to criticality, but through different mechanisms, see Fig. S8). For instance, if 𝑇 max ≠ 1, this means that the system can be re-scaled to increase the complexity of the model dynamics (e.g., if

𝑇 max < 1, 𝑇 needs to decrease to reach the maximum heat capacity, indicating supercritical dynamics).

Furthermore it can be shown that the entropy 𝑆 of the system is equal the integral of the function 𝐶 ℎ (𝑇)/𝑇 from 𝑇 = 0 to 𝑇 = 1 [START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF][START_REF] Tkačik | Searching for collective behavior in a large network of sensory neurons[END_REF]. This is a useful strategy to calculate the entropy that we used to compute the free energy when direct access to the partition function was not feasible (i.e., for the pairwise-MEM). In our study, we calculated 𝐶 ℎ (𝑇) by estimating the variance of the energies through Monte Carlo simulations of the pairwise-MEM for different 𝑇 (five simulations of 5.10 6 steps for each 𝑇).

Finally, note that the energies are equal to the patterns' minus log probabilities, or "surprise", minus the free energy (a constant), i.e., 𝐸(𝝈) = -ln𝑃(𝝈) + ln𝑍. Thus, the variance of the energy (heat capacity) measures the range of surprises of the different collective states. A large heat capacity allows the system to represent sensory events that occur with a wide range of likelihoods (energy states that are distributed, numerous, and separable) [START_REF] Tkacik | Ising models for networks of real neurons[END_REF].

Fisher information matrix and free energy

We were interested in detecting which parameters have the strongest effect on the collective activity.

To measure how distinguishable two models, with parameters Ω and Ω + 𝛿Ω, are based on their predictions, we used the Fisher information matrix (FIM). Indeed, the Kullback-Leibler divergence between the two models can be written as:

𝐷 𝐾𝐿 (𝛀; 𝛀 + 𝛿𝛀) = 𝐺 𝑘𝑙 𝛿Ω 𝑘 𝛿Ω 𝑙 + 𝒪(𝛿Ω 3 ), (S5)

where 1 ≤ 𝑘, 𝑙 ≤ 𝐿, where 𝐿 is the number of parameters, and the matrix 𝑮 is the FIM matrix given by:

𝐺 𝑘𝑙 = -∑ 𝑃(𝝈|𝛀) 𝜕 2 ln 𝑃(𝝈|𝛀) 𝜕Ω 𝑘 𝜕Ω 𝑙 {𝝈} . ( S6
)
The FIM represents the curvature of the log-likelihood of the model, log 𝑃(𝝈|𝛀), with respect to the model parameters, i.e., it quantifies the sensitivity of the model to changes in parameters.

Note that the FIM relates to the free energy 𝐹 = -ln(Z). Indeed, using the Boltzmann distribution, we have:

𝜕 2 ln𝑃 𝜕Ω 𝑘 𝜕Ω 𝑙 = 𝜕 𝜕Ω 𝑘 𝜕Ω 𝑙 [-𝐸(𝝈) -ln (𝑍)]. (S7)
In MEMs the energy is given by the parameters (i.e., Lagrange multipliers) and the target observables as: 𝐸(𝝈) = ∑ Ω 𝑖 𝑖 𝒪 𝑖 (𝝈). Thus, all second derivatives of the energy with respect to the parameters are zero, 𝜕 2 𝐸 𝜕Ω 𝑘 𝜕Ω 𝑙 = 0, and we have:

𝜕 2 ln 𝑃 𝜕Ω 𝑘 𝜕Ω 𝑙 = - 𝜕 2 ln 𝑍 𝜕Ω 𝑘 𝜕Ω 𝑙 . ( S8
)
Since the right-hand term of equation 17 does not depend on 𝝈, equation 15 gives:

𝐺 𝑘𝑙 = 𝜕 2 ln 𝑍 𝜕Ω 𝑘 𝜕Ω 𝑙 = - 𝜕 2 𝐹 𝜕Ω 𝑘 𝜕Ω 𝑙 . ( S9
)
Estimation of maximum entropy models.

Linear and non-linear coupling-MEMs were estimated using the method described in [START_REF] Gardella | A Tractable Method for Describing Complex Couplings between Neurons and Population Rate[END_REF]. Briefly, Newton's method was used to maximize the log-likelihood 𝐿𝐿, by iteratively updating the parameters as:

𝛀 (𝑘+1) = 𝛀 (𝑘) -𝑎𝑮 -𝟏 𝛁𝑳𝑳, (S10)
Where 𝛀 (𝑘) is the vector of parameters at iteration 𝑘, 𝑎 is a learning rate, 𝑮 is the FIM with parameter 𝛀 (𝑘) , and 𝛁𝑳𝑳 is the gradient of 𝐿𝐿 with respect to the model parameters 𝛀 . This method is feasible because, for coupling-MEM, the partition function can be analytically obtained from the model parameters and thus providing 𝑮 and 𝛁𝑳𝑳; see details in [START_REF] Gardella | A Tractable Method for Describing Complex Couplings between Neurons and Population Rate[END_REF]. The algorithm was stopped when the estimation errors for all constraint observables became lower than 10 -6 . We found that this algorithm correctly converged for the sparse binarization. The code for learning the coupling-MEMs is available at:

https://github.com/ChrisGll/MaxEnt_Model_Population_Coupling.

The pairwise-MEM was learned using pseudo-likelihood maximization [START_REF] Ezaki | Energy landscape analysis of neuroimaging data[END_REF]. This method approximates the likelihood function:

𝐿(𝛀) ≈ ∏ ∏ 𝑃 ̃(𝜎 𝑖 |𝛀, 𝜎 ≠𝑖 ), 𝑁 𝑖=1 𝑛 𝑡=1 (S11)
where 𝑛 is the number of time points of the data and 𝑃 ̃(𝜎 𝑖 |𝛀, 𝜎 ≠𝑖 ) is the conditional Boltzmann distribution for a unit 𝜎 𝑖 .

This was done by updating the biases and couplings as:

ℎ 𝑖 (𝑘+1) = ℎ 𝑖 (𝑘) + 𝑎[〈𝜎 𝑖 〉 data -〈𝜎 𝑖 〉 mf ], (S12) 𝐽 𝑖𝑗 (𝑘+1) = 𝐽 𝑖𝑗 (𝑘) + 𝑎[〈𝜎 𝑖 𝜎 𝑗 〉 data -〈𝜎 𝑖 𝜎 𝑗 〉 mf ], (S13)
Where 𝑘 denotes the updating iteration (up to 5.10 3 ) and 〈. 〉 mf are the expected values using the mean-field approximation:

〈𝜎 𝑖 〉 mf = 1 𝑛 ∑ tanh (ℎ 𝑖 + ∑ 𝐽 𝑖𝑗 𝜎 𝑗 (𝑡) 𝑗≠𝑖 ) , 𝑛 𝑡=1 (S14) 〈𝜎 𝑖 𝜎 𝑗 〉 mf = 1 𝑛 ∑ 𝜎 𝑗 (𝑡)tanh (ℎ 𝑖 + ∑ 𝐽 𝑖𝑘 𝜎 𝑘 (𝑡) 𝑘≠𝑖 ) , 𝑛 𝑡=1 (S15)
where 𝜎 𝑖 (𝑡) is the activity of ROI 𝑖 (taking values 1 or -1). The estimator of the maximum pseudolikelihood approaches the maximum-likelihood estimator for 𝑛 → ∞ (56). For this reason, pairwise-MEMs were fitted to concatenated data for each of the six experimental conditions. Also, since this method uses a mean-field approximation, it is accurate when the nodes of the network receive many inputs, thus we used the dense binarization scheme (see Methods). The code for learning the pairwise-MEMs is available at: https://royalsocietypublishing.org/doi/suppl/10.1098/rsta.2016.0287. If the model learned from the data (i.e., 𝑇 = 1) maximized the heat capacity, then its dynamics are critical. If the temperature that maximizes the heat capacity is larger than 1, the dynamics of the learn model are subcritical (the temperature must increase to reach the maximum). Finally, if the temperature that maximizes the heat capacity is lower than 1, the dynamics of the learn model are supercritical (the temperature must decrease to reach the maximum). B) Observed energies 𝐸 for subcritical, critical, and supercritical dynamics. 𝐸 0 correspond to the energy of the state for which all units are silent. For subcritical dynamics, the few visited energy levels are sparsely distributed, and the variance of energies (i.e., heat capacity) is relatively low. For supercritical dynamics, many energy levels are densely distributed, and their variance of energy is also relatively low. For critical dynamics, the energy levels are numerous and separable, leading to a maximal variance of energies.

Figure 1 .

 1 Figure 1. Binarization and statistics. (A) MION fMRI signals were z-scored and binarized by imposing a threshold equal to the standard deviation, for each signal. (B) In each time bin of 2.4 s, the state of signal of ROI 𝑖, noted 𝜎 𝑖 (𝑡), was equal to 1 if the MION signal for this ROI was lower than minus its standard deviation, or equal to 0 otherwise. The average activity of ROI 𝑖 was 〈𝜎 𝑖 〉 = ∑ 𝜎 𝑖 (𝑡) 𝑛 𝑡=1

Figure 2 .

 2 Figure 2. Coupling to population predicts the functional connectivity and is consistent within experimental conditions. (A) Top and left insets: ROIs were ordered according to 𝑧 𝑖 averaged for each scan within a given experimental condition (here for the awake condition). Color map: The average functional connectivity (FC) is shown after ordering the ROIs according to 𝒛. (B) Same as (A) but for the deep propofol (DPP) anesthesia condition. (C)The elements of the FC and the corresponding products of coupling to population (𝑧 𝑖 𝑧 𝑗 ) highly correlated. (D) We tested whether 〈𝜎〉, 𝑧, and FC were similar across scans within the same experimental condition. For example, for the statistic 𝑧, we calculate the correlation of this 𝑁-dimensional variable for all pairs of scans belonging to the same experimental condition and computed the distribution of correlation coefficients (blue distribution). High correlation coefficients indicate that, within experimental conditions, scans yielded similar vector 𝒛. The same can be done for the 𝑁-dimensional variable 〈𝜎〉 (green distribution) and the vector of FC elements (𝑁(𝑁 -1)/2 dimensions; red distribution).

Figure 3 .

 3 Figure 3. Coupling to population predicts the state of the brain.(A-C) Correlation matrix comparing the statistics 𝑧, 〈𝜎〉, and FC among all scans. For example, in panel (A), the element (𝑘, 𝑙) of the matrix represents the correlation between the coupling to population vector 𝒛 of scans 𝑘 and 𝑙. Coupling to population clearly separated awake and anesthesia data. Using k-means, we evaluated how well the different statistics could be used to classify the awake and anesthetized conditions (chance level: 50%). The classification performance using the coupling to population statistic was 96.64%, that was significantly higher than using the mean activity (73.95%) or the functional connectivity (74.79%). Classification of the six experimental conditions was generally lower, but higher for 𝑧 than for 〈𝜎〉 and FC (38% vs. 33% and 29%, chance level: 16.67%). (D) PCA analysis showed that 𝑧 vectors separated the awake and anesthetized conditions along the first principal component (PC1). Each dot represents a scan. (E) The absolute coefficient of PC1 associated to each ROI. (F) During anesthesia, 𝑧 was reduced compare to wakefulness for most of the ROIs. Changes from awake baseline, ∆𝑧(𝑖) = 𝑧(Awake) -𝑧(𝑖), where 𝑧 was averaged over scans, were highly correlated for the different anesthetics (with correlation coefficients ranging from 0.85-0.93). The panel shows the comparison between ∆𝑧(DPP) and ∆𝑧(KET).

Figure 4 .

 4 Figure 4. Maximum entropy models indicate higher free energy and heat capacity during wakefulness than during anesthesia. (A) Fitting of 𝑃(𝐾) using the non-linear coupling-MEM. Data and predictions from all scans from the awake condition. (B-C) Fitting of covariances using the linear (B) and nonlinear (C) coupling-MEM for the awake condition. (D-F) The free energy, the susceptibility, and the heat capacity were derived using non-linear coupling-MEMs for the different conditions. Similar results were obtained using the linear model (see Fig. S4). Squares and error bars indicate means and standard deviations across scans, respectively, and the asterisks indicate significantly different values for the awake condition (p < 0.001 one-way ANOVA followed by Tukey's post hoc analysis). (G-L) same as (A-F) but using pairwise-MEMs. Error bars indicate standard errors across Monte Carlo simulations of the models. Asterisks indicate significantly different values for the awake condition (p < 0.001, oneway ANOVA followed by Tukey's post hoc analysis). (M) Heat capacity as a function of temperature. The peak of heat capacity for 𝑇 = 1 indicates critical dynamics during wakefulness. The heat capacity peaked at 𝑇 < 1 for the anesthetized conditions, indicating supercritical dynamics during anesthetized states.

  The FIM measures the change in the model log-likelihood 𝑃(𝝈|𝛀) with respect to changes in the model parameters 𝛀. As demonstrated in the Appendix, the FIM relates to the second derivatives of the free energy with respect to the model parameters, i.e., 𝐺 𝑖𝑗 = -𝜕 2 𝐹 𝜕Ω 𝑖 𝜕Ω 𝑖 . This relation provides a direct link between a macroscopic quantity (the free energy) describing the collective dynamics of the different brain states and the underlying model parameters. For the linear model, the parameters that contributed the most to the FIM were the parameters 𝛾 𝑖 (Fig5A). This explains how changes in couplings to population, as observed between awake and anesthetized states, effectively change the collective state of the system, leading to the observed shift from critical to supercritical dynamics.

Figure 5 .

 5 Figure 5. Model sensitivity to the different parameters. (A) Fisher information matrix (FIM) calculated using the linear coupling-MEM built using the data concatenated across scans from the awake condition. The FIM measures how much the model log-likelihood changes with respect to changes in the model parameters, i.e., 𝛀 = {𝜶, 𝜷, 𝜸} for the linear coupling-MEM. (B) For each scan, we decomposed the FIM in eigenvectors and the mean contribution to first eigenvector of parameters 𝛾 was represented on the brain's image. This represents the sensitivity of collective activity on the parameters 𝛾 associated to the different ROIs. (C) The ROIs with larger reduction of parameter 𝛾 in anesthetized states with respect to the awake state were those with strongest sensitivity. (D) Change of pairwise correlations between awake and anesthetized states: ∆𝐶 = 𝐶 awake -〈𝐶 anesth 〉. ROIs were ordered according to their contribution to the first eigenvector of ∆𝐶 (top and right insets). Two groups of ROIs were detected according to their positive or negative contribution to this eigenvector, respectively, with both groups reducing the correlation between awake and anesthesia. (E) First eigenvector of ∆𝐶 represented in the brain. (F-G) The two groups of ROIs had significantly different associated sensitivity (p < 0.001, Wilcoxon rank sum test), as measured by the FIM values associated to parameters 𝛾.

  𝛀 = [Ω 1 , … , Ω 𝐿 ] are the Lagrange multipliers enforcing the constraints, 𝐸(𝝈) = -∑ Ω 𝑖 𝐿 𝑖=1 𝒪 𝑖 (𝝈) represents the energy of the pattern, and the normalizing factor 𝑍 = ∑ exp(-𝐸(𝝈)) {𝝈} is the partition function (see Appendix). We estimated different MEMs built on different constrained data observables. Linear coupling-MEM. First, we considered the MEM that is consistent with the probability distribution 𝑃(𝐾), the average activations 〈𝜎 𝑖 〉, and the linear coupling between 𝜎 𝑖 and 𝐾, i.e., 〈𝜎 𝑖 𝐾〉 (which relates to 𝑧 𝑖 ). As shown in Gardella et al. (29) the resulting energy function is given as: 𝐸(𝝈) = -∑ (𝛼 𝑖 + 𝛽 𝐾 + 𝛾 𝑖 )𝜎 𝑖 𝑁 𝑖=1 . The model parameters 𝜶 = [𝛼 1 , … , 𝛼 𝑁 ], 𝜷 = [𝛽 0 , … , 𝛽 𝑁 ], and 𝜸 = [𝛾 1 , … , 𝛾 𝑁 ] are Lagrange multipliers associated to the constrained observables 〈𝜎 𝑖 〉, 𝑃(𝐾), and 〈𝜎 𝑖 𝐾〉, respectively.

  , where 𝐾(𝝈) is the number of active ROIs in pattern 𝝈 and the parameters 𝐻 𝑖,𝐾(𝝈) are associated to the constrained observables 〈𝜎 𝑖 𝛿 𝐾,𝑘 〉, where 𝑃(𝐾 = 𝑘) = 〈𝛿 𝐾,𝑘 〉 and 𝛿 𝐾,𝑘 is the Kronecker's delta. The linear model is a special case of this model with 𝐻 𝑖,𝐾(𝝈) = 𝛼 𝑖 + 𝛽 𝐾 + 𝛾 𝑖 .

  𝐹 = ∑ 𝑃(𝝈)[𝐸(𝝈) + ln𝑃(𝝈)]

  i.e., 𝜒 = 𝜕〈𝐾〉 𝜕𝐵 ⁄ (see Fig S7). Here, we were interested in the so-called zero-field susceptibility obtained for 𝐵 = 0. Using the Boltzmann distribution and noting that 〈𝐾〉 = be shown that 𝜕〈𝐾〉 𝜕𝐵 ⁄ is equal to the variance of 𝐾, and thus: 2 〉 -〈𝐾〉 2 = var(𝐾). (S4)

  the work 𝑊 produced by changing the external stimulus from 𝐵 1 to 𝐵 2 , i.e., 𝑊 = -∫ 〈𝐾〉 𝐵 2 𝐵 1 𝑑𝐵, relates to the variation of the free energy ∆𝐹.

Figure S2 .

 S2 Figure S2. Goodness-of-fit and prediction accuracy of MEMs. We calculated the goodness-of-fit and prediction accuracy of the three types MEMs. For the linear and non-linear coupling-MEMs, the goodness-of-fit was given by the fit of the 𝑃(𝐾) distribution, which was the target observable of the MEMs. For these models, the fit of covariances represent predictions of the models, because covariances were not used to construct them. In contrast, for the pairwise-MEM, the goodness-of-fit was given by the covariance fit (since covariances were the target of the model's learning step) and the fit of 𝑃(𝐾) was a model prediction. A) Jensen-Shannon divergence (𝐷 𝐽𝑆 ) between the model and data 𝑃(𝐾) distributions, for the three types of MEMs (lower values of 𝐷 𝐽𝑆 indicate better approximation of 𝑃(𝐾)). ANOVA: linear coupling-MEM: F(5,118) = 1.3, p = 0.289; non-linear coupling-MEM: F(5,118) = 1.2, p = 0.298. For the pairwise-MEM, ANOVA was not applicable since models were learned using concatenated data. B) Covariance fit for the three types of MEMs. ANOVA: linear coupling-MEM: F(5,118) = 11.4, p < 0.001; non-linear coupling-MEM: F(5,118) = 16.6, p < 0.001. Error bars indicate SEM.

Figure S3 .

 S3 Figure S3. Parameters 𝜸 predict the state of the brain. A) PCA analysis showed that the parameters 𝛾 = [𝛾 1 , … , 𝛾 𝑁 ] from the linear coupling-MEM separated the awake and anesthetized conditions along the first principal component. Each point represents a scan. B) Using k-means, we evaluated how well the different parameters 𝛾 and 𝛼 could be used to classify the awake and anesthetized conditions (chance level: 50%). The classification performance based on 𝛾 and 𝛼 parameters was 85.6 ± 4.4%, and 75.8 ± 0.0%, respectively. Classification of the six experimental conditions was 45.0 ± 4.0%, and 27.9 ± 0.9% based on 𝛾 and 𝛼 parameters, respectively.

Figure S4 .

 S4 Figure S4. Macroscopic quantities using the linear coupling-MEM. A) Free energy 𝐹 for the different conditions, 𝐹 = 〈𝐸〉 -𝑆 = -ln 𝑍. B) Susceptibility 𝜒 for the different conditions, 𝜒 = var(𝐾). C) Heat capacity 𝐶 for the different conditions, 𝐶 ℎ = var(𝐸). In panels (A) to (B), we the linear model was used. Squares and error bars indicate means and standard deviations across scans, respectively, and the asterisks indicate significantly different values for the awake condition (p < 0.001 one-way ANOVA followed by Tukey's post hoc analysis).

Figure S5 .

 S5 Figure S5. Change of the pairwise-MEM parameters 𝒉 and 𝑱 across brain states. A) Change in bias model parameters (ℎ) between the awake and deep propofol (DPP) conditions for each ROI. B) Distribution of the change in couplings (𝐽) between the awake and deep propofol (DPP) conditions for all pairs of ROIs. Inset: the matrix represents the change of the couplings between awake and DPP, i.e., Δ𝐽 𝑖𝑗 = 𝐽 𝑖𝑗 awake -𝐽 𝑖𝑗 DPP . C) Distribution of couplings 𝐽 in the awake and DPP conditions. D) Variance of couplings 𝐽 for each experimental condition. Error bars indicate bootstrap uncertainties (500 repetitions). E) Pearson correlation between the coupling matrix 𝐽 and the structural connectivity, for each experimental condition. The error bars indicate the correlation coefficient's 95% confidence interval.

Figure S6 .

 S6 Figure S6. Mean energy and entropy for the pairwise-MEM. The mean energy 〈𝐸〉 and the entropy 𝑆 of the learned pairwise-MEMs were calculated using Monte Carlo simulations (10 6 steps). To calculate the entropy, we used the heat capacity as a function of the scaled parameter 𝑇, analogous to temperature (see Methods).

Figure S7 .

 S7 Figure S7. Population response to an external stimulus predicted by the pairwise-MEMs. A) An external stimulus was applied to the learned pairwise-MEMs. The stimulus 𝐵 added a term to the energy as: 𝐸(𝜎) + ∑ 𝐵𝜎 𝑖 𝑁 𝑖=1

Figure S8 .

 S8 Figure S8. Heat capacity and energy levels. A) Heat capacity 𝐶 ℎ and entropy 𝑆 of the pairwise-MEM, as a function of the scalling parameter 𝑇, analogous to temperature in statistical physics. The temperature 𝑇 controls the disorder of the system (i.e., the entropy 𝑆(𝑇) increases as a function of 𝑇).

  [𝜎 1 , … , 𝜎 𝑁 ]. Estimating the distribution 𝑃(𝝈) over the 2 𝑁 possible binary patterns from the data is impractical with limited amount of observations, since for N = 82 there are more than 10 24 possible patterns. A useful technique to estimate 𝑃(𝝈) relies on the maximum entropy principle. Maximum entropy models (MEMs) find 𝑃(𝝈) by maximizing its entropy under the constraint that some empirical statistics are preserved (see Methods). As shown above, an interesting statistic for the present study is the coupling between the state of each binary signal, 𝜎 𝑖 , and the population activity 𝐾. The maximum entropy distribution that is consistent with the probability distribution 𝑃(𝐾), the average activations 〈𝜎 𝑖 〉, and the linear coupling between 𝜎 𝑖 and 𝐾, i.e., 〈𝜎 𝑖 𝐾〉 (which relates to 𝑧 𝑖 ),

	is given by the Boltzmann distribution 𝑃(𝝈) = 𝑒 -𝐸(𝝈) /𝑍, where 𝐸(𝝈) represents the energy of the
	pattern 𝝈, given as (29):
	𝑁
	𝐸(𝝈) = -∑(𝛼 𝑖 + 𝛽 𝐾 + 𝛾 𝑖 )𝜎 𝑖
	𝑖=1

Acknowledgements: APA, BJ, and GD received funding from the FLAG-ERA JTC (PCI2018-092891). GD acknowledges funding from the European Union's Horizon 2020 FET Flagship Human Brain Project under Grant Agreement 785907 HBP SGA1, SGA2 and SGA3, the Spanish Ministry Research Project PSI2016-75688-P (AEI/FEDER), the Catalan Research Group Support 2017 SGR 1545, and AWAKENING (PID2019-105772GB-I00, AEI FEDER EU) funded by the Spanish Ministry of Science, Innovation and Universities (MCIU), State Research Agency (AEI) and European Regional Development Funds (FEDER). MLK is supported by the ERC Consolidator Grant: CAREGIVING (n. 615539), Center for Music in the Brain, funded by the Danish National Research Foundation (DNRF117), and Centre for Eudaimonia and Human Flourishing funded by the Pettit and Carlsberg Foundations. BJ received funding from Fondation Bettencourt-Schueller,

Competing interests: No competing interests declared.

Author Contributions: APA designed research, analyzed the data, studied the models, and wrote the manuscript; LU performed the experiments and curated the data; ND studied the implementation of the models; CMS curated data; MK analyzed the data and provided data visualization codes; BJ designed and supervised the experiments; GD designed and supervised research. All authors discussed the results and contributed to the editing of the manuscript.

Supplementary Materials

Includes:

Appendix Supplementary Figures S1 to S8

Appendix Anesthesia Protocols

Monkeys received anesthesia either with propofol, ketamine, or sevoflurane [START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. Levels of anesthesia were defined using a clinical score and continuous electroencephalography monitoring.

Under ketamine, deep propofol anesthesia and deep sevoflurane anesthesia, monkeys stopped responding to all stimuli, reaching a state of general anesthesia [START_REF] Uhrig | Resting-state Dynamics as a Cortical Signature of Anesthesia in Monkeys[END_REF]. The 3 animals scanned at a deep level of ketamine anesthesia received an intramuscular (i.m.) injection of ketamine (20 mg/kg i.m., Virbac, France), followed by a continuous intravenous infusion of ketamine (15-16 mg/kg/h i.v.) to maintain anesthesia. For propofol anesthesia, monkeys were scanned under moderate propofol sedation and deep propofol anesthesia. The awake monkeys were injected with an intravenous (i.v.) propofol bolus (5-7.5 mg/kg i.v.; Fresenius Kabi, France) to induce anesthesia, followed by targetcontrolled infusion (Alaris PK Syringe pump, CareFusion, CA, USA) of propofol (moderate propofol sedation: 3.7-4.0 microg/ml; deep propofol anesthesia 5.6-7,2 microg/ml). During the ketamine and moderate propofol sedation, a muscle-blocking agent was administered (cisatracrium, 0.15 mg/kg bolus i.v. followed by continuous i.v. infusion at a rate of 0.18 mg/kg/h, GlaxoSmithKline, France)

to avoid artifacts related to potential movements during magnetic resonance imaging (MRI) acquisition. For sevoflurane anesthesia, monkeys were scanned under moderate and deep sevoflurane anesthesia. Monkeys received an intramuscular (i.m.) injection of ketamine (20 mg/kg i.m., Virbac, France) for induction of anesthesia, followed by sevoflurane anesthesia (moderate sevoflurane anesthesia: I/E: 2,2/2,1 vol% or deep sevoflurane anesthesia: I/E: 4,4/4,0 vol%) (Abbott, France). For all the anesthesia experiments, monkeys were intubated and ventilated. Heart rate, non-invasive blood pressure (systolic/diastolic/mean), oxygen saturation (SpO2), respiratory rate, end-tidal CO2 (EtCO2), and cutaneous temperature were monitored (Maglife, Schiller, France) and recorded online (Schiller, France).

fMRI Data Acquisition

Monkeys were scanned on a 3-T horizontal scanner (Siemens Tim Trio; TR, 2,400 ms; TE, 20 ms; and 1.5-mm 3 voxel size; 500 brain volumes per scan session). Before each scanning session, a contrast agent monocrystalline iron oxide nanoparticle (MION, Feraheme; AMAG Pharmaceuticals; 10 mg/kg, i.v.) was injected into the monkey's saphenous vein [START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF]. Preprocessing of functional images followed the standard steps described in [START_REF] Barttfeld | Signature of consciousness in the dynamics of resting-state brain activity[END_REF], normalized to the anatomical template of the monkey Supplementary Figures For example, in panel (A), the element (𝑘, 𝑙) of the matrix represents the correlation between the coupling to population vector 𝒛 of scans 𝑘 and 𝑙. Coupling to population clearly separated awake and anesthesia data. Using k-means, we evaluated how well the different statistics could be used to classify the awake and anesthetized conditions (chance level: 50%). The classification performance using the coupling to population statistic was 97.48%, that was significantly higher than using the functional connectivity (84.87%). Classification of the six experimental conditions was generally lower, but higher for 𝑧 than for functional connectivity (48% vs. 38.11%, chance level: 16.67%).