# Spontaneous cortical activity is transiently poised close to criticality

Gerald Hahn, Adrián Ponce-Alvarez, Cyril Monier, Giacomo Benvenuti, Arvind Kumar, Frédéric Chavane, Gustavo Deco, Yves Frégnac

► To cite this version:

Gerald Hahn, Adrián Ponce-Alvarez, Cyril Monier, Giacomo Benvenuti, Arvind Kumar, et al.. Spontaneous cortical activity is transiently poised close to criticality. PLoS Computational Biology, 2017, 13 (5), pp.e1005543. 10.1371/journal.pcbi.1005543. hal-02008357

HAL Id: hal-02008357

https://amu.hal.science/hal-02008357

Submitted on 21 May 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Spontaneous cortical activity is transiently poised close to criticality

Gerald Hahn<sup>1,2,*</sup>, Adrian Ponce-Alvarez<sup>2,*</sup>, Cyril Monier<sup>1</sup>, Giacomo Benvenuti<sup>3</sup>, Arvind Kumar<sup>4,5</sup>, Frédéric Chavane<sup>3</sup>, Gustavo Deco<sup>2,6,7,8</sup>, Yves Frégnac<sup>1</sup>

1 Unité de Neuroscience, Information et Complexité (UNIC), CNRS, Gif-sur-Yvette, France, 2 Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain, 3 Institut des Neurosciences de la Timone, CNRS, Marseille, France, 4 Bernstein Center for Computational Neuroscience, Freiburg, Germany, 5 Dept. of Computational Science and Technology, School of Computer Science and Communication, KTH, Royal Institute of Technology, Stockholm, Sweden, 6 Institució Catalana de la Recerca i Estudis Avançats, Universitat Pompeu Fabra, Barcelona, Spain, 7 Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany, 8 School of Psychological Sciences, Monash University, Melbourne, Clayton, Victoria, Australia

© These authors contributed equally to this work.

* gerald.hahn@upf.edu

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Hahn G, Ponce-Alvarez A, Monier C, Benvenuti G, Kumar A, Chavane F, et al. (2017) Spontaneous cortical activity is transiently poised close to criticality. PLoS Comput Biol 13(5): e1005543. https://doi.org/10.1371/journal. pcbi.1005543

**Editor:** Claus C. Hilgetag, Hamburg University, GERMANY

**Received:** December 8, 2016

**Accepted:** April 26, 2017

**Published:** May 24, 2017

**Copyright:** © 2017 Hahn et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** GH, CM and YF were supported by the CNRS, the Agence Nationale de la Recherche (ANR: V1-Complex) https://www.cnrs.fr. GH and GB were financed by the initial training network program FACETS-ITN (PITN-GA-2009-237955) http://facets.kip.uni-heidelberg.de/ITN/. APA was supported by SEMAINE ERA-Net NEURON Project and by a Juan de la Cierva fellowship (IJCI-2014-

## Abstract

Brain activity displays a large repertoire of dynamics across the sleep-wake cycle and even during anesthesia. It was suggested that criticality could serve as a unifying principle under- lying the diversity of dynamics. This view has been supported by the observation of sponta- neous bursts of cortical activity with scale-invariant sizes and durations, known as neuronal avalanches, in recordings of mesoscopic cortical signals. However, the existence of neuro- nal avalanches in spiking activity has been equivocal with studies reporting both its presence and absence. Here, we show that signs of criticality in spiking activity can change between synchronized and desynchronized cortical states. We analyzed the spontaneous activity in the primary visual cortex of the anesthetized cat and the awake monkey, and found that neu- ronal avalanches and thermodynamic indicators of criticality strongly depend on collective synchrony among neurons, LFP fluctuations, and behavioral state. We found that synchro- nized states are associated to criticality, large dynamical repertoire and prolonged epochs of eye closure, while desynchronized states are associated to sub-criticality, reduced dynam- ical repertoire, and eyes open conditions. Our results show that criticality in cortical dynam- ics is not stationary, but fluctuates during anesthesia and between different vigilance states.

## Author summary

Cortical activity spontaneously displays a large diversity of dynamics in different behav- ioral states and at multiple spatiotemporal scales. In the last decades, a unifying principle has been proposed to underlie the diversity and scale-invariance of brain dynamics: Criti- cality, a particular state of complex systems in which order and disorder coexist. On the other hand, the cortex can exhibit a continuum of states with different levels of collective synchrony and LFP fluctuations. Here, we ask how criticality measures vary as the cortex
21066) from the Spanish Ministry of Economy and Competitiveness. CM, GD and YF received funding from the EC grants BrainScales (FP7-2010-IST-FETPI 269921) and the flagship Human Brain Project (n.604102) https://www.humanbrainproject.eu/. The Utah array recordings were made possible through a loan by S.Grün (Research Center Jülich, INM6, Germany) and were part of a collaborative work with S. Grün and A. Riehle (INT, Marseille). AK received funding from the German Federal Ministry of Education and Research (BMBF 01GQ0420 to BCCN Freiburg and 01GQ0830 to BFNT Freiburg/Tübingen) https://www.bmbf.de/en/. GD is supported by the ERC Advanced Grant: DYSTRUCTURE (n. 295129), by the Spanish Research Project PSI2016-75688-P and by the European Union's Horizon 2020 research and innovation programme under grant agreement n. 720270 (HBP SGA1). GD obtained support from the ERC Advanced Grant DYSTRUCTURE (n. 295129) gustavodecolab.com/dystructure/ and the Spanish Research Project PSI2013-42091- P. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

spontaneously fluctuates across different states. Using recordings from the primary visual cortex of the anaesthetized cat and awake monkey, we show that spiking patterns of cortical states differ in their proximity to criticality and their relation to the vigilance state.

## Introduction

The cortex continuously generates coordinated ongoing patterns of activity during varying behavioral states such as restful wakefulness [1], focused attention [2], sleep and anesthesia [3]. Patterned spontaneous activity has been reported at multiple scales as revealed by multi-electrode arrays [4], optical imaging [5,6] magnetencephalography [7] and functional MRI [1]. In spite of early findings that the statistics of cortical activity are highly variable [8], interest arose to find an overarching physical principle that explains in a unified way the dynamical organization of the brain. Theoretical reasoning made criticality, a concept borrowed from statistical physics, a plausible candidate for such a unifying concept as it has been proposed to account for the brain's inherent complexity necessary to process and represent its environment [9]. In the critical regime, cascading dynamics are expected to produce a great diversity of transient co-activations or "neuronal avalanches" [10,11]. Within this framework, power laws and high levels of long-range correlations are two strong indicators of critical dynamics [11,12].

Even though first reported in *in vitro* preparations [13–18], evidence for criticality in the cortex was found during anesthesia [19,20], sleep [21], resting state in animals and humans [22–27], task-related conditions [28] and stimulus evoked activity [29,30]. However, other studies failed to find the typical signs of criticality, either refuting the neuronal avalanche hypothesis [31–33] or explaining the negative result by subsampling of neuronal activity due to fundamental measurement limitations of current multi-electrode recording techniques [21,28,34]. Alternatively, criticality or its absence in the cortex may not be a constant feature of brain activity, but fluctuates with the ever-varying dynamics of the cortex and associated arousal states. Indeed, the statistics of neuronal population spike activity in the cortex substantially vary over time, continuously moving between highly synchronized burst dynamics and desynchronized, irregular activity [35]. Synchronized cortical states are the hallmark of *in vitro* preparations [13], but also appear during slow wave sleep [36], anesthesia [21,37,38] and also in awake animals during quiet waking and drowsiness [39,40]. In contrast, desynchronized activity is associated with active behavior [41,42], responses to visual stimuli [43] and suggested to be a neuronal correlate of attention [35]. In line with this hypothesis, putative transitions between critical, supercritical and subcritical dynamics have been observed during anesthesia [20,44], shifts between anesthesia and the physiological sleep-wake cycle [21,22,45–48], prolonged waking [49], and in a rest versus task setting [50]. Here, we tested whether the signature of criticality at the level of collective spiking activity remains robust across different cortical states or fluctuates together with changing cortical dynamics recorded in the primary visual cortex of the anesthetized cat and awake monkey using classical criticality markers and data modeling. We found time varying criticality properties which were associated with the local cortical state, both during anesthesia and in the awake state, and changed with prolonged epochs of eye closure in the awake condition.

## Results

In this study, we recorded spontaneous spiking activity and local field potentials (LFP) from the primary visual cortex of four anesthetized cats (~1 hr/animal) and one awake monkey using 32-channel silicon probes (cats) and a chronically implanted 96-channel Utah array
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Fig 1.** Separation of cortical states in spontaneous activity of anesthetized cat and awake monkey. (A) LFP spectrograms of two 100s segments computed with non-overlapping windows of 1s. Bottom: colored bars indicate cortical state as defined in the main text. (B) Coefficients for first three principal components as a function of power spectrum frequency. (C) Principal component space for two entire datasets (cat: 6000s, monkey: 600s). Each circle represents a data segment of 1s duration. Colors indicate different cortical states. (D) Dunn index as a function of the cluster number extracted by k-means. (E) Average (+SD) duration of different states across all datasets of a species. (F) Average power spectrum of different cortical states for all cat and monkey datasets. Dashed lines indicate standard deviation (±SD). Inset: same as in main figure, but in log-log coordinates to show peak in alpha band.

https://doi.org/10.1371/journal.pcbi.1005543.g001

(monkey, 4x10 minutes, total of 40 minutes). To avoid interference of visually evoked activity with the ongoing dynamics all recordings were performed in the dark. Visual inspection of the LFP spectrogram of both cat and monkey datasets showed non-stationary dynamics with prominent power fluctuations in a frequency band between 1 and 15 Hz and stronger alpha oscillations (9–13 Hz) in the monkey recordings (Fig 1A).

## State separation

In order to account for the variability of cortical dynamics, we attempted to pool short segments of the data into sets which were characterized by a similar frequency composition of the
LFP power spectrum. Each set is referred to as a different cortical (or dynamical) state in accordance with previous studies [35,38,41].

We first used a combination of principal component analysis (PCA) and subsequent k-means clustering to find frequency bands that showed the highest variance across an entire dataset and subsequently cluster data epochs with similar power in these bands (see Materials and methods). The first three principal components (PC) explained 97% (cats) and 99% (monkey) of the variance and corresponded to different bands of the frequency spectrum within a range between 1 and 15 Hz (Fig 1B). Each data segment was represented in a three dimensional PC space and we observed a characteristic pattern of segment clusters across cat and monkey datasets. Most segments were concentrated at the origin of the PC space around which the remaining ones were scattered along the three dimensions (Fig 1C). We performed a k-means cluster analysis to extract different clusters in each dataset. The optimal number of clusters was estimated using a validation procedure based on the Dunn index (DI). On average, the optimal number of clusters, as indicated by a peak in the DI distribution, was found to be five in the cat datasets and four in the monkey recordings (Fig 1D), each of which corresponded to different cortical states. In order to compare states across different datasets and species, we coded clusters at the same location in the PC space with the same color (see color code in Fig 1C). Analysis of cluster duration revealed that the cortex spent most of the time in the cluster centered around the origin of the PC space (58 ± 0.07% and 75 ± 0.17%, mean ± SD of the data in the cat and monkey, respectively), while cortical dynamics represented in the other clusters was much less frequent (Fig 1E).

Next, we computed average power spectra across all data segments in one cluster and the analysis results confirmed that the separation algorithm extracted power spectra with a similar frequency composition across different time epochs (Fig 1F). The black clusters consistently corresponded to dynamics with little power in lower frequencies (1–5 Hz) as opposed to the red and green clusters in which the power in slower (red) or faster frequencies (5–15 Hz, green) was high. The blue and magenta clusters had intermediate power at slower frequencies approaching either the black or red cluster, respectively. In the subsequent analysis we will refer to states with relatively lower power in slow frequencies as “desynchronized states”—denoted as Desyn I (black) and Desyn II (blue, only for cat)—, with higher power in slow frequencies as “slow synchronized states”—denoted as SynSlow I (magenta) and SynSlow II (red)—and with high power in faster frequencies as “fast synchronized state”—denoted as SynFast (green).

## Characteristics of spiking activity across cortical states

Spiking activity consistently varied with the cortical states as defined by the LFP (see Fig 2A for an example taken from a cat recording). During epochs of desynchronized LFP with small amplitude fluctuations, spiking activity was characterized by continuous and irregular firing of neurons without apparent synchronization across the different channels. When the LFP shifted to dynamics with larger amplitude fluctuations, the underlying firing pattern of neurons was characterized by bursts of spiking activity that occurred synchronously across most of the channels and was followed by periods of silence.

