Criticality of neuronal avalanches in human sleep and their relationship with sleep macro- and micro-architecture

Summary Sleep plays a key role in preserving brain function, keeping brain networks in a state that ensures optimal computation. Empirical evidence indicates that this state is consistent with criticality, where scale-free neuronal avalanches emerge. However, the connection between sleep architecture and brain tuning to criticality remains poorly understood. Here, we characterize the critical behavior of avalanches and study their relationship with sleep macro- and micro-architectures, in particular, the cyclic alternating pattern (CAP). We show that avalanches exhibit robust scaling behaviors, with exponents obeying scaling relations consistent with the mean-field directed percolation universality class. We demonstrate that avalanche dynamics is modulated by the NREM-REM cycles and that, within NREM sleep, avalanche occurrence correlates with CAP activation phases—indicating a potential link between CAP and brain tuning to criticality. The results open new perspectives on the collective dynamics underlying CAP function, and on the relationship between sleep architecture, avalanches, and self-organization to criticality.


Highlights
During sleep avalanches follow scaling laws in line with a critical branching process The measured critical exponents satisfy theoretically predicted scaling relations

Critical exponents and scaling relations for neuronal avalanches during sleep
We analyze brain activity recorded from 10 healthy subjects during night sleep (EEG: 16 channels in three subjects, 19 channels in four subjects and 25 in the remaining three; STAR methods).To measure neuronal avalanches across EEG signals, we first identify prominent deflections in the continuous signal collected from each sensor (Figure 1).We observe that the signal amplitude distributions differ from a Gaussian fit for values larger than 2.0 SD (Figure 1B).Deviations of the EEG amplitude distributions from a Gaussian behavior indicate presence of spatiotemporal correlations that signal the occurrence of clustered, synchronized neural activity. 33,54,59,60Thus, for further analysis, we identify large positive and negative signal deflections using an amplitude threshold q = 2 SD (STAR methods; see supplementary information for robustness of avalanche statistics with respect to q values).
We define an avalanche as a continuous time interval in which there is at least one excursion beyond threshold in at least one EEG channel (Figure 1C, shaded regions).Avalanches are preceded and followed by time intervals with no excursions beyond threshold on any EEG channel. 30,54The size of an avalanche, s, is defined as the sum over all channels of the absolute values of the signals exceeding the threshold (Figure 1C, bottom).
To characterize cortical dynamics underlying sleep macro-and micro-architecture, we identify neuronal avalanches and investigate signatures of criticality across the entire sleep period.To this end, we compute the distribution of avalanche sizes, PðsÞ, and avalanche durations, PðTÞ.In Figure 2 we show the distributions PðsÞ and PðTÞ for all subjects (pooled).We find that both the size and duration distributions are well described by a power law, PðsÞfs À t and PðTÞfT À a , respectively.In both distributions the power law regime is followed by an exponential cutoff (Figure 2).
Power laws are the hallmark of criticality, and imply absence of characteristic scales in the underlying dynamics. 58In this context, the observed power law distributions indicate that neuronal avalanches have no characteristic size and duration, namely they are scale-free.Our analysis shows that the exponent t for the size distribution is close to 3/2 (t = 1:438G0:001G0:0414; fit G error on the fit G SE), while the exponent a for the duration distribution is close to 2 (1:973G0:002G0:0452; fit G error on the fit G SE) (Figure 2).The power-law fits were performed over about three decades on the size distributions and two decades on the duration distributions.The limited range of the powerlaw regime has been associated with finite size effects. 30,33To account for the uncertainty due to the limited fit range, we added a systematic error (SE) to the power law exponent estimates (STAR methods).
We compared the power law with an exponential fit by evaluating the log likelihood ratio R = ln Lp Le between the likelihood L p for the power law and L e for the exponential fit (STAR methods).We found R = 295 for the size and R = 95 (p value < 10 À 5 ; see STAR methods) for the duration distribution, indicating that the respective power laws better describe the empirical distributions.Importantly, we observe that the power law exponents t and a are robust and weakly depend on the scale of analysis (supplementary information, Figure S3)-e.g. the threshold q used to identify avalanches--, and are consistent across subjects (Figure 3).In Figure 3 we show the avalanche size and duration distributions for each individual subject.Both distributions show little variability across subjects, and follow a power law with exponents t = 1:45G 0:03 and a = 1:96G0:05 (mean G SEM).These values are consistent with the values predicted within the mean-field directed percolation (MF-DP) universality class-3/2 and 2, respectively. 62We note that the cutoff following the power-law regime appears to deviate from an exponential behavior, particularly for very large avalanche durations (Figure 3B).However, in this range of values, distributions are rather noisy (see error (n = 10 subjects; for each subject we pooled all individual EEG channels).The black curve is the grand average over all subjects.The red curve is the best Gaussian fit for the grand average.We notice that the empirical probability density starts deviating from the Gaussian fit around G2 SD. (C) A neuronal avalanche is defined as a continuous sequence of signal excursions beyond threshold (red thick line) on one or more EEG channels (upper panel).An avalanche is preceded and followed by periods in which EEG signal are below the threshold in all channels.The size of an avalanche is defined as the sum over all channels of the absolute values of the signals exceeding the threshold (bottom panel).
bars in Figure S3), making it difficult to reliably assess the nature of this regime.The analysis of finite size effects shows that the onset of the cutoff depends on the size of the electrodes array (Figure S4), and that the tail shortens significantly for small array sizes, getting progressively closer to an exponential behavior.We cannot exclude that the deviation of the tails from the exponential behavior, observed at large scale with very low probability density (< 10 À 7 for the sizes and < 10 À 4 for the durations), may be related to intrinsic properties of the data, e.g., dominant, slow delta oscillations.
Moreover, we find that the avalanche branching parameter 22 is very close to 1, as predicted for critical branching processes (also in the MF-DP universality class), and weakly depends on q (supplementary information, Figure S3).
Next, we analyze the relationship between avalanche sizes and durations.Near criticality the average avalanche size CsD is expected to scale as a power of the duration T, namely CsDfT k 62 .We find that such a power law relationship holds during sleep (Figure 4).In particular, we observe that, for T's smaller than the duration corresponding to the onset of the cutoff regime in the distribution PðTÞ (Figures 2 and 3), the average size scales as CsDfT k with kz2 (Figure 4).For larger durations, we observe a crossover to a power law relationship with a smaller exponent kz1:3 (Figure 4).Importantly, the exponent k is robust and independent of the threshold q used to detect neuronal avalanches (supplementary information, Figure S3).Moreover, we observe that the relation CsDfT k is consistent across individual subjects (Figure 4B), the exponent k showing little variability across subjects.Specifically, we find k = 1:96G0:13 (mean G SD) for T's smaller than the duration corresponding to the onset of the exponential cutoff in the distribution PðTÞ, and k = 1:32G0:19 (mean G SD) for larger T's (Figure 4B).
Notably, we find that the exponent k measured in Figure 4 is in agreement, within errors, with the value predicted by the scaling relation which needs to hold at criticality.Indeed, we have that a À 1 t À 1 = 2:13G0:26 (mean G SEM), and k = 1:96G0:04 (mean G SEM).The scaling relation in Equation 1 has a general validity in avalanche dynamics, as shown in, 63,64 where Equation 1 was derived with the only hypothesis that PðsÞfs À t and PðTÞfT À a , and that the size fluctuations for fixed durations are small and can be neglected.
In sum, during sleep, the values of the critical exponents t, a and k are very close to the ones predicted for the mean field directed percolation (MF-DP) universality class, with exponents t = 3=2 for the size and a = 2 for the duration distribution, and k = 2 62 , and the avalanche branching parameter is very close to the critical value s = 1.

