A new method for detecting state changes in the ... - Wiley Online Library

Paper 30 Disc J. Sleep Res. (1998) 7, Suppl. 1, 48±56

A new method for detecting state changes in the EEG: exploratory application to sleep data MARTIN J. MCKEOWN 1 , COLIN HUMPHRIES 1 , PETER ACHERMANN 2 , ALEXANDER A. BORBEÂ LY 2 and TERRENCE J. SEJNOWSKI 1 , 3 1

Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, CA 92186±5800, USA, 2Institute of Pharmacology, University of Zurich, Winterthurerstr. 190, Zurich CH-8057, Switzerland, 3Department of Biology, University of California, San Diego, La Jolla, CA USA 92093, USA Accepted ???????

SUMMARY

A new statistical method is described for detecting state changes in the electroencephalogram (EEG), based on the ongoing relationships between electrode voltages at different scalp locations. An EEG sleep recording from one NREM-REM sleep cycle from a healthy subject was used for exploratory analysis. A dimensionless function defined at discrete times ti, u(ti), was calculated by determining the log-likelihood of observing all scalp electrode voltages under the assumption that the data can be modeled by linear combinations of stationary relationships between derivations. The u(ti), calculated by using independent component analysis, provided a sensitive, but non-specific measure of changes in the global pattern of the EEG. In stage 2, abrupt increases in u(ti) corresponded to sleep spindles. In stages 3 and 4, low frequency (& 0.6 Hz) oscillations occurred in u(ti) which may correspond to slow oscillations described in cellular recordings and the EEG of sleeping cats. In stage 4 sleep, additional irregular very low frequency (& 0.05±0.2 Hz) oscillations were observed in u(ti) consistent with possible cyclic changes in cerebral blood flow or changes of vigilance and muscle tone. These preliminary results suggest that the new method can detect subtle changes in the overall pattern of the EEG without the necessity of making tenuous assumptions about stationarity. KEYWORDS

EEG, state change, Independent Component Analysis, Sleep.

INTRODUCTION The scalp electroencephalogram (EEG) represents the summed synaptic activity from large groups of neurons primarily from the cerebral cortex (Hubbard et al. 1969), making it well suited for the non-invasive monitoring of widespread oscillations observed during sleep. Many EEG features recorded during sleep tend to be distributed widely over the scalp electrodes (Amzica and Steriade 1997). The classification of sleep stages is largely based on the dominant rhythms of the widespread EEG activity, without taking into consideration the relationships between different derivations. In the human scalp EEG of slow-wave sleep (SWS), oscillations are present in the frequency below 1 Hz, in the delta band (1±4 Hz) and in the spindle frequency band (11±15 Hz) (Achermann and BorbeÂly 1997). Correspondence: Dr M.J. McKeown, Computational Neurobiology Laboratory, Salk Institute for Biological Studies, 10010 North Torrey Pines Road, La Jolla, CA 92037±1099, USA. Tel.: 619 453 4100 (ext. 1561); fax: 619 5870417; e-mail: [email protected]

48

Sleep spindles, which consist of waxing and waning synchronized oscillations and are a hallmark of stage 2, reflect the activity in a reciprocal thalamo-cortical network that has been extensively investigated in animals in vivo (Contreras et al. 1997; Steriade 1997) and modeled computationally (Destexhe et al. 1994). Cells in the reticular nucleus of the thalamus play an important role in the synchronization of the cortical network, and demonstrate spindling properties even after removal of the cortex (Steriade et al. 1987). The waxing part of spindling activity observed in the scalp EEG is the result of progressive recruitment of neurons into a synchronized network oscillation, while the waning part may be due to the activation of a non-inactivating cation current (Ih) in thalamocortical cells (Bal and McCormick 1996). Oscillations in the delta band, which dominate in stages 3 and 4, correspond to widespread cellular activity generated by at least two mechanisms. The cortex itself is capable of rhythmic oscillations in the delta range, even after removal of the thalamus (Steriade et al. 1993a). Thalamocortical cells can #1998 European Sleep Research Society

Ahed Bhed Ched Dhed Ref marker Fig marker Table marker Ref end Ref start

Paper 30 Disc

