Signal processing and computational modeling for interpretation of SEEG-recorded interictal epileptiform discharges in epileptogenic and non-epileptogenic zones

During the pre-surgical evaluation of patients with drug-resistant epilepsy, stereo-electroencephalographic (SEEG, depth electrodes) recordings are performed to identify regions involved in epileptic seizures. Identification of epileptogenic zones (EZs) and non-epileptogenic zones (NEZs, also referred to as propagation, irritative or non-involved zones) is usually based on pre-ictal and ictal patterns and ictal propagation. This categorization is used for deciding the best therapeutic strategy.

Numerous studies have been dedicated to interictal paroxysmal events in humans [1]. Interictal epileptiform discharges (IEDs) observed in SEEG recordings demonstrate a wide range of morphology, ranging from simple monophasic epileptic spikes to more complex multiphasic transient events. Generation of IEDs in partial epilepsies is commonly ascribed to enhanced excitatory interactions within glutamatergic neural networks. However, both human data and animal models have supported the view that gamma-aminobutyric acidergic (GABAergic) signaling does play a significant role [2, 3]. The presence of heterogeneous neural firing patterns during IEDs in epileptic and non-epileptic zones suggests a multitude of agents taking part in IED generation [4], including paroxysmal depolarization shifts (PDSs) [5, 6] and depolarizing GABAergic actions [7, 8]. So far, mechanisms leading to the generation and propagation of SEEG-recorded IEDs remain unclear, and their clinical value for predicting the epileptogenicity level of brain regions, is uncertain.

This paper focuses on a specific IED class called biphasic epileptic spikes, typically observed in SEEG recordings. These events, often referred to as spike-waves (SWs), are characterized by a short-duration spike followed by a longer duration wave, both with the same polarity [9]. Figure 1 provides an example of interictal activity recorded by SEEG electrodes from a patient suffering from temporal lope epilepsy (TLE). As part of a general clinical procedure, the electrode implantation is specific to the patient and it depends on the EZ location provided by non-invasive methods: e.g. electroencephalogram (EEG), MRI, and computed tomography scan. Figure 1 displays the selected bipolar channels where the IEDs are less affected by volume conduction. We recall that a bipolar channel captures the local potential (defined as the electric potential difference between adjacent contacts of an implanted electrode). After this initial selection, SEEG experts labeled the bipolar channels located in the EZ and NEZ. The red and green markers correspond to the SWs in the EZ and NEZ regions, respectively. By visual inspection, we remark that spike and wave amplitudes are higher in the EZ than in the NEZ with narrower spikes and wider waves. These findings are supported by literature by Alarçon et al [10] and more recently by Serafini [11], both investigating the origin and propagation of IEDs in the acute electrocorticogram. In addition to these morphological differences, SWs are more frequent in the EZ, and they can appear as clusters and be preceded and/or followed by spiking activity.

Figure 1. Example SEEG interictal activity recorded in a patient. The visualized SEEG recording is reduced to a small subset of bipolar channels for simplicity. In practice, such SEEG recordings are performed on a more significant number of recording channels (up to 128–256). The red and green markers correspond to spike-waves in the epileptogenic zone (EZ) and non-epileptogenic zone (NEZ). Electrodes are located as follows. H: Heschl gyrus; T: superior temporal gyrus (anterior part); A: amygdala (mesial contacts), middle temporal gyrus (lateral contacts); B: anterior hippocampus (mesial contacts), middle temporal gyrus (lateral contacts); C: posterior hippocampus (mesial contacts), middle temporal gyrus (lateral contacts); TP: temporal pole; TB: temporo-basal region; OF: orbito-frontal region.

Standard image High-resolution image

In the present manuscript, biphasic IEDs are analyzed in-depth and the recorded IEDs are classified according to their morphologies, durations, and generation sites. Their role in interictal epileptogenic networks is discussed accordingly to identify the patient’s EZ versus NEZs. The discriminative role of IEDs’ spatio-temporal features in the assessment of the EZ and NEZ is compared to SEEG experts’ analysis. Subsequently, the physiological hypotheses behind IEDs morphology are formulated using laminar neural mass models (NMMs).

2.1. Patients and SEEG recordings

