Paper—Empirical Characterization of the Temporal Dynamics of EEG Spectral Components Empirical Characterization of the Temporal Dynamics of EEG Spectral Components

The properties of time-domain electroencephalographic data have been studied extensively. There has however been no attempt to characterize the temporal evolution of resulting spectral components when successive segments of electroencephalographic data are decomposed. We analyzed resting-state scalp electroencephalographic data from 23 subjects, acquired at 256 Hz, and transformed using 64-point Fast Fourier Transform with a Hamming window. KPSS and Nason tests were administered to study the trendand wide sense stationarity respectively of the spectral components. Thereafter, the Rosenstein algorithm for dynamic evolution was applied to determine the largest Lyapunov exponents of each component’s temporal evolution. We found that the evolutions were wide sense stationary for time scales up to 8 s, and had significant interactions, especially between spectral series in the frequency ranges 0-4 Hz, 12-24 Hz, and 32128 Hz. The spectral series were generally non-chaotic, with average largest Lyapunov exponent of 0. The results show that significant information is contained in all frequency bands, and that the interactions between bands are complicated and time-varying. Keywords—Trend stationarity, Wide sense stationarity, Chaoticity, Largest Lyapunov Exponent


Introduction
A typical EEG data acquisition pipeline can be modelled as a time-domain operation, with data acquired from multiple electrodes at a fixed sampling rate leading to a time series, [ ] = 1 , 2 The assumption that the J sub-sequences are contiguous and mutually exclusive is common in practice, and simplifies analysis. In many cases however, it can be relaxed, and neighboring sub-sequences of [ ] are allowed to overlap.
The transform [⋅] is most commonly the Fast Fourier Transform (or similar variants) for which = 2 . The ℎ row of [ , ] represents the m spectral components for the ℎ sub-sequence of [ ] , while the ℎ column represents a temporal sequence comprised of the ℎ spectral components from all J sub-sequences of [ ]. Further processing using the whole of [ , ] can be computationally exacting, especially considering the fact that most problems call for data from multiple EEG electrodes to be used. Consequently, there is often interest in dimensionality reduction in some form or the other, even after spectral decomposition of temporal EEG data.
One common approach is to reduce the elements in each row of [ , ] to some number << of elements, by dividing the row into K subset, and replacing each subset with some representative statistic. Another common method of dimensionality reduction is to apply a band-pass filter to extract a single sub-sequence from each row of [ , ], which results in discarding some elements in the row [33], [39]. Neither of the foregoing approaches is optimal. Representing a band with a single statistic is likely to result in failure to detect intra-band dynamics, which may be significant, depending on the task at hand. On the other hand, applying band-pass filtering greatly increases the risk that information-rich spectral components will be discarded. For instance, a significant number of frequency domain approaches to epileptic seizure detectors apply low-pass filters to eliminate components above 25 Hz, in order to limit electromyographic (EMG) artefacts [13]. However, current thinking on the roles of ripple (80 -250 Hz) and fast ripple (250 -500 Hz) high frequency oscillations in the generation of epileptic waveforms Wang et al. [40] strongly suggests that such low pass filtering eliminates important information along with the EMG artifacts.
Specialized dimensionality reduction techniques such as principal component analysis (PCA) ought to be superior to either of the preceding approaches. Before those more sophisticated techniques are applied to reduce [ , ] however, the nature of the dynamics underlying [ , ] must be understood. This study examined the characteristics and evolution of spectral components of EEG time series after spectral decomposition. To the best of our knowledge, there has not been a previous study to characterize the individual and collective temporal evolution of EEG data after spectral decomposition. The rest of the paper is organized as follows. Section 2 provides an overview of time series methods and testing, identifying an organizing framework within which a limited set of measures might be used as a basis for examining the characteristics of EEG time series spectral components. In Section 3, three time series measures were adopted for a number of tests on a selected EEG dataset. Results are presented and discussed in Section 4, with concluding remarks following in the last Section.

