## Abstract

Sleep encompasses approximately a third of our lifetime, yet its purpose and biological function are not well understood. Without sleep optimal brain functioning such as responsiveness to stimuli, information processing, or learning may be impaired. Such observations suggest that sleep plays a crucial role in organizing or reorganizing neuronal networks of the brain toward states where information processing is optimized.

Increasing evidence suggests that cortical neuronal networks operate near a critical state characterized by balanced activity patterns, which supports optimal information processing. However, it remains unknown whether critical dynamics is affected in the course of wake and sleep, which would also impact information processing. Here, we show that signatures of criticality are progressively disturbed during wake and restored by sleep. We demonstrate that the precise power-laws governing the cascading activity of neuronal avalanches and the distribution of phase-lock intervals in human electroencephalographic recordings are increasingly disarranged during sustained wakefulness. These changes are accompanied by a decrease in variability of synchronization. Interpreted in the context of a critical branching process, these seemingly different findings indicate a decline of balanced activity and progressive distance from criticality toward states characterized by an imbalance toward excitation where larger events prevail dynamics. Conversely, sleep restores the critical state resulting in recovered power-law characteristics in activity and variability of synchronization. These findings support the intriguing hypothesis that sleep may be important to reorganize cortical network dynamics to a critical state thereby assuring optimal computational capabilities for the following time awake.

## Introduction

Sleep is crucial for daytime functioning and well being. Although a vital part of life, its purpose and biological function are not yet well understood (Mignot, 2008). The importance of sleep is illustrated by the deteriorating effects of chronic sleep restriction or total sleep deprivation on human performance (Banks and Dinges, 2007). Without sleep optimal brain functioning such as responsiveness to stimuli, information processing, or learning may be impaired. Such observations suggest that sleep may play an important role in organizing or reorganizing neuronal networks in the brain toward states where information processing is optimized.

The general idea that both the computational capabilities of a system and its complexity are maximized at or nearby critical states related to phase transitions or bifurcations (Langton, 1990) led to the hypothesis that neuronal networks in the brain operate at or close to a critical state. The observation of neuronal activity patterns consistently following power-law distributions, a hallmark of systems at a continuous phase transition, further raised the interest in the hypothesis of critical brain dynamics (Linkenkaer-Hansen et al., 2001; Worrell et al., 2002; Beggs and Plenz, 2003; Fraiman et al., 2009; Benayoun et al., 2010; Chialvo, 2010; Poil et al., 2012). Spatiotemporal cascades of activity termed neuronal avalanches obeying a power-law size distribution were observed *in vitro* (Beggs and Plenz, 2003, 2004), *in vivo* (Gireesh and Plenz, 2008; Petermann et al., 2009; Ribeiro et al., 2010), and human magnetoencephalogram (MEG; Palva et al., 2013; Shriki et al., 2013). Recently, the spatiotemporal patterns of coherence potentials, i.e., large-amplitude negative deflections with high similarity, were reported to express neuronal avalanches in local field potentials (LFP; Thiagarajan et al., 2010; Plenz, 2012). Neuronal avalanches have been regarded an indication of balanced dynamics, i.e., avoiding regimes of overexcitation or underexcitation. The balance in activity is captured by the branching parameter σ = 1 indicating that one event on average leads to one future event resulting in the corresponding cascade size distribution to follow a power-law with exponent −3/2 (Zapperi et al., 1995; de Carvalho and Prado, 2000; Haldeman and Beggs, 2005).

The idea of a balanced regime of activity in cortical networks also extends to properties of synchrony between neuronal groups. Recent insights from *in vitro* experiments and computational modeling indicated that a balanced critical state with neuronal avalanches is also characterized by moderate mean and maximal variability of neuronal synchrony (Yang et al., 2012). As an alternative synchronization metric, the durations of transient synchronization events between channel pairs in electroencephalography (EEG) were first reported to follow a power-law density distribution in Gong et al. (2003). Later on, the distribution of phase-lock intervals (PLIs) was observed to follow a power-law function in functional magnetic resonance imaging and electrocorticographic recordings, too (Kitzbichler et al., 2009; Meisel et al., 2012). Although not limited to it (Botcharova et al., 2012), power-law scaling of PLIs also arises in computational models at criticality (Kitzbichler et al., 2009), which led to the hypothesis that its observation in neurophysiological data is indicative of a critical state of brain dynamics.

