The Children’s Hospital of Philadelphia (CHOP) Institutional Review Board approved this study. Each patient’s legal guardian signed written consent in accordance with the Declaration of Helsinki. Full-duration IEEG recordings were obtained from a database of Phase II presurgical evaluations performed at our institution between the years of 2002 and 2009. Each patient required intracranial EEG monitoring following unsatisfactory non-invasive localization of the epileptic foci. Of the 30 patients available for study, 18 patients (12 male, 6 female, mean age = 10.9 years, range = 3–20 years) met the inclusion criteria: (1) availability of detailed intraoperative photos; (2) unambiguous seizure markings; and (3) availability of ≥24 h of recordings. Patients were not screened for a particular clinical history, seizure semiology, seizure onset location, or electrographic onset pattern. All implants were clinically motivated without influence from this study. Retrospective review of patient charts provided detailed information regarding implantation site, pathology, etiology, and MRI description. Postsurgical outcomes were assessed by primary neurologists upon last patient contact (minimum post-operation = 2 years) and classified using Engel’s modified scale (28): Class 1 = seizure free; Class 2 = significant improvement; Class 3 = worthwhile improvement; and Class 4 = no improvement. Patients with complete postsurgical seizure freedom (Engel Score = 1, n = 9) were assigned to the “Sz-Free” group while patients experiencing persistent seizures (Engel Score ≥ 2, n = 9) were assigned to the “Sz-Persist” group.

Intracranial recordings were obtained using a Telefactor Beehive 32-128 channel CTE digital synchronized video-EEG system with 16-bit amplifiers and 200 Hz sampling rate (Astro Med Corp., West Warwick, RI, USA). Subdural grids and strips with embedded platinum electrodes (exposure diameter = 2.3 mm, interelectrode distance = 10 mm) were used for all patients (Adtech Medical Instrument Corporations, Racine, WI, USA). An outward-facing epidural electrode strip served as the online electrical reference. Data were subjected online to an analog antialiasing notch filter (60 Hz). De-identified intraoperative photos and diagrams provided by the surgical staff were used to construct 2D electrode maps for each patient. Cartesian coordinates were assigned to each channel using the interelectrode distance of 10 mm. Distances separating adjacent grids and strips were estimated using convenient anatomical landmarks. Electrode mapping was performed by an author (Samuel B. Tomlinson) blind to patient identity and confirmed by the senior author (Eric D. Marsh) prior to analysis. We acknowledge that this 2D approach does not capture the precise physical distance between electrode pairs, which is influenced by the 3D structure of the brain and the sulcal/gyral patterns beneath the arrays (29). However, approximating spatial relationships among channels using 2D schematics is common practice in the literature when coregistered MRI and CT imaging studies are unavailable.

Two experienced pediatric epileptologists (Brenda E. Porter and Eric D. Marsh) independently inspected each patient’s full IEEG records in referential montage and identified all seizures, marking the following parameters: (i) time of earliest electrical change (EEC) at ictal onset; (ii) time of unequivocal electrographic offset (UEO); and (iii) the seizure onset channels. To isolate interictal segments from full-duration recordings, seizure intervals (i.e., activity between EEC and UEO) were clipped and discarded, leaving only the non-seizure periods for analysis.

In order to detect interictal spikes, we utilized an automated spike detector developed previously by our group and validated against human markings (30). The spike detection algorithm and the procedure for determining patient-specific thresholds are described thoroughly in the original publication. Briefly, spikes are identified based on stereotypic morphological features such as peak amplitude, spike duration, and the presence of a characteristic post-spike slow-wave (7). A previous validation study demonstrated that the spike detector performs comparably to expert epileptologists, with human verification of detected spikes exceeding an average of 75% (30). The detector is designed to process long EEG segments and output a two-column matrix encoding the channel index (column 1) and time (column 2) corresponding to each detected spike. For each patient, full interictal recordings were submitted to the spike detector. Spike detection was performed using MATLAB 2016a (Mathworks, Natick, MA, USA).

Across patients, the duration of interictal EEG and the number of detected spikes were highly variable. In order to construct consistent spike datasets amenable to cross-subject comparisons, we divided each patient’s spike detector output into non-overlapping segments each containing 10,000 interictal spikes. Then, for each patient, 10 segments were randomly selected and concatenated to yield the final interictal spike dataset (spikes/patient = 100,000). Segments were chosen without regard for factors such as time-of-day, proximity to nearest seizure, and physiological state.