Brief Review on Time Series Analysis of EEG Data
Time series arise across all scientific disciplines. Hundreds of measures for time series analysis have been developed over the last few decades with few notable efforts to develop an organizing framework. To the best of our knowledge, there has only been one large, systematic organizing framework for cross-disciplinary classification of time series datasets and methods [11], in which applied over 8651 distinct time series analysis operations to a set of interdisciplinary time series, and elicited a framework for organizing them into small groups. A minimum of four clusters were identified into which all time series measures can be organized. The clusters corresponded to measures of stationarity and scaling, correlation and predictability, complexity, and measures of non-linear or dynamical evolution.
Measures representative of each of Fulcher's four clusters have been applied towards characterizing EEG time series over the years. Some of the earliest studies on EEG time series analysis focused on stationarity. In the early days, the sign-based RUN test and trend test [5], [18] were used to test varying sample segments for their stationarity. More sophisticated tests for testing stationarity are now available. Unit root tests have been explored to test EEG timeseries [22]. Tests such as the Augmented Dickey-Fuller (ADF) test [6] [30] and Kwiatkowski-Phillips-Schmidt-Schin (KPSS) test for trendstationarity or mean-stationarity [24] are the most common unit root employed in this domain.
With the rise in interest in connectivity across brain regions, predictability have largely been replaced by connectivity [25]; and multivariate measures have become more important, allowing activities in different parts of the brain to be compared [3]. The initial use of cross-correlation has however been put aside [14] so that more focus is now enjoyed by methods based on coherence, a measure of the correlation between two signals as a function of frequency [41], [34] or phase synchrony [19], [2]. Coherence can be contaminated by volume conduction and reference electrodes, and effectively only measures linear relationships between two time series [21] and so, phase synchrony is generally considered superior.
The link between entropy and complexity can be made by considering the fact that Kolmogorov complexity and the Shannon entropy are equal up to some constant limit for any recursive probability distribution [38]. In an information theoretic sense therefore, complexity is often synonymous with uncertainty or disorder, and is most often estimated using one of a large number of mathematical methods bearing the name "entropy". Complexity measures have been applied to EEG time series as a way of determining the diversity of cooperating underlying cortical populations and the extent and complexity of their interactions.
EEG data have been shown to be non-linear, and a lot of work has been done on characterizing the nature of their time evolution. Several methods have been employed in testing dynamical systems for sensitivity to initial conditions i.e. chaoticity. Methods such as the largest lyapunov exponent (LLE), whose sign indicates chaos has been studied extensively. Recently, methods that do not analyze time series using the phase state, but rather analyze the time series directly in the time domain have proven successful in characterizing various chaotic systems.

Overview
Data generated by successive spectral decomposition of sub-sequences of an original time series [ ] remain time-domain data in a sense, as they can be organized into new time-indexed series. Hereafter, each of the J-element-long time series generated after successive spectral decomposition of original EEG time series data will be termed a spectral component time series (SCTS). Also, for ease of reference, the SCTS generated from the ℎ frequency bin of successive spectral decompositions of data from electrode q will be represented with the notation . The goal of this study was to characterize the temporal evolution of these spectral component time series, leading to an understanding which would, among other things, allow correct assumptions to be made during further processing of frequency-domain data. Fulcher's taxonomy of time series methods provided a useful basis for selecting a limited number of measures with which to characterize the emergent spectral component time series.

EEG dataset and pre-processing
The Childrens' Hospital, Boston -Massachusetts Institute of Technology (CHB-MIT) scalp EEG dataset [33,12] was used for this study. CHB-MIT subject data are organized as sessions, with each session typically one or four hours long, and stored as a single European data format (edf) file. All but one subjects have a single set of recordings composed of multiple sessions, while the last subject has two sets, with the second set acquired 1.5 years after the first. In all, the records include over 800 hours of recording containing 182 seizures. Seizure start and stop times are included in bundled text files along with montage information. We only considered normal EEG data in this study.

