Despiking SEEG signals reveals dynamics of gamma band preictal activity

We constructed simulations that reproduce a preictal series of epileptiform discharges, inspired by real SEEG traces, with a sampling rate set to 512 Hz. In the real data (see section 2.2), we observed oscillations at different frequencies in the preictal period. Within the gamma band, oscillatory bursts were observed around 55 Hz and around 85 Hz (supplementary figure 1 (stacks.iop.org/PM/38/N42/mmedia)). Simulated data consisted of four components, corresponding to spikes, oscillations, seizure and noise. 100 realizations were constructed, one realization is shown on figure 1.

Zoom In Zoom Out Reset image size

Spikes were simulated as a sum of two gamma functions, using the gampdf(x, a, b) function in Matlab. With time variable x in samples, we set (a, b) to (9.8, 3.9) in channel 1, (7.8, 2) in channel 2, (11.7, 3.9) in channel 3 and (11.7, 2) in channel 4. Amplitudes were set to 10, with fluctuation ranging from 10% to 25% (uniform distribution), except for the first gamma function which presented no fluctuation. Time of occurrence of spikes was every 0.8 s, with jitter ranging between 40 and 0 ms, except for the first gamma function.

Oscillations were simulated as oscillatory pulses, using the gausspuls(t, fc, bw) function of matlab, with t time, fc frequency of oscillation set to 55 Hz in channel 3 and 85 Hz in channel 4 and bw fractional bandwidth was set to 0.15 (corresponding to 7 periods of oscillations in the segment corresponding to the full-width half maximum). Amplitude fluctuations were set to 25% and 15% in channel 3 and 4 respectively, maximum temporal jitter to 50 ms and 40 ms in channel 3 and 4 respectively. The seizure part was simulated as two sine waves at 55 Hz in channel 3 and 85 Hz in channel 4, with envelope set to a Tukey window (tukeywin function of Matlab).

Noise was obtained by smoothing random Gaussian noise (function randn of Matlab) with a moving average of 10 points. Noise was set to an amplitude of 2 during the spike period, and 1 during the seizure period. A smooth transition between the two periods was ensured thanks to a half Tukey window, in order to avoid filter ringing arising from abrupt transitions.

We selected presurgical intracerebral signals (stereo-electroencephalography, SEEG) in one patient with drug-resistant focal epilepsy from the Clinical Neurophysiology Department of the Timone hospital in Marseille. This epilepsy was symptomatic of a focal cortical dysplasia localized in the right occipito-temporal junction. The placement of depth electrodes was based upon available non-invasive information providing hypotheses about the localization of the epileptogenic zone. This patient presented a focal cortical dysplasia in the right lateral occipito-temporal junction. The SEEG data was acquired on a Deltamed System, sampled at 256 Hz, with anti-aliasing low-pass analog filter set to 85 Hz. This particular SEEG recording was selected for the current study because it presented clear preictal patterns with regular spiking and visible epileptic oscillations. SEEG was using eight multi-contact depth electrodes (diameter 0.8 mm, 10–15 contacts, each 2 mm long with an inter-contact distance of 1.5 mm), positioned in the right hemisphere, implanted according to Talairach's stereotaxic method (Talairach and Bancaud 1973). Patients signed informed consent, and the study was approved by the Institutional Review board (IRB00003888) of INSERM (IORG0003254, FWA00005831).

We first detected spikes as events with high amplitude (in absolute value). In order to emphasize the transient part of the spikes, signals were first filtered in the 5 Hz–40 Hz band. Amplitude thresholds (on the signed values) were found based on data quartiles for each channel:

$$\begin{eqnarray}&&\text{th}{{\text{r}}_{\text{high}}}={{Q}_{0.5}}+d\left({{Q}_{0.75}}-{{Q}_{0.25}}\right)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\text{th}{{\text{r}}_{\text{low}}}={{Q}_{0.5}}-d\left({{Q}_{0.75}}-{{Q}_{0.25}}\right)\end{eqnarray} \tag{ 2 }$$

