Sleep spindles detection from human sleep EEG signals using autoregressive (AR) model: a surrogate data approach

A new algorithm for the detection of sleep spindles from human sleep EEG with surrogate data approach is presented. Surrogate data approach is the state of the art technique for nonlinear spectral analysis. In this paper, by developing autoregressive (AR) models on short segment of the EEG is described as a superposition of harmonic oscillating with damping and frequency in time. Sleep spindle events are detected, whenever the damping of one or more frequencies falls below a predefined threshold. Based on a surrogate data, a method was proposed to test the hypothesis that the original data were generated by a linear Gaussian process. This method was tested on human sleep EEG signal. The algorithm work well for the detection of sleep spindles and in addition the analysis reveals the alpha and beta band activities in EEG. The rigorous statistical framework proves essential in establishing these results.


INTRODUCTION
Oscillatory signal activities are ubiquitous in the biomedical signals [1].Multielectrode recordings provide the opportunity to study signal oscillations from a network perspective.To assess signal interactions in the frequency domain, one often applies methods, such as ordinary coherence and Granger causality spectra [2] that are formulated within the frame work of linear stochastic process.Electroencephalogram (EEG) is one of the most important electrophysiological techniques used in human clinical and basic sleep research.In 1979 Bar-low proposed linear modeling system which has a longlasting history in EEG analysis [3].The models are mainly considered as a mathematical description of the signal and less as a biophysical model of the underlying neuronal mechanisms.
In 1985, Frannaszczua et al. [4] proposed a model to interpret linear models as damped harmonic oscillators generating EEG activity based on the equivalence between stochastically driven harmonic oscillators and autoregressive (AR) models.There is a unique transformation between the AR coefficients and the frequencies and damping coefficients of the corresponding oscillators.In particular at times when the EEG is dominated by a certain rhythmic activity e.g. in the case of sleep spindles or alpha activity, on might expect, that this activity will be rejected by a pole with a corresponding frequency and low damping.This idea was the staring point of our analysis [5].
The sleep EEG is always not stationary.However, we demonstrated that the effects of non stationary become relevant only with scales longer than 1s [6].Therefore, short segments with duration of around 1s are sufficiently described by linear models.The non stationary in longer time scales might be rejected by the variation of the AR-coefficients and thus by the corresponding frequencies and damping coefficients.Based on the above considerations we propose an easy way to define oscillatory events.They are detected, whenever the damping of one of the poles of a 1s AR model is below a predefined threshold.
The method of surrogate data is a tool to test whether data were generated by some class of model.In 1992 the method of surrogate data proposed by Theiler et al. [7] is a general procedure to test whether data are consistent with some class of models.In order to test the hypothesis that the data are consistent with being generated by a linear system, the Fourier Transform (FT) algorithm is applied.Based on a example and the theory of linear stochastic systems we will show that this algorithm produce correct The statistical properties of a time series generated by a linear process are specified by the autocovariance function (ACF) or equivalently by its Fourier Transform, the power spectrum X (ω).The purpose of this study is described as even without a rigorous mathematical foundation it is possible to detect sleep spindles as well as alpha and beta activities from human sleep EEG.The theoretical framework of AR model and the procedure to generate the surrogate data are given in Section 2. The method is then tested on one simulation data set in Section 3. Section 4 describes sleep spindles detection using AR model with surrogate data approach.Results are discussed and summarized in Section 5.

METHODS
In this section, we start with AR model of order p and then proceed to outline the procedure for generating surrogate data.