While there is growing evidence for the existence of critical states in cortical networks it is still an unresolved question how this relates to modifications of these networks in the course of wake and sleep (Pearlmutter and Houghton, 2009; Ribeiro et al., 2010; Priesemann et al., 2013). In the present work, we investigated the hypothesis that wakefulness moves cortical networks away from a critical state, which is restored by sleep. We focused on the detection of neuronal avalanches and measures of synchrony (mean, variability, the distribution of PLIs) as means to detect critical dynamics in the EEG during a period of sustained wakefulness (sleep deprivation).

## Materials and Methods

##### EEG recordings during prolonged wakefulness.

We analyzed wake EEG recordings of eight healthy young right-handed males (23.0 ± 0.46 years; mean ± SEM) during 40 h of sustained wakefulness (data from a previous study; Finelli et al., 2000). During this time, participants were under constant surveillance. The waking EEG was recorded in 14 sessions at 3 h intervals starting at 07:00. Another waking EEG was recorded after a recovery night following the sleep deprivation period amounting to a total of 15 EEG sessions. Sessions consisted of a first 5 min eyes-open period, followed by a 4–5 min eyes-closed period, and a second 5 min eyes-open period. Twenty-seven EEG derivations (extended 10–20 system; reference electrode 5% rostral to Cz) were sampled with 256 Hz (high-pass filter at 0.16 Hz; anti-aliasing low-pass filter at 70 Hz). Artifacts including eye blinks were marked by visual inspection.

##### Rating of alertness.

Alertness was self-rated on a pseudo-analog scale, which consisted of a 20 point subjective rating scale using a palmtop computer. Alertness, date, and time of rating were recorded from each subject at each EEG session (Finelli et al., 2000).

##### Detection of cascades of coherence potentials.

We used artifact-free EEG segments during the eyes-open condition. The length of segments was chosen to be 5000 samples long (corresponding to 19.53 s) to strike a compromise between, on one hand, including as many artifact-free segments as possible in the analysis and, on the other hand, having these segments long enough to provide a sufficiently large number of events. An average of 18 segments corresponding to a total of ∼6 min of EEG data was analyzed in each subject and EEG session.

Large-amplitude activity patterns in the LFP, termed coherence potentials, have previously been shown to occur without distortion at other cortex sites via fast synaptic transmission extending up to hundreds of milliseconds (Thiagarajan et al., 2010). In this work, we used coherence potentials to detect cascades of neuronal activity in space and time in EEG. Coherence potentials were defined as EEG segments across time and different EEG channels characterized by crossing a positive or negative threshold at some point and exhibiting high mutual similarity. It should be noted, that although we followed a similar approach as described by Thiagarajan et al. (2010) in identifying coherence potentials, the method applied here to the EEG is different to previous studies where coherence potentials were derived from LFP. LFP allows for detection of neuronal activity with high spatial resolution. In the EEG, the determination of the spatial location of the signal is not possible to the same extent. In this respect, the detection of cascading events in neuronal activity in our work differs from the original works on neuronal avalanches derived from LFP (Beggs and Plenz, 2003; Petermann et al., 2009; Thiagarajan et al., 2010). Furthermore, coherence potentials were originally only defined for negative potentials due to their close correlation to neuronal spiking. This justification, however, does not necessarily hold for the EEG. In view of these differences, we considered both positive and negative potentials to identify coherence potentials.

First, large positive or negative excursions beyond a threshold were identified for each EEG channel after normalization of the signal by subtraction of the mean and dividing by the SD. The comparison of the signal distribution to the best Gaussian fit indicated that the two distributions start to deviate from one another around ±4 SD (see Fig. 1*A*). We therefore chose to use a threshold of ±4 SD throughout the paper. We systematically explored other thresholds around ±4 SD to verify that results were robust for a parameter range around the chosen threshold (see Results). Once a positive or negative excursion beyond the threshold had been detected in one of the channels, we then defined the period of interest on this channel as the continuous positive or negative excursion until the signal returned the baseline (0 SD) on both sides of the threshold crossing similar to Thiagarajan et al. (2010) and termed those segments *F _{i}* trigger events (Fig. 1

*B*). The time of the signal's baseline crossing on either side of the excursion therefore determined the length of the trigger event.