First, we aimed to characterize the spatial distribution of spikes across the recording field (i.e., assessing how frequently spikes occurred at each electrode). Analyzing spikes in this manner closely parallels the traditional clinical interpretation of spike activity, with regions of high spike frequency typically marked as important areas for resection (2). For each patient, interictal spike frequency maps were generated by dividing each channel’s spike count by the total analyzed EEG duration (mins) (Figure 1). Frequency maps were coded with warmer colors representing higher spike frequency (spikes per minute) and cooler colors representing lower spike frequency (Figure 1C, right). In addition to spike frequency maps, spike distributions were represented using the Lorenz curve (31), which visualizes the distribution of assets (spikes) across the population (electrodes) (Figure 1C, left). The Lorenz curve was quantified using the Gini coefficient, which encodes the difference between the Lorenz curve and the linear “Equality Line” (i.e., a uniform distribution in which all channels exhibit the same percentage of the total spikes).

After visualizing the spatial density of mono-channel spike discharges, we next aimed to characterize the spatiotemporal propagation of spikes through the recording field. To do so, we developed a novel MATLAB algorithm that extracts multichannel sequences from the spike detector output (Figure 2). The algorithm initializes by defining the first detected spike as the “leader” of a candidate spike sequence (23). Successive spikes occurring within 50 ms of the leader were then appended to the candidate sequence. Fifty milliseconds was chosen based on previous studies showing that interictal spikes consistently propagate with 10–50 ms latencies (15, 23, 25, 32). Once a spike outside of 50 ms from the leader was encountered, this spike time was compared to the previous spike in the sequence. If the latency between the two events was ≤15 ms, the spike was added to the sequence. This additional parameter afforded the temporal flexibility needed to capture unexpectedly long spike sequences. A “terminating spike” was designated when the algorithm encountered a spike violating both temporal conditions (i.e., occurring ≥50 ms from leader and ≥15 ms from the previous spike). The terminating spike was then initialized as the leader of a new potential sequence, and the procedure restarted, continuing in this way until all spikes in the dataset were processed. To scrutinize the performance of the algorithm, two authors (Samuel B. Tomlinson and Eric D. Marsh) visually inspected the spike sequences extracted from random 60-min EEG segments for five patients. This preliminary review suggested that the majority of true-positive sequences included at least five mono-channel spikes while most false-positive detections included fewer than five spikes (data not shown). Therefore, only sequences containing ≥5 spikes were preserved for analysis.

Two constraints were applied to the extracted spike sequences. First, due to the modest acquisition rate used in this study (200 Hz), spikes were sometimes detected with the same peak time. By default, the spike detector ordered these “tied” events according to increasing channel number, which does not reflect the actual path of the trajectory. To address this issue, tied spikes were rearranged based on spatial proximity (Euclidian distance) to the nearest non-tied event under the assumption that the probability of spike transmission decreases with distance. Second, propagation velocity is presumably limited by the neural conduction speed and the synaptic delay. It is therefore possible that spikes with minimal latency differences appearing at very remote sites belong to separate (but temporally overlapping) spike sequences. To account for this possibility, channels were grouped into partitions containing about four to eight neighboring electrodes. Each successive spike in a given sequence was required to occur within the same partition or a partition immediately adjacent to the previous spike. This procedure minimized the grouping of distinct spike sequences and limited the intrusion of unrelated but frequently spiking channels into the sequences. Because long-range neural connections are possible, a condition was included that allowed for connections traversing multiple partitions, if the connection was frequent. To identify frequent connections, we constructed a channel-by-channel matrix, C, such that C(i,j) encoded the number of times that a spike propagated from channel i to channel j. Then, when assessing whether the (i,j) connection was “frequent,” we divided C(i,j) by the sum of matrix row i. If this fraction exceeded 0.05 (i.e., when a spike is transmitted from channel i, it propagates directly to channel j at least 5% of the time), then the C(i,j) connection was deemed “frequent” and the spike was preserved. Partition assignments were made via inspection of de-identified intraoperative photos by Samuel B. Tomlinson and confirmed by Eric D. Marsh prior to the study.