AR Model
The Parametric description of the EEG signal by means of the AR model makes possible estimation of the transfer function of the system in the straight forward way.From the transfer function it is easy to find the differential equation describing the investigated process.Our detection algorithm is based on modeling 1s segments of the EEG time series using autoregressive (AR) models of order p.From the AR (p)-model Applying the Z-transform to Eq.1 we obtain: where X(z)-the Z-transform of the signal x E(z)-the Z-transform of the noise.
If the system is stable, there exists A -1 (z) and we get In the z domain this filter is expressed by the Eq.4 where A -1 (z) is the transfer function.Denoting it by H(z) and writing if explicitly we obtain: Multiplying numerator and denominator p z we get Factorizing the denominator gives the formula where Using the above formula to estimate the frequencies and damping coefficients ( We assume that there are only single poles of H (Z) which can be written in the form For the single pole coefficients can be found according to the formula: By means of the inverse transform -z -1 from the Eq.8 the impulse response of the system: h (n) can be found.Since from the properties of the z -1 transform we know that we obtain: If the sampling interval ∆t was chosen according to the Nyquist theorem we can express the impulse response as continuous function and write it in the form: where Laplace transform of the Eq.11 which corresponds to the transfer function H(s) of the continuous system is given by the formula: The above expression can be obtained directly from Eq.8

SciRes Copyright © 2009 JBiSE
by means of the integral transform z.Eq.12 can be written as a ratio of two polynomials of the order p-1 and p The polynomial coefficients can be readily calculated from and .This form of the transfer function was found also by Freeman 1975.It leads directly to the differential equations describing the system.The transfer function is the ratio of the Laplace transform of input y(t) and output x(t) functions: Since variable s corresponds to the operator dt d , from Eq.13 and 14 we get: In this way we have obtained the differential equation describing the system which is free of the arbitrary parameters.Its order is determined by the characteristic of the signal and may be found from the criteria based on the principle of the maximum of entropy.

Surrogate Data Method
Generally, a surrogate data testing method involves three ingredients: 1) a null hypothesis; 2) a method to generate surrogate data; and 3) testing statistics for significance evaluation.The null hypothesis in the present study is that investigated data from a linear Gaussian process.If the null hypothesis is rejected, then we conclude that the data are either non-Gaussian or come from nonlinear process.The surrogate data are generated in such a way that it is Gaussian distributed but has the same second order spectral properties (in the bivariate, auto spectra) as the original data.The testing statistics the amplitude of the peridogram we will give the steps for generation the surrogate data and the theoretical rationale behind the steps [8].
Consider two zero mean stationary random processes x(t) and y(t).These processes may or may not be linear Gaussian processes.Their one-sided auto spectra S xx and S yy can be estimated [9].
where T is the duration of the data, the expectation is taken over multiple realizations and X k and Y k are the Fourier Transform of x(t) and y(t).Now we consider how to generate two zero-mean linear Gaussian processes ) (t x and ) (t y which have the same second order statistical properties as x(t) and y (t).Let k X and k Y  expressed in real and imaginary parts, ,  [10].(18) According to these relations the covariance matrix (∑) for the real and imaginary parts of and

Linear Bivariate AR Model Driven by Gaussian White Noise
For easy implementation, we summarize the earlier analysis as follows.
Step 1) Estimate and for original data according to (16) xx S yy

S
The model is written as Step 2) Calculate Covariance matrix ∑ according to (19) Step 3) Draw values (realizations) of R , I X  X  , R Y  and I from the Gaussian processes with zero mean and covariance matrix (∑) for each frequency.This can be done for example Matlab function based on the ziggurat method [11].The one sided power spectra xx and yy for both original (solid curve) and one set of surrogate data (dotted curve) are shown in Figure 1.It is seen that the surrogate data's spectra nearly match well with that of the original data.This is an expected result.The maximum amplitude for x(t) original power spectra is 30.1891[Figure 1(a)] and that for surrogate data x′(t) power spectra is 32.6881.Similarly the maximum amplitude for y(t) original power spectra is 28.5673 [Figure 1(b)] and that for surrogate data y′(t) power spectra is 29.8921.Notice that the power spectra amplitude is higher between 0 to 30 Hz than other regions.Namely, the estimation variances are proportional to the power spectral amplitudes as those frequencies.Considering that the surrogate data contour plots are computed based on a randomly selected data set among P=1500 available, these maximum value comparisons suggests the known fact that there are no nonlinear or non-Gaussian components in the original data.The histogram and Gaussian fit of input data x(t) and y(t) and surrogate data x′(t) and y′(t) are generated by ziggurat algorithm shown in Figure 2. The PDFs are obtained from these parameters