Testing for stationarity
There are many tests for mean or wide-sense (WS) stationarity. Certain tests for mean stationarity work even when a deterministic trend is superimposed on the time series. This is called trend stationarity. A common test for trend stationarity is the KPSS procedure, which tests a null hypothesis that a time series is trend stationary against the alternative hypothesis that it possesses a unit root. The KPSS test procedure assumes the framework of a random walk process with time-dependent drift component. That is, if the time series can be expressed as where the first term on the right hand side of Equation 6 represents a deterministic trend, r[t] is a random walk and u[t] a stationary error process. A test of the null hypothesis ( 0 : = 0) is now evaluated against an alternative hypothesis that the time series is non-stationary due to the presence of a unit root process.
Some of the distributions over SC times series in this study were heavy tailed, and consequently, the Nason test was adopted [29]. We partitioned 2 1 from Subject 6's first electrode (FP1) into 16-sample long blocks. Since 64 -bit FFT was caried out on EEG data acquired at 256 samples per second, each sample in 2 1 resulted from the FFT of a 0.25 s long EEG time-domain data. Each 16-sample block therefore corresponded to 4 s worth of spectral data. Both KPSS and Nason tests were administered to every block of the data, after which , the ratios of number of stationary blocks versus total number of blocks were obtained. This procedure was repeated for each of the time series, , ∈ {1,2, ⋯ ,16}, ∈ {1,2, ⋯ ,33} comprising Subject 1''s non-seizure data. This resulted in a 16 x 33 array containing the proportion of stationary 4-second blocks for each spectral component per electrode for Subject 1. Thereafter, the procedure was repeated for other block sizes. Based on previous studies on EEG stationarity, block sizes of 16, 32, 64, 128, 256, 512, 1024, 2048, and 4096 samples were used, and in each case, a 16x33 array was generated.
Previous studies on the stationarity of time-domain EEG data, reported as much as 98% agreement between results on data from different scalp sites. Consequently, they typically based their methodology on single-electrode data series [5], [36]. For this study, a similar agreement was anticipated for SCTS from different electrodes. To verify this, we tested a hypothesis that the distributions over the combined time series generated from any two electrodes were identical. For each of the 16  Finally, a key concern of this study was investigated, namely, the extent of similarities between the stationarity of different time series 1 , 2 , . . . , 3 for any particular electrode . Data acquired from an EEG electrode result from a superposition of electrical activities of large populations of cortical neurons directly under, or adjacent to the electrode. Different activities of brain populations are believed to give rise to electrical activity with signature frequencies that tend to cluster in bands. It might therefore be anticipated that the patterns of the stationarity of the data from different spectral component might demonstrate correlations. To test this, we used the 33 SCTS generated from the P7 electrode of Subject 10 (subject and electrode were both randomly selected). Dividing the entirely of each SCTS into a number, of blocks, with the length of each block corresponding to samples, we determined the stationarity of each block. This was done for both trend and WS stationarity. By assigning a "1" to each stationary block and a "0" to every non-stationary block, the 33 SCTS gave rise to 33 Bernoulli processes. The task of evaluating the similarities between any two SCTS consequently became a test of the proportions of two binomial distributions.