The mean channel firing rates consistently varied across different states (Fig 2B). Neuronal activity in the anesthetized cat showed a trend to increased firing rates as LFP entered into synchronized states, but differences across states were not significant (one-way rm-ANOVA: $p > 0.11$). In contrast, neurons in the awake monkey fired significantly more vigorously in desynchronized states than during synchronized activity (one-way rm-ANOVA: $F_{3,9} = 39.03$, $p = 0.0005$, $\epsilon = 0.62$).
A

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

B

<figure><img src="image_12.png" /></figure>

C

<figure><img src="image_13.png" /></figure>

D

<figure><img src="image_14.png" /></figure>

E

<figure><img src="image_15.png" /></figure>

F

<figure><img src="image_16.png" /></figure>

**Fig 2. Characteristics of LFP-spike relationship and spiking activity for different cortical states.** (A) LFPs (top) and spike rasters (bottom) for 1s segments of a desynchronized and synchronized cortical state from one cat dataset. Spike counts were computed with a Gaussian kernel (20ms window-size). (B) Average (+SD) population firing rate for all recordings of a species. (C) Examples of spike-triggered averages (STA) for different cortical states. (D) Average (+SD) STA area across all cat and monkey datasets. The area is normalized to the state with the largest value within a dataset. (E) Autocorrelation histograms of the population spike trains (pACH) for all cortical states in one cat and one monkey recording. (F) Average (+SD) area and of ACH peaks computed for all cat and monkey datasets. The values were normalized to the states with maximum area within a dataset.
Horizontal bars: significant Bonferroni multiple comparison test (p<0.05).

https://doi.org/10.1371/journal.pcbi.1005543.g002

To further quantify the relationship between spikes and LFP we computed spike-triggered averages (STA) for each dataset and cortical state (see Fig 2C for two examples of a cat and a monkey dataset). In all cases, as expected, spikes were associated with the negative trough of LFP deflections, whose area between the negative deflection and the zero baseline (one-way rm-ANOVA: $F_{4,12} = 7.26$, p = 0.019, $\epsilon = 0.56$ (cat) and $F_{3,9} = 6.63$, p = 0.027, $\epsilon = 0.71$ (monkey)) varied between states. In desynchronized states the LFP peaks were small and thin, while they increased in amplitude and width in more synchronized states. The synchronized cluster with faster deflections displayed peaks with the highest amplitude, but its width was smaller than synchronized states with lower frequency fluctuations. These results were confirmed across all datasets and both species (Fig 2D).

To assess the level of correlation in the population spiking activity, we computed the population auto-correlation histograms (pACH) of the population spike trains for each state (Fig 2E) which captures all possible correlations either generated intrinsically by neurons (i.e. bursting behavior) or through neuronal interactions. The area of pACH peaks was significantly state dependent (one-way rm-ANOVA: $F_{4,12} = 15.91$, p = 0.002, $\epsilon = 0.6$ (cat) and $F_{3,12} = 4.27$, p = 0.018, $\epsilon = 1$ (monkey)). The pACH of the Desyn I state had lower amplitude and area
compared to the other states, while the pACHs of SynSlow II and SynFast states had the largest peaks. The pACHs of the slow synchronized states had the largest spread, while the temporal correlations were shortened in the other states. In addition, spike synchronization in faster frequency bands (alpha range) was clearly visible in the pACH of the monkey, especially in the SynFast state. These results were consistent across all datasets and species (Fig 2F).

## Cortical states and criticality

To study criticality in spiking activity, we first applied a classical analysis [13] that is based on binning point processes like spike trains and extracting clusters defined as consecutive bins with a number of spikes above a predefined threshold (Fig 3A, threshold ≥ 1). Bin-sizes were given by the average inter-spike interval of each state (see Table 1 for detailed values). Sufficiently sampled critical dynamics underlying the generated spikes predicts power laws in the cluster size distributions, i.e. a straight line in log-log plots. Such analysis was consequently performed for each cortical state separately and the full datasets for comparison (Fig 3B). The distributions in the most desynchronized states (black solid lines) did not follow a power law, albeit not being as curved as expected from shuffled data (Fig 3B, dashed lines). However, moving from desynchronized to more synchronized states, the tail of the cluster size distribution increased and came closest to a power law in the most synchronized states (red lines). This pattern was consistent in both cat and monkey recordings. Analysis of the full datasets revealed clearly curved distributions similar to the desynchronized state in accordance with previous studies in spiking data [21, 31, 32, 51] (Fig 3B, solid black line).

To quantify these tails, we fitted both lognormal and power law distributions to the cluster size distributions. Lognormal functions are more flexible and support different degrees of hea-vytaildness as opposed to power laws with extreme tails. Thus, comparing the goodness of fit between both distributions provides a relative quantification of how closely a distribution approaches a power law. We performed comparisons using loglikelihood ratios (LLR) with positive values indicating superior power law fit and negative values giving more evidence for lognormal fits with reduced tails [22,52]. In our recordings, LLRs were in general negative across all states in both cat and monkey datasets (all datasets: p < 0.001), indicating that tails were less pronounced than expected by a power law. However, how much each state deviated from a power law varied between states as indicated by the absolute values of the LLR. These were more negative during desynchronized states and when analyzing the full datasets, while the discrepancy between lognormal and power law fit was less pronounced during states with more synchronization (Fig 3C). We repeated this analysis using different thresholds for ava-lanche detection (see Methods and material) [22,53]. In the cat, significant state differences (rm-ANOVA: p < 0.05) were found for all but the lowest and highest thresholds, which only yielded a trend (p < 0.1) or were non-significant, respectively. In the monkey, state differences were found for all thresholds (see example in Fig 3D and S1 Fig for a detailed analysis). We also determined the exponent of the fitted power law distributions for all states which yielded values for each species close to the theoretically expected and experimentally found value [13,54] of -1.5 (cats: -1.62 ± 0.05, monkey: -1.58 ± 0.03), desynchronized states included.

Even though power laws can be predictive of criticality, a failure to detect straight lines does not entirely rule out absence of criticality [21,28]. Another test for criticality is the existence of universal scaling functions that capture the system dynamics at different scales. In cortical cultures, it has been observed that the time-courses of long avalanches are scaled versions of short avalanches as predicted close to a critical point [15]. Such invariance across scales, or self-similarity, known as “shape collapse”, was also found in spiking data of rats undergoing all stages of the sleep-wake cycle and resting state EEG recordings of humans [21,50]. Here, we used the
A

<figure><img src="image_17.png" /></figure>

B

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

C

<figure><img src="image_20.png" /></figure>

D

<figure><img src="image_21.png" /></figure>

**Fig 3. Neuronal avalanche analysis of spiking activity.** (A) Spike clusters were defined as a sequence of bins containing spikes ≥1 threshold (in this example threshold = 1). The size of a spike cluster is given by the total number of spikes within a cluster. The lifetime of a cluster is defined as the number of bins. Δt was chosen as the average ISI<sub>pop</sub> interval of the population spike train for each state. (B) Cluster size distributions of different states. Dotted black lines indicate power law with exponent = -1.5. Dashed lines with gray squares represent an inhomogenous Poisson process created from the desynchronized I state by spike time randomization within all 1s segments. Dashed lines with empty squares indicate homogenous Poisson process with the same rate and duration as the entire desynchronized I state. (C) Loglikelihood ratios for power law and lognormal fits to cluster size distributions of different cortical states across all cat (threshold: 3 spikes) and monkey data (threshold: 2 spikes). Negative values indicate a better lognormal fit. Horizontal bars: significant Bonferroni multiple comparison test. (D) Significance between different states expressed as the p-value of an rm-ANOVA test for different thresholds defining a cluster.

https://doi.org/10.1371/journal.pcbi.1005543.g003
<table><caption>Table 1. Average inter-spike intervals (in ms) of the population spike train for all cat and monkey datasets.</caption><thead><tr><th rowspan="2">Cortical State:</th><th colspan="4">Cat</th><th colspan="4">Monkey</th></tr><tr><th>Dataset 1</th><th>Dataset 2</th><th>Dataset 3</th><th>Dataset 4</th><th>Dataset 1</th><th>Dataset 2</th><th>Dataset 3</th><th>Dataset 4</th></tr></thead><tbody><tr><td>Desyn I</td><td>14</td><td>43</td><td>28</td><td>30</td><td>27</td><td>14</td><td>21</td><td>16</td></tr><tr><td>Desyn II</td><td>26</td><td>31</td><td>23</td><td>43</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>SynFast</td><td>16</td><td>47</td><td>26</td><td>55</td><td>21</td><td>24</td><td>19</td><td>29</td></tr><tr><td>SynSlow I</td><td>56</td><td>38</td><td>27</td><td>22</td><td>25</td><td>22</td><td>30</td><td>19</td></tr><tr><td>SynSlow II</td><td>16</td><td>54</td><td>18</td><td>37</td><td>24</td><td>15</td><td>31</td><td>28</td></tr><tr><td>Full</td><td>28</td><td>48</td><td>26</td><td>45</td><td>23</td><td>15</td><td>21</td><td>16</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1005543.t001

method of Marshall et al. [55] to examine whether avalanches time-courses are self-similar for the different states and show shape collapse (for details see Supplementary Material). Since this analysis requires a large amount of data, we restricted our analysis to the cat data. The goodness of collapse was quantified by the collapse index (CI), defined as the variance across rescaled avalanche profiles. We found that the empirical avalanche profiles collapse significantly (p = 0.027) more for the synchronized states (SynSlow I/II) than for the desynchronized states (Desyn I/II) (S3A-S3F Fig). These results are consistent with the power law analysis showing that the presence of curved avalanche size distributions is accompanied by a reduction in self-similarity and indicate stronger deviations from a critical state in desynchronized cortical dynamics.

## Predicting criticality from LFP fluctuations and spike correlations

Next, we attempted to synthesize the above findings by showing how different measures of LFP and spiking activity are linked and relate to criticality in neuronal discharges. To this end, we quantified the correlation between the loglikelihood ratio fitted to all spike cluster size distributions, and the mean autocorrelation between neurons (i.e. area of the population ACH, representing the total amount of correlation in the neuronal population) as well as the amount of spike-triggered LFP fluctuations (i.e. the average STA area across channels). As the STA and ACH peaks became larger, the tails of the cluster size distributions and the associated loglikelihood ratios significantly increased or followed a trend (Fig 4). This linear dependence indicates that the size of the STA and the level of population autocorrelation in spiking activity are directly translated into the degree of tailedness of the spike cluster size distributions. As the LFP varied with cortical state, the results of the neuronal avalanche analysis changed from more power law like to more curved distributions indicating the absence of correlated population activity.

## Lifetime distributions

Another signature of criticality is a power law in the distribution of lifetimes, i.e. the duration of neuronal event clusters (Fig 5A), which decays with a characteristic exponent of -2. This theoretically obtained result [54] was indeed found in neuronal networks in vitro and in vivo using LFP events and spiking activity [13,22]. Here, we also calculated lifetime distributions and assessed their closeness to power laws by comparing LLRs between power laws and lognormals across different states. Similar to the cluster size distributions, LLRs indicated a generally better lognormal fit (all datasets: p < 0.001). However, LLR values during synchronized states were significantly less negative during synchronized states compared to desynchronized states and the full datasets (Fig 5B). These differences were significant (rm-ANOVA: p<0.05) for thresholds one and three in the cat, and for all thresholds in the monkey (Fig 5C and
<figure><img src="image_22.png" /><figcaption>**Fig 4.** Correlation between LFP, spike synchronization and the level of criticality measured by the loglikelihood ratio for power law and lognormal fits. Top: LLR as a function of STA area for all cat and monkey datasets. Bottom: LLR as a function of ACH area.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005543.g004

S2 Fig). Exponents of fit power law distributions remained unchanged across different states and approached a value of -2 (cats: -1.84 ± 0.06, monkey: -1.74 ± 0.04), as predicted by theory and other experiments [13,54].

## Inter-spike interval distributions

