Dual mechanisms of ictal high frequency oscillations in human rhythmic onset seizures

High frequency oscillations (HFOs) are bursts of neural activity in the range of 80 Hz or higher, recorded from intracranial electrodes during epileptiform discharges. HFOs are a proposed biomarker of epileptic brain tissue and may also be useful for seizure forecasting. Despite such clinical utility of HFOs, the spatial context and neuronal activity underlying these local field potential (LFP) events remains unclear. We sought to further understand the neuronal correlates of ictal high frequency LFPs using multielectrode array recordings in the human neocortex and mesial temporal lobe during rhythmic onset seizures. These multiscale recordings capture single cell, multiunit, and LFP activity from the human brain. We compare features of multiunit firing and high frequency LFP from microelectrodes and macroelectrodes during ictal discharges in both the seizure core and penumbra (spatial seizure domains defined by multiunit activity patterns). We report differences in spectral features, unit-local field potential coupling, and information theoretic characteristics of high frequency LFP before and after local seizure invasion. Furthermore, we tie these time-domain differences to spatial domains of seizures, showing that penumbral discharges are more broadly distributed and less useful for seizure localization. These results describe the neuronal and synaptic correlates of two types of pathological HFOs in humans and have important implications for clinical interpretation of rhythmic onset seizures.

frequency ranges of HFOs, ripples (~ 80-200 Hz) and fast ripples (~ 200-800 Hz) 13 , have been studied extensively in the rodent and human medial temporal lobe [14][15][16][17] . In hippocampal slices with impaired inhibition, a single cell can induce population bursting with pathological fast ripples 17 , whereas when inhibition is maintained, the same cells induce oscillations that restrain population activity in time 18,19 . These studies have shown that pathological bursting is characterized by reduced spike-timing reliability across a population during a burst, which results in a broadening of the high frequency LFP spectrum and increased prevalence of fast ripples 15 . To date, there have been few studies aimed at identifying the cellular correlates of pathological high frequency LFPs in the human neocortex, and the context of these correlates within the framework of seizure dynamics.
Previously, our group focused on ictal high frequency LFP as a method of identifying cortex that has been recruited into a seizure, correlating these high-amplitude, phase-locked, broadband ripples with multiunit firing from nearby microelectrode arrays that exhibited signatures of recruitment into the seizure. That is, sites demonstrating a well-defined interictal to ictal transition (tonic firing of the ictal wavefront, followed by burst firing synchronized to each epileptic discharge; Fig. 1) as opposed to sites demonstrating heterogeneous, nonphase-locked firing. Since pathological phase-locked firing is more likely to produce increased high-frequency LFP in macroelectrode recordings than heterogeneous non-phase-locked firing 20,21 , we focused on ictal highfrequency LFP beginning several seconds after seizure onset as a clinically accessible correlate of the specific neuronal firing pattern defining ictal recruitment 4 and as a predictor of postoperative seizure control 22 . In these studies, we also showed that non-phase locked high-frequency LFP at seizure onset correlates poorly with evidence of seizure recruitment, and is a relatively poor predictor of seizure outcome when compared to current clinical standard assessments.
This reasoning appears to be at odds with the long-held principle that electrophysiological changes at seizure onset are more likely to be indicators of epileptic tissue. To resolve this contradiction, we propose that pathological high-frequency LFP may arise from two contrasting mechanisms: paroxysmal depolarizing shifts (PDS) resulting in spatiotemporally jittered population bursting that results in dramatically increased broadband high frequency (BHF) LFP (broadband shift in LFP power, often appears at approximately 70-200 Hz), and oscillations arising from a combination of increased excitation and strong inhibitory interneuron firing, resulting in a pathologically exaggerated narrowband gamma rhythm (true oscillation; approximately 40-80 Hz). Such inhibitory-interneuron sculpted rhythms would be expected to occur in the ictal penumbra, i.e. at sites receiving strong synaptic currents generated from the seizure, yet exhibit the strong inhibitory restraint from surround inhibition preventing that area from being recruited into the ictal core. This pathologically high-amplitude narrowband gamma oscillation may occur throughout a seizure, yet is most apparent at seizure initiation, when the ictal core is small relative to the area of surround inhibition that comprises the penumbra 23 . To test this hypothesis, we examined seizures characterized by phase-locked high-frequency LFPs that occurred before and after the ictal wavefront, and contrast features of these signals in the context of multiunit activity recorded on neocortical and mesial temporal microelectrodes. We then translate these to clinically applicable measures, and leverage the resulting classification of HFOs to demonstrate differences in the ability of each type of HFO to localize seizures.