Testing for dynamic evolution
There is no universally accepted definition for systems or data that are chaotic, however, [35] gives a practical working definition of chaos as an aperiodic long-term behaviour in a deterministic system that exhibits sensitive dependence on initial conditions. Lyapunov exponents provide a direct measure of a systems sensitivity to initial conditions by quantifying the exponential rates at which neighboring orbits on an attractor diverge or converge as the system evolves in time.The sign of the LLE indicates chaos [17].
Adopting the pipeline of methods established by Krakovská et al. [23] which successfully identifies the optimal embedding parameters for calculating the largest lyapunov exponents (LLE), we use the average mutual information algorithm introduced by Fraser and Swinney [9] in finding the time delay and the false nearest neighbour algorithm [20] for identifying the optimal embedding dimension and finally the direct method developed by Rosenstein et al. [32] is used to estimate the LLE.
In the reconstruction of the phase space for the estimation of the largest Lyapunov exponent (LLE), it is necessary to identify the optimal reconstruction delay and embedding dimension . The embedding theorem introduced by Takens [37], Mañé [27] enables the 'reconstruction' from a univariate time series ( ) and its delayed variants* to form a -dimensional vector described below: We used the average mutual information [9], [10], [16] to compute the reconstruction's time delay . Although, there have been several recent methods like the neural networks method [26] and the recurrence networks technique [15] to select the optimal , these methods become more computationally expensive as the number and length of the time series being analysed is increased. False nearest neighbour (FNN) was used in determining the optimal embedding dimension [20]. For each time series, the % of false nearest neighbours was calculated iteratively for each dimension till the preset maximum dimension is reached. The optimal dimension is then chosen to be the least possible dimension with acceptable %FNN.
Due to the computational complexities attached to estimating the full Lyapunov spectrum, many algorithms have been introduced to estimate the LLE [42], [32], a key measure of chaoticity. With its robustness to short, noisy and non-stationary data [8], [4] the technique proposed by Rosenstein et al. [32] outperforms all other usable techniques for this particular case and was thus utilised in this study.

4
Results and Discussion Table 1 presents the average for all 33 SCTS for Subject 6 as block size is increased from 4 s to 1024 s. The values for each SCTS appear to start from low values for 4 s block sizes before rising above 0.9 for 16 s block sizes. Thereafter, there is a steady drop in as block sizes increase. A known deficiency of the KPSS test is that it may have high rate of Type I errors, particularly for time series lengths of 25 and below for certain values of a test parameter. Its power is generally good above 50 samples [7]. Block sizes of 8 s correspond to 32-element long time series, while 16 s block sizes imply 64-element long series. Results lower than 16 s may therefore be discarded. In that case, the KPSS test results suggest that trend-stationarity blocks become less frequent as larger block sizes are considered. Table 2 presents the average following the Nason WS stationarity tests on Subject 6's non-seizure data. There are fewer WS stationary blocks than there were trend-stationary blocks at almost every block size. Approximately half of all 32 s blocks were likely WS stationary, as seen in Table 2. The results of the wide sense stationarity tests are consistent with results from previous studies on time domain EEG signals. For example, they match very well with some of the results of Cohen and Sances [5], especially its Run test of half periods. For that test, at 5% level of significance, the null hypothesis of stationarity was not rejected for 91% of 12 s blocks, and 17 % of 128 s blocks. McEwen and Anderson [28] found out that up to 90% of 4 s blocks were likely WS stationary, while 70-80% of 16 s blocks were likely to be WS stationary. Two-sample KS tests were conducted to assess the extent of similarity between time series generated from different electrodes. For this test, , the proportion of WS stationary blocks for the Subject 6's data at different block sizes were used. The pvalues of the respective KS tests between all pairs of electrodes { , } , ∈ {1,2, . . . ,16}, ≠ are presented in Table 3. The table reveals that for this subject, the distribution over the data from any one electrode is essentially indistinguishable from the distribution over data from any other electrode 87.5% of the time, since the null hypothesis of the KS test could not be rejected for 105 out of 120 pairs of electrodes at 5% significance level. The 87.5% agreement across electrode data is comparable to previous studies such as [5], which found 98% agreement between electrode data. This provided a justification for utilizing single electrode data for additional tests where warranted.

Stationarity tests
The brain is a very complex system with many highly specialized subsystems. It is to be expected then that the instantaneous electrical activities under different electrodes might differ. Over longer time spans however, previous studies discovered similarities in the stationarity of data of the same length acquired from different electrodes. The level of agreement found in the current study (87.5%, as seen from Table 3) was however lower than previously reported (for example, in [5]). This may be due to intersubject differences, or to errors introduced by windowing, spectral decomposition, or differences in statistical tests.  While it has to be noted that the results above are subject to the low spectral resolution caused by use of 64-point FFT as well as computational errors, it should still be apparent that the interactions depicted by the charts are the banding of spectral correlates of neural activity that have been recognized in neuroscience as , , , and bands. This fact in turn leads to an important observation. Repeating the same series of tests with trend-stationarity rather than WS stationarity fails to reveal any of those interactions. We believe that this is a confirmation that in using N-th order stationarity for EEG and possibly other physiological signals, N should be at least 2 (WS stationarity). While mean-or trend-stationarity might have some value from a purely statistical standpoint, these results suggest that they may not correlate to any significant degree with underlying physiological dynamics. Studies based on first order stationary tests such as the ADF and KPSS are most affected by this observation.