Previous studies reported heavy tails in the inter-spike interval distributions of individual neurons [56,57] and spike trains of neuronal populations [20,58]. In this study, we also compared inter-spike interval distributions of population spike trains (ISI<sub>pop</sub>) for each state quantifying differences through the coefficient of variation (CV = std(ISI<sub>pop</sub>) / mean(ISI<sub>pop</sub>)). In both cat and monkey recordings the distributions of desynchronized states and the entire datasets were clearly curved (Fig 5D) and approached a CV of ~1 (Fig 5E), which is expected for independent Poisson processes. In contrast, synchronized states were associated with heavier tails as expressed by significantly larger CVs (one-way rm-ANOVA: F<sub>4,12</sub> = 9.66, p = 0.03, ε = 0.35 (cats) and F<sub>4,12</sub> = 18.18, p < 0.001, ε = 1 (monkey)).

## Temporal evolution of cortical states

Having identified different cortical states and their statistical properties, we next studied how synchronization in the visual cortex evolved over time in the cat, where we had acquired continuous recordings up to one-hour duration. In general, the cortex not only fluctuated between the different states at a time scale of seconds, but the probability to find synchronized or desynchronized states also changed at much longer time scales (Fig 6A). These fluctuations appeared at very slow timescale (<0.01 Hz) with a peak at ~ 0.001–0.002 Hz (Fig 6B), indicating the presence of ultra-slow dynamics during anesthesia with long periods of predominant synchronization or desynchronization.
A

<figure><img src="image_23.png" /></figure>

B

<figure><img src="image_24.png" /></figure>

C

<figure><img src="image_25.png" /></figure>

D

<figure><img src="image_26.png" /></figure>

E

<figure><img src="image_27.png" /></figure>

**Fig 5. Lifetime and ISI<sub>pop</sub> distributions of spiking activity.** (A) Lifetime distributions of different states. Dotted black lines indicate power law with exponent = -2. (B) Loglikelihood ratios for power law and lognormal fits to lifetime distributions of different cortical states across all cat (threshold: 3 spikes) and monkey data (threshold: 2 spikes). Negative values indicate a better lognormal fit. Horizontal bars: significant Bonferroni multiple comparison test. (C) Significance between different states calculated as the p-value of an rm-ANOVA test for varying thresholds constituting a cluster. (D) ISI<sub>pop</sub> distributions of population spike trains for one cat
and one monkey recording. (E) Coefficient of variation for ISI<sub>pop</sub> distributions of different states across all cat and monkey datasets.

https://doi.org/10.1371/journal.pcbi.1005543.g005

In the monkey, short periods of synchronization (several seconds) alternated with longer temporal epochs of desynchronization (Fig 1A). As the monkey was resting in the dark, we hypothesized that state changes coincided with the monkey becoming drowsy and closing the eyes. Eye tracking analysis revealed that the animal indeed closed the eyes during extended periods of time, on average 8 ± 4.8s (Fig 6C). These epochs of eye closure were associated with the cortex shifting to synchronized states (Fig 6D). In particular the occurrence of fast (p < 0.001) and slower synchronized (SynSlow II) states (p < 0.001) was highly correlated with the closed eye condition, while desynchronization showed significant anticorrelation (p < 0.001) (Fig 6E). This confirms that the level of spike synchronization and hence criticality is closely related to the vigilance state of the animal.

## Statistics of the joint neuronal activity patterns

Next, we studied the patterns of the ensemble spiking activity from *N* recording sites during the different cortical states and examined their statistics in the framework of statistical

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

**Fig 6.** Evolution of cortical states over time. (A) LFP spectrogram of an entire recording in one cat calculated with non overlapping windows of 1s duration (top) and the probability of synchronized states (synchronized slow I, II and synchronized fast states) computed with a sliding window of size 100s and an overlap of 1s (bottom). (B) Average power spectrum (SEM: dashed lines) of the time courses of synchronized states across all four cats. (C) Probability to find periods of eye closure with a given duration across all monkey datasets. (D) Time course of synchronized states (synchronized slow II and synchronized fast) (top) and periods of full eye closure (bottom) in one monkey dataset. Vertical bars represent 1s time segments. (E) Pearson correlation coefficient calculated between the time course of eyes closure and the time courses of different cortical states.

https://doi.org/10.1371/journal.pcbi.1005543.g006
mechanisms. This enabled us to have further physical characterization of the different cortical states in terms of criticality. Time was discretized in non-overlapping bins. In a time bin $\Delta t$, the state of the neural ensemble is described by a binary pattern $\vec{\sigma} = [\sigma_1, \sigma_2, ..., \sigma_N]$, where $\sigma_i = +1$ if the $i$-th electrode site recorded one or more spikes (activation) and $\sigma_i = -1$ otherwise (quiescence). $\Delta t$ was chosen to obtain homogeneous averaged activation rates across datasets and cortical states (Fig 7A, see also Table 2). We aimed at obtaining the probability distribution $P(\vec{\sigma})$ over the $2^N$ possible binary patterns. For this, we used a Maximum entropy model (MEM) to find $P(\vec{\sigma})$ by maximizing its entropy under the constraint that the activation rates ($\langle\sigma_i\rangle$) and the pairwise correlations ($\langle\sigma_i\sigma_j\rangle$) are preserved (see Methods). The resulting distribution is $P(\vec{\sigma}) \propto e^{-E(\vec{\sigma})}$, where $E(\vec{\sigma})$ is the energy of the pattern $\vec{\sigma}$, and $E(\vec{\sigma}) = -\sum_i^N [h_i\sigma_i + \frac{1}{2}\sum_j^N J_{ij}\sigma_i\sigma_j]$, where $h_i$ represents the intrinsic tendency of site $i$ towards activation and $J_{ij}$ represents the effective interaction between sites $i$ and $j$ (see Methods). Note that energies are proportional to the patterns' minus log probabilities, or "surprise". The model parameters $\Omega = \{h; J\}$ were estimated from the data using a gradient descent algorithm

<figure><img src="image_34.png" /><figcaption>**Fig 7.** Maximum entropy models (MEMs) of different cortical state. (A) Probability that an electrode site $i$ has $\sigma_i = +1$ (emission probability) using the bin sizes indicated in Table 2 (error bars indicate SEM). The emission probabilities do not significantly depend on the cortical state ($p = 0.092$ for cat data, $p = 0.058$ for monkey data, rm-ANOVA). (B) Goodness-of-fit ($1/D_{JS}$) of pairwise-MEMs (filled bars) and independent-MEMs (open bars) for each cortical state (averaged over the 10 groups of $N$ signals and all datasets; error bars indicate SEM). (C) Prediction of the cortical state using pairwise- and independent-MEMs. The percentage of correct classifications is shown for each cat (left) and monkey (right) dataset. Squares indicate the medians and error bars delineate the 5–95th percentiles of the classification performance. Dash lines: mean and 95th percentile of the number of correct classifications expected by chance. \*: significant classification performances ($p < 0.05$). (D–E) Effect of changing the temperature parameter $T$ on the model activity. The MEM in this example was estimated using the activity of a neural ensemble of the cat in the SynSlow II state. **D**: model activity (1000 steps are shown out of $10^6$ steps). At low temperature ($T = 0.5$) the activity is sparse and correlations are low ($\langle r_c \rangle = 0.07$); at high temperature ($T = 2$) the activity is dense and random and correlations are low ($\langle r_c \rangle = 0.06$); for an intermediate temperature ($T = 1$) the activity is more patterned and correlations are higher ($\langle r_c \rangle = 0.10$). **E**: Left, occupied energy levels. The size of the horizontal lines is proportional to $\log(n_E)$, where $n_E$ is the number of patterns that have the energy $E$. Right, the entropy $S(T)$ increases with $T$ and the heat capacity $C(T)$ peaks at a given temperature $T_{max}$ (for this particular example neural ensemble $T_{max} = 1$). **F–G**: Heat capacity as a function of the temperature parameter ($T$), for each cortical state, for two example anesthetized cat datasets (**F**) and two example awake monkey datasets (**G**) (10 random choices of groups of $N$ signals were used; trace thickness indicates SEM). **Inset**: Peak temperature ($T_{max}$) for each neuronal ensemble and for each cortical state. $T_{max}$ significantly depends on cortical state (p: p-value, rm-ANOVA). (H–I) $T_{max}$ averaged over all cat datasets (**H**) and over all monkey datasets (**I**) for each cortical state (rm-ANOVA: H: $F_{4,156} = 25.03$, $p < 0.001$, $\varepsilon = 0.34$; I: $F_{3,117} = 59.60$, $p < 0.001$, $\varepsilon = 0.90$). Error bars indicate SEM. Horizontal bars: significant differences ($p < 0.05$) of subsequent Bonferroni's test for multiple comparisons. In A–I the size of the models was $N = 6$.</figcaption></figure>
<table><caption>Table 2. Total number of active electrode sites (n), number of observed binary patterns (L), and bin size (b) for each cortical state and each dataset.</caption><thead><tr><th rowspan="3">Cortical State:</th><th colspan="12">Cat</th><th colspan="12">Monkey</th></tr><tr><th colspan="2">Dataset 1<br>(n = 27)</th><th colspan="2">Dataset 2<br>(n = 25)</th><th colspan="2">Dataset 3<br>(n = 22)</th><th colspan="2">Dataset 4<br>(n = 12)</th><th colspan="2">Dataset 1<br>(n = 6)</th><th colspan="2">Dataset 2<br>(n = 9)</th><th colspan="2">Dataset 3<br>(n = 10)</th><th colspan="2">Dataset 4<br>(n = 10)</th></tr><tr><th>L</th><th>b</th><th>L</th><th>b</th><th>L</th><th>b</th><th>L</th><th>b</th><th>L</th><th>B</th><th>L</th><th>b</th><th>L</th><th>b</th><th>L</th><th>b</th></tr></thead><tbody><tr><td>Desyn I</td><td>38,379</td><td>50</td><td>70,248</td><td>50</td><td>77,395</td><td>50</td><td>77,395</td><td>50</td><td>34,700</td><td>10</td><td>51,499</td><td>10</td><td>36,996</td><td>10</td><td>54,199</td><td>10</td></tr><tr><td>Desyn II</td><td>21,259</td><td>50</td><td>25,080</td><td>50</td><td>25,080</td><td>50</td><td>10,440</td><td>50</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td><td>--</td></tr><tr><td>SynFast</td><td>6,860</td><td>50</td><td>14,360</td><td>50</td><td>14,360</td><td>50</td><td>5,020</td><td>50</td><td>6,191</td><td>20</td><td>1,300</td><td>20</td><td>4,350</td><td>20</td><td>650</td><td>20</td></tr><tr><td>SynSlow I</td><td>7,640</td><td>50</td><td>2,880</td><td>50</td><td>2,880</td><td>50</td><td>4,280</td><td>50</td><td>3,850</td><td>20</td><td>2,200</td><td>20</td><td>2,150</td><td>20</td><td>1,750</td><td>20</td></tr><tr><td>SynSlow II</td><td>3,860</td><td>50</td><td>7,380</td><td>50</td><td>7,380</td><td>50</td><td>22,860</td><td>50</td><td>2,600</td><td>20</td><td>750</td><td>20</td><td>5,000</td><td>20</td><td>500</td><td>20</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1005543.t002

(see Methods). We constructed pairwise-MEMs for each dataset and each cortical state separately and analyzed the obtained models. Here, we present analyses of the joint activity of $N = 6$ electrodes for both cat and monkey datasets—but, as shown below, we found consistent results for models of larger size for the cat datasets, for which larger ensembles were recorded (see Table 2).

We first evaluated the importance of correlations by comparing the performance of the pairwise-MEMs to the performance of independent-MEMs for which only the activation rates ($<\sigma_i>$) are preserved (i.e., only *h* is optimized). Specifically, we compared the prediction and the empirical estimation of the probabilities of spiking patterns in each cortical state. We found that taking into account the pairwise correlations improves the agreement with the data by 1–2 orders of magnitude (Fig 7B). Second, we investigated how specific the different MEM parameters are to the different cortical states by measuring the ability of the model to predict the cortical state from neural ensemble activity. We used a jackknife cross-validation procedure, consisting of first estimating the MEM parameters on a subset of the available data for each cortical state and then using the learned MEM parameters to classify the remaining data (Methods). We found that, on average, pairwise-MEMs predicted the cortical state of the cat data with 75.8–98.7% and 74.8–99.7% accuracy for $N = 6$ and $N = 12$, respectively, and of the monkey data with 73.1–83.4% accuracy for $N = 6$ (Fig 7C). However, when using independent-MEMs the performance hardly exceeded chance level. These results confirm that the MEM parameters, specifically the interaction couplings $J_{ij}$, tightly relate to the cortical state.