Avalanche dynamics and sleep macro-architecture
We have shown that, during sleep, neuronal avalanches are characterized by a robust scaling behavior in their size and duration distributions (Figures 2 and 3), and that avalanche size and duration are linked by precise scaling relationships (Equation 1; Figure 4).These observations are robust and consistent across subjects, and indicate underlying tuning to criticality during sleep.Next, we investigate the relationship between critical avalanche dynamics, sleep stages, and sleep stage transitions.
We first characterize sleep macro-architecture across all subjects.The main sleep parameters are described in Table 1 (macro-structural measures).The average total sleep time (TST) across the 10 subjects was 423.9 min, with a mean SE of 88.92%.Around 56% of TST was spent in light sleep (N1 = 7.23%, N2 = 48.47%),23.99% in deep sleep (N3 = 23.99%), and 20.30% in REM.TST for individual subjects is reported in the supplementary information (Table S1).
To study the interplay between sleep macro-architecture and avalanche dynamics, we introduce the avalanche density, F av ðtÞ, defined as the amount of time occupied by avalanches in a sliding window of length u 0 (STAR methods), and study the temporal evolution of F av ðtÞ in relation to the sleep macro-architecture.In the following we fix u 0 = 10 s, which approximately corresponds to the largest avalanche duration A B Figure 2. Avalanche size and duration distributions exhibit a robust power law behavior during sleep periods (A) The distribution of avalanche sizes (red circles) follows a power law with exponent t = 1:438G0:001G0:0414 (fit G error on the fit G SE; pooled data, 10 subjects).The power law regime is followed by an exponential cut off.The Kolmogorov-Smirnov distance between data and fit is D = 0:1, while the log likelihood ratio between the power law and the exponential fit is R = 295 (p < 10 À 5 ).
(B) The distribution of avalanche duration follows a power law with exponent a = 1:973G0:002G0:0452 (fit G error on the fit G SE), followed by an exponentiallike cutoff (pooled data, 10 subjects).The Kolmogorov-Smirnov distance between the data and the fit is D = 0:07.The log likelihood ratio between the power-law and the exponential fit is R = 95 (p < 10 À 5 ).Maximum likelihood estimation of the power law exponents was performed using the Powerlaw Python package 61 over the range of values indicated by the thick black lines.SE is the systematic error on the fit (STAR methods).
we observed (Figures 2 and 3).In Figure 5A we show the avalanche density F av ðtÞ as a function of time for an individual subject, together with the corresponding hypnogram.We observe that F av ðtÞ gradually increases in parallel with sleep deepening, i.e., going from REM to N1, N2, and finally N3: F av is very small during stage N1, reaches an intermediate value during stage N2, and increases substantially during stage N3, where it peaks slightly before the following transition back to N2 and REM (Figure 5A).Although the avalanche density tends to decrease across the night and is, on average, much smaller at the end of the night, we find that this trend repeats throughout the night in correspondence to the descending REM / N3 of the NREM-REM sleep cycle.In contrast to this gradually increasing trend, we observe that the avalanche density decreases rather abruptly with transitions from N3 to N2 and N1-the ascending phase of the NREM-REM sleep cycle.
In sum, we find that the avalanche density gradually increases during the descending slope of each sleep cycle, whilst it rapidly decreases in the ascending slope of the same cycles that precedes the onset of REM sleep (Figure 5A).
Our analysis shows that the density of avalanches is significantly higher during N3 as compared to N1 and REM (Figures 5B and 5C), while the pairwise comparison between N3 and N2 shows non-significant differences (p = 0:06; STAR methods).The analysis of the Pearson correlation coefficient r x;y (STAR methods, Equation 2) shows that avalanche occurrence, on average, is positively correlated with N3, the deepest sleep stage, while it is either weakly or slightly negatively correlated with other sleep stages (Figure 5D).Finally, we observe that, during N3, the avalanche density tends to increase with time (Figure 5A).
Interestingly, we observe that, for relatively low densities (F av < 0:4.), both the number of avalanches, N av , and the mean duration of avalanches, CTD, in the sliding window u 0 , tend to increase with the density F av (supplementary information, Figure S6).On the other hand, larger densities do not correspond to a consistent increase in N av , but are rather associated to longer avalanches.Because F av > 0:4 mostly occur during N3 sleep stage (see Figures 5A and 5B), this indicates the presence of longer avalanches during N3, which could be related to dominance of slow delta oscillations in this stage (see definition of avalanches; Figure 1 and STAR methods).On the other hand, the gradual increase of density from the low values in N1 to the intermediate values in N2 is accompanied by a gradual increase of the number N av of avalanches.
Importantly, we notice that the avalanche density peak-typically located within N3 periods-is higher in the first half of the night, d progressively decreases during the second half of the night.To quantify the significance of this behavior with respect to the characteristics of neuronal avalanches, we compare the avalanche density, as well as avalanche size and duration distributions, in the first and last N3 stage of the sleep recordings.We find that avalanche size and duration distributions in the first N3 are comparable to the distributions calculated in the last N3 (SI, Figure S5).Furthermore, the scaling relation CsDfT k between avalanche size and duration is satisfied both in the first and last N3, with the same values of the exponent k (SI, Figure S5).On the other hand, we observe that the avalanche density is significantly higher during the first N3 as compared to the last N3 (Figures 5E and 5F) (t-test: p = 0:04).This is consistent across subjects (Figure 5F), with the exception of subject #6 for which we observe that the density is higher in the last N3.Such deviation from the average behavior may be related to general differences we observed in sleep of subject #6.For instance, this subject presented an unusually short duration of the N3 stage at the beginning of the night, followed by a gradual increase of N3 in the second half of the sleep.

