- Email: [email protected]
Adaptive time–frequency parametrization in pharmaco EEG
Journal of Neuroscience Methods 117 (2002) 65 /71 www.elsevier.com/locate/jneumeth
Adaptive time frequency parametrization in pharmaco EEG /
Piotr J. Durka a,*,1, Waldemar Szelenberger b, Katarzyna J. Blinowska a, Wojciech Androsiuk b, Maciej Myszka b a
Laboratory of Medical Physics, Institute of Experimental Physics, Warsaw University, ul. Hoz˙a 69, 00-681 Warszawa, Poland b Department of Psychiatry, Medical University of Warsaw, ul. Nowowiejska 27, 00-665 Warszawa, Poland Received 18 January 2002; received in revised form 15 March 2002; accepted 18 March 2002
Abstract Adaptive time /frequency approximations offer description of the local structures of a signal in terms of their time and frequency coordinates, widths and amplitudes. These parameters can then be used to select and study electroencephalogram (EEG) structures like sleep spindles or slow wave activity (SWA) with high resolution. Such a detailed description of relevant structures improves on the sensitivity of the traditionally used spectral power estimates and opens new possibilities of investigation. These advantages are illustrated using a double-blind test of the influence of zolpidem and midazolam on sleep EEG, and the results are compared with the traditional approach. The observed decrease of frequency of the SWA under the influence of sleep-inducing drugs gives an example of an effect elusive to classical methodology. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Matching pursuit; Time /frequency parametrization; Pharmaco EEG; Sleep spindles; Slow wave activity; Zolpidem; Midazolam
1. Introduction The estimation of spectral power in fixed frequency bands by means of the Fourier transform is one of the very few signal processing methods that is widely adopted in the practice of analyzing the electroencephalogram (EEG). Its general acceptance is due to a long record of applications (since the first application of Fourier transform to EEG (Dietsch, 1932)), fast numerical implementations (since the introduction of the fast Fourier transform by Cooley and Tukey (1965)) and the relatively simple interpretation of results in terms of oscillations and cycles per second. In contrast, none of the advanced signal processing methods proposed in the last few decades have gained any comparable acceptance.
* Corresponding author. Tel.: /48-22-553-2126; fax: /48-22-6226154; http://brain.fuw.edu.pl. E-mail addresses: [email protected] (P.J. Durka), [email protected] (W. Szelenberger), [email protected] (K.J. Blinowska), [email protected] (W. Androsiuk), maciek@ psych.waw.pl (M. Myszka). 1 Partially supported by Committee for Scientific Research (Poland).
The Matching pursuit algorithm (MP, Mallat and Zhang, 1993b) opened the way for the use of adaptive time/frequency approximations. This approach possess several features that are desirable in the field of EEG analysis, like the easy interpretation of results, a correspondence to the traditional analysis, and the local description of both transient and oscillatory structures with high time/frequency resolution. These advantages have been presented in the context of event-related ˙ ygierewicz et al., 1998), potentials (Durka et al., 2001b; Z ˙ ygierewicz et al., 1999) and electrosleep spindles (Z graphic seizure patterns (Franaszczuk et al., 1998). This study now presents the first application of adaptive time/frequency parametrization to a standard problem of clinical EEG analysis */the quantification of the influence of drugs on sleep EEG.
2. Methods 2.1. Matching pursuit Given a set of functions (dictionary) D /{g1, g2, . . ., gn } such that jjgi jj /1, we can define an optimal M approximation as an expansion, minimizing the error o
0165-0270/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 0 2 7 0 ( 0 2 ) 0 0 0 7 5 - 4
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
66
of an approximation of signal f (t ) by M waveforms:
k
o f (t)
M X
wi ggi (t)
i1
k
(1)
where {gi }i1. . .M represents the indices of the chosen functions ggi . Finding such an optimal approximation is an NP-hard problem (Davis et al., 1994). A sub-optimal expansion can be found by means of an iterative procedure, such as the matching pursuit algorithm (MP, Mallat and Zhang, 1993b). In the first step of MP, the waveform gg0, which best matches the signal f(t), is chosen. In each of the consecutive steps, the waveform ggn is matched to the signal Rn f, which is the residual left after subtracting results of previous iterations: 8 0
where arg maxgg / /D means the ggi giving the largest2 i value of the product jRn f, ggi j. Orthogonality of Rn1f and ggn in each step implies energy conservation: 2
kf k
m1 X
2
jRn f ; ggn j kRm f k
2
(3)
n0
For a complete dictionary the procedure converges to f: f
X
Rn f ; ggn ggn
(4)
n0
Gabor functions (sine-modulated Gaussians) provide optimal joint time /frequency localization. A real Gabor function can be expressed as: p
gg (t) K(g)e
tu s
2
v sin 2p (tu)f N
(5)
N is the size of the signal for which the dictionary is constructed, K (g ) is such that jjgg jj /1. g/{u , v , s , f} denotes parameters of the dictionary’s functions (time / frequency atoms). The length of signal (N points) suggests ranges for the parameters of Gabor functions, which can be reasonably fitted to such an epoch. However, within these ranges no particular sampling is a priori defined, and we are confronted with a threedimensional3 continuous space; this results in an infinite dictionary size. Therefore, in practical implementations, we use subsets of the possible dictionary functions. 2 Or nearly the largest, as explained in detail in Mallat and Zhang (1993b). 3 The phase f is usually optimized in each step of the MP procedure separately for each fitted atom.
In the dictionary implemented originally by Mallat and Zhang (1993a), the parameters {u , v , s } of dictionary’s atoms are chosen from dyadic sequences of integers. This, as any other fixed scheme of subsampling the parameter space, introduces statistical bias in the resulting parametrization. This bias becomes apparent only when statistics comes into play, like in parametrization of large amounts of data */thousands of structures in hours of sleep EEG are a good example. Illustrations of this bias are presented in Durka et al. (2001a). In the same paper, we proposed a solution in terms of MP with stochastic dictionaries, where the parameters of a dictionary’s atoms are randomized before each decomposition, or drawn from flat distributions. Results presented hereby were obtained with such a bias-free implementation.
2.2. Experimental data The effects of zolpidem (10.0 mg), midazolam (7.5 mg) and placebo, administered at bed time, were studied in seven paid volunteers, three males and four females (mean age: 33.6, S.D.9/5.4, range 25/41 years) with no psychiatric, medical or sleep disorders. The volunteers were recruited from hospital staff and university students. All participants were required to be free of any prescription or nonprescription drugs for at least 2 weeks before the study. Sleep was recorded during three blocks of three consecutive nights, with at least a 2 week interval. The first night was an adaptation night. Zolpidem, midazolam or placebo was administered before the second night, in random order, according to a double-blind, crossover design. Data obtained during the third night is not presented in the current study. Polygraphic monitoring consisted of EEG activity collected from 21 electrodes, according to the international 10 /20 system, electrooculogram (EOG) from two channels, and a submentalis electromyogram. Filters were set between 0.15 and 30.0 Hz. The impedance at each electrode was below 5000 V. Sleep stages and artifacts were identified visually (sleep stages according to the standardized manual (Rechtschaffen and Kales, 1968)). Informed consent was obtained from all subjects. The study was approved by the University Ethics Committee. From each of the 21 all-night recordings (three nights for each of seven subjects), we extracted artifact-free epochs from derivation C3 /A2 of the 10 /20 system. For the analysis of slow wave activity (SWA), we used data from stages III and IV, for the u band from all sleep stages, and for sleep spindles only data from sleep stage II.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
67
Table 1 Criteria chosen for the automatic selection of EEG structures based upon the time /frequency parameters
d (SWA) u s (sleep spindles)
Frequency (Hz)
Time duration (s)
Minimum amplitude (mV)
0.75 /4 4 /8 12 /14
0.5 / / 0.5 /2.5
75 30 15
Table 2 Values of spectral estimates of power in mV2 for subjects 1 /7 #
1 2 3 4 5 6 7
d (SWA, 0.75 /4 Hz)
u (4 /8 Hz)
s (12 /14 Hz)
Placebo
Zolpidem
Midazolam
Placebo
Zolpidem
Midazolam
Placebo
Zolpidem
Midazolam
960.33 652.37 911.54 897.11 847.35 383.77 324.43
637.61 618.80 719.53 615.04 816.65 352.96 291.56
601.03 548.15 850.67 675.22 633.66 422.33 218.77
30.38 35.84 42.87 30.86 44.92 32.37 18.51
29.32 32.25 33.87 21.73 35.57 31.39 14.97
25.64 26.70 29.76 21.79 31.86 33.82 14.29
4.07 8.41 6.83 5.12 10.66 8.53 2.42
6.14 10.02 7.38 6.49 11.64 10.40 2.85
5.69 11.38 9.07 9.58 11.94 12.17 2.33
2.3. Criteria of choice of relevant structures The distinction between SWA (d ), u activity and sleep spindles (s) is traditionally based upon division of the spectrum into fixed bands. We adopted the following definitions: d , 0.75 /4 Hz; u , 4 /8 Hz; s, 12 /14 Hz. These choices are complete and sufficient for the traditional pharmaco-EEG, relying on estimating the spectral power for each band. However, consulting a classical definition of sleep spindles4 from Rechtschaffen and Kales (1968) reveals that frequency is not the only relevant parameter. To tune the automatic selection of a relevant signal’s structures to the generally accepted criteria used in visual analysis, we can also use time width (in this study the width of Gaussian envelope in half-height, or parameter s from Eq. (5)) and amplitude. In the absence of clear-cut criteria as in the definition of sleep spindles, the choice of parameters ranges for defining relevant structures remains an open problem. For this study, we chose to adopt structures for d and u bands with minimum restrictions on time widths, namely above half a second for SWA (to contain a full period for lower frequencies) and no restrictions for structures in the u band. The minimum amplitude of a structure is a separate and complicated issue. For example, it can be related e.g. to the notion of structure’s ‘visibility’ (signal-tonoise ratio) as pertaining to the visual analysis, or 4 ‘The presence of sleep spindle should not be defined unless it is of at least 0.5 s duration, i.e. one should be able to count 6 or 7 distinct waves within the half-second period (. . .). The term should be used only to describe activity between 12 and 14 cps’.
assessed separately for each recording. However, such adaptive criteria are not generally accepted. Therefore, to concentrate this study on its main issue, without introducing more free parameters and methods on which results might depend, we chose fixed thresholds for a structure’s minimum amplitude according to the available literature and our experience.
3. Calculations For the calculation of the spectral estimates of energy, artifact-free data from relevant sleep stages was further divided into 4-s epochs. For each of these epochs, the spectral density of energy was estimated as the squared magnitude of discrete Fourier transform of the epoch multiplied by a 4 s Hamming window (MATLAB function psd(×/, 512, 128)). From this estimate, a spectral integral in the relevant range (Table 1) was extracted separately for each epoch. The mean values of these integrals are presented in Table 2. Fig. 1 presents shapes of the spectra in the s band. Time /frequency parametrization was performed on 20-s epochs selected as above from artifact-free EEG in relevant sleep stages. These epochs were subjected to MP decomposition in dictionaries containing 106 Gabor functions and complete Dirac and Fourier bases, as described in the Section 2.1. Waveforms corresponding to relevant EEG structures were chosen according to the criteria given in Table 1. Table 3 gives estimates of the power of selected structures, the number of their occurrences per min and their average frequency. Fig. 2 presents each of the
68
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
Fig. 1. Spectra in the s band (sleep spindles, stage II) for the 21 recordings analyzed in this study. Columns correspond, from the left, to nights after placebo, zolpidem and midazolam. Each row represents one subject (code above each box), identified in the tables with numerical results by numbers 1 /7 increasing from top to bottom.
Fig. 2. Each of the black dots represents one structure classified as a sleep spindle according to criteria from Table 1 in the frequency (horizontal, Hz) vs. amplitude (vertical, mV) coordinates. Each box corresponds to one overnight recording as in Fig. 1.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
69
Table 3 Estimates of power (mV2), occurrences per minute and average (weighted by amplitudes) frequency (Hz), based upon selective time /frequency parametrization of relevant structures #
1
2
3
4
5
6
7
d (SWA, 0.75 /4 Hz)
Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz)
u (4 /8 Hz)
s (12 /14 Hz)
Plac.
Zolp.
Midaz.
Plac.
Zolp.
Midaz.
Plac.
Zolp.
Midaz.
439.15 9.8 1.17 180.06 6.6 1.19 416.82 11.7 1.28 394.90 10.0 1.26 373.02 9.8 1.28 106.92 4.8 1.19 84.06 3.5 1.11
214.20 6.5 1.17 136.53 5.2 1.12 290.39 9.2 1.24 168.65 5.9 1.18 294.82 8.3 1.25 65.22 3.7 1.14 61.29 2.6 1.02
164.32 4.9 1.12 126.54 4.8 1.10 351.54 10.8 1.25 220.39 6.4 1.18 209.19 6.1 1.14 92.00 4.4 1.16 31.77 1.4 0.98
6.74 6.9 5.40 6.93 8.7 5.44 11.98 11.4 5.63 7.11 7.6 5.58 13.52 12.2 5.62 6.53 8.2 5.68 2.64 3.5 5.41
5.60 6.3 5.35 5.24 7.2 5.58 8.11 8.5 5.64 3.59 4.6 5.53 8.05 9.0 5.55 5.72 7.8 5.56 1.43 2.2 5.15
5.08 5.6 5.29 3.98 6.0 5.58 6.92 7.8 5.67 3.24 4.3 5.63 6.96 8.0 5.56 7.47 9.2 5.64 1.19 1.9 5.09
0.81 1.5 13.36 2.86 4.2 13.12 1.63 2.7 13.22 1.19 2.0 13.41 3.72 4.6 13.05 1.92 3.1 13.03 0.16 0.4 13.60
1.75 2.8 13.46 3.92 5.3 13.23 1.95 3.1 13.23 2.32 3.6 13.37 5.23 5.8 13.02 2.82 4.2 13.16 0.41 1.0 13.48
1.91 3.0 13.36 4.46 5.6 13.25 2.85 4.0 13.21 4.06 5.0 13.28 5.60 5.9 12.83 3.45 5.0 13.12 0.27 0.6 13.53
Each column gives these parameters for nights after placebo (plac.), zolpidem (zolp.) or midazolam (midaz.) for each of seven subjects.
structures, classified according to the Table 1 as sleep spindles (s), in frequency /amplitude coordinates for each of the analyzed recordings.
4. Results and discussion The first issue to be considered in the evaluation of a newly proposed tool is a concordance of the results with
those of previously established methods. This is addressed in Table 4, which gives percentage changes of power in recordings after the administration of zolpidem or midazolam (as related to nights after placebo). We observed that both the spectral power in relevant frequency bands, as well as power estimated for selected structures, exhibit the same general trends: an increase of power for sleep spindles and a decrease for SWA and u activities. This is in line with previous studies
Table 4 Percentage changes in spectral (upper, ‘FT’) and selective (lower row of each cell, ‘MP’) estimates of energy after zolpidem (zolp.) or midazolam (midaz.) #
1 2 3 4 5 6 7
d (0.75 /4 Hz)
FT MP FT MP FT MP FT MP FT MP FT MP FT MP
u (4 /8 Hz)
s (12 /14 Hz)
Zolp.
Midaz.
Zolp.
Midaz.
Zolp.
Midaz.
34 51 5 24 21 30 31 57 4 21 8 39 10 27
37 63 16 30 7 16 25 44 25 44 10 14 33 62
3 16 10 17 21 30 30 47 21 40 3 10 19 41
16 23 26 37 31 39 29 48 29 46 4 14 23 52
51 115 19 37 8 20 27 95 9 41 22 47 18 154
40 134 35 56 33 76 87 242 12 51 43 80 4 68
Percentages related to the night after administering placebo.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
70
Table 5 Values of P from Friedmann’s nonparametric two-way ANOVA for no differences between the nights after placebo, zolpidem and midazolam (F) and Wilcoxon signed rank test for nights after placebo and zolpidem (pla-zol) and placebo and midazolam (pla-mid), presented for spectral estimates of energy (FT) and MP-based power (power), number of occurrences per minute (n /min) and frequency (frequency) d (SWA, 0.75 /4 Hz)
FT Power n /min Frequency
u (4 /8 Hz)
s (12 /14 Hz)
F
Pla-zol
Pla-mid
F
Pla-zol
Pla-mid
F
Pla-zol
Pla-mid
0.021 0.005 0.005 0.006
0.031
0.031
0.018 0.012 0.012 0.45
0.3
0.031 0.031 0.031 0.58
0.018 0.021 0.002 0.46
0.69
0.031 0.56
‘’ marks the value of P 0.0156, minimum in Wilcoxon test for paired vectors of size 7.
reporting a decrease of spectral power for frequencies below 10 Hz and an increase in the spindle band (12 /14 Hz). Such behavior is characteristic of benzodiazepines and similar for different benzodiazepine receptor agonists (Brunner et al., 1991; Feige et al., 1999; Trachsel et al., 1990). It should be noted that a concordance with spectral estimates is insufficient to justify the introduction of a new method. However, we observed that selective MPbased estimates of power are more sensitive to the influence of both drugs than spectral integrals. In fact, percentage changes are about twice as high. This sensitivity becomes crucial in two of the analyzed cases: 1) In a recording of patient #7, who displayed a very low spindling activity (cf. last row in Figs. 1 and 2), fluctuations in the background EEG masked the power of sleep spindles to such an extent, that spectral integrals indicated a partially reversed trend, i.e. a decrease in the power of sleep spindles under the influence of midazolam. However, selective estimates indicated the expected increase. 2) For patient #6, spectral integrals indicated an increase in the power of SWA (d ), while MP estimates revealed a decrease consistent with all the other cases. These two cases provide examples of situations where the classical analysis has too low a sensitivity and selectivity to indicate changes clearly revealed by the proposed approach. Application of selective estimates in the u band does not bring such a significant improvement, since in this case we used little information specific to the structures of interest (Table 1). Finally, observed differences were tested statistically, treating averages for each subject under given condition (placebo, zolpidem and midazolam) as three matched vectors with seven values (seven subjects) each using: 1) Friedmann’s nonparametric two-way analysis of variance (ANOVA; MATLAB software package documentation, 2000) for all three vectors, and
2) Wilcoxon signed rank test (MATLAB software package documentation, 2000) for differences between pairs of vectors. Basically no differences were found between the nights after zolpidem and midazolam. Table 5 gives values of probability for the null hypothesis of no differences between all three conditions from Friedman’s ANOVA, and Wilcoxon tests between placebo/ zolpidem and placebo/midazolam. The first two rows of Table 5 confirm what is easily seen in Table 4: the power of structures selected in respective bands reveals changes consistent with spectral integrals, but usually with better significance. The third row indicates that these changes are also very well expressed by the number of relevant structures counted per min. Finally, the last row reveals a new effect: a significant change of the average frequency (weighted by amplitudes) of the SWA. From Table 3, we read that this change is a drug-induced decrease of frequency.
5. Conclusion The proposed approach links the advantages of visual and spectral analysis of the EEG with advanced time / frequency signal processing in one unified formalism, that is suitable for clinical applications. It offers a detailed high-resolution description of signal structures, that is compatible with classical definitions of some of the EEG transients. This allows the investigation of the evolution of properties of selected structures with a specifity and resolution unmatched by any other currently available method. Although not limited to pharmaco-EEG, this methodology improves the sensitivity and objectivity of the evaluation of the influence of drugs and opens new possibilities for investigating their effects on the EEG. As an example of an effect elusive to classical methodology, we illustrate a decrease in frequency of SWA under the influene of sleepinducing drugs.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
6. Reproducible research Software for matching pursuit with stochastic dictionaries is available from http://brain.fuw.edu.pl/mp.
Acknowledgements We gratefully acknowledge the assistance of our colleagues Dobieslaw Ircha and Dominik Blacha who wrote the software for MP decomposition used in this research. The EEG was recorded on equipment donated by the AJUS&KAJUS Foundation, in memory of Professor Andrzej Jus, the pioneer of Polish Clinical EEG and the first person to introduce polygraphic studies of sleep in Poland. This work was partially supported by the grant of Committee for Scientific Research (Poland) to the Institute of Experimental Physics, Warsaw University.
References Brunner DP, Dijk DJ, Mu¨nch M, Borbe´ly A. Effect of zolpidem on sleep and sleep EEG spectra in healthy young men. Psychopharmacology 1991;104:1 /5. Cooley JW, Tukey JW. An algorithm for the machine calculation of complex fourier series. Math Comput 1965;19:297 /301. Davis G. Adaptive nonlinear approximations. Ph.D. Thesis, New York University, 1994; ftp://cs.nyu.edu/pub/wave/report/DissertationGDavis.ps.Z
71
Dietsch G. Fourier analyse von elektroenzephalogrammen des menschen. Pflugers Arch Ges Physiol 1932;230:106 /12. Durka PJ, Ircha D, Blinowska KJ. Stochastic time /frequency dictionaries for matching pursuit. IEEE Trans Signal Process 2001;49(3):507 /10. Durka PJ, Ircha D, Neuper C, Pfurtscheller G. Time /frequency microstructure of event-related desynchronization and synchronization. Med Biol Eng Comput 2001;39(3):315 /21. Feige B, Voderholzer U, Riemann D, Hohagen F, Berger M. Independent sleep EEG slow-wave and spindle band dynamics associated with 4 weeks of continuous application of short half-life hypnotics in healthy subjects. Clin Neurophysiol 1999;110:1965 / 74. Franaszczuk PJ, Bergey G, Durka PJ, Eisenberg H. Time /frequency analysis using the matching pursuit algorithm applied to seizures originating from the mesial temporal lobe. Electroencephalogr Clin Neurophysiol 1998;106:513 /21. Mallat S, Zhang Z. Matching pursuit software package (mpp), 1993; Available online at ftp://cs.nyu.edu/pub/wave/software/mpp.tar-Z Mallat S, Zhang Z. Matching pursuit with time /frequency dictionaries. IEEE Trans Signal Process 1993;41:3397 /415. MATLAB software package documentation. Statistic Toolbox User’s Guide (Version 3). The MathWorks Inc. 2000. Rechtschaffen A, Kales A. A manual of standardized terminology, techniques and scoring system for sleep stages in human subjects. In: Number 204 in National Institutes of Health Publications. Washington, DC: US Government Printing Office, 1968. Trachsel L, Dijk DJ, Brunner DP, Klene C, Borbely AA. Effect of zopiclone and midazolam on sleep and EEC spectra in a phaseadvanced sleep schedule. Neuropsychopharmacology 1990;3:11 /8. ˙ ygierewicz J, Kelly EF, Blinowska KJ, Durka PJ, Folger S. Time / .Z frequency analysis of vibrotactile driving responses by matching pursuit. J Neurosci Methods 1998;81:121 /9. ˙ ygierewicz J, Blinowska KJ, Durka PJ, Szelenberger W, Niemcewicz .Z S, Androsiuk W. High resolution study of sleep spindles. Clin Neurophysiol 1999;110(12):2136 /47.
Adaptive time frequency parametrization in pharmaco EEG /
Piotr J. Durka a,*,1, Waldemar Szelenberger b, Katarzyna J. Blinowska a, Wojciech Androsiuk b, Maciej Myszka b a
Laboratory of Medical Physics, Institute of Experimental Physics, Warsaw University, ul. Hoz˙a 69, 00-681 Warszawa, Poland b Department of Psychiatry, Medical University of Warsaw, ul. Nowowiejska 27, 00-665 Warszawa, Poland Received 18 January 2002; received in revised form 15 March 2002; accepted 18 March 2002
Abstract Adaptive time /frequency approximations offer description of the local structures of a signal in terms of their time and frequency coordinates, widths and amplitudes. These parameters can then be used to select and study electroencephalogram (EEG) structures like sleep spindles or slow wave activity (SWA) with high resolution. Such a detailed description of relevant structures improves on the sensitivity of the traditionally used spectral power estimates and opens new possibilities of investigation. These advantages are illustrated using a double-blind test of the influence of zolpidem and midazolam on sleep EEG, and the results are compared with the traditional approach. The observed decrease of frequency of the SWA under the influence of sleep-inducing drugs gives an example of an effect elusive to classical methodology. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Matching pursuit; Time /frequency parametrization; Pharmaco EEG; Sleep spindles; Slow wave activity; Zolpidem; Midazolam
1. Introduction The estimation of spectral power in fixed frequency bands by means of the Fourier transform is one of the very few signal processing methods that is widely adopted in the practice of analyzing the electroencephalogram (EEG). Its general acceptance is due to a long record of applications (since the first application of Fourier transform to EEG (Dietsch, 1932)), fast numerical implementations (since the introduction of the fast Fourier transform by Cooley and Tukey (1965)) and the relatively simple interpretation of results in terms of oscillations and cycles per second. In contrast, none of the advanced signal processing methods proposed in the last few decades have gained any comparable acceptance.
* Corresponding author. Tel.: /48-22-553-2126; fax: /48-22-6226154; http://brain.fuw.edu.pl. E-mail addresses: [email protected] (P.J. Durka), [email protected] (W. Szelenberger), [email protected] (K.J. Blinowska), [email protected] (W. Androsiuk), maciek@ psych.waw.pl (M. Myszka). 1 Partially supported by Committee for Scientific Research (Poland).
The Matching pursuit algorithm (MP, Mallat and Zhang, 1993b) opened the way for the use of adaptive time/frequency approximations. This approach possess several features that are desirable in the field of EEG analysis, like the easy interpretation of results, a correspondence to the traditional analysis, and the local description of both transient and oscillatory structures with high time/frequency resolution. These advantages have been presented in the context of event-related ˙ ygierewicz et al., 1998), potentials (Durka et al., 2001b; Z ˙ ygierewicz et al., 1999) and electrosleep spindles (Z graphic seizure patterns (Franaszczuk et al., 1998). This study now presents the first application of adaptive time/frequency parametrization to a standard problem of clinical EEG analysis */the quantification of the influence of drugs on sleep EEG.
2. Methods 2.1. Matching pursuit Given a set of functions (dictionary) D /{g1, g2, . . ., gn } such that jjgi jj /1, we can define an optimal M approximation as an expansion, minimizing the error o
0165-0270/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 0 2 7 0 ( 0 2 ) 0 0 0 7 5 - 4
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
66
of an approximation of signal f (t ) by M waveforms:
k
o f (t)
M X
wi ggi (t)
i1
k
(1)
where {gi }i1. . .M represents the indices of the chosen functions ggi . Finding such an optimal approximation is an NP-hard problem (Davis et al., 1994). A sub-optimal expansion can be found by means of an iterative procedure, such as the matching pursuit algorithm (MP, Mallat and Zhang, 1993b). In the first step of MP, the waveform gg0, which best matches the signal f(t), is chosen. In each of the consecutive steps, the waveform ggn is matched to the signal Rn f, which is the residual left after subtracting results of previous iterations: 8 0
where arg maxgg / /D means the ggi giving the largest2 i value of the product jRn f, ggi j. Orthogonality of Rn1f and ggn in each step implies energy conservation: 2
kf k
m1 X
2
jRn f ; ggn j kRm f k
2
(3)
n0
For a complete dictionary the procedure converges to f: f
X
Rn f ; ggn ggn
(4)
n0
Gabor functions (sine-modulated Gaussians) provide optimal joint time /frequency localization. A real Gabor function can be expressed as: p
gg (t) K(g)e
tu s
2
v sin 2p (tu)f N
(5)
N is the size of the signal for which the dictionary is constructed, K (g ) is such that jjgg jj /1. g/{u , v , s , f} denotes parameters of the dictionary’s functions (time / frequency atoms). The length of signal (N points) suggests ranges for the parameters of Gabor functions, which can be reasonably fitted to such an epoch. However, within these ranges no particular sampling is a priori defined, and we are confronted with a threedimensional3 continuous space; this results in an infinite dictionary size. Therefore, in practical implementations, we use subsets of the possible dictionary functions. 2 Or nearly the largest, as explained in detail in Mallat and Zhang (1993b). 3 The phase f is usually optimized in each step of the MP procedure separately for each fitted atom.
In the dictionary implemented originally by Mallat and Zhang (1993a), the parameters {u , v , s } of dictionary’s atoms are chosen from dyadic sequences of integers. This, as any other fixed scheme of subsampling the parameter space, introduces statistical bias in the resulting parametrization. This bias becomes apparent only when statistics comes into play, like in parametrization of large amounts of data */thousands of structures in hours of sleep EEG are a good example. Illustrations of this bias are presented in Durka et al. (2001a). In the same paper, we proposed a solution in terms of MP with stochastic dictionaries, where the parameters of a dictionary’s atoms are randomized before each decomposition, or drawn from flat distributions. Results presented hereby were obtained with such a bias-free implementation.
2.2. Experimental data The effects of zolpidem (10.0 mg), midazolam (7.5 mg) and placebo, administered at bed time, were studied in seven paid volunteers, three males and four females (mean age: 33.6, S.D.9/5.4, range 25/41 years) with no psychiatric, medical or sleep disorders. The volunteers were recruited from hospital staff and university students. All participants were required to be free of any prescription or nonprescription drugs for at least 2 weeks before the study. Sleep was recorded during three blocks of three consecutive nights, with at least a 2 week interval. The first night was an adaptation night. Zolpidem, midazolam or placebo was administered before the second night, in random order, according to a double-blind, crossover design. Data obtained during the third night is not presented in the current study. Polygraphic monitoring consisted of EEG activity collected from 21 electrodes, according to the international 10 /20 system, electrooculogram (EOG) from two channels, and a submentalis electromyogram. Filters were set between 0.15 and 30.0 Hz. The impedance at each electrode was below 5000 V. Sleep stages and artifacts were identified visually (sleep stages according to the standardized manual (Rechtschaffen and Kales, 1968)). Informed consent was obtained from all subjects. The study was approved by the University Ethics Committee. From each of the 21 all-night recordings (three nights for each of seven subjects), we extracted artifact-free epochs from derivation C3 /A2 of the 10 /20 system. For the analysis of slow wave activity (SWA), we used data from stages III and IV, for the u band from all sleep stages, and for sleep spindles only data from sleep stage II.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
67
Table 1 Criteria chosen for the automatic selection of EEG structures based upon the time /frequency parameters
d (SWA) u s (sleep spindles)
Frequency (Hz)
Time duration (s)
Minimum amplitude (mV)
0.75 /4 4 /8 12 /14
0.5 / / 0.5 /2.5
75 30 15
Table 2 Values of spectral estimates of power in mV2 for subjects 1 /7 #
1 2 3 4 5 6 7
d (SWA, 0.75 /4 Hz)
u (4 /8 Hz)
s (12 /14 Hz)
Placebo
Zolpidem
Midazolam
Placebo
Zolpidem
Midazolam
Placebo
Zolpidem
Midazolam
960.33 652.37 911.54 897.11 847.35 383.77 324.43
637.61 618.80 719.53 615.04 816.65 352.96 291.56
601.03 548.15 850.67 675.22 633.66 422.33 218.77
30.38 35.84 42.87 30.86 44.92 32.37 18.51
29.32 32.25 33.87 21.73 35.57 31.39 14.97
25.64 26.70 29.76 21.79 31.86 33.82 14.29
4.07 8.41 6.83 5.12 10.66 8.53 2.42
6.14 10.02 7.38 6.49 11.64 10.40 2.85
5.69 11.38 9.07 9.58 11.94 12.17 2.33
2.3. Criteria of choice of relevant structures The distinction between SWA (d ), u activity and sleep spindles (s) is traditionally based upon division of the spectrum into fixed bands. We adopted the following definitions: d , 0.75 /4 Hz; u , 4 /8 Hz; s, 12 /14 Hz. These choices are complete and sufficient for the traditional pharmaco-EEG, relying on estimating the spectral power for each band. However, consulting a classical definition of sleep spindles4 from Rechtschaffen and Kales (1968) reveals that frequency is not the only relevant parameter. To tune the automatic selection of a relevant signal’s structures to the generally accepted criteria used in visual analysis, we can also use time width (in this study the width of Gaussian envelope in half-height, or parameter s from Eq. (5)) and amplitude. In the absence of clear-cut criteria as in the definition of sleep spindles, the choice of parameters ranges for defining relevant structures remains an open problem. For this study, we chose to adopt structures for d and u bands with minimum restrictions on time widths, namely above half a second for SWA (to contain a full period for lower frequencies) and no restrictions for structures in the u band. The minimum amplitude of a structure is a separate and complicated issue. For example, it can be related e.g. to the notion of structure’s ‘visibility’ (signal-tonoise ratio) as pertaining to the visual analysis, or 4 ‘The presence of sleep spindle should not be defined unless it is of at least 0.5 s duration, i.e. one should be able to count 6 or 7 distinct waves within the half-second period (. . .). The term should be used only to describe activity between 12 and 14 cps’.
assessed separately for each recording. However, such adaptive criteria are not generally accepted. Therefore, to concentrate this study on its main issue, without introducing more free parameters and methods on which results might depend, we chose fixed thresholds for a structure’s minimum amplitude according to the available literature and our experience.
3. Calculations For the calculation of the spectral estimates of energy, artifact-free data from relevant sleep stages was further divided into 4-s epochs. For each of these epochs, the spectral density of energy was estimated as the squared magnitude of discrete Fourier transform of the epoch multiplied by a 4 s Hamming window (MATLAB function psd(×/, 512, 128)). From this estimate, a spectral integral in the relevant range (Table 1) was extracted separately for each epoch. The mean values of these integrals are presented in Table 2. Fig. 1 presents shapes of the spectra in the s band. Time /frequency parametrization was performed on 20-s epochs selected as above from artifact-free EEG in relevant sleep stages. These epochs were subjected to MP decomposition in dictionaries containing 106 Gabor functions and complete Dirac and Fourier bases, as described in the Section 2.1. Waveforms corresponding to relevant EEG structures were chosen according to the criteria given in Table 1. Table 3 gives estimates of the power of selected structures, the number of their occurrences per min and their average frequency. Fig. 2 presents each of the
68
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
Fig. 1. Spectra in the s band (sleep spindles, stage II) for the 21 recordings analyzed in this study. Columns correspond, from the left, to nights after placebo, zolpidem and midazolam. Each row represents one subject (code above each box), identified in the tables with numerical results by numbers 1 /7 increasing from top to bottom.
Fig. 2. Each of the black dots represents one structure classified as a sleep spindle according to criteria from Table 1 in the frequency (horizontal, Hz) vs. amplitude (vertical, mV) coordinates. Each box corresponds to one overnight recording as in Fig. 1.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
69
Table 3 Estimates of power (mV2), occurrences per minute and average (weighted by amplitudes) frequency (Hz), based upon selective time /frequency parametrization of relevant structures #
1
2
3
4
5
6
7
d (SWA, 0.75 /4 Hz)
Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz) Power (mV2) Occurrences (/min) Frequency (Hz)
u (4 /8 Hz)
s (12 /14 Hz)
Plac.
Zolp.
Midaz.
Plac.
Zolp.
Midaz.
Plac.
Zolp.
Midaz.
439.15 9.8 1.17 180.06 6.6 1.19 416.82 11.7 1.28 394.90 10.0 1.26 373.02 9.8 1.28 106.92 4.8 1.19 84.06 3.5 1.11
214.20 6.5 1.17 136.53 5.2 1.12 290.39 9.2 1.24 168.65 5.9 1.18 294.82 8.3 1.25 65.22 3.7 1.14 61.29 2.6 1.02
164.32 4.9 1.12 126.54 4.8 1.10 351.54 10.8 1.25 220.39 6.4 1.18 209.19 6.1 1.14 92.00 4.4 1.16 31.77 1.4 0.98
6.74 6.9 5.40 6.93 8.7 5.44 11.98 11.4 5.63 7.11 7.6 5.58 13.52 12.2 5.62 6.53 8.2 5.68 2.64 3.5 5.41
5.60 6.3 5.35 5.24 7.2 5.58 8.11 8.5 5.64 3.59 4.6 5.53 8.05 9.0 5.55 5.72 7.8 5.56 1.43 2.2 5.15
5.08 5.6 5.29 3.98 6.0 5.58 6.92 7.8 5.67 3.24 4.3 5.63 6.96 8.0 5.56 7.47 9.2 5.64 1.19 1.9 5.09
0.81 1.5 13.36 2.86 4.2 13.12 1.63 2.7 13.22 1.19 2.0 13.41 3.72 4.6 13.05 1.92 3.1 13.03 0.16 0.4 13.60
1.75 2.8 13.46 3.92 5.3 13.23 1.95 3.1 13.23 2.32 3.6 13.37 5.23 5.8 13.02 2.82 4.2 13.16 0.41 1.0 13.48
1.91 3.0 13.36 4.46 5.6 13.25 2.85 4.0 13.21 4.06 5.0 13.28 5.60 5.9 12.83 3.45 5.0 13.12 0.27 0.6 13.53
Each column gives these parameters for nights after placebo (plac.), zolpidem (zolp.) or midazolam (midaz.) for each of seven subjects.
structures, classified according to the Table 1 as sleep spindles (s), in frequency /amplitude coordinates for each of the analyzed recordings.
4. Results and discussion The first issue to be considered in the evaluation of a newly proposed tool is a concordance of the results with
those of previously established methods. This is addressed in Table 4, which gives percentage changes of power in recordings after the administration of zolpidem or midazolam (as related to nights after placebo). We observed that both the spectral power in relevant frequency bands, as well as power estimated for selected structures, exhibit the same general trends: an increase of power for sleep spindles and a decrease for SWA and u activities. This is in line with previous studies
Table 4 Percentage changes in spectral (upper, ‘FT’) and selective (lower row of each cell, ‘MP’) estimates of energy after zolpidem (zolp.) or midazolam (midaz.) #
1 2 3 4 5 6 7
d (0.75 /4 Hz)
FT MP FT MP FT MP FT MP FT MP FT MP FT MP
u (4 /8 Hz)
s (12 /14 Hz)
Zolp.
Midaz.
Zolp.
Midaz.
Zolp.
Midaz.
34 51 5 24 21 30 31 57 4 21 8 39 10 27
37 63 16 30 7 16 25 44 25 44 10 14 33 62
3 16 10 17 21 30 30 47 21 40 3 10 19 41
16 23 26 37 31 39 29 48 29 46 4 14 23 52
51 115 19 37 8 20 27 95 9 41 22 47 18 154
40 134 35 56 33 76 87 242 12 51 43 80 4 68
Percentages related to the night after administering placebo.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
70
Table 5 Values of P from Friedmann’s nonparametric two-way ANOVA for no differences between the nights after placebo, zolpidem and midazolam (F) and Wilcoxon signed rank test for nights after placebo and zolpidem (pla-zol) and placebo and midazolam (pla-mid), presented for spectral estimates of energy (FT) and MP-based power (power), number of occurrences per minute (n /min) and frequency (frequency) d (SWA, 0.75 /4 Hz)
FT Power n /min Frequency
u (4 /8 Hz)
s (12 /14 Hz)
F
Pla-zol
Pla-mid
F
Pla-zol
Pla-mid
F
Pla-zol
Pla-mid
0.021 0.005 0.005 0.006
0.031
0.031
0.018 0.012 0.012 0.45
0.3
0.031 0.031 0.031 0.58
0.018 0.021 0.002 0.46
0.69
0.031 0.56
‘’ marks the value of P 0.0156, minimum in Wilcoxon test for paired vectors of size 7.
reporting a decrease of spectral power for frequencies below 10 Hz and an increase in the spindle band (12 /14 Hz). Such behavior is characteristic of benzodiazepines and similar for different benzodiazepine receptor agonists (Brunner et al., 1991; Feige et al., 1999; Trachsel et al., 1990). It should be noted that a concordance with spectral estimates is insufficient to justify the introduction of a new method. However, we observed that selective MPbased estimates of power are more sensitive to the influence of both drugs than spectral integrals. In fact, percentage changes are about twice as high. This sensitivity becomes crucial in two of the analyzed cases: 1) In a recording of patient #7, who displayed a very low spindling activity (cf. last row in Figs. 1 and 2), fluctuations in the background EEG masked the power of sleep spindles to such an extent, that spectral integrals indicated a partially reversed trend, i.e. a decrease in the power of sleep spindles under the influence of midazolam. However, selective estimates indicated the expected increase. 2) For patient #6, spectral integrals indicated an increase in the power of SWA (d ), while MP estimates revealed a decrease consistent with all the other cases. These two cases provide examples of situations where the classical analysis has too low a sensitivity and selectivity to indicate changes clearly revealed by the proposed approach. Application of selective estimates in the u band does not bring such a significant improvement, since in this case we used little information specific to the structures of interest (Table 1). Finally, observed differences were tested statistically, treating averages for each subject under given condition (placebo, zolpidem and midazolam) as three matched vectors with seven values (seven subjects) each using: 1) Friedmann’s nonparametric two-way analysis of variance (ANOVA; MATLAB software package documentation, 2000) for all three vectors, and
2) Wilcoxon signed rank test (MATLAB software package documentation, 2000) for differences between pairs of vectors. Basically no differences were found between the nights after zolpidem and midazolam. Table 5 gives values of probability for the null hypothesis of no differences between all three conditions from Friedman’s ANOVA, and Wilcoxon tests between placebo/ zolpidem and placebo/midazolam. The first two rows of Table 5 confirm what is easily seen in Table 4: the power of structures selected in respective bands reveals changes consistent with spectral integrals, but usually with better significance. The third row indicates that these changes are also very well expressed by the number of relevant structures counted per min. Finally, the last row reveals a new effect: a significant change of the average frequency (weighted by amplitudes) of the SWA. From Table 3, we read that this change is a drug-induced decrease of frequency.
5. Conclusion The proposed approach links the advantages of visual and spectral analysis of the EEG with advanced time / frequency signal processing in one unified formalism, that is suitable for clinical applications. It offers a detailed high-resolution description of signal structures, that is compatible with classical definitions of some of the EEG transients. This allows the investigation of the evolution of properties of selected structures with a specifity and resolution unmatched by any other currently available method. Although not limited to pharmaco-EEG, this methodology improves the sensitivity and objectivity of the evaluation of the influence of drugs and opens new possibilities for investigating their effects on the EEG. As an example of an effect elusive to classical methodology, we illustrate a decrease in frequency of SWA under the influene of sleepinducing drugs.
P.J. Durka et al. / Journal of Neuroscience Methods 117 (2002) 65 /71
6. Reproducible research Software for matching pursuit with stochastic dictionaries is available from http://brain.fuw.edu.pl/mp.
Acknowledgements We gratefully acknowledge the assistance of our colleagues Dobieslaw Ircha and Dominik Blacha who wrote the software for MP decomposition used in this research. The EEG was recorded on equipment donated by the AJUS&KAJUS Foundation, in memory of Professor Andrzej Jus, the pioneer of Polish Clinical EEG and the first person to introduce polygraphic studies of sleep in Poland. This work was partially supported by the grant of Committee for Scientific Research (Poland) to the Institute of Experimental Physics, Warsaw University.
References Brunner DP, Dijk DJ, Mu¨nch M, Borbe´ly A. Effect of zolpidem on sleep and sleep EEG spectra in healthy young men. Psychopharmacology 1991;104:1 /5. Cooley JW, Tukey JW. An algorithm for the machine calculation of complex fourier series. Math Comput 1965;19:297 /301. Davis G. Adaptive nonlinear approximations. Ph.D. Thesis, New York University, 1994; ftp://cs.nyu.edu/pub/wave/report/DissertationGDavis.ps.Z
71
Dietsch G. Fourier analyse von elektroenzephalogrammen des menschen. Pflugers Arch Ges Physiol 1932;230:106 /12. Durka PJ, Ircha D, Blinowska KJ. Stochastic time /frequency dictionaries for matching pursuit. IEEE Trans Signal Process 2001;49(3):507 /10. Durka PJ, Ircha D, Neuper C, Pfurtscheller G. Time /frequency microstructure of event-related desynchronization and synchronization. Med Biol Eng Comput 2001;39(3):315 /21. Feige B, Voderholzer U, Riemann D, Hohagen F, Berger M. Independent sleep EEG slow-wave and spindle band dynamics associated with 4 weeks of continuous application of short half-life hypnotics in healthy subjects. Clin Neurophysiol 1999;110:1965 / 74. Franaszczuk PJ, Bergey G, Durka PJ, Eisenberg H. Time /frequency analysis using the matching pursuit algorithm applied to seizures originating from the mesial temporal lobe. Electroencephalogr Clin Neurophysiol 1998;106:513 /21. Mallat S, Zhang Z. Matching pursuit software package (mpp), 1993; Available online at ftp://cs.nyu.edu/pub/wave/software/mpp.tar-Z Mallat S, Zhang Z. Matching pursuit with time /frequency dictionaries. IEEE Trans Signal Process 1993;41:3397 /415. MATLAB software package documentation. Statistic Toolbox User’s Guide (Version 3). The MathWorks Inc. 2000. Rechtschaffen A, Kales A. A manual of standardized terminology, techniques and scoring system for sleep stages in human subjects. In: Number 204 in National Institutes of Health Publications. Washington, DC: US Government Printing Office, 1968. Trachsel L, Dijk DJ, Brunner DP, Klene C, Borbely AA. Effect of zopiclone and midazolam on sleep and EEC spectra in a phaseadvanced sleep schedule. Neuropsychopharmacology 1990;3:11 /8. ˙ ygierewicz J, Kelly EF, Blinowska KJ, Durka PJ, Folger S. Time / .Z frequency analysis of vibrotactile driving responses by matching pursuit. J Neurosci Methods 1998;81:121 /9. ˙ ygierewicz J, Blinowska KJ, Durka PJ, Szelenberger W, Niemcewicz .Z S, Androsiuk W. High resolution study of sleep spindles. Clin Neurophysiol 1999;110(12):2136 /47.