Delta-gamma phase-amplitude coupling as a biomarker of postictal generalized EEG suppression.

Postictal generalized EEG suppression is the state of suppression of electrical activity at the end of a seizure. Prolongation of this state has been associated with increased risk of sudden unexpected death in epilepsy, making characterization of underlying electrical rhythmic activity during postictal suppression an important step in improving epilepsy treatment. Phase-amplitude coupling in EEG reflects cognitive coding within brain networks and some of those codes highlight epileptic activity; therefore, we hypothesized that there are distinct phase-amplitude coupling features in the postictal suppression state that can provide an improved estimate of this state in the context of patient risk for sudden unexpected death in epilepsy. We used both intracranial and scalp EEG data from eleven patients (six male, five female; age range 21-41 years) containing 25 seizures, to identify frequency dynamics, both in the ictal and postictal EEG suppression states. Cross-frequency coupling analysis identified that during seizures there was a gradual decrease of phase frequency in the coupling between delta (0.5-4 Hz) and gamma (30+ Hz), which was followed by an increased coupling between the phase of 0.5-1.5 Hz signal and amplitude of 30-50 Hz signal in the postictal state as compared to the pre-seizure baseline. This marker was consistent across patients. Then, using these postictal-specific features, an unsupervised state classifier-a hidden Markov model-was able to reliably classify four distinct states of seizure episodes, including a postictal suppression state. Furthermore, a connectome analysis of the postictal suppression states showed increased information flow within the network during postictal suppression states as compared to the pre-seizure baseline, suggesting enhanced network communication. When the same tools were applied to the EEG of an epilepsy patient who died unexpectedly, ictal coupling dynamics disappeared and postictal phase-amplitude coupling remained constant throughout. Overall, our findings suggest that there are active postictal networks, as defined through coupling dynamics that can be used to objectively classify the postictal suppression state; furthermore, in a case study of sudden unexpected death in epilepsy, the network does not show ictal-like phase-amplitude coupling features despite the presence of convulsive seizures, and instead demonstrates activity similar to postictal. The postictal suppression state is a period of elevated network activity as compared to the baseline activity which can provide key insights into the epileptic pathology.

Introduction

Frequently, after an epileptic seizure, the brain enters a postictal state—a transition period to a following baseline interictal state. This postictal state has been correlated with impaired cognition, prolonged confusion and various other co-morbidities (Fisher and Schachter, 2000; Theodore 2010). One particular postictal state frequently follows convulsive seizures—postictal generalized EEG suppression (PGES), which manifests as a reduction or even total cessation of apparent background activity (Lhatoo ; Surges ), and is assumed to represent underlying markedly decreased electrophysiological activity. PGES tends to be present in patients with a history of uncontrollable generalized tonic-clonic seizures (Moseley and DeGiorgio, 2015), particularly those at greatest risk for sudden unexpected death in epilepsy (SUDEP) (DeGiorgio ). In recent years, studies on postictal states focused primarily on PGES since several studies have connected its incidence and duration with SUDEP (Lhatoo ; Kang ).

Since the presence and duration of PGES is correlated with the risk of SUDEP, the ability to predict PGES duration is an important step in identifying at-risk patients. While the relationship between PGES and SUDEP is apparent, conflicting reports have emerged as to how these two phenomena are connected. While one study reported that longer PGES states lead to a higher risk of SUDEP (Lhatoo ), another has shown that SUDEP patients have shorter PGES duration as compared to other epilepsy patients (Kang ). Hence, the standard amplitude and duration definitions of PGES may not be sufficient to address the relationship between the postictal state and SUDEP risk.

One way to analyse epileptiform discharges in brain networks is to decompose the original signal into its constituent oscillations. The coupling between specific neuronal oscillations can be assessed using phase-amplitude cross-frequency coupling (PAC; Tort ) and is an important part of cognition (Lisman and Jensen, 2013). It can also signify how brain communication changes under pathological conditions like Parkinson’s disease (de Hemptinne ), and epilepsy (Guirgis ). In particular, the phase of the delta (0.5–4 Hz) oscillation and amplitude of the high-frequency oscillation (>30 Hz) were previously successful in localizing the epileptogenic zone of patients (Guirgis ), in predicting pre-clinical seizure onset (Jacobs ), and for differentiating seizure types (Amiri ). Therefore, we hypothesized that PAC can better characterize PGES. The objectives were to identify coupled neuronal oscillations specific to postictal states, and to use those PAC features to highlight the active network state during PGES in refractory epilepsy patients and in a case study of a patient with confirmed SUDEP.

Materials and methods

Human EEG recordings