Second, coherent segments between a trigger event *F _{i}* and periods

*F*of equal duration in all channels were identified when correlations

_{j}*R*as follows: were

*R*≥ 0.75. A correlation threshold in this range was identified by Thiagarajan et al. (2010) to robustly detect coherence potentials. We systematically changed the correlation threshold to verify that results were robust for a broad range around this value (see Results). Together with the trigger event, these coherent segments form the set of coherence potentials.

Next, coherence potentials were discretized in a raster plot to identify cascades in space and time. The time point of a single event was defined as the time point of the occurrence of the minimum (or maximum, depending on the trigger event) for each coherence potential (Fig. 1*C*). For each trigger event, the time series of events from each channel was discretized with time bins of duration Δ*t*. Along the line of previous work (Beggs and Plenz, 2003; Gireesh and Plenz, 2008; Ribeiro et al., 2010; Friedman et al., 2012; Priesemann et al., 2013; Shriki et al., 2013), a continuous sequence of time bins with events in any channel ending with a time bin with no events in any channel was defined as a cascade of events. The cascade size *S* was defined as the number of events in all channels in a cascade. To prevent double counting of events and cascades we required trigger events not to have been an event in a previous cascade.

For control, we generated phase-shuffled surrogate data for each channel, which shuffles the phase of the original signal but keeps the amplitude spectra unchanged in the frequency domain (Theiler et al., 1992).

As an alternative approach to coherence potentials, we also derived neuronal avalanches by identifying large positive or negative deflections above a threshold at single EEG channels and combined them into spatiotemporal cascades by discretizing them in time by their most extreme excursion as described previously (Palva et al., 2013; Shriki et al., 2013). This method was recently shown to allow identification of neuronal avalanches in EEG and MEG.

##### Branching parameter σ.

The branching parameter σ quantifies the ratio of future events to ancestor events. Previously, it has been used to characterize cascading events on different scales, from activity in neuronal cultures to events in the MEG. We calculated σ as the ratio of the number of events in the second time bin of a cascade to the one in the first time bin. The ratio was averaged over all cascades for each subject and each EEG session.

##### Synchrony and entropy measurements.

We followed the approach described previously (Yang et al., 2012) to quantify mean and variability of synchronization. The calculation was performed on the same artifact-free segments used for the derivation of coherence potentials. After filtering the data in the alpha band (8–12 Hz), we first obtained a phase trace θ* _{i}*(

*t*) from each EEG trace

*F*(

_{i}*t*) using its Hilbert transform

*H*[

*F*(

_{i}*t*)]: Next, we quantified the mean synchrony in each EEG segment by where

*L*= 5000 is the length of our EEG segments in samples and

*r*(

*t*) is the Kuramoto order parameter: which was used as a time-dependent measure of phase synchrony with

*n*= 27 being the number of EEG channels in our data.

As a measure for the variability of synchronization we derived the entropy of *r*(*t*) in each EEG segment by
where we estimated a probability distribution of *r*(*t*) by binning values into intervals. *p _{i}* is then the probability that

*r*(

*t*) falls into a range

*b*<

_{i}*r*(

*t*) ≤

*b*

_{i}_{+1}. Similar to Yang et al. (2012), we found results to be robust over a broad range for the number of bins

*B*used. We applied

*B*= 24 in the current analysis.

##### Derivation of the distribution of PLIs.

PLIs were calculated for all possible pairs of derivations of artifact-free EEG segments of 19.53 s duration (5000 samples, same segments as for the analysis of coherence potentials and synchronization measures). The analysis was performed for scale 4 (8–16 Hz; alpha band; see below) and scale 5 (4–8 Hz; theta band).

##### Hilbert wavelet transforms.

To derive a scale-dependent estimate of the phase difference between two EEG channels, we follow the approach described previously (Kitzbichler et al., 2009) using Hilbert transform derived pairs of wavelet coefficients (Whitcher et al., 2005). We define the instantaneous complex phase vector for two signals *F _{i}*(

*t*) and

*F*(

_{j}*t*) as follows: where

*W*denotes the

_{k}*k*th scale of a Hilbert wavelet transform and † its complex conjugate. Here

*F*(

_{i}*t*) and