Five patients (labeled from P1 to P5) undergoing pre-surgical evaluation of drug-resistant partial epilepsy were selected for the study. Their clinical, anatomical, and electrophysiological data showed TLE in which parahippocampal and neocortical structures were involved. The group of patients in whom we analyzed SWs was homogeneous in terms of underlying pathology. The selection was based on precise criteria. First, patients were MRI-negative showing no focal cortical dysplasia. Second, very frequent SWs were recorded at electrode contacts located in mesial temporal lobe structures as well as in temporal neocortical structures.

SEEG recordings were performed using intracerebral multiple contact electrodes (10–15 contacts, length: 2 mm, diameter: 0.8 mm, 1.5 mm apart) placed following Talairach’s stereotactic method [12, 13]. The positioning of electrodes was determined based on the available non-invasive information and clinical hypotheses about the localization of the EZ. The SEEG signals were recorded on a Deltamed system with a maximum of 128 channels. They were sampled at 256 Hz using no digital filter. For each patient, a mean of five SEEG recordings were analyzed (see table 1 for details). For each patient, interictal and ictal (where at least one seizure was recorded) periods were present. Seizures recorded during presurgical video-EEG monitoring were relatively reproducible although drug withdrawal had an impact on seizure patterns, as expected. The SEEG experts used the ictal periods to assess the epileptogenicity level of the selected bipolar channels. The interictal periods were exploited for the SWs detection and the SWs features extraction in order to affiliate the bipolar channel to the EZ versus the NEZ.

Table 1. Information of ictal and interictal recordings per patient.

PatientP1P2P3P4P5 Ictal records24141Average ictal records duration (min)2538602230Inter-ictal records52965Average inter-ictal records duration (min)5260565360SEEG contacts931169511662Detected SWs per bipolar channel over all interictal records3850555710339874269Total detected SWs on interictal records358 02464 402674 809462 518264 6832.2. SEEG signals analysis2.2.1. SWs detection

A bipolar montage was adopted for the SEEG analysis, and prior to the investigation, a notch finite impulse response (FIR) filter was applied to remove the 50 Hz power supply component from signals. Among the IEDs detection approaches [14], the one proposed in [15] was implemented in the present work. This approach was adopted because it compares the input signal to a family of wavelets whose shape resembles the SWs morphology (complex Mexican hat wavelet). This comparison allows for enhancing the detection of IEDs in the records. In more detail, this method uses the mean value of the squared modulus of a wavelet filter banks output to increase the signal-to-noise ratio, as in [16]. The obtained signal was then used as input to the Page-Hinkley algorithm [17, 18]. This method builds a monotonically decreasing statistic, sample by sample from the input signal. In correspondence to an SWs, the behavior of the statistic changes because the signal is above its mean: in fact, the metric starts to increase. It is at this point that the algorithm decides whether a spike is present, by comparing the increase of the statistic with a pre-defined threshold. The algorithm parameters were adjusted by visual inspection to minimize false-negative and false-positive detections.

2.2.2. Bipolar channels selection

Generally speaking, SEEG recordings make use of multiple electrodes (usually between five and ten), and each electrode contains up to fifteen contacts. Among all the bipolar channel recordings, the objective of this step is to retain only those pairs of electrodes contact revealing a large number of IEDs. Some contacts may record the same IEDs generators: to avoid this problem the same strategy was adopted for all five patients considered in the study. It consisted in discarding bipolar channels containing IEDs generated by the same generator (usually adjacent plots referring to similar brain structures) or generated by volume conduction, and retaining only bipolar channels with a maximum number of IEDs. Table 1 provides the average number of SWs detected per interictal record and patient. On average, 365 000 SWs per patient were detected and exploited in the first screening. It is worth mentioning that the SWs distribution across the bipolar channels was not uniform. For this reason, the first screening of bipolar channels was performed by selecting only those showing the most significant number of IEDs with respect to the adjacent SEEG contacts.

2.2.3. EZ vs. NEZ bipolar channels discrimination