For analysis of PGES and its biomarkers, EEG data were obtained from 12 patients (see Table 1) at Toronto Western Hospital (TWH—Toronto, Canada) and Phramongkutklao Hospital (PH—Bangkok, Thailand). All patients suffered from refractory epilepsy, and one patient died from SUDEP within 2 years of EEG collection. The majority of patients underwent surgical resection to remove a suspected epileptogenic zone. Epileptic surgical intervention outcomes were classified using Engel Class system (Durnford ), with Engel Class 1 (EC1) indicating that the patient was free from seizures following the surgery, Engel Class 3 (EC3)—a worthwhile improvement without abolishing seizures, and Engel Class 4 (EC4)—no worthwhile improvement. To choose an electrode closest to the seizure onset zone, a combination of the clinician-defined region of focal epileptic seizure generation and analyses done previously by our lab (Guirgis ) was used. In total, 26 substantially artefact-free seizure recordings were identified in patients—ten seizures with intracranial EEG (iEEG) recordings, nine seizures with scalp EEG recordings (using the standard 10–20 system), and seven seizures with simultaneous intracranial-scalp recordings. Three patients with five seizures (Patients 9–11) had a postictal state classified as PGES by an expert electroencephalographer. One patient (Patient 12, referred to as SUDEP patient) was an individual who died of SUDEP around 2 years after the recording and did not show any postictal suppression. The other patients had a postictal EEG suppression state (PES) that did not generalize across all channels. The acronym P(G)ES was used in this study to refer to both generalized and non-generalized postictal suppression states. Referential recordings were used in all analyses, with the exception of network connectivity where a bipolar montage was computed to remove common elements for spatial analysis. Recordings used an acquisition reference at FCz, grounded at Fpz; with artefactual channels removed before data analysis. For the SUDEP patient, one iEEG and one scalp EEG recording were obtained. A further summary of patient data is in Table 1.

Table 1

Summary of patient data

Patient no.	Age/ sex	Surgical outcome	MRI findings	Electrode placement	Intracranial EEG		Scalp EEG		Average P(G)ES duration (s)
					No. of recordings	Sampling rate	No. of recordings	Sampling rate
1	36/F	N/A	Abnormal intensity lesion	Left FT	1	2 kHz	0		10
2	28/M	EC3	Non-lesional	Left FT	3	2 kHz	1	200 Hz	19.5
3	21/M	EC1	Cortical dysplasia	Left DF	3	2 kHz	2	256 Hz	21.5
4	41/M	EC1	Mesial temporal sclerosis	Right T	0		2	512 Hz	4.2
5	22/M	EC1	Normal	Left T	1	500 Hz	0		13
6	31/M	EC1	Mesial temporal sclerosis	Right T	1	500 Hz	2	200 Hz	9.8
7	26/F	EC1	Normal	Right F	2	500 Hz	3b	500/256 Hz	11.7
8	36/F	EC1	N/A	N/A	0		1	256 Hz	29.5
9	29/F	EC3	Normal	Bilateral F	1	1000 Hz	0		11
10	32/F	EC4	Normal	Left T	2	500 Hz	2b	500 Hz	39
11	29/M	N/A	Normal	Bilateral F	2	250 Hz	2b	250 Hz	18
12a	30/M	EC1	Normal	Left FT	1	1 kHz	1b	1 kHz

Patient with intractable epilepsy which later died from SUDEP.

Scalp EEG seizures are simultaneous recordings with iEEG with the same sampling rate.

D = dorsolateral; F = frontal; T = temporal.

Seizure onset and termination were marked by experienced clinical electroencephalographers (J.M.C., Y.C. and R.W.), with occasional differences resolved by consensus. A formal comparison for the inter- and intra-rater reliability was not performed. The durations of postictal suppression states were determined visually (with confirmation from at least two EEG specialists), and seizures were selected that had P(G)ES durations of at least 4 s long. Previous studies have allowed for PGES of as short as one second (Surges ); however, we are looking at very low delta frequencies (0.5–1.5 Hz), requiring at least four seconds of P(G)ES.

Patient data were originally filtered with a 0.1 Hz high-pass filter during acquisition and were later pre-processed by removing power line interference using an FIR notch filter at 50 Hz or 60 Hz (depending on the data source) and associated harmonics.

Phase-amplitude CFC analysis

To investigate PAC dynamics between baseline (defined as the quiet period before the seizure without visible artefacts, spikes or other activity at least 60 s before seizure onset) and the postictal states, wavelet coefficients corresponding to frequency bands of interest were extracted from patient data using the continuous wavelet transform. We used the complex Morlet wavelet with a mother wavelet of 0.8125 Hz centre frequency and 5 Hz bandwidth, as these parameters have been previously found to work best with EEG data for epilepsy (Guirgis ).

The strength of PAC is measured by the cross-frequency coupling (CFC) index following the protocol established by Tort and adapted by our lab (Guirgis ; Colic ). PAC is a normalized measure of how much preference the amplitude of the higher frequency component has for particular phases of the lower frequency component. In short, wavelet coefficients obtained from continuous wavelet transform of the form were used to recreate these two time-series (Colic ), with the amplitude envelope signal F covering the range f of 30–100 Hz in 1 Hz increments (except in select iEEG cases, where higher sampling rate allowed to expand the range to 30–200 Hz), and the instantaneous phase signal F covering the range f of 0.5–5 Hz in 0.1 Hz increments:

F was binned into 20° windows—for a total of 18 bins—and F within each bin was averaged. A discrete probability density function was then calculated for each individual averaged amplitude envelope (denoted as : and the PAC was effectively the normalized Kullback–Leibler distance of this distribution from the uniform distribution (where N = 18):

While in the case of iEEG recordings of EC1 patients it is fairly easy to establish the location of the seizure onset electrodes (since they are captured by the electrode grid and lie within the successful resection area), the task becomes significantly harder in scalp EEG recordings. In order to reduce the effects of potential artefacts and analyse phenomena in a more spatially generalized manner, we used a global PAC measure, inspired by Jacobs . To calculate global PAC with a reduced influence of potential artefactual channels, a median for each phase-amplitude combination was obtained. For comparison, in addition to scalp EEG, global PAC was also calculated for iEEG.

To identify changes in PAC between the baseline state and the P(G)ES state, we binned the PAC values into 63 bins, with each bin having a value of the median PAC across 0.5 Hz of phase and 10 Hz of amplitude. The absolute difference between baseline state and P(G)ES state was calculated for each bin, and regions of increased PAC in both states were identified as those having PAC greater than 3 dB (∼0.707) of the maximal coupling during P(G)ES (a common ‘half-maximum-power’ threshold).

Brain network connectivity

To assess changes in the network connectivity from baseline to P(G)ES, network analysis was performed on the iEEG recordings. Noisy or artefactual channels were removed from the analysis and all iEEG traces were z-score normalized in order to reduce the influence of amplitude differences on the network connectivity estimation (van Mierlo ).

A time-varying directed transfer function (DTF) (Astolfi ) was computed following the implementation described in detail by Omidvarnia et al. (Omidvarnia ) using a multivariate adaptive autoregressive model (MVAAR) based on the Kalman filter approach (Vidaurre ). Model parameters were found in the literature of similar studies (Staljanssens ; van Mierlo ): model order of 8, an update coefficient of 10−3, and a downsampling factor of 10. The time-varying DTF was computed in the frequency ranges corresponding to the 3 dB peak global PAC specific to each window, corresponding to f ∈ (0.5, 1.5 Hz) in 0.1 Hz intervals, and various subsets of f ∈ (30, 70 Hz) in 1 Hz intervals depending on the window (ranges shown in Supplementary Fig. 3 and summarized in Supplementary Table 1). Thus, the DTF for a given window can be represented as 4D matrix (NChannel × NChannel × NTime × NFrequency), where the first channel dimensions are information sinks and the second dimension are information sources, and NFrequency = NfH + NfL depending on the global PAC for the window in question. To better represent PAC effects, we took the mean of the DTF coefficients at each (f, f), which results in the 5D DTF matrix DTFCC = (NChannel × NChannel × NTime × NfH × NfL).

In order to properly compare the DTFCC across states and subjects, network characteristics were summarized using the graph theory measure of closeness centrality (CC). Here, CC is computed at each (f) pairing for each event as the inverse of the sum of the shortest distances from the jth source channel to all other channel sinks, using the Floyd–Warshall algorithm (Rubinov and Sporns, 2010). For each baseline versus P(G)ES state comparison, we control for the differences in window duration, frequency range and number of channels by normalizing to the baseline state median.

Assessing postictal state duration with hidden Markov models

To classify the post-ictal state, a first-order hidden Markov model (HMM) was applied as per Breton . We adapted the features for application to a 4-state model with two Gaussians used to fit each state (as depicted in Supplementary Fig. 10). To make the feature set used for the classifier, first, the continuous wavelet transform was computed using 30–50 Hz (in 1 Hz intervals) for the high-frequency amplitude signal and 0.5–1.5 Hz for the low-frequency phase signal (in 0.5 Hz intervals). Second, the phase (φ) signal was combined with the amplitude (A) signal over time (t) using the following formula (Li ): with n1 = 1 … N1, n2 = 1 … N2, where N1 = 21 and N2 = 3 for this study. The envelope of the resultant nth feature signal was averaged on two second intervals, with no overlap. Then, the envelope signal was normalized to the maximum and minimum for the entire trace for each feature independently. Training was performed on the ictal and post-ictal period of Patient 4 seizure 1 data. Testing was performed on the entire signal using six second windows with two second moving windows for all patients.

Association of P(G)ES duration with ictal PAC

In order to establish the correlation between P(G)ES duration and ictal PAC dynamics, PAC was extracted from onset iEEG electrodes following the algorithm described in the earlier section [see eqs (1)–(5)]. Specifically, PAC was calculated on consecutive four second windows, with no overlap, throughout the seizure. Maximal PAC measures were then tracked throughout the seizure, and the phase of the PAC measure was used with a maximal low frequency of 5 Hz to obtain the dominant PAC curve. To allow for better comparison across patients, phase frequency was normalized by subtracting the minimal value of the phase of maximal PAC, setting the baseline to 0. To quantify the rate of dominant PAC rise and decay, a single measure, β, was obtained by fitting the equation below to the extracted dominant PAC curve: where H was fitted to the height of the dominant PAC curve. The β-fit has been proposed in Zalay ) to characterize system excitability, and has been favourably compared to other excitability measures—such as Teager energy—and used in evaluating spontaneous epileptiform discharges (Grigorovsky and Bardakjian, 2018).

Statistical analysis