Identifying EEG spikes is an imperfect process, even for experienced human reviewers (7, 30). In this study, false-positive spike detections frequently occurred during periods of EEG artifact (e.g., during motor movements or ambient electrical noise), which were often falsely identified as propagating sequences. In order to discard erroneous sequence detections in an efficient and unbiased manner, we developed an unsupervised sequence “cleaning” procedure using a recent trajectory clustering algorithm (33). Briefly, the clustering algorithm computes a similarity score (min = 0, max = 1) between pairs of spike sequences based on their spatiotemporal overlap (Figure 3). For each patient, we identified all interictal spike sequences and submitted them to the comparison algorithm. Spike sequences were represented as a series of points, each taking the form (x-location, y-location, latency from lead spike). Then, a similarity score was computed between each pair of sequences in the set according to the following procedure. First, one sequence was designated as the “reference” sequence and the other was designated as the “test” sequence (Figure 3A, right). For each point in the reference sequence, the algorithm searched for data points (or “matches”) in the test sequence that satisfied predefined spatial and temporal thresholds (ss_thresh and tt_thresh, respectively). Specifically, the algorithm found all points in the test sequence that (i) fell within 1.5 cm (Euclidean distance) from the reference point, and (ii) occurred within 15 ms of the reference point (i.e., occurred with similar latency from lead spike). Then, the “matched” point with the shortest Euclidean distance (d_match) from the reference point was used to compute the following score: sim_point = 1 − d_match/ss_thresh, which ranged from min = 0 to max = 1. When no “matches” were found in the test sequence, sim_point = 0. This procedure was repeated for the remaining points in the reference sequence, and the average of the set of sim_point values constituted the similarity score between the reference and test sequences (again, min = 0 and max = 1). By iterating this process through all possible sequence pairs in the dataset, we constructed a non-symmetrical sequence-by-sequence similarity matrix, S, in which the S(i,j) cell contained the similarity score between sequences i and j (Figure 3B, right). Using this matrix, we computed the degree centrality (34) of each sequence to identify sequences sharing minimal spatiotemporal overlap with other discharges in the dataset. Based on the degree distribution (Figure 3B, bottom left), sequences were classified as “Low Degree,” “Mid Degree,” or “High Degree” using the MATLAB k-means clustering algorithm with a fixed cluster (k) count of k = 3, and all “Low Degree” sequences were eliminated as suspected outliers.

After extracting multichannel spike sequences, we aimed to characterize global patterns of spike propagation. To do so, we created “recruitment latency” maps for each patient (Figure 4). Recruitment latency maps were constructed by determining each electrode’s average temporal position (ms) within spike sequences (i.e., the average latency from the first spike in the discharge). This information was first represented as a cumulative probability distribution (CPD) encoding the tendency for each channel to occur at a given latency across all spike sequences (Figure 4B, middle). Channels that consistently appeared early in spike sequences thus exhibited a left-shifted CPD. Then, to visualize recruitment latencies across electrodes, recruitment latency maps were constructed (Figure 4B, right), with warmer colors representing shorter mean recruitment times (29). This approach allowed us to visualize the global movement of spikes from “source” regions with lower recruitment latencies (warmer colors) to “sink” regions with higher recruitment latencies (cooler colors).

As described in the previous sections, two spatial maps were constructed per patient in order to describe interictal spike behavior: the spike frequency map (see Construction of Spike Frequency Maps), and the recruitment latency map (see Construction of Recruitment Latency Maps). Using these distinct approaches, we aimed to characterize the spatial organization of spike patterns. To examine the spatial organization of spike frequency maps, we asked: do regions of high (low) spike frequency cluster together in space? Similarly, for the recruitment latency maps, we assessed: do regions appearing early (late) in propagating discharges cluster together in space? Both questions could be readily examined using the Moran Index (35), a spatial autocorrelation measure that determines whether channels close in spatial proximity exhibit similar patterns of behavior (or in this case, similar spike frequencies/recruitment latencies). For the spike frequency maps, the Moran Index (I) was computed using the formula: I=N∑i∑jwij∑i∑jwij(Fi−Favg)(Fj−Favg)∑i(Fi−Favg)2,

where N is the number of channels, F_i is the spike frequency of channel i (spikes per minute), F_avg is the mean spike frequency across all channels, and w_ij is the spatial “weight” between channels i and j (29). For recruitment latency maps, the same equation was used, but F_i corresponded to the mean recruitment latency of channel i (ms) and F_avg was the mean recruitment latency across channels. In this study, if the Euclidean distance (d) between channels i and j (d_ij) was ≤1.5 cm (i.e., the channels were immediately adjacent to one another, including diagonally), then w_ij was set to 1/d_ij. When d_ij exceeded 1.5 cm, w_ij was set to 0. An illustration of the Moran Index applied to various spatial maps is provided in Figure 5 using simulated data. As the simulation demonstrates the Moran Index ranges from −1 (perfect anti-autocorrelation between neighboring channels) to +1 (perfect autocorrelation between neighboring channels), with 0 indicating no pattern of spatial autocorrelation (29).