A new method for detecting state changes in the EEG also produce oscillations in this range by an interaction of voltage-gated currents, It and Ih (Pape and McCormick 1995). Thalamocortical interactions organize globally the spatiotemporal coherence of delta oscillations (Contreras et al. 1996). Recently, Achermann and BorbeÂly (1997) have demonstrated a low-frequency component in the human sleep EEG. In anesthetized and sleeping cats, low-frequency rhythmic depolarizations followed by prolonged hyperpolarizations have been demonstrated in widespread areas of the cortex (Steriade et al. 1993b). A close temporal relationship exists between low-frequency components of the cat sleep EEG and components in the cortical and thalamic cells, with the EEG oscillations corresponding to synchronized oscillatory activity in cortical networks. The hyperpolarizations are associated with prolonged depth positive waves in the cortical EEG and the slow oscillations may group the thalamically generated spindles and the delta waves in sequences that recur every 1.1± 3.3 s (Steriade et al. 1993b, Contreras et al. 1995). The low frequencies present in the sleep EEG have not been fully appreciated until relatively recently. The high pass filters typically used in standard EEG recordings may have prevented the precise specification of these low frequency oscillations. Achermann and BorbeÂly (1997) circumvented this problem by compensating the attenuation of the EEG signal by the high pass filters and then measuring the EEG power spectrum. However, as Amzica et al. (1997) point out, higher frequency waveforms recurring at a low frequency are relatively robust to filtering at low frequencies. Nevertheless, slow (5 1 Hz) dynamic changes in the EEG may be difficult to characterize because spectral techniques, such as the fast Fourier transform (FFT), require the assumption of stationarity of the EEG signal (Press et al. 1992), a prerequisite that may not apply to the large window size required to capture these low frequencies (Geering et al. 1993). Since the EEG reflects the ever-changing cortical activity which varies spatially as well as temporally, a complementary way to monitor the fluctuating widespread patterns in EEG activity is to measure the relationships between different derivations rather than to monitor the frequencies recorded from a single derivation. Measuring the changes from interactions between derivations would be relatively robust to the filtering of any specific frequency, and would preclude the problem of varying spectral characteristics in a data set. Here we describe a sensitive technique for detecting changes in the interactions in the EEG by employing an indirect method: calculating the (±log) likelihood of observing the EEG data under the (probably incorrect) assumption that the interactions between EEG electrodes are constant over a given period of time or epoch. We apply this method to demonstrate and characterize low-frequency interactions present in the human sleep EEG. A PROBABILISTIC INTERPRETATION TO THE EEG The digitized EEG data for a given segment of time, or epoch, can be stored in an n by p matrix X, where n is the number of

49

channels and p is the number of time points in the epoch. The time waveforms recorded from each of the EEG electrodes during the epoch can be transformed into a new set of waveforms, C, by taking linear combinations, defined by an n by n matrix Z, of the recordings from each channel. Thus, C = ZX

(1)

where C is an n by p matrix of transformed waveforms, Z is an n by n matrix containing combinations of channels, and X is the n by p data matrix. If Z is an invertable matrix, we can define A = Z71, so that: X=AC

(2)

Equation (2) implies that the recorded data can be accurately modeled as component waveforms, C, linearly combined as specified in the matrix A. Using Bayes' theorem, we can determine the probability of observing the instantaneous voltages across all scalp electrodes at a given point of time, P(Xi), under the model specified by equation (2), i.e. P(Xi|A) (see Technical Appendix). By estimating -log(P(Xi|A)) for each time point i in a given epoch, we define a sequence of scalar values through time, which we define as u(ti). As this is an estimate of the log-likelihood of observing the data under the model specified in equation (2), the u(ti) function is a dimensionless sequence of numbers, with the ti's ranging over the number of time points in the epoch. In order to detect the second-to-second changes in the EEG, the u(ti) function can be smoothed using a running average technique to create a smoothed function, us(ti). As discussed in the Technical Appendix, calculation of u(ti) is greatly facilitated by judiciously selecting the combinations of channels represented by Z in (1) so that the component waveforms contained in the rows of C are maximally statistically independent. This is possible with the independent component analysis (ICA) algorithm of Bell and Sejnowski (1995, 1997). The general ICA algorithm assumes that the data consists of n independent, linearly mixed underlying components, and tries to blindly ùnmix' the data to recover the original components. When applied to EEG data (Makeig et al. 1996), ICA determines the optimum combination of channels via an iterative training process that separates the EEG into component waveforms that are as statistically independent as possible so that the component waveforms sum together to the original EEG (Fig. 1 A). Separation of the sleep EEG in this manner allows the specification of the time course and spatial distribution of components which reflect many features in the sleep EEG, including spindles, K-complexes and electrocardiographic (ECG) artifacts (Fig. 1 A). The ICA algorithm is related to principal component analysis (PCA) techniques that have been applied to EEG (Lagerlund et al. 1997; Skrandies 1993). PCA requires that the scalp distributions are uncorrelated, a requirement which is unnecessary with ICA. Also, PCA separates the EEG into waveforms that are merely orthogonal, a weaker criterion for statistical independence than the one used by ICA. These differences can be attributed to the fact that the goals of PCA and ICA differ. PCA tries to summarize the variability seen in

#1998 European Sleep Research Society, J. Sleep Res., 7 (Suppl. 1), 48±56.