S S
Step 5) Repeat Steps 3) and 4) to generate multiple realizations of the surrogate data.
R and I , could be drawn from Gaussian distribution with zero mean and covariance matrix ∑ x.This simplified method is similar to the method proposed by Timmer [12].The probability density function for the extreme value distribution (type I) with location parameter µ and scale parameter is σ [13].
The exact PDF is determined once µ and σ are known.From this distribution, one can determine the threshold for any desired significance level (i.e.-p value).

SIMULATION EXAMPLE
We have performed simulation studies to test the effectiveness of the method proposed before.Matlab "Signal processing and Spectral analysis toolbox" is used in our analysis.

SciRes Copyright © 2009
JBiSE according to Eq.20 and are plotted along with the histogram in blue color in Figure 3. Notice that the PDF here are multiplied by the total number of surrogate data sets (1500) to match the histogram.From the PDFs, the threshold for a significance level of p < 0.005.By com-paring the original data's maximum power spectra values with the thresholds, it can be seen from the null hypothesis that the original data coming from a linear Gaussian processes cannot be rejected, a theoretically expected result.We have also performed a study to examine whether fitting an extreme value distribution to the empirical histogram is a viable approach.The fitted model parameters µ and σ versus the number of surrogate data sets used are plotted in Figure 4.It can be seen after around 500 surrogate data sets, the estimation becomes linear.

