Bursts with High and Low Load of Epileptiform Spikes Show Context-Dependent Correlations in Epileptic Mice

Visual Abstract


Introduction
Epilepsy is a disease of hypersynchronized brain activity that can be monitored electrophysiologically as epileptiform activity (EA) in local field potentials (LFP). Aside from EA accompanied by behavioral symptoms, e.g., as in clinical seizures, the EA patterns that have received the most attention thus far are epileptiform spikes and compound events termed "electrographic seizures." The exact definition and even the existence of purely electrographic seizures, however, is being debated. There seems to be no clear boundary between electrographic seizures and other compound events, which can differ widely in spike rate and duration (Gotman, 2011). Although the evolution of spikes within larger bursts has been analyzed (Chauvière et al., 2012), a systematic quantitative time series analysis of compound events is lacking. In fact, clustering of seizures in patients (Osorio et al., 2009;Karoly et al., 2017) suggests that the incidence of other compound events might also fluctuate over time. This raises the questions whether information about the state of epileptic networks could be read from the time series of compound events and whether the occurrence of one type of EA pattern influences the occurrence of other patterns. The aim of this study was to identify different types of compound events, investigate their temporal succession and analyze interactions between different types of events in an animal model for mesial temporal lobe epilepsy (MTLE).
In MTLE, the most common form of focal epilepsy, partial seizures are frequent, while generalized convulsive seizures are rare (Engel, 2001). The intrahippocampal kainate mouse model is a well-accepted animal model to investigate MTLE because it exhibits several phenomena similar to the human condition, e.g., spatially localized EA and histologic changes (Riban et al., 2002;Guillemain et al., 2012;Häussler et al., 2012;Krook-Magnuson et al., 2013;Twele et al., 2016;Janz et al., 2017). While large behavioral seizures are as rare in this as in other rodent models of epilepsy (Cavalheiro et al., 1991;Rattka et al., 2013), events that have been labeled "electrographic seizures," "hyperparoxysmal discharges," or "high voltage sharp waves" recur particularly frequently in kainateinjected mice (Klee et al., 2017;Kilias et al., 2018). As detailed in Twele et al. (2016), previous studies using this model identified and classified these large EA events based on heuristic criteria and experts applied these criteria with highly divergent results. The same study demonstrated that even slight changes in event definition can have a critical impact on whether a candidate drug is found to be anti-epileptic or not. This highlights the need for a more data-driven and reproducible approach to identify different patterns of EA and their interaction during epileptiform dynamics.
Interactions between EA event types have mostly been addressed based on the incidence of individual epileptiform spikes and seizures. While specific sub-types of spikes are thought to be involved in seizure onset (De Curtis and Avoli, 2016), several studies interpreted interictal spiking as preventing seizure generation (Engel and Ackermann, 1980;Librizzi and De Curtis, 2003;Goncharova et al., 2016). Other reports proposed either that changes in interictal spike rate are governed by seizure dynamics (Jensen and Yaari, 1988;Gotman, 1991) or found the two phenomena to be independent from each other (Selvitelli et al., 2010). Recently, Chang et al. (2018) suggested that epileptiform spikes stabilize network activity in states of low excitability while the same events precipitate seizures if they occur during transition phases of excitability. Hence, temporal context should be considered when interpreting the role of EA events.
In this study, we systematically describe patterns of EA in the kainate mouse model of MTLE using a data-driven approach and investigate how such patterns interact. Motivated by the observation that spikes preferentially occurred in bursts, we developed a method that automatically transformed hippocampal LFPs into a time series of solitary spikes and bursts. These bursts were classified using machine learning techniques according to a composite measure that could best be summarized as spike load. In these time series, bursts of different categories aggregated and thus formed phases reflecting slow fluctuations of EA dynamics. We delimited these phases and found that the relationship between low-load and highload bursting depended on temporal context. Our results suggest that the burst structure of EA plays an important role in epileptiform dynamics, that different burst patterns of EA interact and that these interactions vary over time.

Animals and kainate injection
Data originated from 18 adult male C57BL/6N mice (age 9 -12 weeks, Charles River, Sulzfeld, Germany) that were part of previous studies and in which a detailed description of the surgical and electrophysiological procedures can be found (dataset A with N ϭ 7 from Janz et al., 2017;dataset B with N ϭ 12 from Froriep et al., 2012). Briefly, mice were anesthetized (ketamine hydrochloride 100 mg/kg, xylazine 5 mg/kg, atropine 0.1 mg/kg body weight, i.p.) and kainate solution (50 nl, 20 mM in 0.9% NaCl, Tocris) was injected unilaterally into the dorsal right hippocampus. After kainate injection, all mice showed symptoms of status epilepticus, such as convulsions, chewing and immobility (Riban et al., 2002;Heinrich et al., 2006). All animal procedures were in accordance with the German Animal Welfare Act and were approved by the regional council in Freiburg, Germany.

Electrophysiological recordings
Starting two weeks after kainate injection, when recurrent EA occurred reliably, mice were re-anesthetized to implant platinum-iridium wire electrodes into both dorsal hippocampi (dataset A) or steel electrodes at the ipsilateral dorsal hippocampus and medial entorhinal cortex (dataset B). Signals were amplified 1000ϫ, filtered at 1-5000 Hz (MPA8I and PGA32; Multi Channel Systems), and sampled at 10 kHz (Power 1401 ADC; Spike2 software, RRID: SCR_000903, Cambridge Electronic Design). Data were recorded in the molecular or granule cell layer of the dentate gyrus at 14 -39 d after kainate injection from awake, freely behaving, chronically epileptic animals (three to nine sessions per mouse, each lasting 1.5-3 h, 105 sessions total). Mice from dataset A were recorded on consecutive days, mice from dataset B every other day. Large, putative behavioral seizures (compare Zeidler et al., 2018) were identified visually based on morphology and correlated EA across at least two different brain regions (dataset A: ipsilateral dentate gyrus and contralateral dentate gyrus; dataset B: ipsilateral dentate gyrus and entorhinal cortex).