Paper 30 Disc

50

M. McKeown et al.

Figure 1. The ICA algorithm applied to the EEG. (A) The ICA algorithm assumes that the multichannel EEG recording (top half of figure) is composed of a linear mixture of underlying components (lower half of figure). The goal of the algorithm is to blindly ùnmix' the data to determine the underlying components. Each component consists of a waveform, which is as independent as possible from all other waveforms, and a specific scalp distribution. The scalp distribution, shown here as a gray-scale, specifies the relative contribution of the associated waveform to the raw EEG at each channel location. Components that can be separated out from the sleep EEG may be related spindles (second row, left), K-complexes (second row, right), EKG artifacts (bottom row, left) and slow rhythms (bottom row, right). The components shown are illustrative and did not all come from the same epoch. Note that the ICA algorithm assumes that the scalp patterns remain constant throughout the epoch. (B) Procedure used for kernel density estimation. Once the component waveforms were derived using the ICA algorithm (A), a histogram of values was determined for each waveform. This histogram was then convolved (smoothed) with a gaussian to provide an estimate of the underlying probability distribution. Under the independence assumption, the probability of observing all voltages at a given point will be the product of the individual probabilities from each of the component waveforms.

all channels of the EEG in as few components as possible. In contrast, ICA assumes that the EEG is composed of mixed, underlying components and tries to determine how to ùnmix' the data to recover the underlying components. The ICA algorithm determines the average relationships between channels during all the time points in an epoch, and hence assumes that the particular combinations of channels required to separate the data into independent component waveforms remain constant throughout the epoch. If there is a violation of any of the assumptions of ICA, i.e. there is not exactly the same number of independent components as scalp electrodes, the mixing of the components changes during the epoch, or there is non-linear mixing of the components, then the ICA algorithm will be less able to separate components which are truly statistically independent. Models of the data derived under the assumption that the data can be separated into independent components will therefore appear less likely. Thus, the derived u(ti) will increase (since it estimates log(P(Xi|A))) when parts of an EEG epoch are not an accurate reflection of any of these assumptions.

NREM-REM sleep cycle. The subject was free of sleep disorders and was instructed to abstain from caffeine, alcohol or any other pharmacologically active substances. The data were collected with a portable 32-channel system (Neurofile II, Nihon Kohden, Japan). Twenty-nine EEG channels (including mastoids) with electrode positions according to an extended 10±20 system, an EMG and a bipolar EOG were recorded. Data were sampled at 256 Hz, digitally filtered (low-pass FIR filter, 73 dB at 40 Hz) and stored with a resolution of 128 Hz. Twenty-second artifact-free epochs were separated with the ICA algorithm of Bell and Sejnowski (1995) and sleep stages were scored according to conventional criteria. To detect very slow oscillations in stages 3 and 4, 50-s epochs were used. Kernel density estimation (Fig. 1B) was performed on each row of C in equation 2, Ci, by first calculating a histogram, using 100 equally spaced bins, of scalp voltages from that channel. This histogram was then convolved (smoothed) with a Gaussian function (Baxter et al. 1997) whose standard deviation was, 1.0592* n* s71/5

EXPERIMENTAL METHODS The data were recorded from one healthy male from 23.00 to 07.00 hours on one night. Analysis was restricted to the first

where n was the number of time points in the epoch, and s the standard deviation of the values recorded from the ith channel during the epoch, i.e. Ci (Fig. 1B).


Paper 30 Disc

A new method for detecting state changes in the EEG The u(ti) function was calculated for each epoch. The smoothed u(ti) function, us(ti), was determined by calculating a running average using a 250 point Hanning window (Press et al. 1992) centered on the time point of interest. This window was shifted one time point (1/128 s) at a time to calculate us(ti) at each time point ti. To determine whether the u(ti) function was related to specific frequency bands contained in the raw EEG, we correlated the u(ti) with bandpass-filtered portions of the EEG. For all stage 2 and 4 epochs, the EEG was bandpass-filtered at delta (0±4), theta (4±8), alpha (8±13) and beta (13±64 Hz) frequencies. The u(ti) was calculated for each 20-s stage 2 and 4 epochs as described above. The u(ti)'s from each epoch were then concatenated together to form a long sequence which was then correlated separately with each channel in the corresponding band-pass filtered EEG. To ascertain if the u(ti) function tracked the relative amounts or modulations of the different frequency rhythms contained in the EEG, we correlated the u(ti) with an estimate of the instantaneous power contained in each of the four frequency bands (delta, theta, alpha and beta) for each channel. The EEG was divided into 1-s, Hanning windowed (Press et al. 1992) sections, which were then zero padded to form overlapping, 2-s segments. A discrete Fourier transform (n = 256) was used to estimate the spectral power density at the center of the window, which was then shifted 0.5 s to obtain the next estimate, and so on. This procedure provided estimates for the instantaneous power of each frequency band every 0.5 s For comparisons of the instantaneous power estimates to the u(ti), the u(ti) was calculated for all 20-s stage 2 and 4 epochs separately and then concatenated together. This long sequence was resampled using a Chebyshev type I lowpass filter with a cutoff frequency of 0.8 Hz to also give values every 0.5 s (64 points) which was then correlated with the instantaneous power estimate of each frequency band at each channel. To monitor how the u(ti) function varied across sleep stages, the mean of the u(ti) function for epoch, Z(Tj), was calculated for each epoch: Tj