SLEEP SPINDLES DETECTION
A seven minutes recording of 9 channels of EEG (C3, A2, O1, O2, C4, P3, P4, F3, and F4) was used in our test (16,207  The EEG derivation C3A2 and C4P4 was analyzed shown in Figure 6.Oscillatory events are detected, if the damping coefficients r k falls above a pre defined threshold and hence r k , exceeds the corresponding threshold.In practice, we use two thresholds, a lower one, r a , to detect candidate events scanning the EEG with non-overlapping 1-s segments.When r a is crossed we go back to the previous segment and use a smaller step size of 1/16s (overlapping 1-s segments).If r k exceeds a second threshold r b > r a the beginning of an oscillatory events is detected.The oscillatory events are terminated by the time r k is lower than r b for the last time before it falls below r a .The frequency and time at the position of the maximal value r max are considered as the frequency and occurrence of the event respectively.Spindles or Sigma waves are a poor indicator for sleep onset.Some subjects do not have discernable sigma waves.These spindles are more prominent in Stage 2 than Stage 1.In the present analysis the order of the AR model was set to p=8 and the threshold \value set to r a =0.75 and r b =0.85 parameter were chosen in such a way, that clearly visible sleep spindles were reliable detected by the algorithm.Sleep stages were visually scored according to standard criteria in Figure 7.

SciRes Copyright © 2009 JBiSE
In Figure 7 shows the detected events from 2 recordings.The distributions show modes in four stages: in the delta (±1:0-2 Hz), in the fast delta (±2: 2-4 Hz), in the alpha (®: 8-12 Hz) and sigma (sleep spindles: 11.5-16 Hz) bands.These modes are also evident in the distribution of the events in the different sleep stages (Figure 7).Alpha waves predominate in waking, spindles in Stage 2 while the occurrence of delta waves decreases from Stage 2 to Stage 4 (deepening of sleep).Delta waves are the prevalent events in REM sleep and Stage 1.Note that they also occur in Stage 2 with almost the same incidence.Events in the alpha frequency range correspond to continuous alpha activity during waking and to small amplitude, sometimes spindle-like activity in NREM sleep.However, in particular during NREM sleep Stages 3 and 4, slow waves are often not detected as oscillatory events because they yield a relaxatory pole (frequency zero) in the AR-model.Table 2 provides detection of different frequency bands and their amplitude calculated from the Fast Fourier Transform.The two most important characteristics of EEG elements are frequency and amplitude.Frequency is inverse value of duration of EEG segment.The frequency range is divided in to four bands: beta (12Hz and higher), alpha (8-12Hz), theta (4-8Hz) and delta (0, 1-4Hz).Amplitude of EEG is taken as the peak to peak value.After the magnitude of amplitude EEG signal is divided into low-voltage, middle and high-voltage EEG.For delta, theta and alpha bands these values are following: 10-30 µV for low voltage, 30-70 µV for middle voltage and above 70 µV for high voltage EEG.For beta band, the values are lower: below 10 µV low voltage, 10-25 µV middle voltage, and above 25 µV high voltage EEG.
Beta activity is typical during wakefulness.Alpha occurs in relaxed state with eyes closed.Theta waves are dominant in normal wake state in children; in adult people theta activity appears only in small amount especially in drowsiness.Delta waves are present in deep sleep; in healthy people do not exist in wake state.Sleep spindles are rhythmic activity in 12-14Hz frequency range and amplitude below 50 µV.The power spectra for the original (solid line) and the surrogate (dotted line) data are shown in Figure 8. Three peaks around 1-14Hz can be seen, indicating the presence of synchronized activities at these frequencies.According to the tradi-tional classification first peak belongs to the delta band (1-4Hz) and the second peak belongs to the alpha band (8-12Hz) and third peak belongs to the beta band (>12Hz).Past research has examined whether such oscillatory neural activities self-couple and couple each other in a nonlinear fashion [17].We study this issue with the new surrogate test method.Table 3 provides comparison between E. Olbrich et al. and our study.However, a few spindles are not detected by the algorithm with the chosen parameters, i.e. the maximum of r k remains smaller than the threshold r b .Therefore, in our method proposed is using a surrogate data approach needs to be further optimized for more sensitive spindle detection.

. DISCUSSION AND SUMMARY
In this study sleep spindles were detected in 1s overlapped segments, wide 0.5-16Hz band containing delta and theta was used for better eye movement separation.Alpha activity was excluded as it would increase auto covariance.With consensus scoring we would expect the agreement to improve [5].The cost benefit of this laborious task would have been low.There was no manual artifact handling and it is possible that some misclassifications are the results of the artifacts undetected by the corresponding threshold r k and were generated surrogate data with ziggurat method.In this algorithm there are two different thresholds but as it can be seen from Table 3, the amplitude criterion was most important (i.e.mostly amplitude in REM sleep is less than 50 µV and NREM sleep is greater than 100 µV).Schwilden et al. [14] reported that 90% of the EEG did not show any nonlinearity.They suggested that only under some pathological conditions, such as epilepsy, the brain signals manifest nonlinear effects.Other studies [15,16], have shown that non linear characteristics exist between theta and gamma EEG signals during short term memory processing.To settle these debates, one needs carefully constructed statistical tests.The method proposed here represents our effort in this direction.We will contrast our method with other often applied techniques in this area.

Comparison with other Surrogate Data Method
Multiple ways are currently in use to generate surrogate data sets to test the significance of the sleep spindle detection.In one method, the Fourier Transform (FT) is estimated for each segment and the phase of the Fourier Transform is randomized without changing the Fourier amplitude.The result is then inverse Fourier Trans-formed to generate a surrogate time series [14,17].While the phase randomization destroyed nonlinearity and leads to linear Gaussian process [18], the Fourier amplitudes are random variables for a stationary process as well.By keeping them constant, one loses an important degree of freedom [12].In contrast, our method, in which both amplitudes and phases of the Fourier Transform are randomly generated, produces surrogate data sets that are explicitly Gaussian, linear and share the same second order statistics with the original data.Another method, called amplitude adjusted Fourier Transform (AAFT), has been used in recent studies [19].
AAFT was designed to test the null hypothesis that the observed time series is a monotonic nonlinear transformation of a linear Gaussian process.This method has the same problems as the previous phase randomization method.In addition, the generated surrogate data sets usually do not have the same power spectrum as the original data, leading to false rejections when the discriminating statistics are sensitive to second-order statistical properties [20].

Final Remarks
We make several additional remarks regarding our method.First, generating the null hypothesis distribution for the test statistic is a very time-consuming process.
The application of the extreme value theory mitigates this problem.Our examples show that the probability distribution function based on the extreme value theory fits the empirical distribution well.In agrees with that from the empirical histogram.Secondly, when the original data with length L need to be zero padded to length N (N > L), the surrogate data we generated have length N. To make sure that the surrogate data still have the same power spectra with the original data, the surrogate data need to be normalized by L, instead of the real data length N. Third, the periodogram method is used to estimate the power spectra.Consequently, all of the limitations associated with the periodogram method were inherent in our method, including poor frequency resolution for short data, fourth when we generated surrogate data set using ziggurat method does not generate random numbers effectively in tail regions because Matlab does not execute greater than 500 function calls [11].Finally, This can be determined empirically by repeating the test many times on different realization for the data [20].We will consider this issue as part of our future research.In summary, detection of sleep spindles using the surrogate method proposed in this paper is shown to give accurate results when applied to test the significance of power spectral amplitudes.It is based on solid statistical principles and overcomes some weaknesses in previous methods for the same purpose.It is expected to become a useful addition to the repertoire of nonlinear analysis methods for neuroscience and other biomedical signal processing applications.
value of the sampled signal at the moment n n  -Zero mean uncorrelated white noise process.

Y
 and I Y  are then generated by sampling from a Gaussian distribution with zero mean and covariance ma-trix∑.Once R X  , I X  , and are generated for each frequency, the surrogate data and

Figure 1 .
Figure 1.One sided power spectra of a) x(t) and x′(t) b) y(t) and y′(t).The solid curve indicates the result from the original data and the dotted curve indicates the result from one of the 1500 surrogate data sets.