A Pearson coefficient was calculated to measure the degree of linear relationship between the HMM-defined PGES-like states and visually identified P(G)ES durations. HMM-derived ictal and postictal state duration distributions were also fitted using a gamma function (gamfit in the MATLAB package): where is the shape parameter and is the rate parameter. When the shape parameter is equal to one, gamma distribution collapses to an exponential distribution, suggesting that the underlying process is a Poisson process. When, however, the shape parameter is larger than one, the process is non-Poisson in nature and instead is time-dependent. When comparing the SUDEP case study patient to the population of other patients, S2 (ictal-like) state durations were compared using both Vysochanskij–Petunin inequality and Chebyshev inequality. Chebyshev inequality is a conservative approach to a non-parametric estimation of confidence interval, while Vysochanskij–Petunin inequality assumes a unimodal distribution. This established lower bounds for the probability that the SUDEP patient’s S2 duration belonged to the general patient population:

To assess whether the area of the regions of increased PAC (calculated as the number of phase-amplitude combination above the 3 dB threshold, or pixels) across the channels and patients varies significantly, a Shapiro–Wilk test was first performed to determine whether the data were distributed normally or not. If it was, a two-sample t-test was used to test for significant change in the average area; if, however, the data were not normally distributed, Wilcoxon rank sum test was used to test for significant change in the median area. A similar approach for significance calculation was used for connectivity analysis. The method of surrogate data testing was adapted from Guirgis and Breton . A 200-surrogate test with 5% significance threshold was used on a sample patients’ iEEG data, to determine whether the PAC observed was significant.

Data availability

The data that support the findings of this study are available on reasonable request from the corresponding author. The data are not publicly available due to institutional restrictions associated with original data acquisition protocols.

Results

Phase-amplitude coupling in P(G)ES

To identify the coupled neuronal oscillations specific to postictal suppression states, we first looked at the PAC features that are present during PGES. As PGES is a generalized phenomenon, a global PAC analysis was used to compare the PGES state to the baseline state (Fig. 1A). Increased coupling between the phase of 0.5–1 Hz and amplitude of 30–40 Hz is observed in both scalp and iEEG (Fig. 1B). Both iEEG and scalp EEG showed a significant (P < 0.0005) increase in coupling during PGES as compared to pre-seizure baseline, when measured by the median area of increased coupling across channels (i.e. the number of individual phase-amplitude combinations over the 3 dB threshold; Fig. 1C). A computation of PAC using surrogate data confirmed that this identified PGES marker was significant, and not due to noise present in the recording (Supplementary Fig. 2). While sampling rates vary between patients based on the data source, we found no distinguishable difference in PAC analyses and elected to use the original data sampling.

Figure 1

Global PAC of PGES patient. (A) Scalp EEG traces of a baseline, seizure onset, seizure offset and a following PGES state in patient 10 seizure 2—seizure duration of 88 s, baseline 100 s before onset, PGES state begins 15 s after the seizure. (B) Global PAC comparison between baseline and PGES states in both intracranial and scalp EEGs, as well as the area of increased (>3 dB) PAC. (C) Median area of increased coupling—in number of phase-amplitude combinations—during baseline and PGES states in both intracranial and scalp EEG recordings for patient 10 (*Wilcoxon rank sum test, P < 0.0005). Red lines show 25% and 75% quartiles. Reference electrode—FCz.

A similar comparison between the baseline and non-generalized PES was done for Patient 7, Seizure 2 (Fig. 2). In both scalp and iEEG, we observed an increase in PAC strength during the PES state within the same ranges as in the PGES case (Fig. 2C), which was confirmed by the Wilcoxon rank test of the median area (Fig. 2D).

Figure 2

Global PAC of PES patient. (A) Scalp EEG traces of a baseline, seizure onset, seizure offset and a following PES state in Patient 7 Seizure 2—seizure duration of 65 s, baseline 900 s before onset. (B) Zoomed in traces of the scalp EEG showing the actual PES state. Note lack of generalization of the state across all of the electrodes. (C) Global PAC comparison between baseline and PES states in both intracranial and scalp EEGs, as well as the area of increased (>3 dB) PAC. (D) Median area of increased coupling—in number of phase-amplitude combinations—during baseline and PES states in both intracranial and scalp EEG recordings for Patient 7 (*Wilcoxon rank sum test, P < 0.005). Red lines show 25% and 75% quartiles. Reference electrode—FCz.

These postictal PAC features are not patient-specific—across patients with both PGES and PES, P(G)ES states demonstrated similar ranges of increased PAC (Fig. 3A and C) and a similar increase in the median area of increased PAC (Fig. 3B and D). Since the ranges and the areas were comparable for both PGES and PES patient groups, results were pooled into an overall comparison for P(G)ES (Fig. 3E and F). The sum of the area was significantly larger in the P(G)ES state (P < 0.005, Wilcoxon test). Note that, while the area of increased PAC shows the binned ranges of interest, the median area represented the number of individual phase-amplitude frequency combinations—without binning—above the 3 dB threshold (Fig. 3B, D and F). Overall, these PAC results show that there is a significant increase in global PAC in P(G)ES, regardless of whether postictal suppression is generalized or not.

Figure 3

Comparison of global PAC between baseline and P(G)ES states across patients. (A) Median CFC across all of the patients with PGES (n = 3) within iEEG and scalp EEG, highlighting the increased coupling in the PGES state using binned frequency ranges. (B) Median area of increased coupling—in number of phase-amplitude combinations—during baseline and PGES states in both intracranial and scalp EEG recordings for patients with PGES. Red lines show 25% and 75% quartiles. (C) Similar to part (A), median CFC across all of the patients with PES (n = 8) in iEEG and scalp EEG. (D) Median area of increased coupling during baseline and PES states in both intracranial and scalp EEG across patients with PES. Red lines show 25% and 75% quartiles. (E) Median CFC of both (A) and (C) pools joined together (n = 11). (F) Median area of increased coupling during baseline and postictal states across all patients. In parts (B), (D), (F), postictal states have significant increase in mean area—Wilcoxon rank sum test, P < 0.005.