n 1X uti n i1

51

for example, yet sensitive to abrupt overall changes in the EEG pattern, such as spindling at 6, 17 and 19 s (Fig. 2 A). Increase of the u(ti) in response to spindling was a consistent feature in all stage 2 epochs. In order to verify that broad peaks in the us(ti) were secondary to changes in the overall relationships between the channels of the EEG and not random fluctuations, we differentially permutated the time points from each channel of the EEG for a given epoch to create a new data set, and then re-calculated the us(ti). The slow fluctuations were no longer apparent after the shuffling procedure (Fig. 2C). Recomputing the u(ti) and us(ti) using only 9 representative channels (F3, F4, C3, C4, P3, P4, O1, O2, and Cz) for this epoch instead of the full 29 channels was qualitatively similar (data not shown), with a correlation of r = 0.63 for u(ti) and r = 0.89 for us(ti). We performed a number of simulations to investigate the characteristics of us(ti). To assess the effect of altered withinepoch linear mixing of the presumed underlying sources of variability in the data, a simulated 29-channel data set was created by taking 29 linear combinations (i.e. linearly mixing) of 29 20-s random time series. Each random time series was

3

where n is the number of time points in the epoch and Tj represents the epoch for which Z() is calculated.

RESULTS Calculation of the channel combinations to optimally separate the EEG into independent waveforms by the ICA algorithm (Z in equation 1) for a 20-s epoch took & 80 s on a DEC Alpha computational server. Subsequent calculation of the u(ti) for the entire epoch took an additional 7 s. In stage 2, after the ICA algorithm was trained on the entire duration of the epoch, the us(ti) function (Fig. 2B) appeared to be relatively insensitive to large amplitude changes from eye movements,

Figure 2. The u(ti) function for a stage 2 sleep epoch. (A) A stage 2 epoch includes intermittent sleep spindles. (B) The u(ti) (thin line) and us(ti) (thick line) derived from this epoch showed irregular oscillations, with peaks corresponding to the spindling changes. (C) Shuffling the order of the time points from each electrode and recalculating the us(ti) based on the shuffled data resulted in removal of the oscillations (thick line).


Paper 30 Disc

52

M. McKeown et al.

sampled from a kurtotic distribution (i.e. having more values near zero and longer tails than a Gaussian distribution) at 128 Hz, the same frequency as for the EEG. For the first 10% of the epoch, a random set of linear combinations of the time series was chosen, and then the mixing was changed suddenly to another random set of combinations for the remaining 90% of the epoch (Fig. 3). Note the sudden change in us(t) at the time the mixing which shows that the interactions between channels abruptly changed (Fig. 3). To determine which increases in us(ti) during spindling were due to within-epoch changes in the linear mixing of the assumed underlying sources in the EEG and which were due to violations of the other assumptions of ICA such as linearity, we calculated a new us(ti) based on a subset of time points. We first selected time points from an epoch which displayed the highest us(ti) (Fig. 4 A, hatched areas), and then calculated new combinations of channels with the ICA algorithm based on this reduced data set. The new combinations were then applied over the whole epoch to create a new us(ti). Many peaks in the new us(ti) corresponded to the original us(ti), especially at t = 6, 13, 17, and 19 s (Fig. 4B). In stage 3, the us(ti) demonstrated irregular low-frequency oscillations of & 0.6 Hz (Fig. 5B). K-complexes occurred near the peaks in the us(ti) (Fig. 5B). As progressively more delta appeared signaling the onset of stage 4, the slow (& 0.6 Hz) oscillations became more regular, and were superimposed on much slower irregular oscillations of & 0.05 Hz (Fig. 6B). Other sections of the EEG in stage 4 demonstrated similar irregular low-frequency oscillations of & 0.05±0.2 Hz (data not shown). To test if the us(ti) was merely following the amplitude changes in the EEG, we correlated the us(ti) with individual