*F*(

_{j}*t*) are different EEG derivations. A local mean phase difference in the frequency interval defined by the

*k*th wavelet scale is then given by with being a less noisy estimate of

*C*(

_{i,j}*t*) where 〈 · 〉 indicates the temporal average over a time window Δ

*t*= 2

*8 in sampling steps (Kitzbichler et al., 2009).*

^{k}Intervals of phase locking can then be identified as periods when |Δφ* _{i,j}*(

*t*)| is smaller than some arbitrary threshold, which we set to π/4. We also require the modulus squared of the complex time average,

##### Estimation of the goodness of power-law fit.

To determine the quality of the fit of a power-law function to the observed distribution we performed a goodness-of-fit test, which leads to a *p* value that quantifies the plausibility of the hypothesis that the distribution is power-law like (Clauset et al., 2009). Upon fitting the empirical distribution, power-law distributed synthetic data were generated with parameters derived from the power-law fit and individually fit to their own power-law model. The *p* values were then defined as the fraction of times the resulting Kolmogorov–Smirnov statistics for each synthetic data (*n* = 1000) relative to its own model was larger than the value of the empirical data. Larger *p* values therefore indicate a higher probability that the observed distribution could be explained by an underlying power-law characteristic. Conversely, low *p* values make a power-law distribution a less likely distribution to describe the empirical data.

As a measure quantifying the difference of a given empirical distribution of some quantity *X* (in our case *S* or PLI) from a power-law distribution with exponent α, we defined
for which the cumulative probability distribution *P*(*X*) of the distribution's tail, i.e., the *N _{X}* number of

*X*values with values larger than some minimal value

*X*

_{min}, was used. For distributions of cascade sizes the minimum

*S*

_{min}was set to 1 and cascade sizes up to system size (i.e., number of channels 27) were included in the sum. For the distribution of PLI the minimal value

*PLI*

_{min}given by the fitting algorithm was applied.

##### Power density spectra.

EEG power density spectra were computed for derivation C3A2 of artifact-free 20 s epochs (same starting points as segments used for the derivation of coherence potentials, synchronization measures, and PLIs; fast Fourier transform routine; Hanning window; frequency resolution 0.25 Hz). Power in the theta (5–8 Hz) and alpha (8.25–12 Hz) range was determined and averaged for eyes-open segments per session and subject.

## Results

### Coherence potentials organize as neuronal avalanches in the EEG

We investigated the spatiotemporal organization of coherence potentials in artifact-free EEG intervals in the eyes-open condition. Artifact-free EEG intervals were analyzed in segments of length 5000 samples (corresponding to 19.53 s; see Materials and Methods). Multiple segments were analyzed for each subject and EEG session. On each of these segments, we first determined potentials with either a positive or a negative deflection larger than a certain threshold, which we termed trigger events and second identified segments with high similarity to these trigger events (see Materials and Methods). Amplitude distributions from EEG signals start to deviate from a Gaussian distribution for deflections larger than 4 SD (Fig. 1*A*). For the detection of coherence potentials we therefore focused on a threshold of ±4 SD throughout the paper. Systematic exploration of other thresholds around ±4 SD verified that results were independent of the exact choice of threshold (Fig. 2*C*,*D*). Figure 1*B* shows some exemplary trigger events, their mean duration in all subjects was 468 ± 5 ms. Next, segments of the same length and with high similarity to trigger events were located. Potentials were considered similar to a trigger event when their correlation *R* was equal to or greater than 0.75 (see Materials and Methods). Again, systematic analysis with different *R* values confirmed that results were robust over a range of values (Fig. 2*C*,*D*).

We observed these coherence potentials to organize in cascades of continuous events in time and space (*n* = 59258 cascades; Fig. 1*C*). Cascade sizes exhibited a high degree of variability. At the beginning of the sleep deprivation period, the probability distribution of cascade sizes *S* closely followed a power-law function *p*(*S*) ∼ *S*^{α} with α close to −3/2 (Fig. 1*D*). The total number of EEG channels determined the observed cutoff of the power-law function at ∼27. In line with reports of neuronal avalanches observed in other model and experimental systems (Beggs and Plenz, 2003; Priesemann et al., 2013; Shriki et al., 2013), the distribution of cascade sizes remained a power-law distribution for different bin sizes Δ*t* with shallower slopes for larger Δ*t*. Conversely, phase shuffling of the data destroyed the power-law distribution indicating that the long-range spatiotemporal correlations were captured by the power-law distribution of cascade sizes (Fig. 1*D*). Along with a power-law exponent of ∼−3/2, we observed the branching parameter, i.e., the ratio of descendant events to ancestor events, to be σ = 1.17 ± 0.03 (Fig. 2*B*).