Brain network connectivity using P(G)ES-specific PAC features

After establishing the P(G)ES-specific PAC frequencies, we wanted to investigate what effect these frequency ranges have on network connectivity. Network-level connectivity analysis showed significantly elevated CC during the postictal state within the same PAC frequency bands as those present during P(G)ES. To identify network connectivity, the shortest-path weight between iEEG channel pairs was measured, and from this weight, the CC of the networks was computed using 0.5–1.5 Hz coupled to 30–50 Hz features (Fig. 4, Supplementary Fig. 5), as the average of all nodes. In the case of PGES, we show a representative shortest-path network between the channel with the median-level CC during the PGES state (channel 3) and all other channels (Patient 10, seizure 2, Fig. 4A; iEEG channel locations are shown in Supplementary Fig. 4). There was a decrease in the connection path length between channel pairs during the postictal state as compared to baseline. This decrease summed to an overall greater CC during PGES. In comparison, in a case of PES (Patient 7, Seizure 2; Fig. 4C channel 13), the same increase in CC occurred (Fig. 4D; electrode montage shown in Supplementary Fig. 4). Across all nine subjects with iEEG recordings (Fig. 4E), the P(G)ES CC was found to be significantly elevated above the baseline (P < 0.0005, Wilcoxon rank sum test; for individual patient comparison, see Supplementary Fig. 5). Hence, the connectivity of the particular network discussed above is enhanced in the post-ictal state, and this occurs both in the generalized and non-generalized EEG suppression cases.

Figure 4

CC of iEEG networks is elevated for P(G)ES state versus baseline. (A) Example network shown for baseline versus PGES state for Patient 10, Seizure 2. Network consists of outward shortest-path connections for the PES median channel (3), normalized to the PGES state. Edge length is scaled by physical distance in the iEEG bipolar montage, and colour and line width is scaled by the connection strength. (B) Summary of networks for all channel pairings for P10, Seizure 2, using cross-channel mean of CC. Here, it is shown as a function of PAC frequency ranges and normalized to the median of the baseline. (C, D) Same as (A) and (B) but for Patient 7, Seizure 2. (E) CC is generally elevated for the PES state when considering the lumped normalized data from all subjects (Nsubjects = 9, Nevents = 16). (*Wilcoxon rank sum test, P < 0.0005).

Assessing postictal state duration with HMMs

Investigating whether identified PAC frequency ranges can be used to classify different states within a seizure episode, HMM were used with iEEG recordings to better understand if, and how, postictal suppression states compare to the statistical properties of ictal states. Classification resulted in a clear delineation of the iEEG into four states (including seizure state S2 and post-seizure state S3) (Fig. 5A). To validate this classification, the model output was compared to visually identified epochs of the iEEG. First, the durations of the HMM-classified PGES-like state were compared with visually determined P(G)ES (Fig. 5B). We found a significant positive correlation of 0.77 (P = 0.009) between the HMM classified and the visually identified P(G)ES, suggesting that the P(G)ES state was well classified using the HMM. Then, classification was performed on visually identified interictal data from an iEEG trace, which designated state S4 as interictal (Supplementary Fig. 6).

Figure 5

State classification of iEEG recordings using a four-state HMM. (A) State classification of an example iEEG trace from Patient 2 (reference electrode—FCz), where S2 is a seizure-like state, and S3 is PGES-like state. (B) Comparison of HMM-driven S3 state classification to the visually estimated PGES durations, with the overall correlation of 0.77 (P = 0.009). (C) Histogram of S2 durations, with a gamma function fitted to the data. The shape parameter of the gamma function—alpha—was calculated to be 3.54 (95% CI 1.59–7.89), which is consistent with the average alpha for seizure duration distributions from Suffczynski of 3.03 and from Bauer of 2.66. (D) Histogram of S3 durations, with gamma function fitted to the data with alpha parameter calculated to be 2.21 (95% CI 1.01–4.82) which is consistent with alpha for PGES duration distributions from Bauer of 1.54. This can also be compared to alpha value of 1.83 (95% CI 1.04–3.20) obtained from distribution of visually estimated PGES durations (not shown). These shape parameters suggest that the transitions (to and from seizure states) occur not according to Poisson process, but rather that the probability of transition varies with time. For both sections (C) and (D), n = 11.

Then, the statistical properties of the ictal and PGES states were compared by measuring the shape parameters of the duration of the HMM-identified states. HMM analysis was performed on all patient iEEG data, and collated results are shown in Fig. 5C and D for S2 and S3 durations, respectively. A gamma function was fitted to each duration distribution; the shape parameter, alpha, was calculated as 3.54 (95% confidence interval 1.59–7.89) for S2 (Fig. 5B), and 2.21 (CI: 1.01–4.82) for S3 (Fig. 5C). These results validated the existence of the P(G)ES state, and identified similar statistical processes underlying the ictal and postictal state.