Figure 3. Effects on u(ti) of abruptly changing the linear mixing of synthetic sources during an epoch. A 29 channel synthesized data set was generated by randomly mixing 29 random times series. For the first 10% of the epoch, the random time series were mixed with using 29 fixed, randomly chosen linear combinations. The combinations were abruptly changed to a different mixing for the remaining 90% of the epoch. Note the abrupt reduction in the u(ti) corresponding to the time where the mixing in the synthesized data set changed.

Figure 4. Possible mechanisms for increases in the us (ti). (A) The u(ti) (thin line) and us(ti) (thick line) from a stage 2 epoch demonstrated time periods where the us(ti) was greater (hatched areas), and hence the data were less likely to be observed under an assumed ICA model. A reduced EEG data-set incorporating these time periods in the epoch was used by ICA to calculate a new set of electrode combinations. (B) The new electrode combinations were then applied to the whole epoch to create a new us(ti). Note the similar peaks observed in both the original us(ti) (dotted line) and the new us(ti) (solid line) at t = 6, 13, 17, and 19 s.

channels from band-pass filtered EEG (Fig. 7 A). There was minimal correlation between u(ti) and the individual channels in both stage 2 and stage 4 for all frequency bands. As the us(ti) appeared to follow the modulation of the different frequency rhythms in EEG in stage 2±4, we estimated the instantaneous power in each of the frequency bands in the EEG by spectrogram techniques. Over all stage 2 epochs, the resampled u(ti) correlated best (r = 0.47±0.48) with the instantaneous power estimates of the EEG in the theta and alpha bands. In stage 4 epochs, the u(ti) correlated better with instantaneous power estimates in the delta and theta bands (Fig. 7B). The mean of the u(ti) for a whole epoch, Z(Tj), appeared to be consistent within sleep stages and sensitively detected transitions between sleep stages (Fig. 8). Stage 4 sleep appeared to have consistently higher values for Z(Tj), compared with stage 2 or REM sleep. The highest values for Z(Tj) were present when movement or other artifacts contaminated the data as these represented sudden changes in the overall EEG pattern.


Paper 30 Disc

A new method for detecting state changes in the EEG

Figure 5. The us(ti) function for a stage 3 epoch. (A) The EEG displays a prominent K-complex in the center of the epoch. (B) The us(ti) for this epoch demonstrated slow, non periodic oscillations roughly every 2 s, with the highest peak in the us(ti) corresponded to the K-complex. (C) A power spectrum of the us(ti) based on the Fast Fourier Transform (Hanning Window of 512 points, zero-padded to 1024 points, with overlap of 448 points). The lower frequencies are cut off to more clearly demonstrate the peak at about 0.6 Hz.

DISCUSSION The standard tool for quantitative EEG analysis is the FFT. The FFT makes the assumption that the EEG is stationary, i.e. the spectral content of the EEG is constant during a prespecified window (Press et al. 1992). In order to capture very low frequencies in the EEG, large window sizes are required, making stationarity a tenuous assumption (Geering et al. 1993). The log-likelihood of observing the EEG under the assumption that the relationships between EEG channels were stationary or constant, u(ti), sensitively detected changes in the overall structure of the EEG as it appeared across all scalp electrodes without the need to make stationarity assumptions. The irregularity of the slow oscillations detected in stages 2, 3 and 4 (Figs 2, 5, 6) suggests that the assumption of stationarity may be problematic for other methods estimating low frequency oscillations in the sleep EEG. The u(ti) appeared to increase when there were changes in the global pattern of the EEG. Thus, the u(ti) correlated poorly with the different frequency components of the raw EEG (Fig. 7A), yet correlated relatively well with the modulation of the

53

Figure 6. The us(t) function for a stage 4 epoch. (A) The EEG is dominated by delta activity in this stage 4 epoch, of which 5 representative channels are shown. (B) The calculated us(ti) demonstrated regular oscillations & 0.6 Hz, superimposed on much slower irregular oscillations of periods & 20 s (C) A power spectrum of the us(ti) based on the Fast Fourier Transform (1024 points, Hanning Window of 512 points, overlap of 448 points), again shows a peak at about 0.6 Hz.

different frequency bands over time (Fig. 7B). The similarity between the correlation of u(ti) with the instantaneous estimate of delta power and the delta band-passed filtered EEG in stage 4 (Fig. 7A-B) might reflect the different frequencies of oscillations present in the us(ti). Very low-frequency oscillations (& 0.05 Hz) observed in the us(ti) corresponded to the increased content of delta in the EEG noticeable by inspection (Fig. 6), consistent with the computed correlation between the us(ti) and the estimate of instantaneous delta power (Fig. 7B). Oscillations in the EEG at slightly higher frequencies (e.g. in the delta range) may still be considered a state change for the size of smoothing window (n = 250) used to calculate the running average, explaining the correlation between the us(ti) and the delta band-passed EEG (Fig. 7A). When spindles appeared within an epoch in stage 2 sleep, there were abrupt increases in the us(ti) (Fig. 2), suggesting that time points during the spindle were less likely to be observed under a model where the data were considered as the sum of linearly mixed, temporally independent components. Spindles occur when the cortex is functionally isolated from afferent signals normally passing through the thalamus and oscillations are generated in the reciprocal thalamocortical network