Approaching spike activity from a spatial density perspective (i.e., spike frequency maps) versus a spatiotemporal trajectory perspective (i.e., recruitment latency maps) leads us to ask: does either approach provide insights into the likely surgical outcome of the patient? This question follows from our original hypothesis that spike propagation would reliably differentiate seizure-free and seizure-persistent patients whereas the more traditional spike frequency analysis would not capture such differences. To assess whether spike measures were related to surgical outcome, we conducted two comparisons between the study’s outcome groups (i.e., Sz-Free, n = 9; Sz-Persist, n = 9). In the first comparison, Moran Indices corresponding to spike frequency maps were compared across groups. In the second comparison, Moran Indices from the recruitment latency maps were compared across groups. Group comparisons were performed using the non-parametric Wilcoxon rank sum test with a Bonferroni-adjusted significance threshold of p = 0.05/2 = 0.025.

Unless otherwise noted, group-level statistics are reported as mean (μ) ± SD (σ), and p-values correspond to Wilcoxon rank sum tests.

Clinical data from the patient cohort are presented in Table 1. Across patients, surgical implants ranged from 64 to 126 electrodes and did not differ in size between Sz-Free (102.4 ± 16.7 electrodes) and Sz-Persist (100.1 ± 16.2 electrodes) patients (p = 0.779). Implant locations were highly variable across patients and did not noticeably differ between groups. The percentage of channels included in the SOZ varied across patients (range = 4.0–100%) but did not differ between outcome groups (p = 0.436). Other variables including average follow-up duration, years since surgery, seizure frequency, and implant duration did not differentiate groups (data not shown).

Seizures were clipped from full-duration EEG datasets, leaving only the non-seizure windows for analysis (Table 2). For each patient, 10 segments of 10,000 interictal spikes were randomly extracted from the interictal spike detector output. Across patients, the duration of EEG (minutes) encompassed in the spike dataset varied dramatically but did not differ between Sz-Free (804.5 ± 702.8 min) and Sz-Persist (529.8 ± 262.8 min) groups (p = 0.667). Note from Table 1 that two patients (Patient 01, Sz-Free; Patient 08, Sz-Persist) had only five available interictal segments, yielding 50,000 spikes for analysis. Reanalyzing the data using five segments for each patient did not appreciably change the results presented below (Figure S1 in Supplementary Material).

Interictal spike frequency maps (spikes per minute) were constructed for each patient to visualize the distribution of spike events across channels. The spatial organization of spike frequency maps (quantified by the Moran Index) varied substantially across patients (Table 2). When surgical outcome groups were compared (Figure 6, top), Moran Indices (min = −1, max = 1) were similar for Sz-Free (0.407 ± 0.161) and Sz-Persist (0.389 ± 0.136) groups (p = 0.863). This result suggests that the spatial organization of spike frequency maps was not associated with surgical outcome.

We applied our novel spike propagation trajectory methods to the interictal detector output. Consistent with previous studies (17, 18), the frequency of multichannel spike sequences (sequences per minute) varied tremendously across patients (Table 2). The average sequence frequency across patients was 4.32 ± 6.10 sequences/min, with no differences between Sz-Free (5.59 ± 8.39 sequences/min) and Sz-Persist (3.05 ± 2.23 sequences/min) groups (p = 0.863). To characterize the global flow of propagating spikes across the recording field, recruitment latency maps were constructed by determining each channel’s average temporal position within spike sequences. When recruitment latency maps were compared across outcome groups (Figure 6, bottom), we observed significantly increased spatial organization (Moran Index) among Sz-Free (0.447 ± 0.160) compared with Sz-Persist (0.275 ± 0.088) patients (p = 0.003).