Avalanche dynamics and sleep micro-architecture
The analysis of the avalanche density across sleep stages has shown that neuronal avalanches tend to occur with higher frequency during NREM sleep.However, NREM sleep has a complex micro-architecture that is characterized by the CAP (see supplementary information, Figures S1 and S2 for a representative example of CAP oscillations and NCAP during NREM sleep). 3In our data, the mean CAP rate was 49.19% with the following distribution across NREM stages: N1 = 41.69%,N2 = 48.36%,and N3 = 53.37%( stage N3 (50.21%) as compared to N2 (5.72%) and N1 (1.49%), in agreement with previous studies. 65On the other hand, subtypes A2 and A3 predominated in stage N1 (particularly A3, 37.77%) and N2 (14.39% for A2 and 17.46% for A3) (Table 2).
To dissect the relationship between CAP and occurrence of neuronal avalanches during NREM sleep, we compare the time course of the avalanche density with the density of distinct CAP phases (Figure 6A) defined as F X ðtÞ = ðu X ðtÞÞ=u 0 , where X denotes the specific CAP phase-A, A1, A2, A3, B-and u X ðtÞ the time occupied by the specific CAP phase in a window of length u 0 = 10 s.We observe a remarkable time correspondence between the temporal profile of the density of avalanches F av ðtÞ and the density of CAP, with the peaks in avalanche density corresponding to high density of CAP-in particular phase A and A1 (Figure 6A).Specifically, we notice that, with sleep deepening, the progressive increase of CAP density is accompanied by a parallel increase in avalanche density.We find that the percentage of phase A occupied by neuronal avalanches is about 42.16%, while the percentage of sleep time occupied by avalanches is 19.21% (STAR methods).Interestingly, CAP phase A1 is even richer in avalanches compared to CAP A phases A2 and A3 (53.32% versus 43.84% and 27.72%, respectively).
The physiological increase of CAP cycles during N2 and N3, indirectly leads to a reduction of time occupied by NCAP sleep.Furthermore, during the deepest stages of NREM sleep, CAP's typically present shorter phases B. These changes in the sleep micro-dynamics lastly sustain the observed increase of avalanche density.
Next, we measure the Pearson correlation coefficients between occurrence of neuronal avalanches and different CAP phases (see STAR methods, Equation 2).We find positive correlations between occurrence of avalanches and CAP phase A, in particular CAP phase A1 (Figure 6B).On the contrary, we observe negative correlations between occurrence of avalanches, CAP phase B, NCAP periods.This indicates that the occurrence of avalanches during NREM sleep is strictly related to occurrence of CAP, and in particular CAP phase A1.These results are consistent across subjects, as shown in Table 3.

DISCUSSION
In this paper we analyzed the scaling properties of neuronal avalanches during sleep in healthy volunteers, and investigated the relationship between avalanche dynamics and sleep macro-and micro-architecture, with a particular focus on the CAPs.We showed that the scaling exponents characterizing neuronal avalanches are consistent with the MF-DP universality class, and obey the theoretically predicted scaling relations.This indicates that, during physiological sleep, brain dynamics is consistent with criticality and is satisfactorily described by the MF-DP universality class.Furthermore, we introduced a measure-the density of avalanches-to quantify the relationship between avalanche dynamics and sleep macro-and micro-architecture.Our analysis showed that the distribution of avalanches in time is not random but closely follows the descending and ascending phase of the NREM-REM cycles.Within such cycles, the presence of neuronal avalanches is linked to the occurrence of CAP during NREM sleep.Specifically, we found that the density of avalanches is higher during NREM, and, within NREM sleep, avalanche occurrence is positively correlated with the phase A of the CAP, in particular the phase A1.This suggests a close relationship between modulation and control of brain criticality, sleep macro-and micro-architecture, and brain function, which we discuss in turn.

Brain dynamics and criticality during sleep
][52] In particular, recent studies suggest that criticality plays a key role in determining the temporal organization of sleep stage and arousal transitions, 11,12 and sleep deprivation progressively disrupts signatures of criticality 54 and alters brain connectivity. 66 easily accessible otherwise, or even expected, opening the way to further experimental and theoretical development.This is made possible by a set of scaling relationships connecting the critical exponents.Specifically, for the directed percolation (DP) universality class there are only three independent exponents, and all others can be derived from known scaling relationships. 67For systems out of equilibrium, as the brain, the DP universality class is expected to describe any absorbing phase transition-i.e., a transition between an active and an inactive (absorbing) state-with just one absorbing state (or more but non-equivalent).Systems with more than one equivalent absorbing state falls in different universality classes, and thus are described by a different set of exponents. 68o the best of our knowledge, this is the first study exploring the scaling relations among critical exponents of neuronal avalanches during sleep.We reported a picture that is consistent with the MF-DP universality class.Indeed, we have shown that (1) the critical exponents for the avalanche size and duration distributions are very close to the prediction of MF-DP universality class, i.e., t = 3=2, a = 2, respectively; (2) the exponent k connecting sizes and durations is very close to 2, as predicted; (3) the exponents t, a, and k correctly satisfy the expected scaling relation.The fact that sleep criticality in healthy subjects seems to be in line with the MF-DP class is important to assess deviations from criticality in pathological sleep, and a key step toward sleep biomarkers based on EEG criticality measures.
0][71][72] In line with our findings, recent works in awake nonhuman primates 69 and awake mice 44 found that kx2 in the range corresponding to the power law regime of the size and duration distributions, and k ˛½1; 1:5 only in the region that corresponds to the cut-off of the distributions-where we found kx1:3.Similar results were found in zebrafish. 37Deviation from the MF-DP value k = 2 was observed in the resting-state of the human brain, 71 in ex vivo turtle visual cortex, 73 in the barrel cortex of anesthetized rats, 42 in cortex slice cultures, 74 and in freely behaving and anesthetized rats. 70Sub-sampling in brain activity recordings has been suggested as a potential origin of such deviation from the theoretically predicted value, 75,76 and could affect scaling exponents in our analysis (we have between 16 and 25 electrodes).Alternatively, a recent work has shown that, in a 2D Wilson-Cowan neural network, the value of the exponent k is related to the network connectivity, with kx 1:3 for a 2D connectivity, and k = 2 when the mean-field approximation holds. 77This may suggest that, in our case, the mean-field approximation is justified for relatively small avalanches involving few electrodes, but not for large-scale EEG avalanches.However, to verify whether this is due to the structure of the underlying, collective neural activity, simultaneous multi-scale recordings would be needed.In this respect, an important open question is the relationship between avalanches across scales (from spike avalanches to EEG/MEG)-a key question to correctly interpret discrepancies in scaling exponents across experiments.Alternatively, a recent work has shown that, in a 2D Wilson-Cowan neural network, the value of the exponent k is related to the network connectivity, with kx1:3 for a 2D connectivity, and k = 2 when the mean-field approximation holds. 77This may suggest that, in our case, the mean-field approximation is justified for relatively small avalanches involving few electrodes, but not for large-scale EEG avalanches.However, to verify whether this is due to the structure of the underlying, collective neural activity, simultaneous multi-scale recordings would be needed.In this respect, an important open question is the relationship between avalanches across scales (from spike avalanches to EEG/MEG)-a key question to correctly interpret discrepancies in scaling exponents across experiments.Moreover, the presence of a large-avalanche regime with kx1 after a regime with kx2 has been also observed at the critical point, near a bistability regime in a stochastic Wilson-Cowan model whose function of activation mimics cooperative effects. 78n this model, the regime with kx1 is related to closeness of the system operating point to an underlying self-sustained regime.In sum, together with subsampling and coarse-graining effects, further investigation is needed to understand the crossover from kx2 to kx 1. Analyses of scalp EEG and human intracranial depth recordings have shown that avalanche size and duration distributions follow a similar power law behavior across the sleep-wake cycle, with exponents in line with our observations. 50,52Similarly, the analysis of whole-brain fMRI data have confirmed a critical (or near-critical) behavior from wakefulness to deep sleep, with little differences in the power-law exponent of the avalanche size distribution (in particular between wakefulness and stage N2). 51n the other hand, here we have shown that, although the static properties remain fairly stable across different sleep stages, [50][51][52] avalanche dynamics is modulated by the ascending and descending slope of the NREM-REM sleep cycles.By analyzing the temporal evolution of the avalanche density, we found that avalanche occurrence markedly and progressively increases with NREM sleep stages N2 and N3 and, specifically, during periods of sleep deepening (descending slope of sleep cycles), in parallel with the increase of SWA.On the contrary, the abrupt decrease in avalanche density during the ascending slope of sleep cycles suggests a negative influence from REM-on/wakefulness circuits with respect to their appearance.The different behavior of avalanche density during the descending and ascending slopes of the sleep cycles was not previously observed, despite the crucial role of such dynamics for sleep regulation.In terms of sleep physiology, the descending and ascending slopes of sleep cycles are markedly different: during the descending slope, sleep-promoting forces, i.e., homeostatic sleep propensity, are stronger, 79,80 the thalamocortical system works in the burst-firing mode and brainstem cholinergic pathways are tonically repressed.Conversely, during the ascending slope, the NREM driving forces become weaker, sleep is more vulnerable toward pro-arousal intrusions and REM-promoting outputs prevail. 65Taking this into account, our results suggest that avalanche occurrence is not random across the sleep cycles, but instead contributes to define the dynamical interplay between sleep-wake promoting networks.
The reported scale-free avalanche activity within the NREM sleep coexists with SWA which, in contrast, has a characteristic timescale.The coexistence of scale-free and scale-specific brain activity, and specifically the relationship between neuronal avalanches and oscillations, is still poorly understood and deserves further investigation.In a recent human study, it has been shown that the dynamics of avalanches depends on the dominant, coexisting brain oscillations, and differs between awake resting-state and NREM sleep. 60During resting wake, alpha oscillations induce attenuation of avalanches (i.e., consecutive sizes tend to be smaller) within a single alpha cycle, and corresponding amplification over several alpha cycles (i.e., consecutive sizes tend to be smaller).In contrast, during NREM sleep, only the attenuation regime has been observed, both at short and long timescales. 602][83] In rodents, nested q and b/g oscillations were found to be embedded in neuronal avalanches, 81 while in cortex slice cultures a hierarchical organization of avalanches in relation to q and g oscillations was described. 82rom a theoretical perspective, such a coexistence has been studied in a few models, 71,[84][85][86][87] and quantitatively captured in the resting human brain by an adaptive Ising model indicating that the coexistence of oscillations and avalanches in brain activity occurs close to a non-equilibrium critical point at the onset of self-sustained oscillations. 71,88Yet, the precise relationship between collective oscillations and avalanches remains an open question.

