Spatiotemporal Mapping of Interictal Spike Propagation: A Novel Methodology Applied to Pediatric Intracranial EEG Recordings

Abstract

Synchronized cortical activity is implicated in both normative cognitive functioning and many neurologic disorders. For epilepsy patients with intractable seizures, irregular synchronization within the epileptogenic zone (EZ) is believed to provide the network substrate through which seizures initiate and propagate. Mapping the EZ prior to epilepsy surgery is critical for detecting seizure networks in order to achieve postsurgical seizure control. However, automated techniques for characterizing epileptic networks have yet to gain traction in the clinical setting. Recent advances in signal processing and spike detection have made it possible to examine the spatiotemporal propagation of interictal spike discharges across the epileptic cortex. In this study, we present a novel methodology for detecting, extracting, and visualizing spike propagation and demonstrate its potential utility as a biomarker for the EZ. Eighteen presurgical intracranial EEG recordings were obtained from pediatric patients ultimately experiencing favorable (i.e., seizure-free, n = 9) or unfavorable (i.e., seizure-persistent, n = 9) surgical outcomes. Novel algorithms were applied to extract multichannel spike discharges and visualize their spatiotemporal propagation. Quantitative analysis of spike propagation was performed using trajectory clustering and spatial autocorrelation techniques. Comparison of interictal propagation patterns revealed an increase in trajectory organization (i.e., spatial autocorrelation) among Sz-Free patients compared with Sz-Persist patients. The pathophysiological basis and clinical implications of these findings are considered.

Keywords: epilepsy surgery; epileptogenic zone; interictal spike propagation; pediatric epilepsy; surgical outcome.

Figure 1

Spike discharges were detected and encoded in spike frequency maps. (A) For each patient, interictal spike discharges (red asterisk) were identified using an automated MATLAB detector. Spikes were identified based on characteristic morphological features and patient-specific thresholds. (B) Raster plot demonstrates the spike detector output from one representative 60-min EEG segment (Patient 18). Black dots correspond to the channel and time of each detected spike. The raster is sorted by spike count. (C) Left: spikes are unevenly distributed across channels. For illustrative purposes, the Lorenz curve (dotted) is used to describe the uniformity of the spike distribution. The Lorenz curve quantifies: “X percent of channels account for Y percent of the total spikes.” Deviation from the “Equality Line” (solid) confirms that the spike distribution is not uniform. Right: spike frequency maps were calculated by dividing each channel’s spike count by the total analyzed duration (spikes per minute).

Figure 2

Interictal spikes are sometimes observed at multiple electrode sites with discernible temporal latency. Left: electrode schematic for Patient 12, a 7-year-old male with a temporal lobe implant. The neurologist-defined seizure onset zone is shown in red. Right: 1-s IEEG epoch containing a multichannel interictal spike discharge. Detailed inspection of discharge reveals subtle latency differences within the spike sequence. Here, the spike on channel 6 (red) precedes the spike on channel 47 (blue). Dotted gray line marks the peak time of the lead spike (red channel).

Figure 3

Outlier spike sequences were eliminated using a sequence “cleaning” procedure. (A) Left: four spike trajectories (labeled 1–4) are spatially mapped on the 2D electrode rendering. Right: a trajectory comparison algorithm was used to compare pairs of spike trajectories. Similarity scores were based on the extent of spatiotemporal overlap between trajectory pairs, ranging from min = 0 to max = 1. (B) Pairwise similarity scores were encoded in a sequence-by-sequence similarity matrix, S. The inset (left) shows the similarity mapping between trajectories 1 and 4 from part A. As expected, similarity scores between sequences 1 and 2 are high while sequences 3 and 4 share considerable spatiotemporal overlap. Bottom left: the degree distribution of matrix S was calculated in order to identify outlier sequences. A standard k-means clustering algorithm with fixed cluster count (k) of k = 3 was used to classify sequences as “Low Degree,” “Medium Degree,” and “High Degree.” Low Degree sequences shared minimal spatiotemporal overlap with the rest of the sequence set and were discarded as suspected outliers.

Figure 4

Multichannel spike sequences were extracted using automated techniques. (A) A 2.5 s epoch (Patient 17) containing a multichannel spike sequence (dotted box, nine channels shown). The time of the first peak in the sequence (red channel) is shown as a solid vertical line (B) Left: spike overlay (duration = 625 ms) reveals latency differences between spikes in the sequence. Here, channel 17 (red) is the “leader” (i.e., first spike) of the multichannel discharge. Channel 28 (blue) peaks considerably later in the sequence. Middle: cumulative probability distribution (CPD) encodes the tendency for each channel to occur at a given latency across all spike sequences. The heterogeneity of the CPD demonstrates the preference for channels to appear at different recruitment latencies. Right: the mean recruitment latency was computed for each channel (red = early recruitment, blue = late recruitment). Again, channel 17 tends to be recruited earlier in spike sequences than channel 28. Spatial organization of recruitment latency maps was characterized using the Moran Index. For this patient, channels appearing early in spike discharges (warm colors, “source” regions) cluster together while “sink” regions (cool colors) cluster together, resulting in a high Moran Index.

Figure 5

Spatial organization of heat maps was quantified using the Moran Index. Spike frequency maps and recruitment latency maps were characterized using the Moran Index, a spatial autocorrelation technique. Moran Indices ranged from −1 (perfect spatial anti-autocorrelation) to +1 (perfect spatial autocorrelation), with 0 corresponding to no patterns of spatial autocorrelation. To illustrate the Moran Index technique, simulated data are used to construct spatial maps (electrodes = 128) of varying degrees of spatial organization.

Figure 6

Recruitment latency maps (but not spike frequency maps) differentiate Sz-Free and Sz-Persist groups. The spatial organization of spike frequency maps and recruitment latency maps were assessed using the Moran Index. Box plots are used to display group results (box = 25th–75th percentile, horizontal line = median, whiskers = 5th–95th percentile). The spatial organization of spike frequency maps (top) does not differentiate clinical groups (Wilcoxon rank sum test, p = 0.863). When recruitment latency maps were examined (lower), we observed significantly increased spatial organization (Moran Index) among Sz-Free (0.447 ± 0.160) compared with Sz-Persist (0.275 ± 0.088) patients (Wilcoxon rank sum test, **p = 0.003). n.s., not significant.