A team of SEEG experts classified the selected bipolar channels of SEEG recordings performed in each patient as being located in the epileptogenic (EZ) [19] versus propagation or not involved zone (NEZ), by visual inspection. Ictal data was used to determine the degree of ‘epileptogenicity’ of a given region and its subsequent contribution to the EZ. The experts considered three criteria: the capability of a given region to generate high-frequency oscillations (rapid discharges) at seizure onset [18, 2022], the delay of involvement of the region to the seizure onset and the seizure-related signal amplitude. SEEG experts worked independently and afterward, an agreement was reached to choose two bipolar channels per patient: one located in the EZ and another one in the NEZ.

Table 2 summarizes the selection of bipolar channels by SEEG experts. Except for A1–A2 and B1–B2, all channels are located in neocortical structures. Table 2 also provides the number of SWs detected per patient and bipolar channel. It is possible to observe the occurrence variability of interictal events for each patient and a higher number of SWs in the EZ than in the NEZ (except for P3).

Table 2. SEEG experts selected one bipolar channel located in the EZ and one positioned in the NEZ for each patient. Electrode locations: TB, temporal basalis area; TP, temporal pole; A, amygdala complex; PH: parahippocampal gyrus; B, anterior part of Hippocampus. The number SWs detected per bipolar channel is reported together with the total of SWs detected per patient in both EZ and NEZ.

PatientP1P2P3P4P5 ZoneEZn-EZEZn-EZEZn-EZEZn-EZEZn-EZChannelTB1–TB2TB9–TB10TP2–TP3TB5–TB6A1–A2TB6–TB7TP9–TP10PH3–PH4B1–B2A9–A10SWs491823862679607782111 0069810319111 33710 562SWs Sum7304328618 82713 00121 8992.2.4. SWs mean waveform computation

The last step consisted of computing the mean waveforms of the detected SWs per selected bipolar channel per record. This task aimed to identify qualitative morphological differences among SWs. For each record, segments of 1.5 s centered on the SWs detections were extracted and normalized using the z-score. A first SW mean waveform was obtained by averaging all channel detections. Then, each detection was aligned to the maximum of the cross-correlation function, computed pairwise between the detection and the preliminary average waveform. The final waveform was computed by averaging all the aligned detections. As shown in table 1 several thousands of SWs (3953 on average) were exploited to compute the mean SW waveform per bipolar channel.

By averaging waveforms of IEDs over time, we assumed that the morphology of these is stationary over time. We provided some evidence of this assumption by also visualizing a small subset of SWs (n = 30), randomly chosen across the time range for the same bipolar channels.

2.2.5. SWs features extraction and statistical tests

For each patient, we randomly selected 100 SWs at the EZ bipolar channel and 100 SWs at the NEZ bipolar channel to extract the SW features. If more than one interictal period was recorded in a patient, we randomly extracted the same amount of SWs from each period. Bipolar channels were visually inspected to avoid IEDs polarity inversion: when it was the case, the signal was inverted for SWs features extraction. An expert visually verified each SW in the set of 100 samples. False detections were discarded and replaced by the correct ones. For a selected SW, several features were extracted: spike and wave amplitude (calculated from the baseline depicted in red, in figure 2), spike and wave full width at half maximum (FWHM, related to the sharpness of the spike and duration of the wave), spike peak to wave peak time delay (SW Delay), spike FWHM to wave FWHM time delay (FWHM Delay), the ratio between FWHM wave and FWHM spike, the ratio between Spike Amplitude and Wave Amplitude and the ratio between FWHM wave and FHWM Delay. The selected features are shown in figure 2 and enumerated in table 3.

Figure 2. Extracted featured per SW: spike and wave amplitude, spike and wave full width at half maximum (FWHM), spike peak to wave peak time delay (SW Delay), spike FWHM to wave FWHM time delay (FWHM Delay), the ratio between FWHM wave and FWHM spike, the ratio between Spike Amplitude and Wave Amplitude and the ratio between FWHM wave and FHWM Delay.

Standard image High-resolution image

Table 3. Spike and wave extracted features.

FeaturesUnitAbbreviation Spike AmplitudemicroVolt (µV)Spike AmpWave AmplitudemicroVolt (µV)Wave AmpSpike-Wave Delaysecond (s)SW DelaySpike Full Width at Half Maximumsecond (s)FWHM SpikeWave Full Width at Half Maximumsecond (s)FWHM WaveDelay between FWHM Spike and FWHM Wavesecond (s)FWHM DelayRatio between Spike Amp and Wave Amp—Spike Amp/Wave AmpRatio between FWHM Wave and FWHM Spike—FWHM Wave/FWHM SpikeRatio between FWHM Wave and FWHM Delay—FWHM Wave/FWHM Delay