Figure 2 .Figure 3 .
Figure 2. Histogram and Gaussian fit of the a) original and b) surrogate data.

Figure 4 .
Figure 4. µ and σ versus the number of surrogate data sets.
trials) shown in Figure 5. Precise numbers of sleep stages and standard characteristics derived form hypnograms are in Table 1.Number of sleep stages (17s records) are in the first column; data in the second column are the average values over 5 subjects (C3A2, C4P4, F3C3, P3O1, P4O2) related to total sleep time, the beginning of sleep is set as the first appear-ance of the sleep Stage 2. Sleep efficiencies is the ratio of time spent in Stages 1-4 or REM sleep to the whole sleep time, where also movement time and some awakenings during the night are included; in sleep medicine this parameter discriminates some sleep disorders.Surrogate data were generated following the procedure in Section 2.

Figure 8 .
Figure 8.One sided power spectra estimated recording from P3O1 electrode on 1s period.The solid line is for the original data and dotted line is for surrogate data.

SciRes Copyright © 2009 JBiSE distribution
of time series and therefore might not generally yield the correct distribution of the test statistics.The surrogate data method differs from a simple Monte Carlo implementation of a hypothesis test in that it tests not against a single model, but a class of models, i.e. linear systems driven by Gaussian noise.The idea is to select single model from the class on the basis of the measured data x, and then do a Monte Carlo hypothesis test for the selected model and the original data.
In practice, S xx and S yy estimated from the data.The question is how to draw Gaussian variables for R X  , I X  , R Y  and I Y  such that and R X  , I X  , R Y  and I Y  and xx S and are as follows

Table 1 .
Number of sleep stages (first column) and sleep efficiency and average percentage of sleep stages during the sleep (mean ± standard deviation).

Table 2 .
Different frequency band detection and the corresponding true value determined from one sided Power spectra.
able 3. Comparison between AR model and AR model with surrogate data approach.JBiSE the discrimination power of our statistical test is unclear.