Scalp EEG recordings for 14 non-SUDEP (with 77 peri-ictal segments) and 8 SUDEP (with 25 peri-ictal segments) patients were provided through a consortium formed by the Toronto Western Hospital, the NYU Comprehensive Epilepsy Center, and the Phramongkutklao Royal Army Hospital (Table 1). Non-SUDEP patients had focal (temporal or extratemporal lobe) epilepsy, were resistant to anti-seizure medications and were undergoing presurgical evaluation. Ictal segments were marked by board-certified neurologists and electroencephalographers. The institutional review boards of the consortium approved the study protocol and all patients gave informed consent.

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 a finite impulse response (FIR) notch filter at 50 Hz or 60 Hz (data centre location dependent) and associated harmonics. Recordings used an acquisition reference at FCz, grounded at Fpz. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto (Loken et al., 2010; Ponce et al., 2019).

Wavelet phase coherence (WPC) was computed between pairs of scalp EEG electrodes. The phase of different frequency bands was extracted through the complex wavelet transform (Cotic et al., 2015). The Morlet complex wavelet transform was used with a mother wavelet of central frequency of 0.8125 Hz and bandwidth of 5 Hz, as previously used on EEG data (Grigorovsky et al., 2020). The relative phase difference was obtained using Δ ϕ ( s ,   τ ) = tan − 1 ( W 1 ∗ ( s , τ ) W 2 ( s , τ ) − W 1 ( s , τ ) W 2 ∗ ( s , τ ) W 1 ( s , τ ) W 2 ( s , τ ) −   W 1 ∗ ( s , τ ) W 2 ∗ ( s , τ ) ) (1) where " W " is wavelet coefficient, " W ∗ " is the complex conjugate, “s” is the scaling coefficient and " τ " is the time shift. The phase coherence between two electrodes is computed over a time window ( N · Δ t ) expressed as an integer multiple N of the sampling period Δ t . The phase coherence is defined as: ρ ( s ,   τ ) = | 〈 e j Δ ϕ ( s , τ ) 〉 | = 1 N + 1 ∑ k = − N / 2 N / 2 e j Δ ϕ ( s , τ + k Δ t ) (2)

WPC was applied to each wavelet central frequency in the delta (0.5–2 Hz) and HFO (80–120 Hz) ranges, in increments on a logarithmic base two scale. Wavelet frequency scales took the form of 2 k where " k " ranged from − 2.0 to 1.0 for the delta range and 6.3 to 6.9 for the HFO range in 0.1 incremental steps. The WPC window size was proportional to 8 cycles for each frequency.

The undirected connectivity matrices were computed at 1 s intervals by assigning each edge corresponding to a pair of scalp electrodes to the WPC averaged over the frequency range (delta or HFO) and over 1 s temporal windows. The edges between electrodes are undirected, yielding to symmetrical connectivity matrices. e i j ( t ) = e j i ( t ) = 〈 ρ i j ( s , τ ) 〉   |     e i j ∈ E ,   ( i , j ) ∈ V 2 (3) where " e i j ( t ) " and " ρ i j ( s , τ ) " is the edge and WPC between node " i " and " j " , " E " and " V " is the set of all edges and the set of vertices respectively in the graph.

The undirected connectivity matrices were derived from delta–HFO WPC. The degree of each node was computed for both the delta and HFO WPC connectivity adjacency matrices to construct the directed delta-HFO network. Edges between pairs of nodes i and j were computed from the product of the delta degree of node i and HFO degree of node j . This measure represents information flow across the brain pertaining to the coexistence of simultaneous cohered delta and cohered HFO regions, which are two frequency ranges associated with seizure activity (Guirgis et al., 2015; Grigorovsky et al., 2020). The nature of this simultaneous coexistence may or may not be associated (Figure 2) with classical cross-frequency coupling (Tort et al., 2008; Canolty and Knight 2010; Stankovski et al., 2017). Additionally, this measure can identify phase-amplitude cross-frequency coupling (Supplementary Figure S1), where low frequency oscillations coordinate large scale activity driving local HFOs.

Global coherence was computed as a measure for the undirected delta and HFO connectivity networks. First the eigenvalues of the connectivity matrix were computed and sorted. The global coherence is a ratio between the largest eigenvalue to the sum of all eigenvalues and has been previously used to analyze spatiotemporal EEG dynamics (Cimenser et al., 2011). C G l o b a l = λ max   ∑ ​ ⁡ λ (4) The temporal mean of the global coherence was used to assess group differences.

Two network measures were computed: the clustering coefficient and the global efficiency. Together, these measures provide a description of the connection topology in the brain network.

The clustering coefficient for each node/electrode " C → i " in the directed graph is defined as the fraction of directed edges between adjacent nodes of node i over the maximum amount of directed edges C → i = 1 2 ∑ j , h ∈ V ( e i j + e j i ) ( e i h + e h i ) ( e j h + e h j ) ( k i out + k i i n ) ( k i out + k i i n − 1 ) − 2 ∑ j ∈ V e i j e j i (5) where " k i " is the degree of node i (Watts and Strogatz 1998; Fagiolo 2007; Liao et al., 2011). The clustering coefficient of the network is the mean clustering coefficient of all nodes. C = 1 n ∑ i ∈ V C → i (6)