With Q_p value in the distribution of amplitudes corresponding to the p percentile, and d distance parameter to be tuned heuristically. Typically, we found a value of 3 to be useful in most cases. Detections were found as the local peaks in the signal passing the amplitude thresholds, with a minimum distance between consecutive peaks set to 10 ms (function 'findpeaks' of Matlab). Then, for each detection, a model consisting of a modified Gaussian function was fitted non-linearly to the data:

$$\begin{eqnarray}&&G(t,a,b,\gamma )=\left\{\begin{array}{c} -A{{\text{e}}^{-\frac{{{\left[(t-a)\ast \gamma \right]}^{2}}}{b}}},~t<a \\ A{{\text{e}}^{-\frac{{{\left[(t-a)/\gamma \right]}^{2}}}{b}}},~t\geqslant a \end{array}\right. \end{eqnarray} \tag{ 3 }$$

With A amplitude, t time, a shift parameter (corresponding to zero of time), b scale parameters (width of the model function), γ asymmetry parameter (1 corresponds to the original—symmetric—Gaussian function). In order to estimate the parameters (a, b, γ) corresponding to the best fit, we used the non-linear simplex method (fminsearch from Matlab). Parameter A was fitted linearly by a least square method, with all other parameters fixed. The whole procedure results in a model of the spike components of the signal.

Time-frequency analysis was performed based on a Morlet wavelet transform with oscillation parameter ξ  =  5 (Bénar et al 2009). Data were then presented under the form of spatio-temporal maps (#channels x #time points). In a first representation, data were filtered in the frequency band of the oscillations (David et al 2011), using the EEGLAB matlab function 'eegfilt' (Delorme et al 2004). A Hilbert transform was applied in order to estimate the envelope of fluctuations, followed by a smoothing of 256 ms. In a second representation this transform was normalized by dividing the map obtained at a lower frequency band, as in Bartolomei et al (2008).

In order to compute connectivity graphs, we used the nonlinear correlation h² (van der Heyden et al 1999, Wendling et al 2000). Window size was set to 2 s, the maximum delay to 100 ms, with a heuristic threshold set to 0.35 on the original signals and 0.15 on the despiked signals.

The robustness of the procedure was validated on a series of 100 realizations, in the two channels with a mixture of oscillations and spikes (channels 3 and 4). For each realization, we computed the correlation coefficient between the simulated data filtered in the band (35–95 Hz) on the one hand and the simulated oscillations in the preictal phase on the other hand. This was done both on the original data and on the despiked data.

We first tested the impact of the despiking procedure on the filtered data in the temporal domain, on one channel of simulated data (channel 4, presenting a mixture of spikes and oscillations at 85 Hz). Results are presented in figure 2. On the original signals (figure 2(a)), the (sharp) spikes produced filter ringing (i.e. spurious oscillations) that arise from the impulse response of the filter (figure 2(b)), of shorter duration but higher amplitude than the actual oscillation. Such oscillations are often referred to as 'Gibbs effect'—although originally the Gibbs phenomenon was obtained from Fourier' Series—and have been called 'false ripples' (Bénar et al 2010) in the context of epilepsy. This implies that energy will be detected in the gamma band whenever spikes are present, whether accompanied by actual oscillations or not. In contrast, for the despiked signal (figure 2(c)), the impact of the high-frequency part of the spike has been drastically reduced (figure 2(d)). This implies that the filtering procedure will be much more specific for detecting actual oscillations, reducing the false detections produced by filter ringing. In supplementary figure 2, we illustrate this effect for simulated oscillations with higher frequency (in the fast ripple band), in both a channel with a mixture of spikes and oscillation and a channel with no true oscillation.

Zoom In Zoom Out Reset image size

We computed the ratio between the amplitude of the despiked signal versus that of the original signal, at the peak of the first (i.e. sharper) wave. The median ratio across channel, events and realizations was 4.6% (quartiles 1.7%–7.5%), which means that most of the time of reduction of the spike of more than 90% was obtained. We did the same computation for the envelope of signals filtered in the 35–95 Hz band. In channel 1, the median ratio was 48.9% (quartiles 47.8%–50%). In channel 2, with sharper spikes, the median ratio was 13.7 (quartiles 12.6%–15.3%).

As the distinction between spikes and oscillations is often better done in the time-frequency plane (van der Heyden et al 1999, Wendling et al 2000, Bénar et al 2010, Jmail et al 2011), we studied the impact of despiking signals on the time-frequency representation of signals. An example is shown in figure 3, corresponding also to channel 4. On the original signals (figure 3(a)), spikes result in a wide-band pattern, and one can easily picture than depending on the time lag between spike and oscillations, the spike part could cover the 'blob' corresponding to the oscillation (Jmail et al 2011). On the despiked signals (figure 3(b)), the oscillatory part is more visible on the time-frequency representation.

Zoom In Zoom Out Reset image size

A useful representation in clinical practice is based on a spatio-temporal map, where signals are filtered in a frequency band of interest (David et al 2011), with possible division by energy in low frequencies (Bartolomei et al 2008). This help mapping channels where the seizure starts, where typically there is a sudden increase in high frequency content and lowering of low frequencies (Bartolomei et al 2008). In figure 4, we illustrate such representation on the simulated data, with or without despiking. On the original signal, when filtered in the frequency band that contains oscillations (35–95 Hz), a high energy level is visible on all channels presenting spikes, whether presenting actual oscillations or not (figure 4(a)). After despiking, there is better discrimination of channels with actual oscillations (figure 4(b))—only little residuals from filtered spikes remain. When normalizing by low frequencies, all channels present low values, even when there are actual oscillations. Indeed, spikes present energy at all frequencies, and the division by energy in low frequencies thus prevents detecting high frequency activity (figure 4(c)). The best representations (in terms of specificity) are obtained for maps after despiking (figures 4(b) and (d)).

Zoom In Zoom Out Reset image size

Finally, we were interested in quantifying the gain in accuracy of reconstructing for the oscillatory part of the signal obtained by the despiking procedure. We computed the correlation of both original signal and the despiked signal with the actual simulated oscillatory component, across 100 realizations and for the channels with oscillations (channels 3 and 4). The median correlation coefficient was 0.70 for the original signal (quartiles 0.57–0.82) and 0.93 for the despiked signal (quartiles 0.92–0.94), showing an improvement which is both highly significant, and robust to different realizations of noise (figure 5).

Zoom In Zoom Out Reset image size

Real data consisted of a section of intracerebral EEG data, in the period immediately preceding the epileptic seizure (i.e. the 'preictal' period). In figure 6, we show the spatio-temporal distribution of automatic detections. A large number of spikes were detected (range 0–82, median 34), mostly in the preictal phase and in most channels.

Zoom In Zoom Out Reset image size

Figure 7 shows a time-frequency analysis in one channel presenting both spikes and oscillations, on a subset of the time window. One can observe that in the despiked signals there almost no spikes are left; only oscillations with some slow waves (which do not impact time-frequency representation in the high frequency range) remain. In the time frequency representation of the despiked signal (figure 7(d)), clear bursts of high frequency activity are now visible (indicated by white arrows), which were masked in the original representation (figure 7(d)) due to the high frequency part of the spikes. To be noted, some real spikes can present larger slow waves than what we decided to include in our simulations (figures 1 and 2). However, this has no impact in the frequency bands of interest (above 15 Hz).

Zoom In Zoom Out Reset image size

Finally, we computed the spatio-temporal maps on all channels, filtered in the 65–85 Hz band after normalization by the 8–30 Hz band (figure 8). Interestingly, a buildup of gamma oscillations can be observed during the preictal period only in the despiked signals (figure 8(b)), whereas in the original signals this activity is masked by the high energy of the low frequency part of the spikes. Interestingly, several channels present energy in high frequency band during the seizure onset, but no preictal oscillations, which suggest that seizure onset and interictal oscillations may be different biomarkers of the epileptogenic zone.

Zoom In Zoom Out Reset image size