Paper 30 Disc

54

M. McKeown et al.

Figure 7. Comparison of u(ti) with bandpass-filtered EEG and spectrogram estimates. (A) Correlation of u(ti) with channels from bandpass-filtered EEG. The amounts shown are the mean of the correlation between each channel and the concatenation of the u(ti)'s calculated for all 20-s stage 2 (dark bars) and stage 4 epochs (light bars), with the error bars signifying one standard deviation. Correlation values were first transformed with a Fisher's z-transform. Note the relatively poor correlation over both stages and all frequency bands (compare with (B)). (B) The correspondence of u(ti) with estimates of instantaneous power contained in different frequency bands across stages 2 and 4. The values show the mean of the correlation of the instantaneous power estimate for all channels with the resampled u(ti), with error bars specifying one standard deviation. Correlation values were first transformed with a Fisher's z-transform. Note the much higher correlations compared to (A).

Figure 8. The variability of Z(Tj) across sleep stages for the first nonREM-REM sleep cycle divided into 20 s epochs (1± 4: stages 1±4 of nonREM sleep, r: REM sleep; m: movement time; w: waking, a: artifact). The Z(Tj) was consistent within sleep stages. The highest values of Z(Tj) corresponded to artifact. Stage 4 sleep corresponded to consistently higher values of Z(Tj) compared to stages 2±3.

(Steriade 1995). Interactions between the EEG channels are therefore likely to change at the onset and completion of the spindles, violating the assumption of within-epoch constant electrode relationships underlying the ICA. The shuffling of the data removed the oscillations present in the us(ti) (Fig. 2C), suggesting that there is a distinct temporal structure to the periods within an epoch where the ICA assumptions are violated. Changes in linear mixing of presumed underlying components during an epoch were shown via simulation to be one potential mechanism for increases in the u(ti) (Fig. 3). In this simulation (Fig. 3), the ICA algorithm, in an attempt to find the average channel combinations for unmixing over the entire epoch, was biased towards determining the unmixing suitable for 90% of the epoch. The different mixing used for the first 10% of the epoch was therefore less likely to be observed under the model which assumes that the mixing is constant over the entire epoch, reflected by the higher u(ti). Another simulation where the two different mixings were split 50%±50% (not shown) revealed no abrupt changes in the u(ti), as the algorithm determined the average mixing over the entire epoch, making the ICA-derived model equally improbable during both halves of the epoch. In order to determine the relative importance of changes in linear mixing within an epoch for increases in the u(ti), we calculated new channel combinations for unmixing based on periods during the epoch that the linear model fit less well (Fig. 4B). The re-calculated us(ti), based on the new calculated electrode relationships, demonstrated increases similar to the original us(ti) during parts of the epoch (Fig. 4B). This suggests that these peaks in the original us(ti) occurred not because the linear mixing of the presumed sources of signal changed, but due to a possible non linear mixing of the presumed underlying sources, and/or because the number of sources contributing to


Paper 30 Disc

A new method for detecting state changes in the EEG the scalp EEG changed during the epoch. The changing interactions between electrodes in the sleep EEG may therefore be the result of several mechanisms. These results may have implications for techniques used for EEG artifact removal (Jung et al. 1998; Lagerlund et al. 1997) and evoked response potential (ERP) analysis (Makeig et al. 1997) which make the assumption that the relationship between channels remains constant throughout an epoch. In stages 3 and 4, the & 0.6 Hz oscillations in the us(ti) (Figs 5±6) are consistent with low-frequency oscillations detected in sleeping and anesthetized cats (Steriade et al. 1993a; Steriade et al. 1993b) and human subjects (Achermann and BorbeÂly 1997). Under urethane anesthesia, feline cortical cells in layers II-VI show spontaneous low-frequency oscillations (Steriade et al. 1993a; Steriade et al. 1993b). As the oscillation spreads cortically, the interactions between the EEG channels would be expected to differ, increasing the us(ti). The EEG correlate emerging from these slow oscillations recently has been shown to be the K-complexes (Amzica et al. 1997). During the early part of the K-complex, the slow oscillations are widely spread throughout the cortex (Amzica et al. 1997) via intracortical linkages. The peaks in the us(ti) occurring at & 1.5±2-s intervals may therefore reflect subtle K-complexes not obviously discernible in the raw EEG. The very slow oscillations of & 0.05±0.2 Hz observed in the us(ti) in stage 4 sleep (Fig. 6) have been recently documented for human sleep (Achermann and BorbeÂly 1997). They probably reflect a different type of state changes in the EEG. The cyclic alternating pattern, a correlate of sleep associated with changes in vigilance and muscle tone, typically has periods of & 20±40 s (Terzano and Parrino 1993; Terzano et al. 1988). Slow cyclical variations in cerebral blood flow velocity have been observed during sleep in normal newborns with a frequency of between 2 and 6 cycles/min (Ferrarri et al. 1994). In spinalized rats, Golanov and Reis (1995) explored the relationship between changes in regional cerebral blood flow and electrocortigram (ECoG) activity. Spontaneous cerebrovascular waves with a frequency of 6/min were always preceded by high-amplitude bursts of ECoG activity. Any of these very low frequency mechanisms could affect the overall pattern of the EEG, and thus could be detected by the us(ti). The EEG in stages 3 and 4 is increasingly dominated by a single cortical rhythm, making the assumption that there are as many independent temporal components as there are channels, and that their contributions sum linearly, less valid. This probably explains why the us(ti), while still reflecting withinepoch state changes, tended to increase in stages 3 and 4 when averaged over an entire epoch (Z(Tj)) (Fig. 8).