Given that the gamma function hinted at similar statistical properties between HMM-classified ictal- and P(G)ES-like states, our aim was to find an association between PAC features and the P(G)ES state. Some of the CFC changes used by the HMM appear not only between P(G)ES and baseline states, but within the ictal state as well. Figure 6B highlights phase-amplitude coupling dynamics during the ictal and PGES states, compared to the baseline state in patient 10, seizure 2 (Fig. 6A). Similar ictal PAC dynamics were also observed in patients with PES (Fig. 6C and D). In both cases, PAC features during P(G)ES and baseline states are not clearly visible on the scale of ictal coupling since ictal PAC is significantly elevated as compared to P(G)ES PAC (baseline and P(G)ES co-modulograms obtained from iEEG data are shown separately in Supplementary Fig. 1 with PAC). These data showed variability in PAC, shifting from ∼4–5 Hz to ∼1–2 Hz of the maximal coupled frequency of the low-frequency phase signal, over the course of the ictal state.

Figure 6

Ictal phase-amplitude coupling dynamics in the onset electrode of P(G)ES patients. (A) iEEG trace from PGES patient 10 (electrode LAT 2) with seizure occurring at 9845–10 033 s. (B) Four second samples of EEG trace with concomitant CFC co-modulograms of different points in the ictal state, baseline state and PGES state, plotted on the same scale. Clear ictal PAC dynamics appear, showing a phase slowing from theta to low delta frequency ranges throughout the seizure. (C) iEEG trace from PES Patient 7 (electrode RPF 6) with seizure occurring at 1038–1103. (D) Similar analysis to part (B) showing the same theta to low delta shift in frequencies. Reference electrode—FCz.

To further establish a connection between the ictal PAC variability and the P(G)ES duration, we applied a β-fit to the frequency of the phase part of the maximal ictal PAC. From empirical observations of the results, for a proper β-fit, the variability in the frequency of the phase signal during the seizure had to be at least 2 Hz and postictal state needed to be <20% of the duration of the preceding seizure. Conversely, if the duration of the postictal state exceeded 40% of the seizure duration, ictal PAC variability flattened and no reliable β-fit could be achieved. Under these constraints, there was a strong correlation between the inverse of the β value and the HMM-defined postictal state duration with Pearson coefficient of 0.98 (P < 0.001; see Fig. 7B; for sensitivity analysis see Supplementary Fig. 7), suggesting a link between relatively long postictal state durations and low ictal PAC variability.

Figure 7

Correlation between frequency of the phase signal of PAC during seizure and postictal state duration. (A) An example patient showing the evolution of normalized phase frequency of the dominant PAC during and right after the seizure, with the red bar showing seizure offset and the blue line showing the fitted beta function. (B) Correlation between the inverse of the beta value and HMM-based postictal state durations, with the dotted line showing trendline.

Note that for both Figs 2C and 6D, the broader f range in the PAC did highlight additional activity in the 100–200 Hz range; however, it did not impact the final analyses (significant PAC during the P(G)ES and the frequency of the phase of the dominant PAC, respectively), and so was omitted from other patients.

A SUDEP patient case study

To see how SUDEP patients differ from our cohort (Patients 1–11), one SUDEP patient with concurrent scalp and iEEG recordings was used (Fig. 8). In this SUDEP patient, there was an absence of the PAC variability as previously observed in the other patients. Specifically, the PGES-like PAC features appeared throughout the recording, even though a seizure had been identified by an electroencephalographer (Fig. 8A and B; Supplementary Fig. 9 for zoomed in intracranial traces). When the same HMM from the previous section was tested on the iEEG trace of the SUDEP patient, it identified primarily S3 (PGES-like) and S4 (interictal-like) states throughout the recording outside of the seizure and primarily S1 state during the seizure. Meanwhile, the S2 (ictal-like) state was largely absent (Fig. 8C) with only a short duration at the beginning of the seizure.

Figure 8

