- Email: [email protected]
Simultaneous multi-channel spikes and inverted spikes in focal epileptic ECoG are more after offset than during the seizure
Biomedical Signal Processing and Control 10 (2014) 58–64
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Short Communication
Simultaneous multi-channel spikes and inverted spikes in focal epileptic ECoG are more after offset than during the seizure Kaushik Majumdar ∗ Indian Statistical Institute, 8th Mile, Mysore Road, Bangalore 560059, India
a r t i c l e
i n f o
Article history: Received 9 May 2013 Received in revised form 31 December 2013 Accepted 2 January 2014 Available online 1 February 2014 Keywords: Electrocorticogram (ECoG) Focal epilepsy Differential operator Noise suppression Spike detection Seizure offset
a b s t r a c t A simple, real time, differential operator based spike detection algorithm has been described, which can efficiently detect seizure spikes in noisy ECoG signals. Simultaneous spikes and inverted spikes have been detected across all focal ECoG channels during preictal, ictal and postictal periods. Out of 79 seizures recorded from 21 patients, about 80% showed more occurrence of simultaneous spikes and inverted spikes across focal channels after the seizure offset than during the seizure in 0–40 Hz range, where the duration studied after the offset is equal to the duration of the seizure. This is an important finding because this goes contrary to the prevailing wisdom that epileptic seizure is a hyper-synchronous phenomenon. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Epileptic seizures are marked by spikes and sharp waves in the electrophysiological signals both on the scalp (Latka et al. [1]; Godoy et al. [2]) and in depth recordings (Godoy et al. [2]). Also epileptic seizures are characterized by hyper-synchronous neuronal activities (Fisher et al. [3]). Single cell studies of human focal seizures in and around the focal region have revealed a complex pattern of synchronous activity in the neurons. In the beginning of the seizure neurons fire heterogeneously. The firing becomes more and more homogeneous as the seizure progresses toward the termination (Truccolo et al. [4]). In other words, single cell recordings show more and more simultaneous or correlated spiking activities toward the end of seizure than in the beginning. Automatic spike detection in electrophysiological signals is an active area of research (Latka et al. [1]; Pan and Tompkin [5]; Obeid and Wolf [6]; Quiroga et al. [7]; Shahid et al. [8]). It has been shown in Majumdar [9] that first and second order differential operators can augment the epileptic spikes relative to the background electrocorticogram (ECoG) signals. Seizure focal ECoG is a convenient signal for analyzing spiking behavior of the epileptic focus. Where multi-channel focal ECoG signal is
∗ Tel.: +91 80 2848 3002x332. E-mail address: [email protected] 1746-8094/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.bspc.2014.01.002
available, studies on how synchronously or asynchronously spikes occur in those channels before, during and after seizure can offer us an insight into the seizure (termination) dynamics at neuronal ensemble level consisting of several thousand cells (Buzsaki et al. [10]). Efficacy of differential operators in accentuating epileptic spikes in the ECoG signals was demonstrated in Majumdar [9] and Majumdar and Vardhan [11]. For a focal ECoG signal x(t), the seizure spikes gets enhanced at x (t) with respect to the background (Majumdar [9]; Majumdar and Vardhan [11]), where x (t) is the first order difference with respect to t. This particular property of differentiation of the seizure EEG has been observed in White et al. [12]. The same is also true for x (t) (Majumdar [9]). Both x (t) and x (t) have been combined to devise a robust method to detect epileptic spikes and inverted spikes in focal ECoG signals, which works with high accuracy even in presence of a good amount of noise. With this novel method simultaneous occurrence of spikes and inverted spikes across all the focal channels has been measured during and after the offset of seizures. Duration after the offset has been taken equal to the duration of the seizure. Out of 79 seizures (there have to be adequate data after the offset, for which 8 other seizures could not be tested) of 21 patients, 82.28% have shown more simultaneous spikes occurred across all 3 focal channels after the offset than during the seizure and 77.22% have shown occurrence of more simultaneous inverted spikes. This happened even for patients whose seizure foci were significantly apart.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64 -14
2
x 10
E 1
B
Amplitude
A 0
F -1
C
-2
D
-3 1
2
3
4
5
6
Sampled time Fig. 1. Plot of six consecutive sample points of a real ECoG signal. ABC is a spike and CDE is an inverted spike. For a spike the signal is an increasing function on the left of the spike reaching maximum at B and then it is a decreasing function on the right of B. For an inverted spike the signal is a decreasing function on the left of the spike reaching minimum at D and then it is an increasing function on the right of D.
2. Method 2.1. Spike detection The spikes and inverted spikes in a digital signal will appear as in Fig. 1. Notice that for a spike ABC, the signal is an increasing function of time on the AB segment reaching maximum at B and then decreasing function on the BC segment. Let the coordinate of B be (m, s(m)), that of A be (m − 1, s(m − 1)) and of C be (m + 1, s(m + 1)). Clearly, s(m) − s(m − 1) > 0 (for increasing function), s(m + 1) − s(m) < 0 (for decreasing function). Note that s (m) = s(m + 1) − s(m) − (s(m) − s(m − 1)) < 0 (for maximum value at B). Let us define a spike-inverted spike operator P as P(x) = s (x) * s (x). Clearly, P(m −) = s (m)(s(m) − s(m − 1)) < 0 and P(m +) = s (m)(s(m + 1) − s(m)) > 0. Therefore a spike at m can be identified by the property P(m−) < 0 & P(m+) > 0 & s (m) < 0,
(1)
where P(m −) signifies value of P at a point less than m, but close to m and P(m +) signifies value of P at a point greater than m, but close to m. Similarly, the inverted spike at D in Fig. 1 with coordinate (k, s(k)) can be identified by the property P(k−) < 0 & P(k+) > 0 & s (k) > 0.
(2)
Difference operation is a linear time operation. Therefore the algorithm described above for spike detection is a linear time or real time operation. This is a big computational advantage over wavelet or independent component analysis (ICA) or principal component analysis (PCA) based algorithms. A wavelet based algorithm, which always involves convolution operation, cannot be executed in less than n log n time, where n is the input size. ICA and PCA based algorithms will take more than O(n2 ) time, for an input size n. Note that a difference operation works as a high-pass filter (Figs. 2(B) and 3(B)). The operator P will therefore work as a high-pass filter leading to allowing the high frequency noise to pass through and suppression of the low frequency information. To avoid the former a low-pass filter should be used before applying P. Regarding the latter, note that the signal morphology is not being used after being operated by P. Only the sign of P and s are being used to identify
59
a spike (condition (1)) or an inverted spike (condition (2)), which remains true for spikes and inverted spikes irrespective of their shape, amplitude or any other morphological properties. This is one major advantage of the proposed algorithm over any other morphology based algorithm. Therefore P acting as a high-pass filter is not affecting the detection performance. There are a plethora of methods to identify spikes and inverted spikes. The easiest among them is possibly to check in every three successive points whether the midpoint value is higher or lower than the values at other two points. But P operator has an advantage of noise reduction before a spike identification (by suppressing the background as can be seen comparing Fig. 3(B) with (A)), which leads to a more accurate identification of simultaneous occurrence of seizure spikes in focal ECoG before, during and after a seizure (Fig. 2(C) and (D)). Seizure spikes usually have higher amplitude than the background spikes and therefore gets enhanced more compared to the background spikes under difference operation as explained above. Notice that before running the P operator for spike detection in all 3 focal ECoG channels the signals were band pass equiripple-FIR filtered between 0 and 40 Hz with stop band attenuation of 40 dB (Fig. 2(C) and (D)). It is clear from Fig. 2 that the number of simultaneous spikes is much higher after the seizure offset than during the seizure. Clearly, P operator is executable in real time. Another advantage of this method is difference operation is inbuilt in most of the scientific computation packages like MATLAB. To further test the efficacy of P operator in spike detection in noisy signals, particularly contaminated ECoG with (electrode scratching) artifacts (average of 3 different hours of ECoG recorded from one focal channel of patient 4) quite akin to seizure spikes was chosen. Each hour contained preictal + ictal + postictal ECoG. In all three hours seizures started at different time points and ended at different time points without any overlap among seizure duration in three different hours. Next, the mean one hour ECoG of the three hour’s recording was taken (Fig. 3(A)). Fig. 3(B) shows the P operated mean. The mean ECoG was low-pass filtered with cutoff frequency 40 Hz with equiripple FIR filter and the P operator was run on the filtered signal (Fig. 3(C)). Notice, in Fig. 3(C) how all the three seizure spike clusters were identified clearly amidst a highly noisy background (Fig. 3(A)). In case x is a focal ECoG signal of an epileptic patient, the ictal spikes of x gets enhanced considerably with respect to the background (which is suppressed to a good extent). The reason is very simple. Let a, b and c are successive time points. If a spike occurs at b then statistically x(b) − x(a) and x(c) − x(b) both have high numerical value, which are two successive points in x . On the other hand if all the three points belong to normal background signal then x(b) − x(a) and x(c) − x(b) both will have small values (actually smaller than the average background signal amplitude). Thus x will enhance the spikes of x, but will suppress its background. This particular property of differentiation of the seizure EEG has already been observed in White et al. [12]. For the same reason x will accentuate the seizure spikes even further. It has been observed that subsequent higher derivatives have only marginal effects on spike enhancement. The ratio of maximum of the P(x) and maximum of x is greater than 1 (in many cases of the order of 103 –104 ), x is one hour long preictal + ictal + postictal ECoG of a focal or nonfocal channel. The ratio of number of zeros in P(x) and in x is also greater than 1 (often of the order of 10). This supports the assertion that the P operator enhances the higher amplitude seizure spikes, but suppresses the lower amplitude background activities. Both the aforementioned ratios are higher in the focal channels than in the nonfocal channels. Higher the values of the ratios are in the focal channels more accurate is the automatic seizure detection (Majumdar and Vardhan [11]).
60
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64 4
(A) Raw single channel focal ECoG signal s
Amplitude
x 10 5
0
-5 0
1
2
3
4
5
8
Amplitude
x 10 2
# simultaneous peaks
7
8
9
10 5
x 10
0
-2 0
# simultaneous troughs
6
Time (B) s'' * s'
1
2
3
4
5
6
7
8
9
Time (C) Number of simultaneous peaks in all focal channels
10 5
x 10
60 40 20 0 0
100
200
300
400
500
600
700
800
900
800
900
Time: 1 time point = 1000 time point window (D) Number of simultaneous troughs in all focal channels 60 40 20 0 0
100
200
300
400
500
600
700
Time: 1 time point = 1000 time point window
Fig. 2. (A) Raw focal ECoG signal s from an epileptic patient with electrode scratching artifacts. The recording is for one hour of preictal, ictal and postictal signals. (B) s * s suppressed the artifacts to a good extent. (C) Number of simultaneous occurrence of spikes (peaks) in all three focal channels (blue) and in all three nonfocal channels (red), and (D) number of simultaneous inverted spikes (troughs) in focal channels (blue) and nonfocal channels (red). The horizontal line in (C) and (D) indicates the value of statistical significance. Parallel vertical lines indicate seizure onset and offset times. In (C) and (D) one time point = a 1000 time point window. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Boolean product of the ith entries of Ur and Vr for all i ∈ {1, . . ., r}. Let it be expressed as Ur | ∗ |Vr = Wr , where |*| denotes the Boolean vector product. It is obvious that the Boolean vector product is valid not just between a pair of vectors, but it can be performed among any N ≥ 2
2.2. Multi-channel extension Definition 1. Let Ur and Vr are two r dimensional vectors with only binary entries, i.e., 0 and 1. Boolean vector product of Ur and Vr means the r dimensional vector Wr , where the ith entry of Wr is the 4
Amplitude
5
(A) Raw average ECoG s
x 10
0
-5 0
1
2
3
4
5
6
7
8
9
10 5
x 10
7
(B) P(s) = s''*s'
x 10 4
Amplitude
2 0 -2 -4
0
1
2
3
4
5
6
7
8
9
10 5
x 10
6
(C) P(band pass filtered s)
x 10 10
Amplitude
5
0
-5
0
1
2
3
4
5
Time
6
7
8
9
10 5
x 10
Fig. 3. (A) One hour average of raw ECoG of a focal channel across three different hours, each containing preictal + ictal + postictal data, denoted by s. (B) P operated s. (C) P operated band pass filtered s between 0 and 40 Hz with 40 dB stop band attenuation. Parallel vertical line-pairs denote onset and offset time points for each of the three seizures.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Spike/inverted spike detection in channel 1
Matching simultaneous spikes/inverted spikes across all channels
Spike/inverted spike detection in channel 2
Counting simultaneous spikes/inverted spikes ni non -overlapping windows
Spike/inverted spike detection in channel 3 Fig. 4. Block diagram of simultaneous spikes/inverted spikes detection and counting across three channels, which can be extended to any number of channels in a straightforward manner.
number of vectors and the result will be the same irrespective of the order in which the vectors have been taken. Simultaneous spikes and inverted spikes detection method across two or more ECoG channels have succinctly been put in the block diagram in Fig. 4. 2.3. Statistical significance What is the guarantee that the simultaneous spikes and inverted spikes detection in a multi-channel recording environment is not merely a random phenomenon? To address this issue a statistical significance test has been devised the way it is usually done for pairwise phase synchronization measure (Tass et al. [13]). To justify the statistical significance of an N channel simultaneous spikes
61
detection task, N shifted surrogate signals are to be generated by randomly shuffling the values (amplitudes) of each of the N signals across the time points. This completely destroys the temporal structure of each signal. Now the spike detection is to be performed (after necessary filtering, exactly the way it has been done for the N signals of interest) on this N surrogate signals. This has to be done for a large number of times (at least 100 times). Then the values have to be sorted in increasing order. The 95th value represents a level of statistical significance that has been used in this paper as has been done in (Tass et al. [13]). This means that the null hypothesis H0 = “Simultaneous occurrence of spikes across all focal channels is only due to chance,” has 0.05 probability to be true. Any number of times of simultaneous occurrence of spikes across the N signals of interest that falls above the level of statistical significance value (only 3 time points across 3 channels for a 1000 time point long window, surrogate data were generated out of the same ECoG data that had to be tested in order to keep the amplitude same) may be accepted as statistically significant. 3. Data The ECoG data (publicly available from Freiburg Seizure Prediction Project [14]) were collected using Neurofile NT digital video EEG system (It-med, Usingen, Germany) with 128 channels, 256 Hz sampling rate, and a 16 bit analog to digital converter. The ECoG from only six sites were made accessible to researchers all over the world that were analyzed. Three of them from the focal areas and the other three from outside the focal areas (this is the configuration available to us). Data for each channel has been provided in a one hour long segment consisting of preictal, ictal and postictal parts. See Table 1 for the patient details. The patient population had earlier been studied in Aschenbrenner-Scheibe [15], Maiwald et al. [16], and Winterhalder et al. [17], where further details can be found. Although there are a total of 87 seizures, only 79 of them could be tested for simultaneous occurrence of spikes and inverted spikes, because 8 of the recordings didn’t have enough
Patient 2, 15th hour, simultaneous peaks focal channels 200
nonfocal channels
# simultaneous peaks
line of statistical significance (horizontal) 150
vertical line 1: seizure onset vertical line 2: seizure offset
100
50
0 0
100
200
300
400
500
600
700
800
900
700
800
900
Patient 2, 15th hour, simultaneous troughs
# simultaneous troughs
200
150
100
50
0 0
100
200
300
400
500
600
Time: 1 time point = 1000 time point window
Fig. 5. (Top) Occurrence of simultaneous spikes (peaks) in all three focal channels (in blue) and all three nonfocal channels (red). (Bottom) Simultaneous inverted spikes (troughs) in all three focal (in blue) and all three nonfocal channels. Parallel vertical lines demarcate seizure onset and offset times. Horizontal line in each plot indicates the statistically significant value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
62
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Table 1 Patient detail. Patient
Gender
Age
Seizure type
H/NC
Electrode
Origin
No of seizures
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
F M M F F F F F M M F F F F M F M F F M M
15 38 14 26 16 31 42 32 44 47 10 42 22 41 31 50 28 25 28 33 13
SP,CP SP,CP,GTC SP,CP SP,CP,GTC SP,CP,GTC CP,GTC SP,CP,GTC SP,CP CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC SP,CP SP,CP,GTC SP,CP,GTC SP,CP
NC H NC H NC H H NC NC H NC H H H and NC H and NC H NC NC NC NC NC
g,s D g,s D,g,s g,s D,g,s D g,s g,s D g,s D,g,s d,s d,s d,s d,s S S S D,g,s g,s
Frontal Temporal Frontal Temporal Frontal Temporo/Occipital Temporal Frontal Temporo/Occipital Temporal Parietal Temporal Temporo/Occipital Fronto/Temporal Temporal Temporal Temporal Frontal Frontal Temporo/Parietal Temporal
4 3 5 5 5 3 3 2 5 5 4 4 2 4 4 5 5 5 4 5 5
SP = simple parietal, CP = complex parietal, GTC = generalized tonic-clonic, H = hippocampal, NC = neocortical. Electrode: grid (g), strip (s), depth (d). Seizure frequency varies between 0.1 and 6.8 per day (Table 1 of [15]).
Data were FIR low-pass filtered at 40 Hz with phase distortion correction. FIR filters have linear phase distortion, to rectify which the filter should be run once from left to right of the signal and then again from right to left. This is accomplished by the filtfilt command in MATLAB. Stop-band attenuation was 40 dB. From Fig. 5 it is clear that simultaneous occurrence of spikes in all three focal channels is more after the offset than during the seizure. The same is true for inverted spikes. Table 2 gives the complete statistics for all the 79 seizures recorded from 21 patents. The trend that is reflected in Table 2 is a robust one in the sense that it was not influenced by random noise. In order to ensure this, the test was repeated by replacing the preictal + ictal + postictal ECoG of one randomly chosen focal channel with a purely interictal ECoG of same duration and from the same channel, but from an arbitrarily chosen recording hour (Fig. 6). In almost all cases number of simultaneous spikes and inverted spikes came down significantly (Fig. 6), yet the trend remained the same (i.e., number of simultaneous spikes is more after the seizure offset than during the seizure, where the duration after the offset is same as the seizure duration, etc.). This is happening because two out of the three channels are still focal channels. On the other hand if the time duration was altered this trend was easily changed. Same behavior was observed in nonfocal channels as well. In Fig. 6(A) simultaneous occurrence of spikes is shown in three seizure focal channels during a one hour recording of Table 2 Percentage of the cases of simultaneous occurrence of more spikes and inverted spikes after the seizure offset across all focal and nonfocal channels in the 79 cases that have been studied.
Focal Nonfocal
Spikes
Inverted spikes
82.28% 55.7%
77.22% 49.37%
A 30
# simultaneous events
4. Results
preictal + ictal + postictal ECoG. Fig. 6(B) shows the same except for a channel the recording was replaced by a one hour pure interictal recording. In Fig. 6(B) the number of simultaneous occurrence of spikes has sharply gone down compared to Fig. 6(A), which shows simultaneous occurrence of spikes goes up around the time of occurrence of an epileptic seizure (Wyler et al. [18]). Fig. 7 shows the focal channel locations in the left hippocampus of the ECoG signals of patient 7 (Table 1) that were analyzed to generate Fig. 6. The artifacts are much less in the deep brain region (Ball et al. [19]), but the volume conduction is likely to be high due to proximity of the focal channels. Still the recording hour plays a significant role in generating simultaneous spikes, which supports the claim that the simultaneous spikes are due to seizure, not due to random noise. Number of spikes and inverted spikes in a signal also depends on the filter used on the signal. However using equiripple FIR band
20
10
0 0
100
200
300
400
500
600
700
800
900
500
600
700
800
900
B # simultaneous events
signal after the seizure offset so that a duration equal to the seizure duration could be tested after the offset. Apart from 2 to 5 h of preictal + ictal + postictal recordings for each patient at least 24 h of purely interictal recordings were also provided. There was no interictal recording available for patient 2.
30 20 10 0 0
100
200
300
400
Time: 1 point = 1000 time point window Fig. 6. (A) Simultaneous spikes in all three focal channels (number of spikes counted in every 1000 time point window) during a one hour recording of preictal + ictal + postictal ECoG. (B) Simultaneous spikes during one hour recording as in (A), but one channel’s ECoG is replaced by one hour purely interictal recording in the same channel. Discontinuation of signals between 500 and 600 time points is due to ECoG amplifier disconnection.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Fig. 7. Focal channel locations (red) of patient 7 (Table 1) for the ECoG analyzed in Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
pass filter for 0–40 Hz band and using Gaussian low pass filter with cutoff frequency 40 Hz the trends that have been shown in Table 2 remain the same. For both the filters stop band attenuation was same. Simultaneous occurrence of spikes and inverted spikes measure in more than two channels is generally different from correlation measure as shown in Fig. 8. Here the correlation among three focal channels was measured as the highest eigen value of the
63
symmetric matrix of the correlation coefficients of channel pairs. Each time point is a 1000 time point window in the filtered data. The correlation matrix is for the entire 1000 time point window. Highest eigen value of each matrix was chosen and plotted over time (Fig. 8(A)), which gives the dominant correlation trend among the channels. Multi-channel correlation of focal epileptic ECoG signals was calculated the same way by Schindler et al. [20]. Simultaneous inverted spikes (Fig. 8(B)) and spikes (Fig. 8(C)) have been calculated in the same window. It was observed that correlation among the focal ECoG channels go up toward the end of seizures, rather than in the beginning (Schindler et al. [20]), which is compatible with the trends of simultaneous occurrence of spikes and inverted spikes in Table 2 despite the measures are different. Since ECoG spikes are largely due to simultaneous depolarization of a large number of neurons and ECoG inverted spikes are mainly due to simultaneous hyperpolarization of a large number of neurons, it is natural that number of spikes and inverted spikes will more or less match with each other. This is evident in Fig. 8(B) and (C). Research on seizure termination is relatively less compared to seizure initiation. In this paper our concern is entirely on seizure offset rather than on the onset. From Table 1 it is clear that out of a total of 21 patients there were 11 with neocortical origin of seizure, 8 with hippocampal origin and 2 with both hippocampal and neocortical origin. 6 out of 11 patients with neocortical origin of seizure had larger number of spikes and inverted spikes after offset than during the seizure in all the seizures that were tested. 74% of all the seizures that were tested in the 11 patients showed the trend. 5 out of 8 patients with hippocampal origin of seizure showed the trend in all the seizures that were tested. 86% of the seizures that were tested in the 8 patients had larger number of spikes and inverted spikes after the offset than during the seizure. 2 patients with both neocortical and hippocampal origin of seizure showed a weaker trend of having more spikes and inverted spikes
Correlation
(A) Correlation among three focal channles
2 1 0 0
100
200
300
400
500
600
(B) Simultaneous troughs across three focal channels
# troughs
80 60 40 20 0 0
100
0
100
200 300 400 500 (C) Simultaneous peaks across three focal channels
600
# peaks
80 60 40 20 0 200
300
400
500
600
Time: 1 time point = 1000 time point window Fig. 8. (A) Plot of highest eigen value of correlation matrix of pairwise correlation coefficients of focal channels over time. (B) Simultaneous occurrence of inverted spikes (troughs) across focal channels. (C) Simultaneous occurrence of spikes (peaks) across focal channels. All the three graphs are cut between time points between 580 and 630 because of discontinuation of signals due to ECoG amplifier disconnection.
64
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
after the offset than during the seizure, where duration for measure of spikes and inverted spikes after the offset is same as the duration of seizure. None of the two patients showed the trend for all seizures that were tested. Only 4 out of all 7 seizures tested from these two patients, i.e., 57% showed the trend of having more spikes and inverted spikes after the offset than during the seizure. More seizures from more patients with both hippocampal and neocortical origin of seizures need to be tested to ascertain if the trend of having more spikes and inverted spikes after offset than during the seizure is weaker in those patients. 5. Conclusion It is widely accepted that epileptic seizure is a hypersynchronous phenomenon (Fisher et al. [3]), which implies simultaneous occurrence of spikes and inverted spikes should be more during the seizure than after its termination. However in this work it has been shown that simultaneous occurrence of spikes and inverted spikes in focal ECoG channels is more after the seizure offset than during the seizure for an overwhelming majority of the cases across wide ranging focal onset seizures (see Table 1 for the diversity of seizure types), where the duration taken after the offset is equal to the seizure duration. This is somewhat similar to the finding that amplitude correlation across the focal ECoG channels goes up toward the seizure offset rather than soon after the onset (Schindler et al. [20]). This happens despite the fact that amplitude correlation and simultaneous occurrence of spikes and inverted spikes are different measures. This is also compatible with the recent finding that in epileptogenic foci single neurons fire heterogeneously in the beginning of the seizure, but more homogeneously toward the seizure offset (Truccolo et al. [4]). In patients with frequent spontaneous focal seizures, firing of single neurons is reported to be correlated with local ECoG spikes (Wyler et al. [18]). Also high degree of homogeneity in neuronal firing is necessary within the epileptogenic focus for having spikes in the focal ECoG (Wyler et al. [18]). From Fig. 2 it is clear that even after the post offset spikes are largely suppressed due to filtering the number of spikes remains higher after the offset than during the seizure. Figs. 2 and 3 show that the post offset spikes are rather diminished (otherwise they would have been enhanced by the P operator), and therefore may be due to homogeneous firing of a relatively smaller number of neurons in the epileptogenic foci. Possibility of existence of a small network responsible for termination of the seizures can be a potential target for intervention by new drugs in hitherto pharmacologically intractable epilepsies. Single cell studies like Truccolo et al. [4] in this network may offer new insights into the seizure dynamics with potential benefit for the patients with intractable epilepsy. Simultaneous neuronal activities are captured more clearly in simultaneous occurrence of spikes and inverted spikes across multiple channels than phase synchronization and multi-channel correlation even in the ECoG data where seizure characteristics are not clearly visible and phase synchronization among focal ECoG channels is weak or statistically insignificant during seizure (for example, in case of patient 5).
Acknowledgement Freiburg Seizure Prediction Project is acknowledged for providing the ECoG data. References [1] M. Latka, Z. Was, A. Kozik, B.J. West, Wavelet analysis of epileptic spikes, Phys. Rev. E 67 (2003) 052902. [2] J. Godoy, H. Luders, D.S. Dinner, H.H. Morris, E. Wyllie, D. Murphy, Significance of sharp waves in routine EEGs after epilepsy surgery, Epilepsia 33 (2) (1992) 285–288. [3] R.S. Fisher, W. van Emde Boas, W. Blume, C. Elger, P. Genton, P. Lee, J. Engel Jr., Epileptic seizures and epilepsy: definitions proposed by the International League Against Epilepsy (ILAE) and the International Bureau for Epilepsy (IBE), Epilepsia 46 (4) (2005) 470–472. [4] W. Truccolo, J.A. Donoghue, L.R. Hochberg, E.N. Eskandar, J.R. Madsen, W.A. Anderson, E.N. Brown, E. Halgren, S.S. Cash, Single-neuron dynamics in human focal epilepsy, Nature Neuroscience 14 (2011) 635–641. [5] J. Pan, W.J. Tompkins, A real-time QRS detection algorithm, IEEE Transactions on Biomedical Engineering 32 (3) (1985) 230–236. [6] I. Obeid, P.D. Wolf, Evaluation of spike-detection algorithms for a brainmachine interface application, IEEE Transactions on Biomedical Engineering 51 (6) (2004) 905–911. [7] R.Q. Quiroga, Z. Nadasdy, Y. Ben-Shaul, Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering, Neural Computation 16 (2004) 1661–1687. [8] S. Shahid, J. Walker, L.S. Smith, A new spike detection algorithm for extracellular neural recordings, IEEE Transactions on Biomedical Engineering 57 (4) (2010) 853–866. [9] K. Majumdar, Differential operator in seizure detection, Computers in Biology and Medicine 42 (1) (2012) 70–74. [10] G. Buzsaki, C.A. Anastassiou, C. Koch, The origin of extracellular fields and currents – EEG, ECoG, LFP and spikes, Nature Reviews Neuroscience 13 (2012) 407–420. [11] K. Majumdar, P. Vardhan, Automatic seizure detection in ECoG by differential operator and windowed variance, IEEE Transactions on Neural Systems and Rehabilitation Engineering 19 (4) (2011) 356–365. [12] A.D. White, P.A. Williams, D.J. Ferraro, S. Clark, S.D. Kadam, F.E. Dudek, K.J. Staley, Efficient unsupervised algorithms for the detection of seizures in continuous recordings from rats after brain injury, Journal of Neuroscience Methods 152 (2006) 256–266. [13] P. Tass, M.G. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volkmann, A. Schnitzler, H.-J. Freund, Detection of n:m phase locking from noisy data: application to magnetoencephalography, Physical Review Letters 81 (15) (1998) 3291–3294. [14] Freiburg Seizure Prediction Project, Freiburg, Germany, 2008 https:// epilepsy.uni-freiburg.de/ [15] R. Aschenbrenner-Scheibe, T. Maiwald, M. Winterhalder, H.U. Voss, J. Timmer, A. Schulze-Bonhage, How well can epileptic seizures be predicted? An evaluation of a nonlinear method, Brain 126 (2003) 2616–2626. [16] T. Maiwald, M. Winterhalder, R. Aschenbrenner-Scheibe, H.U. Voss, A. SchulzeBonhage, J. Timmer, Comparison of three nonlinear seizure prediction methods by means of the seizure prediction characteristics, Physica D 194 (2004) 357–368. [17] M. Winterhalder, B. Schelter, T. Maiwald, A. Brandt, A. Schad, A. SchulzeBonhage, J. Timmer, Spatio-temporal patient-individual assessment of synchronization changes for epileptic seizure prediction, Clinical Neurophysiology 117 (2006) 2399–2413. [18] A.R. Wyler, G.A. Ojemann, A.A. Ward, Neurons in human epileptic cortex: correlation between unit and EEG activity, Annals of Neurology 11 (3) (2004) 301–308. [19] T. Ball, M. Kern, I. Mutschler, A. Aertsen, A. Schulze-Bonhage, Signal quality of simultaneously recorded invasive and non-invasive EEG, Neuroimage 46 (2009) 708–719. [20] K. Schindler, H. Leung, C.E. Elger, K. Lehnertz, Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG, Brain 130 (2007) 65–77.
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc
Short Communication
Simultaneous multi-channel spikes and inverted spikes in focal epileptic ECoG are more after offset than during the seizure Kaushik Majumdar ∗ Indian Statistical Institute, 8th Mile, Mysore Road, Bangalore 560059, India
a r t i c l e
i n f o
Article history: Received 9 May 2013 Received in revised form 31 December 2013 Accepted 2 January 2014 Available online 1 February 2014 Keywords: Electrocorticogram (ECoG) Focal epilepsy Differential operator Noise suppression Spike detection Seizure offset
a b s t r a c t A simple, real time, differential operator based spike detection algorithm has been described, which can efficiently detect seizure spikes in noisy ECoG signals. Simultaneous spikes and inverted spikes have been detected across all focal ECoG channels during preictal, ictal and postictal periods. Out of 79 seizures recorded from 21 patients, about 80% showed more occurrence of simultaneous spikes and inverted spikes across focal channels after the seizure offset than during the seizure in 0–40 Hz range, where the duration studied after the offset is equal to the duration of the seizure. This is an important finding because this goes contrary to the prevailing wisdom that epileptic seizure is a hyper-synchronous phenomenon. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction Epileptic seizures are marked by spikes and sharp waves in the electrophysiological signals both on the scalp (Latka et al. [1]; Godoy et al. [2]) and in depth recordings (Godoy et al. [2]). Also epileptic seizures are characterized by hyper-synchronous neuronal activities (Fisher et al. [3]). Single cell studies of human focal seizures in and around the focal region have revealed a complex pattern of synchronous activity in the neurons. In the beginning of the seizure neurons fire heterogeneously. The firing becomes more and more homogeneous as the seizure progresses toward the termination (Truccolo et al. [4]). In other words, single cell recordings show more and more simultaneous or correlated spiking activities toward the end of seizure than in the beginning. Automatic spike detection in electrophysiological signals is an active area of research (Latka et al. [1]; Pan and Tompkin [5]; Obeid and Wolf [6]; Quiroga et al. [7]; Shahid et al. [8]). It has been shown in Majumdar [9] that first and second order differential operators can augment the epileptic spikes relative to the background electrocorticogram (ECoG) signals. Seizure focal ECoG is a convenient signal for analyzing spiking behavior of the epileptic focus. Where multi-channel focal ECoG signal is
∗ Tel.: +91 80 2848 3002x332. E-mail address: [email protected] 1746-8094/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.bspc.2014.01.002
available, studies on how synchronously or asynchronously spikes occur in those channels before, during and after seizure can offer us an insight into the seizure (termination) dynamics at neuronal ensemble level consisting of several thousand cells (Buzsaki et al. [10]). Efficacy of differential operators in accentuating epileptic spikes in the ECoG signals was demonstrated in Majumdar [9] and Majumdar and Vardhan [11]. For a focal ECoG signal x(t), the seizure spikes gets enhanced at x (t) with respect to the background (Majumdar [9]; Majumdar and Vardhan [11]), where x (t) is the first order difference with respect to t. This particular property of differentiation of the seizure EEG has been observed in White et al. [12]. The same is also true for x (t) (Majumdar [9]). Both x (t) and x (t) have been combined to devise a robust method to detect epileptic spikes and inverted spikes in focal ECoG signals, which works with high accuracy even in presence of a good amount of noise. With this novel method simultaneous occurrence of spikes and inverted spikes across all the focal channels has been measured during and after the offset of seizures. Duration after the offset has been taken equal to the duration of the seizure. Out of 79 seizures (there have to be adequate data after the offset, for which 8 other seizures could not be tested) of 21 patients, 82.28% have shown more simultaneous spikes occurred across all 3 focal channels after the offset than during the seizure and 77.22% have shown occurrence of more simultaneous inverted spikes. This happened even for patients whose seizure foci were significantly apart.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64 -14
2
x 10
E 1
B
Amplitude
A 0
F -1
C
-2
D
-3 1
2
3
4
5
6
Sampled time Fig. 1. Plot of six consecutive sample points of a real ECoG signal. ABC is a spike and CDE is an inverted spike. For a spike the signal is an increasing function on the left of the spike reaching maximum at B and then it is a decreasing function on the right of B. For an inverted spike the signal is a decreasing function on the left of the spike reaching minimum at D and then it is an increasing function on the right of D.
2. Method 2.1. Spike detection The spikes and inverted spikes in a digital signal will appear as in Fig. 1. Notice that for a spike ABC, the signal is an increasing function of time on the AB segment reaching maximum at B and then decreasing function on the BC segment. Let the coordinate of B be (m, s(m)), that of A be (m − 1, s(m − 1)) and of C be (m + 1, s(m + 1)). Clearly, s(m) − s(m − 1) > 0 (for increasing function), s(m + 1) − s(m) < 0 (for decreasing function). Note that s (m) = s(m + 1) − s(m) − (s(m) − s(m − 1)) < 0 (for maximum value at B). Let us define a spike-inverted spike operator P as P(x) = s (x) * s (x). Clearly, P(m −) = s (m)(s(m) − s(m − 1)) < 0 and P(m +) = s (m)(s(m + 1) − s(m)) > 0. Therefore a spike at m can be identified by the property P(m−) < 0 & P(m+) > 0 & s (m) < 0,
(1)
where P(m −) signifies value of P at a point less than m, but close to m and P(m +) signifies value of P at a point greater than m, but close to m. Similarly, the inverted spike at D in Fig. 1 with coordinate (k, s(k)) can be identified by the property P(k−) < 0 & P(k+) > 0 & s (k) > 0.
(2)
Difference operation is a linear time operation. Therefore the algorithm described above for spike detection is a linear time or real time operation. This is a big computational advantage over wavelet or independent component analysis (ICA) or principal component analysis (PCA) based algorithms. A wavelet based algorithm, which always involves convolution operation, cannot be executed in less than n log n time, where n is the input size. ICA and PCA based algorithms will take more than O(n2 ) time, for an input size n. Note that a difference operation works as a high-pass filter (Figs. 2(B) and 3(B)). The operator P will therefore work as a high-pass filter leading to allowing the high frequency noise to pass through and suppression of the low frequency information. To avoid the former a low-pass filter should be used before applying P. Regarding the latter, note that the signal morphology is not being used after being operated by P. Only the sign of P and s are being used to identify
59
a spike (condition (1)) or an inverted spike (condition (2)), which remains true for spikes and inverted spikes irrespective of their shape, amplitude or any other morphological properties. This is one major advantage of the proposed algorithm over any other morphology based algorithm. Therefore P acting as a high-pass filter is not affecting the detection performance. There are a plethora of methods to identify spikes and inverted spikes. The easiest among them is possibly to check in every three successive points whether the midpoint value is higher or lower than the values at other two points. But P operator has an advantage of noise reduction before a spike identification (by suppressing the background as can be seen comparing Fig. 3(B) with (A)), which leads to a more accurate identification of simultaneous occurrence of seizure spikes in focal ECoG before, during and after a seizure (Fig. 2(C) and (D)). Seizure spikes usually have higher amplitude than the background spikes and therefore gets enhanced more compared to the background spikes under difference operation as explained above. Notice that before running the P operator for spike detection in all 3 focal ECoG channels the signals were band pass equiripple-FIR filtered between 0 and 40 Hz with stop band attenuation of 40 dB (Fig. 2(C) and (D)). It is clear from Fig. 2 that the number of simultaneous spikes is much higher after the seizure offset than during the seizure. Clearly, P operator is executable in real time. Another advantage of this method is difference operation is inbuilt in most of the scientific computation packages like MATLAB. To further test the efficacy of P operator in spike detection in noisy signals, particularly contaminated ECoG with (electrode scratching) artifacts (average of 3 different hours of ECoG recorded from one focal channel of patient 4) quite akin to seizure spikes was chosen. Each hour contained preictal + ictal + postictal ECoG. In all three hours seizures started at different time points and ended at different time points without any overlap among seizure duration in three different hours. Next, the mean one hour ECoG of the three hour’s recording was taken (Fig. 3(A)). Fig. 3(B) shows the P operated mean. The mean ECoG was low-pass filtered with cutoff frequency 40 Hz with equiripple FIR filter and the P operator was run on the filtered signal (Fig. 3(C)). Notice, in Fig. 3(C) how all the three seizure spike clusters were identified clearly amidst a highly noisy background (Fig. 3(A)). In case x is a focal ECoG signal of an epileptic patient, the ictal spikes of x gets enhanced considerably with respect to the background (which is suppressed to a good extent). The reason is very simple. Let a, b and c are successive time points. If a spike occurs at b then statistically x(b) − x(a) and x(c) − x(b) both have high numerical value, which are two successive points in x . On the other hand if all the three points belong to normal background signal then x(b) − x(a) and x(c) − x(b) both will have small values (actually smaller than the average background signal amplitude). Thus x will enhance the spikes of x, but will suppress its background. This particular property of differentiation of the seizure EEG has already been observed in White et al. [12]. For the same reason x will accentuate the seizure spikes even further. It has been observed that subsequent higher derivatives have only marginal effects on spike enhancement. The ratio of maximum of the P(x) and maximum of x is greater than 1 (in many cases of the order of 103 –104 ), x is one hour long preictal + ictal + postictal ECoG of a focal or nonfocal channel. The ratio of number of zeros in P(x) and in x is also greater than 1 (often of the order of 10). This supports the assertion that the P operator enhances the higher amplitude seizure spikes, but suppresses the lower amplitude background activities. Both the aforementioned ratios are higher in the focal channels than in the nonfocal channels. Higher the values of the ratios are in the focal channels more accurate is the automatic seizure detection (Majumdar and Vardhan [11]).
60
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64 4
(A) Raw single channel focal ECoG signal s
Amplitude
x 10 5
0
-5 0
1
2
3
4
5
8
Amplitude
x 10 2
# simultaneous peaks
7
8
9
10 5
x 10
0
-2 0
# simultaneous troughs
6
Time (B) s'' * s'
1
2
3
4
5
6
7
8
9
Time (C) Number of simultaneous peaks in all focal channels
10 5
x 10
60 40 20 0 0
100
200
300
400
500
600
700
800
900
800
900
Time: 1 time point = 1000 time point window (D) Number of simultaneous troughs in all focal channels 60 40 20 0 0
100
200
300
400
500
600
700
Time: 1 time point = 1000 time point window
Fig. 2. (A) Raw focal ECoG signal s from an epileptic patient with electrode scratching artifacts. The recording is for one hour of preictal, ictal and postictal signals. (B) s * s suppressed the artifacts to a good extent. (C) Number of simultaneous occurrence of spikes (peaks) in all three focal channels (blue) and in all three nonfocal channels (red), and (D) number of simultaneous inverted spikes (troughs) in focal channels (blue) and nonfocal channels (red). The horizontal line in (C) and (D) indicates the value of statistical significance. Parallel vertical lines indicate seizure onset and offset times. In (C) and (D) one time point = a 1000 time point window. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Boolean product of the ith entries of Ur and Vr for all i ∈ {1, . . ., r}. Let it be expressed as Ur | ∗ |Vr = Wr , where |*| denotes the Boolean vector product. It is obvious that the Boolean vector product is valid not just between a pair of vectors, but it can be performed among any N ≥ 2
2.2. Multi-channel extension Definition 1. Let Ur and Vr are two r dimensional vectors with only binary entries, i.e., 0 and 1. Boolean vector product of Ur and Vr means the r dimensional vector Wr , where the ith entry of Wr is the 4
Amplitude
5
(A) Raw average ECoG s
x 10
0
-5 0
1
2
3
4
5
6
7
8
9
10 5
x 10
7
(B) P(s) = s''*s'
x 10 4
Amplitude
2 0 -2 -4
0
1
2
3
4
5
6
7
8
9
10 5
x 10
6
(C) P(band pass filtered s)
x 10 10
Amplitude
5
0
-5
0
1
2
3
4
5
Time
6
7
8
9
10 5
x 10
Fig. 3. (A) One hour average of raw ECoG of a focal channel across three different hours, each containing preictal + ictal + postictal data, denoted by s. (B) P operated s. (C) P operated band pass filtered s between 0 and 40 Hz with 40 dB stop band attenuation. Parallel vertical line-pairs denote onset and offset time points for each of the three seizures.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Spike/inverted spike detection in channel 1
Matching simultaneous spikes/inverted spikes across all channels
Spike/inverted spike detection in channel 2
Counting simultaneous spikes/inverted spikes ni non -overlapping windows
Spike/inverted spike detection in channel 3 Fig. 4. Block diagram of simultaneous spikes/inverted spikes detection and counting across three channels, which can be extended to any number of channels in a straightforward manner.
number of vectors and the result will be the same irrespective of the order in which the vectors have been taken. Simultaneous spikes and inverted spikes detection method across two or more ECoG channels have succinctly been put in the block diagram in Fig. 4. 2.3. Statistical significance What is the guarantee that the simultaneous spikes and inverted spikes detection in a multi-channel recording environment is not merely a random phenomenon? To address this issue a statistical significance test has been devised the way it is usually done for pairwise phase synchronization measure (Tass et al. [13]). To justify the statistical significance of an N channel simultaneous spikes
61
detection task, N shifted surrogate signals are to be generated by randomly shuffling the values (amplitudes) of each of the N signals across the time points. This completely destroys the temporal structure of each signal. Now the spike detection is to be performed (after necessary filtering, exactly the way it has been done for the N signals of interest) on this N surrogate signals. This has to be done for a large number of times (at least 100 times). Then the values have to be sorted in increasing order. The 95th value represents a level of statistical significance that has been used in this paper as has been done in (Tass et al. [13]). This means that the null hypothesis H0 = “Simultaneous occurrence of spikes across all focal channels is only due to chance,” has 0.05 probability to be true. Any number of times of simultaneous occurrence of spikes across the N signals of interest that falls above the level of statistical significance value (only 3 time points across 3 channels for a 1000 time point long window, surrogate data were generated out of the same ECoG data that had to be tested in order to keep the amplitude same) may be accepted as statistically significant. 3. Data The ECoG data (publicly available from Freiburg Seizure Prediction Project [14]) were collected using Neurofile NT digital video EEG system (It-med, Usingen, Germany) with 128 channels, 256 Hz sampling rate, and a 16 bit analog to digital converter. The ECoG from only six sites were made accessible to researchers all over the world that were analyzed. Three of them from the focal areas and the other three from outside the focal areas (this is the configuration available to us). Data for each channel has been provided in a one hour long segment consisting of preictal, ictal and postictal parts. See Table 1 for the patient details. The patient population had earlier been studied in Aschenbrenner-Scheibe [15], Maiwald et al. [16], and Winterhalder et al. [17], where further details can be found. Although there are a total of 87 seizures, only 79 of them could be tested for simultaneous occurrence of spikes and inverted spikes, because 8 of the recordings didn’t have enough
Patient 2, 15th hour, simultaneous peaks focal channels 200
nonfocal channels
# simultaneous peaks
line of statistical significance (horizontal) 150
vertical line 1: seizure onset vertical line 2: seizure offset
100
50
0 0
100
200
300
400
500
600
700
800
900
700
800
900
Patient 2, 15th hour, simultaneous troughs
# simultaneous troughs
200
150
100
50
0 0
100
200
300
400
500
600
Time: 1 time point = 1000 time point window
Fig. 5. (Top) Occurrence of simultaneous spikes (peaks) in all three focal channels (in blue) and all three nonfocal channels (red). (Bottom) Simultaneous inverted spikes (troughs) in all three focal (in blue) and all three nonfocal channels. Parallel vertical lines demarcate seizure onset and offset times. Horizontal line in each plot indicates the statistically significant value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
62
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Table 1 Patient detail. Patient
Gender
Age
Seizure type
H/NC
Electrode
Origin
No of seizures
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
F M M F F F F F M M F F F F M F M F F M M
15 38 14 26 16 31 42 32 44 47 10 42 22 41 31 50 28 25 28 33 13
SP,CP SP,CP,GTC SP,CP SP,CP,GTC SP,CP,GTC CP,GTC SP,CP,GTC SP,CP CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC CP,GTC SP,CP,GTC SP,CP,GTC SP,CP,GTC SP,CP SP,CP,GTC SP,CP,GTC SP,CP
NC H NC H NC H H NC NC H NC H H H and NC H and NC H NC NC NC NC NC
g,s D g,s D,g,s g,s D,g,s D g,s g,s D g,s D,g,s d,s d,s d,s d,s S S S D,g,s g,s
Frontal Temporal Frontal Temporal Frontal Temporo/Occipital Temporal Frontal Temporo/Occipital Temporal Parietal Temporal Temporo/Occipital Fronto/Temporal Temporal Temporal Temporal Frontal Frontal Temporo/Parietal Temporal
4 3 5 5 5 3 3 2 5 5 4 4 2 4 4 5 5 5 4 5 5
SP = simple parietal, CP = complex parietal, GTC = generalized tonic-clonic, H = hippocampal, NC = neocortical. Electrode: grid (g), strip (s), depth (d). Seizure frequency varies between 0.1 and 6.8 per day (Table 1 of [15]).
Data were FIR low-pass filtered at 40 Hz with phase distortion correction. FIR filters have linear phase distortion, to rectify which the filter should be run once from left to right of the signal and then again from right to left. This is accomplished by the filtfilt command in MATLAB. Stop-band attenuation was 40 dB. From Fig. 5 it is clear that simultaneous occurrence of spikes in all three focal channels is more after the offset than during the seizure. The same is true for inverted spikes. Table 2 gives the complete statistics for all the 79 seizures recorded from 21 patents. The trend that is reflected in Table 2 is a robust one in the sense that it was not influenced by random noise. In order to ensure this, the test was repeated by replacing the preictal + ictal + postictal ECoG of one randomly chosen focal channel with a purely interictal ECoG of same duration and from the same channel, but from an arbitrarily chosen recording hour (Fig. 6). In almost all cases number of simultaneous spikes and inverted spikes came down significantly (Fig. 6), yet the trend remained the same (i.e., number of simultaneous spikes is more after the seizure offset than during the seizure, where the duration after the offset is same as the seizure duration, etc.). This is happening because two out of the three channels are still focal channels. On the other hand if the time duration was altered this trend was easily changed. Same behavior was observed in nonfocal channels as well. In Fig. 6(A) simultaneous occurrence of spikes is shown in three seizure focal channels during a one hour recording of Table 2 Percentage of the cases of simultaneous occurrence of more spikes and inverted spikes after the seizure offset across all focal and nonfocal channels in the 79 cases that have been studied.
Focal Nonfocal
Spikes
Inverted spikes
82.28% 55.7%
77.22% 49.37%
A 30
# simultaneous events
4. Results
preictal + ictal + postictal ECoG. Fig. 6(B) shows the same except for a channel the recording was replaced by a one hour pure interictal recording. In Fig. 6(B) the number of simultaneous occurrence of spikes has sharply gone down compared to Fig. 6(A), which shows simultaneous occurrence of spikes goes up around the time of occurrence of an epileptic seizure (Wyler et al. [18]). Fig. 7 shows the focal channel locations in the left hippocampus of the ECoG signals of patient 7 (Table 1) that were analyzed to generate Fig. 6. The artifacts are much less in the deep brain region (Ball et al. [19]), but the volume conduction is likely to be high due to proximity of the focal channels. Still the recording hour plays a significant role in generating simultaneous spikes, which supports the claim that the simultaneous spikes are due to seizure, not due to random noise. Number of spikes and inverted spikes in a signal also depends on the filter used on the signal. However using equiripple FIR band
20
10
0 0
100
200
300
400
500
600
700
800
900
500
600
700
800
900
B # simultaneous events
signal after the seizure offset so that a duration equal to the seizure duration could be tested after the offset. Apart from 2 to 5 h of preictal + ictal + postictal recordings for each patient at least 24 h of purely interictal recordings were also provided. There was no interictal recording available for patient 2.
30 20 10 0 0
100
200
300
400
Time: 1 point = 1000 time point window Fig. 6. (A) Simultaneous spikes in all three focal channels (number of spikes counted in every 1000 time point window) during a one hour recording of preictal + ictal + postictal ECoG. (B) Simultaneous spikes during one hour recording as in (A), but one channel’s ECoG is replaced by one hour purely interictal recording in the same channel. Discontinuation of signals between 500 and 600 time points is due to ECoG amplifier disconnection.
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
Fig. 7. Focal channel locations (red) of patient 7 (Table 1) for the ECoG analyzed in Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
pass filter for 0–40 Hz band and using Gaussian low pass filter with cutoff frequency 40 Hz the trends that have been shown in Table 2 remain the same. For both the filters stop band attenuation was same. Simultaneous occurrence of spikes and inverted spikes measure in more than two channels is generally different from correlation measure as shown in Fig. 8. Here the correlation among three focal channels was measured as the highest eigen value of the
63
symmetric matrix of the correlation coefficients of channel pairs. Each time point is a 1000 time point window in the filtered data. The correlation matrix is for the entire 1000 time point window. Highest eigen value of each matrix was chosen and plotted over time (Fig. 8(A)), which gives the dominant correlation trend among the channels. Multi-channel correlation of focal epileptic ECoG signals was calculated the same way by Schindler et al. [20]. Simultaneous inverted spikes (Fig. 8(B)) and spikes (Fig. 8(C)) have been calculated in the same window. It was observed that correlation among the focal ECoG channels go up toward the end of seizures, rather than in the beginning (Schindler et al. [20]), which is compatible with the trends of simultaneous occurrence of spikes and inverted spikes in Table 2 despite the measures are different. Since ECoG spikes are largely due to simultaneous depolarization of a large number of neurons and ECoG inverted spikes are mainly due to simultaneous hyperpolarization of a large number of neurons, it is natural that number of spikes and inverted spikes will more or less match with each other. This is evident in Fig. 8(B) and (C). Research on seizure termination is relatively less compared to seizure initiation. In this paper our concern is entirely on seizure offset rather than on the onset. From Table 1 it is clear that out of a total of 21 patients there were 11 with neocortical origin of seizure, 8 with hippocampal origin and 2 with both hippocampal and neocortical origin. 6 out of 11 patients with neocortical origin of seizure had larger number of spikes and inverted spikes after offset than during the seizure in all the seizures that were tested. 74% of all the seizures that were tested in the 11 patients showed the trend. 5 out of 8 patients with hippocampal origin of seizure showed the trend in all the seizures that were tested. 86% of the seizures that were tested in the 8 patients had larger number of spikes and inverted spikes after the offset than during the seizure. 2 patients with both neocortical and hippocampal origin of seizure showed a weaker trend of having more spikes and inverted spikes
Correlation
(A) Correlation among three focal channles
2 1 0 0
100
200
300
400
500
600
(B) Simultaneous troughs across three focal channels
# troughs
80 60 40 20 0 0
100
0
100
200 300 400 500 (C) Simultaneous peaks across three focal channels
600
# peaks
80 60 40 20 0 200
300
400
500
600
Time: 1 time point = 1000 time point window Fig. 8. (A) Plot of highest eigen value of correlation matrix of pairwise correlation coefficients of focal channels over time. (B) Simultaneous occurrence of inverted spikes (troughs) across focal channels. (C) Simultaneous occurrence of spikes (peaks) across focal channels. All the three graphs are cut between time points between 580 and 630 because of discontinuation of signals due to ECoG amplifier disconnection.
64
K. Majumdar / Biomedical Signal Processing and Control 10 (2014) 58–64
after the offset than during the seizure, where duration for measure of spikes and inverted spikes after the offset is same as the duration of seizure. None of the two patients showed the trend for all seizures that were tested. Only 4 out of all 7 seizures tested from these two patients, i.e., 57% showed the trend of having more spikes and inverted spikes after the offset than during the seizure. More seizures from more patients with both hippocampal and neocortical origin of seizures need to be tested to ascertain if the trend of having more spikes and inverted spikes after offset than during the seizure is weaker in those patients. 5. Conclusion It is widely accepted that epileptic seizure is a hypersynchronous phenomenon (Fisher et al. [3]), which implies simultaneous occurrence of spikes and inverted spikes should be more during the seizure than after its termination. However in this work it has been shown that simultaneous occurrence of spikes and inverted spikes in focal ECoG channels is more after the seizure offset than during the seizure for an overwhelming majority of the cases across wide ranging focal onset seizures (see Table 1 for the diversity of seizure types), where the duration taken after the offset is equal to the seizure duration. This is somewhat similar to the finding that amplitude correlation across the focal ECoG channels goes up toward the seizure offset rather than soon after the onset (Schindler et al. [20]). This happens despite the fact that amplitude correlation and simultaneous occurrence of spikes and inverted spikes are different measures. This is also compatible with the recent finding that in epileptogenic foci single neurons fire heterogeneously in the beginning of the seizure, but more homogeneously toward the seizure offset (Truccolo et al. [4]). In patients with frequent spontaneous focal seizures, firing of single neurons is reported to be correlated with local ECoG spikes (Wyler et al. [18]). Also high degree of homogeneity in neuronal firing is necessary within the epileptogenic focus for having spikes in the focal ECoG (Wyler et al. [18]). From Fig. 2 it is clear that even after the post offset spikes are largely suppressed due to filtering the number of spikes remains higher after the offset than during the seizure. Figs. 2 and 3 show that the post offset spikes are rather diminished (otherwise they would have been enhanced by the P operator), and therefore may be due to homogeneous firing of a relatively smaller number of neurons in the epileptogenic foci. Possibility of existence of a small network responsible for termination of the seizures can be a potential target for intervention by new drugs in hitherto pharmacologically intractable epilepsies. Single cell studies like Truccolo et al. [4] in this network may offer new insights into the seizure dynamics with potential benefit for the patients with intractable epilepsy. Simultaneous neuronal activities are captured more clearly in simultaneous occurrence of spikes and inverted spikes across multiple channels than phase synchronization and multi-channel correlation even in the ECoG data where seizure characteristics are not clearly visible and phase synchronization among focal ECoG channels is weak or statistically insignificant during seizure (for example, in case of patient 5).
Acknowledgement Freiburg Seizure Prediction Project is acknowledged for providing the ECoG data. References [1] M. Latka, Z. Was, A. Kozik, B.J. West, Wavelet analysis of epileptic spikes, Phys. Rev. E 67 (2003) 052902. [2] J. Godoy, H. Luders, D.S. Dinner, H.H. Morris, E. Wyllie, D. Murphy, Significance of sharp waves in routine EEGs after epilepsy surgery, Epilepsia 33 (2) (1992) 285–288. [3] R.S. Fisher, W. van Emde Boas, W. Blume, C. Elger, P. Genton, P. Lee, J. Engel Jr., Epileptic seizures and epilepsy: definitions proposed by the International League Against Epilepsy (ILAE) and the International Bureau for Epilepsy (IBE), Epilepsia 46 (4) (2005) 470–472. [4] W. Truccolo, J.A. Donoghue, L.R. Hochberg, E.N. Eskandar, J.R. Madsen, W.A. Anderson, E.N. Brown, E. Halgren, S.S. Cash, Single-neuron dynamics in human focal epilepsy, Nature Neuroscience 14 (2011) 635–641. [5] J. Pan, W.J. Tompkins, A real-time QRS detection algorithm, IEEE Transactions on Biomedical Engineering 32 (3) (1985) 230–236. [6] I. Obeid, P.D. Wolf, Evaluation of spike-detection algorithms for a brainmachine interface application, IEEE Transactions on Biomedical Engineering 51 (6) (2004) 905–911. [7] R.Q. Quiroga, Z. Nadasdy, Y. Ben-Shaul, Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering, Neural Computation 16 (2004) 1661–1687. [8] S. Shahid, J. Walker, L.S. Smith, A new spike detection algorithm for extracellular neural recordings, IEEE Transactions on Biomedical Engineering 57 (4) (2010) 853–866. [9] K. Majumdar, Differential operator in seizure detection, Computers in Biology and Medicine 42 (1) (2012) 70–74. [10] G. Buzsaki, C.A. Anastassiou, C. Koch, The origin of extracellular fields and currents – EEG, ECoG, LFP and spikes, Nature Reviews Neuroscience 13 (2012) 407–420. [11] K. Majumdar, P. Vardhan, Automatic seizure detection in ECoG by differential operator and windowed variance, IEEE Transactions on Neural Systems and Rehabilitation Engineering 19 (4) (2011) 356–365. [12] A.D. White, P.A. Williams, D.J. Ferraro, S. Clark, S.D. Kadam, F.E. Dudek, K.J. Staley, Efficient unsupervised algorithms for the detection of seizures in continuous recordings from rats after brain injury, Journal of Neuroscience Methods 152 (2006) 256–266. [13] P. Tass, M.G. Rosenblum, J. Weule, J. Kurths, A. Pikovsky, J. Volkmann, A. Schnitzler, H.-J. Freund, Detection of n:m phase locking from noisy data: application to magnetoencephalography, Physical Review Letters 81 (15) (1998) 3291–3294. [14] Freiburg Seizure Prediction Project, Freiburg, Germany, 2008 https:// epilepsy.uni-freiburg.de/ [15] R. Aschenbrenner-Scheibe, T. Maiwald, M. Winterhalder, H.U. Voss, J. Timmer, A. Schulze-Bonhage, How well can epileptic seizures be predicted? An evaluation of a nonlinear method, Brain 126 (2003) 2616–2626. [16] T. Maiwald, M. Winterhalder, R. Aschenbrenner-Scheibe, H.U. Voss, A. SchulzeBonhage, J. Timmer, Comparison of three nonlinear seizure prediction methods by means of the seizure prediction characteristics, Physica D 194 (2004) 357–368. [17] M. Winterhalder, B. Schelter, T. Maiwald, A. Brandt, A. Schad, A. SchulzeBonhage, J. Timmer, Spatio-temporal patient-individual assessment of synchronization changes for epileptic seizure prediction, Clinical Neurophysiology 117 (2006) 2399–2413. [18] A.R. Wyler, G.A. Ojemann, A.A. Ward, Neurons in human epileptic cortex: correlation between unit and EEG activity, Annals of Neurology 11 (3) (2004) 301–308. [19] T. Ball, M. Kern, I. Mutschler, A. Aertsen, A. Schulze-Bonhage, Signal quality of simultaneously recorded invasive and non-invasive EEG, Neuroimage 46 (2009) 708–719. [20] K. Schindler, H. Leung, C.E. Elger, K. Lehnertz, Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG, Brain 130 (2007) 65–77.