### Sustained wakefulness leads to changes in the distribution of cascade sizes and branching parameter σ

With increasing duration of sustained wakefulness, we observed changes in the distributions of cascade sizes. During prolonged wakefulness the probability for larger cascade sizes increased giving rise to a shallower power-law slope with a heavier tail in the distribution (Fig. 2*A*). We quantified the deviation of the distribution's tail from a power-law function with exponent α = −3/2 and calculated Δ*D* by combining cascade sizes into one distribution for each subject and EEG session (see Materials and Methods). To account for differences in absolute values between subjects, Δ*D* and σ values were first transformed into *z*-scores, i.e., subtraction of the mean and division by the SD, for each individual before averaging over subjects. Regardless of the underlying distribution, *z*-scores reflect performance relative to some group, rather than relative to an absolute standard. We denote the normalized data by subscript *n* throughout the paper. With growing time awake Δ*D* progressively increased. Similarly, the branching parameter σ increased during sustained wakefulness. Figure 2*B* illustrates the evolution of the two metrics as a function of time awake.

To statistically evaluate the changes of Δ*D* and σ in the course of sleep deprivation we averaged the values of the first four recordings (0–9 h awake) and compared them to the corresponding averages over the last four recordings (30–39 h awake). The increase in both metrics was significant for a wide range of correlations *R* and thresholds (Fig. 2*C*; bar heights reflect the difference of average metrics at 30–39 h and at 0–9 h, two-tailed paired *t* test). The return of σ and Δ*D* to lower values after recovery sleep was similarly observed across a broad range of parameters (significant decrease of values after recovery sleep compared with the value after 39 h of wakefulness; two-tailed paired *t* test; Fig. 2*D*). In absolute values, σ increased from 1.17 ± 0.03 during the first four EEG recordings (0–9 h awake) to 1.45 ± 0.09 at the end of the sleep deprivation period (30–39 h awake; Fig. 2*B*, middle).

We also derived cascades of activity using only events with a positive or negative excursion larger than a certain threshold as previously reported (Palva et al., 2013; Shriki et al., 2013). Using only the threshold as an event detection criterion substantially lowered the number of detected events and cascade sizes when compared with the approach involving coherence potentials (5992 cascades vs 59,258 cascades using coherence potentials). This alternative approach similarly produced power-laws of cascade sizes and a significant increase of σ with increasing duration of sleep deprivation (Fig. 3).

### Decreasing variability and increase of mean synchrony during sustained wakefulness

The extent and variability of synchrony has been demonstrated to be sensitive to the balance between excitation and inhibition similar to cascades of activity (Yang et al., 2012). Given the changes observed in the distribution of cascade sizes and σ in our data during sustained wakefulness, we next investigated whether these were accompanied by corresponding alterations in synchronization metrics. We calculated mean and variability of synchrony from the same artifact-free segments used for the analysis of coherence potentials during the eyes-open condition. In our analysis we focused on the alpha and theta frequency bands since power changes in these bands during wake were found to be associated with sleep propensity (Torsvall and Akerstedt, 1987; Cajochen et al., 1995; Aeschbach et al., 1999; Finelli et al., 2000; Strijkstra et al., 2003). During sustained wakefulness we observed an increase in mean synchrony as a function of time awake in the alpha band while the variability of synchrony quantified by its entropy decreased significantly (Fig. 4). After consecutive recovery sleep (*p*s) both metrics recovered in the direction of initial values. Interestingly, no significant changes were observed in the theta band for the eyes-open condition. Extending our analysis to the eyes-closed condition in the theta band revealed similar, albeit weaker changes than in the alpha band in both synchronization metrics indicating the observed effects to occur predominantly in the alpha band.

The initially low mean and high variability of synchrony in our data closely resemble the findings in Yang et al. (2012) where a maximum of phase synchrony was found to be associated with a balanced regime characterized by neuronal avalanches and the onset of synchrony.