Avalanches and sleep micro-architecture
Sleep architecture is composed of numerous oscillatory patterns, including, above all, the CAP. 89CAP's occur on time scales of seconds or minutes, accompany sleep stage shifts, and contribute to the organization of sleep cycles.Our analyses demonstrated positive correlations between CAP and avalanche occurrence, and negative correlations for NCAP sleep.Such link suggests a close relationship between CAP and brain tuning to criticality during sleep, a key aspect that should be further investigated in future work.Although the definition of avalanches (large, collective non-gaussian fluctuations of brain activity) is not related to the definition of CAP phase A, our results show that neuronal avalanches are correlated with the occurrence of CAP phase A. In particular, we observed stronger correlations between avalanche occurrence and the CAP A1 subtype, and weaker positive correlation with subtypes A2 and A3.Interestingly, the correlation between avalanches and the phase A of the CAP is more prominent than the correlation with the CAP itself-phase A and phase B together.We speculate that this could be due to the opposite significance of CAP phase A and B with respect to sleep dynamics.Electrophysiologically the phase B is characterized by the rebound of background EEG activity after the strong ''activation'' driven by the phase A. Compared to phase A, the phase B could be described as a ''lower arousal reaction'' or mechanism of deactivation. 90Phases B reflect a period of transient inhibition and have been associated with nocturnal sleep apnea in patient with sleep-breathing disorder as well as with inhibitory effect on nocturnal epileptiform discharges in patients with generalized epilepsies. 89,91Importantly, we did not observe significant correlation between avalanche occurrence and phase B, corroborating our assumption about the relationship between CAP phase A and avalanches.The prominent correlation between avalanche occurrence and CAP ''activation phase'' A1 may suggest that neuronal avalanches emerge at the edge of a synchronization phase transition, as recent numerical studies indicate. 85,86,92inally, we note that CAP-A1 physiologically prevail in the first half of the night and during the descending slope of each sleep cycle, boosting or maintaining SWS.Similarly, the avalanche density decreases moving from the first to the last sleep cycle.Hence, both CAP phase A and neuronal avalanches follow a physiological, homeostatic decay throughout the night.