Next, we further studied the learned MEMs to obtained relevant features of the collective dynamics in each cortical state. One important quantity is the heat capacity, $C$, that relates to how the distribution of energies changes as a function of a parameter $T$, analogous to the temperature in statistical physics and equivalent to scaling all model parameters as $\Omega \to \Omega/T$ (see Methods). Specifically, $C$ is given by $C(\Omega,T) = \text{var}[E]/T^2$. The “temperature” $T$ controls the level of disorder and its effect can be understood by examining the energy levels $\{E\}$ that are accessible to the system. At low temperatures the system is predominantly silent, it accesses few and separated excited states, it is relatively ordered (its entropy $S$ is low), and, since interactions dominate, the system scarcely fluctuates, leading to weak correlations (Fig 7D and 7E). At high temperatures the system has a high probability of occupying the excited states, which are barely separated, making it easy to fluctuate among them, thus increasing the disorder and decreasing the correlations (fluctuations dominate over interactions, effectively decorrelating the system). As for low temperatures, high temperatures lead to a low $C$. However, for a particular temperature, $T_{max}$, a large range of energies is accessible to the system, leading to a maximal $C$, and, moreover, the balance of fluctuations and interactions leads to a maximal mean correlation. This is the expected behavior close to a critical point where both order and
disorder coexist [59]. In conclusion, a maximal heat capacity close to the operating point (i.e., $T_{max}$ = 1, corresponding to the model of the real data) suggests that the system is likely to be close to a critical state.

We calculated the heat capacity $C$ as a function of $T$ for each cortical state (Fig 7F-7I). We found that the peak of $C$ is located at different temperatures for different cortical states: in all 4 cat datasets and in all 4 monkey datasets the temperature $T_{max}$ at which $C$ is maximal is significantly cortical-state-dependent ($p < 0.01$, one-way rm-ANOVA, Fig 7F and 7G) and is the highest for the most desynchronized cortical state and the lowest (approaching $T_{max} \approx 1$) for the most synchronized cortical state (Fig 7H and 7I).

We further studied the collective dynamics for different ensemble sizes $N$. For increasing $N$ the heat capacity of the data grows and the peak temperature decreases for all cortical states, except for the most desynchronized state for which $T_{max}$ increases (Fig 8A and 8B). Notably, the heat capacity of the data significantly grows even when it is normalized by $N$, as expected close to a critical point [59], for all but the desynchronized states (Fig 8B). Lastly, as a consequence of the relation between entropy and $C$ (see Methods), the entropy of the spiking patterns significantly depends on the cortical state ($p < 0.05$, rm-ANOVA) and it is higher in synchronized states (Fig 8C). The mean entropy differences between SynSlow II and Desyn I is $\Delta S = 0.428$ bits (average over all cat and monkey datasets, for $N = 6$); this is a substantial change, since this implies that the pattern repertoire is reduced $2^{NAS} \approx 6$-fold from the synchronized to desynchronized states. Altogether these results show that synchronized cortical states are closer to a maximum possible range of surprise (energies) and further suggest that synchronized cortical states are close to criticality.