### Distribution of PLIs

To further test for changes in synchronization in the EEG we computed the distribution of PLIs. We calculated PLIs between all pairs of EEG derivations of the same artifact-free EEG segments during the eyes-open condition that were used for the derivation of coherence potentials and the analysis of the mean and variability of synchrony. The probability distribution of PLIs closely followed a power-law function *p*(*PLI*) ∼ *PLI*^{−α} as reported previously (Gong et al., 2003; Kitzbichler et al., 2009) with estimated exponents α in the range of 2 to 2.5 (Fig. 5*A*). The distribution's *p* value based on the Kolmogorov–Smirnov statistic provides a measure of the plausibility of a power-law fit to the data (Clauset et al., 2009). Statistical analysis revealed high *p* values in many time intervals suggesting that the power-law hypothesis cannot be rejected. A recent comprehensive analysis of various fitting functions applied to PLI distributions had revealed a power-law function to be the most likely fit (Kitzbichler et al., 2009).

During prolonged wakefulness, we observed changes in the PLI distributions similar to the ones for cascades of coherence potentials. With increasing duration of sleep deprivation, the probability for longer PLIs increased (Fig. 5*C*) giving rise to an increasing deviation from a power-law distribution. Δ*D* (see Materials and Methods) again quantifies the deviation of the distribution's tail from a power-law function with exponent α. Both *p* and Δ*D* values can be regarded as complementary measures characterizing a power-law distribution with on average larger *p* values corresponding to smaller Δ*D* values (Fig. 5*B*). Our analysis was performed on cumulative distributions *P*(*PLI*), with α and *PLI*_{min} being averaged values of the power-law fit of each participant's first EEG recording (0 h of wakefulness). To compare *p* and Δ*D* values and to account for differences in absolute values between subjects, values were first transformed to *z*-scores for each individual before averaging over subjects. Figure 6 shows normalized *p* values (left), Δ*D* (middle), and spectral power (right) averaged over all eight participants in the 8–16 Hz frequency band (alpha band; corresponding to scale 4; see Materials and Methods) and the 4–8 Hz frequency band (theta band; corresponding to scale 5). Along with time course at 3 h intervals, averages of values of the first (blue, 0–9 h of wakefulness) and last four recording sessions (red, 30–39 h of wakefulness) are illustrated. With increasing duration of sleep deprivation the PLI distributions deviated stronger from a power-law distribution (as quantified by the decrease in *p* values) showing a larger tail (captured by the increase in Δ*D* values). The change in both measures was significant for the alpha band (two-tailed paired *t* test; Fig. 6*A*) and qualitatively also observed in the theta band (Fig. 6*B*). Similar to the cascade size of coherence potentials, the observed changes correspond to the increasing incidence of larger events compared with scale-free activity observed earlier in the course of sustained wakefulness.

Previously, changes in EEG spectral power have been described in the alpha and theta bands during sustained wakefulness (Torsvall and Akerstedt, 1987; Cajochen et al., 1995; Aeschbach et al., 1999; Finelli et al., 2000; Strijkstra et al., 2003). The time-dependent changes in the metrics of PLI distributions as well as mean and variability of synchronization in these two frequency bands (scales) differed from prior observations in the time course of spectral power in the corresponding frequency bands during sleep deprivation. Normalized power calculated from derivation C3A2 for the same EEG segments showed a predominantly circadian pattern for the alpha band (Fig. 6*A*, right plot; note the period of ∼24 h in spectral power) and a circadian component superimposed on an increasing trend for the theta band (Fig. 6*B*, right plot) similar to previously reported data (Aeschbach et al., 1999; Finelli et al., 2000) while the indices *p* and Δ*D* quantifying the PLI and also coherence potential distribution as well as the other metrics σ, <*r*(*t*)> and *H*(*r*(*t*)) exhibited a more monotonic trend.

### Correlation of EEG indices to subjective alertness