Neuronal avalanches, CAP, and learning mechanisms: An intriguing hypothesis
Sleep is crucial to renormalize synaptic weight, ensure an optimal and effective network state for information processing, and preserve cognition. 9Renormalization of synaptic weights taking place during sleep may serve to keep the network close to criticality. 93In line with this view, the here reported higher concentration of avalanches during SWS and CAP-A1 indicate that these states may exert a pivotal role in modulating and restoring brain criticality.Furthermore, because CAP-A1 has been proposed to play a role in the sleep-dependent learning processes, 8 our observations point to a functional link between critical avalanche dynamics and sleep-dependent learning processes, as shown in recent numerical studies. 92,94Specifically, it has been demonstrated that, within the alternation of up-and down-states observed during SWS, the sequence of avalanches occurring in the up-states correspond to an intermittent reactivation of stored spatiotemporal patterns, a mechanism that is key for memory consolidation. 95

Conclusions
Overall, our findings open a novel perspective on the relationship between critical brain dynamics and physiological sleep.We provided a comprehensive account of the critical exponents and scaling relations for neuronal avalanches, demonstrating that brain dynamics during sleep is consistent with the MF-DP universality class.This sets the bases for future investigation of neural collective behaviors occurring during sleep, including their functional role in relation to criticality.As a first step in this direction, our study provides evidence of a functional link between avalanche occurrence, slow-wave sleep dynamics, sleep stage transitions and occurrence of CAP phase A during NREM sleep.As CAP is considered one of the major guardians of NREM sleep 96 that allows the brain to react dynamically to any external perturbation and contributes to the cognitive consolidation processes occurring in sleep, we speculate that neuronal avalanches at criticality might be associated with flexible response to external inputs and to cognitive processes-a key assumption of the critical brain hypothesis.1][12] Moreover, based on our results, one could speculate that a relationship between occurrence of neuronal avalanches and physiological sleep measures exists.To address this point, additional studies in pathological sleep conditions where both CAP and criticality-based metrics show a deviation from the physiological parameters are needed. 96,97