55

oscillations detectable in slow wave sleep. However, to demonstrate the validity of this new method, it needs to be tested by a systematic application to larger data sets from several subjects and by a statistical analysis.

TECHNICAL APPENDIX The digitized EEG data, X, for an epoch can be modeled as: X=AC

(2)

where X is an n by p matrix, n is the number of channels and p is the number of time points in the epoch, C is an n by p matrix of component waveforms, and A is an n by n matrix containing linear combinations of component waveforms. Determining the likelihood of observing the all EEG channel voltages at the ith time point, Xi, under the model specified in equation 2 by A and C gives, 1

PXi jA

PXi jA; CPCdC

4

ÿ1

This can be simplified by making a change of variables (Olshausen and Field 1997). From (2): X=AC then, C Aÿ1 X; dC

1

PXi jA

dx jdetAj

5

Xi ÿ XPAÿ1 X

ÿ1

dX jdetAj

PAÿ1 Xi PCi jdetAj jdetAj

6

where Ci is the ith column of the matrix in (1±2). If it is possible to choose Z (= A71) in (1) such that C is separated into rows that are approximately statistically independent, then computation of (6) can be estimated by: n Q

PXi jA

PkCki PCi k1 jdetAj jdetAj

7

Manipulation of (7) is simplified by taking the logarithm of both sides of the equation and multiplying by 71: ÿlogPXi jA logdetA ÿ

n X

logPk Ck i uti

8

k1

CONCLUSION We have introduced a new method for monitoring overall changes in the EEG, by providing a statistical interpretation to the ongoing interactions between derivations. Our exploratory results suggest that this method can detect abrupt changes in the sleep EEG, such as spindles, and also the lower frequency

The matrices A and C were determined by applying ICA to the data matrix X as described in (Bell and Sejnowski 1995, 1997). Performing kernel density estimation (Baxter et al. 1997), or more crudely, taking a smoothed histogram of each row of C (over the whole epoch) provides an estimate of the probability distribution function for each of the Pk(Ck), from


Paper 30 Disc

56

M. McKeown et al.

which the probability of the ith point in the row, Pk(Cki), can be estimated (Fig. 1B). ACKNOWLEDGEMENTS Dr McKeown is supported by the Heart and Stroke Foundation of Ontario, Canada. Drs Achermann and BorbeÂly acknowledge the support by the Swiss National Science Foundation. The study was supported by the Human Frontiers Science Program project RG-81/96. The authors are grateful for many useful discussions with Drs Mike Lewicki, Bruno Olshausen, Te-Won Lee, Scott Makeig and Tzyy-Ping Jung. We also thank Mr Michiel van Burik from the Institute of Biomedical Engineering, Graz University of Technology, for providing us with the mesh data required to create the head images in Figure 1. REFERENCES Achermann, P. and Borbely, A. Low-frequency (51Hz) oscillations in the human sleep EEG. Neuroscience, 1997, 81: 213±222. Amzica, F. and Steriade, M. The K-complex: Its slow (51Hz) rhythmicity and relation to delta waves. Neurology, 1997, 49: 952±959. Bal, T. and McCormick, D. A. What stops synchronized thalamocortical oscillations? Neuron, 1996, 17: 297±308. Baxter, M., Beardah, C. and Wright, R. Some archaeological examples of kernel density estimates. J. Archaeol Sci., 1997, 24: 347±354. Bell, A. J. and Sejnowski, T. J. An information-maximization approach to blind separation and blind deconvolution. Neural Comput., 1995, 7: 1129±11 59. Bell, A. J. and Sejnowski, T. J. The Ìndependent Components' of natural scenes are edge filters. Vision Res., 1997, 37: 3327±3338. Contreras, D., Destexhe, A., Sejnowski, T. J. and Steriade, M. Control. of spatiotemporal coherence of a thalamic oscillation by corticothalamic feedback. Science, 1996, 274: 771±774. Contreras, D., Destexhe, A., Sejnowski, T. J. and Steriade, M. Spatiotemporal patterns of spindle oscillations in cortex and thalamus. J. Neurosci., 1997, 17: 1179±1196. Contreras, D. and Steriade, M. Cellular basis of EEG slow rhythms: a study of dynamic corticothalamic relationships. J. Neurosci., 1995, 15: 604±622. Destexhe, A., Contreras, D., Sejnowski, T. J. and Steriade, M. A model of spindle rhythmicity in the isolated thalamic reticular nucleus. J. Neurophysiol., 1994, 72: 803±818. Ferrarri, F., Kelsall, A. W., Rennie, J. M. and Evans, D. H. The relationship between cerebral blood flow velocity fluctuations and sleep state in normal newborns. Pediatr. Res., 1994, 35: 50±54.