EEG recording analysis from a SUDEP patient. (A) iEEG recordings from the SUDEP patient with seizure onset identified by a red line and seizure offset by a blue line. Below are the corresponding global PAC throughout the recording. These results highlight the relatively low variability in the PAC, compared to coupling dynamics shown in Figs 1 and 2 and the presence of P(G)ES-like PAC features throughout the recording. (B) Scalp EEG recordings (simultaneous with iEEG in (A). Below are the corresponding global PAC throughout the recording. Note the predominance of the same PAC marker throughout the recording as in (A). (C) Single iEEG trace (from LOF1 electrode) centred around the seizure. This iEEG trace is classified using an HMM. Reference electrode—FCz.

To measure the probability that the S2 state of this SUDEP patient belonged to the general patient population, a Chebyshev inequality was used, establishing a 95% interval for ictal state durations: 28.2–77.6 s, while Vysochanskij–Petunin inequality gave a 99% confidence interval of 11.1–90.8 s. Hence, the SUDEP patient S2 duration of 6 s most likely did not belong to that of the general population of epilepsy patients (P = 0.009). Overall, this case study suggested that P(G)ES-like PAC features dominate the network activity of patient at risk of SUDEP, unlike the activity analysed in the general patient cohort (Patients 1–11).

Discussion

We used phase-amplitude CFC features from ictal and postictal states to show that (i) the P(G)ES state has a distinct PAC signature, compared to the baseline state, regardless of whether it is generalized or not; (ii) this PAC signature can be used with HMM to identify a similar postictal state as seen by visual inspection, and with graph-theoretic measures to highlight changing network connectivity between P(G)ES state and the baseline; and (iii) in the case of the SUDEP patient, PAC features and HMM analysis suggest that the network largely stays in a postictal PGES-like state, regardless of the EEG amplitude measure.

The role of delta-gamma PAC in epileptic state classification

Epilepsy has been described as a pathological manifestation of a multistate oscillatory system (Kalitzin ), and the coupling between neural oscillations was an obvious target for investigation. We identified 0.5–1.5 Hz and 30–50 Hz coupled signals as markedly elevated in the PGES state, as compared to baseline. Furthermore, this marker was largely consistent in the postictal suppression state regardless of whether or not the suppression was generalized across channels. Compared to a previous study, which highlighted delta-gamma coupling during seizure termination (Guirgis ), we have identified a narrower range of delta-gamma features that persist within the P(G)ES state. Furthermore, this range was able to resolve four states of the epileptic cortex, as was previously not obtainable using a broader PAC frequency range (Guirgis ). We have also demonstrated that the PAC features observed in the iEEG can be extended to analysis of scalp EEG signals. Given that as few as 5% of patients with intractable epilepsy are eligible for surgery (Ryvlin and Rheims, 2008), and scalp EEG is much less invasive than iEEG, these PAC features have the potential for use in a broader epilepsy population. Moreover, these PAC features can be used for identification and classification of PGES and other postictal states, which can improve tasks such as automated detection of PGES using machine learning algorithms (Theeranaew ).

Activity during P(G)ES

It is important to distinguish ictal and non-ictal PAC. While PAC during P(G)ES is significantly smaller than during a seizure, it is larger than a state before the seizure (see Figs 1–3). Network connectivity showed increased CC in the P(G)ES state compared to baseline suggesting that, despite having decreased voltage amplitude, PGES state is more active than baseline. Yet, despite these results an earlier study argued that there exists a ‘neuronal emergency brake’ mechanism that prevents status epilepticus and when ‘activated too strongly or persistently, PGES occurs’ (Bauer ), suggesting that the PGES is a result of a dampened network activity. In line with our results, a recent modelling study, demonstrated that seizure cessation and an overall suppressed state was achieved by increasing the coupling of oscillatory units, rather than reducing their activity (Farah ). These data provide evidence that the PGES state is not characterized by reduced activity, but rather by activity that is different from either ictal or interictal states.

Pathophysiological mechanisms of PGES

While the mechanisms behind PGES generation are unclear, one hypothesis is that PGES relates to increased neuronal inhibition at and after seizure termination, and that longer PGES occurs after a more severe seizure (Rajakulendran and Nashef, 2015). This is supported by a recent modelling study which showed that the excitation/inhibition balance drove the ictal frequency and coupling dynamics (Liu ). Several modelling papers also showed that microglial effects were responsible for increasing the time between spontaneous epileptiform discharges (Grigorovsky and Bardakjian, 2018; Farah ), suggesting that microglial pruning could be involved in longer PGES durations. This is supported by experimental evidence that ictal activity can enhance motility and drive alterations in microglial physiology (Sepulveda-Rodriguez ), and that seizure activity increases microglial process numbers, suggesting that under certain conditions microglia have neuroprotective effects in the epileptic brain (Eyo ).

Furthermore, similar PAC features appeared in a computational model of neuroglial networks (Grigorovsky and Bardakjian, 2018), suggesting that glial factors play an important role in hyperexcitability, and dysfunction in excitability pathways leads to an emergence of PAC features similar to those seen in epileptic patients.

Ictal PAC features and P(G)ES

Both HMM and β-fit were used to investigate the P(G)ES state features and duration, and how they relate to ictal PAC. Using HMM, S2 (ictal-like state) and S3 (PGES-like state) were classified, and alpha values from gamma distributions of these state durations were consistent with previous studies (Suffczynski ; Bauer ). For ictal-like states, previous studies found alpha in seizure duration distributions in humans to be, on average, 3.03 (Suffczynski ) and 2.66 (Bauer ). For S3, our alpha value (2.21) was qualitatively similar to the earlier study which used a larger patient cohort with 48 seizures and arrived at the alpha of 1.54 (CI 1.01–2.33) (Bauer ). These values can be further compared to the alpha obtained from the distribution of visually estimated P(G)ES durations in our study—1.83 (CI 1.04–3.20; not shown). These findings suggest that both transitions to and from S2 and S3 states (ictal- and postictal-like states, respectively) occur not according to a Poisson process, but rather the probability of transition varies with time (Suffczynski ). Despite the lower number of seizures in our study than in previous studies, the general trend in ictal and postictal state distributions are consistent.

Furthermore, the β-fit of the phase frequency of dominant PAC during seizures showed a strong correlation to the P(G)ES duration (Pearson coefficient of 0.98, P < 0.001). While this P(G)ES duration marker does not appear to apply universally—as some of the patients did not show PAC dynamics necessary for a β-fit—in the subset of patients where it does apply, these PAC features show a higher correlation with P(G)ES duration compared to the more broadly applicable feature developed by Bauer .

These dominant PAC features (as seen in Fig. 6), and their association with the P(G)ES state duration suggest that rather than two distinct and unrelated states, ictal and postictal states show a continuation of activity, but with a reduction in the frequency of phase—from high-delta to low-delta. This link between ictal and postictal states is further supported by the gamma distribution analysis which suggested similar statistical features of S2 and S3 states.

Association of PGES duration with SUDEP

The nature of the relationship between PGES and SUDEP is still controversial—as mentioned earlier, while (Lhatoo ) reported that increased PGES duration (after a 50-s threshold) correlated with higher chance of SUDEP, (Kang ) found that SUDEP patients had statistically shorter PGES duration, and (Surges ) concluded that neither presence nor duration of PGES state was associated with SUDEP. Moreover, another study found that the presence of a long PGES state was not consistent in patients with multiple seizures, suggesting that it might not be a reliable clinical tool for SUDEP risk assessment (Lamberts ). It is possible that this discrepancy in research findings is due to different aetiologies of seizures and associated PGES states. Another reason could be, that the conventional way of identifying PGES—by visual assessment or using the EEG amplitude threshold—highlights a subset of postictal states. This is supported by our SUDEP case study, where despite the EEG being above the PGES threshold, PAC features identified in P(G)ES states dominated the recording.

These networks might be dynamic—in a few instances, dominant PAC regions differed between baseline and PGES states—however, they have equally strong activity. These states were further distinguished using a four-state HMM, which clearly separated the PGES-like S3 state from ictal and other states. Furthermore, when the same model was applied to the iEEG recording of the SUDEP patient, it identified S3 state to be present both before and after the seizure, even when no suppression was identified. Because the S3 state is confined to after a seizure in the general patient cohort, whereas in the SUDEP patient, the S3 state is present during seemingly interictal periods, we suggest that when measuring ‘PGES’ it is critical to consider that the markers that distinguish this state may begin to invade the ‘normal’ brain state when epilepsy severity increases. Further support to this proposed network state change can be seen in Supplementary Fig. 8, where a separate SUDEP patient with generalized tonic-clonic seizure and identified PGES state is briefly analysed and the same PGES-like biomarker is present throughout most of the recording (see Supplementary Fig. 8). This idea, then, may lead to future tests of this potential novel biomarker for susceptibility to SUDEP.

Limitations and future directions

Our work is limited by the sample size of our patient population and is prone to artefacts in EEG or outliers among the patients. To account for that, we relied on median measures as opposed to averages; however, a larger collection of patient data would potentially reduce variability in the results. Our analysis—especially with scalp EEG data—was also constrained by the range of usable frequency bands. On one hand, due to the high-pass filter during data acquisition, extremely low frequencies (on the order of 0.01–0.3 Hz) could not be reliably used for coupling analysis. Given that 0.3 Hz coupling has been shown to increase in the postictal state (Lundstrom ), and 0.2 Hz is a marker in paediatric epilepsy (Nariai ), using lower high-pass filtering cut-off can shed more light on PGES and SUDEP features. On the other hand, sampling frequency provided an upper bound on the HFO range that could be used for PAC calculation, whereas very high-frequency oscillations (VHFOs; 600–2000 Hz) (Jacobs ; Brázdil ) could be combined with low-delta as another biomarker for postictal states.

Furthermore, while the data Though Lundstrom found that the coupling strength between 0.3–1 Hz and 20–50 Hz increased in the postictal state—providing more support to our results (Lundstrom ) —they also saw some overlap in the broadband delta and HFO power and CFC strength between awake state, postictal state and sleep state. Also, given that previously obtained distributions of the ictal epochs are different depending on whether the EEG was obtained during the day or night (Suffczynski ), future studies may consider the interaction between sleep and postictal states.

It is prudent to highlight the fact that a long time after EEG collection, patient 10 died of what could have been possible SUDEP. Since more than nine years have elapsed between the EEG recordings and death—unlike the SUDEP case study, which was within 2 years—the patient was treated as part of the overall cohort population (compare, for example, with Fig. 8 and Supplementary Fig. 8). However, this situation highlights the need for longitudinal SUDEP analysis to understand at which point and at what rate PGES-like biomarkers begin to dominate EEG recordings.

While our SUDEP case study is an important comparison to the ictal and postictal PAC dynamics in the rest of the patients, a more in-depth SUDEP patient analysis is necessary to draw more general conclusions. This includes not only having more SUDEP patients but also incorporating EKG and respiration data into analysis, as recently identified potential contributors to SUDEP susceptibility (Rajakulendran and Nashef, 2015; Lertwittayanon ). Furthermore, since the exponential fit used by Bauer et al. has been used as an input (together with video recording data) to PGES and, indirectly, to a SUDEP classification algorithm using support vector machines (van Beurden ), PAC markers identified in this study can be used as alternate features for a similar machine learning approach.

Supplementary material

Supplementary material is available at Brain Communications online.

Funding

This study was funded in part by the Natural Sciences and Engineering Research Council (NSERC), the Canadian Institutes of Health Research (CIHR), and the Ontario Brain Institute (OBI).

Competing interests

The authors report no competing interests.

Click here for additional data file.