Limitations of the study
The question of the universality class of neuronal avalanches, both in sleep and wake, is still debated. 44,70,76,86Although our results show that neuronal avalanches during sleep are described by scaling exponents and scaling relationships consistent with the MF-DP universality class, this is not sufficient to prove that during sleep neuronal avalanches belong to the MF-DP universality class.The main reasons are: (1) although most non-equilibrium processes with an absorbing phase transition belong to the DP universality class, we cannot exclude, based on our estimates of the exponents, that neuronal avalanches do not belong to a different universality class with similar mean-field exponents, e.g., the mean-field dynamical percolation. 56,67Further investigation of this point would benefit from multi-scale recordings to measure other critical exponents that could confirm (or reject) the MF-DP as the universality class for neuronal avalanches; (2) we define neuronal avalanches as spatiotemporal sequences of threshold-crossing events across EEG sensor array-a sort of coarse-grained measure-but we do not have access to the actual activity of underlying neural populations and, thus, we cannot prove that the same results hold for spike avalanches.Simultaneous multiscale recordings would be needed to demonstrate consistency with the MF-DP across scales, also taking into account potential subsampling effects. 75Indeed, in the awake state, deviation from the MF-DP exponents in spike avalanches in mice 70 and in zebrafish 37 has been reported.However, this deviation may arise from subsampling issues. 98inally, the relatively small number of subjects requires a word of caution about the presented results.Although our individual recordings are rather long (424G63 min, Table 1) and provide a fairly robust individual statistics for avalanches (> 10 5 avalanches per subject), sleep stages (NREM: 335G42 min; mean G SD), and CAP (CAP time: 162G42 min; mean G SD), our cohort comprises 10 subjects with a relatively wide age range-from 28 to 53 years.This limits the statistical power of our analysis and, in particular, does not allow stratification of results according to age, which may be relevant in this context.We observed that results for avalanche dynamics were robust across subjects (Figures 3 and 4; Figure S7), with small variability in the scaling exponents.Furthermore, we found that avalanche size and duration distributions of the youngest and oldest subjects (28 vs. 53 years) are very close to each other, showing no larger differences than those observed between the two youngest subjects (28 years) (Figure S7).Similarly, we observed that the relationship between avalanches and sleep architectures is consistent across subjects (Figures 5B and 5C; Table 3), with no apparent age-related effects when comparing the youngest and oldest subjects (Figure S8).However, because of the limited number of subjects, we cannot exclude that avalanche dynamics and its relationship with sleep architecture depends on aging.Future work on larger cohorts-including a significant number of healthy subjects in each age group-is needed to address this key point, as well as to avoid potential small-sample biases (e.g., due to the intrinsic nature of the EEG) and thus confirm the presented results.