Tests of dynamic evolution
The LLE for each SCTS per electrode for non-seizure data for all subjects were computed. For each subject, the edf file corresponding to a 1-hour session was processed. A random selection of the results is presented in Figure 1. Data series from most electrodes exhibit positive Lyapunov exponents in at least one spectral component. However, the superimposition of LLE for all electrodes allows a pattern to emerge wherein the Lyapunov exponents collectively cluster about zero. The immediate implication of these results is that the individual spectral components of decomposed EEG data do not appear to evolve chaotically. There has recently been some controversy about whether EEG time series are chaotic. While studies such as [1] had found evidence of chaotic dynamics, the introduction of surrogate data methods suggested that at least part of the "chaotic" behavior might be more correctly adduced to noise. EEG data are known to be noisy, with contamination arising from both the instrumentation pipeline and the the physiological processes that give rise to the scalp potentials themselves.Interpret the Lyapunov exponents of EEG time series in the context of chaotic nonlinear evolution has therefore attracted some controversy in recent times. In fact, it has been said that in the case of highly noisy and possibly highdimensional EEG data, the chaotic measures probably do not measure chaos any more, but reflect macroscopic statistical properties of the studied data [31]. Positive Lyapunov exponents may therefore not convincingly establish chaoticity. Nevertheless, a nonpositive LLE should result in the hypothesis of chaotic dynamics being rejected with some confidence. In Figure 1, some individual SCTS exhibit positive LLE for certain electrodes. This can be attributed to the effects of noise. The collective clustering about zero evens out the random effects of noise in all spectral bands, and allows a true picture the chaotic dynamic evolution of the EEG data to emerge.
A key takeaway from the series of tests in this study is that at some scales, EEG spectral components evolve distinctly, while at others, there are interactions and correlations with other spectral components. For example, from the WS stationarity tests, choice of different block lengths made the interactions between spectral components more or less apparent. The interactions between spectral component time series and the fact that there appear to be changing interactions within each band is an important finding from this study. For example, it suggests that studies in which dimensionality reduction is predicated on representing bands of frequencies with single statistics are likely to lead to information loss. Likewise, the results from tests of complexity show that all spectral components hold information, and in fact, there are interactions between the higher spectral components. This indicates the error of chopping off higher spectral components for dimensionality reduction or in a bid to increase signal to noise ratio.
We also found no evidence that the evolution of the individual spectral components in EEG data are chaotic. Many studies which studied time domain EEG data have found positive LLE, and concluded that this indicated chaotic dynamics. Even with the question of the effect of noise on the validity of equating positive LLE with chaotic evolution, it is obvious that they measure "something", because classifiers based on LLE in EEG have shown statistically significant effects. This finding however shows that for the individual spectral components after EEG data is decomposed, the question of interpreting LLE is moot; the series do not appear to be chaotic.

Conclusion
Using tests of stationarity, and largest Lyapunov exponents, this study investigated the temporal evolution of each spectral component of EEG data after spectral decomposition with a transform such as the Fourier transform. Some of our findings have implications for the selection of appropriate frequency domain methods that can be affected by the characteristics of the spectral time series. For example, we found that such data can be assumed to be WS stationary for segments up to 8 s. Also, higher spectral components contain useful information, and should not be discarded. The results of this study will lead to improved frequency domain methods for processing EEG data.

Acknowledgement
The workstation used for this research was donated by the Saffioti Family Foundation, Wollongong. It contained a Titan Xp GPU donated by the NVIDIA Corporation.