Results
Temporal profile of seizure recruitment. Simultaneous microelectrode and standard clinical macroelectrode recordings were acquired in patients undergoing intracranial EEG monitoring as part of neurosurgical treatment for pharmacoresistant focal epilepsy. There were 9 subjects with clinical grids and strip electrocorticography (ECoG) and "Utah" microelectrode arrays (MEAs) 24 , and 16 subjects with stereo-electroencephalography (sEEG) and Behnke-Fried depth microwire bundles 25 . Recording was carried out continuously throughout the duration of the patients' hospital stays in order to capture typical seizures. From this large cohort, we included a subset of subjects based on the following criteria: (1) rhythmic slow oscillations at seizure onset, in a region that included at least one microelectrode recording site, (2) high signal-to-noise multiunit action potentials (see Methods) recorded in the time periods leading into these seizures, and (3) multiunit activity (MUA) demonstrating the classic pattern of ictal recruitment, i.e. a period of tonic firing followed by repetitive burst firing 23,26 . Seizures meeting these criteria allowed us to study epileptiform discharges locally before and after the passage of the ictal wavefront, that is before and after the failure of strong feedforward inhibition in the local tissue. Six spontaneous seizures in four (three female) patients met these criteria, four seizures in the MEA group (two patients) and two seizures in the Behnke-Fried group (two patients), with stereotactically placed depth electrodes and microwires in mesial temporal lobe (MTL) structures. Clinical characteristics of these patients are provided in Supplementary  Fig. 2a-c. The timing of ictal invasion of each microelectrode site was determined from the stereotypical barrage of asynchronous multiunit activity that signals this transition (Fig. 2d) 23 . The ictal wavefront functionally and temporally divides the seizure into two stages: a pre-recruitment stage, in which the seizure has not yet invaded the microelectrode recording site, and a post-recruitment stage following seizure invasion 5 . The pre-recruitment stage is the temporal analog of the penumbra, as described in our earlier reports ( Fig. 1) 23 . Consistent with our prior studies, the transition between the two stages could not be directly visualized in low frequency (< 50 Hz) recordings, which is the typical range monitored with standard clinical EEG. Ictal discharges were detected by identifying peaks in the high frequency LFP (Fig. 2e,f) 8 and MUA (Fig. 2d). Each detected discharge was matched between the two recording modalities and classified as either pre-recruitment or post-recruitment depending on its timing relative to the ictal wavefront recorded on the nearest microelectrode. In each seizure, 287.5 ± 51.4 (mean ± s.d.) discharges were detected. The same number of pre-and post-recruitment discharges were retained for further analysis (N = 154.8 ± 116.2 (mean ± s.d. per seizure)) and differences between neurophysiological features of these discharges were tested for in a statistically balanced fashion.
Comparing the normalized amplitude of MUA firing rate and broadband high frequency (BHF) amplitude during each discharge yielded two discriminable axes for differentiating pre-vs. post-recruitment discharges. 74.4% of pre-recruitment discharges had higher normalized MUA firing rates than normalized BHF amplitudes, and 94.5% of post-recruitment discharges had higher normalized BHF amplitudes than normalized MUA firing rates. Accordingly, and as previously observed leading up to seizure termination 5 , normalized MUA firing rates and BHF amplitudes across all discharges were significantly anticorrelated (  Fig. 1), despite overall shifts in the relative amplitudes of MUA firing rate and BHF power. Although MUA is not available from clinical ECoG or sEEG recordings, BHF amplitude differences still provided an informative axis of discriminability between preand post-recruitment discharges. For macroelectrodes, the peak LFP amplitude and BHF amplitudes during discharges were significantly correlated ( Fig. 2h; Pearson's r = 0.42, N = 400, p < 10 -6 ).
Differential MUA coupling to narrow-band HFOs before and after recruitment. We hypothesized that maintenance of inhibition is critical for the generation of pre-recruitment HFOs. That is, an exaggerated expression of the feed-forward inhibition mechanism is associated with narrowband gamma oscillations, which have been associated with parvalbumin-positive interneurons [27][28][29] . We began by examining proposed multiunit indicators of maintained inhibition 30 . Specifically, we hypothesized that if inhibition were maintained during an epileptiform discharge, coupling of multiunit firing to the dominant LFP frequency (Supplementary Fig. 2; see methods for operational definition) of the surrounding population should exhibit phase symmetry and reliability, which would reflect a true narrowband oscillation of excitatory and inhibitory drive. Alternatively, if inhibition has failed, multiunit events should couple to more asymmetrically distributed phases of the dominant high frequency activity.
Our hypothesis was confirmed by the multiunit event coupling relationship between pre-and post-recruitment discharges, examples of which are shown in Fig. 3a,b. Spike coupling distributions changed significantly before and after recruitment into a seizure (permutation Kuiper tests, all V 2000 > 7.34, all p < 0.001; See methods, supplementary Fig. 2). More specifically, neocortical unit firing during pre-recruitment discharges coupled to two specific and opposing phases of narrowband gamma oscillations, in a bimodal circular distribution (Fig. 3c). However, neocortical unit firing during post-recruitment discharges coupled more broadly to a single phase of gamma oscillations in a unimodal circular distribution (Fig. 3d). The variance of this unimodal distribution was significantly greater than that of the angle-doubled bimodal distribution (two-sample tests for equal concentration parameters on the interval [0 π), all p < 0.007, all F > 1.03). Critically, these differences were only determined by examining the patient-specific, dominant high-frequency bands. If predefined narrowband gamma  and ripple (70-200 Hz) frequency ranges were used, no significant differences in unit coupling between pre-and     Figure 3e shows multiunit event coupling from a representative MTL unit that was recruited into the seizure, demonstrating a shift in phase-locking between pre-and post-recruitment events (permutation Kuiper test, V 200 = 6.8, p < 0.001). However MTL units that were not recruited into the seizure (n = 6) maintained the bimodal spike-timing distribution during post-recruitment discharges ( Fig. 3f; permutation Kuiper test, V 200 < 2.16 , p > 0.61). Not only do these data provide evidence for differential unit coupling in different spatial domains in the same seizure, these results show that both neocortical and MTL multiunit events couple to the dominant narrowband gamma frequency of pre-recruitment discharges, with a shift in multiunit event coupling occurring in post-recruitment discharges. This provides support for the hypothesis that pre-recruitment discharges occur in tissue in which feedforward inhibition is maintained, and that this mechanism is similar in both hippocampus and neocortex.
Spatiotemporal profile of seizure recruitment across clinical recordings. As schematized in Fig. 1a, and described in the previous sections, neuronal data from microelectrodes indicate a spatiotemporal structure of seizure expansion, where adjacent groups of cells are successively recruited into the ictal core. We use the terms pre-and post-recruitment to define temporal seizure domains before and after passage of the ictal wavefront. We use the terms "ictal penumbra" and "ictal core" to describe the corresponding spatial domains (Fig. 1a). Even as this successive recruitment occurs across groups of neurons, the effects of the seizure core extend beyond the ictal wavefront, to the penumbra, where inhibition is likely preserved, yet neural ensembles are assailed by discharges from the adjacent, already recruited tissue. In the previous section we showed data indicating disparate coupling between neuronal firing and high frequency LFP during discharges in these two territories. In the following two sections, we examined correlates of these two types of discharges on a larger, more clinically relevant scale, relating changes in the high frequency LFP on clinical macroelectrodes to intense, disorganized, phase-locked action potential firing in the seizure core.
In order to extend our analysis to HFOs recorded on macroelectrodes throughout the brain, we utilized phase-locked BHF, a previously defined ECoG-based measure of recruitment 4,22 , to operationally define the spatial extent of ictal recruitment across the entire coverage area. We maintained the same temporal delineation of recruitment defined by passage of the ictal wavefront on simultaneously recorded microelectrodes. Examples of recruited and unrecruited macroelectrodes are shown in Fig. 4a-d (MTL seizure shown in Supplementary  Fig. 4). Phase-locked BHF yielded clearly bimodal distributions of channels during the post-recruitment phase ( Fig. 4e; t 363 = 2.8, p = 0.006), confirming that not all implanted macroelectrodes were recruited.
We therefore classified discharges recorded on macroelectrodes into 4 mutually exclusive categories that correspond to different spatiotemporal seizure domains: pre-recruitment discharges on core macroelectrodes, post-recruitment discharges on core macroelectrodes, pre-recruitment discharges on penumbral macroelectrodes, and post-recruitment discharges on penumbral macroelectrodes (i.e. discharges after seizure recruitment on the microelectrode array in contacts that are never recruited into the ictal core). That is, in the temporal domain, pre-and post-recruitment discharges were distinguished based on the MUA-based temporal designation described above. In the spatial domain, core and penumbra contacts were distinguished by clustering the bimodally-distributed post-recruitment phase-locked BHF values into the two distinct unimodal components of the overall distributions using Gaussian Mixture Models (see methods). The mean spectra for these four distributions are shown in Fig. 4f. The pink spectrum clearly shows how the post-recruitment core is distinguished from the other spatiotemporal categories by its increased BHF.
As the tissue under the penumbra electrodes is not recruited, we compared penumbral discharges to discharges located within the ictal core, occurring both before and after the ictal wavefront, in order to understand large-scale properties of the penumbra. We found that during the pre-recruitment period, phase-locked narrowband gamma did not differ between core and penumbral macroelectrodes ( Fig. 4g; t 363 = 0.2, p = 0.82). Examining the Bayes Factor from this model suggested that narrowband gamma oscillations occurring in the pre-recruitment period do not distinguish the core and penumbra (Bayes Factor = 1.2*10 -32 ). Unsurprisingly, this similarity was also evident between the phase-locked dominant frequency LFP in the core and penumbra during the pre-recruitment period (t 363 = 1.9, p = 0.05). However, phase-locked dominant frequency power during post-recruitment discharges differed between core and penumbral electrodes ( Fig. 4h; t 363 = 8.5, p < 0.01). While dominant frequencies of pre-recruitment discharges were lower than those from the post-recruitment epoch (prerecruitment mean± std: 57.94 ± 11.38; post-recruitment mean± std: 65.01 ± 16.95; t 728 = − 3.3, p = 0.001), they both represented a strong narrow band gamma component that increased in frequency as the seizure approached recruitment and remained consistent in a high frequency range following recruitment.
Together these results show that in rhythmic onset seizures, the pre-recruitment time period and penumbral spatial domain are marked by a narrowband gamma oscillation whose dominant frequency increases until the point of seizure recruitment. Importantly, this narrowband gamma oscillation does not differentiate between core and penumbra territories, as it can appear across widely distributed brain areas including those distant from the seizure onset zone. Given that all but one of the patients were seizure free after resection surgery (apart from the patient who received responsive neurostimulation, with more than 33 months follow-up at time of writing), this provides additional evidence that pre-recruitment HFOs are not reliable biomarkers of epileptic brain in rhythmic onset seizures. Following recruitment, BHF amplitude increases specifically in the slowly expanding seizure core.
Gauging seizure recruitment with discharge complexity. We next sought to further clarify the features underlying differences between pre-and post-recruitment discharges, and relate these differences to the pathological nature of HFOs previously reported in animal models 15  www.nature.com/scientificreports/ between HFOs during pre-and post-recruitment discharges we observed was that HFOs during pre-recruitment discharges appeared smoother and less complicated than those during post-recruitment discharges. We sought to quantify this feature by measuring entropy of the relative distributions of recorded voltages through the duration of each discharge (differential entropy) and of the relative distribution of LFP frequencies (spectral entropy) during each discharge. These two quantities may alternatively be thought of as respective measures of timedomain and spectral-domain complexity of each epileptiform discharge. As in the aforementioned studies, we compare these measures to the proportion of the spectrum in the ripple and fast ripple ranges (ripple and fast ripple indices; see methods). Figure 5a,b show the distributions of the fast ripple indices and spectral entropies for all discharges and MEA electrodes for a single seizure. After accounting for any random variance that may have occurred across patients and seizures using linear mixed-effects models, we found that pre-and post-recruitment discharges recorded on microelectrodes had no significant difference in fast ripple indices ( Fig. 5c; generalized linear mixed effects model: index ~ 1 + epoch + (1 + epoch | seizure); t 87,723 = − 0.59, p = 0.55; Bayes factor = 0) or differential entropy ( Fig. 5d; generalized linear mixed effects model: entropy ~ 1 + epoch + (1 + epoch | seizure), t 87,723 = 0.005, p = 0.99; Bayes factor = 3.5*10 -120 ). However, spectral entropies were significantly higher for post-recruitment discharges ( Fig. 5e; generalized linear mixed effects model: entropy ~ 1 + epoch + (1 + epoch | seizure), t 87,723 = − 3.09, p < 0.002).
In order to determine whether these measures translated to clinical ECoG recordings, we also calculated similar entropy measures for discharges recorded on macroelectrodes, again employing the four spatiotemporal seizure domains described above. The sampling rates typically used to record clinical data, and that were used to record these data, did not allow for the detection of fast ripples, but a wealth of evidence, both previously presented, and presented herein, suggest that BHF is informative about the intense firing of recruited units 4 . We therefore calculated a "ripple index" (see Methods for operational definition) analogous to fast ripple indices derived from microelecrode recordings to compare with spectral and differential entropy measures derived from the signals recorded on macroelectrodes. This index differs from standard HFO definitions by its basis on the area under the normalized spectrum, rather than amplitude of particular spectral bands. The "ripple index" therefore incorporates broadening of the spectrum, as well as the normalized amplitude of the ripple band, relative to the rest of the spectrum. Figure 5f,g show the distributions of the "ripple index" and spectral entropy for each discharge in a representative seizure. This measure was significantly higher for post-recruitment discharges in the putative ictal core ( These results further confirm, that increased complexity and significant broadening of the LFP spectrum in the ictal core is related to large phase-locked increases in BHF, and disorganized, intense unit firing that localizes specifically to the ictal core following recruitment.

Discussion
Here we provide evidence that early discharges in spontaneous human seizures with rhythmic onset morphology occur in brain tissue in which feed-forward inhibition is maintained. More specifically, we show that discharges occurring before the neuronal signature of recruitment exhibit pathologically high-amplitude phase-locked narrowband gamma-range activity (approximately 40-80 Hz) across a broad cortical area. Multiunit firing couples to these narrowband gamma oscillations in a bimodal pattern during pre-recruitment discharges. Following recruitment, a more restricted cortical area that corresponds to the ictal core exhibits high-amplitude BHF activity that is phase locked to the seizure's dominant rhythm. This BHF activity is tightly temporally linked to phase-locked multiunit activity, reflective of actively seizing brain tissue. These post-recruitment discharges exhibit increased complexity on both microelectrodes and clinical electrodes. We therefore translate our understanding of the dynamics of multiunit activity during seizures to clinical recordings with important implications for seizure localization in the clinical setting, and demonstrate that early discharges in rhythmic onset seizures exhibit marked narrowband gamma, and that these are present outside of seizing brain tissue.
Together these results provide further means to decode the neuronal and synaptic processes underlying the discharges that make up focal seizures, and therefore have important clinical implications. While pre-recruitment discharges indicate seizure activity somewhere, they are likely representative of a response to topologically adjacent seizing tissue. That is, inhibitory firing sculpts excitatory bursts into a fast oscillatory structure, resulting in a high-amplitude narrowband gamma oscillation. In the seizure core, which has been shown to be dominated by paroxysmal depolarizing shifts 26,32 , this sculpting effect is absent, and the appearance of a broadband pseudooscillation, termed BHF here, likely arises from summated, jittered postsynaptic activity 21 . Importantly, early phase-locked narrowband gamma oscillations do not necessarily precede post-recruitment phase-locked BHF. These features of epileptiform discharges could however be extremely informative if interpreted as follows: early discharges with increased phase-locked narrowband oscillatory gamma power are likely not yet recruited into the seizure core, yet they are receiving strong excitatory synaptic input from a connected cortical territory. www.nature.com/scientificreports/ The current results also inform our basic understanding of seizure pathophysiology by replicating results from the animal literature in which destructive, pharmacological, or optogenetic manipulations have been performed. Previous studies in animal models have shown that neuronal firing is less precise when feed-forward inhibition is altered 15,16,33,34 . This lack of action potential timing reliability has been related to several cellular processes. Changes in coupling to the high frequency LFP, as we report here, have been associated with cellular chloride loading in early ictal discharges 35 . However, in Alfonsa et al. such chloride loading was insufficient to induce seizure-like activity alone 35 . Both AMPA and acetylcholine agonists were able to rescue spike timing reliability, and blocking potassium channels increased spike timing reliability in hippocampal slices 15 . The current results support a combination of these mechanisms, suggesting that strong feedforward inhibition, likely mediated by parvalbumin-expressing GABAergic interneurons 19,[36][37][38] , is reflected in narrowband gamma range neuronal synchrony, which is unrelated, both spatially and temporally as we show here, to the BHF increases that correlate with unconstrained pyramidal cell firing and PDS 39 .

Scientific Reports
As with physiological high frequency LFP, the generating mechanisms of HFOs in human epilepsy remain an open area of investigation. Results from animal studies suggest jittered spike timing as a mechanism for converting ripples to fast ripples, or broadening the high frequency spectrum in general, both in hippocampus 7,15,16 and neocortex 35 . Accordingly, we propose that HFOs associated with PDS, while meeting commonly accepted HFO criteria, are not true oscillations, but rather are generated from summated spatiotemporally jittered fast potentials, which are far stronger during PDS than during non-pathological neuronal bursting activity 21 . Such strong, temporally disorganized neuronal firing results in an aperiodic shift in the high frequency spectrum, resulting in broadband high frequency activity (BHF) 2,40 . BHF contrasts with HFOs (e.g. narrowband gamma oscillations reported here) that may be regarded as true oscillations, such as those produced from phasic population firing that exhibit clear, periodic excursions from the overall ( ≈ f −n ) shape of the neural spectrum 1,41 . It is hypothesized that synchronized fast-spiking inhibitory interneuronal activity plays the crucial role of delineating temporal windows for pyramidal cells to fire 30 . Ironically, such true oscillatory narrowband gamma HFOs may therefore occur in brain sites that have not yet been invaded by seizures due to an intact inhibitory restraint mechanism 42,43 . The data reported here therefore support the existence of two distinct mechanisms for pathological high frequency LFP with differing biophysical mechanisms and utilities as biomarkers.

Materials and methods
Subjects and ethics statement. Four human patients with pharmacoresistant epilepsy were included in this study. Despite these relatively small sample sizes of patients and seizures, our statistics are based on the numbers of units, electrodes, and discharges observed in these seizures and are therefore sufficiently powered. The included patients were implanted with standard clinical electrocorticography (ECoG) electrodes (2 patients) or stereo-electroencephalographic (sEEG) electrodes as part of the monitoring prescribed for surgical treatment of their epilepsy. In addition, all four patients were implanted with microelectrodes. ECoG patients received a "Utah-style" microelectrode array (MEA), and stereo EEG patients had Behnke-Fried (B-F) style microelectrodes. MEAs are 4 × 4 mm square arrays with 96 recording microelectrodes (10 × 10 electrodes, with the corner electrodes disconnected) that penetrate 1 mm into layer 4/5 of human association cortex 42 . B-F microelectrodes consist of a bundle of eight microelectrodes that protruded approximately 4 mm from the tips of depth electrodes. Detailed methodology for surgical implantation of these arrays are described in 24 for MEAs and 25 for B-Fs. The MEAs were all implanted into the clinically defined seizure onset zone. Both the New York University and Columbia University Medical Center Institutional Review Boards approved this study and each patient provided informed consent prior to surgery. Research was performed in accordance with the relevant regulations specified by these review boards.
Data collection and preprocessing. Neural data were acquired from each microelectrode at 30 kilosamples per second (0.3 Hz-7.5 kHz bandpass, 16-bit precision, range ± 8 mV) using a Neuroport Neural Signal Processor (Blackrock Microsystems, LLC). ECoG and sEEG signals were acquired using the clinical amplifier (Natus Medical, Inc.) at either 500 or 2000 samples per second (0.5 Hz high pass, low pass set to ¼ sampling rate, 24-bit precision). All ECoG data with sampling rates higher than 500 samples per second were subsequently filtered (4 th order Butterworth) and downsampled to 500 samples per second in order to facilitate comparisons across patients and clinically-relevant sampling rates.
Spectral analysis. Spectrograms of full seizures and of individual discharges were each generated using 5-cycle Morlet wavelet scaleograms using 107 logarithmically-spaced frequency bins from 1 to 200 Hz for ECoG data, and 1 to 850 Hz for MEA data. The absolute values of these scaleograms were normalized by a theoretical spectrum, S, such that.
where f represents the frequency range of the spectrum. These methods were used for both the full duration of the seizure and for each discharge and for both ECoG and MEA recordings.

Dominant frequency analysis.
In order to detect the dominant frequencies of seizure-associated LFP through the duration of the seizure, the aforementioned spectral analysis was carried out on each channel of electrophysiological data from 5 s before any discharge detection to 5 s after the final discharge detection. The maximum frequency in each time bin was then smoothed by a moving average with a width of one eighth of the Per-discharge analyses. In order to understand differences among ictal discharges before recruitment into the ictal core, we first detected each ictal discharge using a thresholding and clustering process. First, all local maxima above one standard deviation in the mean LFP power between 50 and 200 Hz were detected. This method yielded a few false positives due to the noisy nature of LFP, particularly on ECoG. In order to remove false positives, the differences in amplitude and timing of the aforementioned detections were clustered into seven groups using k-means and outlying clusters were removed in a supervised manor until aberrant detections were no longer present. The remaining detection times were used for all further ECoG analyses. We detected discharges on the MEA using a different procedure, which capitalized on the lack of noise in the MEA recordings and match firing rate and LFP measurements across discharges. Specifically, we detected discharges on the Utah array by cross-referencing the normalized maxima of firing rate and instantaneous high gamma amplitude over channels. First, both instantaneous high gamma and firing rate were estimated for each channel and convolved with a Gaussian kernel with 10 ms standard deviation. These signals were normalized by subtraction of their absolute minimum and division by their absolute maximum. Local maxima were detected in both high gamma and firing rate using differentiation. Each maximum exceeding one standard deviation of the firing rate or high gamma amplitude during the seizure was retained as an ictal discharge. Any discharge that occurred within 50 ms of the discharge peak was excluded from further analysis. Discharges that were detected in the BHF signal, but not in the firing rate signal, or vice versa, were also discarded. The resulting data are therefore comprised of co-occuring ictal discharges in the HFA and firing rate signals.

MUA detection and time of recruitment.
Since we previously showed that sorting single units during the recruited period is impossible, here we examine multiunit activity (MUA) recorded on MEAs and microwires. MUA on these electrodes was determined by first zero-phase-lag filtering the data between 300 and 3000 Hz and then finding peaks greater than 3.5 times two thirds of the median absolute value of the filtered data 26 . We refer to the set of times of these peaks on each channel as a multiunit, and firing rate calculated from these times, via convolution with a gaussian kernel, as MUA 44 .
We determined that a multi-unit was recruited into a seizure from its pattern of firing, as previously described 42 . If the multi-unit exhibited a period of intense tonic firing that was flanked in time by burst firing that was phase locked to the dominant seizure rhythm, we considered that multi-unit as recruited after that period of tonic firing 42 . All of the multi-units from the MEA recordings were recruited into the seizure core. Only a single MTL unit from the microwires was recruited into the core. The time of recruitment was operationally defined using methods from 5 . Specifically, a Gaussian kernel with 500 ms width was convolved with the times MUA firing on each channel in order to emphasize slow elements of the signal. The mean peak of this slow firing rate estimate across channels defined the time of seizure recruitment.
Action potential phase coupling. The timing of multiunit action potentials relative to the LFP recorded on the nearest macroelectrode was examined. MUA coupling was determined by extracting the phase from the complex elements of the analytical signal of the bandpass filtered LFP (96th order fir filter) at ±5 Hz above and below the dominant frequency of the LFP, calculated as aforementioned, for each MUA event time. Histograms of action potential times relative to LFP phase (36 bins) were then constructed for each unit and across units for each seizure. We tested for differences in coupling for each seizure using a modified Kuiper test in which a random partition of 2000 action potential times was taken from the cumulative density functions calculated above for both the pre and post-recruitment periods. The test statistic, V, was defined as.
where D + and D − represent the maximum two-tailed distances between the cumulative density functions, as in a standard Kuiper or Kolmogorov-Smirnov test. This procedure was carried out 10,000 times and the test statistic calculated on data from all units in each seizure was compared to the permutation distribution in order to derive a p-value. Kuiper tests in which the test statistic controlled for sample size were carried out on a per-unit basis for MTL neurons using a permutation size of 200 45 . Defining recruited clinical electrodes using phase-locked measures. As a formal method for defining recruited vs unrecruited electrodes, Gaussian Mixture Models with two univariate normal density components were considered. These models were fit using the iterative expectation maximization algorithm. Each of the models converged and exhibited clearly-separated posterior probability distributions.
Per-discharge entropy. In order to quantify the temporal complexity of voltages measured during ictal oscillations, V(t), differential entropy of the signals, h(V), for each discharge was calculated as where p(v) is the probability density function of voltage values for each discharge and n is the number of samples examined for each discharge, where n was 98 samples for macroelectrodes and 296 samples for microelectrodes.
The entropy of the signals in the spectral domain, as in 15,16,31 , were also examined for each ictal discharge recorded on the MEA. The spectral entropy was calculated as www.nature.com/scientificreports/ where p(s) is the distribution of power in each frequency between 100 and 800 Hz. This measure has a maximum at the log of the number of frequency bins used in the spectral estimation. The current study used 50 frequency bins for ECoG and sEEG electrodes and 107 frequency bins for microelectrodes, yielding maximal entropies of 5.64 and 6.71 bits, respectively. This entropy measure was compared with spectral mode and fast ripple index, as calculated in 16 . The spectral mode is the most prevalent frequency in the spectrum between 100 and 800 Hz and the fast ripple index is the proportion of the normalized spectrum above 250 Hz. We added an additional index to these analyses in order to observe similar relationships in macroelectrode recordings. This "ripple index" was defined as the proportion of the normalized spectrum between 80 and 200 Hz in spectra spanning frequencies from 1 to 250 Hz.
Bayesian statistics. Several of the linear mixed effects models (LMM) implemented did not exhibit significant slopes. We therefore sought to quantify the extent to which our beliefs about the hypotheses tested in these models were updated by using Bayesian analysis. We calculated Bayes factors for the LMM fixed effects from the models' Bayesian Information Criteria (BIC) 46 . The Bayes Factors were calculated as where BIC refers to the difference between the BIC for the alternative hypothesis, BIC H1 , and the BIC for the null hypothesis, BIC H0 , defined as where k is the number of model parameters and n is the total sample size. We set k equal to two for the models analyzed in this study, comprising the two LMM parameters: intercept and slope.

Data availability
Matlab code for the measurements and analyses performed herein can be found at https ://githu b.com/ellio thsmi th/seizu reCod es. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.