Figure 1 .
Figure 1.Identification of neuronal avalanches and definition of avalanche size and duration (A) Segments (2 h) of Z score normalized EEG signal traces for an individual subject.Each trace corresponds to an EEG channel.(B) Probability density of the Z score normalized EEG signal amplitude.The cyan curves in the background are the probability densities for all individual subjects(n = 10 subjects; for each subject we pooled all individual EEG channels).The black curve is the grand average over all subjects.The red curve is the best Gaussian fit for the grand average.We notice that the empirical probability density starts deviating from the Gaussian fit around G2 SD. (C) A neuronal avalanche is defined as a continuous sequence of signal excursions beyond threshold (red thick line) on one or more EEG channels (upper panel).An avalanche is preceded and followed by periods in which EEG signal are below the threshold in all channels.The size of an avalanche is defined as the sum over all channels of the absolute values of the signals exceeding the threshold (bottom panel).

Figure 4 .
Figure 4. Avalanche sizes and durations are connected by the scaling relationship CsDfT k consistent with underlying criticality (A) Average avalanche size as a function of the avalanche duration T (red stars; pooled data, 10 subjects).The average avalanche size scales as CsDðTÞf T k with k = 1:89 for T's within the scaling regime of the distribution PðTÞ.This power-law regime is followed by a crossover to a power-law with a significantly smaller exponent k = 1:4 for larger T's.The thick black line is a power law fit for 0:01 < T < 0:4 s; dashed black line is a power law fit for 0:4 < T % 5. Blue dots: (s,T) scatterplot.(B) Average avalanche size as a function of the avalanche duration T for all individual subjects.The relationship between avalanche sizes and durations is consistent across subjects, showing a crossover from an exponent k = 1:96G0:13, and k 1 = 1:32G0:19 (mean G SD).