Data analysis Overview
We delimited bursts of epileptiform spikes as EA events (Fig. 1) and classified them according to spike load (Fig. 2) using a machine learning algorithm. From the succession of classified EA events we derived phases of larger temporal scale (Figs. 3,4;. LFPs obtained from ipsilateral septal recording sites were analyzed in detail. The first 10 min of each recording were discarded to avoid residual influence of the brief isoflurane anesthesia applied to connect the headstage. Periods with movement or other recording artifacts were identified visually and masked before further analysis.

Detecting epileptiform spikes in the LFP
LFPs were downsampled to 500 Hz. We then calculated spectrograms (fast Fourier transform, 256 ms sliding windows, normalized per frequency bin) and averaged the spectrograms in the 4 -40 Hz band. We detected peaks in this activity band as epileptiform spikes. The dead time before a new spike could be detected was 0.0833 s. In addition, we identified slow high-amplitude deflections that escaped the spectrogram-based detection. Such peaks were detected if they crossed a threshold of Ϯ4.5ϫ the SD of the LFP baseline, which consisted of episodes where no peaks had been detected in the spectrogram. Cutouts around detections were then sorted according to their waveform by principal component analysis (five clusters identified by a Gaussian mixture model on the first three components, separately for positive and negative spikes). Events in the cluster with the lowest average peak-to-peak amplitude did not conform to the typical morphology of epileptiform spikes and were discarded. Overlapping peaks within dense bursts were directly included into the pool of spikes without prior sorting.
The spike detection algorithm was validated by a human expert who manually marked spikes in 197 randomly selected 20-s epochs (1550 spikes in total). Automatically detected spikes that occurred within a 150-ms tolerance window around the manually selected spikes were considered as true positives (TPs) and spikes outside this window as false positives (FPs). Manually selected spikes without corresponding automatically detected spikes within the tolerance window were considered false negatives (FNs). The spike detection algorithm had a sensitivity score [TP/(TPϩFN)] of 0.87 and a precision score [TP/ (TPϩFP)] of 0.90. The average FP rate was 1.9 FP/min.

Delimiting bursts
A popular method to delimit bursts is based on interspike interval distributions: When these distributions are bimodal, the valley can be used to separate short intraburst intervals from long interburst intervals (Selinger et al., 2007). The interspike interval distributions of our recording sessions showed valleys or plateaus at around 2.5 s (Fig. 1D). We therefore used this value as a threshold New Research to delimit spike bursts. Spikes closer than 2.5 s were grouped into the same event. Bursts closer than 3.5 s were joined. Our dataset contained 22,578 bursts and 17,814 solitary spikes.

Classification of bursts
To classify bursts and derive a measure of spike load we used a self-organizing map (SOM; Kohonen, 2013). A SOM can transform a high-dimensional feature space into a two-dimensional map of nodes such that neighboring nodes represent similar feature combinations. The map thus creates a, not necessarily linear, continuous representation of the dataset in feature space that can be readily visualized and interpreted.
Feature selection and SOM training To develop a measure of spike load, we exclusively used features associated with the time series of spikes contained in a burst. Since subtle differences in electrode placement and recording quality can affect signal amplitude, we did not use measures related to LFP amplitude. The feature vectors contained values of the following weighted quantities for each burst: [log(number of spikes) ϫ 2; log(median interspike interval) ϫ 2; (SD of interspike interval) ϫ 1]. Feature values were z-scored across all bursts from all animals and sessions (N ϭ 12,373) and then used to train the SOM with the batch algorithm provided by Carrasco Kind and Brunner (2014) Spike load index The node representing the highest spike count and lowest median interspike interval (H‫)ء‬ was located at the top left. For each node, we calculated the Euclidean distance between its feature vector and the vector of H‫.ء‬ We then converted these distances to values between 0 and 1 to define a spike load index, with spike load index ϭ 1 assigned to H‫.ء‬ Categories To identify groups of nodes with similar features we hierarchically clustered the feature vectors corresponding to the nodes (Fig. 2C) and selected an appropriate number of categories according to the Thorndike criterion (Thorndike, 1953;Extended Data Fig. 2-1). The dendrogram was constructed using Ward's method (Ward, 1963) implemented in the Python package scipy.cluster.hierarchy (version: 0.17.0). We thus derived three categories and labeled them according to the average spike load index of their nodes "high-load," "mediumload," and "low-load."

Re-mapping spike load index and categories to
bursts Each burst inherited cluster affiliation and spike load index from its best matching node, i.e., the node of the SOM whose features were most similar to the features of the burst (criterion: smallest Euclidean distance in feature space). High-load and low-load bursts gathered in opposite regions of the SOM ( Fig. 2A,B).
Bursts with less than five spikes Bursts with less than five spikes were not used to train the SOM because preliminary tests had shown that they would disrupt the continuous representation of the map due to their high internal variability with respect to the features used. Across sessions, the rates of burst with less than five spikes strongly correlated with the rates of low-load bursts (N ϭ 105, Kendall's ϭ 0.43, p ϭ 6.6 ϫ 10 Ϫ11 ) and they were therefore assigned to the lowload category (Fig. 2B,E). Burstiness scores can adopt values ranging from -1 (indicating completely regular interevent intervals) over 0 (Poissonian interval distribution) to 1 (extreme burstiness of events). Spike series in the recording sessions (black, N ϭ 105) have distinctly higher burstiness than their Poissonian surrogates (gray, N ϭ 105,000). D, Defining bursts. Interspike intervals from three sessions representing the 5th (deep purple), 50th (green), and 95th (yellow) percentile of spike rates. Spike rates are displayed next to the percentiles. We interpreted the location of the valley/plateau in these distributions as separating short intraburst intervals from long interburst intervals and grouped spikes closer than 2.5 s (red vertical line) into the same burst event. Interspike interval distributions of surrogates generated with matching spike rates are shown as dotted lines in the color of their reference recording. E, We thus converted continuous LFP patterns into time series of discrete EA events: bursts (red) and solitary spikes (black ticks). Extended Data Figure 1-1 compares the succession of spikes and bursts to a rate profile generated from spike timings and shows that bursts match distinct peaks in the rate profile.

Defining phases
Transforming the LFP into a sequence of classified EA events and silent periods allowed us to discern their patterns on a larger scale. In a series of processing steps, we identified and delimited three main temporal phases : (1) clusters of high-load bursts, (2) transition phases around these high-load clusters, and finally (3) background phases. The burst categories (high-load, medium-load, and continued into categories by a dendrogram analysis (see C). B, Examples of EA events after classification. Vertical ticks above the LFP traces mark individual spikes found by our detection algorithm. Horizontal bars delimit the extent of spike bursts, with color indicating the category assigned to the burst. Lower case letters (a-i) refer to their placement on the SOM. Bursts with less than five spikes were not classified on the SOM and later added to the category low-load bursts. The inset shows how many bursts of each category were present in the dataset. Participation of the contralateral DG and ipsilateral EC varied for bursts of different categories (Extended Data Figure 2-2). C, Nodes of the SOM were hierarchically clustered using Ward's dendrogram method (Ward, 1963). To create the dendrogram, nodes were successively merged according to their proximity in feature space. The y-axis shows the distance in feature space between merges. The gray dotted line indicates the threshold used to obtain the three categories as suggested by Thorndike's criterion (Extended Data Figure 2-1). D, Spike rate against duration of bursts. Each dot marks one burst and is colored according to its category. Bursts of the same duration or spike rate can belong to different categories. Purple crosses mark visually identified, putative behavioral seizures. E, Overview of terminology and event rates. The population of EA events had been split into spike bursts and solitary spikes (so.sp.) by burst detection (Fig. 1). Bursts with more than four spikes were classified using a SOM (A) and subsequent hierarchical clustering (C), resulting in three main categories. The category high-load bursts included events with the highest spike load index (idx. ϭ 1; inset). Seizures that had been identified visually beforehand all scored spike load index ϭ 1 in the SOM. Numbers indicate median and 10th-90th percentile range of event counts per hour across sessions (N ϭ 105). Black bars below the time series indicate high-load clusters defined according to D. Our definition allowed for high-load clusters to consist of a single high-load burst. B, Burstiness calculated from the intervals between all high-load bursts in the full dataset (value marked by red triangle) and for each session separately (distribution shown in black, N ϭ 87, including only sessions with at least five intervals between high-load burst). The distribution of surrogate burstiness scores (gray) was derived from surrogate series of high-load bursts matching the rate of high-load bursts in the whole dataset (N ϭ 1000). C, Number of sessions scoring burstiness values significantly below ("regular"), significantly above ("clustered"), or not significantly different from surrogate interval distributions. For each recording, we generated 1000 surrogate interval distributions of matching rate and compared their burstiness distribution to the original data. Significance level: ␣ Յ 0.05. D, The distributions of intervals between high-load bursts were used to define clusters of high-load bursts. Colored lines indicate distributions for individual animals; the average distribution across animals is shown in black. High-load bursts separated by Ͻ3 min (vertical line) were grouped into the same high-load cluster. E, Distribution of cluster size, i.e., the number of high-load bursts within high-load clusters. F, Distribution of duration of high-load clusters. Gray, high-load clusters consisting of a single high-load burst; black, clusters containing more than one high-load burst. G, Fraction of time spent in high-load clusters for all sessions (circles, colors as in D). Horizontal dashes indicate medians across recording sessions of the same mouse (separated by vertical lines). low-load bursts) used to define these phases are described in detail in the results section.
High-load clusters As high-load bursts clustered temporally ( Fig. 3A-C), we defined the temporal extent of these aggregations as high-load clusters. To delimit highload clusters, we analyzed the distribution of intervals between high-load bursts. At 3 min, the peaks of these distributions had leveled off (Fig. 3D). Hence, high-load bursts closer than 3 min were grouped into the same high-load cluster. A high-load cluster started with the onset of its first high-load burst and ended with the offset of its last high-load burst. High-load clusters could consist of individual high-load bursts.
Transition phases with increased burst rate We observed that the probabilities for low-load and medium-load bursts were increased around high-load clusters. To detect putative transition phases we z-scored the intervals between bursts and applied the CUSUM algorithm (Gustafsson, 2001) to identify change points in the interburst interval sequence between two high-load clusters (threshold: cumulative sum of change larger than 1.5; drift compensation: 0.1). We then selected the last change point toward shorter intervals before a high-load cluster to delimit the beginning of the pre-phase. Similarly, the first change point toward longer intervals following a high-load cluster was defined as the end of the post-phase. Transition phases were not assigned when (1) the CUSUM algorithm did not detect a change point or when (2) the fraction of time spent in bursts during a candidate transition phase was less than the fraction of time occupied by such events in the remaining unassigned portion of the recording session. We thereby ensured that bursting in transition phases was above the average of the session excluding high-load clusters.

Transitions with depression of bursting after high-load
clusters Some high-load clusters were followed by extended burst-free periods, reminiscent of post-ictal depression. These periods could be components of the transitions from high-load clusters. To identify depression periods, we pooled interburst intervals from unassigned periods of all recordings. We assigned burst-free periods longer than the 95th percentile (112 s) of the interburst interval distribution to post-phases (Extended Data Fig.  4-1B).

Exclusion of session borders and regions around arti-
facts To avoid potential contributions of ongoing postphases at the beginning of a recording to the analyses of background phases, we excluded the first 3 min of each Post-phases could also consist of a depression of EA that was sometimes followed by a rebound of EA (Extended Data Fig. 4-2). B, Distributions of the start of pre-phases (left) and end of post-phases (right) with respect to the onset (left), respectively, offset (right) of high-load clusters. These distributions are equivalent to the duration of the phases. Note that the y-axis is scaled logarithmically. C, Spike rates for pre-phases (N ϭ 347), post-phases (N ϭ 311), and background phases (N ϭ 340). b Within each box the middle bar indicates the median, while the lower and upper boundaries mark the 25th and 75th percentiles, respectively. Whiskers extend 1.5 times the range from the 25th to 75th percentiles. Crosses show data points outside of the whisker range. D, Burstiness of spikes in all phases combined, i.e., the whole dataset ("whole") and calculated for each type of phase separately. Significance was derived by comparison to 1000 spike trains with interspike intervals generated from a Poissonian distribution and matching rate. All phases scored burstiness values higher than any of their Poissonian surrogates ‫‪p‬ءءء(‬ Ͻ 0.001). E, Rates of medium-load (left) and low-load (right) bursts in pre-phases, post-phases, and background phases. c,d Statistics: Kruskal-Wallis test for differences across groups (p Ͻ 0.001; C, E), followed by pairwise Mann-Whitney test with Bonferroni correction ‫‪p‬ءءء(‬ Ͻ 0.001; C, E). recording (following the 10 min excluded to avoid the effects of anesthesia) unless a high-load cluster ended in this window. To avoid contributions of putative prephases ending after the recording, we likewise excluded the last 5 min of a session from background phases. These settings were derived from the distributions of pre-phase and post-phase durations (Fig. 4B). We further excluded phases rendered ambiguous by movementinduced artifacts.
Background phases To avoid misidentification of short lapses between high-load dynamics as background activity, we excluded periods Ͻ90 s ("short") from background phases (Extended Data Fig. 4-1C). Since mediumload bursts were strongly associated with transition dynamics, we further excluded periods with Ͼ15% of their time occupied by medium-load bursts (Extended Data Fig. 4-1C). All remaining periods not classified as any of the above were considered background phases. Recording time was spent mostly in background phases (43%) and high-load clusters (40%; Extended Data Fig.  4-1E). Only complete high-load clusters were used to calculate the statistics of high-load clusters (data presented in Fig. 3E,F). Likewise, only data of background phases of known duration is shown in Figure 5B-G.

Code accessibility
Data analysis was conducted using a custom algorithm developed in Python (optimized for versions 2.7 and 3.4) 3). Note that burstiness was calculated from interspike intervals of all background phases pertaining to an animal. Hence background phases with more intervals had a greater influence on overall burstiness of an animal. Because longer background phases, i.e., those with more intervals and/or higher spike rates, scored higher burstiness, the center of mass of the burstiness distribution per animal is shifted toward higher burstiness values compared to the center of mass of the distribution per background phase shown in B. E, Correlation between duration of background phase and rate of low-load bursts during background phase. Only background phases of known duration were included (N ϭ 240, ϭ 0.28, p surrogate Ͻ 0.001). F, Relation between background rate of low-load bursts and duration of background phase per animal. Each line shows the least-squares regression across background phases for each animal individually (same data as in E but analyzed per animal). G, Correlation between average duration of background phase and average background rate of low-load bursts per animal (N ϭ 18, ϭ 0.20, p surrogate ϭ 0.5). The slopes were positive in 15 out of 18 animals. The fraction of time spent in background phases per session/animal and the survival function of background phases are shown in Extended Data New Research on Linux and Windows machines. The code for detecting and classifying EA and delimiting phases is available on reasonable request. The Python code for simulating ratemodulated spike trains (Extended Data Fig. 6-3) is freely available at doi:10.5281/zenodo.3376778.

Statistics
Statistical tests and data analyses were performed with Python 2.7 (RRID: SCR_008394, Python Software Foundation). Since most variables were not normally distributed, statistical dependence was assessed using the nonparametric correlation coefficient Kendall's tau (, p value indicated as p ) in all cases. Correspondingly, regression lines show the Theil-Sen estimator (median of slopes through all pairs of sample points), unless stated otherwise. Burstiness was defined as standarddeviation Ϫ mean / standarddeviation ϩ mean of interevent intervals (Goh and Barabási, 2008). Survival functions were calculated according to the Kaplan-Meier method (Lifelines Python package, version: 0.9.4, doi:10.5281/zenodo.805993). Differences between groups were tested with the Kruskal-Wallis test. Differences between pairs of independent groups were tested by a Mann-Whitney U test. For pairwise comparison of dependent samples, we used Wilcoxon signed-rank test; p values were corrected for multiple testing with the Bonferroni method. A one-sample t-test was used to assess whether the distribution of fitted regression slopes significantly deviated from a random population with mean ϭ 0. Statistical power was assessed with G‫ء‬Power 3.1 (Faul et al., 2009). Observed power and exact p values are reported in Table 1. Datasets A and B yielded comparable results throughout this study and were therefore pooled. All value ranges reported in this study are given as the 10th-90th percentile. Per-session background rate of events was calculated by counting events within all background phases of a session and dividing by the total time spent in background phases in a session. We defined the rate of high-load bursts as the count of all high-load bursts in a session divided by the total duration of that session excluding session borders and periods spent with movement artifacts. Averages per animal were computed as weighted averages, with duration of background phase (Fig. 5D,G) or time spent in background phases in a session (Fig. 6C,F) as a weighting factor.

Surrogate tests
Significance of burstiness scores (Figs. 1C, 3B,C, 4D) and of the correlations shown in Figures 5, 6 was assessed using surrogate tests. Surrogate event series were generated by concatenating intervals derived from Poisson processes. The rate parameter rate adjusted of the Poisson process was set to yield surrogate series matching the event rates in the original data (rate observed ) after adjustment for dead time (T dead ) and event duration (T event ) as rate adjusted ϭ rate observed / 1 Ϫ rate observed ·͑T dead ϩ T event ͒ , with T dead ϭ 0.0833s and T event ϭ 0 s for surrogate spike series. For surrogate burst series we used T dead ϭ 3.5 s and defined T event as the mean duration of burst events occurring during all episodes for which the surrogate was generated. rate observed was a constant calculated as a duration-weighted average of event rates from all episodes considered, e.g., when assessing burstiness of spikes in all background phases rate observed was defined as the average rate of spikes across all background phases, weighted by background phase duration. The duration of surrogate bursts was randomly drawn from the pool of burst events in the respective episodes. Surrogate tests based on event series are shown and further described in Extended Data Figures 5-2A-C, 6 -2A-C. Significance of correlations across the set of animals was assessed with a permutation test. To generate surrogate animals, we randomly permuted episodes, i.e., background phases (Extended Data Fig. 5-2B,D) or sessions (Extended Data Fig. 6-2B,D), from the pool of all episodes such that each surrogate animal was randomly assigned as many episodes as it originally contributed to the dataset.
To test significance, the original values were compared to distributions of 1000 surrogate values. The p value p surrogate was defined as the fraction of surrogate values above (right-tailed) or below (left-tailed) the original value.

Results
This study addresses the question what patterns of EA occur in LFPs recorded from the hippocampus of kainateinjected mice and how these patterns interact.

Epileptiform spikes come in bursts
As is typical for this animal model, epileptiform spikes were frequent in our recordings. These spikes occurred mainly as components of larger bursts or aggregated loosely into smaller bursts that contained fewer spikes. In contrast, large, putative behavioral seizures were rare (estimated range: 0.0 -9.6/d).
To describe burst patterns systematically, we first detected spikes using a custom algorithm and transformed LFPs into spike trains (Fig. 1A,B). Spike trains in all sessions had burstiness scores (0.37-0.70) significantly higher than Poissonian surrogate trains (burstiness -0.04 to 0.00; p Ͻ 0.001 for all recordings; Fig. 1C) and nonnormal interspike interval distributions (Fig. 1D). This corroborated our observation that spikes preferentially occurred in bursts. We then delimited bursts by grouping spikes (Fig. 1D), thereby creating time series of two general types of EA events: Bursts and solitary spikes (Fig.  1E). Such a time series approach captures the discrete nature of the process more directly and requires fewer parameters compared to an analysis of continuous rate profiles (Extended Data Fig. 1-1).

Categories of bursts
We observed that bursts did not simply differ in a single property but rather in a composite of multiple features. To classify bursts according to features associated with spike load we used a machine learning algorithm (SOM; Figure 6 High background rates of low-load bursts, reduced susceptibility to high-load bursts. A, Correlation across sessions between the average rate of low-load bursts during background phases and the rate of high-load bursts. Only sessions with background phases were included (N ϭ 91, ϭ -0.43, p surrogate Ͻ 0.001). B, Relation between background rate of low-load bursts and rate of high-load bursts per animal. Each line shows the least-squares regression across all sessions of each animal individually (same data as in A but analyzed per animal). The slopes were negative in 15 out of 17 animals. Line width indicates the number of background phases (thin: N Ͻ 5, thick: N Ն 5). C, Correlation between the average background rate of low-load bursts and the average rate of high-load bursts per animal (N ϭ 18, ϭ -0.33, p surrogate ϭ 0.31). The gray dot shows data for an animal that had only one session containing background phases. D, Correlation across sessions between the average rate of low-load bursts during background phases and the percentage of time spent without EA (N ϭ 91, ϭ 0.39, p surrogate Ͻ 0.001). EA-free episodes were defined as periods within the background phase lacking EA events. Long EA-free episodes occurred throughout the extent of interictal phases (Extended Data Fig. 6-1). j E, Relation between background rate of low-load bursts and percentage of time spent without EA per animal. Each line shows the least-squares regression across all sessions of an animal (same data as in D but analyzed per animal). The slopes were positive in 14 out of 17 animals. F, Correlation between the average background rate of low-load bursts and the total percentage of time spent without EA per animal (N ϭ 18, ϭ 0.29, p surrogate ϭ 0.32). Surrogate tests based on Poissonian event trains for the relationships across sessions (A, D) and permutation tests for the relationships across the set of animals (C, F) are shown in Extended Data Figure 6-2. In Extended Data Figure 6-3, we compare our results to simulations of rate-modulated spiking processes. Fig. 2A,B) and assigned a spike load index ranging from 0 (loose aggregation of spikes) to 1 (long spike burst with high rate) to each burst. A hierarchical clustering analysis ( Fig. 2C; Extended Data Fig. 2-1) allowed us to distinguish three distinct burst categories: low-load, medium-load and high-load bursts. Spike load was loosely linked to spike rate and duration in that high-load bursts showed elevated values for both. Short low-load bursts could also have high spike rates (Fig. 2D), whereas long low-load bursts tended to have low spike rates. Medium-load bursts covered intermediate ranges of duration and spike rate in highly variable combinations. Bursts corresponding to visually identified behavioral seizures had some of the highest spike rates and longest durations.
Although this classification was based specifically on the spike statistics of the bursts, bursts within a category showed consistent morphologic features ( Fig. 2B; Extended Data Fig. 2-2). Bursts corresponding to visually identified large behavioral seizures had spike load index ϭ 1 (N ϭ 18). During such bursts, we observed highly synchronous activity on the contralateral side and in the entorhinal cortex (Extended Data Fig. 2-2A). The remaining bursts with spike load index ϭ 1 (N ϭ 125) and the other high-load bursts with spike load index Ͻ1 were accompanied by intense EA in the contralateral dentate gyrus but only sparse spiking in the entorhinal cortex (Extended Data Fig. 2-2B,C). In low-load bursts, the spike component was typically followed by a clearly discernible wave of opposite polarity, smaller amplitude, and longer duration, similar to a spike-wave discharge (comparable to type 1 spikes in Chauvière et al., 2012). During low-load bursts, contralateral and entorhinal activity was sparse or absent (Extended Data Fig. 2-1D). Medium-load bursts often contained periods in which spikes were more densely packed than in low-load bursts and where the wave component was less salient or absent. During medium-load bursts, contralateral activity was typically more pronounced than during low-load bursts.

High-load bursts form clusters
The temporal succession of the three burst types suggested patterns of higher order: Periods dominated by high-load bursts appeared to alternate with extended phases containing only low-load bursts (Fig. 3A). In line with this observation, we found that high-load bursts clustered significantly across all sessions pooled (p Ͻ 0.001 burstiness; Fig. 3B) as well as in most sessions individually (64 sessions significantly clustered with p Ͻ 0.05, 20 not clustered; Fig. 3C). We delimited clusters of high-load bursts as the temporal extent of high-load bursts connected by intervals shorter than 3 min (Fig. 3D). Note that within such clusters intervals between high-load bursts were typically not devoid of EA but could contain additional bursts of lower load and solitary spikes (Fig.  3A). For consistency, we treated isolated high-load bursts (36%; Fig. 3E) as high-load clusters. High-load clusters with more than one high-load burst could last from 20 s up to ϳ 50 min (Fig. 3F). The fraction of time spent in high-load clusters fluctuated considerably from session to session (Fig. 3G) and was independent of the time elapsed since injection of kainate ( ϭ 0.01, p ϭ 0.8) and of the daytime of the recording (ANOVA a contrasting the morning and afternoon group).

Transitions are rich in low-load and medium-load bursts
High-load clusters were often preceded and followed by aggregations of low-load and medium-load bursts (Fig.  4A). We identified these pre-phases and post-phases through change points in the interburst interval series (Extended Data Fig. 4-1A). Depressions following highload clusters were also considered as post-phases and could be concluded by rebound bursts (Extended Data Figs. 4-1B, 4 -2). Of all high-load clusters, 86% were preceded by pre-phases and 81% followed by postphases.
After delimiting high-load clusters and transition phases, we defined background phases by exclusion (Extended Data Fig. 4-1A). Median spike rate was more than twice as high in the pre-phase and post-phase than in the background phase (Fig. 4C). b Spikes in each of the phase types had a strong tendency to occur in bursts (Fig. 4D), corroborating that even during background phases spiking was not random but structured in bursts.
Median burst rate was significantly higher in the prephase than in the post-phase (Fig. 4E). c,d During background phases, low-load bursts occurred at substantially higher rates than medium-load bursts (median rate and range for low-load bursts: 0.84/min, 0.00 -1.90/min; for medium-load bursts: 0.06/min, 0.00 -0.53/min). e Transition phases thus appeared to be more specifically associated with medium-load bursts than with low-load bursts.

Background phase duration is positively correlated with burstiness and the rate of low-load bursts
Duration and EA content of background phases varied (Fig. 5A, Extended Data Fig. 5-1). We therefore asked, whether there was a connection between the duration of background phases and the rate and structure of EA within them. Across the pool of background phases, the burstiness of spikes within a background phase correlated positively with its duration ( ϭ 0.27, p surrogate ϭ 0.009; Fig. 5B; Extended Data Fig. 5-2A), i.e., burstier background phases lasted longer. To test whether this relationship held for each animal individually, we calculated least-squares regressions between burstiness and duration across all background phases per animal (Fig.  5C). The distribution of regression slopes significantly differed from random (t-test f , p ϭ 7.8 ϫ 10 Ϫ3 ) and most slopes were positive (15 out of 18 mice), demonstrating a consistently positive relationship between burstiness and duration of background phases. However, animals with on average longer background phases did not score significantly higher overall burstiness ( ϭ 0.07, p surrogate ϭ 0.3; Fig. 5D; Extended Data Fig. 5-2B). Hence, while overall burstiness of background spiking in a mouse did not indicate the average duration of its background phases, burstiness and background phase duration were positively correlated across background phases of most mice individually.
Background phases also lasted longer if they had higher rates of low-load bursts ( ϭ 0.28, p surrogate Ͻ 0.001; Fig. 5E; Extended Data Fig. 5-2C). This relationship was reproduced for most animals individually (15 out of 18; Fig. 5F). g However, the average background rate of low-load bursts of a mouse did not indicate the average duration of its background phases ( ϭ 0.20, p surrogate ϭ 0.5; Fig. 5G; Extended Data Fig. 5-2D).

Sessions with higher rates of low-load bursts in background phases have lower rates of high-load bursts
Longer background phases in a recording session would leave less time available for high-load bursts. Alternatively, the fact that long background phases were associated with increased background rates of low-load bursts could be due to higher excitability, which in turn might lead to higher rates of high-load bursts. The former situation would lead to an anti-correlation between the background rates of low-load bursts and the rates of high-load bursts across sessions, whereas the latter would result in a positive correlation. In addition, the clustering of high-load bursts makes it difficult to predict the specific relation between these EA patterns. We found that across sessions, the average background rate of low-load bursts was significantly anti-correlated to the rate of high-load bursts ( ϭ -0.43, p surrogate Ͻ 0.001; Fig.  6A; Extended Data Fig. 6-2A). Regressions across sessions had negative slopes in 15 out of 17 mice (t-test h against random distribution of slopes: p ϭ 2.8 ϫ 10 Ϫ5 ; Fig. 6B). This indicates that the inverse relationship between the background rate of low-load bursts and the rate of high-load bursts is robust across mice. Again, there was no significant (anti-)correlation between the average rate of high-load bursts and the average background rate of low-load bursts of a mouse ( ϭ -0.33, p surrogate ϭ 0.05; Fig. 6C; Extended Data Fig. 6-2B). The background rate of low-load bursts thus more closely reflected the state of a mouse during a session than its overall condition.
Although elevated background rates of low-load bursts were linked to fewer high-load bursts in a session, lowload bursts themselves are also pathologic activity. We therefore tested whether the increased rates of low-load bursts in longer background phases would increase the overall time spent with EA. Across sessions, the percentage of time spent without EA in background phases was larger in sessions with higher background rates of lowload bursts ( ϭ 0.39, p surrogate Ͻ 0.001; Fig. 6D; per animal: Fig. 6E i ; across animals: Fig. 6F; survival probability of EA-free episodes: Extended Data Fig. 6-1A). In addition, the probability for long EA-free episodes did not decrease with increasing background phase duration (Extended Data Fig. 6-1B). We thus rejected the hypothesis that elevated background rates of low-load bursts curtail the overall time spent without EA.
In summary, the background rate of low-load bursts was significantly correlated to (1) the duration of the background phase and (2) to the percentage of time spent without EA, but was (3) significantly anti-correlated to rate of high-load bursts. We found that the background rate of low-load bursts more strongly correlated with these signatures of susceptibility than the background rate of spikes or solitary spikes and that these results were robust against changes in event definition (Extended Data Fig. 5-3).

Discussion
We analyzed the temporal structure of epileptiform spiking in LFPs recorded in a mouse model of MTLE. Since large parts of EA consisted of well-defined spike bursts, we chose an approach based on time series. Time series are well suited for correlation analyses to investigate statistical interactions between event types.
In previous studies, large EA events have been described by features such as duration, spike rate, spike amplitude, spike wave form, and the evolution of these parameters over time (Twele et al., 2016). To avoid effects resulting from variations in recording conditions across animals and time, such as electrode properties, glial scarring, etc., we classified bursts based on a set of features best summarized as spike load. Importantly, spike load was only loosely correlated to burst duration and spike rate within a burst; it could differ considerably for bursts of the same duration or spike rate.
Although they were classified according to statistical measures of their spike structure, events in a burst category were morphologically similar. Low-load bursts typically consisted of spikes followed by longer deflections of opposite polarity with lower amplitudes. These bursts were similar to trains or groups of spike-wave discharges reported previously for this animal model (Riban et al., 2002;Twele et al., 2016). Medium-load bursts typically comprised spike trains with higher rates, where wave components were less apparent or obfuscated by subsequent spikes.
In high-load bursts, most spikes were densely packed. High-load bursts often started with loosely spaced spikewave discharges and evolved into dense spiking with progressively lower amplitudes, like type A bursts described by Chauvière et al. (2012) for kainate injected rats. In terms of both morphology and incidence, high-load bursts were similar to what other studies using our animal model labeled as hyperparoxysmal discharges, high voltage sharp waves, or electrographic seizures based on visual inspection and manual classification (Riban et al., 2002;Twele et al., 2016;Zeidler et al., 2018). With respect to long-term dynamics, the clustering of high-load bursts appeared similar to seizure clusters commonly described in rodent models (Lim et al., 2018) and in patients (Karoly et al., 2017). A small fraction of high-load bursts were similar to large behavioral seizures described elsewhere based on their characteristic, stereotypical wave form and highly synchronous activity across recording sites (Riban et al., 2002;Zeidler et al., 2018). Their incidence was comparable to other reports on various rodent models of epilepsy (Rattka et al., 2013;Klee et al., 2017). We did not investigate these in further detail because they were too rare for a meaningful correlation analysis.
Differences between burst categories also manifested in their temporal occurrence. Clusters of high-load bursts were often preceded by aggregations of low-load and medium-load bursts and alternated with phases populated almost exclusively by low-load bursts. This temporal cohesion of morphologically similar events substantiates our classifications and suggests distinct temporal dynamics dominated by different types of bursts. Segregating such temporal phases of distinct EA dynamics is crucial to evaluate relationships between different types of EA. Separating transition phases rich in low-load and mediumload bursts enabled us to specifically address EA in background phases and to investigate how background activity relates to the rate of high-load bursts.
Our findings indicate that the interpretation of low-load bursts should depend on temporal and functional context: While the rates of low-load and medium-load bursts increased before the onset of high-load dynamics, background rates of low-load bursts were anti-correlated to rates of high-load bursts. We further found that higher rates of low-load bursts during background phases were associated with an increased duration of the background phase and a larger relative amount of time spent without any EA. This antagonistic relationship between low-load bursts and susceptibility to high-load bursts is in agreement with the hypothesis that epileptiform spikes could reduce seizure susceptibility (Barbarosie and Avoli, 1997;De Curtis and Avanzini, 2001;Khosravani et al., 2003;Muldoon et al., 2015;Goncharova et al., 2016). On a mechanistic level, epileptiform spikes, and especially the wave component following the sharp spike during spikewave discharges, have been proposed to reflect GABAergic inhibitory processes (De Curtis and Avanzini, 2001;Muldoon et al., 2015). In our recordings, low-load bursts mostly consisted of spike-wave discharges and thus might be interpreted as a signature of boosted inhibition. Conversely, in medium-load bursts, the wave was often less prominent and spikes were more densely packed. While this could be due to an overlap on the signal level, the increase of medium-load bursts preceding high-load clusters could also indicate a gradual breakdown of inhibition. Along the same line, high rates of low-load and medium-load bursts during transition phases could promote high-load dynamics through excessive GABAergic input: high rates of spikes might increase the intracellular chloride concentration and hence make GABA act depolarizing (Pallud et al., 2014;Magloire et al., 2019), thereby enhancing excitability and facilitating high-load bursts. Alternatively, an increase of intracellular chloride concentration could trigger potassium extrusion into the extracellular space and thus enhance excitability (De Curtis and Avoli, 2016;Librizzi et al., 2017).
One could argue that the existence of different categories of bursts with lower or higher spike load and the alternation we observe between background, transition phases and high-load clusters is governed by modulations of the spike rate. Such modulations and an accompanying buildup in susceptibility have been suggested to be due to slowly changing variables such as extracellular ion concentrations or metabolic factors. For example, a slow "permittivity variable" has been used in computational models of seizure generation (Jirsa et al., 2014).
While it is possible that underlying slow modulations contributed to transitions from background to high-load dynamics observed in our data, such modulations alone cannot explain the antagonistic relationship between the rate of low-load bursts in background phases and the susceptibility to high-load bursts. When simulated as rate-modulated Poisson processes (Extended Data Fig.  6-3), nested modulations indeed reproduced high-load bursts and clusters of high-load bursts (Extended Data Fig. 6-3A-F) but they did not result in an anti-correlation between the rate of high-load bursts and the background rate of low-load bursts found in our data (compare Fig.  6A, Extended Data Fig. 6-3H). Neither did this spike rate modulation replicate the observed positive correlation between the rate of low-load bursts in the background phase and the duration of background phases (compare Fig. 5E, Extended Data Fig. 6-3I). To reproduce the correlation statistics we observed would require an additional interaction between the rate of low-load bursts in background phases and high-load dynamics.
The ideas of a slow process modulating excitability and a context-dependent interaction between high-load dynamics and background bursting, however, may be reconciled. Based on experiments and computational modeling, Chang et al. (2018) recently proposed that epileptiform spikes interact with slowly increasing network excitability in a context-dependent manner: When the excitability of the network was low, epileptiform spikes further reduced excitability and thus lengthened interictal states. In a transition range of excitability, sufficiently strong, synchronized firing of a population of neurons could prematurely trigger seizures. In analogy to Chang et al. (2018), high rates of low-load bursts could stabilize background phases in our animals, counteracting a continuous increase of network excitability. Low rates of low-load bursts, however, would allow a net increase of excitability. In the transition range of excitability, increased rates of low-load events could then eventually facilitate high-load dynamics. Since bursts of spikes were a more powerful correlate to highload susceptibility, we hypothesize that spikes grouped into bursts might be particularly effective in keeping excitability in check and thus prolong background dynamics.
In summary, we characterized the nested structure of EA patterns in the kainate mouse model of MTLE using time series analyses and machine learning techniques. Systematically segregating periods of intense high-load bursting and transition phases from background phases allowed us to address the relationship between different types of EA specific to temporal context. Our study corroborates the hypothesis that the role of low-level EA depends on the current state of network dynamics and suggests that the rate of bursts with low spike load could be a powerful correlate of susceptibility dynamics.