The novel methodology introduced here immediately advances the study of interictal spike propagation in several important ways. First, our method is the only published study that combines a logical automated spike propagation algorithm with a sophisticated trajectory “cleaning” procedure, allowing us to examine spike trajectories across the entire interictal recording without introducing overwhelming artifact to the dataset. In stark contrast, most previous studies have limited their analyses to brief, manually selected EEG segments with prototypic spike activity, which prohibits analysis of spike patterns across the full recording period and leads to sampling biases. Indeed, recent work finds that spikes originating from the same lead region follow highly variable trajectories (22), suggesting that limited samples of spikes may provide misleading information about the overall propagation patterns. Next, this study is among the first to examine the relationship between presurgical spike propagation and postsurgical seizure outcome. Demonstrating that spike propagation analyses can disentangle surgical outcome groups is an important step for translating this paradigm to the clinical setting. Additionally, although several groups have developed source reconstruction models of propagating spikes in EEG (17, 25, 27, 32, 36) and MEG recordings (37–39), few have attempted to capture the state of the entire recording field with a single global metric. By applying the Moran Index to the recruitment latency maps constructed for each patient, we reduced the dimensionality of the multivariate electrode grid to a global statistic with an intuitive interpretation (i.e., spatial organization). This approach is advantageous compared to more complex analyses because it may allow clinicians to efficiently incorporate spike propagation into the workflow of the presurgical evaluation. Finally, our study is unique in demonstrating the superiority of the spike propagation approach over a more traditional spike analysis based on spatial density. While the methodology of this study immediately advances the field in these regards, the clinical efficacy of our technique must be validated in a larger study with prospectively collected IEEG data from both adult and pediatric patients.

In our pediatric cohort, patients with seizure-free outcomes exhibited increased spatial organization of spike trajectories relative to patients with persistent seizures. This finding aligns with recent work by Martinet et al. (29), who mapped the spatiotemporal propagation of neocortical seizures and demonstrated that more consistent and organized ictal recruitment patterns were associated with improved surgical outcomes. Critically, these authors relied on a similar application of the Moran Index to analyze ictal recruitment latency maps, maximizing the comparability of the findings. The concordance of these studies is somewhat surprising given that the relationship between spikes, seizures, and surgical outcome is highly inconsistent across studies (10). However, researchers have yet to determine whether spike propagation and seizure propagation networks overlap in space, which may reconcile the findings of these two studies. Critically, the analysis of spike propagation has practical benefits over seizure propagation, as interictal spikes are more readily acquired and may reduce the necessary implant time in patients undergoing Phase II surgical evaluation. Nonetheless, the findings bolster one another in suggesting that epileptiform propagation elucidates patterns of EZ connectivity that may offer insights into the eventual surgical outcome of the patient.

We propose that spike propagation analyses differ from more traditional spike density analyses because the former captures the network structure and organization of the EZ in a manner inaccessible to the latter. Indeed, in recent years, network analysis of EEG recordings has emerged as a major area of research (40, 41). Researchers have used principles from network science and graph theory to localize the EZ (42, 43), examine the events surrounding ictal onset (44) and predict surgical outcome (45) with increasing levels of success. Spike propagation is an intuitive method for visualizing connectivity patterns within the EZ, allowing researchers to quantify both the strength and directionality of pairwise connections (i.e., how often is a spike transmitted from channel A to B, and vice versa). The results from this study suggest that spike propagation is a valuable and potentially under-utilized tool for assessing the network configuration of the EZ.

This study provides evidence for the relevance of spike propagation in the surgical management of epilepsy patients. Although the methodological and clinical advances of this study are significant, future steps must be taken to bolster our results. First, a major limitation of this study is the small number of patients included in the analyses. While small samples are common in both adult and pediatric IEEG studies, validating the clinical findings of our study requires a prospective study with considerably larger group sizes. Second, the reliance on two-dimensional spatial electrode maps is a potential limitation of the design. While recent procedures have been introduced to localize electrodes using coregistered presurgical CT imaging, many IEEG studies continue to rely on 2D electrode schematics. We posit that representing the electrode array in two dimensions is a reasonable simplification, given our interest in mapping spike propagation across the cortical surface. However, future studies using stereo-EEG electrodes or synthetic depth electrodes (MEG) may reveal important contributions from deep neural sources, and validating electrode positions using co-registration technologies may strengthen the present results. Finally, as with all IEEG studies, the study is limited by the clinically chosen electrode locations and may not sample the full extent of cortical spiking.

This study presents a novel methodology for detecting, extracting, and characterizing presurgical spike propagation patterns. Our results demonstrate the potential clinical utility of spike propagation in assessing surgical candidacy and characterizing the EZ in pediatric epilepsy surgery. The findings of this study represent a positive step toward a presurgical biomarker of seizure freedom and highlight the need for improved understanding of interictal spike activity in the clinical work-up for epilepsy surgery. This work justifies a large-scale prospective study that will further elucidate the role of spike propagation in the surgical management of refractory neocortical seizures.