Figure 5 .
Figure 5. Overnight sleep macro-architecture is associated with strong modulation of avalanche dynamics (A) The density of avalanches (blue dots), F av ðtÞ, is shown as a function of time, together with the corresponding sleep stages and sleep stage transitions (REM, N1, N2, N3 black line) for an individual subject.F av ðtÞ increases gradually in N2 and N3, and then abruptly decreases when transitioning from N3 to either N2, N1 or REM.Waking periods during sleep have been removed.(B) Mean avalanche density for each sleep stage (REM, N1, N2, N3) averaged across subjects.The density F av ðtÞ is highest in N3 and gradually decreases for N2, N1, and REM (one-way ANOVA on ranks: p = 5,10 À 6 ; significant pairwise comparison: N3 vs. REM (p = 7,10 À 5 ), N3 vs. N1 (p = 3,10 À 5 ).(C) Mean avalanche density for each sleep stage and for each individual subject.The behavior observed for the group average is consistent across individual subject, N3 being the sleep stage with the highest density of avalanches.(D) The mean Pearson correlation coefficients r x;y (see Equation 2) between avalanche occurrence and sleep macro-architecture (namely REM, N1, N2, N3) shows that avalanches tend to occur mostly during N3 (One-way ANOVA group comparison: p = 3,10 À 23 ).Pairwise differences between N3 and all other sleep stages are significant (N3 versus N2: p = 10 À 21 ; N3 versus N1: p = 3,10 À 19 ; N3 versus REM: p = 10 À 21 ).N2 is significantly different from N1 (p = 0:006), and N1 is significantly different from REM (p = 0:013).(E) Mean density of avalanches in the first and last N3 stage of the recordings averaged over all subjects.We observe that the density is significantly higher during the first N3 (p = 0:04).(F) Avalanche density in the first N3 (blue) and last N3 (red) for each individual subject.The density is higher in the first N3 for all subjects but the subject #6.Error bars indicate the SE.Significance legend: *** for p < 0:001; ** for p < 0:01; * for p < 0:05.The *** in panel D refers to the pairwise comparison between N3 and all the other sleep stages.Differences are not significant where no stars are reported.

Figure 6 .
Figure 6.Occurrence of neuronal avalanches is coupled with the occurrence of the CAP (A) Density of avalanches versus density of CAP phases as function of time for an individual subject.Density of avalanches in blue, density of phase A in red, density of phase A1 in green, density of CAP in black.(B) The mean Pearson correlation coefficients rðx; yÞ (average over subjects; see STAR methods, Equation 2) between avalanche occurrence and microarchitecture features (NCAP, CAP, B, A, and A subtypes A3, A2, A1).Error bars indicate the SE.A one-way ANOVA test shows that group differences are significant (p = 9,10 À 16 ).Pairwise comparisons show significant differences for all couples (p < 0:02) but av-NCAP vs. av-B, av-CAP vs. av-A3, av-CAP vs. av-A2, av-B vs. av-A3, av-A2 vs. av-A3, and av-A1 vs. av-A.Significance legend: *** for p < 0:001; ** for p < 0:01; * for p < 0:05.The ** in panel B refers to the pairwise comparison between av-A1 and all the other bars but av-A.Differences are not significant where no stars are reported.

Table 2 )
61On average, subjects presented 37.1 CAP sequences per night, with a mean duration of 4.55 min.With respect to CAP subtypes distribution, 206 were A1 Figure 3. Avalanche size and duration distributions consistently follow a power law behavior across individual subject (A) The distributions of avalanche sizes for individual subjects follow a power law with little variability across subjects (t = 1:45G0:09; mean over subjects G SD).(B) The distributions of avalanche durations for individual subjects consistently follow a power law, with little variability across subjects (a = 1:96G 0:16; mean over subjects G SD).For each individual subject, maximum likelihood estimation of the power law exponents was performed over the range of values corresponding to the thick black line using the Powerlaw Python package61.
(25.7% of the CAP time); 67.2 were A2 (9.2% of the CAP time), and 83.8 were A3 (14.19% of the CAP time).A1's were more present during A B

Table 1 .
Average characteristics of sleep macro-architecture across the analyzed subjects (n = 10)For each measure mean and standard deviation (SD) are reported.SE, sleep efficiency; TST, total sleep time; WASO, wake after sleep onset.
iScience Article Neuronal avalanches and sleep macro-architecture

Table 2 .
Average characteristics of sleep micro-architecture across the analyzed subjects (n = 10)

Table 3 .
Pearson correlation coefficient between avalanche occurrence and CAP subtypes for the analyzed individual subjects (n = 10) (STAR methods)