The global efficiency of the network is a measure which quantifies how information flows throughout the network. Graphs with a high global efficiency have on average shorter paths connecting any two nodes within the network. The global efficiency was used instead of the characteristic path length as the shortest path length is not defined when a network contains two nodes that are not connected by any path. The global efficiency is the average efficiency, defined as the inverse of the shortest path between two nodes, over all electrode pairs (Latora and Marchiori 2001; Mitsis et al., 2020). E = 1 n ( n − 1 ) ∑ i , j ∈ N ,   i ≠ j 1 d i j (7) where " d i j " is the shortest path between nodes i and j in the directed graph. Graph theory measures were computed using the Brain Connectivity Toolbox in Python (Rubinov and Sporns 2010). Graph visualization was created using the circular layout graph of the MNE-Python software package (Gramfort et al., 2013).

DABEST Python toolbox was used for two group comparisons of graph measures generating Gardner-Altman estimation plots for independent group mean differences (Ho et al., 2019). Bootstrapping was used to obtain distribution and confidence intervals for difference in groups. The Wilcoxon rank sum test was also used to test for significance between the two groups, as a separate test from the bootstrapping confidence intervals.

FCN gives insight into disease induced changes in synaptic plasticity and efficiency of communication within neural networks in the brain (Bettus et al., 2008). WPC has previously been used as a measure for FCN in the brain, depicting coupling of different brain regions by way of specific brain rhythms (Cotic et al., 2015). Figure 1 shows scalp EEG traces of representative seizures in a non-SUDEP (P2) and SUDEP patient (P19), and examples of the corresponding WPC between two electrodes. The chosen electrodes showed the highest closeness centrality during seizure. Differences in the WPC distributions reaffirmed the choice of the two frequency ranges. The analysis was repeated for each pair of nodes and averaged over the frequency range and temporal windows to provide the connectivity strength between the two nodes.

Complex information flow involves multi-frequency large-scale organization in the brain (Buzsaki 2006; Jirsa and Müller 2013). We used directed networks to provide information about the coupling directionality between low and high frequency rhythms. The directed graph is constructed using the low frequency and high frequency with edge weights corresponding to the product of the low frequency graph degree of one electrode to the high frequency graph degree of another electrode. To validate our cross-frequency directed network, we simulated EEG rhythms and low/high frequency hubs. Starting with an empty set of EEG signals, we began populating specific channels with chosen delta and HFO rhythms. All channels contained Gaussian white noise. We aimed to create one low frequency hub and one high frequency hub and confirm that the directed connectivity network showed a directed edge between them. The low frequency hub was chosen to be electrode P7. This hub was designed to contain two different delta rhythms which then would each spread to another electrode (P3 and O1). The high frequency hub was chosen to be electrode F8. This hub was designed to contain two different HFO rhythms which then would each spread to a nearby electrode (Fp2 and C4). As expected, the low and high frequency networks shown in Figures 2C,D showed the desired connections. The directed connectivity network correctly identified the connection between the low and high frequency hubs as per our network design. Furthermore, the network analysis proved to be consistent when applied to recorded EEG data from a SUDEP (P19) and non-SUDEP patient (P1) (Figure 3). Taken together, these results validate and demonstrate how the directed graph captures cross-frequency interactions within the brain.

We constructed the delta-HFO directed network for consecutive 10-s duration time windows during peri-ictal regions in representative SUDEP (P19) and non-SUDEP (P1) patients (Figure 4). The graph visualizations indicated differences in network seizure dynamics in a non-SUDEP patient that were not observed in a SUDEP patient (Figure 4). The directed graphs (Figures 4D–G) further highlighted the seemingly unchanging network of the non-SUDEP patient during the ictus. Topological measurements of the directed graphs were used to compare the networks between the two patients. The pre-ictal, ictal, and post-ictal mean clustering coefficient and global efficiency were higher in the SUDEP than in the non-SUDEP patient (Figure 4D).

The clustering coefficient and global efficiency measures were used to compare group differences between non-SUDEP and SUDEP epileptic patients (Figure 5). The clustering coefficient is an average measure of how node triples are connected within the network and specifies the tendency for nodes to cluster together. The clustering coefficient was significantly higher in the SUDEP group during the pre-ictal, ictal, and post-ictal segments (pre-ictal: p = 0.00100, ictal: p = 0.00001, post-ictal: p = 0.00100). The global efficiency, which indicates how efficiently information is transferred between nodes, was shown to be significantly higher in the SUDEP group pre-ictally, ictally, and post-ictally (pre-ictal: p = 0.00012, ictal: p = 0.00001, post-ictal: p = 0.00071). These results show that network topology is a potential biomarker in assessing SUDEP risk.