Geering, B., Achermann, P., Eggimann, F. and Borbely, A. Periodamplitude analysis and power spectral analysis: a comparison based on all-night sleep EEG recordings. J. Sleep Res., 1993, 2: 121±129. Golanov, E. V. and Reis, D. J. Vasodilation evoked from medulla and cerebellum is coupled to bursts of cortical EEG activity in rats. Am. J. Physiol., 1995, 268: R454±R4 67. Hubbard, J., Llinas, R. and Quastel, M. Electrophysiological Analysis of Synaptic Transmission. Arnold, London, 1969. Jung, T.-P., Humphries, C., Lee, T.-W., Makeig, S., McKeown, M., Iragui, V. and Sejnowski, T. Extended ICA removes artifacts from electroencephalographic recordings. In: D. Touretzky, M. Mozer and M. Hasselmo (eds). Advances in Neural Information Processing Systems, 10. MIT Press, Cambridge, MA, 1998. Lagerlund, T. D., Sharbrough, F. W. and Busacker, N. E. Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition. J. Clin. Neurophysiol., 1997, 14: 73±82. Makeig, S., Bell, A. J., Jung, T.-P. and Sejnowski, T. J. Independent component analysis of electroencephalographic data. Adv. Neural Information Processing Systems, 1996, 8: 145±151. Makeig, S., Jung, T.-P., Bell, A., Ghahremani, D. and Sejnowski, T. Blind separation of auditory event-related brain responses into independent components. Proc. Nat Acad. Sci. (USA) 1997, 94: 10979±10984. Olshausen, B. A. and Field, D. J. Sparse coding with an overcomplete basis set: a strategy employed by V1? Vision Res., 1997, 37: 3311±3325. Pape, H. C. and McCormick, D. A. Electrophysiological and pharmacological properties of interneurons in the cat dorsal lateral geniculate nucleus. Neuroscience, 1995, 68: 1105±1125. Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. Numerical Recipes in C. The Art of Scientific Computing. Cambridge University Press, Cambridge, 1992 (second edition). Skrandies, W. EEG EP: new techniques. Brain Topogr., 1993, 5: 347±350. Steriade, M. Thalamic origin of sleep spindles: Morison and Bassett (1945). J. Neurophysiol., 1995, 73: 921±922. Steriade, M. Synchronized activities of coupled oscillators in the cerebral cortex and thalamus at different levels of vigilance. Cereb Cortex, 1997, 7: 583±604. Steriade, M., Domich, L., Oakson, G. and Deschenes, M. The deafferented reticular thalamic nucleus generates spindle rhythmicity. J. Neurophysiol., 1987, 57: 260±273. Steriade, M., Nunez, A. and Amzica, F. A novel slow (5 1 Hz) oscillation of neocortical neurons in vivo: depolarizing and hyperpolarizing components. J. Neurosci, 1993a, 13: 3252±3265. Steriade, M., Nunez, A. and Amzica, F. Intracellular analysis of relations between the slow (5 1 Hz) neocortical oscillation and other sleep rhythms of the electroencephalogram. J. Neurosci., 1993b, 13: 3266±3283. Terzano, M. G. and Parrino, L. Clinical applications of cyclic alternating pattern. Physiol. Behav., 1993, 54: 807±813. Terzano, M. G., Parrino, L. and Spaggiari, M. C. The cyclic alternating pattern sequences in the dynamic organization of sleep. Electroencephalogr. Clin. Neurophysiol., 1988, 69: 437±447.
