Discrete Scale Invariance of Human Large EEG Voltage Deflections is More Prominent in Waking than Sleep Stage 2

Electroencephalography (EEG) is typically viewed through the lens of spectral analysis. Recently, multiple lines of evidence have demonstrated that the underlying neuronal dynamics are characterized by scale-free avalanches. These results suggest that techniques from statistical physics may be used to analyze EEG signals. We utilized a publicly available database of fourteen subjects with waking and sleep stage 2 EEG tracings per subject, and observe that power-law dynamics of critical-state neuronal avalanches are not sufficient to fully describe essential features of EEG signals. We hypothesized that this could reflect the phenomenon of discrete scale invariance (DSI) in EEG large voltage deflections (LVDs) as being more prominent in waking consciousness. We isolated LVDs, and analyzed logarithmically transformed LVD size probability density functions (PDF) to assess for DSI. We find evidence of increased DSI in waking, as opposed to sleep stage 2 consciousness. We also show that the signatures of DSI are specific for EEG LVDs, and not a general feature of fractal simulations with similar statistical properties to EEG. Removing only LVDs from waking EEG produces a reduction in power in the alpha and beta frequency bands. These findings may represent a new insight into the understanding of the cortical dynamics underlying consciousness.


INTRODUCTION
Although rarely explicitly stated, the dominant model for the analysis of human EEG signals for more than 50 years has been spectral analysis, implicitly viewing the time-dependent changes in cortical local field potentials as a set of dynamic and standing electrical waves originating from the top 5 mm of the cerebral cortex (Nunez and Srinivasan, 2006;Nunez, 2010). Such techniques have proven spectacularly successful in both research and clinical medicine (Schwilden, 2006;Arciniegas, 2011;Schiff et al., 2014).
However, many lines of evidence originating in multielectrode array recordings of rodent cortex (Plenz and Thiagarajan, 2007;Gireesh and Plenz, 2008;Klaus et al., 2011), extending through human electrocorticography (Priesemann et al., 2013), electroencephalography (EEG; Poil et al., 2012;Palva et al., 2013), magnetoencephalography (MEG; Shriki et al., 2013;Yu et al., 2013), and magnetic resonance imaging (MRI; Kitzbichler et al., 2009) have resounding demonstrated that a fundamental organizing principle of neuronal activity in the cerebral cortex is via dynamic distributions of scale-free ''neuronal avalanches''. Power-lawdistributed neuronal avalanches have been shown to demonstrate many hallmarks of dynamic systems in a critical state Yu et al., 2013), which simulations have shown are likely to lead maximal information transmission throughout the neuronal ensemble (reviewed in Shew and Plenz, 2013). A priori, such non-stationary dynamical systems would not be best analyzed via spectral analysis (Weisstein, 2014).
EEG has also been linked to neuronal avalanche activity. The developmental time course of neuronal avalanche maturation was studied in a rodent model, where the onset of powerlaw distributed neuronal avalanches was correlated with the maturation of beta and gamma band EEG signals (Gireesh and Plenz, 2008). Neural network models with avalanche dynamics can quite readily produce a signal that approximates human EEG, however with reduced power in the 20-50 Hz (''beta'' and ''gamma'' bands;de Arcangelis and Herrmann, 2012). Many convincing demonstrations of neuronal avalanche activity in rodents, non-human primates, and humans have been reported, showing that functional mammalian cortical neuronal ensembles organize into scale-free, critical-state dynamics, with a characteristic power law exponent and branch ratio (Shew and Plenz, 2013;Shriki et al., 2013). Additionally, EEG and MEG signals themselves have been shown to be highly non-stationary, albeit amenable to analysis as ''quasi-stationary'' (Kaplan et al., 2005).
While compelling, these observations beg the question as to what role, if any, neuronal avalanches may play in the origin of conscious mental states, especially given the clear observations that rodent brain tissue in vitro exhibits similar power law behavior and critical state dynamics via branch ratio analysis of avalanches to that of waking human MEG signals (Plenz and Thiagarajan, 2007;Shriki et al., 2013). In fact, cortical neuronal avalanche size has been found to follow a power law with an exponent of −1.5 in rodent cortex cultures in vitro, in vivo anesthetized rodent cortex, and awake primates (Klaus et al., 2011). Consciousness is thought to be mediated at the neuronal level via cortical-thalamic feedback loops capable of global informational binding without strict localization (Edelman and Tonioni, 2000;Baars et al., 2013;Boly et al., 2013), and many investigators have pointed out the theoretical advantages of critical-state neuronal dynamics in terms of rapid, global cortical information transfer (Poil et al., 2012;Palva et al., 2013;Shew and Plenz, 2013). Cortical avalanche dynamics have been explored using power law techniques in awake vs. sleeping rodents, and found to have a very similar structure (Ribeiro et al., 2010). By contrast, avalanche distributions were found to differ somewhat among sleep stages in a group of human epileptic patients implanted with hippocampal electrodes (Priesemann et al., 2013).
Investigators have also studied changes in neuronal physiology with conscious vs. unconscious mental states using methods looking for Long Range Temporal Correlation (LRTC) in neuronal physiology in many different systems (Li et al., 2008;Nikulin et al., 2012;Blythe et al., 2014). Using multiscale entropy and detrended fluctuation analysis (DFA) of renal sympathetic nerve activity, anesthetized rats were found to have decreased LRTC compared to awake rats (Li et al., 2008). Using Hilbert transformed alpha and gamma band information, subsequent application of DFA demonstrated that subjects with schizophrenia had decreased alpha and beta band LRTC compared to healthy controls (Nikulin et al., 2012). A further application of LRTC in EEG demonstrated that an improved source localization method was able to improve LRTC estimates using DFA (Blythe et al., 2014).
With an eye to these data, we make the fairly obvious observation that critical-state neuronal avalanche dynamics are unable to explain much of the regularity in structure seen in waking human EEG signals, both when analyzed by spectral analysis and autocorrelation function (ACF) analysis. Spectral analysis and ACF analysis yield markedly diferent results for waking and sleep 2, whereas cortical avalanche dynamics have been shown to be very similar across different states of cortical function, both in vivo and in vitro (Plenz and Thiagarajan, 2007;Klaus et al., 2011;Shriki et al., 2013). We utilized a publicly available EEG dataset, taking advantage of continuous somnogram information to get waking and sleep stage EEG tracings from the same individuals, in order to directly compare changes in EEG large voltage deflection (LVD) dynamics within subjects. Taking the viewpoint that ensembles of neuronal avalanches contribute to EEG, we hypothesized that discrete scale invariance (DSI) might be superimposed on criticalstate dynamics in waking more than in sleep stage 2 LVDs. This is because log-periodic oscillations in the observables of natural processes can be caused by a partial breaking of scale invariance, resulting in a ''periodicity'' in the data (cf. Sornette, 1998Sornette, , 2006Zhou et al., 2005). This hypothesis could represent a new insight into our understanding of the cortical neuronal dynamics of brain function, including consciousness.
For the following brief description of DSI, we adapt from several expositions by Sornette (1998Sornette ( , 2006. This report is too brief to provide a comprehensive review of DSI, for which we refer the interested reader to Sornette (1998) and references therein. Power-law-distributed systems are typically characterized by continuous scale invariance, where an observable function O(x) is scale invariant under the arbitrary change x → λx (where λ → 1 + ) when there is a function µ(λ) such that: The solution to which is the basic power law formula: where C is a constant and Studies in complex systems prone to catastrophes, including financial markets, ruptures of pressurized containers, and earthquakes have demonstrated that certain dynamical systems can exhibit a weaker form of DSI where the control parameter λ acts as a scaling ratio, and can only exhibit discrete values, e.g., λ 1 , λ 2 , λ 3 . . . , λ n (Sornette, 1998(Sornette, , 2006. As a result, in DSI, the power law behavior of the underlying system is also controlled by complex exponents decorating the power laws, leading to periodicity (Sornette, 1998(Sornette, , 2006: where n is an integer and i the imaginary unit. In the case of n = 0, the imaginary part drops out and continuous scale invariance is recovered as in (3).

EEG Data
We utilized a publicly available database of single-channel EEG recordings from the MIT-BIH polysomnographic database recorded using the 10-20 EEG reference system (slpdb; http://www.physionet.org; Ichimaru and Moody, 1999;Goldberger et al., 2000). We utilized the same dataset previously reported, and details of the subjects and EEG data utilized are as described (Zorick and Mandelkern, 2013). Briefly, our database consists of eight different non-contiguous 1 min EEG tracings (recorded at 250 Hz) from 14 subjects in both waking and sleep stage 2 consciousness, one lead per subject. Approval for this study was provided by the local VA West Los Angeles Institutional Review Board.

Simulations
Our choice of simulations was motivated by the fact that these simulations exhibit similar fractal (and in some cases, multifractal) properties to that of human waking EEG (Zorick and Mandelkern, 2013). The log normal sigma 0.1 multifractal series (32,768 data points) was downloaded from http://www.physionet.org/physiotools/multifractal/, made from the log-normal wavelet cascade algorithm with parameters v = ln(2)/4 and σ = 0.1 as described (Muzy et al., 1993). The fractional Brownian motion monofractal series was generated with Hurst exponent (H) of 0.2 using the ''dvfBm'' R package (120,000 data points; version 1.0; Coeurjolly, 2001). The binomial multifractal series (BMS) was constructed with 120,000 data points after Kantelhardt et al. (2002); Zorick and Mandelkern (2013). The BMS is a series of N = 2 n max numbers with index k = 1, . . ., N, defined by x k = a n(k−1) (1 − a) n max −n(k−1) For this series, a is a user-defined parameter which can take values 0.5 < a < 1. We chose the parameter a = 0.6 such that the resulting multifractal spectrum roughly matches that of the human waking EEG samples (Zorick and Mandelkern, 2013). Here n(k) is the sum of digits equal to 1 in the binary representation of the index k (120,000 data points). As an example, choosing an index value of k = 13, n (13) = 3, as the binary representation of the decimal number 13 is 1101.

EEG Large Voltage Deflections
Based upon data from MEG showing that only spikes >3 SD from the mean in size are able to be differentiated from the Gaussian distribution , we identified LVDs of size >2.5 SD from the mean of the EEG segment, separately both for positive and negative voltages. The results obtained are quite similar when identifying LVDs of size >2 SD from the mean. However, the relatively short length of available EEG traces and the relatively coarse sampling rate (250 Hz) limited analysis in our dataset to thresholds <3 SD from the mean, as thresholds higher than this would not generate enough LVDs to reliably analyze in subsequent steps. We defined LVD size s as the sum of voltage readings between times at which the EEG voltage is zero (zero crossings). We utilized LVD size as an appropriate measure for further analysis, as it led to robust PDF estimations given adequate sequence length. Each EEG segment was analyzed separately for LVD PDF estimation and periodicity.

Probability Density Function Analysis
We utilized the R ''density'' function (RCoreTeam, 2012) to obtain the normalized probability density function (PDF) f(s) for each LVD segment. Segments with fewer than 50 avalanches per segment did not produce reliable PDF estimates, so these were excluded. After PDF estimation, the scale sizes f(s) were natural log transformed for further analysis.

Lomb-Scargle Periodogram Analysis
To assess for DSI in EEG LVD sizes, we utilized the Lomb-Scargle nonparametric periodogram method from the log-scale PDF f(s) of the neuronal avalanches. The Lomb-Scargle method is motivated as the necessary logarithmic transformation step of the frequencies presents non-uniformly spaced samples (Zhou and Sornette, 2002;Zhou et al., 2005;Press et al., 2007). The R function ''spec.ls'' in the ''cts'' package (written by Zhu Wang) was used, which was designed to follow (Press et al., 2007). Normalized Lomb-Scargle periodograms were constructed by dividing by the variance of the frequency data f(s) to permit statistical testing (Press et al., 2007). We utilized normalized Lomb-Scargle periodograms to assess for statistical significance using the probability threshold of p = 0.01, with P(>z) Me −z , where z is the threshold value at the desired significance level, and M is a maximum estimate of the number of independent frequencies in the periodogram (which we heuristically estimate at twice the number of data points in the input series (512 data points; Press et al., 2007).

Statistical Analysis
All statistical analyses were performed using R (RCoreTeam, 2012) and SPSS (IBM, Chicago, 2014). Statistically significant Lomb-Scargle periodogram peaks were broadly defined as any point with a positive slope immediately preceding, and a negative slope immediately following. For each segment we have the number of significant (p < 0.01) Lomb-Scargle peaks and the area of the Lomb-Scargle periodogram above the p = 0.01 probability threshold, for both positive-and negative-LVD periodograms. We analyze these four data sets separately using a linear mixed model, regressing the number of peaks (or area) against consciousness condition (waking or sleep stage 2). The covariance model for the within-subject measurements is taken as compound symmetry.

Waking EEG Differs from Sleep Stage 2 EEG by Both Spectral Analysis and Autocorrelation Function Analysis
Using a dataset with 8 min of EEG from both waking and sleep stage 2 in 14 individuals, we performed spectral analysis via FFT and ACF analysis on each of eight separate non-contiguous 1 min tracings for each subject (Figure 1). While by no means a novel finding, we simply make the observation that waking EEG is readily distinguished from sleep stage 2 EEG via spectral analysis ( Figure 1A). We furthermore note the characteristic oscillatory patterns in the ACF that differ markedly between waking and sleep stage 2 EEGs ( Figure 1B).

LVD PDFs Appear to have a Similar
Structure for Waking and Sleep Stage 2, though more Avalanches are found in Waking We isolated 2.5 SD LVDs from the 14 subjects' waking and sleep stage 2 EEGs, and show a semilog plot of the PDF vs. LVD size for one subject (Figure 2). Positive (Figure 2A) and negative ( Figure 2B) LVD PDFs appear somewhat similar for both waking and sleep stage 2. By contrast, we observe a clear distinction between waking and sleep stage 2 in terms of both the mean number of LVDs per segment, and the percentage of segments containing more than 50 LVDs, both shown by statistical comparison to be highly significant ( Table 1). The number of segments we were able to utilize per state of consciousness for both positive and negative LVDs is listed in Table 1.

Waking Differs from Sleep Stage 2 in the Extent of DSI Seen in the EEG LVD PDF Fluctuations
LVD PDFs were subjected to log transformation for Lomb-Scargle periodogram analysis (Figure 3) to look for evidence FIGURE 1 | Mean spectral analysis and autocorrelation function analysis for waking vs. sleep stage 2 EEG. Fourteen subjects had 8 min each of both waking and sleep stage 2 EEG data, in non-continuous 1 min segments. (A) Fast Fourier Transform (FFT) was performed on 8 × 1 min segments of both waking (black) and sleep stage 2 (blue) EEG for each subject, and the spectra were averaged across segments and subjects to generate the mean spectrum shown. Waking EEG has more power in the alpha (8-13 Hz), beta, and gamma band (∼20-40 Hz) ranges. (B) Autocorrelation function (ACF) with a maximum lag of 500 data points was performed on 8 × 1 min EEG segments of both waking (black) and sleep stage 2 (blue) EEG for each subject. Note characteristic oscillatory pattern for waking EEG ACF.
Frontiers in Human Neuroscience | www.frontiersin.org of DSI in the PDF as a function of LVD size. While the total mean area of statistically significant periodogram signals does not differ between waking and sleep 2 positive LVDs, there is a modest difference in total mean area for negative LVDs (Figure 3; Table 2). However, for both positive and negative LVDs, there were clearly more statistically significant periodogram peaks for waking compared to sleep stage 2 (Figure 3; Table 2). As a further test of this hypothesis, in order to confirm that this effect persisted for longer segment lengths, we isolated additional 5 min long continuous EEG tracings from the same subjects for both waking and sleep stage 2 (all 14 subjects, one 5 min segment for each state of consciousness). These were then processed using the same techniques of isolation of 2.5 SD LVDs, followed by PDF estimation, log transformation, and Lomb-Scargle periodogram analysis ( Table 3). Given that there is less data per subject in this analysis, it is unsurprising that there is now only a strong trend to significance for positive LVD areas (p = 0.11) and peaks (p = 0.06), with no difference for negative LVD areas (p = 0.28) and peaks (p = 0.25; Table 3). However, in all cases the average values for waking EEGs are greater than those for sleep 2 EEGs, supporting the patterns seen for the 1 min EEG tracing analysis (Tables 2, 3;  Figure 3).

Simulation Data and LVD DSI
In order to evaluate whether DSI is specific to waking EEG avalanches, we utilized mono-and multifractal simulation data which have similar fractal dynamics (as assessed by Hölder exponents) to waking EEGs (cf. Zorick and Mandelkern, 2013). Positive 2.5 SD amplitude deflections were isolated from BMS, fractional Brownian Motion with H = 0.2 (fBM 0.2) and log normal sigma 0.1 (LNS 0.1) multifractal series, and subjected to PDF estimation and log transformation followed by normalized Lomb-Scargle periodogram analysis as for EEGs (Figure 4). For perspective, these were plotted together with mean waking FIGURE 3 | Mean significant Lomb-Scargle normalized periodogram areas from PDF plots of log scale frequencies for positive (A) and negative (B) LVDs. Data represent mean normalized Lomb-Scargle periodogram plots of power > p = 0.01 across 14 subjects for 8 min of EEG for waking (black) and sleep stage 2 (blue) LVDs. * * p < 0.001 for mean number of significant peaks for waking vs. sleep 2 by linear mixed effect modeling. # p < 0.05 for total significant area for waking vs. sleep 2 by linear mixed effect modeling; see Table 2. and sleep stage 2 positive LVD normalized Lomb-Scargle periodogram data derived from a natural log scale (Figure 4). Significance level for normalized Lomb-Scargle periodogram peaks are shown with a horizontal line for the p = 0.01 significance level (Figure 4). As can be seen, normalized Lomb-Scargle periodogram data of 2.5 SD positive deflections from BMS, LNS 0.1, and fBM 0.2 do not approach significance, whereas the mean waking and sleep stage 2 LVDs clearly exceed the significance threshold (Figure 4). This indicates that DSI is not a general feature of fractal series with mean Hölder exponents of ∼0.2, but rather specific for EEG-derived LVDs (Figure 4).

Removal of EEG LVDs Reduces Alpha and Beta Spectral Power
In order to assess whether >2.5 SD EEG LVDs contribute to the spectral structure of EEG, we removed all EEG values containing 2.5 SD positive and negative LVDs from the 14  subjects' waking consciousness tracings, or about 4.2% of voltage readings per segment on average (Figure 5). For comparison, we included both complete waking EEG tracings, and waking EEG tracings where 5% of the voltage values were removed randomly from each segment, and waking EEG tracings where 5% of contiguous EEG data was excised from the middle of the segment (Figure 5). FFT was performed on these data, and spectral power in the alpha (8-13 Hz), beta (16-31 Hz), and gamma (35-45 Hz) frequency bands was summed for each of eight waking EEG segments from the 14 participants ( Figure 5; Table 4). Compared to complete waking EEG, random deletion of either 5% of the voltage values or a segment of 5% of contiguous values from the middle of the EEG tracing had minimal effect on the resulting mean spectral power in the frequency bands, whereas removing both positive and negative 2.5 SD LVDs significantly reduced mean spectral power in the alpha and beta bands ( Figure 5; Table 4). Therefore, 2.5 SD LVDs characterized by DSI contribute substantially to various spectral features in EEG.

EEG-Derived LVDs Exhibit DSI
It has been well-established that power-law dynamics of cortical neuronal avalanches differ little between in vitro and in vivo animal models and human MEG (i.e., Plenz and Thiagarajan, 2007;Ribeiro et al., 2010;Shriki et al., 2013). Given that  the EEG is likely determined by time-dependent dynamical neuronal avalanches, we reasoned that such structure in the waking ACF could indicate the loss of some measure of scale invariance in the neuronal avalanche distribution. The central finding of our study is that DSI is characteristic of the PDFs of EEG-derived LVDs in both waking and sleep stage 2 consciousness (Figure 3). Therefore, there is likely to be widespread coordination of the sizes of cortical neuronal avalanches in cortical brain activity, indicating a characteristic ''lacunarity'' in the power-law distribution (Sornette, 1998(Sornette, , 2006.

DSI is more Prevalent in Waking than Sleep Stage 2
Our study showed that the total number of significant Lomb-Scargle peaks of positive and negative 2.5 SD LVDs derived from 1 min EEG segments was greater in waking than in sleep stage 2, and the total significant Lomb-Scargle power was greater in negative LVDs for waking than sleep stage 2 (Figure 3, Table 2). These results strongly suggest that improved EEG datasets would demonstrate a stronger difference in many measures of DSI between waking and sleep 2. Given that sleep stage 2 EEG had less overall variance than waking, fewer segments of sleep stage 2 consciousness gave >50 LVDs per segment than in waking ( Table 1); for the purposes of this study, these segments were eliminated from further analysis, but this may have led to overestimating the level of DSI in sleep stage 2 compared to waking. While this is not proof, it is a strong inference that DSI is more prominent in waking, compared to sleep stage 2.

Limitations
We utilized a small, publicly-available database with limited demographic and clinical information, so certainly these results should be repeated on larger, more complete datasets. While we demonstrated the presence of DSI in EEG-derived LVDs, we cannot demonstrate at this time that the degree of DSI can completely account for regular structures in seen in EEG FFT or ACF analysis (Figures 1, 5). Our definition of EEG-derived LVDs (derived from a single lead) lacks a spatial component, which may limit the applicability of these results to other studies of neuronal avalanche dynamics. While we have observed a difference in DSI derived from EEG LVDs between waking and sleep stage 2, these differences may be characteristic only of sleep, and not characteristic of other unconscious brain states (i.e., anesthesia, coma). More definitive evidence that DSI of EEG-derived LVDs is attenuated in states of reduced conscious awareness will have to await analysis with better datasets.

DSI of EEG LVDs: a New Clue in the Study of Consciousness?
The study of consciousness has long fascinated scientists, and much has been learned about the phenomenology (Boly et al., 2013) and neurobiology (Edelman and Tonioni, 2000;Baars et al., 2013) of conscious mental states. Additionally, much recent work has gone into the role of LRTC in states of consciousness and brain pathology (Li et al., 2008;Nikulin et al., 2012;Blythe et al., 2014). However, a neural dynamical understanding of consciousness has been elusive. We believe that the finding that DSI is more prominent in conscious mental states could be an important observation that will lead a paradigm shift in the understanding of the neurobiology of consciousness, as it would be able to link the behavioral manifestations of consciousness with a specific statistical physical pattern of neuronal activation via EEG. These provocative findings, while preliminary, do raise the hypothesis that the degree of DSI in EEG LVDs is intimately linked to cortical functioning in consciousness. In this hypothesis, conscious brain states would be characterized by increased DSI as compared to unconscious brain states. However, this excitement should be tempered by the preliminary nature of these findings. Certainly many more studies in both in vitro and in vivo systems are needed to confirm and expand upon these observations in order to have a strong indication that LVD DSI is indeed important for the understanding of consciousness. Nonetheless, we find these early results compelling, and worthy of future investigation.

AUTHOR CONTRIBUTIONS
TZ came up with the initial ideas for this work. TZ and MM jointly performed all the experiments listed, discussed findings and interpretations, and wrote the manuscript.