Features were z-scored per patient: each feature value was normalized by subtracting its mean and dividing by its standard deviation, both computed for every patient. For example, a Spike-Wave Delay value equal to −0.9 indicates that the time delay between the spike and wave peaks is almost −1 standard deviation less with respect to the mean value of the delay variable computed on the patient data. This step was introduced to compare SWs features extracted from different patients. Thanks to the z-score, we did not compare the absolute value of features but normalized deviations from their mean. This allowed considering patient-specific physiology. In other words, we were able, for example, to assess similarity in the deviation from the mean for the feature FWHM Spike, in the EZ and NEZ across different patients. The same comparison, but using features absolute values would have been more difficult because of inter-patient SWs morphology variability.

We used the Wilcoxon signed-rank test and the Kolmogorov–Smirnov test to compare the EZ versus the NEZ features distributions. The Kolmogorov–Smirnov test verifies the null hypothesis that the EZ (and NEZ) samples distribution follows a Gaussian law for a given feature. The Wilcoxon signed-rank test verifies the null hypothesis that two related paired samples come from the same distribution. A p-value less than 0.05 gives confidence of 95% to reject the null hypothesis. We tested if feature distributions followed a Gaussian law or not (Kolmogorov Smirnov test), and selected a nonparametric test (Wilcoxon signed-rank test) instead of Student’s t-test in the latter case. Tests were performed in both ways: per patient and by merging all patients’ SWs features data together.

2.2.6. SWs features clustering

We used k-means unsupervised clustering to discern if the SW morphology could be discriminative of the affiliation of SW to the EZ or NEZ. The clustering problem was performed on combined data from all patients and on data per patient to understand patient-specific features variations. The k-means clustering algorithm was applied with default parameters [23] by setting the number of clusters to two. In order to reduce the grouping dimensionality of the proposed features, we performed the clustering task by using tuples of three and two features.

Our interest was to establish if interictal activity can reveal concordant information on the EZ and NEZ to partition the data in a manner similar to the SEEG experts’ analysis. In fact, the choice of EZ and NEZ was performed a priori on seizure data and then compared to the analysis performed on interictal data.

2.3. Computational model

NMMs are mathematical representations of the dynamics of the mean activity of synaptically connected neuronal subpopulations [24, 25]. The neural mass formalism assumes that these populations are almost homogeneous and ‘synchronized’, therefore, their behaviors can be reflected by average firing rates and membrane perturbations. NMMs have been extensively used to study both physiological [26] and pathological brain rhythms [2729]. In the NMM considered in the present study, the activity of each subpopulation is given by a sigmoid ‘pulse to wave’ function that transforms the synaptic inputs v into a firing rate:

Equation (1)

A second-order differential equation relates the average presynaptic firing rate to post-synaptic potential (PSP):

Equation (2)

where W represents the average synaptic gain, and τw is the synaptic time constant lumping the rise and decay times (assumed equal), and axonal delays. The time constant parameter τw reflects the kinetics of glutamatergic and GABAergic PSPs of cortical neurons. Given that an NMM formulates the ‘synchronized’ average temporal activity of a homogeneous population, single cell synaptic time constants can be conveyed (to some extent) to the model. On the other hand, parameters like average synaptic gain and connectivity coefficients are rather kept at an abstract level (they are ‘lumped’ parameters collecting various aspects not explicitly described by the model). Appendix D provides the system equations and parameters used in the manuscript.

