Non-equilibrium critical dynamics of bursts in θ and δ rhythms as fundamental characteristic of sleep and wake micro-architecture

Origin and functions of intermittent transitions among sleep stages, including short awakenings and arousals, constitute a challenge to the current homeostatic framework for sleep regulation, focusing on factors modulating sleep over large time scales. Here we propose that the complex micro-architecture characterizing the sleep-wake cycle results from an underlying non-equilibrium critical dynamics, bridging collective behaviors across spatio-temporal scales. We investigate θ and δ wave dynamics in control rats and in rats with lesions of sleep-promoting neurons in the parafacial zone. We demonstrate that intermittent bursts in θ and δ rhythms exhibit a complex temporal organization, with long-range power-law correlations and a robust duality of power law (θ-bursts, active phase) and exponential-like (δ-bursts, quiescent phase) duration distributions, typical features of non-equilibrium systems self-organizing at criticality. Crucially, such temporal organization relates to anti-correlated coupling between θ- and δ-bursts, and is independent of the dominant physiologic state and lesions, a solid indication of a basic principle in sleep dynamics.

Introduction The brain's ability to adapt and perform complex functions crucially depends on the cooperation of assemblies of neurons across multiple spatial and temporal scales. Cortical rhythms represent one of the most fascinating collective phenomena emerging from the self-organized synchronous activity of large neuronal populations, and are consistently associated with complex brain functions and distinct physiologic states such as sleep and wake [1,2]. Different brain rhythms characterize distinct phases of the sleep-wake cycle. During deep NREM sleep, brain dynamics are generally dominated by δ rhythm, low-frequency high-amplitude oscilla-have successfully accounted for consolidated sleep and wakefulness over time scales of hours and days, reproducing homeostatic, ultradian and circadian influences [9]. Further, a flip-flop switching mechanism involving mutually inhibitory interactions between sleep-and wakepromoting neuronal assemblies distributed across brain areas [19][20][21], and regulated by slow homeostatic factors such as adenosine and nitric oxide [9], has been proposed to explain the transition from stable sleep to stable wakefulness. However, the mechanism responsible for turning the switch on and off, the dynamic characteristics and temporal organization of these transitions remain unclear. Moreover, existing homeostatic models of sleep regulation (i) do not address empirical observations of transient behaviors at scales of seconds and minutesan intrinsic sleep micro-architecture on time scales much shorter than consolidated sleep and wake states which last several hours, and sleep-stage episodes lasting from many minutes to an hour; and (ii) do not account for the emergent complex structure of sleep stage and arousal/ brief wake transitions and the related micro-architecture of bursts in cortical brain waves activity. The intrinsic fluctuations in activity of brain rhythms in response to nonlinear feedback interactions among multi-component sleep-and wake-promoting neuronal pathways, the high susceptibility to abrupt transitions and the resulting complex temporal organization of sleep micro-architecture at scales of seconds and minutes rise the hypothesis that non-equilibrium critical dynamics may underlie sleep regulation at short time scales, in coexistence with the well-established homeostatic behavior at large time scales.
The traditional description of sleep stages in terms of dominant cortical rhythms indeed provides only an 'average', coarse-grained phenomenological picture that does not reflect the complex dynamics of sleep microstates and brief arousals. Current guidelines for sleep-stage classification do not consider the role of dynamical interactions among brain rhythms [3], although these interactions are potentially an integral feature of the collective neuronal activity driving physiologic states and related sleep-stage transitions [22][23][24][25]. Furthermore, the relation between emergent cortical brain-waves dynamics and neuronal activity in sleep-and wakepromoting areas remains largely unknown. Although many brain areas involved in sleep control have been identified [9,10], the complex nonlinear dynamics of extended neuronal networks (comprised of diverse mutually connected neuronal populations) continues to pose severe limitations to a mechanistic understanding of the collective neuronal behavior underlying observed sleep-stage transition patterns. This is a general problem that arises in systems where the emerging behavior at the system level originates from collective dynamics of a large number of interacting units. While elegant theories bridging micro-and macro-scales are sometimes available for such systems at thermodynamic equilibrium, emergent collective phenomena out of equilibrium often require a top-down approach, where basic mechanisms are inferred from the statistical characteristics of emergent dynamics.
Following this approach, here we investigate the temporal organization in bursting activity of brain rhythms across the sleep-wake cycle, with the aim to understand the basic physical principles bridging neuronal interactions at the microscopic level and physiologic state at the system level. The present study is motivated by recent works showing that the complex dynamics of sleep stage transitions give rise to power-law probability distributions for the durations of brief awakenings and arousals, a robust scale-invariant organization observed in both humans and animal models [13,[26][27][28][29][30][31]. Power-law distributions P(x) = Nx −α , where α is the scaling exponent and N is a normalization constant, are the statistical hallmark of scale invariance, i.e., they are not altered by a change of scale from x to Lx, hence lacking a relevant characteristic scale. Power laws are typical features of physical systems at the critical point of a second order phase transition in equilibrium thermodynamics [32,33]. At criticality systems exhibit high susceptibility and sensitivity to interactions among elements, leading to emergent collective behavior across scales, and thus, power laws. The critical point is located at the border between an ordered and a disordered phase, and can be reached by fine tuning external parameters. In contrast to this scenario, in non-equilibrium systems the dynamics can be spontaneously driven at criticality, where an active phase characterized by bursts/avalanches with power-law distributed sizes and durations coexists with a quiescent phase with exponential-like statistics [15,17,34,35]. We note that such avalanche criticality does not necessarily co-occur with edge-of-chaos criticality [36].
Since brief awakenings/arousals can be viewed as 'active' states of the brain that interrupt the 'inactive' phase represented by sleep periods, the scale-invariant organization of arousals has been interpreted as a fingerprint of criticality in sleep dynamics [13,26,31,37]. Contrary to the established interpretation of arousals as random and detrimental disruptions of sleep [38,39], caused by external stimuli or as part of the patho-physiology of sleep disorders [12,14,[40][41][42][43][44], such robust scale-invariant temporal organization indicates that arousals are an integral part of sleep regulation, resulting from a single dynamic principle responsible for the emergent temporal organization of sleep stages and sleep-stage transitions.
The hypothesis that sleep phenomenology and micro-architecture may reflect underlying scale-invariant dynamics resulting from self-tuning of a system at criticality [34] would imply presence of power-law distributions and long-range correlations as basic characteristics for sleep-related neuronal activity across spatial and temporal scales. Although power-law distributions [45][46][47] and long-range correlations [48] have been previously reported in neuronal and brain dynamics, and some attempts to correlate them with behavioral scaling laws has been made [49], the connection between (i) dynamics of sleep-related brain waves activity, (ii) neuronal circuitry and pathways related to sleep regulation, and (iii) emergent scale-invariant organization of brief arousals/awakenings, remains not understood.
To test this hypothesis, we consider the dominant brain waves in the sleep-wake cycle of rats, and study their dynamics and coupling in relation to the neuronal circuitry responsible for wake and sleep control. We focus on sleep-promoting PZ [50], a region of the rostral medulla, that plays a significant role in the regulation of slow-wave sleep [6][7][8]51]. Both PZ cell body specific lesion and disruption of GABAergic transmission in PZ result in insomnia, indicating that within the PZ, GABAergic neurons are primarily in sleep control. Similar to other NREM sleep promoting neurons, PZ GABAergic neurons project to and inhibit wake promoting systems, such as the parabrachial neurons that project to basal forebrain neurons, which in turn project to the cortex and are critical for cortical activation and wakefulness [6]. It has been shown that lesions of all PZ neurons result in a significant decrease of total NREM sleep during the light period (when rats are predominantly asleep), whereas chemogenetic activation of PZ neurons results in an increase of NREM sleep and slow-wave activity with dominant δ rhythm [6]. However, the influence of PZ sleep-promoting neurons on cortical brain waves activity remains not clear. In particular, the role played by PZ neurons in the dynamics and temporal organization of δ waves, as well as in the emergence of transient θ-bursts arousal activations during the sleep-wake cycle, has not been investigated.
To this end, we analyze long-term continuous EEG recordings in control rats and rats where the PZ brain area is lesioned, and we investigate the complex dynamics of θ-and δbursts in relation to PZ neuronal integrity, during both light and dark periods. In particular, we focus on the emergent scale-invariant features in the temporal organization of θ-and δ-bursts, and study whether alterations in the sleep-wake cycle are mirrored by a reorganization of dominant brain rhythms, a question that has not been addressed so far in sleep research. Our aim is to establish basic mechanisms that underlie cortical dynamics and sleep micro-architecture, and whether these dynamics exhibit characteristics of a system at criticality. Confirming our hypothesis, that cortical activations underlying sleep micro-architecture exhibit critical dynamics, would lay the foundation for a novel unified framework of the sleep-wake cycle, where sleep and arousals/wake, consolidated and transient states, originate from the same fundamental principles and a common mechanism that bridge collective behaviors across spatio-temporal scales, from neuronal assemblies to brain rhythms and emerging physiologic states.

Transient dynamics in bursting activity of θ and δ rhythms
To dissect the temporal organization of δ-and θ-bursts in the broadband brain activity during sleep and wake, we analyzed the time course of the EEG signal by evaluating the spectral power in several frequency bands on non-overlapping windows of length w (Materials and methods, Data analysis). In Fig 1a we show a typical spectrogram S(f) as a function of time for a 2 h recording of a rat in the control group. We notice that, in each window, the spectral power is primarily concentrated in either the δ-wave frequency range (0 − 4Hz) or in the θwave band (4 − 8Hz), and we observe sharp transitions from periods with dominant δ to periods with dominant θ waves. Such dynamics can be understood as the temporal evolution of the ratio R θδ = S(θ)/S(δ) between θ and δ spectral power in association with different physiological states-NREM, REM and arousals/wake ( The ratio R θδ (t) exhibits irregular, intermittent fluctuations between values larger and smaller than a threshold Th-a typical characteristic of non-equilibrium dynamics: R θδ > Th = 1 indicates that the spectral power in the θ-wave band is dominant; vice-versa, for R θδ < Th = 1 the spectral power is dominated by the δ-wave. We define bursts in θ and δ brain rhythms as sequences of consecutive time windows where R θδ > Th = 1 and R θδ < Th = 1, respectively (Fig 1b). The focus of this study is to establish the dynamical features and underlying mechanisms of bursts in cortical activity across the sleep-wake cycle by investigating the temporal organization of such bursts and their coupling. To this aim, we associate a duration d = n � w to each burst (Materials and methods, Data analysis), where n is the number of consecutive windows belonging to a given burst and w is the window length (Fig 1b).

Distinct functional forms of θ-and δ-burst duration distributions
We next study the probability distribution of the durations of θ-and δ-bursts over a 24h period for control and PZ-lesioned rats (Materials and methods, experimental setup). We notice that θ and δ bursts follow very different statistics. The probability density P θ of the θ-burst durations exhibits a power-law behavior, followed by an exponential cut-off (Fig 2a), where α denotes the scaling exponent of the power law. Power laws are the fingerprint of scale invariance and, depending on the context, they imply that events of any size, length or duration are likely to occur with some finite probability that is larger than expected in random or short-range correlated processes. Presence of power law indicates absence of characteristic time scales in the underlying control mechanism, which is a typical feature of physical systems at the critical point of continuous phase transition-a highly sensitive state where cooperative behavior spontaneously emerges over a range of time scales characterized by long-range correlations (presence of such correlations in the duration of consecutive δ-and θ-bursts is shown in Fig 7 and is discussed later in the manuscript). Notably, such scale-invariant power-law behavior is not influenced by lesions of the PZ neurons: in both control and PZ-lesioned group we find an exponent of α ' 2.35 (Fig 2a), indicating that lesions of the PZ do not alter the dynamical micro-architecture of θ-bursts across the 24h sleep-wake cycle. Smoothed ratio R θδ of the spectral power in the θ and δ band during a 30 min segment of 12-hour dark (lights-off) period for a control rat (top panel) and a PZ-lesioned rat (bottom panel). R θδ is calculated in non-overlapping windows w = 5 s; smoothing is performed using a 5 point moving average. θ-and δ-bursts are defined as sequences of consecutive windows w where either the power in θ or δ band is dominant.
In contrast to the power-law characteristic of θ-bursts, δ-burst durations follow a distinct behavior that is described by a Weibull distribution (Fig 2b), where λ is the characteristic time scale and β is the shape parameter. Further, we find that the distribution of δ-burst durations follows the same Weibull functional form for both control and PZ-lesioned groups, with similar values of the Weibull parameters λ and β. The functional forms established for the distributions of θ-and δ-burst durations in Eqs 1 and 2 and Fig 2, indicate a very different temporal organization of θ-and δ-bursts. A surrogate test based on randomizing the sequence of windows w in the EEG spectrogram (Fig 1a) leads to a different profile of θ-and δ-burst durations with an exponential functional form for their distributions (Fig A in S1 File), indicating that the observed temporal organization in bursting activity of brain rhythms is physiologically relevant and relates to underlying regulation.
The results demonstrate a remarkable duality of power-law scale-invariant dynamics for θbursts and Weibull dynamics with characteristic time scale for δ-bursts. Such duality of two unlikely processes appears to be a general feature of brain cortical activity of individual subjects in each group across the entire sleep-wake cycle. Moreover, this feature is robust as it remains unchanged in both control and PZ-lesioned rats. Coexistence of scale-invariant and exponential type behaviors is a hallmark of non-equilibrium systems at criticality characterized by alternating active and inactive states, and exhibiting self-organization (i.e., maintaining critical behavior without external tuning) [35,52,53]. Thus, our observations indicate an intrinsic common mechanism that underlies the temporal organization of bursting activity in both δ and θ cortical waves across the distinct physiological states of wake/arousals and sleep.
To achieve a more detailed understanding of the temporal organization of bursting activity of θ-and δ-waves during sleep and wake, and the role of PZ neurons in these dynamics, we next analyze the probability distributions of θ-and δ-burst durations separately during 12-hour dark (nighttime) and light (daytime) periods. In contrast to humans, rats are predominantly awake during the dark period and asleep during the light period. Although sleep and wake are characterized by different dominant brain rhythms with distinct dynamics, synchronization and coupling patterns across cortical areas [25,54], our analyses show that the duality of power-law for the θ-bursts and Weibull distribution for the δ-bursts durations is robust and remains independent of the dominant physiologic state (sleep or wake). Such coexistence of scale-invariant and scale-specific temporal organization in θ-and δ-bursts appears to be a basic characteristic of cortical activity during both dark and light periods (Fig 3). Importantly, the scaling exponent α of the power law as well as the Weibull parameters λ and β concurrently change comparing dark to light periods, indicating a coordinated modulation of θ and δ brain dynamics across sleep and wake. This rises the hypothesis of a coupling mechanism controlling the durations of consecutive θ-and δ-bursts (confirmed by empirical results shown in Fig 8). The observations are consistent for both control and PZ-lesioned group, demonstrating a remarkable robustness of the power-law and Weibull functional form, which remain present across the dark and light periods, and even after lesioning the PZ neurons (Fig 3).
Further, considering separately dark and light periods, we observe that the power-law scaling exponent α is higher for the dark period (Fig 3a and 3c), indicating that θ-bursts of longer duration are more likely during the light period when rats are predominantly asleep. Higher probability for longer lasting θ-bursts during light periods could be associated with the presence of longer episodes of REM sleep, where θ-wave oscillations are dominant. This leads to an average increase of spectral power in the θ band, and thus, to a higher likelihood for longer θ-bursts (smaller α) during the light period. For both control and PZ-lesioned rats, the scaling exponent α � 2.45 is higher during the dark period, and lower α � 2.25 during the light period (Fig 3a and 3c), in concurrence with change in Weibull parameters λ and β (Fig 3b and 3d). Notably, comparing the control vs the PZ-lesioned group within a given dark or light period, we find no significant difference for the distributions of θ-and δ-burst durations, indicating that the temporal organization of these fundamental brain rhythms across the sleep-wake cycle is not influenced by neuronal assemblies in the PZ area.

Robust organization of θ-and δ-bursts across time scales
Our findings show that bursts associated with θ and δ rhythms exhibit a distinct temporal organization described by specific duration distributions: power law for θ-bursts indicating absence of a characteristic time scale (i.e. scale-invariant behavior), and Weibull distribution for δ-bursts with a characteristic time scale λ (Figs 2 and 3). These findings are obtained based on a particular choice for the observational window size w and threshold Th utilized to analyze bursting dynamics (Fig 1; Materials and methods, Data analysis). We note that in our analyses, θ-and δ-bursts are defined as time periods for which the ratio R θδ = S(θ)/S(δ) is above or below the threshold Th = 1, where R θδ is evaluated over consecutive windows of w = 5 s. To demonstrate that our results are indeed independent of the particular choice of (a) Probability distributions of θ-burst durations for control (circles) and PZ-lesioned (triangles) rats over the 12 h dark lights-off period (pooled data) follow a power-law with an exponent α ctrl = 2.44 ± 0.06 and α PZ = 2.40 ± 0.06 (higher than for the 24 h sleep-wake cycle, Fig 2), where the line shows a power-law fit for the control group. (b) Probability distributions of δ-burst durations for control and PZ-lesioned rats over 12 h dark period (pooled data) follow a Weibull form, with no significant differences in the fitting parameters (β ctrl = 0.59, λ ctrl = 0.16; β PZ = 0.54, λ PZ = 0.14), where the line shows a Weibull fit for the control group. (c) Probability distributions of θ-burst durations for control and PZ-lesioned rats over the 12h lights-on period (pooled data) also follow a power-law but with smaller exponent α ctrl = 2.28 ± 0.07 and α PZ = 2.24 ± 0.08 compared to the dark period, indicating higher probability for longer durations. (d) Probability distributions of δ-burst durations for control and PZ-lesioned rats over the 12 h light period (pooled data) follow Weibull behavior for both groups with no significant differences in the fitting parameters (β ctrl = 0.54, λ ctrl = 0.12; β PZ = 0.56, λ PZ = 0.14). Lines in (c) and (d) show fits for the distributions of the control group. All durations are calculated using window size w = 5 s for the spectrogram and threshold Th = 1 on the ratio R θδ (Fig 1). Error bars are calculated for each value and where not shown are smaller than the symbol size. Error bars calculation and binning procedure are described in Materials and methods, Data analysis. The power-law exponent value for control and PZ rats do not show significant difference for both light and dark periods (t-test, p > 0.46).
parameters Th and w, we repeat the analyses for a range of parameter values. We find that universal scaling functions describe the dynamics of burst durations across the 24-hour sleep-wake cycle.
We first examine the duration distributions of θ-and δ-bursts for different threshold values Th, keeping the window size w fixed. By increasing the threshold on the ratio R θδ from Th = 1 to Th = 2, we find that the scaling exponent α characterizing the power-law distribution of θburst durations remains stable (data collapse on to a single curve, Fig 4a and 4c). The scaling behavior is followed by an exponential cut-off that, with increasing Th values, shifts to shorter burst durations d θ . Non-equilibrium critical dynamics in brain waves define sleep and wake micro-architecture For both control and PZ-lesioned rats this behavior is captured by the relation where f θ (d θ /Th −� ) is a universal scaling function, α is the power-law scaling exponent and � expresses the dependence of the cut-off on Th. This universal scaling function is confirmed by the data collapse obtained by plotting d a y P y ðdÞ versus d θ /Th −� for a range of Th values (insets in Fig 4a and 4c). In addition to the universal scaling function, such exponential cut-off is reminiscent of finite size effects observed in systems at criticality.
Similarly, we demonstrate that a unique function f δ describes the distribution of δ-burst durations that is independent of the threshold Th (Fig 4b and 4d). Since a δ-burst is defined as a time period of consecutive windows w where R θδ < Th = 1 (Fig 1b), to properly explore the behavior of the duration distribution for states with increasingly dominant δ power, we repeat the analyses for different values Th < 1. We observe that as Th decreases the probability for long δ-bursts decreases, while short δ-bursts become more likely (insets in Fig 4b and 4d). However, when distributions are rescaled by their respective mean δ-burst duration hd δ i, they all collapse onto a unique function f δ , which is the same for both control and PZ-lesioned rats (shown in main panels of Fig 4b and 4d respectively). This universal function f δ is defined by the scaling relation, where η = 1.2 for both rat groups, and is well fitted by a Weibull functional form as we find by rescaling d δ and P δ . Repeating the analyses for 12-hour dark and light periods separately, we find that the universal scaling forms in Eqs 3 and 4 consistently describe the dynamics of δ-and θ-bursts in both control and PZ-lesioned groups ( Fig B and Fig C in S1 File). Thus, our results indicate that the duality of power-law and Weibull distribution, as well as the scaling properties summarized in Eqs 3 and 4, are robust features of the bursting dynamics of θ and δ rhythms across the sleep-wake cycle, and do not depend on the particular threshold Th utilized to study the time evolution and intermittent dynamics of θ-and δ-bursts embedded in the EEG spectral power.
Moreover, we also find that the functional behavior of the distributions of θ-and δ-burst durations is to a large extend independent of the window size w, used to investigate the time course of the EEG spectral power (Fig 1a). Intuitively, larger w would tend to fail in identifying short bursts and merge them together, thus causing an increase in the probability of observing longer durations. In that regard, considering the power-law temporal organization of θ-bursts, larger window sizes w mainly influence the tail of the distribution with longer θ-bursts durations, thus leading to decrease of the scaling exponent α (insets in Fig 5a and 5c). We note that the window size effect becomes visible only for extremely large w compared to the average θburst duration hd θ i. However, when the θ-bursts distributions (curves in insets of Fig 5a and  5c) are rescaled by their respective window size w, all data collapse onto a single power law (Fig 5a and 5c), confirming the robustness of the results obtained in Fig 2a for both control and PZ-lesioned rats. This rescaling is represented by the following relation P y ðdÞ � w À 1 � f y ðd=wÞ: ð5Þ Separate analyses of 12-hour dark and light periods for different window sizes w (shown in Fig D a and c, Fig E a and c in S1 File) further confirm the robustness of the established scaleinvariant power-law form for the θ-burst durations (Fig 3a and 3c).
A similar data collapse characterizes the dependence of the δ-burst duration distribution on window size w. Generally, we observe that for increasing w the probability for long δ-bursts increases, while short δ-bursts become less likely (insets in Fig 5b and 5d). We find that the relation between w and the average duration of δ-bursts hd δ i is given by w � hd δ i ξ , with exponent ξ = 1.2 for both control and PZ-lesioned rats. When δ-burst duration distributions corresponding to different window sizes w, are rescaled by their respective mean duration hd δ i, we find that all distributions collapse onto a unique function f δ of a Weibull form (Fig 5b and 5d), Distributions of θ-and δ-burst durations for different scales of observation defined by the window size w (Fig 1). (a) Rescaled distribution P θ of θ-burst durations for control rats over a 24 h period (pooled data). Distributions obtained for different scales of observation are rescaled by the window size w and consistently show the same power-law behavior with α = 2.35, as proven by the data collapse (red line). Small deviations observed on the tail of P θ for w > 6 are due to coarse-graining effect at large-window sizes. Inset: Distributions P θ for different window sizes w (not rescaled). (b) Rescaled distribution of δ-burst durations for control rats over a 24 h period (pooled data). Distributions obtained for different window sizes w are rescaled by hd δ i ξ , where hd δ i is the mean δ-burst duration and ξ = 1.2, and collapse onto a single function that is well described by a Weibull distribution f(d; λ, β) with λ = 0.55 and β = 0.59 (black line). Inset: Distributions P δ for different window sizes (not rescaled). (c) Rescaled distribution of θ-burst durations for PZ-lesioned rats over a 24h period (pooled data). Distributions obtained for different window sizes w are rescaled by the corresponding window size, and consistently show the same power-law behavior with α = 2.35, as proven by the data collapse (red line). Small deviation observed on the tail of P θ for w > 6 are due to large-window effects. Inset: Distributions P θ for different window sizes (not rescaled). (d) Rescaled distribution of δ-burst durations for PZ-lesioned rats over a 24 h period (pooled data). Distributions are rescaled by hd δ i ξ , with ξ = 1.2, and collapse onto a single function following a Weibull behavior f(d; λ, β) (black line) with λ = 0.44 and β = 0.54. Inset: Distributions P δ for different window sizes (not rescaled). Results in all panels are obtained for fixed threshold Th = 1 on the ratio R θδ (Fig 1). Results are consistent when considering separately light and dark periods (Fig D and Fig E in S1 File).
Non-equilibrium critical dynamics in brain waves define sleep and wake micro-architecture expressed by the following scaling relation As in the case for θ-burst dynamics, we find that the temporal organization of δ-bursts is robust and characterized by a scaling function (Eq 6), which is universal for the control and PZ-lesioned groups, and remains stable during light and dark periods (Fig D b and d, Fig E b and d in S1 File).
The existence of universal scaling functions (Eqs 3-6) not only demonstrates that duration distributions are independent of the specific set of parameters used to identify θ-and δ-bursts, but also constitutes a striking evidence of scale invariance, a property associated with systems operating at criticality.

Self-similar micro-architecture of active and quiet states in the sleep-wake cycle
Our investigations reveal a scale-invariant power-law structure for the durations of θ-bursts and a homogeneous functional form of Weibull type with a characteristic time scale for δburst durations, that are remarkably robust during light and dark periods across the sleepwake cycle, and are consistent for both the control and PZ-lesioned group (Figs 2 and 3). The coexistence of these two distinct types of dynamics draws a strong parallel with far-from-equilibrium physical phenomena that are characterized by bursting dynamics and abrupt transitions between active and quiet states, such as avalanches and earthquakes [15,35,[55][56][57]. For instance, the energy released during avalanches/earthquakes (active states) is also distributed according to a power law, while the distribution of time intervals between consecutive avalanches/earthquakes (quiet states) is described by a generalized Gamma distribution with a characteristic time scale (exponential tail). Gamma is a universal scaling function that is independent of spatial scales and minimum magnitude thresholds, and is consistently observed for a broad range of conditions despite the large variability associated with phenomena such as earthquakes and avalanches [55,[58][59][60].
In the context of sleep dynamics, wake and brief arousals during sleep can be considered as active states that, in rodents, are characterized by bursts in θ rhythms. Thus, we focus on the organization of θ rhythms in time. Specifically, we investigate the relationship between the duration of θ-bursts and their occurrence in time (Fig 6), and we hypothesize that a selfsimilar structure, invariant across time scales, may also characterize the occurrence of θbursts. In analogy with non-equilibrium phenomena, where quiet times are associated with the magnitude of active states, presence of a self-similar structure between occurrence and duration of θ-bursts (active states) would provide additional evidence for criticality in sleep micro-architecture.
To this end, we consider the time sequence of θ-bursts, and we investigate the statistical features of the quiet times Δt separating consecutive bursts, taking into account the duration d θ of each θ-burst (Fig 6a). Since θ-bursts vary in duration, we impose a threshold D 0 representing the time scale of analysis, and we define quiet time Δt i as the period from the end of θ i -burst to the beginning θ i+1 -burst. Thus, the statistical characteristics of Δt i depend on the threshold value D 0 . We obtain the probability distribution P(Δt; D 0 ) of quiet times Δt i for different values of D 0 (insets in Fig 6c and 6d). With increasing threshold (scale of observation) D 0 , the probability of longer Δt i increases, while the probability of short quiet times decreases, leading to different curves for the distributions P(Δt; D 0 ).
Visual inspection of the complex profile formed by the time sequence of θ-bursts and their respective durations shows an apparent similarity when comparing short segments of the profile with the entire sequence above a given threshold D 0 (Fig 6b). Indeed, the transformation illustrated in Fig 6b resembles a renormalization-group transformation, analogous to those studied in the context of critical phenomena. In particular, the picture in the bottom panel of (Fig 6b) can be seen as a result of contracting the time axis by a factor a with respect to the middle panel, and successively decimating the number of events by the same factor. To demonstrate statistical self-similarity in the sequence of θ-bursts, we systematically analyze the functional form of the probability distributions P(Δt; D 0 ) for different thresholds D 0 by rescaling each distribution on the average quiet time hDti D 0 . Remarkably, we find that all distribution curves collapse onto a single functional form G (Fig 6c and 6d), defined by the following scaling relation The scaling relation in Eq 7 represents a quantitative, mathematical expression of the statistical self-similarity in the profile formed by the quiet times and θ-burst durations shown in Fig  6b. We find that the functional form G is well approximated by the generalized Gamma distribution GðDt=hDti; b; n; pÞ ¼ ðp=b n ÞðDt=hDtiÞ nÀ 1 e À ðDt=hDtibÞ p =Gðn=pÞ [61], where in our analysis Δt/hΔti is a dimensionless quiet time. The Gamma functional form is homogeneous [62], i.e. rescaling the variable leaves the functional form unchanged. The existence of such scaling function indicates that the distribution of quiet times between consecutive θ-bursts is independent on the scale of observation D 0 . In the limit of D 0 = 0, the quiet time distribution P(Δt; D 0 ) coincides with the distribution of δ-burst durations P δ (Eq 2 and Figs 2b, 3b and 3d)-a Weibull functional form that belongs to the same class of homogeneous functions as the generalized Gamma.
Our data analysis shows that the scaling relation in Eq 7 and the associated Gamma functional form for the quiet times are robust: (i) we find it across the 24-hour sleep-wake cycle ( Fig  6) as well as separately during light and dark periods (Fig F in S1 File), and (ii) it does not significantly change with lesion of the sleep-promoting neurons in the PZ brain area (compare Fig 6c and Fig 6d). Further, the presence of self-similar structure in quiet time indicates specific temporal order in the occurrence of θ-bursts. To explicitly verify this, we randomly reshuffle the data sequence of θ-burst durations, while preserving the δ-burst durations corresponding to quiet times at D 0 = 0, and we perform the analysis on the reshuffled sequence to obtain quiet time distributions P rand (Δt; D 0 ) for different thresholds D 0 . After rescaling the distributions P rand (Δt; D 0 ) by the average quiet time hDti D 0 , their curves collapse onto an exponential distribution (dashed lines in Fig 6c and 6d)-a hallmark of temporal independence between consecutive events [63]. This clearly demonstrates that temporal correlations are intimately related to the existence of universal non-exponential scaling functions (Eqs 2 and 7) [63,64].
Notably, a similar temporal organization characterized by coexistence of power law and generalized Gamma distribution has been reported for active states and quiet times between them in a range of non-equilibrium systems self-tuning at criticality [15,53,55,58]. Thus, our findings are a strong evidence in support of the hypothesis that bursting activity of fundamental brain rhythms and the associated sleep micro-architecture exhibit critical non-equilibrium behavior.

Long-range scale-invariant correlations in θ and δ bursts
Physical systems at criticality exhibit high sensitivity to interactions among components [32,65]. This leads to the emergence of collective cooperative behavior, where interactions span the entire system across space and time scales [32], leading to long-range correlations. Indeed, scaling features in such systems often arise in conjunction with long-range spatio-temporal correlations of power-law (scale-invariant) type, as observed at the critical point of continuous phase transitions [32]. Notably, physiological systems under neuroautonomic regulation also exhibit dynamics characterized by long-range power-law correlations-a scale-invariant structure that undergoes a phase transition with transitions from sleep to wake [66][67][68][69], with circadian rhythms [70][71][72][73] and under clinical conditions [74][75][76]. Further, the randomization procedure in the previous subsection (Fig 6) clearly demonstrates that a self-similar structure in quiet times, characterized by a Gamma scaling function (Eq 7), can arise only in the presence of a certain temporal order in θ-bursts. Thus, we next perform correlation analysis to quantify long-range features in the temporal organization of δ-and θ-burst durations.
To this end, we utilize the detrended fluctuation analysis (DFA)-a random walk based method, specially tailored to quantify long-range power-law correlations embedded in nonstationary signals with various polynomial trends and bursting dynamics [77,78]. The DFA method is based on evaluation of the root mean square (r.m.s.) fluctuation function F(n) (Materials and methods, Data analysis), where n is the scale of analysis expressed in number of consecutive bursts (Fig 7). A scaling relationship of the form FðnÞ / n a d indicates presence of long-range power-law correlations in the time series of burst durations. Correlation exponent α d 2 [0, 0.5) indicates anti-correlations (where short burst durations tend to be followed by longer burst durations), while α d 2 (0.5, 1] indicates positive persistent correlations (long bursts tend to be followed by longer bursts)-a scale-invariant behavior that is consistent over several decades of time scale n; α d = 0.5 corresponds to a random walk and absence of correlations.
We performed DFA on sequences of θ-and δ-burst durations separately, distinguishing between dark and light periods as well as between control and PZ-lesioned rats (Fig 7). We find that both θ-and δ-bursts exhibit long-range power-law correlations, with an exponent α d = 0.60 ± 0.02 and α d = 0.65 ± 0.02, respectively. These exponents are robust and do not Non-equilibrium critical dynamics in brain waves define sleep and wake micro-architecture change during dark or light periods, in line with our observations for the duration distributions (Fig 3). Moreover, our analyses indicate that PZ lesions do not affect the nature and strength of temporal correlations, and consistently shows that θ-and δ-bursts are long-range correlated across the sleep-wake cycle (Fig 7).

Anti-correlated coupling between δ-and θ-bursts
Majority of physical and biological systems at equilibrium (homoeostasis) are controlled by mechanisms that either lead to dynamics with specific space or time scales characterized by exponential behaviors or to scale-invariant dynamics without characteristic scales following power laws. Non-equilibrium systems self-organizing at criticality are unique in the sense that they combine two distinct processes-a scale-invariant process related to the dynamics of active states and an exponential process related to quiet states-both of which emerge out of a single regulatory mechanism [13]. In that context, our findings of (i) power-law distribution for theta-burst durations (scale-invariant dynamics of active states) in coexistence with Weibull distribution for δ-burst durations (characteristic time scale for the dynamics of the quiet states) shown in Figs 2 and 3, (ii) universal Gamma distribution characterizing the temporal organization of quiet times across a range of scales (Fig 6), and (iii) long-range power-law correlations in the time sequence of θ-and δ-burst durations (Fig 7)-all these characteristics are typical for systems at criticality-indicate a common sleep regulatory mechanism for the bursting activity of both θ and δ rhythms and associated sleep micro-architecture. Presence of such common mechanism would imply coupling between θ-and δ-bursts. Indeed balanced excitation and inhibition in neuronal networks is essential to maintain critical-state dynamics [79][80][81]. Further, the concurrent change we find in both power-law and Weibull distribution parameters with transition from dark to light periods (Fig 3) is an additional indication of possible coupling between θ-and δ-bursts. Thus, we ask whether there is a cross-correlation between the durations of consecutive δ-and θ-bursts.
To further understand the temporal organization of bursting dynamics in relation to neuronal integrity in the PZ, we next investigate the coupling between consecutive δ-and θ-bursts, and the role of such coupling in the emergent scaling behavior of duration distributions in control and PZ-lesioned rats.
We first focus on the relationship between ranks of consecutive δ-and θ-burst durations, d δ and d θ . We rank burst durations in ascending order, with the shortest duration corresponding to the smallest rank, and examine the scatter plots between the ranks of consecutive d δ and d θ (Fig 8a and 8b). We find that δ-bursts of high ranks (i.e. long durations) tend to be followed by θ-bursts of low ranks (i.e. short durations). This anti-correlated coupling is consistently present in both control (Fig 8a) and PZ-lesioned rats (Fig 8b), and appears to be a basic characteristic of dynamics as it is observed throughout the entire sleep-wake cycle in both dark and light periods.
To quantify the coupling between consecutive δ-and θ-burst durations we utilize Spearman's correlation coefficient, which assesses monotonic relationships between two variables (Materials and methods, Data analysis). The Spearman's cross-correlation is positive when observations of two variables have similar ranks, and negative if observations of two variables have opposite ranks. Our analyses show that the cross-correlation coefficient calculated for consecutive δ-and θ-burst durations is always (24h, dark, light period) significantly negative (Fig 8c)-a clear sign of anti-correlated coupling. This is verified by a surrogate test where the sequence of consecutive δ-and θ-burst durations is randomized (Fig 8c; Materials and methods, Data analysis). We find that the δ-θ anti-correlated coupling is more pronounced for PZ-lesioned rats, and this is consistently observed in both light and dark periods. Comparing light vs dark period, our results show an increase in the anti-correlated coupling during the light period within each group (Fig 8c).
The presence of anti-correlated coupling between consecutive δ-and θ-bursts is further supported by conditional probability analysis (Fig G in S1 File). Specifically, we ask how the conditional probability distribution Pðd d i jd y iÀ 1 > d � Þ of δ-burst durations d d i depends on the length of preceding θ-burst duration d y iÀ 1 above a given threshold d � -i.e., we consider the probability distribution of the subset of δ-bursts which follow θ-burst durations longer than d � . Note that when no condition is imposed on d y iÀ 1 , i.e. d � = 0, the conditional probability is equivalent to the Weibull distribution of δ-burst durations, namely Pðd d i jd y iÀ 1 > 0Þ ¼ Pðd d i Þ (Figs 2 and 3). On the other hand, for threshold d � > 0, Pðd d i jd y iÀ 1 > d � Þ 6 ¼ Pðd d i Þ implies that d d i and d y iÀ 1 are not independent of each other. As discussed in further details in the Supplementary Information, we find that increasing the conditional threshold for d y iÀ 1 the probability

Fig 8. Coupling between δ-and θ-burst durations indicates a common mechanism regulating the activity of these rhythms in relation to sleep micro-architecture.
Scatter plots and rank correlation analysis demonstrate coupling between consecutive δ-and θ-burst durations. (a) Scatter plot of δ-burst ranks vs following θ-burst ranks in the 24h period for control rats. Each dot represents a pair formed by a δ-burst and the following θ-burst, with burst durations separately ranked among the δ-bursts and the θ-bursts (longest duration corresponding to highest rank). (b) Scatter plot of δ-burst ranks vs following θburst ranks in the 24h period for PZ lesioned rats. For each rat group, ranks are calculated separately for each rat and then plotted together. (c) Average Spearman's cross-correlation coefficient for control and PZ-lesioned rats in dark, light and 24h periods. Anti-correlations between consecutive θ-and δ-bursts are stronger during light than during dark periods in each of the two rat groups. Comparing dark vs light periods, the Student's t-test gives p = 0.651 for control rats and p = 0.461 for PZlesioned rats. Importantly, PZ-lesioned rats generally exhibit stronger anti-correlations than the control group, in particular during dark periods, where the strength of anti-correlated coupling increases with � 30% compared to control rats (control vs PZ t-test: 24h, p = 0.158; dark, p = 0.121; light, p = 0.064). All correlation coefficients calculated in both groups are significantly different from the corresponding values obtained in the surrogates (red bars) after randomly reshuffling the original order of θand δ-bursts (t-test: p < 0.001). All durations are calculated using a window w = 5 s and threshold Th = 1 on the ratio R θδ (as in Fig 1). This finding of anti-correlated coupling between θ-and δ-bursts durations is further supported by an independent analysis based on conditional probability (Fig G in S1 File). https://doi.org/10.1371/journal.pcbi.1007268.g008 Non-equilibrium critical dynamics in brain waves define sleep and wake micro-architecture for longer d d i significantly drops, while the probability for shorter d d i significantly increases (Fig G insets in S1 File). This dependence is observed during light and dark periods for both control and PZ-lesioned groups, and clearly confirms the anti-correlated coupling between the durations of consecutive δ-and θ-bursts.

Phenomenological model of coupling and criticality in θ-and δ-bursts dynamics
We next test whether the established anti-correlated coupling between consecutive δ-and θburst durations is essential for the emergent duality of power-law and Weibull distribution as a basic characteristic of systems at criticality. We develop a phenomenological model (Fig 9) based on anti-correlated pairing of θ and δ durations randomly drawn from the empirical distributions of the δ-and θ-burst durations established in our study (Fig 2). The model allows to control the degree of anti-correlations between consecutive θ-and δ-burst durations, and thus to examine if and how anti-correlations affect the emerging power-law and Weibull distributions of burst durations, and the related scale-invariant temporal structure of θ-and δ-bursts.
The basic steps to generate sequences of alternating θ-and δ-burst durations with the desired degree of correlations are schematically outlined in Fig 9 (Materials and methods, model of anti-correlated burst coupling.). Specifically: (1) we randomly draw durations d θ and d δ from their respective power-law and Weibull distributions obtained from our empirical data analyses (Fig 9a); (2) d θ and d δ are next separately sorted in ascending order, from shortest to longest, and are assigned a unique rank (Fig 9b); (3) the durations d θ and d δ are then paired based on their ranks with a certain degree of anti-correlation (defined by and monotonically depending on a single parameter) to generate an artificial time series of alternating θ-and δburst durations (Fig 9b); (4) the obtained artificial time series is then binarized in windows w corresponding to the window size used in our EEG spectral power analysis of the original data (Fig 1a), where '+' is assigned for θ-bursts and '−' for δ-bursts. Finally, the binary series is coarse-grained at larger time scale Δ (Fig 9c).
We utilize this model to test our hypothesis that coupling between δ-and θ-bursts is essential for the emergent duality of power-law and Weibull behavior across time scales as a hallmark of system at criticality. We test to what extent δ-θ coupling strength plays role in the emergent scale-invariant organization of sleep micro-architecture. Our simulations show that the generated distributions of δ-and θ-burst durations depend on the degree of anti-correlation introduced in the model. When the Spearman's cross-correlation coefficient of burst durations generated by the model corresponds to the empirical values found in real data, the distributions obtained from the model approximate the empirical distributions (Fig 10), and scale-invariant temporal organization in burst durations emerges over a range of coarse-graining scales Δ. In contrast, absence of anti-correlated coupling in our model (i.e. random pairing of δ and θ durations in the generated time series) leads to exponential distribution for both δand θ-bursts, significantly different from real data (Fig 10).

Discussion
We show that transient bursts in both the θ and δ cortical rhythms continuously occur during the sleep-wake cycle and exhibit a complex temporal organization that is invariant across a range of time scales, from seconds to minutes, and underlies sleep micro-architecture. We discover a remarkable duality of scale-invariant power-law distribution for θ-burst durations (active states) and a Weibull distribution with a exponential characteristic time scale for the δburst durations (quiet states) (Fig 2)-a behavior which is typically observed in non-equilibrium systems self-organizing at criticality, where a quiescent phase with exponential dynamics coexists with active events following power-law distributed sizes and durations [15,34,35,57]. Further, we identify presence of coupling between θ-and δ-bursts dynamics which is characterized by significant anti-correlation in consecutive θ-and δ-burst durations (Fig 8), and we demonstrate through both empirical and modeling approaches that this anti-correlated coupling is essential part of the mechanism responsible for the emergent duality of power-law and Weibull behavior across time scales (Figs 9 and 10). Importantly, we find that sequences of consecutive θ-or δ-burst durations are long-range power-law correlated, indicating a scaleinvariant organization in the temporal order of burst durations and a unique underlying process with persistent 'memory' spanning over a wide range of scales that statistically couples the Non-equilibrium critical dynamics in brain waves define sleep and wake micro-architecture  θ-and δ-bursts, Fig 9). (a) Distributions P θ (d) of θ-burst durations for: (i) 24 h control rats data (red diamonds), (ii) model-generated time series of θ-and δ-bursts durations with anti-correlations (green circles), and (iii) model-generated time series with random pairing of θ-and δ-bursts durations (magenta dashed line). Inset shows results from same analysis on P θ (d) for the group of PZ-lesioned rats. (b) Distribution P δ (d) of δ-burst durations for: (i) 24h control rats data (blue diamonds), (ii) model-generated time series with anti-correlations (green circles), and (iii) model-generated time series with random pairing of θ-and δ-bursts durations (magenta dashed line). Inset shows results from same analysis on P δ (d) for the group of PZ-lesioned rats. In both (a) and (b), durations are in units of Δ, which is the window size used to coarse grain the sequences of θ-and δbursts durations. The distributions obtained from the model using anti-correlated d θ and d δ pairing (green circles) closely match the duration distributions for the original data (diamonds)-power law for P θ (d) and Weibull for P δ (d) -for both control and PZ-lesioned rats. In contrast, a random pairing of d θ and d δ produces duration distributions following the Poisson functional form (magenta dashed lines) that significantly deviates from the original data. duration of a given burst with the durations of hundreds of following bursts (Fig 7). Presence of complex temporal organization and coupling in cortical rhythms is also manifested through the self-similar structure we uncover in the quiet times separating consecutive θ-busts (active events) above a given duration (Fig 6)-described by a homogeneous generalized Gamma distribution [61]. This self-similar structure links, across time scales, the duration of a given θburst with the time of its occurrence.
Our empirical analyses show that the characteristics of θ-and δ-bursts dynamics do not depend on the scale of observation or on the threshold used to separate θ-from δ-bursts, and remain continuously present during dark and light periods (Figs 4, 5 and Figs B, C, D, E in S1 File), under different dominant physiologic states, and both in control rats and rats where all PZ neurons are lesioned. Thus, our findings indicate that the discovered scale-invariant organization in the bursting activity of θ and δ cortical rhythms is independent of behavioral factors across the sleep-wake cycle, and is not affected by the loss of neuronal inputs from the PZ, a major component of the sleep-promoting circuitry. Further, the presence of multiple scaleinvariant characteristics related to distributions, correlations, coupling and timing of bursting events is a strong evidence of critical self-organization in θ and δ cortical rhythms-a signature of collective behavior over a range of time scales that emerges from neuronal interactions across brain areas which is essential to maintain system's susceptibility and flexibility for abrupt sleep-stage and arousal transitions. The reported self-organization at the integrated cortical level, with dynamics spontaneously driven at criticality by active and quiet states, indicates that a non-equilibrium mechanism regulates sleep micro-architecture on time scales from seconds to minutes-a behavior which is in contrast to the homeostatic (equilibrium) regulation that controls consolidated sleep and wake, ultradian and circadian rhythms at large horizons of hours and days.
Our observations provide new insights about the role of PZ. Previous studies have demonstrated that PZ lesions cause a significant decrease of NREM sleep during the light period, without affecting the distribution of EEG spectral power among different brain rhythms [6,51]. Our results suggest that such stability of the spectral power distribution results from the robust temporal organization and cross-talk of δ and θ rhythms, which seems not to be significantly altered by lesion of the PZ neurons. Furthermore, our findings show that the decrease in NREM sleep is not associated with significant changes in the dynamic characteristics and temporal organization of δ rhythms.
Sleep and wake are under control of complex regulatory circuitry involving multiple neuronal assemblies and specialized brain areas. GABAergic neurons from the PZ, together with other NREM sleep-promoting neurons, inhibit wake-promoting populations and arousal pathway which are crucial for cortical activation and wakefulness. During the past few years multiple NREM sleep promoting neuronal populations have been described [82][83][84][85]. Understanding the specific role of PZ sleep-promoting neurons and how they interact with other sleep-and wake-promoting brain areas to influence cortical activity is crucial to develop an integrated picture of sleep-wake control. For example, it is likely that other neuronal assemblies compensate for the loss of PZ neurons, e.g. the sleep-promoting neurons located in the ventrolateral preoptic area (VLPO) [86]. Our findings of scale-invariant features of θ-and δbursts dynamics which remain present after lesioning the PZ neurons, while at the same time total NREM sleep declines [6,51], are a strong indication that the function of PZ neurons relates to sleep initiation, and may not be actively involved in regulating the dynamics once a sleep episode is initiated-i.e., lesioning the PZ reduces sleep initiation, and thus leading to decline in total NREM sleep; however, once sleep is initiated the dynamics and micro-architecture of θ-and δ-bursts are not significantly altered, raising the hypothesis that another sleeppromoting center (possibly the VLPO) may be involved in the dynamic aspects of sleep regulation. Moreover, our empirical findings of increased anti-correlation between durations of consecutive θ-and δ-bursts under PZ lesion (Fig 8) and model simulations demonstrating that such θ-δ coupling is essential for the emerging scale-invariant temporal organization in these cortical rhythms (Figs 9 and 10), indicate that PZ neurons may have dual role for both sleep and arousal/brief wake activation.
Importantly, the uncovered complex dynamics in θ and δ cortical rhythms do not depend on the time scales of observation, and share striking similarities with the dynamical characteristics of natural phenomena exhibiting non-equilibrium behavior of active and quiet states, and self-organization at criticality. Our analyses and findings are based on the interpretation of θ-bursts as active states and δ-bursts as quiet states, and on the hypothesis that the observed scale-invariant temporal organization of θ and δ rhythms emerges out of a common sleep regulatory mechanism, in analogy to the cooperative mechanisms underlying the dynamics of non-equilibrium systems at criticality. This approach is consistent with the basic neurophysiological understanding of δ rhythm as the quiet state cortical default mode-indeed, lesion and transection experiments have reported that interruption of sensory inputs to the cortex results in a cortical EEG similar to that in NREM sleep [87][88][89]. In contrast, oscillations in the θ band are associated with REM, arousals and wakefulness [10,11]. Due to the respective amount of wakefulness and REM sleep in our data [50], most of the analyzed θ-bursts are likely associated with arousals and wake. Further, during wakefulness cortical θ activations are desynchronized with a wide range of burst durations, and correspondingly our analyses show that collective cortical neuronal activity follows robust universal scaling laws: power-law distribution for burst durations, long-range power-law correlations in the sequence of burst durations, and a self-similar structure in quiet times between bursts, all reminiscent of the dynamic characteristics and temporal organization of active states found in avalanche and earthquake dynamics [59,60].
The presented findings provide a general picture unifying previous empirical observations of criticality in spontaneous brain dynamics at different levels-from networks of dissociated cortical neurons [90] and local field potentials (LFP) in cortex slice cultures [45], human EEG and fMRI resting state dynamics [48,49,91,92], awake monkeys [93] and resting magnetoencephalography of the human brain [46], to the dynamics of sleep-stage and arousal transitions across species [13,26,27,31,94]-where either distributions or temporal correlations of active events have been studied and discussed in the context of self-organization at criticality. Several models based on statistical physics have also shown that the functional properties of the healthy brain resemble those of systems at criticality, and that altered physiologic states may correspond to super-critical or sub-critical states [95,96]. In particular, models including disorder and frustration [47,97]-essential features of spin glasses as well as of multi-layer cortical networks where neurons may be frustrated by receiving both excitatory and inhibitory inputs-provide critical exponents in good agreement with those observed experimentally depending on the balance between excitatory and inhibitory neurons and their distributions among cortical layers [81]. Crucially, here we show that important physiologic functions may benefit from underlying critical dynamics, and demonstrate presence of the full spectrum of scaling characteristics typical for non-equilibrium systems self-organizing at criticality. Furthermore, we link our observations to the collective behavior of a key sleep-promoting neuronal population leading to emerging cortical rhythms in relation to physiological alternation of sleep and wake. We find that the power-law scaling exponent of the distribution for θ-burst durations is α ' 2.35, a robust temporal organization in bursting activity across the sleep-wake cycle in both control and PZ lesioned rats. Notably, this scaling exponent is close to the power-law exponent α ' 2 for the distribution of neuronal avalanches in cortex slice cultures [45] and dissociated cortical neurons [90], to α ' 2.2 reported for arousal/wake episodes in coarse-grained sleep-stage recordings in humans [13,31] and other species [26,27,30,94].
In summary, our findings of scaling features for a full spectrum of dynamic characteristics in the bursting activity of cortical θ and δ rhythms strongly support the hypothesis of an underlying critical dynamics for sleep regulation [26,65]. In systems far from equilibrium, emerging bursting activity described by power laws and exhibiting long-range spatio-temporal correlations has been proposed as an indication of self-organized criticality (SOC) [16,17,45]. In this context, bursts do not have a characteristic duration, and short as well as long bursts are expression of the same underlying dynamics [17]. Consecutive bursts are separated by quiescence periods whose distribution depends on the details of the system and generally exhibit an exponential tail [55,58,98,99], and is an exponential for the paradigmatic sandpile model of SOC [17,35]. Thus, in systems exhibiting self-organized criticality power-law and exponential dynamics for active and quiet states coexists, and emerge out of the same regulatory mechanism. The robust duality of power-law (scale invariant) and Weibull (exponential tail) distribution for the bursting dynamics of θ and δ rhythms is closely reminiscent of this scenario, where scale-free θ-bursts in cortical activity can be seen as avalanches or earthquakes [17,55], while δ-bursts can be interpreted as the quiet periods between active states. Notably, the Weibull functional form for the distribution of quiet periods represents an extreme (minimal) event statistics. In the proposed here criticality framework of cortical rhythm dynamics and sleep micro-architecture, one can interpret the Weibull distribution as realization of shortest quiet periods between spontaneous initiation of active states (θ-bursts). Since such activations are spontaneously initiated at multiple brain locations, what one measures at the EEG probe is the closest (fastest arriving) activation. Following the analogy with SOC systems, we further demonstrated that the organization (occurrence in time) of θ-bursts is coupled with their durations, forming a scale-invariant structure for the quiet times between consecutive θ-burst above a given duration described by a universal Gamma distribution (also observed in earthquake dynamics [55]).
Overall, the combined empirical observations and modeling simulations reported here lay the foundation for a new paradigm for the investigation of sleep dynamics and sleep-stage transitions mechanisms, considering sleep micro-architecture as result of a non-equilibrium process and self-organization among neuronal assemblies to maintain a critical state-a behavior which is in stark contrast to the traditional homeostasis paradigm of sleep regulation at large time scales. Within this criticality paradigm arousals/brief wake and sleep-stage dynamics emerge out of a common regulatory mechanism, where arousals play an essential part in maintaining a non-equilibrium behavior at criticality, as evidenced by our observations of coupling between consecutive θ-and δ-bursts durations. Systems at criticality exhibit high susceptibility and sensitivity to interactions, leading to cooperative behaviors over a range of space/time scales, and thus, maintaining a critical state through self-organization is important for system's flexibility and for generating spontaneous transitions [100,101]. Such transitions can not occur intrinsically in a homeostatic (equilibrium) systems. In the context of sleep micro-architecture, the proposed criticality-based paradigm may provide new insights on the origin and mechanisms underlying the dynamics of sleep-stage and arousal transitions, and offers a unifying picture of sleep and wake.

Experimental setup
Animals. Pathogen-free adult male Sprague Dawley rats (Harlan; 275-300 g; n = 20) were used in this study. Care of these animals in the experiments met National Institutes of Health standards, as set forth in the Guide for the Care and Use of Laboratory Animals, and all protocols were approved by the Beth Israel Deaconess Medical Center and Harvard Medical School Institutional Animal Care and Use Committees.
Surgery. For details on brain injection and implantation for polysomnographic recording in rats, see Anaclet et al., 2012 [50]. To perform cell-specific lesions, 10 rats received brain microinjections of 0.1% anti-orexin-B IgG saporin (OX-SAP; 130-330 nl; Advanced Targeting Systems) within the PZ [AP, -10.3 mm; L, ±2.1 mm; DV, -6.6 mm, as per the rat atlas of Paxinos and Watson (2005)]. The rat control group included five rats that received saline brain injections and five rats without brain injection. After the brain injections, the rats were implanted with polysomnographic electrodes. Four electroencephalogram (EEG) screw electrodes (Plastics One, stock # E363/20/4.8/SP) were implanted into the skull, in the frontal (2) and parietal bones (2) of each side. Two flexible electromyogram (EMG) wire electrodes (Plastics One, stock # E363/76) were placed in the neck muscles. After EEG/EMG implantation, rats were single housed until the end of the experiments.
Sleep-wake recording. Ten days after surgery, the rats were moved in an insulated soundproofed recording chamber maintained at an ambient temperature of 22 ± 1 C and on a 12-h light/dark cycle (lights-on at 7 a.m.) with food and water available ad libitum. They were connected via flexible recording cables and a commutator (Plastics One, stock # 363-363) to an analog amplifier (A-M Systems, Model 3500) and computer, with an analog-to-digital converter card and running Vital Recorder (KISSEI COMTEC CO., LTD., Japan). After 3-5 day habituation period, EEG/EMG were recorded for 48 h, beginning at 7:00 P.M. Cortical EEG and EMG signals were amplified and digitalized with a resolution of 256 Hz. Recordings used in this study have already been subjected to sleep-wake analysis and published [50].
Histology. At the end of the experiments, rats were perfused under deep anesthesia (200 mg per kg of chloral hydrate) with 50-ml saline, followed by 200-ml of neutral phosphate-buffered formalin (4%, vol/vol, Fischer Scientific). After perfusion, the brains were removed, postfixed in neutral phosphate-buffered formalin for 2 hr, equilibrated in PBS containing sodium azide (0.02%) (PBS-azide) and sucrose (20%) for at least 1 day, and then sectioned at 40 μm on a freezing microtome into four series. For verification of OX-SAP-induced brain lesions, one series of tissue was processed for Nissl staining as done previously (Lu et al., 2000). PZ lesion for each rat has previously been published (see Fig. 2 in [50]).

Data recording and analysis
Ten control rats and ten rats with PZ lesion were used for this study. Cortical EEG signals were recorded continuously for 48 h from the left and right hemisphere with a resolution of 256 Hz. Two electrodes, one frontal and one parietal, were placed on each hemisphere. Cortical EEG is the differential potential between a frontal electrode (for δ frequencies) and a parietal electrode (for θ frequencies). The parietal electrode picks up θ rhythm from the hippocampus during REM sleep; both the frontal and parietal electrodes pick up θ rhythm during wakefulness. Hippocampal θ rhythm can be also present during wakefulness, mainly during cognitive wakefulness. The signal analyzed in this study is the difference between frontal and parietal EEG electrode potentials (frontal − parietal EEG) from one hemisphere (ipsilateral).

Data pre-processing
EEG recordings were first normalized to zero mean, μ = 0, and unit standard deviation, σ = 1. For each rat, EEG signals were visually inspected and noisy segments were discarded. Only clean segments were finally included in the analysis. Data filtering. Data were bandpass filtered in the range 0.5 − 25 Hz using a FIR (Finite Impulse Response) filter designed in Matlab.

Data analysis
Spectral analysis. The clean EEG signal is divided in N non-overlapping windows of size w and the spectral power in the δ band (0.5 − 4 Hz), S δ , and in the θ band (4 − 8 Hz), S θ , is estimated in each window using Welch's method [102]. The analysis is performed for several values of the window size w, from 2 s to 10 s. Results are generally independent of w, as shown in Fig 5, as well as Fig D and Fig E in S1 File, and extensively discussed in the main text.
θ-and δ-burst detection and definition. The ratio R θδ = S θ /S δ between θ and δ power is calculated in each window k, with k = 1, 2, . . ., N, and a time series R θδ (k) is obtained. Given a threshold Th � 1, a θ-burst is defined as a sequence of n consecutive windows where R θδ > Th, while a δ-burst consists in a sequence of n consecutive windows where R θδ < 1/Th (Fig 1). The duration of a burst is given by d = n � w. Durations of θ (δ) bursts are denoted by d θ (d δ ). The threshold Th is set equal to 1 throughout the analysis. Results are independent of Th, as shown in Surrogate test for θ-and δ-burst duration distributions (Fig A in S1 File). For each rat, the time series R θδ (k) is randomly reshuffled to obtain a surrogate R � yd ðk 0 Þ. Surrogate θ-and δbursts durations are then calculated from R � yd ðk 0 Þ following the procedure illustrated in the previous paragraph. The corresponding θ-and δ-burst duration distributions are shown in the where p = P(D)dD is the probability to observe a duration D in the range [D, D + dD] and N is the total number of bursts and dD is the corresponding bin size.
Spearman's correlation. Given to variables X and Y, the Spearman's correlation coefficient is defined as where rg X and rg Y are the tied rankings of X and Y [104], respectively, s rg X and s rg Y their standard deviations, and cov(rg X , rg Y ) indicates the covariance between rg X and rg Y . Surrogate test for correlations between consecutive θ-and δ-burst duration. To test significance of correlations between consecutive θ-and δ-burst durations, a surrogate sequence of burst durations is generated for each rat by randomly reshuffling the original order of θ-and δ-bursts. The Spearman's correlation coefficient ρ s between consecutive θand δ-bursts is calculated for each surrogate. The average Spearman's correlation coefficient obtained from all surrogates is then compared with the average correlation coefficient calculated from the original sequences of burst durations via t − test (Section Results, Fig 8). Correlation coefficients for surrogate data for both control and PZ-lesioned rats during dark, light and 24h are all with value |ρ s | < 10 −3 .
Detrented Fluctuations Analysis (DFA). The DFA is a method based on random walk [105]. It improves the classical fluctuation analysis (FA) for non-stationary signals where embedded polynomial trends mask the intrinsic correlation properties in the fluctuations [105]. The performance of DFA for signals with different types of non-stationarities and artifacts has been extensively studied and compared to other methods of correlation analysis [78,[106][107][108]. The DFA method is briefly described by the following steps [105]: 1. A given signal u i (i = 1, . . ., N, where N is the length of the signal) is integrated to obtain yðkÞ � P k i¼1 ½uðiÞ À hui�, where hui is the mean of u i . 2. The integrated signal y(k) is divided into boxes of equal length n.
3. In each box of length n we fit y(k) using a first order polynomial function which represents the trend in that box. The y coordinate of the fit curve in each box is denoted by y n (k).
6. Repeat the above computation over a broad range of box lengths n, where n represents a specific space or time scale, to obtain a functional relationship between F(n) and n.
For a power-law correlated time series, the average r.m.s. fluctuation function F(n) and the box size n are connected by a power-law relation, that is FðnÞ � n a d . The exponent α d is a parameter which quantifies the long-range power-law correlation properties of the signal. Values of α d < 0.5 indicate the presence of anti-correlations in the time series, α d = 0.5 absence of correlations (white noise), and α d > 0.5 indicates the presence of positive correlations in the time series.

Conditional probability distributions
The conditional probability of an event H for a given event X is defined as where P(H \ X) is the probability that H and X jointly occur, and P(X) > 0 is the probability of the event X. The condition X reduces the statistics and increases the fluctuations of the distribution P(H|X) as compared to P(H). For this reason one associates the following error to each bin of the densities [56]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi pð1 À pÞ N r ; ð13Þ where p = P(H)dH is the probability to observe a H in the range [H, H + dH] and N is the total number of events. It is a basic result of probability theory that P(H|X) = P(X) if and only if H does not depend on X. On the contrary, P(H|X) 6 ¼ P(X) implies that H and X are not independent of each other, and their relation can be quantified by a suitable correlation measure. In the analysis of burst coupling, H and X are considered significantly correlated if P(H|X) − P(X) > � H .
To obtain a detailed picture of temporal correlations and coupling between consecutive bursts, we investigated how the individual δ-burst durations are influenced by the preceding θbursts. To this end, we evaluated the conditional probability density Pðd i d jXði À 1ÞÞ for different conditions X(i − 1) on preceding burst durations d iÀ 1 y (Fig G in S1 File) (Fig G in S1 File). Selecting only longer preceding θbursts, the probabilities for long lasting δ-bursts systematically decrease, while very short δbursts become more and more likely (Fig G in S1 File). Importantly, the analysis of conditional probabilities shows that increasing the threshold d � , i.e. progressively retaining only longer and longer preceding θ-bursts, produces a selective increase only in the probabilities of very short δ-bursts, a feature that could not be captured by the Spearman's correlation coefficient.

Model of anti-correlated burst coupling
The model consists of the following steps.
Random drawing and ranking. N durations d θ and d δ are randomly drawn from the empirical distributions previously obtained using a specific window size w. d θ and d δ are separately sorted in ascending order, i.e. from shortest to longest, and get a distinct ordinal numbers from k = 1, 2, . . ., N, which corresponds to their rank. This procedure ensures that each duration has a unique rank. The ranked d θ and d δ are then paired with a tunable degree of anti-correlation and a new time series of alternating θ-and δ-burst durations is thus generated. The coarse-grained properties of the resulting time series depends on the degree of anti-correlations used in the pairing.
Correlated pairing. Once d θ and d δ are ranked and a distinct, unique ordinal number is associated to them, one randomly choose a d θ with rank k 1 between 1 and N. To choose the following d δ , one draws a random number k 2 from a Gaussian distribution with mean μ = 1 + N − k 1 and standard deviation σ, and takes d δ as the duration corresponding to rank k 2 . This procedure is iterated N times, and at each iteration i the mean of the Gaussian from which one draws the next random rank, k i , depends on k i−1 , i.e. μ = 1 + N − k i−1 . At each iteration, k i will correspond to a duration d δ from the sorted δ-burst durations if the preceding burst was a θ-burst with duration d θ , whereas k i will select a duration d θ from the sorted θburst durations if the preceding burst was a δ-burst with duration d δ . As a result one obtains a sequence of d θ and d δ whose degree of anti-correlations is controlled by a single parameter, σ. The smaller σ, the stronger anti-correlations are.
Binary series and coarse-graining. To characterize the coarse-grained properties, the time series is first converted in a binary sequence, namely a sequence of '+' and '-'. Since each duration is by definition a multiple n of the unit window w, namely d = nw, the n windows belonging to a d θ are populated with '+', while the n windows belonging to a d δ with '-' (Fig 9c). As a results one has a sequence of windows populated with '+' and '-'. This binary sequence is then coarse-grained grouping a given number Δ of consecutive windows, with Δ odd number, and assigning '+' or '-' to the new windows of size Δ according to a majority rule, i.e. one assign '+' ('-') if the number of '+' is larger (smaller) than the number of '-' (Fig 9c). A Coarse-Grained Binary Sequence (CGBS) is thus obtained, and d CG coarse-grained durations are calculated as shown in Fig 10. Supporting information S1 File. Supporting information presenting tests to validate reported empirical results for recordings during 12-hour light and 12-hour dark periods and for different parameter values: (i) surrogate test to confirm physiological origin of the power-law distribution for θ-burst and Weibull distribution for δ-burst durations (Fig. A); (ii) range of threshold values Th for the ratio R θδ (t) between δ and θ spectral power (Fig. B and