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

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.

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.

INtRodUCtIoN
Interictal spikes are transient, high-amplitude discharges resolvable using both non-invasive (i.e., scalp) and invasive EEG electrodes. The presence of EEG spikes is strongly associated with epilepsy and contributes to the diagnosis; about 90% of epilepsy patients show EEG spiking compared to <1% of people without epilepsy (1). In the presurgical evaluation of intractable epilepsy patients, neurologists traditionally interpret spikes as focal biomarkers of pathological tissue that should be resected to ensure seizure freedom (2). However, evidence supporting this interpretation of spike activity is inconclusive (3). Resection of the spiking region has been associated with improved surgical outcomes (4,5), though for many patients, persistent postoperative spiking is compatible with seizure freedom (6). Spikes are routinely observed in regions beyond the neurologist-defined seizure onset zone (SOZ), and the localization of spikes to the SOZ varies tremendously between patients and across time (7). Although spikes have been implicated in the formation of epileptic foci (8)(9)(10)(11), several groups have suggested that periods of increased spikes actually protect against seizure occurrences (12,13), challenging the assumption that "spikes beget seizures" (14).
Interictal spikes are sometimes recorded at multiple sites with discernible latency, suggesting that spikes can propagate rapidly across the cortex (15)(16)(17)(18). In recent years, several groups have used spike propagation to visualize patterns of connectivity between remote neural populations (18,19). This approach aligns with the emerging view of the epileptogenic zone (EZ) as a densely interconnected network capable of generating and sustaining seizures (20,21). By mapping spike propagation through the cortex, several studies have correlated propagation trajectories with important clinical variables such as epileptic pathology (22) and SOZ localization (23)(24)(25)(26)(27). However, previous studies have been limited by important methodological concerns, e.g., the construction and overinterpretation of small, visually selected prototype spike datasets. Such methods require a priori information about spike events and have limited use in the surgical setting. While evidence is mounting for the clinical relevance of spike propagation, further study and methodological refinement are needed.
To assess the role of spike propagation in the presurgical evaluation, we developed a novel methodology for detecting, visualizing, and characterizing spike trajectories (or sequences) and applied it to a sample of 18 pediatric intracranial EEG recordings. After constructing large, unbiased spike datasets, we tested the hypothesis that the spatial organization of spike trajectories would reliably differ between patients with favorable (i.e., seizure-free) versus unfavorable (i.e., seizure-persistent) surgical outcomes. Additionally, we estimated the clinical impact of the methodology by comparing our approach to a more traditional spike analysis (i.e., mapping the focal density of spikes rather than their spatiotemporal spread). Our results suggest that patients with seizure-free outcomes exhibit more spatially organized interictal propagation patterns than patients with recurrent postoperative seizures. This work advances our understanding of propagating spike discharges and identifies a potential role for spike propagation as a quantitative surgical candidacy biomarker.

Patient selection
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 postoperation = 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.

eeG Acquisition
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. (i.e., activity between EEC and UEO) were clipped and discarded, leaving only the non-seizure periods for analysis.

Automated spike detection
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).

Interictal segmentation
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-ofday, proximity to nearest seizure, and physiological state.

Construction of spike Frequency Maps
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).

Identifying Multichannel spike sequences
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.

Applying Constraints to spike sequences
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 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). 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.

eliminating outlier spike trajectories
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  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. 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 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 (ssthresh and ttthresh, 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 (dmatch) from the reference point was used to compute the following score: simpoint = 1 − dmatch/ssthresh, which ranged from min = 0 to max = 1. When no "matches" were found in the test sequence, simpoint = 0. This procedure was repeated for the remaining points in the reference sequence, and the average of the set of simpoint 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-bysequence 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.

Construction of Recruitment Latency Maps
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 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. from "source" regions with lower recruitment latencies (warmer colors) to "sink" regions with higher recruitment latencies (cooler colors).

Assessing spatial organization of spike Frequency and Recruitment Latency Maps
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: where N is the number of channels, Fi is the spike frequency of channel i (spikes per minute), Favg is the mean spike frequency across all channels, and wij is the spatial "weight" between channels i and j (29). For recruitment latency maps, the same equation was used, but Fi corresponded to the mean recruitment latency of channel i (ms) and Favg was the mean recruitment latency across channels. In this study, if the Euclidean distance (d) between channels i and j (dij) was ≤1.5 cm (i.e., the channels were immediately adjacent to one another, including diagonally), then wij was set to 1/dij. When dij exceeded 1.5 cm, wij 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).

Comparing spike Frequency and Recruitment Latency Maps across Clinical Groups
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.

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

ResULts
The objective of this study was to assess the clinical utility of spike propagation as a biomarker for EZ organization and eventual seizure outcome in pediatric epilepsy surgery. To accomplish this goal, we developed a novel methodology for detecting, clustering, and visualizing global spike patterns from presurgical IEEG recordings. Then, we conducted a series of spike analyses on 18 patients (mean age = 10.9 ± 4.8 years) divided into two outcome groups: Sz-Free (Engel Score = 1, n = 9) and Sz-Persist (Engel Score, n = 9).

Patient Characteristics
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).

eeG extraction and spike detection
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).
spatial organization of spike Frequency Maps does Not differ by surgical outcome 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.

spatial organization of Recruitment Latency Maps Is Increased Among Patients with Favorable 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 tABLe 2 | the spatial density ("spike frequency analysis") and spatiotemporal propagation ("recruitment latency analysis") of interictal spike discharges were examined. 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).

dIsCUssIoN
Defining and deciphering the network structure of the EZ is critical for improving the surgical management of intractable seizures. In this study, we examined connectivity patterns within the epileptic brain by mapping the propagation of intracranially recorded spike discharges in pediatric epilepsy surgery patients. Using a spatial autocorrelation approach, we found that the global organization of spike trajectories reliably distinguished patients with favorable versus unfavorable surgical outcomes. Further, we found that group differences were not observed using a more traditional (i.e., density-based) approach to spike activity, suggesting that our technique elucidates group differences that are inaccessible to spatial density analyses.

Methodological Advancements
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)(38)(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.
spike Propagation Is More spatially organized in Patients with Favorable outcomes 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.
spike Propagation Reveals Group differences that Are Not observed in the More traditional, density-Based spike Analysis 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.

Future directions and Limitations
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.

Conclusion
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.

AUthoR CoNtRIBUtIoNs
The following contributions were made by the authors: experimental hypothesis and study design (ST, EM, and BP), data extraction and curation (ST, CB, CC, MB, and EM), data analysis (ST and CB), manuscript preparation (ST), manuscript editing (ST, EM, and BP), and final manuscript approval (EM).