The generation of realistic SEEG data to compare with real recordings requires combining the NMM framework with a physical formalism that can account for the location of the SEEG electrode contacts, the biophysical features of the volume conductor (tissue morphology and conductivity), and the physics of electrical current propagation deriving from Maxwell’s equation. The integration of these aspects is crucial for solving the SEEG-forward problem and hence simulating realistic SEEG signals. The present model includes a laminar architecture that represents the cortical layers of the human brain, following the framework in [30, 31] and also used in [32] for the analysis of seizure data. A laminar representation of a human cortical column of six layers with a physiological thickness [33] and uniform conductivity $\sigma = 0.3\times 10^$ S mm−1 [34] was considered. Given the anatomical and spatial characteristics of the pyramidal cells, we assumed that only the synaptic inputs on the pyramidal cells contribute to the local field potential (LFP) signal measured by a virtual SEEG electrode [35]. Synapses on the pyramidal cells were considered as point contacts as an approximation of the spatial distribution of the synaptic locations. They were placed along a one-dimensional fiber passing through cortical layers. Physiological constraints regarding the distribution of the synaptic contacts across the layers were loosened, and average synaptic locations were considered. Nevertheless, this simplification is appropriate given the level of the mesoscopic NMM formulation.

To a large extent, extra-cellular field potential (also known as LFP) is caused by the synaptic interactions between sub-populations of neurons and interneurons in a given neural population [35, 36]. Briefly, activation of an excitatory (inhibitory) synapse leads to a net inward (outward) flow of cations to the neuron, which depolarizes (hyperpolarizes) the membrane potential. Inflow (outflow) of positive charges into the neuron causes a negative (positive) current at the level of synaptic contact (so-called active sink (active source)). Meanwhile, due to the depolarizing (hyperpolarizing) membrane potential, further intra- and extra-cellular ion currents move along the membrane and a positive (negative) current will develop elsewhere along the neuronal element (so-called passive source (passive sink)). Then, there will be an extra-cellular current I(t) flowing from the source to the sink. The field potential measured by a contact with reference $V_ = 0$ is then the difference between the source and the sink,

Equation (3)

where r1 and r2 are the distance to the measurement point of the sink and of the source in a homogeneous infinite media, respectively. The amplitude I(t) equals the post-synaptic current at the level of the synapses, and it is computed from the PSP given by the solution to equation (2),

Equation (4)

where $\eta = 10^$ S is a conversion factor relating the PSP to post-synaptic current [31, 32]. In reality, η is a lumped parameter that reflects neuronal population characteristics, such as cell density and morphology, but these factors were ignored.

We did not model the electrode geometry but considered the voltage at the center of two contact points 2 mm away from each other. We assumed that the virtual SEEG electrode was parallel to the cortical column at a 10 mm distance from the considered neural mass. The system was simulated using the Euler–Murayama numerical method with a step size $dt = 10^$ s and the same noise vector following a normal distribution.

In the present section, results are reported following the methods presentation order.

3.1. SEEG analysis3.1.1. SWs mean waveform computation

The mean waveform is computed for each patient’s EZ and NEZ bipolar channels. This representation is used to qualitatively understand morphological similarities and differences among the SWs and highlight the presence of discriminating characteristics to identify the EZ. In figures 3(a) and (b), the mean waveform of the patient P1’s EZ and NEZ bipolar channels is represented in red; in blue are depicted all the SWs used for the mean waveform computation. The mean waveform was computed from an average of 6432 SWs per bipolar channel.

Figure 3. Detected spike-waves (SWs) and mean waveforms in patient P1. Detected SWs (blue) on the epileptogenic zone (a) and non-epileptogenic-zone (b) bipolar channels are used to compute the corresponding mean waveform (red). In (c) and (d) are reported the same plots in (a) and (b) but for a restrained number of SWs (n = 30) randomly selected over time: the green line is the mean SWs waveform computed on those 30 SWs and its shape is comparable to the red one computed respectively in (a) and (b). The plotted SWs were normalized, resulting in a dimensionless amplitude.

Standard image High-resolution image

To show the IEDs morphology is stationary over time, figures 3(c) and (d) display in green the mean SWs mean waveform computed on a subset of IEDs (n = 30) randomly extracted over time for the same bipolar channels. The green signal is similar to the red one depicted in figures 3(a) and (b): this visually confirms that SWs morphology does not vary so much and that an average is appropriate.

The computation of the SWs mean waveform reveals qualitative differences between interictal events generated in the EZ with respect to those generated in the NEZ. In particular, the FWHM Spike is less pronounced in the EZ and the time delay between spike and wave is shorter. Also, FWHM Wave appears larger in the NEZ than in the EZ.

3.1.2. SWs features distribution and statistical tests