Self-rated alertness declined significantly with increasing time awake (Fig. 7*A*). Before averaging over subjects, alertness scales were transformed to *z*-scores for each subject in the same way as the EEG indices. To compare its evolution with EEG indices, we calculated their correlation coefficients *R* for which we used averaged data over all subjects (Fig. 7*B*; see Materials and Methods). Similar to a previous report (Finelli et al., 2000), theta power and alertness were negatively correlated (*r* = −0.94) and exhibited the largest absolute correlation value of all indices. The avalanche metrics and synchronization measures also exhibited a high correlation with alertness. From these indices, σ correlated best with alertness (*r* = −0.77).

## Discussion

In the present work, we reported several changes in EEG indices during sustained wakefulness: the organization of cascade sizes of coherence potentials, the mean and variability of synchronization, and the distribution of PLIs. These indices changed as a function of time awake and recovered toward baseline values after subsequent recovery sleep. The changes in synchronization measures (mean synchronization, variability of synchronization, and distribution of PLIs) were most predominant in the alpha band. This is in contrast to the changes in spectral EEG power as a marker of sleep propensity, which are most evident in the theta band (Finelli et al., 2000). Thus, the changes in synchronization cannot be explained as a direct consequence of the alterations in spectral power.

At the beginning of the sleep deprivation period coherence potentials in the EEG were organized as neuronal avalanches, i.e., events in space and time characterized by a power-law size distribution with exponent close to −3/2 and branching parameter of ∼1 (σ = 1.17). Neuronal activity patterns following a power-law size distribution with exponent −3/2, a lifetime distribution with power-law exponent −2, and a branching ratio of 1 have previously been interpreted as an indication that cortical neuronal networks operate at criticality (Beggs and Plenz, 2003; Larremore et al., 2011; Friedman et al., 2012; Palva et al., 2013; Shriki et al., 2013). This conclusion is based on insights derived from theoretical models of systems poised at a phase transition exhibiting the same scaling laws. The case of neuronal avalanches refers to a critical branching process giving rise to scale-free avalanches of activity and avoiding activity regimes comprising only either large or small parts of a network (Bak and Paczuski, 1995). This balanced propagation of activity is also reflected in the branching parameter, which in such a case is expected to be exactly 1. Although close to 1, we observed a baseline branching parameter of 1.17, which could be caused by uncertainty in its exact determination or, if taken by itself, indicate a slightly supercritical state.

With increasing duration of wakefulness, both size distributions of coherence potentials and PLIs increased in tail as measured by Δ*D* and *p* values in the case of PLIs. For the neuronal avalanches we observed a coincidental increase in σ, defined as the ratio between descendant events to ancestor events. Furthermore, mean synchronization increased while its variability decreased.

A critical branching process provides an interesting framework connecting these seemingly different observations. It was recently shown that power-law neuronal avalanches are accompanied by a maximum of synchronization variability *in vitro* and in modeling systems (Yang et al., 2012). Conversely, disinhibited networks exhibited cascades of activity with an imbalance toward larger avalanches, which coincided with decreased variability in synchronization and increased mean of synchronization in some metrics. These findings closely resemble the observations in our data. In the context of a critical branching process they can be interpreted in the sense that network activity is increasingly shifted toward a state in which larger events prevail in the dynamics, unlike in the critical state with scale-free avalanches of activity and high variability of synchronization.

The balance between excitation and inhibition depends crucially on synaptic strength in neuronal systems. It was recently demonstrated that a reduction of inhibitory synaptic transmission by application of GABA receptor blockers resulted in an artificial imbalance toward excitation related to a supercritical regime with larger neuronal avalanches than predicted by a power-law characteristic (Beggs and Plenz, 2003; Shew et al., 2009). Although the changes in such pharmacologically disinhibited networks were more drastic than what we observed, they qualitatively well correspond to alterations in the distribution of neuronal avalanches and σ reported here. Indirect evidence for changes in excitability and synaptic strength during wake and sleep comes from observations in the EEG, where power changes in the theta and alpha bands during wake were found to be associated with sleep propensity (Torsvall and Akerstedt, 1987; Cajochen et al., 1995; Aeschbach et al., 1999; Finelli et al., 2000; Strijkstra et al., 2003) and slow-wave activity during sleep to be associated with sleep intensity reflecting a regulatory process termed sleep homeostasis (Achermann and Borbély, 2011). Tononi and Cirelli (2003, 2006) hypothesized that synaptic homeostasis is underlying sleep homeostasis: synaptic strength is high at the beginning of a night, due to plastic processes occurring during waking, and decreases by means of synaptic downscaling during sleep. Recent investigations of changes in synaptic strength in *Drosophila* (Bushey et al., 2011) and neuronal excitability in human cortex (Huber et al., 2012) support the hypothesis by providing evidence of structural changes occurring in neuronal networks during waking and their reorganization during sleep. It is conceivable that the increase in excitability of cortical circuits caused by the buildup in synaptic strength provides the cellular basis for the observed shift away from critical dynamics in the course of wakefulness. Similarly, rebalancing and synaptic downscaling during sleep could tune network dynamics back to criticality. However, a plausible hypothesis for the detailed mechanisms, how the changes in synaptic strength could induce the observed changes both in spectral power and in the signatures of criticality, is still missing. From a clinical perspective, disinhibited networks are reminiscent of epileptic seizures where activity engages most of the network. In fact, for many forms of epilepsy sleep deprivation is known to increase the probability of seizure occurrence (Ellingson et al., 1984).