<figure><img src="image_35.png" /><figcaption>**Fig 8.** Heat capacity and entropy as a function of ensemble size. (A) Peak temperature ($T_{max}$) averaged across anesthetized cat datasets for each cortical state and for ensembles of different sizes $N$. This analysis was restricted to cat datasets, for which more electrodes and longer recordings were available (datasets #1–4 for $N \le 12$ and datasets #1–3 for $N > 12$). Inset: heat capacity functions for Desyn I (gray scale) and SynSlow II (red scale) for an example cat dataset. (B) Magnitude of the heat capacity, $C$ (top), and specific heat capacity, $C/N$ (bottom), evaluated at $T = 1$ for different ensemble sizes $N$. For all tested $N$, the heat capacity was significantly cortical-state-dependent ($p < 0.05$, rm-ANOVA). (C) Entropy (top) and specific entropy (bottom) as a function of $N$. For all tested $N$, the entropy was significantly cortical-state-dependent ($p < 0.05$, rm-ANOVA). In both (B) and (C) lines indicate least-squares linear fit, for which the slopes are showed in the brackets and the asterisks indicate slopes significantly different from zero (*: $p < 0.05$, **: $p < 0.01$), and error bars indicate SEM.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1005543.g008
## Spiking model and sub-sampling analysis

In the following, we evaluated the impact of using a small sample of recorded neurons on our analyses. Sub-sampling has been shown to underestimate neuronal correlations [60] and can lead to a breakdown of power law relationships [21,28,34]. In contrast, absence of power laws may also be due to genuine deviations from critical dynamics as we show above. To test how sub-sampling affects the discrimination between cortical states and their critical properties, we modeled a realistic spiking neuron network (see Methods) which is capable of generating population dynamics that closely resemble those seen in physiological recordings [35,53,61,62]. Depending on the overall drive, the network can exhibit irregular bouts of synchronized bursts (synchronous irregular state [61], SI, low input) or largely desynchronized activity (asynchronous irregular state [61], AI, high input), akin to the physiological states shown in this study (Fig 9A and 9B). In accordance with the predictions from criticality and similar to the presented in vivo data, synchronized activity in the fully sampled network displays power laws in cluster size distributions with an exponent close to -1.5, power laws in lifetime distributions with exponents approaching -2 and CV values of inter-spike interval distributions >>1 (Fig 9C, 9E and 9G). In contrast, desynchronized states in the full model created clearly curved distributions with loglikelihood ratios being much more negative than in synchronized states (Fig 9D and 9F). In addition, CV values of the inter-spike interval distribution were close to 1, as expected from a Poisson process. Finally, we investigated whether the simulated neuronal avalanches produced self-similar dynamics (shape collapse) as predicted by criticality theory (S3G-S3J Fig). We found that, in the synchronized state, the spiking network produces avalanches with a temporal evolution that could be rescaled, while desynchronized model activity patterns did not collapse suggesting a departure from critical dynamics.

Next, we sub-sampled the model by randomly removing neurons (50 times) such that the cut off in power laws seen in the cluster size distributions resembled those seen in the cat (cut-off: ~60 spikes; population rate: 450 Hz) and monkey recording (cutoff: ~20 spikes; population rate: 150 Hz). In both cases, LLRs of cluster size and lifetime distributions as well as CV values of inter-spike interval distributions remained distinguishable from desynchronized activity sub-sampled with the same degree. With more severe sub-sampling (population rate: 50 spikes) cluster size and lifetime distributions could not be delineated by loglikelihood ratios anymore. Only CV values of interspike interval distributions remained sensitive enough to separate synchronized and desynchronized states at this level of sub-sampling. In conclusion, these results demonstrate that if synchronization levels of the population dynamics remain sufficiently strong, cluster size and lifetime distributions remain sensitive to differences in dynamical states despite heavy sub-sampling.

Sub-sampling could be also an issue when analyzing the spiking activity with the MEM. We thus studied the behavior of the heat capacity and entropy functions obtained for different ensemble sizes *N* for the spiking model activity in the subcritical and critical regimes. As in our previous analyses on the empirical data, we constructed MEMs of size *N* = 6–20 by randomly selecting *N* model neurons. The bin sizes were chosen to obtain the same averaged emission probability as in the empirical data (equal to 0.048). We next used the learned models to calculate the heat capacity curves given by estimating the variance of the energies (*E*) from Monte Carlo simulations (5.10⁵ steps) for different temperature parameters *T*. As shown in Fig 9I, the heat capacity behaved as expected for the two different network regimes, even for small neuronal ensembles of *N* = 6. The peak temperature *T*<sub>max</sub> was significantly different between the two dynamical regimes (p < 0.001, t-test), it was the lowest for the critical regime and approached 1 for increasing *N*. Furthermore, the normalized heat capacity (*C*/*N*) significantly increased with *N* for the critical regime (slope = 0.005, p < 0.001) but remained constant in
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

**Fig 9.** Sub-sampling analysis of a spiking neuronal model. (A-B) Model spike rasters of synchronized and desynchronized states. (C) Cluster size distributions of synchronized and desynchronized model activity. (D) Loglikelihood ratios of fitted power law and lognormal distributions fitted to the cluster size distributions for different levels of sub-sampling. (E) Lifetime distributions for synchronized and desynchronized model dynamics. (F) Same as in (D) for lifetime distributions. (G-H) Inter-spike interval distributions of population spike trains (ISI<sub>pop</sub>) and corresponding coefficients of variation for three different degrees of sub-
sampling. (I) Heat capacity as a function of the temperature parameter (T), for each dynamical regime of the spiking network (blue: critical, red: subcritical). Twenty random choices of groups of *N* model neurons were used; error bars indicate SEM. The insets show: the peak temperature (Tmax) as a function of *N*, the amplitude of the heat capacity normalized by *N*, and the entropy per neuron in both dynamical regimes. Lines indicate least-squares linear fits.

https://doi.org/10.1371/journal.pcbi.1005543.g009

the subcritical regime (slope = 0.002, p > 0.05). Finally, the resulting entropy was higher in the critical regime than in the subcritical regime (p < 0.01, t-test). Note also that the values of the heat capacity and the entropy were very similar to the ones obtained using the empirical data. We conclude that MEMs constructed using a limited number of neurons effectively discriminate between the different network regimes.

Lastly, we also investigated whether the statistical differences in our data are due to artificially dividing critical dynamics into different states. To this end we divided the simulated critical model dynamics into different synchronization regimes based on the Fano factor of the population spike train (see Supplementary Material), which shows differences across states in the experimental spiking data (S4A–S4C Fig). In contrast to the recordings, collapse and power laws, albeit with smaller cutoff sizes, were preserved during periods of reduced levels of population synchronization and firing rate in the model SI state (S4D–S4F Fig). These results are expected since the self-similarity in the model creates dynamics with statistics that are immune to changes of observation scale, in this case limiting analysis to less synchronized activity. They also indicate that the differences in power law statistics in the recordings reflect genuine state changes and are not due to artificial separation of a critical dynamical regime into different states.

## Discussion

We have shown that neuronal dynamics in the primary visual cortex during anesthesia and in the awake state transiently unfold in a critical state. In our study, the cortex underwent continuous state transitions with varying levels of LFP power in slow frequency bands and overall spike synchronization. Furthermore, the degree of synchronization was linked to the degree of criticality as measured by power laws and heat capacities calculated using collective spiking activity. Notably, the extent of critical dynamics could be predicted from the amount of fluctuations in local field potentials and the vigilance state of the animal, associating states of prolonged eye closure and putative drowsiness with higher levels of criticality.

These findings are in contrast to previously reported invariant presence [21] or absence of critical features in spiking activity [31,32,51]. In the latter studies, invariant deviations from power law were found across several species, brain areas and stages of the sleep wake cycle after averaging criticality measures across entire datasets. However, after having separated statistics for various cortical states, we have shown that the cortex is indeed capable of transiently generating dynamics that come close to criticality. Moreover, our analyses using the MEM confirm that cortical states have consistently different second-order statistics and demonstrate that cortical states are truly different physical states in the sense of statistical physics. Our results are in line with a previous study which reported fluctuations of criticality measures in spiking activity recorded in the barrel cortex of rats under urethane anesthesia [44]. We extended this finding to the primary visual cortex of anesthetized cat and awake monkey and showed in the latter condition that the vigilance state as measured by prolonged periods of eye closure can predict a switch to dynamics closer to criticality. It has been known for a long time that EEG synchronization in the visual cortex is stronger when eyes are closed as compared to when they are open [63]. Here, we link the “closed eye” state not only with higher synchronization in LFP and spiking activity, but also with dynamical changes towards a critical state.
Moreover, we found that the level of synchronization during isoflurane anesthesia in the cat can vary at very slow time scales which is in accordance with a recent report in opioid anesthetized monkeys [64].

Another recent study, however, reported conflicting results obtained with calcium imaging in the frontal cortex of head-fixed awake and isoflurane anesthetized rats [22]. Firing rates and correlations remained stable within each vigilance state, while power laws were only found during waking. In contrast, we found largely subcritical distributions during the awake state during the eyes open condition, which were replaced by more critical distributions during the closed eye condition. In addition, correlations during anesthesia were not stationary in our study despite having used the same anesthetic agent, similar to other studies [38,44,64]. These discrepancies may originate in different levels of anesthesia (0.8% isoflurane concentration in this study, 1-2% in [22] potentially influencing the appearance of state fluctuations, temporal resolution of spike recordings (fast electrophysiological recordings vs. slow calcium imaging) or recording site (primary visual cortex vs. frontal cortex). Further research needs to be undertaken to clarify these different results.

Our findings of a tight link between state changes are difficult to reconcile with the sub-sampling hypothesis of criticality which attributes subcritical distributions in spiking activity to the inevitable sub-sampling of neuronal discharges in empirical data and has been used to explain absence of criticality [20,21,28,34]. This susceptibility to sub-sampling has rendered spiking activity a less likely candidate for finding robust features of criticality [11], even though spiking activity is the backbone of neuronal avalanches in the brain. In contrast, theoretical studies have demonstrated that neural networks can display different collective behaviors, such as asynchronous irregular dynamics, in which excitatory avalanches are precluded by fast inhibition and strong external drive, or synchronized bursting activity, which can be accompanied by statistics expected from a critical state [32,53,62,65]. As predicted from these models, previous experimental work has indeed found similar dynamical transitions of spiking activity [37-41,66]. The spiking neuronal network model used in the current study was also able to reproduce different states with statistics similar to our *in vivo* results, showing synchronized states with signs of criticality and subcritical desynchronized states as a function of external drive. Importantly, differences between distinct states in the model were detectable in criticality statistics even with strong sub-sampling that mimicked the cutoff of power laws seen in the monkey and cat recordings, providing support that the sampling in our experiment was sufficient to observe deviations from criticality in our data. Moreover, the MEM model revealed robust differences in criticality between states requiring only a small sample of neurons (~10). In conclusion, we provide evidence that criticality statistics vary despite invariant spatial sampling of neuronal discharges and can show near critical features, when the cortex is engaged in synchronized population dynamics.

The existence of different cortical states and corresponding levels of criticality raise the question of their functional significance. Traditionally, the larger number of possible neuronal patterns provided by the complex dynamics in a critical state has been suggested to optimize a number of cortical functions for which some evidence has been found experimentally [9,67].

However, recent experimental data suggest that active processing may actively desynchronize spiking activity in the cortex [35,39-41,43,68] through neuromodulatory or glutamatergic signals [35,69], in a way similar to the desynchronized and subcritical dynamics described in this study. This desynchronization reduces response variability and enables a reliable representation of stimuli, in contrast to synchronized states [70-72]. In addition, desynchronized activity may be necessary to establish precise and unambiguous communication between areas [73], possibly embedded within gamma oscillations [74,75] that are more prevalent during cortical desynchronization [76,77]. Consistent with this line of reasoning, task related focused
attention was found to trigger subcritical dynamics in EEG recordings of large cortical networks of humans [50]. The same study, however, found near critical dynamics during resting state arguing in favor of criticality when cognitive demands require exploration of different network states. Similarly, we observed that collective dynamics during synchronized states approach the maximal range of surprise (heat capacity), indicating that the cortical system might be spontaneously exploring its dynamical repertoire (i.e., the system visits a set of configurations, or patterns, more diverse in terms of associated energies; Fig 7). Another hypothesis considers desynchronized states as quasi-critical, maintaining optimization of computational functions in the absence of strict criticality [78]. Finally, our results are consistent with recent studies showing that hierarchical connectivity of neural networks can enlarge the region where critical-like behavior is observed, the so-called Griffith phase [79–81]. In this phase, rare occurrence of local order can be observed even though the system is globally in the disordered phase, thus making it a candidate mechanism for the observed transitions between cortical states. Further studies are required to clarify the precise cognitive role of critical dynamics in the brain.

In conclusion, we have shown that the cortex can transiently approach criticality during anesthesia and awake states in the primary visual cortex. These fluctuations of criticality can be predicted from more global variables such as the LFP and also reflect different vigilance states. Our results suggest that the cortex evolved mechanisms to synchronize and desynchronize its activity according to computational needs, thereby continuously switching between critical and more subcritical dynamics.

# Materials and methods

## Preparation

Recordings were obtained from four anesthetized (isoflurane) and paralyzed cats (three females, one male), and one awake monkey (male). The isoflurane concentration was kept at a constant value of ~0.8%. Experimental protocols on behaving monkeys have been approved by the Marseille Ethical Committee in Neuroscience (approval #A10/01/13, official national registration #71-French Ministry of Research). The animals used in the experiments were bred in the Central CNRS Animal Care (French Agriculture Ministry Authorization: B91-272-105) under required veterinary and National Ethical Committee supervision.

## Recording

We inserted 32 channel silicon-based micro-electrode arrays (four shanks with eight electrode contacts each, distance between electrode contacts: 400 µm, electrode contact impedance: 0.3–0.5 MΩ at 1000 Hz, shank length: 3mm, Neuronexus Technologies, Ann Arbor, USA) into area 17 of four cats (one array per cat) and chronically implanted a Utah array (Blackrock Microsystems, Salt Lake City, UT, USA) into the primary visual cortex (near-foveal retinotopic region) of an adult macaque monkey (macaca mulatta—10kg). The Utah array was composed of 96, 1 mm long, electrodes arranged in a 10 x 10 matrix with an inter-electrode distance of 400 microns. In both preparations, continuous spontaneous multiunit activity (sampling frequency: 30 kHz,) and local field potentials (LFP, sampling frequency: 1 kHz) were recorded with a Cerebus acquisition system (Blackrock Microsystems, Salt Lake City, UT, USA). For each cat we acquired one dataset lasting between 3900 and 6000 seconds, while four datasets of 600 seconds each were obtained in the monkey on four different days. The monkey was sitting in the complete dark, head fixated. For subsequent analysis the LFP was filtered between 1 and 100 Hz. Spikes were detected by manually setting a one-sided threshold that would result in a signal to noise ratio larger than ~2. Action potential waveforms were sorted and separated
from remaining noise using the T-EM clustering algorithm and manual cluster cutting (Off-line Sorter, Plexon Inc, Dallas, USA). In the cat, simultaneous recordings of spiking activity were extracted from an average of 21.5±6.7 channel sites and in the monkey 8.75±1.9 channels sites. Note that no attention was paid to extract single units, as during later analysis spikes from each electrode were pooled together into one population spike train. The mean firing rate per site ranged between 0.61 Hz and 1.44 Hz (1.01±0.38) in the cat and between 4.95 Hz to 7.32 Hz (6.31±1.08) in the monkey recordings. The rates for the entire recordings were 20.64±7.9 Hz for the cat data and 54.24±10.57 Hz for the monkey. The cat recordings were performed in uniform scotopic illumination level, while any source of light was eliminated during the monkey recordings. During three recording sessions in the monkey we obtained eye tracking data using the Eyelink 1000 system (SR research, sampling frequency 1000Hz).

## Separation of cortical states

Cortical states and degrees of synchronization of neuronal activity can switch at a time scale of seconds [35]. In order to capture this fast dynamics a moving temporal window of one second was chosen to analyze different synchronization levels in cortical activity. We based the separation of cortical states on differences in the frequency composition of the LFP between the one second segments. To this end we computed power spectra for each channel of a given segment with the multitaper method (see http://chronux.org) and averaged across all channels to obtain one spectrum per segment. Next, we split the power spectrum into 1 Hz frequency bins between 1 and 100 Hz and fed these 100 values as variables into a principal component analysis (PCA). As a result, we reduced the power spectrum of each segment to its first three principal components and each segment was represented by its position in a three dimensional PCA state space. In order to find segments with similar power spectra, we applied a k-means clustering algorithm with a different number of clusters (2–10 clusters) after normalization of each principal component to a value of 1. To determine the optimal number of clusters we validated the clustering results with the Dunn index [82] according to the following formula:

$$DI = \min_{1 \le i \le n} \left\{ \min_{1 \le j \le n, i \ne j} \left\{ \frac{d(i, j)}{\max_{1 < k < n} d'(k)} \right\} \right\},$$

with distance $d(i, j)$ between the centers of clusters $i$ and $j$ and distance $d'(k)$ between the center and elements of cluster $k$. This equation evaluates the compactness of clusters by calculating the ratio between the minimal distance between clusters to the maximal distance within clusters. The Dunn index was computed for each dataset and then averaged separately across all cat datasets and all monkey datasets. The Dunn index reaches a maximum at the optimal cluster number and we identified five clusters for each cat dataset and four clusters for each monkey dataset (Fig 1C and 1D). The clusters of each recording were consistently color coded according to their position in the state space, such that the clusters were comparable across datasets and species. The different clusters were associated with different cortical states. All the segments of a cluster were concatenated and used as separate datasets in subsequent analysis.

## Correlation analysis

To investigate the correlation characteristics of the recorded neuronal population, we collated together the spiking activity of all channels into one population spike train. Correlations were assessed by computing the auto-correlation histograms (ACH). An ACH was calculated for each 1 second segment and subsequently averaged across all segments of a cortical state. We also subtracted a shift predictor that was calculated from the average of 100 randomly shuffled
surrogates with the same firing rate to normalize the baseline for each segment. All ACHs were normalized to a center peak of one and quantified based on the peak amplitude (i.e. after removal of the peak at zero time lag) and the integral between the ACH curve and the baseline for time lags between -250 and 250ms. In order to allow a comparison between different datasets and species, these measures were normalized such that the cortical state with the maximum value for each dataset was set arbitrarily to one.

We also estimated the correlation between spiking activity and the LFP by constructing spike triggered averages (STA) for each cortical state of a given dataset. The STAs of all clusters and datasets showed a negative deflection. LFP was z-scored and STAs were quantified by calculating the amplitude of the deflection peak and the integral of the negative deflection with respect to the zero crossings. Like for the ACH, these two measures were normalized to the cortical state with the maximum value which was set to one.

## Neuronal avalanche analysis

The spiking data were binned with bin-size $\Delta t$ and event clusters were defined as groups of consecutive bins each containing a number of spikes above a predefined threshold. This threshold was varied between 1–5 spikes in the cat and 1–3 spikes in the monkey. In both cases the upper bound was given by the threshold that yielded enough avalanches to fit distributions (see below). The first and last avalanches of an analyzed time window were discarded, if they were not preceded or followed by a bin with a number of spikes below the threshold, respectively. The cluster size was equal to the number of events within a cluster. In accordance with previous studies [13,21], $\Delta t$ was chosen as the average inter-spike interval of the population spike trains (also denoted as ISI<sub>pop</sub>, Table 1) which reflects a compromise between spurious concatenation of small clusters (large $\Delta t$) and separation of larger clusters (small $\Delta t$). The lifetime of an avalanche was defined as the number of bins spanning an avalanche. Distributions of cluster sizes and lifetimes were plotted in log-log coordinates and further analyzed.

In order to statistically characterize the distribution of cluster sizes $s$, we first fitted a power law distribution $P(s, \alpha) \sim C s^{\alpha}$, with power law exponent $\alpha$. Since on visual inspection size distributions in desynchronized states clearly deviated from a power law, we also fitted lognormal distributions with the following probability density function:

$$N(s, \mu, \sigma) = \frac{1}{s\sigma\sqrt{2\pi}} \exp\left(-\frac{(\ln s - \mu)^2}{2\sigma^2}\right); \sigma > 0, \mu \ge 0 \text{ with scale parameter } \sigma \text{ and location parameter } \mu. \text{ Depending on } \sigma, \text{ this distribution can assume heavy tails close to a power law or lighter tails. We then compared the fits to both distribution by calculating the loglikelihood ratio and a corresponding p-value determining the significance of model fit differences [22,52]. In our study, positive LLR values indicated superior power law fits, while negative values favored lognormal distributions. All fits were performed using the Python toolbox for analysis of heavy-tailed distributions [83]. Parameters were calculated for the entire distribution. The same fitting procedures were applied to lifetime distributions. Maximum entropy models (MEMs). We further studied the patterns of the ensemble spiking activity from $N$ recording sites during the different cortical states in the framework of statistical mechanics. The ensemble activity was binarized in non-overlapping time bins of 50 ms for anaesthetized cat data, and 10 or 20 ms for awake monkey data (different time bins were chosen to compensate for the different averaged activation rates between datasets and cortical states; see Table 2 and Fig 7A). In a time bin $\Delta t$, a single electrode site $i$ either did ($\sigma_i = +1$) or did not ($\sigma_i = -1$) generate one or more spikes, thus, the state of the neural ensemble is described by a binary pattern $\vec{\sigma} = [\sigma_1, \sigma_2, \dots, \sigma_N]$. We used a Maximum entropy model (MEM) to estimate the distribution $P(\vec{\sigma})$ of each of the $2^N$ possible patterns, i.e., we estimated $P(\vec{\sigma})$ by maximizing its entropy under the constraint that some empirical statistics are
preserved. A pairwise-MEM provides a solution under the constraint that the activation rates (<σᵢ>) and the pairwise correlations (<σᵢσⱼ>) are preserved. It is known that the maximum entropy distribution P that is consistent with these expectation values is given by [59,84]:

$$P(\vec{\sigma}) = \frac{e^{-E(\vec{\sigma})}}{Z},$$

where the normalization factor Z = Σ<sup>2N</sup><sub>s=1</sub> exp [E(σ<sub>s</sub>)] is the partition function and E(σ<sub>s</sub>) is the energy of the pattern s, with s ∈ {1, ..., 2<sup>N</sup>}, and E is given by:

$$E(\vec{\sigma}) = - \sum_{i=1}^{N} h_i \sigma_i - \frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} J_{ij} \sigma_i \sigma_j,$$

where, hᵢ represents the intrinsic tendency of site i towards activation (σᵢ = +1) or silence (σᵢ = -1) and Jᵢⱼ represents the effective interaction between sites i and j. Note that the energy represents the patterns' minus log probabilities, log P, or surprise, plus the constant log Z. The estimation of the model parameters **h** and **J** was achieved through a gradient descent algorithm (see below). Once we know the parameters Ω = {**h**; **J**} the expected probability of any pattern is known, by first calculating the energy associated to the pattern and then computing P. For each recording session, we constructed models for each cortical state separately by, first, randomly selecting N signals, second, concatenating the subset of all binarized patterns belonging to the same cortical state, and, then, running the parameter learning procedure for each subset of patterns (Table 2 summarizes the characteristics of the data). We did this for a total of Q = 10 random choices of ensembles of N signals.

**Estimation of MEM parameters.** The MEM parameters hᵢ and Jᵢⱼ were iteratively adjusted to minimize the absolute difference between the empirical activation rates (<σᵢ>) and correlations (<σᵢσⱼ>) and those (<σᵢ>model, <σᵢσⱼ>model) predicted by the model through Monte Carlo simulations. Specifically, each iteration is given by: h<sup>new</sup><sub>i</sub> = h<sup>old</sup><sub>i</sub> - α(<σᵢ>model - <σᵢ>), and J<sup>new</sup><sub>ij</sub> = J<sup>old</sup><sub>ij</sub> - α(<σᵢσⱼ>model - <σᵢσⱼ>), where α is the learning rate (α = 0.1). In our study, we stopped the re-estimations once the differences between the empirical and model values are less than a tolerance threshold (0.005) or if this tolerance was not reached within a maximum number of iterations (100).

**MEM goodness-of-fit.** The goodness-of-fit of the MEMs was evaluated using the Jensen-Shannon divergence (D<sub>JS</sub>) between the probability distribution of the empirical and model binary patterns [85]. D<sub>JS</sub> is a symmetric version of the Kullback-Leibler divergence (D<sub>KL</sub>) and is given as:

$$D_{JS}(P_{\text{emp}}|P_{\text{model}}) = \frac{1}{2}D_{KL}(P_{\text{emp}}|(P_{\text{model}} + P_{\text{emp}})/2) + \frac{1}{2}D_{KL}(P_{\text{model}}|(P_{\text{model}} + P_{\text{emp}})/2).$$

**MEM decoding.** We classified the cortical state from the spiking ensemble activity using a jackknife cross-validation approach. We used a total of Q = 10 random choices of ensembles of N signals. For each cortical state, 70% of the data was used to learn the parameters of a MEM (train-sets). The remaining data (test-sets) was associated to a cortical state by choosing the model with maximal goodness-of-fit (1/D<sub>JS</sub>) between the probability distribution of the test data and the learned model. The percentage of correct classifications was computed across all neuronal ensembles. The entire procedure was repeated 30 times with randomly selected train-sets to get confidence intervals of the classification performance. To assess statistical significance of the classification performance we calculated the probability of getting k correct classifications (hits) by chance, which is given by: Pr(k) = C<sup>k</sup><sub>n</sub>p<sup>k</sup>(1 - p)<sup>n-k</sup>, where p is the
probability of getting a hit by chance ($p = 1/N_s$, where $N_s$ is the number of cortical states) and $n$ is the number of tests ($n = QN_s$). Significant decoding is reached when the decoding performance exceeds the 95% percentile of $Pr(k)$.

**Parameter re-scaling, heat capacity, and criticality.** In the framework of statistical mechanics, one can characterize the physical state of the ensemble activity through the heat capacity. The heat capacity $C$ is given by $C = \text{var}[E]/T^2$ and it is calculated by introducing a “temperature” parameter $T$ that acts as a scaling factor for all model parameters as $\Omega \to \Omega/T$. In our study, we estimated the variance of the energies ($E$) from a large number of Monte Carlo simulations (5.10⁵ steps) for different $T$. A useful relation is the one linking the entropy and the heat capacity: the entropy $S$ is the integral of the function $C(T)/T$ from $T=0$ to $T=1$ [59]. As shown above (Fig 7D and 7E), a maximum of the heat capacity close to $T=1$ suggests that the system is likely to be close to a critical state. Hence, $C$ can be used as an additional diagnosis tool to assess criticality.

**Neuronal network model.** To test the effect of sub-sampling on the assessment of criticality, we modeled a network of 2500 integrate and fire neurons (70% excitatory, 30% inhibitory) which were arranged on a two-dimensional grid (50x50). Every neuron was assigned a spatial position (x,y) and connected to neighboring with connection probability ($p$) which decayed with distance according to Gaussian distribution: $p = \exp(-\frac{(x_i-x_j)^2+(y_i-y_j)^2}{2\tau_w^2})$, with spatial decay constant $\tau_w = 100$ micrometers and a distance to the nearest neighbor of 25 microns [53,86].

The membrane potential $v_m$ of each neuron evolves according to the following equation:

$$ \frac{dv_m}{dt} = \frac{v_{rest} - v_m}{\tau_m} + \frac{g_{exc}(v_{exc} - v_m) + g_{inh}(v_{inh} - v_m)}{C}, $$

where $v_{rest}$ is the resting membrane potential with $v_{exc}$ and $v_{inh}$ representing the excitatory and inhibitory reversal potentials. The parameters $g_{exc}$ and $g_{inh}$ reflect the excitatory and inhibitory conductances, respectively and decay according to $\frac{dg_{exc}}{dt} = -g_{exc} \frac{1}{\tau_{exc}}$ and $\frac{dg_{inh}}{dt} = -g_{inh} \frac{1}{\tau_{inh}}$ with time constants $\tau_{exc}$ and $\tau_{inh}$ (see Table 3 for detailed values of neuronal parameters used in the simulation). The values of excitation and inhibition were chosen such that synchronized network activity emerged with power law statistics. A spike was emitted when the membrane potential reached a fixed threshold $v_{th}$ of -50 mV and was reset to the resting membrane potential (-60 mV). Each neuronal discharge was followed by an absolute refractory period of 1ms. After a delay of 1ms, a spike triggered excitatory or inhibitory postsynaptic potentials in the receiving neurons, which decayed according to an exponential function.

<table><caption>**Table 3. Neuronal parameters used for the spiking neuronal network model.**</caption><thead><tr><td>Neuronal parameter</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>C</td><td>200 pF</td><td>Membrane capacitance</td></tr><tr><td>τ<sub>m</sub></td><td>20 ms</td><td>Membrane time constant</td></tr><tr><td>V<sub>rest</sub></td><td>-60 mV</td><td>Resting membrane potential</td></tr><tr><td>V<sub>th</sub></td><td>-50 mV</td><td>Firing threshold</td></tr><tr><td>τ<sub>exc</sub></td><td>5 ms</td><td>Time constant excitatory synapses</td></tr><tr><td>τ<sub>inh</sub></td><td>10ms</td><td>Time constant inhibitory synapses</td></tr><tr><td>g<sub>exc</sub></td><td>0.7 nS</td><td>Excitatory peak conductance</td></tr><tr><td>g<sub>inh</sub></td><td>50 nS</td><td>Inhibitory peak conductance</td></tr><tr><td>V<sub>exc</sub></td><td>0 mV</td><td>Reversal potential of excitatory synapses</td></tr><tr><td>V<sub>inh</sub></td><td>-70 mV</td><td>Reversal potential of inhibitory synapses</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1005543.t003
Each neuron received a constant external excitatory Poisson input $I_{ext}$, adding the value $\frac{g_{exc}(v_{exc} - v_{m})}{C}$ to the membrane potential for each Poisson spike. The rate of $I_{ext}$ controls the state of the network: for low input ($I_{ext} = 6,000$ Hz) the network displayed patterned synchronized activity, for high input ($I_{ext} = 9,000$ Hz) the network was desynchronized. For each value of the external input, we performed two long simulations of 10 minutes. Thereafter, neurons were randomly removed from each simulation to test the consequences of sub-sampling for various network dynamic statistics. Each simulation was sub-sampled such that the population firing rate was reduced to 450Hz, 150Hz and 50Hz. For each condition we sub-sampled the network 50 times and subsequently averaged statistics were calculated. Standard deviations between the individual sub-sampling trials were exceedingly small, such that they were not plotted in Fig 6D, 6F and 6H. All simulations were performed using the Brian2 simulation environment [87].

**Statistical analysis.** Statistical differences between different cortical states were assessed using one-way repeated measures (rm) ANOVA followed by multiple comparisons using Bonferroni tests. The threshold for statistical significance was set to p-values<0.05. Where the ANOVA's sphericity assumption was not met (using the Mauchly test), p-values and degrees of freedom were corrected using the Huynh-Feldt estimates of sphericity ($\epsilon$). The confidence intervals for the slope of least-squares linear fit was calculated as: $t_{n-2}(i) \sqrt{S^2/S_{xx}}$, where n is the number of data points, $S^2$ is the unbiased estimator of the variance of the estimation errors, $S_{xx}$ is the sum of squared deviations of observed predictors from their mean, $t_{n-2}(i)$ is the i-th percentile for the Student t distribution with (n-2) degrees of freedom (we used 95% and 99%).

**Supporting information**

**S1 Fig.** Loglikelihood ratios for power law and lognormal fits to size distributions of different cortical states across all cat and monkey data and for different thresholds of spike cluster definition. Negative values indicate a better lognormal fit. State differences were assessed using a one-way rm-ANOVA test (threshold 1: cat: $F_{4,12} = 4.03$, p = 0.1, $\epsilon = 0.38$; monkey: $F_{3,9} = 24.57$, p = 0.002, $\epsilon = 0.61$; threshold 2: cat: $F_{4,12} = 5.6$, p = 0.03, $\epsilon = 0.6$; monkey: $F_{3,9} = 25.31$, p = 0.0005, $\epsilon = 0.78$; threshold 3: cat: $F_{4,12} = 9.48$, p = 0.005, $\epsilon = 1$; monkey: $F_{3,9} = 20.89$, p = 0.001, $\epsilon = 0.73$; threshold 4: cat: $F_{4,12} = 3.68$, p = 0.029, $\epsilon = 1$; threshold 5: cat: $F_{4,12} = 1.9$, p = 0.12, $\epsilon = 1$).

(TIF)

**S2 Fig.** Loglikelihood ratios for power law and lognormal fits to lifetime distributions of different cortical states across all cat and monkey data and for different thresholds of spike cluster definition. Negative values indicate a better lognormal fit. State differences were assessed using a one-way rm-ANOVA test (threshold 1: cat: $F_{4,12} = 7.45$, p = 0.02, $\epsilon = 0.6$; monkey: $F_{3,9} = 25.6$, p = 0.002, $\epsilon = 0.58$; threshold 2: cat: $F_{4,12} = 5.21$, p = 0.05, $\epsilon = 0.46$; monkey: $F_{3,9} = 27.02$, p = 0.0006, $\epsilon = 0.74$; threshold 3: cat: $F_{4,12} = 12.2$, p = 0.00003, $\epsilon = 1$; monkey: $F_{3,9} = 22.34$, p = 0.002, $\epsilon = 0.63$; threshold 4: cat: $F_{4,12} = 2.03$, p = 0.2, $\epsilon = 0.58$; threshold 5: cat: $F_{4,12} = 1.26$, p = 0.34, $\epsilon = 0.98$).

(TIF)

**S3 Fig.** Avalanche shape collapse. (A–B) Averaged temporal profile of avalanche of lifetime $\Delta t$, i.e., $\langle S(t,\Delta t) \rangle$, in the Desyn I cortical state (A) and the SynSlow cortical state (B) for an example cat dataset. (C–D) Scaled avalanche profiles as a function of the scaled time $t/\Delta t$, in the Desyn I cortical state (C) and the SynSlow cortical state (D). Red line: averaged scaled avalanche profile; $\sigma v_z$: best scaling parameter. (E) Collapse index (CI) for each cortical state,
averaged over all cat datasets ($F_{4,12} = 1.53$, $p = 0.254$; $\epsilon = 1$). Error bars indicate SEM. (F) CI calculated by grouping the avalanches of the desynchronized cortical states (Desyn I/II) and, separately, those of the synchronized states (SynSlow I/II) ($p = 0.027$, paired t-test). Error bars indicate SEM. (G-H) Averaged temporal profile of avalanche of lifetime $\Delta t$ of avalanches display by the spiking model in the desynchronized states (G) and in the synchronized state (H). (I-J) Scaled avalanche profiles for the spiking model in the desynchronized state (I) and the synchronized state (J).

(TIF)

### S4 Fig. Separation of different states within the synchronized state of the model. (A-B)

Mean Fano factor (FF) for different states in cat and monkey recordings (bin-size = 100ms). Error bars indicate SEM. (C) Distribution of FFs computed for each one second segment of the modeled synchronized state (bin-size = 50ms). (D) Mean firing rate of model neurons for different parts of the FF distribution. (E) Cluster size distributions of different levels of synchronization in the SI state defined by the FF distribution.

(TIF)

### S1 Dataset. Spiking data from four cat and four monkey datasets. The data include the spike timings, electrode index and state information. Details can be found in the README file.

(ZIP)

### S1 Text. Shape collapse and state analysis of the synchronized model state.

(DOCX)

## Acknowledgments

We thank Aurélie Daret (UNIC-CNRS) for help with the preparation of the cat experiments and Etienne Hugues for helpful discussions.

## Author Contributions

**Conceptualization:** GH APA.

**Formal analysis:** GH APA.

**Investigation:** GH CM GB.

**Methodology:** GH APA.

**Supervision:** AK FC GD YF.

**Writing - original draft:** GH APA CM GB AK FC GD YF.

## References

1. Fox MD, Raichle ME. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat Rev Neurosci. 2007; 8: 700–11. https://doi.org/10.1038/nrn2201 PMID: 17704812

2. Petersen SE, Posner MI. The attention system of the human brain: 20 years after. Annu Rev Neurosci. 2012; 35: 73–89. https://doi.org/10.1146/annurev-neuro-062111-150525 PMID: 22524787

3. Alkire MT, Hudetz AG, Tononi G. Consciousness and Anesthesia. Science (80-). 2008; 7: 876–880. https://doi.org/10.1126/science.1149213 PMID: 18988836

4. Takahashi K, Kim S, Coleman TP, Brown KA, Suminski AJ, Best MD, et al. Large-scale spatiotemporal spike patterning consistent with wave propagation in motor cortex. Nat Commun. 2015; 6: 7169. https://doi.org/10.1038/ncomms8169 PMID: 25994554
5. Kenet T, Bibitchkov D, Tsodyks M, Grinvald A, Arieli A. Spontaneously emerging cortical representations of visual attributes. Nature. 2003; 425: 954–956. https://doi.org/10.1038/nature02078 PMID: 14586468

6. Tsodyks M, Kenet T, Grinvald A, Arieli A. Linking spontaneous activity of single cortical neurons and the underlying functional architecture. Science. 1999; 286: 1943–6. Available: http://www.ncbi.nlm.nih.gov/pubmed/10583955

7. Brookes MJ, Woolrich M, Luckhoo H, Price D, Hale JR, Stephenson MC, et al. Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc Natl Acad Sci U S A. 2011; 108: 16783–8. https://doi.org/10.1073/pnas.1112685108 PMID: 21930901

8. Adrian ED, Matthews BH. The interpretation of potential waves in the cortex. J Physiol. 1934; 81: 440–71. Available: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1394145&tool=pmcentrez&rendertype=abstract PMID: 16994555

9. Shew WL, Plenz D. The functional benefits of criticality in the cortex. Neuroscientist. 2013; 19: 88–100. https://doi.org/10.1177/1073858412445487 PMID: 22627091

10. Bak P, Tang C, Wiesenfeld K. Self-Organized Criticality: An Explanation of 1/f Noise. Phys Rev Lett. APS; 1987; 59: 381–384. Available: http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle: Self-Organized+Criticality:+An+Explanation+of+1/f+Noise#8

11. Plenz D, Thiagarajan T. The organizing principles of neuronal avalanches: cell assemblies in the cortex? Trends Neurosci. 2007; 30: 101–10. https://doi.org/10.1016/j.tins.2007.01.005 PMID: 17275102

12. Chialvo D. Emergent complex neural dynamics: the brain at the edge. Nat Phys. 2010; 6: 744–750.

13. Beggs J, Plenz D. Neuronal avalanches in neocortical circuits. J Neurosci. 2003; 23: 11167–77. Available: http://www.ncbi.nlm.nih.gov/pubmed/14657176 PMID: 14657176

14. Beggs J, Plenz D. Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J Neurosci. 2004; 24: 5216–29. https://doi.org/10.1523/JNEUROSCI.0540-04.2004 PMID: 15175392

15. Friedman N, Ito S, Brinkman B, Shimono M, Deville R, Dahmen K, et al. Universal critical dynamics in high resolution neuronal avalanche data. Phys Rev Lett. 2012; 108: 208102. Available: http://www.ncbi.nlm.nih.gov/pubmed/23003192 https://doi.org/10.1103/PhysRevLett.108.208102

16. Mazzoni A, Broccard F, Garcia-Perez E, Bonifazi P, Ruaro M, Torre V. On the dynamics of the spontaneous activity in neuronal networks. PLoS One. 2007; 2: e439. https://doi.org/10.1371/journal.pone.0000439 PMID: 17502919

17. Pasquale V, Massobrio P, Bologna L, Chiappalone M, Martinoia S. Self-organization and neuronal avalanches in networks of dissociated cortical neurons. Neuroscience. 2008; 153: 1354–69. https://doi.org/10.1016/j.neuroscience.2008.03.050 PMID: 18448256

18. Lombardi F, Herrmann HJ, Plenz D, de Arcangelis L. Temporal correlations in neuronal avalanche occurrence. Sci Rep. 2016; 6: 24690. https://doi.org/10.1038/srep24690 PMID: 27094323

19. Gireesh E, Plenz D. Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3. Proc Natl Acad Sci U S A. 2008; 105: 7576–81. https://doi.org/10.1073/pnas.0800537105 PMID: 18499802

20. Hahn G, Petermann T, Havenith M, Yu S, Singer W, Plenz D, et al. Neuronal avalanches in spontaneous activity in vivo. J Neurophysiol. Am Physiological Soc; 2010; 104: 3312. https://doi.org/10.1152/jn.00953.2009 PMID: 20631221

21. Ribeiro T, Copelli M, Caixeta F, Belchior H, Chialvo D, Nicolelis M, et al. Spike Avalanches Exhibit Universal Dynamics across the Sleep-Wake Cycle. Sporns O, editor. PLoS One. 2010; 5: e14129. https://doi.org/10.1371/journal.pone.0014129 PMID: 21152422

22. Bellay T, Klaus A, Seshadri S, Plenz D. Irregular spiking of pyramidal neurons organizes as scale-invariant neuronal avalanches in the awake state. Elife. 2015; 4: 1–25. https://doi.org/10.7554/eLife.07224 PMID: 26151674

23. Petermann T, Thiagarajan T, Lebedev M, Nicolelis M, Chialvo D, Plenz D. Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc Natl Acad Sci U S A. 2009; 106: 15921–6. https://doi.org/10.1073/pnas.0904089106 PMID: 19717463

24. Yu S, Yang H, Shriki O, Plenz D. Universal organization of resting brain activity at the thermodynamic critical point. Front Syst Neurosci. 2013; 7: 42. https://doi.org/10.3389/fnsys.2013.00042 PMID: 23986660

25. Yu S, Yang H, Nakahara H, Santos G, Nikolić D, Plenz D. Higher-order interactions characterized in cortical activity. J Neurosci. 2011; 31: 17514–26. https://doi.org/10.1523/JNEUROSCI.3127-11.2011 PMID: 22131413
26. Shriki O, Alstott J, Carver F, Holroyd T, Henson RNA, Smith ML, et al. Neuronal avalanches in the resting MEG of the human brain. J Neurosci. 2013; 33: 7079–90. https://doi.org/10.1523/JNEUROSCI.4286-12.2013 PMID: 23595765

27. Tagliazucchi E, Balenzuela P, Fraiman D, Chialvo D. Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front fractal Physiol. 2012; 3: 15. https://doi.org/10.3389/fphys.2012.00015 PMID: 22347863

28. Priesemann V, Munk M, Wibral M. Subsampling effects in neuronal avalanche distributions recorded in vivo. BMC Neurosci. 2009; 10: 40. https://doi.org/10.1186/1471-2202-10-40 PMID: 19400967

29. Arviv O, Goldstein A, Shriki O. Near-Critical Dynamics in Stimulus-Evoked Activity of the Human Brain and Its Relation to Spontaneous Resting-State Activity. J Neurosci. 2015; 35: 13927–13942. https://doi.org/10.1523/JNEUROSCI.0477-15.2015 PMID: 26468194

30. Shew W, Clawson W, Pobst J, Karimipanah Y. Adaptation to sensory input tunes visual cortex to criticality. Nat Phys. 2015; Available: http://www.nature.com/nphys/journal/v11/n8/abs/nphys3370.html

31. Bédard C, Kroeger H, Destexhe A. Does the 1/f frequency scaling of brain signals reflect self-organized critical states? Phys Rev Lett. APS; 2006; 97: 118102. https://doi.org/10.1103/PhysRevLett.97.118102 PMID: 17025932

32. Priesemann V, Wibral M, Valderrama M, Pröpper R, Le Van Quyen M, Geisel T, et al. Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Front Syst Neurosci. 2014; 8: 108. https://doi.org/10.3389/fnsys.2014.00108 PMID: 25009473

33. Touboul J, Destexhe A. Can power-law scaling and neuronal avalanches arise from stochastic dynamics? PLoS One. 2010; 5: e8982. https://doi.org/10.1371/journal.pone.0008982 PMID: 20161798

34. Ribeiro TL, Ribeiro S, Belchior H, Caixeta F, Copelli M. Undersampled critical branching processes on small-world and random networks fail to reproduce the statistics of spike avalanches. PLoS One. 2014; 9. https://doi.org/10.1371/journal.pone.0094992 PMID: 24751599

35. Harris K, Thiele A. Cortical state and attention. Nat Rev Neurosci. Nature Publishing Group; 2011, 12: 509–23. https://doi.org/10.1038/nrm3084 PMID: 21829219

36. Destexhe A, Hughes S, Rudolph M, Crunelli V. Are corticothalamic “up” states fragments of wakefulness? Trends Neurosci. 2007; 30: 334–42. https://doi.org/10.1016/j.tins.2007.04.006 PMID: 17481741

37. Luczak A, Barthó P, Marguet S, Buzsáki G, Harris K. Sequential structure of neocortical spontaneous activity in vivo. Proc Natl Acad Sci. National Acad Sciences; 2007; 104: 347. Available: http://www.pnas.org/content/104/1/347.short

38. Renart A, De la Rocha J, Bartho P, Hollender L, Parga N, Reyes A, et al. The asynchronous state in cortical circuits. Science (80-). American Association for the Advancement of Science; 2010; 327: 587. Available: http://www.sciencemag.org/content/327/5965/587.short

39. Crochet S, Petersen C. Correlating whisker behavior with membrane potential in barrel cortex of awake mice. Nat Neurosci. 2006; 9: 608–10. https://doi.org/10.1038/nn1690 PMID: 16617340

40. Poulet J, Petersen C. Internal brain state regulates membrane potential synchrony in barrel cortex of behaving mice. Nature. 2008; 454: 881–5. https://doi.org/10.1038/nature07150 PMID: 18633351

41. Okun M, Naim A, Lampl I. The subthreshold relation between cortical local field potential and neuronal firing unveiled by intracellular recordings in awake rats. J Neurosci. 2010; 30: 4440–8. https://doi.org/10.1523/JNEUROSCI.5062-09.2010 PMID: 20335480

42. Poulet J, Fernandez L, Crochet S, Petersen C. Thalamic control of cortical states. Nat Neurosci. 2012; 15: 370–2. https://doi.org/10.1038/nn.3035 PMID: 22267163

43. Tan AYY, Chen Y, Scholl B, Seidemann E, Priebe NJ. Sensory stimulation shifts visual cortex from synchronous to asynchronous states. Nature. Nature Publishing Group; 2014; 509: 226–9. https://doi.org/10.1038/nature13159 PMID: 24695217

44. Gautam SH, Hoang TT, McClanahan K, Grady SK, Shew WL. Maximizing Sensory Dynamic Range by Tuning the Cortical State to Criticality. PLoS Comput Biol. 2015; 11: e1004576. https://doi.org/10.1371/journal.pcbi.1004576 PMID: 26623645

45. Scott G, Fagerholm ED, Mutoh H, Leech R, Sharp DJ, Shew WL, et al. Voltage imaging of waking mouse cortex reveals emergence of critical neuronal dynamics. J Neurosci. 2014; 34: 16611–20. https://doi.org/10.1523/JNEUROSCI.3474-14.2014 PMID: 25505314

46. Priesemann V, Valderrama M, Wibral M, Le Van Quyen M. Neuronal Avalanches Differ from Wakefulness to Deep Sleep–Evidence from Intracranial Depth Recordings in Humans. Sporns O, editor. PLoS Comput Biol. 2013; 9: e1002985. https://doi.org/10.1371/journal.pcbi.1002985 PMID: 23555220

47. Fagerholm ED, Scott G, Shew WL, Song C, Leech R, Knöpfel T, et al. Cortical Entropy, Mutual Information and Scale-Free Dynamics in Waking Mice. Cereb Cortex. 2016; 26: 3945–52. https://doi.org/10.1093/cercor/bhw200 PMID: 27384059
48. Tagliazucchi E, Chialvo DR, Siniatchkin M, Amico E, Brichant J-F, Bonhomme V, et al. Large-scale signatures of unconsciousness are consistent with a departure from critical dynamics. J R Soc Interface. The Royal Society; 2016; 13: 20151027. https://doi.org/10.1098/rsif.2015.1027 PMID: 26819336

49. Meisel C, Olbrich E, Shriki O, Achermann P. Fading signatures of critical brain dynamics during sustained wakefulness in humans. J Neurosci. 2013; 33: 17363-72. https://doi.org/10.1523/JNEUROSCI.1516-13.2013 PMID: 24174669

50. Fagerholm ED, Lorenz R, Scott G, Dinov M, Hellyer PJ, Mirzaei N, et al. Cascades and Cognitive State: Focused Attention Incurs Subcritical Dynamics. J Neurosci. 2015; 35: 4626-4634. https://doi.org/10.1523/JNEUROSCI.3694-14.2015 PMID: 25788679

51. Dehghani N, Hatsopoulos N, Haga Z, Parker R, Greger B, Halgren E, et al. Avalanche Analysis from Multielectrode Ensemble Recordings in Cat, Monkey, and Human Cerebral Cortex during Wakefulness and Sleep. Front fractal Physiol. 2012; 3: 302. https://doi.org/10.3389/fphys.2012.00302 PMID: 22934053

52. Klaus A, Yu S, Plenz D. Statistical Analyses Support Power Law Distributions Found in Neuronal Avalanches. Zochowski M, editor. PLoS One. 2011; 6: e19779. https://doi.org/10.1371/journal.pone.0019779 PMID: 21720544

53. Poil S-S, Hardstone R, Mansvelder HD, Linkenkaer-Hansen K. Critical-State Dynamics of Avalanches and Oscillations Jointly Emerge from Balanced Excitation/Inhibition in Neuronal Networks. J Neurosci. 2012; 32: 9817-9823. https://doi.org/10.1523/JNEUROSCI.5990-11.2012 PMID: 22815496

54. Eurich C, Herrmann J, Ernst U. Finite-size effects of avalanche dynamics. Phys Rev E. APS; 2002; 66: 66137. Available: http://link.aps.org/doi/10.1103/PhysRevE.66.066137

55. Marshall N, Timme NM, Bennett N, Ripp M, Lautzenhiser E, Beggs JM. Analysis of Power Laws, Shape Collapses, and Neural Complexity: New Techniques and MATLAB Support via the NCC Toolbox. Front Physiol. 2016; 7: 250. https://doi.org/10.3389/fphys.2016.00250 PMID: 27445842

56. Teich M, Heneghan C, Lowen S, Ozaki B, Kaplan E. Fractal character of the neural spike train in the visual system of the cat. J Opt Soc Am. 1997; 14: 529-46. Available: http://www.opticsinfobase.org/abstract.cfm?URI=josaa-14-3-529

57. Safonov L a., Isomura Y, Kang S, Struzik ZR, Fukai T, Câteau H. Near Scale-Free Dynamics in Neural Population Activity of Waking/Sleeping Rats Revealed by Multiscale Analysis. Sporns O, editor. PLoS One. 2010; 5: e12869. https://doi.org/10.1371/journal.pone.0012869 PMID: 20927400

58. Segev R, Benveniste M, Hulata E, Cohen N, Palevski A, Kapon E, et al. Long term behavior of lithographically prepared in vitro neuronal networks. Phys Rev Lett. 2002; 88: 118102. Available: http://www.ncbi.nlm.nih.gov/pubmed/11909430 https://doi.org/10.1103/PhysRevLett.88.118102

59. Tkačik G, Mora T, Marre O, Amodei D, Palmer SE, Berry MJ, et al. Thermodynamics and signatures of criticality in a network of neurons. Proc Natl Acad Sci. 2015; https://doi.org/10.1073/pnas.1514188112

60. Cohen M, Kohn A. Measuring and interpreting neuronal correlations. Nat Neurosci. Nature Publishing Group; 2011; 14: 811-819. https://doi.org/10.1038/nn.2842 PMID: 21709677

61. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci. 2000; 8: 183-208. Available: http://www.ncbi.nlm.nih.gov/pubmed/10809012

62. Touboul J, Destexhe A. Power-law statistics and universal scaling in the absence of criticality. Phys Rev E. 2017; 95: 012413. https://doi.org/10.1103/PhysRevE.95.012413 PMID: 28208383

63. Klimesch W. Alpha-band oscillations, attention, and controlled access to stored information. Trends Cogn Sci. Elsevier Ltd; 2012; 16: 606-617. https://doi.org/10.1016/j.tics.2012.10.007 PMID: 23141428

64. Ecker AS, Berens P, Cotton RJ, Subramaniyan M, Denfield GH, Cadwell CR, et al. State dependence of noise correlations in macaque primary visual cortex. Neuron. Elsevier Inc.; 2014; 82: 235-248. https://doi.org/10.1016/j.neuron.2014.02.006 PMID: 24698278

65. Benayoun M, Cowan J, van Drongelen W, Wallace E. Avalanches in a Stochastic Model of Spiking Neurons. Friston KJ, editor. PLoS Comput Biol. 2010; 6: e1000846. https://doi.org/10.1371/journal.pcbi.1000846 PMID: 20628615

66. Ecker A, Berens P, Keliris G, Bethge M, Logothetis N, Tolias A. Decorrelated neuronal firing in cortical microcircuits. Science (80-). 2010; 327: 584-7. https://doi.org/10.1126/science.1179867 PMID: 20110506

67. Beggs J, Timme N. Being critical of criticality in the brain. Front fractal Physiol. 2012; 3: 163. https://doi.org/10.3389/fphys.2012.00163 PMID: 22701101

68. Schurger A, Sarigiannidis I, Naccache L, Sitt JD, Dehaene S. Cortical activity is more stable when sensory stimuli are consciously perceived. Proc Natl Acad Sci U S A. 2015; 112: E2083-2092. https://doi.org/10.1073/pnas.1418730112 PMID: 25847997

69. Zagha E, McCormick D a. Neural control of brain state. Curr Opin Neurobiol. Elsevier Ltd; 2014; 29: 178-186. https://doi.org/10.1016/j.conb.2014.09.010 PMID: 25310628
70. Goard M, Dan Y. Basal forebrain activation enhances cortical coding of natural scenes. Nat Neurosci. 2009; 12: 1444–9. https://doi.org/10.1038/nn.2402 PMID: 19801988

71. Marguet SL, Harris KD. State-dependent representation of amplitude-modulated noise stimuli in rat auditory cortex. J Neurosci. 2011; 31: 6414–20. https://doi.org/10.1523/JNEUROSCI.5773-10.2011 PMID: 21525282

72. Pachitariu M, Lyamzin DR, Sahani M, Lesica N. State-Dependent Population Coding in Primary Audi-tory Cortex. 2015; 35: 2058–2073. https://doi.org/10.1523/JNEUROSCI.3318-14.2015 PMID: 25653363

73. Kumar A, Rotter S, Aertsen A. Spiking activity propagation in neuronal networks: reconciling different perspectives on neural coding. Nat Rev Neurosci. Nature Publishing Group; 2010; 11: 615–627. https://doi.org/10.1038/nrn2886 PMID: 20725095

74. Fries P. Rhythms for Cognition: Communication through Coherence. Neuron. Elsevier Inc.; 2015; 88: 220–235. https://doi.org/10.1016/j.neuron.2015.09.034 PMID: 26447583

75. Hahn G, Bujan AF, Frégnac Y, Aertsen A, Kumar A. Communication through Resonance in Spiking Neuronal Networks. PLoS Comput Biol. 2014; 10: e1003811. https://doi.org/10.1371/journal.pcbi.1003811 PMID: 25165853

76. Munk M, Roelfsema P, König P, Engel A, Singer W. Role of reticular activation in the modulation of intra-cortical synchronization. Science (80-). 1996; 272: 271–4. Available: http://www.ncbi.nlm.nih.gov/pubmed/8602512

77. Herculano-Houzel S, Munk M, Neuenschwander S, Singer W. Precisely synchronized oscillatory firing patterns require electroencephalographic activation. J Neurosci. 1999; 19: 3992–4010. Available: http://www.ncbi.nlm.nih.gov/pubmed/10234029 PMID: 10234029

78. Williams-García R V, Moore M, Beggs JM, Ortiz G. Quasicritical brain dynamics on a nonequilibrium Widom line. Phys Rev E Stat Nonlin Soft Matter Phys. 2014; 90: 062714. Available: http://www.ncbi.nlm.nih.gov/pubmed/25615136 https://doi.org/10.1103/PhysRevE.90.062714

79. Moretti P, Muñoz MA. Griffiths phases and the stretching of criticality in brain networks. Nat Commun. 2013; 4: 2521. https://doi.org/10.1038/ncomms3521 PMID: 24088740

80. Ódor G, Dickman R, Ódor G. Griffiths phases and localization in hierarchical modular networks. Sci Rep. 2015; 5: 14451. https://doi.org/10.1038/srep14451 PMID: 26399323

81. Hilgetag CC, Hütt M-T. Hierarchical modular brain connectivity is a stretch for criticality. Trends Cogn Sci. 2014; 18: 114–115. https://doi.org/10.1016/j.tics.2013.10.016 PMID: 24268289

82. Dunn J. A Fuzzy Relative of the ISODATA Process and Its Use in Detecting Compact Well-Separated Clusters. J Cybern. 1973; 3: 32–57.

83. Alstott J, Bullmore E, Plenz D, Michel M, Kirk H, Myers P, et al. powerlaw: A Python Package for Analysis of Heavy-Tailed Distributions. Rapallo F, editor. PLoS One. Public Library of Science; 2014; 9: e85777. https://doi.org/10.1371/journal.pone.0085777 PMID: 24489671

84. Schneidman E, Berry MJ, Segev R, Bialek W. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature. 2006; 440: 1007–12. https://doi.org/10.1038/nature04701 PMID: 16625187

85. Marre O, El Boustani S, Frégnac Y, Destexhe A. Prediction of spatiotemporal patterns of neural activity from pairwise correlations. Phys Rev Lett. 2009; 102: 138101. Available: http://www.ncbi.nlm.nih.gov/pubmed/19392405 https://doi.org/10.1103/PhysRevLett.102.138101

86. Kumar A, Rotter S, Aertsen A. Conditions for propagating synchronous spiking and asynchronous firing rates in a cortical network model. J Neurosci. 2008; 28: 5268–80. https://doi.org/10.1523/JNEUROSCI.2542-07.2008 PMID: 18480283

87. Stimberg M, Goodman DFM, Benichoux V, Brette R. Brian 2—the second coming: spiking neural network simulation in Python with code generation. BMC Neurosci. BioMed Central; 2013; 14: P38.