For each patient, 200 SWs were used for features extraction: 100 to identify EZ SWs morphology and 100 to discover NEZ SWs characteristics. Features extracted from combined patients’ records show different distributions between the EZ and NEZ. Figure 4 summarizes these findings in violin plots (the red area refers to the EZ while the green one refers to the NEZ). Dotted horizontal lines represent the quartiles of the distribution. These results indicate that spike amplitude is lower in the EZ than in the NEZ, while the opposite is observed for wave amplitude. The SW Delay and FWHM Delay are shorter for the EZ SWs than the NEZ SWs. The FWHM Spike distribution demonstrates sharper spikes in the EZ SWs than in the NEZ SWs. Almost no difference exists in the FWHM wave feature in these two regions. The ratio features highlight a higher median value for FWHM Wave over FWHM Spike in the EZ, indicating that the EZ spikes have a sharper form than the NEZ spikes. The ratio of FWHM Wave over FWHM Delay reflects a smaller FWHM Delay in the EZ. Finally, the ratio Spike Amplitude/Wave Amplitude informs about a lower amplitude of spikes with respect to the Wave Amplitude in the EZ, but a higher Spike Amplitude with respect to the Wave Amplitude in the NEZ. This finding could be linked to the fact that a stronger inhibition effect could be found in the EZ due to the high excitability of the region. These general features are considered for the SWs modeling. Inter-patient differences still exist, and those are reported in appendix A.

Figure 4. Features distribution between epileptogenic zone (red) and non-epileptogenic zone (green). Dotted horizontal lines represent the quartiles of the distributions. For abbreviations refer to table 3. Data is normalized according to z-score; for details, refer to section 2.2.5.

Standard image High-resolution image

Table 4 reports the mean value of the features extracted from the SWs. In particular, the reported values correspond to the SWs located in the EZ and NEZ per patient; finally, the mean value of the features computed across patients in the EZ and NEZ groups is reported. From the table it is possible to draw some considerations: Wave Amplitude seems to be higher in the EZ (with respect to the NEZ) for all patients; SW Delay, FWHM Spike, FWHM Delay, and FWHM Wave resulted, instead, to be lower in the EZ for four patients over five.

Table 4. Mean values of features computed per patient in the epileptogenic zone (EZ) and non-epileptogenic zone (NEZ). The lowest feature mean value, between the EZ and NEZ, is highlighted in blue, the highest in orange. For abbreviations see table 3.

 P1P2P3P4P5All dataFeatureEZNEZEZNEZEZNEZEZNEZEZNEZEZNEZ Spike Amp3.5262.6033.1664.8173.1752.3592.0303.8583.1513.7483.0103.477Wave Amp2.8572.1322.3862.1142.9112.8953.4102.5093.2092.8502.9552.500SW Delay0.1570.2160.2840.2780.2610.2610.0940.1740.0890.1820.1770.222FWHM Spike0.0150.0430.0450.0610.0320.0350.0250.0230.0200.0350.0270.040FWHM Delay0.0600.0960.1440.1340.1180.1200.0450.0890.0330.0800.0800.104FWHM Wave0.1990.2500.3110.2380.2530.2650.1140.1440.1460.2010.2040.219FWHM W/FWHM S15.1616.6027.1363.9128.5569.4705.2926.4608.1955.7918.8686.447Spike Amp/Wave Amp1.4401.4221.6362.4521.1600.8790.6091.8671.1061.4851.1901.621FWHM W/FWHM D3.5982.7812.3492.0042.4852.4862.6411.9775.1032.5323.2352.356

The null hypothesis of the Kolmogorov Smirnov test (which assumes no difference between the observed and a normal distribution) is always rejected: p-value smaller than 0.05. The Wilcoxon signed-rank test, instead, shows that the EZ and the NEZ distribution differ for each feature, as the p-value is again smaller than 0.05.

3.1.3. SWs features classification

Overall, our results suggest that the EZ SWs present some subtle morphological differences with respect to those generated in a NEZ. We conducted a k-means classification using both two and three SWs extracted features from the selected patient recordings. The best results for the two features classification are accomplished by coupling the FWHM Delay and the FWHM Spike. The best three-element tuple candidate is composed using the following features: FWHM Delay, the spike FWHM, and spike amplitude. Performances of these classifications are reported in figure 5 and table 5.

留言 (0)

沒有登入
gif