Power-law probability distributions of PLIs between pairs of neurophysiological time series were recently studied (Kitzbichler et al., 2009; Meisel et al., 2012) and interpreted as signatures of criticality in human brain dynamics. Evidence for power-law PLI distributions as an indicator for criticality derives from computational models, which show power-law PLI distributions at the transition between an ordered and a disordered phase, i.e., when they are in a critical state (Kitzbichler et al., 2009; Meisel et al., 2012). A PLI power-law distribution can, however, also arise when systems are away from a phase transition making the PLI distribution a sensitive but not specific indicator for critical dynamics (Botcharova et al., 2012). A clear power-law PLI distribution in neurophysiological data can, as any power-law scaling (Touboul and Destexhe, 2010; Beggs and Timme, 2012), therefore only be seen as an indication for the possibility of a phase transition. In our case, this indication is further supported by the simultaneous detection of neuronal avalanches in coherence potentials following a power-law distribution with exponent −3/2, a σ of 1.17, and high variability in synchronization. Conversely, PLI probability distributions with a large tail and weak power-law statistics along with similar changes in the distribution of neuronal avalanches and a σ much larger than 1 together with reduced variability of synchronization support the conclusion that the system is not directly located at a critical state.

It was previously hypothesized that normal activity is located in a slightly subcritical regime and that sleep could function to establish a safety margin by reorganizing activity toward it while wakefulness drives it closer to criticality (Pearlmutter and Houghton, 2009). A recent analysis of neuronal avalanches from human intracranial depth recordings in fact implied that the human brain might not operate at criticality directly but in a somewhat subcritical regime (Priesemann et al., 2013). While we cannot exclude the possibility that network activity in the brain is poised in a slightly subcritical regime, when interpreted in the context of criticality, the fading power-law statistics of the PLI distribution, the increase of σ to values much larger than 1, and the decrease in variability of synchronization in our opinion provide more support to the notion that cortical networks are near or at a critical state initially and shift toward a supercritical regime during wakefulness. Furthermore, the characteristic changes of synchronization metrics accompanying the alterations in the power-law characteristics of neuronal avalanches suggest a growing deviation from criticality as a function of time awake instead of the system remaining critical with a different exponent as could be concluded by looking solely at the distribution of neuronal avalanches.

Critical dynamics are often regarded to support optimal computational functioning (Langton, 1990; Bertschinger and Natschläger, 2004; Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006; Shew et al., 2009, 2011). The observation of fading signatures of critical dynamics during prolonged wakefulness and their correlation to alertness could provide an interesting link to behavioral observations of impaired cognitive functioning and information processing after sleep deprivation (Banks and Dinges, 2007). Our findings support the intriguing hypothesis that sleep might serve to reorganize cortical network dynamics to a critical state thereby assuring optimal computational capabilities for the time awake.

## Footnotes

The study was supported by the Swiss National Science Foundation Grant 320030-130766 (P.A.) and by the European Community's Seventh Framework Programme (FP7/2007-2013) under Grant agreement no. 258749 (E.O.). We thank Thilo Gross for comments on an earlier version of this manuscript.

- Correspondence should be addressed to Christian Meisel, Department of Neurology, University Clinic Carl Gustav Carus, Fetscherstrasse 74, 01307 Dresden, Germany. christian{at}meisel.de