Table Of ContentWavelet analysis of epileptic spikes
Miroslaw Latka∗ and Ziemowit Was†
Institute of Physics, Wroclaw University of Technology,
Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland
Andrzej Kozik‡
Video EEG Lab, Department of Child Neurology,
3
Regional Medical Center, ul. Traugutta 116, 40-529 Wroclaw, Poland ‡
0
0
2 Bruce J. West§
Mathematics Division, Army Research Office, P.O. Box 12211, Research Triangle, NC 27709-2211, USA
n
(Dated: January 17, 2003)
a
J
Interictal spikes and sharp waves in human EEG are characteristic signatures of epilepsy. These
7 potentialsoriginateasaresultofsynchronous,pathologicaldischargeofmanyneurons. Thereliable
2 detection of such potentials has been the long standing problem in EEG analysis, especially after
long-term monitoring becamecommon in investigation of epilepticpatients. Thetraditional defini-
]
tion of a spike is based on its amplitude, duration, sharpness, and emergence from its background.
h
p However,spikedetectionsystemsbuiltsolelyaroundthisdefinitionarenotreliableduetothepres-
- ence of numerous transients and artifacts. We use wavelet transform to analyze the properties of
o EEG manifestations of epilepsy. We demonstrate that the behavior of wavelet transform of epilep-
i tic spikes across scales can constitute the foundation of a relatively simple yet effective detection
b
algorithm.
.
s
c
PACSnumbers: 87.10. +e,05.45.-a
i
s
y
h Recordings of human brain electrical activity (EEG) growingsteadilyandnowincludesthecharacterizationof
p
havebeenthe fundamentaltoolforstudying the dynam- encephalopaties [4], monitoring of anesthesia depth [5],
[
ics of cortical neurons since 1929. Even though the gist characteristics of seizure activity [6], and prediction of
1 of this technique has essentially remained the same, the epileptic seizures [7]. Several other approaches are also
v methods of EEG data analysis have profoundly evolved used to elucidate the nature of electrical activity of the
5 during the last two decades. In 1985 Babloyantz et al. humanbrainrangingfromcoherencemeasures[8, 9]and
6
demonstrated that certain nonlinear measures, first in- methods of nonequilibrium statistical mechanics [10] to
0
troduced in the context of chaotic dynamical systems, complexity measures [11, 12].
1
changed during slow-wave sleep [1]. The flurry of re-
0 One of the most important challenges of EEG analy-
3 search work that followed this discovery focused on the
sis has been the quantification of the manifestations of
0 application of nonlinear dynamics in quantifying brain
epilepsy. The main goal is to establish a correlation be-
/ electrical activity during different mental states, sleep
s tween the EEG and clinical or pharmacological condi-
c stages, and under the influence of the epileptic process
tions. One of the possible approaches is based on the
i (for a review see for example [2, 3]). It must be em-
s properties of the interictal EEG (electrical activity mea-
y phasized that a straightforward interpretation of neu-
sured between seizures) which typically consists of lin-
h ral dynamics in terms of such nonlinear measures as the
ear stochastic backgroundfluctuations interspersed with
p largest Lyapunov exponent or the correlation dimension
: transient nonlinear spikes or sharp waves. These tran-
v is not possible since most biological time series, such as
sient potentials originate as a result of a simultaneous
i EEG, are nonstationary and consequently do not sat-
X pathological discharge of neurons within a volume of at
isfy the assumptions of the underlying theory. On the
least several mm3.
r
other hand, traditional power spectral methods are also
a
based on quite restrictive assumptions but nevertheless The traditional definition of a spike is based on its
have turned out to be successful in some areas of EEG amplitude, duration, sharpness and emergence from its
analysis. Despite these technical difficulties, the number background [13, 14]. However, automatic epileptic spike
of applications of nonlinear time series analysis has been detection systems based on this direct approach suffer
from false detections in the presence of numerous types
of artifacts and non-epileptic transients. This shortcom-
ing is particularly acute for long-term EEG monitoring
∗Electronicaddress: [email protected] ofepileptic patientswhich became commonin 1980s. To
†Electronicaddress: [email protected] reduce false detections Gotman and Wang [15] made the
‡Electronicaddress: [email protected] process of spike identification dependent upon the state
§Electronicaddress: [email protected]
of EEG (active wakefulness, quiet wakefulness, desyn-
2
chronizedEEG,phasicEEG,andslowEEG).Thismodi- where k = 0,...,N 1 is the frequency index. If one
−
ficationleadstosignificantoverallimprovementprovided notes that the Fourier transform of a function ψ(t/a) is
that state classification is correct. aψˆ(af) then by the convolution theorem
Diambra and Malta [16] adopted nonlinear prediction | |
for epileptic spike detection. They demonstrated that N−1
when the model’s parameters are adjusted during the Wn(a)=√aδt hˆnψ∗(afk)e2πifknδt, (5)
“learning”phase to assure good predictive performance kX=0
for stochastic background fluctuations, the appearance
frequencies f are defined in the conventional way. Us-
of an interictal spike is marked by a very large forecast- k
ing(5)andastandardfastFouriertransform(FFT)rou-
ing error. This novel approach is appealing because it
tine it is possible to efficiently calculate the continuous
makes use of changes in EEG dynamics. One expects
wavelet transform (for a given scale a) at all n simul-
goodnonlinearpredictiveperformancewhenthe dynam-
taneously [20]. It should be emphasized that formally
ics ofthe EEGintervalusedfor building up the model is
equation (5) does not yield the discrete linear convolu-
similar to the dynamics of the interval used for testing.
tion corresponding to (3) but rather a discrete circular
However,it is uncertainatthis pointwhether it is possi-
convolution in which the shift n′ n is taken modulo
ble to develop a robust spike detection algorithm based −
N. However, in the context of this work, this problem
solely on this idea.
does not give rise to any numerical difficulties. This is
AsClarket al. putitsuccinctly,automaticEEGanal-
because, for purely practical reasons, the beginning and
ysisisaformidabletaskbecauseofthelackof“...features
theendoftheanalyzedpartofdatastreamarenottaken
that reflect the relevant information ”[17]. Another dif-
into account during the EEG spike detection.
ficulty is the nonstationary nature of the spikes and the
From a plethora of available mother wavelets, we em-
backgroundin which they are embedded. One technique
ploy the Mexican hat
developed for the treatment of such nonstationary time
series is wavelet analysis [18, 19]. The goal of this paper ψ(t)= 2 π−1/4(1 t2)e−t2/2 (6)
is to characterizethe epileptic spikes and sharp wavesin √3 −
terms of the properties of their wavelet transforms. In
particular, we search for features which could be impor- which is particularly suitable for studying epileptic
tant in the detection of epileptic events. events.
The wavelet transform is an integral transform for In the top panel of Fig. 1 we present two pieces of
which the set of basis functions, known as wavelets, are the EEG recording joined at approximately t=1s. The
well localizedboth in time and frequency. Moreover,the digital 19 channel recording sampled at 240 Hz was ob-
wavelet basis can be constructed from a single function tained from a juvenile epileptic patient according to the
ψ(t) by means of translation and dilation: international 10-20 standard with the reference average
electrode. The epileptic spike in this figure (marked by
t t
ψa;t0 =ψ(cid:18) −a 0(cid:19). (1) tohfeFaigrr.o1w)diisspfloalyloswtehdebcyontwtoouarrmtifaapctosf.tThheeabbosottluomtepvaanlueel
ψ(t) is commonly referred to as the mother function or of Mexican hat wavelet coefficients W(a,t0). It is ap-
analyzing wavelet. The wavelet transform of function parent that the red prominent ridges correspond to the
h(t) is defined as position of either spike or the motion artifacts. What
is most important, for small scales, a, the values of the
W(a,t )= 1 ∞ h(t)ψ∗ dt, (2) wavelet coefficients for the spike’s ridge are much larger
0 √aZ a;t0 than those for the artifacts. The peak value along the
−∞
spike ridge corresponds to a = 7. In sharp contrast, for
where ψ∗(t) denotes the complex conjugate of ψ(t). The
the range of scales used in Fig. 1 the absolute value of
continuous wavelet transform of a discrete time series
coefficients W(a,t )for the artifactsgrowmonotonically
h N−1 of length N and equal spacing δt is defined as 0
{ i}i=0 with a.
δt N−1 (n′ n)δt The question arises as to whether the behavior of the
Wn(a)=ra hn′ψ∗(cid:20) −a (cid:21). (3) wavelet transform as a function of scale can be used to
nX′=0 develop a reliable detection algorithm. The first step in
this direction is to use the normalized wavelet power
The above convolution can be evaluated for any of N
values of the time index n. However, by choosing all N
w(a,t )=W2(a,t )/σ2 (7)
successivetimeindexvalues,theconvolutiontheoremal- 0 0
lowsusto calculateallN convolutionssimultaneouslyin
instead of the wavelet coefficients to reduce the depen-
Fourier space using a discrete Fourier transform (DFT).
dence on the amplitude of the EEG recording. In the
The DFT of {hi}Ni=−01 is aboveformulaσ2 isthevarianceoftheportionofthesig-
nalbeinganalyzed(typicallyweusepiecesoflength1024
N−1
1
hˆ = h e−2πikn/N, (4) for EEG tracings sampled at 240 Hz). In actual numeri-
k n
N nX=0 calcalculations we prefer to use the square ofw(a,t0) to
3
merely increase the range of values analyzed during the (8).
spike detection process. In Fig. 2 w2 for the signal used
For testing purposes, we built up the database of arti-
in Fig. 1 is plotted for three scales A = 3, B = 7 and
facts and spikes. We made available some of these EEG
C =20.
tracings at [21] along with the examples of the numeri-
In the most straightforward approach, we identify an
calcalculations. While the analysisofthe pieces ofEEG
EEG transient potential as a simple or isolated epileptic
recordingssuchasthoseshowninFig. 2and3isessential
spike if and only if:
in determining the generic properties of epileptic events,
the value of w2 at a = 7 is greater than a predeter-
it can hardly reflect the difficulties one can encounter
•
mined threshold value T ,
1 in interpretation of clinical EEG. Therefore we selected
the square of normalized wavelet power decreases
four challenging EEG tracings with 340 epileptic events.
•
from scale a=7 to a=20,
The algorithm described in this work marked 356 events
the value of w2 at a = 3 is greater than a predeter-
out of which 239 turned out to be the epileptic events.
•
mined threshold value T .
2 Thusthe sensitivityofthe algorithmwas70%anditsse-
The threshold values T and T may be considered as
1 2 lectivity was equal to 67%. We then analyzed the same
themodel’sparameterswhichcanbeadjustedtoachieve
tracings with the leading commercial spike detector de-
the desired sensitivity (the ratio of detected epileptic
veloped by the Persyst Development Corporation (In-
events to the total number of epileptic events present
sight 2001.07.12). This software marked 654 events out
in the analyzed EEG tracing) and selectivity (the ratio
of which 268 were epileptic events. Thus slightly better
of epileptic events to the total number of events marked
sensitivityof79%wasachievedattheexpenseofthe low
by the algorithm as epileptic spikes).
41%selectivity. Theperformanceofpreliminarynumeri-
Whilethissimplealgorithmisquiteeffectiveforsimple
cal implementation of the detection algorithm presented
spikes such as one shown in Fig. 1 it fails for the com-
in this work is excellent and allows to process 24 hour
mon case of an epileptic spike accompanied by a slow
EEG recording (19 channels sampled at 240 Hz) in a
wave with comparable amplitude. The example of such
matter of minutes on the average personal computer.
complex is given in Fig. 3(a). The overlap of the neg-
Thegoalofwaveletanalysisofthe twotypesofspikes,
ative tail of the Mexican hat with the slow wave yields
presented in this paper, was to elucidate the approach
theinherentlylowvaluesofw2 atscaleA(panel(b))and
to epileptic events detection which explicitly hinges on
scaleB (panel(c))ascomparedtothosecharacteristicof
the behavior of wavelet power spectrum of EEG signal
the“isolated”spike. Nevertheless,thenormalizedwavelet
across scales and not merely on its values. Thus, this
power does decrease from scale B to C. Consequently,
approach is distinct not only from the detection algo-
in the same vein as the argument we presented above,
rithms based upon discrete multiresolution representa-
we can develop an algorithm which detects the epileptic
tionsofEEGrecordings[22,23,24,25,26]butalsofrom
spike in the vicinity of a slow wave by calculating the
the method developed by Senhadji and Wendling which
following linear combination of wavelet transforms:
employs continuous wavelet transform [27].
W˜(a,t )=c W(a,t )+c W(a ,t +τ) (8) Epilepsy is a common disease which affects 1-2% of
0 1 0 2 s 0
the population and about 4% of children [28]. In some
and checking whether the square of corresponding nor- epilepsy syndromes interictal paroxysmal discharges of
malized power w˜(a,t0) = W˜2(a,t0)/σ2at scales a = 3 cerebral neurons reflect the severity of the epileptic dis-
anda=7exceedsthe thresholdvalueT˜ andT˜ ,respec- order and themselves are believed to contribute to the
1 2
tively. The second term in (8) allows us to detect the progressivedisturbancesincerebralfunctions(eg. speech
slow wave which follows the spike. The parameters a impairment, behavioraldisturbances) [29]. In such cases
s
and τ are chosen to maximize the overlapof the wavelet precise quantitative spike analysis would be extremely
with the slow wave. For the Mexicanhatwe use a =28 important. Theepilepticeventdetectordescribedinthis
s
and τ = 0.125s. By varying the values of coefficients c paper was developed with this particular goal in mind
1
and c , it is possible to control the relative contribution and its applicationto the studies of the Landau-Kleffner
2
of the spike and the slow wave to the linear combination syndrome will be presented elsewhere.
[1] A. Babloyantz, J. M. Salazar, and C. Nicolis, Physics [4] C.J.Stam,E.M.H.vanderLeijR.W.M.Keunen,and
Letters A 111, 152 (1985). D. L. J. Tavy,Theory Biosci. 118, 209 (1999).
[2] B.J. West,M.N.Novaes,and V.Kovcic,in Fractal Ge- [5] G. Widman, T. Schreiber, B. Rehberg, A. Hoeft, and
ometryinBiologicalSystems,editedbyP.M.Iannoccone C. E. Elger, Phys. Rev.E 62, 4898 (2000).
andM.Khokha(CRCPress,BocaRaton,FL,1995),pp. [6] M. C. Casdagli, L. D. Iasemidis, R. S. Savit, R. L.
267–316. Gilmore, S. N. Roper, and J. C. Sackellares, Electroen-
[3] W.S.PritchardandD.W.Duke,Intern.J.Neuroscience cephalogr. Clin. Neurophysiol. 102, 98 (1997).
67, 31 (1992). [7] K. Lehnertz and C. E. Elger, Phys. Rev. Lett. 80, 5019
4
(1998). [18] S. Mallat, A Wavelet Tour of Signal Processing (Aca-
[8] P. L. Nunez, R. Srinivasan, A. F. Westdorp, R. S. Wi- demic Press, San Diego, 1998).
jesinghe, D. M. Tucker, R. B. Silberstein, and P. J. [19] M. Unserand A. Aldroubi, Proc. IEEE 84, 626 (1996).
Cadusch, Electroencephalogr. Clin. Neurophysiol. 103, [20] C. TorrenceandG.P.Compo, Bull.Amer.Meteor.Soc.
499 (1997). 79, 61 (1998).
[9] P. L. Nunez, R. B. Silberstein, Z. P. Shi, M. R. Car- [21] URL http://republika.pl/eegspike/.
penter, R. Srinivasan, D. M. Tucker, S. M. Doran, P. J. [22] S. Blanco, C. E. D’Attellis, S. I. Isaacson, O. A. Rosso,
Cadusch,andR.S.Wijesinghe,Clin.Neurophysiol.110, and R.O. Sirne, PRE 54, 6661 (1996).
469 (1999). [23] C. E. D’Attellis, S. Isaacson, and R. O. Sirne, Ann.
[10] L. Ingber,Physica D 5, 83 (1982). Biomed. Eng. 25, 286 (1997).
[11] I.A.RezekandS.J.Roberts,IEEETrans.Biomed.Eng. [24] F. Sartoretto and M. Ermani, Clin. Neurophysiol. 110,
45, 1186 (1998). 239 (1999).
[12] M.J.A.M.vanPuttenandC.J.Stam,Phys.LettersA [25] M.Calvagno, M.Ermani,R.Rinaldo,andF.Sartoretto,
281, 131 (2001). in Proceedings of the 2000 IEEE International Confer-
[13] J.GotmanandP.Gloor, Electroencephalogr. Clin.Neu- ence onAcoustics, Speech, and Signal Processing (2000),
rophysiol. 41, 513 (1976). vol. 6.
[14] J. Gotman and L. Y. Wang, Electroencephalogr. Clin. [26] J. Guti´errez, R. Alc´antara, and V. Medina, Medical En-
Neurophysiol. 54, 530 (1982). gineering & Physics23, 623 (2001).
[15] J. Gotman and L. Y. Wang, Electroencephalogr. Clin. [27] L.SenhadjiandF.Wendling,Neurophysiol.Clin.32,175
Neurophysiol. 79, 11 (1991). (2002).
[16] L.DiambraandC.P.Malta,Phys.Rev.E59,929(1999). [28] P. Jallon, Epileptic Disord. 4, 1 (2002).
[17] I. Clark, R. Biscay, M. Echeverria, and T. Virues, Com- [29] J. Engel, Epilepsia 42, 796 (2001).
put.Biol. Med. 25, 373 (1995).
5
2
← S
V] 1
m
G [
E
E 0
−1
0 0.5 1 1.5 2
time [s]
30
20
a
10
20 40 60 80 100 120
FIG. 1: Top panel: simple epileptic spike (marked by S) fol-
lowed by two artifacts. Bottom panel: contour map of the
absolute value of the Mexican hat wavelet coefficients (arbi-
trary units)calculated fortheEEG signal shown above. The
shadesofbluecorrespondtolowvaluesandtheshadesofred
to high values.
2
V] (a)
m 1
G [ 0
E
E
−1
0 0.5 1 1.5 2
1000
er scale A (b)
w 500
o
p
velet 00 0.5 1 1.5 2
a10000
w
d scale B (c)
e
z 5000
ali
m
or 0
n 0 0.5 1 1.5 2
of 10000
are scale C (d)
qu 5000
S
0
0 0.5 1 1.5 2
Time [s]
FIG.2: Squareofnormalizedwaveletpowerforthreedifferent
scales A < B < C (Panels (b)-(d)). The EEG signal shown
in panel(a) is the same as the oneused in Fig. 1
6
1
V] (a)
m
G [ 0
E
E
−1
0 0.5 1 1.5 2
50
scale A (b)
er 25
w
o
et p 00 0.5 1 1.5 2
vel 1000
wa scale B (c)
ed 500
z
ali
m 0
or 0 0.5 1 1.5 2
e of n1105000000 scale C (d)
ar
qu 5000
S
pFtsquaIudGnaee.r lose3f:(otbfh()nea-0(o)0sdrlEm)owpfaoillrwiezpatethvdiercwe0ies.sa5pvdciekoifflemeet-prepsanlTorotiawmwbs1eeclr we[asfal]toevorsettAhchoai<smt1sopB.i5gflent<xha.leCTiss.phsiehkoea2.wmTnphliine-
