952
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39. NO. 9, SEPTEMBER 1992
On the Tracking of Rapid Dynamic Changes in Seizure EEG Isak Gath, Claude Feuerstein, Dinh Tuan Pham, and Gerard Rondouin
Abstract-Estimation of autospectra and coherence and phase spectra of seizure EEG, using the FFT technique, will cause “smearing” of the rapid dynamic changes which occur during the seizure. This is inherent to FFT spectral estimation, due to the averaging process which is necessary in order to get consistent spectral estimates. A different approach suggested in the present study is to carry out multivariate autoregressive modeling of the multichannel seizure EEG, combined with adaptive segmentation. In order to obtain good estimates in cases of short record length, the vectorial AR modeling was based on residual energy ratios. The method has been tested on multichannel seizure EEG recordings from rats with focal epilepsy, caused by intracerebral administration of Kainic acid, and in depth EEG recordings in patients with temporal lobe epilepsy.
I. INTRODUCTION
P
OWER spectrum estimation is a useful tool in the quantitative analysis of EEG recorded during epileptic seizures. Autospectra, coherence and phase spectra, might be helpful in the interpretation of the multichannel dynamic changes which occur prior and during the seizures. The ictal EEG have been most often analyzed using the fast fourier transform (FFT) [1]-151. However, the “epileptic” EEG signal recorded either in experimental animals or in humans might show marked nonstationarity, and hence, estimation of the power spectrum using the FFT is not without problems. Periodograms have to be computed for short stationary signal segments, usually not exceeding 2-3 s, which means frequency resolution of 0.33-0.5 Hz. In order to obtain consistent power spectrum estimates, averaging of several periodograms is required to decrease the variability of the estimate. This averaging process inherent in the FFT technique will cause “smearing” of the dynamic changes which might occur during the seizure. An alternative approach to the estimation of the power spectrum and phase and coherence spectra is based on parametric modeling of the multichannel EEG. A time varying multivariate autoregressive (AR) modeling, based Manuscript received September 12, 1990; revised December 5 , 1991. I. Gath is with the Department of Biomedical Engineering, Technion, Israel Institute of Technology, 32000 Haifa, Israel. C. Feuerstein is with the Laboratoire de Physiologie Section Neurophysiologie, INSERM U318 CHU, Grenoble. France. D. T. Pham is with the Laboratoire de Modelisation et CalculilMAG, CNRS, and Universite Joseph Fourier, Grenoble, France. G. Rondouin is with the Laboratoire de Medecine Experimentale, INSERM U 239, University of Montpellier, Montpellier. France. IEEE Log Number 9201968.
on the deterministic-regression method using a family of basis functions have been suggested by Gersch 161. This method overcomes the problem of nonstationarity but suffers from a heavy computational load. The number of parameters fitted to the model increases as the square of the number of EEG channels [6], and hence in practice this method is useful for modeling of only a small number of EEG channels. A still different approach will be to carry out multichannel AR modeling in combination with adaptive segmentation of the EEG [7]-[9]. However, whatever the procedure, in order to trace the rapid dynamic changes inherent in the epileptic EEG, the AR model estimator is required to be efficient in case of short record lengths (2 s or less). The aim of the present paper is to investigate a procedure for parametric multivariate modeling of the epileptic EEG, a procedure aimed to detect rapid dynamic changes in the recorded EEG signal. 11. METHODS A . Data Acquisition Generalized seizures with focal onset were generated in two male Wister rats by intracerebral injection of Kainic acid. A cannula was inserted stereotaxically under anesthesia into the right basolateral amygdala, and at the same time bipolar twisted electrodes were implanted in both amygdala for continuous depth recording. Surface EEG recordings using three bipolar channels (left frontal and left and right occipital leads) were obtained by implanted gold plated screws over the respective cerebral cortex, and the electrode assembly was secured to the skull with a dental acrylic cement. Recordings were started ten days after the implantation, acute seizures were obtained by a microinjection of 0.2 pg Kainic acid into the right amygdala. Data acquisition was carried out on an IBM PClAT computer with a Data Translation A/D card, using a sampling rate of 128 Hz. In addition, recordings were carried out in several patients having temporal lobe epilepsy. Depth electrodes were placed using the stereotactic technique in both mesial and cortical structures, in particular, around the limbic system on both sides. Long-time monitoring of these patients was carried out lasting 2-14 days, until sufficient ictal data was accumulated. Typical generalized seizures having a focal onset were chosen for the
0018-9294/92$03.00
1992 IEEE
953
GATH er a l . : TRACKING OF RAPID DYNAMIC CHANGES IN SEIZURE EEG
analysis, and data acquisition was carried out as described for the rats treated with Kainic acid.
B. Data Analysis Signal segments including both normal EEG and ictal activity were chosen for analysis. Two different strategies were applied to the nonstationary seizure EEG signal. In order to obtain quasi-stationary signal segments from the nonstationary seizure EEG, dual channel segmentation of homologous EEG channels was carried out [9]. The segmentation was based on the difference between the spectral energies in a fixed window (beginning of a segment) and a sliding window, using T as criterion for the segmentation: M
r$lv, (0)r$?.*(0)
-
sm
The partial correlation matrix of order k is then given by
(v,”’)~
Pk
= ~,‘/’~[ef(l)e;~(l)]
(4)
where U:/’ and V k / 2 are the Cholesky factors of the covariance matrices u k and vk of e { ( l ) and e : ( l ) , respectively. Let q k = E { y ( l ) [ ~ , l / ’ e ~ ( l ) l r )It. may be shown that
and hence from (4)
The equations above are the basis of the estimation method. Let n denote the sample size and p the order of the model to be fitted, the estimator g k of (n - P)”’qkr can be obtained by computing the Cholesky factor Lp, of the mp by mp matrix ypy,‘, where
and then solving the equation
(41 . g& = - [ y ( p + 1) *
Also, the estimator
ir, of ( n - p )
*
y ( n ) l ( y p ) T . (8)
is obtained from
with
Finally, Pk is estimated by
k
C Af(j)Y(l-j ) j= I
= -
k
9; (1)
4 (0= Y ( 0
(1)
) the zero-lag autocorrelawhere r{ly,(0)and T - $ ~ ~ ( Oare tions o f the fixed window for the two channels, and r{,>?(k)and rGI7(k) are the k-lag cross correlations for the fixed and sliding windows, respectively. The sliding window was moved on the signals sequentially in steps of one sample, and whenever Texceeded a preset threshold value a segment boundary was defined, and the procedure of segmentation restarted. Typically, the EEG signal was divided into uneven segments of 1-2 s length, and the vectorial AR modeling applied to each segment for computation of the power spectra. A second scheme for computation of the power spectra through the multivariate AR modeling was based on using a sliding window of 1 s length, and estimating the spectra as a function of time, with steps of 0.25-0.5 s . Power spectrum estimates and coherence, and phase spectra (@,(f) = tan-’ [Q,,/C,,(f)I~ where C,,(f> and Q x , . ( f )are the real and imaginary part of the crosscorrelation function, respectively) were calculated on each of the EEG segments, using multichannel (vectorial) autoregressive modeling. For short record length, i.e, for EEG signal segments 3 s long or less, the performance of the commonly used covariance method, applying the vectorial version of the Levinson algorithm to solve the YuleWalker equations [ 111, [ 121 is not good enough [ 101, and a better spectral estimator is required. The Dickinson method [13], [14], based on residual energy ratios, estimating the partial correlation matrices using Cholesky decomposition, was found to be the method of choice. The Dickinson method estimates directly the partial correlation matrices, and then uses the vectorial version of the Levinson algorithm to compute the AR parameters. Let y ( l ) be the rn-vector valued observed process, and denote the forward and backward prediction of y ( l ) by
sm
where A { ( J ) and A,b ( j ) are the forward and backward AR coefficient matrices of the kth order model, respectively. The corresponding prediction errors are given by
= -
CI A ; ( j ) y ( l + j )
J =
(2)
From the above estimated partial correlation matrices, the forward and backward autoregressive matrix coefficients are computed using the vectorial version of the Levinson algorithm.
954
’
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 9, SEPTEMBER 1992
111. RESULTS I5O
Figs. 1-3 demonstrate the superiority of the multivariate AR spectral estimator over the FFT for short record length. Sharp distinct spectral peaks are discerned in the PSD computed from 1 s long signal segments, using the parametric modeling, as compared to the more noisy PSD computed via the FFT on a stretch of 8 s long signal (averaging eight periodograms), Fig. 1. In Fig. 2, a stretch of 1.6 s long EEG signal and it’s spectra are demonstrated. Due to the short record length, averaging of the periodogram was not carried out, and therefore, the spectral estimator (FFT) is inconsistent, and it’s standard deviation is as great as the quantity being measured. The PSD calculated via AR modeling gives a consistent estimate, and it’s standard error is approximately equivalent to that of averaging over 13 periodograms (each derived from 1.6 s long signal). The same can be seen for computation of coherence spectra, Fig. 3 . The coherence spectrum calculated using the FFT on a signal stretch of 2 s (averaging of two periodograms) shows that only two of it’s peaks are significantly different from zero coherence. Therefore, all other peaks are considered to be spurious, and the coherence spectrum noisy. For the two significant peaks, at 10 and 15 Hz, the 99% confidence interval gives low coherence values of 0.43 and 0.56, respectively. In order to obtain by the FFT estimator coherence peaks that are significantly different from zero coherence, averaging of several periodograms is necessary, Fig. 3 (left), and this will cause masking of the short-time dynamic changes which might occur in the epileptic EEG. On the other hand, coherence spectrum calculated through multivariate AR modeling on the 2 s long EEG signal discerns significant peaks at 3.25, 10, and 18.5 Hz. The 99% confidence interval for the main peak of 3.25 Hz gives a coherence of 0.98, reflecting the synchronized activity of the two signals at this frequency. Comparison between the commonly used covariance method for the AR modeling, using the Levinson algorithm for solving the YuleWalker equations, and the Dickinson method (both with a model order 8, chosen empirically as the best one), demonstrates the better efficiency of the latter in case of short record length, resolving sharp spectral peaks, Fig. 4. A stretch of rat multichannel EEG signal showing a generalized epileptic seizure, recorded acutely after the microinjection of Kainic acid to the right amygdala, is given in Fig. 5. It can be seen that the signal is highly nonstationary, and the dynamic changes are to be followed by adaptive segmentation, and computation of the spectra, applying the multivariate AR model to each of the short segments. Three-D plots of the power spectra, Fig. 6, and contour maps, Fig. 7, demonstrate that the epileptic activity which starts with a dominant frequency at around 7 Hz undergoes a slow decrease to below 5 Hz. Coherence spectra, depicted in Fig. 8, show that synchro-
m
FREQUENCY
FREQUENCY Hz 0.5 I
yz
I
FREQUENCY Hz
FREQUENCY HZ
Fig. I . Comparison between power spectra computed via the FFT and via parametric modeling of the EEG. Seizure EEG (left frontal lead) from a rat treated acutely with Kainic acid to the right amygdala. Upper left: power spectrum density of an 8 s long signal stretch computed via FFT. Averaging of the eight periodograms calculated on signal segments 1 s long each. Upper right and lower left and right: power spectrum calculated on 3 different 1 s long segments (included in the 8 s long signal stretch) through AR modeling. Model order 8, sampling rate 128 Hz. The y-axis is in arbitrary units.
-400‘ 0
0.5
1
1.5
2I
TIME SEC. 102
10‘ 100
lo-‘
-
10-2
0
10 20 30 FREQUENCY Hz
Fig. 2. Comparison between power spectra computed via FFT and via parametric modeling of the EEG. Above: 1.6 s long EEG stretch derived by adaptive segmentation, recording conditions as in Fig. 1. Below: Broken line-raw periodogram (FFT) of the EEG signal displayed above. No averaging of the periodogram was possible due to the short length of the signal at hand. Thus, the FFT estimate is inconsistent, and it’s standard deviation is as great as the quantity being measured. Solid line-power spectrum calculated on the EEG signal through AR modeling as in Fig. 1 . The estimate is consistent, and it’s standard error is approximately equivalent to that of averaging over 13 periodograms derived from 1.6 signal segments.
GATH
('I (11
955
TRACKING OF RAPID DYNAMlC CHANGES IN SEIZURE EEG
I
I
L w
w
ti
ti
w
I w
0
0
- 20
0
s
IS
10
-
0.1 0
20
5
10
15
Fig. 3 . Comparison between coherence spectra calculated via the FFT and via parametric modeling of the EEG. Same recording conditions as in Fig. I . First channel: left amygdala. second channel: right occipital lead. Broken line, left: coherence spectrum calculated via the FFT on a signal stretch o f 8 s (average of eight periodograms). Broken line, right: coherence spectrum calculated via the FFT on a signal stretch of 2 s (average of two periodograms). Solid lines left and right: coherence spectrum calculated using multichannel AR modeling, signal segment of 2 s length included in the 8 s long signal stretch. The clear peaks discerned in the coherence spectrum calculated using AR modeling are in agreement with the frequency content of the signals, whereas only two of the coherence peaks (at I O and IS Hz) calculated via the FFT on the 2 s long EEG signal are significantly different from zero coherence. Therefore, all other peaks in the coherence spectra calculated using the FFT are considered to be spurious. The 99% confidence interval for the two significant peaks gives low coherence values.
I!\
8
6
4
14
12
10
16
TIME SEC
Fig. 5 . A 16 s long segment of seizure EEG rat acutely treated with Kainic acid to the right amygdala. The channels from above are: I ) right amygdala. 2) left amygdala. 3) left frontal lead. 4) left occipital lead.
Fig. 6 . 3-D plot of the power spectrum density (multichannel Ar modeling) of the signal in Fig. 5. Upper left: right amygdala. Upper right: left amygdala. Lower left: left frontal lead. Lower right: left occipital lead.
0.2
-
p i
20
FREQUENCY Hz
FREQUENCY H r
I 0
0 0
IO
20
30
40
FREQUENCY Hz
Fig. 4. Comparison between the commonly used covariance method for computation of the AR parameters, and Dickinson's method. Recording conditions as in Figs. 1-3, left frontal lead, record length 182 samples (1.4 5 ) . Solid line: power spectrum density calculated using the Dickinson's method. Broken line: power spectrum density calculated using the Levinson algorithm for solving the Yule-Walker equations. The y-axis is in arbitrary units. The Dickinson method gives sharper spectral peaks. For further explanation see text.
nization between right and left amygdala is leading that between the left frontal and left occi~italleads. Also the frequency for maximum coherence for the fight and left from around HZ at amygdala decreases during Some the beginning of the seizure to around 5 Hz. The decrease
0
5
10
15
20
0
5
FREQUENCY Hz
0' 0
5
IO
15
FREQUENCY Hz
10
15
20
FREQUENCY Hz
I 20
0' 0
5
IO
15
20
FREQUENCY HZ
Fig. 7. Contour maps corresponding to the spectra in Fig. 6 (the signal in Fig. 5). The dominant frequency of the seizure is first seen at about 7 s , and is gradually decreasing to below 5 Hz. A sharp second harmonic is seen at the occipital lead.
956
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 9, SEPTEMBER 1992
herence for right amygdaldleft amygdala, and right 7 amygdala/right hippocampus leads, Fig. 12. All coher-
i
ence values at points of maximum coherence where greater than 0.9, and significant, p < 0.01.
I
01 0
5
I
OL
15
10
0
5
10
15
FREQUENCY Hz
FREQUENCY Hz
Fig. 8. Contour maps corresponding to 3-D plots of coherence spectra for the signal in Figs. 5-7. Left: coherence between right and left amygdala leads. Right: coherence hetween left frontal and left occipital leads. The arrows show the instants where a clear increase of coherence for the dominant seizure frequency is taking place. The synchronization between the rightileft amygdala is leading that between the left frontalioccipital leads.
I
1
I
I
0
2
4
6
8
10
12
14
16
TIYE SEC.
Fig. 9. A sample from depth electrode recording from a patient having generalized epileptic seizure with focal onset at the right temporal lobe. The first 16 s from the start of the seizure. The channels from above are: 1) Right amygdala. 2) Right hippocampus. 3) Left amygdala.
in the frequency of maximum coherence for the left frontallleft occipital leads for the same time period is less conspicuous, Fig. 8. A segment showing the early part of a focal seizure, from depth electrode recording of a patient with temporal lobe epilepsy is given in Fig. 9. Also here the signal is highly nonstationary, and the dynamic changes can be traced by looking at the autospectra, and coherence spectra of short record length, Figs. 10 and 11. The seizure starts with a burst of high frequency, at around 10 Hz, clearly depicted in both amygdala leads, but lacking in the right hippocampus, Fig. 10. The instantaneous spectral peaks corresponding to the dominant frequency of the seizure were found to stabilize at around 5 Hz, after a gradual decrease in frequency from about 7 Hz, during some 25 seconds, Figs. 10 and 11. Coherence between all three channels increased during the early part of the seizure. Complete synchronization between all three channels was obtained some 10 s after start of the seizure, as can be seen by plotting the frequency of maximum co-
IV. CONCLUSIONS Due to the necessity of frequency smoothing or segment averaging inherent in the FFT technique for spectral estimation [ 151, this technique is not suitable for the detection of rapid dynamic changes that might occur during preictal and ictal EEG. Parametric spectral estimation, on the other hand, would give better spectral estimate than the FFT for the same record length. The improvement is related to the number of degrees of freedom of the AR model, which is given by N / p where N is the number of signal samples and p is the model order, and the asymptotic variance of the AR spectral estimate is similar to that of the smoothed periodogram with the same number of degrees of freedom [ 161. Thus, for a record length of 128 samples (1 second) and model order 8, the number of degrees of freedom for the AR spectral estimate will be 16. In that case, the standard error of the AR spectral estimate will approximately be equivalent to that of applying the FFT to 8 EEG segments of 1 s length, and averaging over all eight periodograms. A good frequencykime resolution is a prerequisite for tracking the dynamic ictal changes. Thus, for the multivariate AR modeling of the epileptic EEG, the Dickinson . . method [14] is in pirticular attractive, since for short record length it would give better estimates than the vectorial version of the Levinson algorithm [lo]. Compared to the multivariate time varying AR modeling (MVTVAR) [6], the present method of combining adaptive segmentation [9] with vectorial AR spectral estimation is more efficient from the computational point of view. The inherent problem of the MVTVAR is the tendency to overparametrization, the number of parameters increasing as the square of the number of EEG channels, and thus, this method would be limited to modeling of a small number of EEG channels. A different approach for processing nonstationary (time-varying) signals is based on a conjoint time-frequency representation using the WignerVille distribution [17], [18]. However, it has been shown that the value of the Wigner distribution at a certain point is influenced by the rest of the signal in such a way that very often overwhelms what is happening near that point itself [ 191. Smoothing might alleviate this problem, but will induce decrease of the resolution power for detection of refined structures [18]. The examples of epileptic EEG given in the present paper, both from rats having generalized epileptic seizures with focal onset, and from depth recordings in humans with temporal lobe epilepsy, show that the analysis procedure described above give quantitative information related to the dynamic changes which occur prior and during the seizure. Conclusions as to pathophysiological
957
GATH er a l . : TRACKING OF RAPID D Y N A M I C C H A N G E S IN SEIZURE EEG
16 3
0
Fig. 10. Three-D plots of the power spectrum density corresponding to the signal in Fig. 9 (the first 46.3 s from the start of the seizure). Adaptive segmentation, AR model order 8 . Upper left: right hippocampus, Upper right: right amygdala. Lower left: left amygdala. The seizure starts with a sharp peak at about 10 Hz in both amygdala, lacking at the right hippocampus.
I
''i
"i
35:
30
1
2
-
1" I
35 -I
i
-0 0
2
I
6
I
10
LL.
I.
,6
18
20
FRCPl E Y C Y H:
Fig. 11. Contour maps corresponding to the spectra in Fig. 10. The dominant seizure frequency is at about 7 Hz, decreasing slowly to around 5 Hz.
mechanisms involved in the propagation of focal seizures, based on power spectrum estimation using the FFT technique [2], [3], will have to be reevaluated. The reason is, as has been shown in the present paper, the "smearing" of the rapid dynamic changes that occur during the sei-
zure, when using the FFT technique. It is thus concluded that power spectrum estimation together with estimation of coherence and phase spectra, using parametric multichannel modeling of the EEG, could be an aid to the quantitative analysis of EEG during epileptic seizures.
958
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL 39. NO 9. SEPTEMBER 1992
vol. 35, 217-250, 276-300, 372-389, 1980. [I81 W. Martin and P. Flandrin, “Wigner-Ville spectral analysis of nonstationary process,” IEEE Trans. Acoust. , Sound, Signal Process.,
I2 -
0
2
4
6
8
10
12
14
I6
TIME SEC.
Fig. 12. Frequency for maximum coherence for right amygdalaileft amygdala leads (solid line) and for right amygdalairight hippocampus leads (dotted line) of the signal in Fig. 9. Multichannel AR modeling, sliding window of 1 s length, step size 0.5 s. Synchronization between al three channels (at about 7 Hz) is obtained after some 10 s from the start of the seizure.
Isak G a t h received the M D degree trom the Hadassah Medical School, Hebrew University. Jerusalem, Israel in 1964, the D 1 C degree in biomedical engineering trom Imperial College ot Science and Technology, London, Engldnd, in 1970, and the D Sc degree i n Electrical Engineering from the Technion-Israel Institute 01 Technology, Haifa, Israel, in 1975 In 1975 he joined the Technion-Iwael Institute of Technology, where he is currently Professor of Biomedical Engineering. During 1981 he was a research associate at the Monakow laboratory, Umversity of Zurich, Switzerland, and during 1989- 1990 Professor of Applied Mathematics. Depdrtment of Computer Science and Applied Mathematics, University of Grenoble, France His major areas of research interest are signal processing dnd pattern recognition applied to biological 5igndls
REFERENCES [ I ] M . B. A. Brazier, “Spread of seizure discharge in epilepsy: Anatomical and electrophysiological considerations,” Exp. Neurol., vol. 3 6 , pp.263-272, 1973. [2] J . Gotman, “Measurement of small time differences between EEG channels: Method and application to epileptic seizure propagation,” Electroencephal. clin. Neurophysiol., vol. 56, pp. 501-514, 1983. 131 J. P. Lieb and C. Skomer, “Interhemispheric propagation of limbic seizures: A coherenceiphase analysis,” Electroencephul. clin. Neurophysiol., vol. 61, S235, 1985. 141 T. Darcy and P. Williamson, “Spatio-temporal analysis of human intracranially recorded epileptic seizures: Power spectra and derived measures,” Electroencephal. clin. Neurophysiol. vol. 61, 573-587, 1985. [ 5 ] I. Gotman, “Interhemispheric interactions in seizures of focal onset: data from human intracranial recordings,” Electroencephal. clin. Neurophysiol., vol. 67, 120-133, 1987. [6] W. Gersch, “Non stationary multichannel time series analysis,“in EEG Handbook Revised Series, Vol. 11, A. S . Gevins and A. Remond, Eds. 1987, pp. 261-296. [7] G. Bodenstein and H. M. Praetorius, “Feature extraction from EEG by adaptive segmentation,” Proc. IEEE, vol. 65, pp. 642-652, 1977. [8] I. Gath and E. Bar-On, “Computerized method for scoring of polygraphic sleep recordings,” Comp. Progr. Biomed.. vol. 1 I , pp. 217223, 1980. [9] I. Gath, A. Michaeli, and C. Feuerstein, “A model for dual channel segmentation of the EEG signal,” B i d . Cybrrn., vol. 64, pp. 225230, 1991. [IO] D. T . Pham and D. Q . Tong, “Maximum likelihood estimation for multivariate autoregressive model,” Research Rep. RR 807-M, Laboratoire LMCIIMAG, Universite Joseph Fourier, Grenoble. 1990. [ I I ] P. Whittle, “On the fitting of multivariate autoregressions and the approximate canonical factorization of a spectral density matrix,” Biomrtrika. vol. 50. pp. 129-134, 1963. [I21 U. A. Wiggins and E. A . Robinson, “Recursive solution to the multichannel filtering problem,’’ J . Geophysical R e s . , vol. 70, pp. 18851891, 1965. [ 131 B. W. Dickinson, “Autoregressive estimation using residual energy ratios,” IEEE Trans. Inform Theory, vol. 24, pp. 503-506, 1978. 1141 B. W . Dickinson, “Estimation of partial correlation matrices using Cholesky decomposition,” IEEE Trans. Automat. Contr., vol. 24, pp. 302-305, 1979. [IS] J . S. Bendat and A. G. Piersol, Random Data: Anulysis and Meusurements Procedures. New York: Wiley, 1971. [ 161 K . N. Berk, “Consistent autoregressive spectral estimates,” Ann. s t a t . , V O I . 2, pp. 489-502, 1974. 1171 T . A. C . Classen and W. F. G . Mecklenbauker, “The Wigner dis-
Claude Feuerstein was born in Grenoble, France, in 1951. He received the M.D. degree in 1977 from the University of Grenoble, France. Through a combined M.D.-Ph.D. program in Neurophysiology and short training periods in Sweden, US, and Canada, he received the “Diplbm d‘Etudes et de Recherche en Biologie Humain” at Lyon University in 1982, and the “Diplbm d’Habilitation A Diriger des Recherches” at Grenoble University in 1986. From 1978-1981 he was assistant, received tenure in 1981, and is currently Professor of Physiology at the University of Grenoble since 1990. His research interest focus on dopamine neurons, and on clinical neurophysiology. in particular sleep and presurgical evaluation in epilepsy.
Dinh T u a n P h a m was born in Hanoi, Vietnam, in 1945. He graduated from the School of Applied Mathematics and Computer Science (ENSIMAG) of the National Polytechnic Institute of Grenoble (INPG), and received the Ph.D. degree in statistics from the University of Grenoble in 1975. He was a postdoctoral fellow at the Department of Statistics, Berkeley, Berkeley, CA. in 19771970, and a visiting professor at Indiana University, Bloomington (Dept. Mathematics) in 19791980. He is presently a senior researcher at the CNRS (The French National Center for Scientific Research) Laboratory of Modelling and Computation. His research interests include time series analysis, array processing, signal estimation and modelling. and biomedical signal processing.
G e r a r d Rondnuin was born in 1947. He received the M.D. and Ph.D. degrees at the University of Montpellier. France, in 1978 and 1991, respectively. He is currently Maitre de Conferences at the University of Montpellier. His research interests include on one side experimental epilepsies, in particular kindling models, the role o f neurotransmitters, the effect of drugs, and propagation of physiological and pathological activities, and on the other side clinical neurophysiology.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39. NO. 9, SEPTEMBER 1992
On the Tracking of Rapid Dynamic Changes in Seizure EEG Isak Gath, Claude Feuerstein, Dinh Tuan Pham, and Gerard Rondouin
Abstract-Estimation of autospectra and coherence and phase spectra of seizure EEG, using the FFT technique, will cause “smearing” of the rapid dynamic changes which occur during the seizure. This is inherent to FFT spectral estimation, due to the averaging process which is necessary in order to get consistent spectral estimates. A different approach suggested in the present study is to carry out multivariate autoregressive modeling of the multichannel seizure EEG, combined with adaptive segmentation. In order to obtain good estimates in cases of short record length, the vectorial AR modeling was based on residual energy ratios. The method has been tested on multichannel seizure EEG recordings from rats with focal epilepsy, caused by intracerebral administration of Kainic acid, and in depth EEG recordings in patients with temporal lobe epilepsy.
I. INTRODUCTION
P
OWER spectrum estimation is a useful tool in the quantitative analysis of EEG recorded during epileptic seizures. Autospectra, coherence and phase spectra, might be helpful in the interpretation of the multichannel dynamic changes which occur prior and during the seizures. The ictal EEG have been most often analyzed using the fast fourier transform (FFT) [1]-151. However, the “epileptic” EEG signal recorded either in experimental animals or in humans might show marked nonstationarity, and hence, estimation of the power spectrum using the FFT is not without problems. Periodograms have to be computed for short stationary signal segments, usually not exceeding 2-3 s, which means frequency resolution of 0.33-0.5 Hz. In order to obtain consistent power spectrum estimates, averaging of several periodograms is required to decrease the variability of the estimate. This averaging process inherent in the FFT technique will cause “smearing” of the dynamic changes which might occur during the seizure. An alternative approach to the estimation of the power spectrum and phase and coherence spectra is based on parametric modeling of the multichannel EEG. A time varying multivariate autoregressive (AR) modeling, based Manuscript received September 12, 1990; revised December 5 , 1991. I. Gath is with the Department of Biomedical Engineering, Technion, Israel Institute of Technology, 32000 Haifa, Israel. C. Feuerstein is with the Laboratoire de Physiologie Section Neurophysiologie, INSERM U318 CHU, Grenoble. France. D. T. Pham is with the Laboratoire de Modelisation et CalculilMAG, CNRS, and Universite Joseph Fourier, Grenoble, France. G. Rondouin is with the Laboratoire de Medecine Experimentale, INSERM U 239, University of Montpellier, Montpellier. France. IEEE Log Number 9201968.
on the deterministic-regression method using a family of basis functions have been suggested by Gersch 161. This method overcomes the problem of nonstationarity but suffers from a heavy computational load. The number of parameters fitted to the model increases as the square of the number of EEG channels [6], and hence in practice this method is useful for modeling of only a small number of EEG channels. A still different approach will be to carry out multichannel AR modeling in combination with adaptive segmentation of the EEG [7]-[9]. However, whatever the procedure, in order to trace the rapid dynamic changes inherent in the epileptic EEG, the AR model estimator is required to be efficient in case of short record lengths (2 s or less). The aim of the present paper is to investigate a procedure for parametric multivariate modeling of the epileptic EEG, a procedure aimed to detect rapid dynamic changes in the recorded EEG signal. 11. METHODS A . Data Acquisition Generalized seizures with focal onset were generated in two male Wister rats by intracerebral injection of Kainic acid. A cannula was inserted stereotaxically under anesthesia into the right basolateral amygdala, and at the same time bipolar twisted electrodes were implanted in both amygdala for continuous depth recording. Surface EEG recordings using three bipolar channels (left frontal and left and right occipital leads) were obtained by implanted gold plated screws over the respective cerebral cortex, and the electrode assembly was secured to the skull with a dental acrylic cement. Recordings were started ten days after the implantation, acute seizures were obtained by a microinjection of 0.2 pg Kainic acid into the right amygdala. Data acquisition was carried out on an IBM PClAT computer with a Data Translation A/D card, using a sampling rate of 128 Hz. In addition, recordings were carried out in several patients having temporal lobe epilepsy. Depth electrodes were placed using the stereotactic technique in both mesial and cortical structures, in particular, around the limbic system on both sides. Long-time monitoring of these patients was carried out lasting 2-14 days, until sufficient ictal data was accumulated. Typical generalized seizures having a focal onset were chosen for the
0018-9294/92$03.00
1992 IEEE
953
GATH er a l . : TRACKING OF RAPID DYNAMIC CHANGES IN SEIZURE EEG
analysis, and data acquisition was carried out as described for the rats treated with Kainic acid.
B. Data Analysis Signal segments including both normal EEG and ictal activity were chosen for analysis. Two different strategies were applied to the nonstationary seizure EEG signal. In order to obtain quasi-stationary signal segments from the nonstationary seizure EEG, dual channel segmentation of homologous EEG channels was carried out [9]. The segmentation was based on the difference between the spectral energies in a fixed window (beginning of a segment) and a sliding window, using T as criterion for the segmentation: M
r$lv, (0)r$?.*(0)
-
sm
The partial correlation matrix of order k is then given by
(v,”’)~
Pk
= ~,‘/’~[ef(l)e;~(l)]
(4)
where U:/’ and V k / 2 are the Cholesky factors of the covariance matrices u k and vk of e { ( l ) and e : ( l ) , respectively. Let q k = E { y ( l ) [ ~ , l / ’ e ~ ( l ) l r )It. may be shown that
and hence from (4)
The equations above are the basis of the estimation method. Let n denote the sample size and p the order of the model to be fitted, the estimator g k of (n - P)”’qkr can be obtained by computing the Cholesky factor Lp, of the mp by mp matrix ypy,‘, where
and then solving the equation
(41 . g& = - [ y ( p + 1) *
Also, the estimator
ir, of ( n - p )
*
y ( n ) l ( y p ) T . (8)
is obtained from
with
Finally, Pk is estimated by
k
C Af(j)Y(l-j ) j= I
= -
k
9; (1)
4 (0= Y ( 0
(1)
) the zero-lag autocorrelawhere r{ly,(0)and T - $ ~ ~ ( Oare tions o f the fixed window for the two channels, and r{,>?(k)and rGI7(k) are the k-lag cross correlations for the fixed and sliding windows, respectively. The sliding window was moved on the signals sequentially in steps of one sample, and whenever Texceeded a preset threshold value a segment boundary was defined, and the procedure of segmentation restarted. Typically, the EEG signal was divided into uneven segments of 1-2 s length, and the vectorial AR modeling applied to each segment for computation of the power spectra. A second scheme for computation of the power spectra through the multivariate AR modeling was based on using a sliding window of 1 s length, and estimating the spectra as a function of time, with steps of 0.25-0.5 s . Power spectrum estimates and coherence, and phase spectra (@,(f) = tan-’ [Q,,/C,,(f)I~ where C,,(f> and Q x , . ( f )are the real and imaginary part of the crosscorrelation function, respectively) were calculated on each of the EEG segments, using multichannel (vectorial) autoregressive modeling. For short record length, i.e, for EEG signal segments 3 s long or less, the performance of the commonly used covariance method, applying the vectorial version of the Levinson algorithm to solve the YuleWalker equations [ 111, [ 121 is not good enough [ 101, and a better spectral estimator is required. The Dickinson method [13], [14], based on residual energy ratios, estimating the partial correlation matrices using Cholesky decomposition, was found to be the method of choice. The Dickinson method estimates directly the partial correlation matrices, and then uses the vectorial version of the Levinson algorithm to compute the AR parameters. Let y ( l ) be the rn-vector valued observed process, and denote the forward and backward prediction of y ( l ) by
sm
where A { ( J ) and A,b ( j ) are the forward and backward AR coefficient matrices of the kth order model, respectively. The corresponding prediction errors are given by
= -
CI A ; ( j ) y ( l + j )
J =
(2)
From the above estimated partial correlation matrices, the forward and backward autoregressive matrix coefficients are computed using the vectorial version of the Levinson algorithm.
954
’
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 9, SEPTEMBER 1992
111. RESULTS I5O
Figs. 1-3 demonstrate the superiority of the multivariate AR spectral estimator over the FFT for short record length. Sharp distinct spectral peaks are discerned in the PSD computed from 1 s long signal segments, using the parametric modeling, as compared to the more noisy PSD computed via the FFT on a stretch of 8 s long signal (averaging eight periodograms), Fig. 1. In Fig. 2, a stretch of 1.6 s long EEG signal and it’s spectra are demonstrated. Due to the short record length, averaging of the periodogram was not carried out, and therefore, the spectral estimator (FFT) is inconsistent, and it’s standard deviation is as great as the quantity being measured. The PSD calculated via AR modeling gives a consistent estimate, and it’s standard error is approximately equivalent to that of averaging over 13 periodograms (each derived from 1.6 s long signal). The same can be seen for computation of coherence spectra, Fig. 3 . The coherence spectrum calculated using the FFT on a signal stretch of 2 s (averaging of two periodograms) shows that only two of it’s peaks are significantly different from zero coherence. Therefore, all other peaks are considered to be spurious, and the coherence spectrum noisy. For the two significant peaks, at 10 and 15 Hz, the 99% confidence interval gives low coherence values of 0.43 and 0.56, respectively. In order to obtain by the FFT estimator coherence peaks that are significantly different from zero coherence, averaging of several periodograms is necessary, Fig. 3 (left), and this will cause masking of the short-time dynamic changes which might occur in the epileptic EEG. On the other hand, coherence spectrum calculated through multivariate AR modeling on the 2 s long EEG signal discerns significant peaks at 3.25, 10, and 18.5 Hz. The 99% confidence interval for the main peak of 3.25 Hz gives a coherence of 0.98, reflecting the synchronized activity of the two signals at this frequency. Comparison between the commonly used covariance method for the AR modeling, using the Levinson algorithm for solving the YuleWalker equations, and the Dickinson method (both with a model order 8, chosen empirically as the best one), demonstrates the better efficiency of the latter in case of short record length, resolving sharp spectral peaks, Fig. 4. A stretch of rat multichannel EEG signal showing a generalized epileptic seizure, recorded acutely after the microinjection of Kainic acid to the right amygdala, is given in Fig. 5. It can be seen that the signal is highly nonstationary, and the dynamic changes are to be followed by adaptive segmentation, and computation of the spectra, applying the multivariate AR model to each of the short segments. Three-D plots of the power spectra, Fig. 6, and contour maps, Fig. 7, demonstrate that the epileptic activity which starts with a dominant frequency at around 7 Hz undergoes a slow decrease to below 5 Hz. Coherence spectra, depicted in Fig. 8, show that synchro-
m
FREQUENCY
FREQUENCY Hz 0.5 I
yz
I
FREQUENCY Hz
FREQUENCY HZ
Fig. I . Comparison between power spectra computed via the FFT and via parametric modeling of the EEG. Seizure EEG (left frontal lead) from a rat treated acutely with Kainic acid to the right amygdala. Upper left: power spectrum density of an 8 s long signal stretch computed via FFT. Averaging of the eight periodograms calculated on signal segments 1 s long each. Upper right and lower left and right: power spectrum calculated on 3 different 1 s long segments (included in the 8 s long signal stretch) through AR modeling. Model order 8, sampling rate 128 Hz. The y-axis is in arbitrary units.
-400‘ 0
0.5
1
1.5
2I
TIME SEC. 102
10‘ 100
lo-‘
-
10-2
0
10 20 30 FREQUENCY Hz
Fig. 2. Comparison between power spectra computed via FFT and via parametric modeling of the EEG. Above: 1.6 s long EEG stretch derived by adaptive segmentation, recording conditions as in Fig. 1. Below: Broken line-raw periodogram (FFT) of the EEG signal displayed above. No averaging of the periodogram was possible due to the short length of the signal at hand. Thus, the FFT estimate is inconsistent, and it’s standard deviation is as great as the quantity being measured. Solid line-power spectrum calculated on the EEG signal through AR modeling as in Fig. 1 . The estimate is consistent, and it’s standard error is approximately equivalent to that of averaging over 13 periodograms derived from 1.6 signal segments.
GATH
('I (11
955
TRACKING OF RAPID DYNAMlC CHANGES IN SEIZURE EEG
I
I
L w
w
ti
ti
w
I w
0
0
- 20
0
s
IS
10
-
0.1 0
20
5
10
15
Fig. 3 . Comparison between coherence spectra calculated via the FFT and via parametric modeling of the EEG. Same recording conditions as in Fig. I . First channel: left amygdala. second channel: right occipital lead. Broken line, left: coherence spectrum calculated via the FFT on a signal stretch o f 8 s (average of eight periodograms). Broken line, right: coherence spectrum calculated via the FFT on a signal stretch of 2 s (average of two periodograms). Solid lines left and right: coherence spectrum calculated using multichannel AR modeling, signal segment of 2 s length included in the 8 s long signal stretch. The clear peaks discerned in the coherence spectrum calculated using AR modeling are in agreement with the frequency content of the signals, whereas only two of the coherence peaks (at I O and IS Hz) calculated via the FFT on the 2 s long EEG signal are significantly different from zero coherence. Therefore, all other peaks in the coherence spectra calculated using the FFT are considered to be spurious. The 99% confidence interval for the two significant peaks gives low coherence values.
I!\
8
6
4
14
12
10
16
TIME SEC
Fig. 5 . A 16 s long segment of seizure EEG rat acutely treated with Kainic acid to the right amygdala. The channels from above are: I ) right amygdala. 2) left amygdala. 3) left frontal lead. 4) left occipital lead.
Fig. 6 . 3-D plot of the power spectrum density (multichannel Ar modeling) of the signal in Fig. 5. Upper left: right amygdala. Upper right: left amygdala. Lower left: left frontal lead. Lower right: left occipital lead.
0.2
-
p i
20
FREQUENCY Hz
FREQUENCY H r
I 0
0 0
IO
20
30
40
FREQUENCY Hz
Fig. 4. Comparison between the commonly used covariance method for computation of the AR parameters, and Dickinson's method. Recording conditions as in Figs. 1-3, left frontal lead, record length 182 samples (1.4 5 ) . Solid line: power spectrum density calculated using the Dickinson's method. Broken line: power spectrum density calculated using the Levinson algorithm for solving the Yule-Walker equations. The y-axis is in arbitrary units. The Dickinson method gives sharper spectral peaks. For further explanation see text.
nization between right and left amygdala is leading that between the left frontal and left occi~italleads. Also the frequency for maximum coherence for the fight and left from around HZ at amygdala decreases during Some the beginning of the seizure to around 5 Hz. The decrease
0
5
10
15
20
0
5
FREQUENCY Hz
0' 0
5
IO
15
FREQUENCY Hz
10
15
20
FREQUENCY Hz
I 20
0' 0
5
IO
15
20
FREQUENCY HZ
Fig. 7. Contour maps corresponding to the spectra in Fig. 6 (the signal in Fig. 5). The dominant frequency of the seizure is first seen at about 7 s , and is gradually decreasing to below 5 Hz. A sharp second harmonic is seen at the occipital lead.
956
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 9, SEPTEMBER 1992
herence for right amygdaldleft amygdala, and right 7 amygdala/right hippocampus leads, Fig. 12. All coher-
i
ence values at points of maximum coherence where greater than 0.9, and significant, p < 0.01.
I
01 0
5
I
OL
15
10
0
5
10
15
FREQUENCY Hz
FREQUENCY Hz
Fig. 8. Contour maps corresponding to 3-D plots of coherence spectra for the signal in Figs. 5-7. Left: coherence between right and left amygdala leads. Right: coherence hetween left frontal and left occipital leads. The arrows show the instants where a clear increase of coherence for the dominant seizure frequency is taking place. The synchronization between the rightileft amygdala is leading that between the left frontalioccipital leads.
I
1
I
I
0
2
4
6
8
10
12
14
16
TIYE SEC.
Fig. 9. A sample from depth electrode recording from a patient having generalized epileptic seizure with focal onset at the right temporal lobe. The first 16 s from the start of the seizure. The channels from above are: 1) Right amygdala. 2) Right hippocampus. 3) Left amygdala.
in the frequency of maximum coherence for the left frontallleft occipital leads for the same time period is less conspicuous, Fig. 8. A segment showing the early part of a focal seizure, from depth electrode recording of a patient with temporal lobe epilepsy is given in Fig. 9. Also here the signal is highly nonstationary, and the dynamic changes can be traced by looking at the autospectra, and coherence spectra of short record length, Figs. 10 and 11. The seizure starts with a burst of high frequency, at around 10 Hz, clearly depicted in both amygdala leads, but lacking in the right hippocampus, Fig. 10. The instantaneous spectral peaks corresponding to the dominant frequency of the seizure were found to stabilize at around 5 Hz, after a gradual decrease in frequency from about 7 Hz, during some 25 seconds, Figs. 10 and 11. Coherence between all three channels increased during the early part of the seizure. Complete synchronization between all three channels was obtained some 10 s after start of the seizure, as can be seen by plotting the frequency of maximum co-
IV. CONCLUSIONS Due to the necessity of frequency smoothing or segment averaging inherent in the FFT technique for spectral estimation [ 151, this technique is not suitable for the detection of rapid dynamic changes that might occur during preictal and ictal EEG. Parametric spectral estimation, on the other hand, would give better spectral estimate than the FFT for the same record length. The improvement is related to the number of degrees of freedom of the AR model, which is given by N / p where N is the number of signal samples and p is the model order, and the asymptotic variance of the AR spectral estimate is similar to that of the smoothed periodogram with the same number of degrees of freedom [ 161. Thus, for a record length of 128 samples (1 second) and model order 8, the number of degrees of freedom for the AR spectral estimate will be 16. In that case, the standard error of the AR spectral estimate will approximately be equivalent to that of applying the FFT to 8 EEG segments of 1 s length, and averaging over all eight periodograms. A good frequencykime resolution is a prerequisite for tracking the dynamic ictal changes. Thus, for the multivariate AR modeling of the epileptic EEG, the Dickinson . . method [14] is in pirticular attractive, since for short record length it would give better estimates than the vectorial version of the Levinson algorithm [lo]. Compared to the multivariate time varying AR modeling (MVTVAR) [6], the present method of combining adaptive segmentation [9] with vectorial AR spectral estimation is more efficient from the computational point of view. The inherent problem of the MVTVAR is the tendency to overparametrization, the number of parameters increasing as the square of the number of EEG channels, and thus, this method would be limited to modeling of a small number of EEG channels. A different approach for processing nonstationary (time-varying) signals is based on a conjoint time-frequency representation using the WignerVille distribution [17], [18]. However, it has been shown that the value of the Wigner distribution at a certain point is influenced by the rest of the signal in such a way that very often overwhelms what is happening near that point itself [ 191. Smoothing might alleviate this problem, but will induce decrease of the resolution power for detection of refined structures [18]. The examples of epileptic EEG given in the present paper, both from rats having generalized epileptic seizures with focal onset, and from depth recordings in humans with temporal lobe epilepsy, show that the analysis procedure described above give quantitative information related to the dynamic changes which occur prior and during the seizure. Conclusions as to pathophysiological
957
GATH er a l . : TRACKING OF RAPID D Y N A M I C C H A N G E S IN SEIZURE EEG
16 3
0
Fig. 10. Three-D plots of the power spectrum density corresponding to the signal in Fig. 9 (the first 46.3 s from the start of the seizure). Adaptive segmentation, AR model order 8 . Upper left: right hippocampus, Upper right: right amygdala. Lower left: left amygdala. The seizure starts with a sharp peak at about 10 Hz in both amygdala, lacking at the right hippocampus.
I
''i
"i
35:
30
1
2
-
1" I
35 -I
i
-0 0
2
I
6
I
10
LL.
I.
,6
18
20
FRCPl E Y C Y H:
Fig. 11. Contour maps corresponding to the spectra in Fig. 10. The dominant seizure frequency is at about 7 Hz, decreasing slowly to around 5 Hz.
mechanisms involved in the propagation of focal seizures, based on power spectrum estimation using the FFT technique [2], [3], will have to be reevaluated. The reason is, as has been shown in the present paper, the "smearing" of the rapid dynamic changes that occur during the sei-
zure, when using the FFT technique. It is thus concluded that power spectrum estimation together with estimation of coherence and phase spectra, using parametric multichannel modeling of the EEG, could be an aid to the quantitative analysis of EEG during epileptic seizures.
958
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL 39. NO 9. SEPTEMBER 1992
vol. 35, 217-250, 276-300, 372-389, 1980. [I81 W. Martin and P. Flandrin, “Wigner-Ville spectral analysis of nonstationary process,” IEEE Trans. Acoust. , Sound, Signal Process.,
I2 -
0
2
4
6
8
10
12
14
I6
TIME SEC.
Fig. 12. Frequency for maximum coherence for right amygdalaileft amygdala leads (solid line) and for right amygdalairight hippocampus leads (dotted line) of the signal in Fig. 9. Multichannel AR modeling, sliding window of 1 s length, step size 0.5 s. Synchronization between al three channels (at about 7 Hz) is obtained after some 10 s from the start of the seizure.
Isak G a t h received the M D degree trom the Hadassah Medical School, Hebrew University. Jerusalem, Israel in 1964, the D 1 C degree in biomedical engineering trom Imperial College ot Science and Technology, London, Engldnd, in 1970, and the D Sc degree i n Electrical Engineering from the Technion-Israel Institute 01 Technology, Haifa, Israel, in 1975 In 1975 he joined the Technion-Iwael Institute of Technology, where he is currently Professor of Biomedical Engineering. During 1981 he was a research associate at the Monakow laboratory, Umversity of Zurich, Switzerland, and during 1989- 1990 Professor of Applied Mathematics. Depdrtment of Computer Science and Applied Mathematics, University of Grenoble, France His major areas of research interest are signal processing dnd pattern recognition applied to biological 5igndls
REFERENCES [ I ] M . B. A. Brazier, “Spread of seizure discharge in epilepsy: Anatomical and electrophysiological considerations,” Exp. Neurol., vol. 3 6 , pp.263-272, 1973. [2] J . Gotman, “Measurement of small time differences between EEG channels: Method and application to epileptic seizure propagation,” Electroencephal. clin. Neurophysiol., vol. 56, pp. 501-514, 1983. 131 J. P. Lieb and C. Skomer, “Interhemispheric propagation of limbic seizures: A coherenceiphase analysis,” Electroencephul. clin. Neurophysiol., vol. 61, S235, 1985. 141 T. Darcy and P. Williamson, “Spatio-temporal analysis of human intracranially recorded epileptic seizures: Power spectra and derived measures,” Electroencephal. clin. Neurophysiol. vol. 61, 573-587, 1985. [ 5 ] I. Gotman, “Interhemispheric interactions in seizures of focal onset: data from human intracranial recordings,” Electroencephal. clin. Neurophysiol., vol. 67, 120-133, 1987. [6] W. Gersch, “Non stationary multichannel time series analysis,“in EEG Handbook Revised Series, Vol. 11, A. S . Gevins and A. Remond, Eds. 1987, pp. 261-296. [7] G. Bodenstein and H. M. Praetorius, “Feature extraction from EEG by adaptive segmentation,” Proc. IEEE, vol. 65, pp. 642-652, 1977. [8] I. Gath and E. Bar-On, “Computerized method for scoring of polygraphic sleep recordings,” Comp. Progr. Biomed.. vol. 1 I , pp. 217223, 1980. [9] I. Gath, A. Michaeli, and C. Feuerstein, “A model for dual channel segmentation of the EEG signal,” B i d . Cybrrn., vol. 64, pp. 225230, 1991. [IO] D. T . Pham and D. Q . Tong, “Maximum likelihood estimation for multivariate autoregressive model,” Research Rep. RR 807-M, Laboratoire LMCIIMAG, Universite Joseph Fourier, Grenoble. 1990. [ I I ] P. Whittle, “On the fitting of multivariate autoregressions and the approximate canonical factorization of a spectral density matrix,” Biomrtrika. vol. 50. pp. 129-134, 1963. [I21 U. A. Wiggins and E. A . Robinson, “Recursive solution to the multichannel filtering problem,’’ J . Geophysical R e s . , vol. 70, pp. 18851891, 1965. [ 131 B. W. Dickinson, “Autoregressive estimation using residual energy ratios,” IEEE Trans. Inform Theory, vol. 24, pp. 503-506, 1978. 1141 B. W . Dickinson, “Estimation of partial correlation matrices using Cholesky decomposition,” IEEE Trans. Automat. Contr., vol. 24, pp. 302-305, 1979. [IS] J . S. Bendat and A. G. Piersol, Random Data: Anulysis and Meusurements Procedures. New York: Wiley, 1971. [ 161 K . N. Berk, “Consistent autoregressive spectral estimates,” Ann. s t a t . , V O I . 2, pp. 489-502, 1974. 1171 T . A. C . Classen and W. F. G . Mecklenbauker, “The Wigner dis-
Claude Feuerstein was born in Grenoble, France, in 1951. He received the M.D. degree in 1977 from the University of Grenoble, France. Through a combined M.D.-Ph.D. program in Neurophysiology and short training periods in Sweden, US, and Canada, he received the “Diplbm d‘Etudes et de Recherche en Biologie Humain” at Lyon University in 1982, and the “Diplbm d’Habilitation A Diriger des Recherches” at Grenoble University in 1986. From 1978-1981 he was assistant, received tenure in 1981, and is currently Professor of Physiology at the University of Grenoble since 1990. His research interest focus on dopamine neurons, and on clinical neurophysiology. in particular sleep and presurgical evaluation in epilepsy.
Dinh T u a n P h a m was born in Hanoi, Vietnam, in 1945. He graduated from the School of Applied Mathematics and Computer Science (ENSIMAG) of the National Polytechnic Institute of Grenoble (INPG), and received the Ph.D. degree in statistics from the University of Grenoble in 1975. He was a postdoctoral fellow at the Department of Statistics, Berkeley, Berkeley, CA. in 19771970, and a visiting professor at Indiana University, Bloomington (Dept. Mathematics) in 19791980. He is presently a senior researcher at the CNRS (The French National Center for Scientific Research) Laboratory of Modelling and Computation. His research interests include time series analysis, array processing, signal estimation and modelling. and biomedical signal processing.
G e r a r d Rondnuin was born in 1947. He received the M.D. and Ph.D. degrees at the University of Montpellier. France, in 1978 and 1991, respectively. He is currently Maitre de Conferences at the University of Montpellier. His research interests include on one side experimental epilepsies, in particular kindling models, the role o f neurotransmitters, the effect of drugs, and propagation of physiological and pathological activities, and on the other side clinical neurophysiology.
