Extracting the transition network of epileptic seizure onset

The onset of an epileptic seizure is often accompanied in the electroencephalogram (EEG) by a more or less characteristic transition of the so-called background EEG to abnormal seizure discharges.¹ In invasive recordings of focal-onset seizures from patients during presurgical monitoring, the onset of abnormal discharges is seen in a subset of recording channels.² The aim is to find out which anatomical areas are involved in the cooperative dynamics that represents the transition. However, the interpretation of interactions between channels may not be feasible from visual inspection. As such, new multivariate methods are needed to reduce the dimensionality of the problem and to quantitatively assess the variability of interactions between locations at seizure onset in a patient. Recently, this has been addressed by so-called network approaches, which combine the interpretation of uni- and bivariate measures.³ 

Electrographically, the seizure onset is viewed as the appearance of abnormal rhythmic (but not necessarily periodic) discharges from a background that is assumed to be pseudo-stationary for a period of a few seconds prior to seizure.⁴ Here, we adopt this view and assume a two-state model (the background and the ictal state) with a transition period between them. With that assumption, it is possible to consider the last few seconds of background and the first few seconds of seizure as two pseudo-stationary dynamics on two respective branches of a slow manifold in state space and two (comparatively) fast transitions from one branch to another. This idea was discussed in the context of deterministic flows in low-dimensional state space by Rössler in 1976.⁵ In that example, the deterministic chaos was explained as a transition from a pseudo-fixed point to an unstable oscillatory dynamics by recurrent transitions that were illustrated as fast switchings from one branch of the slow manifold to another [see Fig. 1(a)]. The transition to oscillations was illustrated as a straight line as if falling from a cliff. However, in the context of more complex bursting dynamics, the same kind of illustration reveals that there are many possibilities of how the transition can be organized.⁶ While it is straightforward to illustrate this in deterministic dynamical systems with three variables, the picture becomes more complex in bursting dynamical networks.^7,8 A major problem when trying to compare the computational model output to the clinical EEG of epileptic seizures is that the EEG is too high dimensional for a simple illustration.⁹ Nevertheless, knowledge about the transition dynamics (or at least the deterministic components of it) is crucial to better understand the nature of seizure onset (and equally, the seizure offset).^10–13 

Here, we propose a data-driven workflow to extract uni- and bivariate features of the seizure onset dynamics from the EEG in clinical focal-onset seizures. Changes in recurrence of state-space motifs are used to separate the interictal and the ictal state and the joint evolution of trajectories are illustrated in state space. The quantification allows the composition of the seizure-onset network and quantification of the consistency of contributions from anatomical locations to the onset. The method can be automated and can be applied to any state transition in the EEG.

We use intracranially recorded EEG (iEEG) of patients with recurrent focal-onset seizures. Data were from the Neurophysiology Department, Great Ormond Street Hospital, London, UK. All patients are pediatric and stereo EEG recordings for drug-resistant epilepsy with spontaneous, focal-onset seizures. The aim of the implantation was the identification of the seizure onset zone for subsequent epilepsy surgery. Sensors were invasive depth electrodes (Dixi Medical, Besançon, France) with contacts −1.5 mm apart. For data acquisition, Natus NeuroLink IP amplifiers (Natus, Oakville, Canada) with a sampling rate of 1024 Hz were used. The position of the electrodes in patient 1 is shown in Fig. 2(a).

Seizures included are those showing a focal onset (not generalized) and are of the same clinical type. The seizure onset is determined from the visual inspection of the complete recording by a trained neurophysiologist according to its definition as the “first unequivocal intracranial EEG sign of change from the background that led to a clear seizure discharge.”⁴ The number of channels and implantation varied between patients. Data are imported from EDF format using the PYEDFLIB library in Python. Before processing, data are filtered with a high pass to reduce the slow shifts on time scales longer than the selected segment and with a low pass set to half of the sampling rate (the Nyquist frequency) to exclude components above 512 Hz. Filtering is done with a bidirectional Butterworth filter of order 5. We selected 10 s segments of iEEG centered at the clinically determined seizure onset, see Fig. 2(b). Details of the recordings are summarized in Table I.

The workflow can use any univariate measure to determine the state transition. Here, we tested the frequency of maximal power, band power, and autocorrelation. The results below are obtained for the change in the frequency of maximal power obtained from the Fourier spectrum of the sliding window. It has the advantage that it is independent of band-specific contributions and picks up the onset in seizures with notable seizure rhythm following the transition period.

As with the univariate measure, it is possible to use any bivariate quantification in the workflow. Here, we searched and selected a measure that is not restricted to the similarity between pairs of time series (like correlation, coherence, or phase-locking value). Instead, we implemented a measure that quantifies the recurrence of motifs in two-dimensional state space. A sliding window (referred to as recurrence window) is defined to search for similar recurrences of motifs within each large window. Recurrence is determined as the distance between two pieces of trajectory in two-dimensional state space. Each segment is normalized to minimum and maximum to make this distance independent from changes in amplitude. We calculated the Euclidean distance between normalized pieces of trajectories of the same length. The distance is zero if the two trajectories are identical and are large for unrelated trajectories. Using an overlap (shift) of 20 data points, all possible distances are calculated, and the smallest distance is recorded for each large window, and the results are stored as a time series of similarities for all pairs of channels. Two exemplary recurrence matrices are shown in Fig. 3. They are taken from the beginning and the end of the iEEG segment of 10 s in Fig. 2, respectively, and should thus represent the best recurrences during the background (A) and the ictal state (B), respectively. The matrices are symmetric and thus only one triangular part is filled with similarity values. The pattern of similarities varies between the two large windows and shows, in this case, a decrease in distances (and thus increase of similarities) for certain windows. The subsequent analysis is done using only the lowest score (highest similarity) from each matrix.

After assessing the uni- or bivariate quantity for channels or channel pairs, we fit a sigmoid function to the resulting time series. The function is taken as $f=a1+exp(−b(x−c))+d$ where we set $d=0$ because the data are min–max normalized before fitting, and a, b, and c are free parameters adjusting the separation of states, the location of the inflection point, and the steepness at the inflection point, respectively. The fit returns the optimal parameters of the function for each channel, and the result is scored to find the best candidates for a state transition.

Scoring of the fit is calculated from three contributions: (i) the goodness of fit as expressed in the $R2$ value; (ii) the separation of the states (the more the better, maximum 1, minimum 0); and (iii) the steepness at the inflection point (the steeper the better, normalized to the maximal steepness in the dataset between 0 and 1). The total score is taken as the arithmetic mean of the three partial scores.

We include optimization steps to address the problem of robustness of results. We scan the sliding window length to optimize the sum of the N best scores for the N best channels to be displayed as network nodes. Similarly, we optimize the sum of the M best scores of the recurrence window length for all channel pairs to be displayed as network edges.

In the first round, the scoring is done for all iEEG channels using the univariate quantification to return the channel order according to the fit. From a network perspective, these form the weighted nodes of the iEEG transition network. To reduce the computational load, the N top-scoring channels are then selected for bivariate analysis. A bivariate quantity is calculated from the same signal segment to obtain its temporal change as a function of time. This is done for all $N(N−1)/2$ (bidirectional) channel pairs to obtain the weighted edges of the onset network. Choosing the M top-scoring pairs of channels, we finally obtain the network with N nodes and M edges, along with their achieved score and the individual onset time (inflection point of the fitted sigmoid). The results are calculated for $N=M=10$ for the seizures specified in Table I.

The ranked channels and channel pairs along with their scores are then used to derive an across-seizure scoring. To calculate the across-seizure score, we check which channel or channel pair is found among the 10 highest-scoring pairs of a seizure. Each appearance gets a full mark per seizure. For each candidate, we then add ranking marks, from 1 for rank 1 to 0 for rank 10. The result is then normalized between 0 and 1. All channels or channel pairs that never make it into the first 10 score 0, while a channel or channel pair that scores first in all seizures obtains 1. All other channels or channel pairs lie somewhere in between. As the scoring is based on ranks, it is limited to discrete values.

Figure 2(b) shows the time series of voltages of all recorded iEEG channels during a segment centered on the clinically determined seizure onset. There are transitions from irregular background to large amplitude rhythmic voltage changes as indicated by the darker areas. The visibly discernible onset is not homogeneous but occurs differentially in some groups of channels. The arrows point to groups of channels with early onset. As an example, based on the state separation using the frequency of maximal power in the sliding window of the Fourier spectrum, the following channels for returned for Sz01 of patient 1 with the highest scoring: Q15–Q14, Q14–Q13, J11–J10, I4–I3, F14–F13, E7–E6, P10–P9, J6–J5, Q12–Q11, E2–E1, i.e., channels from different sEEG depth electrodes (cf. Fig. 2). For the results of all seizures, see Table II. These 10 top-scoring channels were selected for bivariate analysis.

For each of the channel pairs (in this case 45), the recurrence matrix is calculated within each large sliding window. Two exemplary recurrence matrices are shown in Fig. 3. They are taken from the beginning and the end of the iEEG segment in Fig. 2(b), respectively, and thus represent an instance during the background (left) and the ictal state (right), respectively. The matrices are symmetric and thus only one triangular part is filled with similarity values. The pattern of similarities varies between the two large windows and shows, in this case, a decrease in distances and thus an overall increase of similarities. Note that the downstream analysis is done only with the lowest score (highest similarity) within the matrix.

Figure 4(a) shows the sliding window of the highest similarity for all channel pairs of the selected iEEG segment. Comparing the time series in Fig. 2(b), the onset of a region of higher similarity (dark blue, right-hand side) seems to occur around the annotated seizure onset. Some temporal structure can be seen within the seizure activity but is qualitatively similar for the channel pairs displayed. Note that the change in color reflects a transition in the dynamics in terms of recurrence in state space and is not a reflection of changes in signal amplitudes. The fact that the recurrence measure changes consistently implies that the recurrence found during the ictal state is not the same as the one before the clinically determined seizure onset.

The seizure onset timings (as determined by the inflection point of the fitted sigmoid) for the 10 pairs with the highest score are displayed in Fig. 4(b). The two-state coloring shows that the onset points are distributed within the interval and, in particular, three pairs show onset before the marked onset point. According to the coloring, all switches are from low to high recurrence similarity.

The bivariate scores for all channel pairs are displayed in Fig. 4(c). The scores are ordered from high to low and demonstrate the general finding that some pairs score significantly higher than the majority of pairs. We take this as an indicator that these pairs show the strongest change in recurrence feature and are thus optimal for the display of the dynamics of the state transition from background to ictal activity. Overall, the scores decrease continuously in this particular case and do not indicate any specific groups of pairs.

From this analysis, we thus obtain the best channel and the best channel pairs with their respective scoring and onset time based on the state transition model. For all seizures studied, we find that the best channels and best pairs show good separation (scoring above 50%), which is consistent with the clearly visible seizure onset.

Figure 5(a) shows the reconstructed network containing the nodes and edges, which display the strongest transition in the period of seizure onset. All nodes are included and the nodes that have already passed their onset point according to the analysis are colored in red. At this point in time, the two nodes with the highest degree are still gray, meaning they have not yet passed the onset point. All edges are included and in this case, all edges have already passed the onset point and are thus colored in blue. Figure 5(b) shows the same network but with node size varying according to the degree of the node (based on M edges). One node is found to not have any of the highest-scoring edges attached to it.

Figure 6 shows a state-space representation of the voltage for the two highest-scoring channel pairs as a time series (from top to bottom), similar to the right-hand side of the sketch in Fig. 1. The coloring for the individual trajectories shows the individual phases, black for background and magenta for ictal dynamics. There are relatively few turns of the trajectories because the transition is sharp. Also, the transition (as determined by this method) does not involve large-amplitude oscillations. Large amplitude oscillations appear when the trajectory reaches the bottom of the display, i.e., after the transition.

Having the unique result of onset properties for a feature, we exemplify the transition dynamics with two attempts to visualize the seizure onset dynamics. Figure 7(a) is a ribbon plot of the voltages of two pairs for the transition period in a single display. Note that one of the axes shares the same variable, but the other has two variables. This is equivalent to plotting two scatterplots in one. The ribbon then highlights the co-evolution of the two trajectories. Figure 7(b) is based on the data for the same channels but the time series were band pass filtered around the frequency for which the optimal state separation was obtained (i.e., the optimization of the recurrence window length).

Figure 8 shows displays of the (Euclidean) distances between the voltages of the three highest-scoring pairs that share at least one channel. Figure 8(a) is the time series of the three distances, which show the co-evolution of the distances in a waxing and waning pattern. Figure 8(b) displays the same distances in a two-state-time space as in Fig. 7 (time from top to bottom). The co-evolution indicates a strongly coordinated evolution (as if in a swarm) between the electric potentials of the three-channel pairs.

Doing the analysis for repeated seizures in the same patient, we calculated the across-seizure score. Results are summarized in Table II. The channels are the findings for each seizure after optimization. The results indicate that different patients may show different degrees of consistency of channel onset. Specifically, patient 2 shows the highest consistency, whereas patients 1 and 3 have lower onset channel consistency which agrees with a visual inspection.

A major unknown concerning the dynamics of epileptic seizures is the components of the EEG signal that reflect genuine contributions to the recurrent onset of seizures in patients with intractable epilepsy.¹⁴ Also, the dynamical nature of the transition is still unclear.¹⁵ It may involve random switching between states in a global state of bistability,¹⁶ change of dynamics following the passing of a control parameter through a bifurcation point,^17,18 or a case of the randomly induced complex transient in an excitable system.¹⁹ However, it is widely assumed that the transition follows some organizing nonrandom components, mostly due to the emergence of pronounced rhythmic activity in some EEG channels.²⁰ Computational modeling has typically settled on one or the other possibility without actually providing evidence for either against the others (e.g., in terms of making testable predictions). Here, we have elaborated a data-driven approach to extract and quantify seizure onset dynamics, which is independent of mechanistic assumptions. The main assumption made is that of the state transition between a nonepileptic and an epileptic state which is, however, the widely accepted two-state model of epilepsy with a transition period separating the two states.¹⁴ There is a confirmed correspondence between EEG abnormalities and the clinical symptoms in many cases even if their exact correspondence cannot always be determined quantitatively. As we have not included clinical symptoms, our analysis focuses exclusively on iEEG signatures of seizure onset. Our main point of reference is the visually determined onset time which is a single point in time for all iEEG channels. The concept of a preictal state is not considered in the present analysis as there does not seem to be a general agreement about generic markers to identify its presence in invasive recordings.

EEG analysis and its application to the characterization of the transition from background to ictal activity has a long history.^21–23 Many studies were done using univariate features (mostly Fourier components or related), but there is also a line of research into multivariate analysis.²⁴ As in previous approaches, we use a sliding window approach that implicitly assumes the existence of (pseudo-) stationary states. However, we exploit the knowledge of the existence of a neurologically important state transition (seizure onset) to optimize the analysis based on the feature of state separation. In an everchanging signal like the EEG, a definition of “state” can only be empirical and will depend on the time scale of interest. The scale chosen in our analysis is consistent with the clinical practice involving the onset of electrographic signs over periods of less than a minute.

The inspiration for looking at a state transition as a (comparatively) fast process happening between two branches of a slow manifold representing two pseudo-stationary states comes from the treatment of deterministic chaos in state space as introduced, e.g., in Ref. 4. There, the dynamical property of mixing of trajectories was illustrated by a (deterministic) flow that switches between two subtypes of dynamics (see Fig. 1 in Ref. 5). Rössler extended this idea to arbitrary dynamics happening on either branch of the slow manifold with autonomous switching between them (see Fig. 2 in Ref. 5). Applying this epileptic dynamics led to the proposal to treat epileptiform rhythms as a type of bursting^25–27 and the subsequent modeling of spontaneous transitions as a slowly modulated transition between a time-varying fixed point and complex bursting.⁷ The nature of the slow process has been a subject of speculation²⁸ but so far no conclusive evidence has been obtained for the case of human data (however, see Ref. 29 for some recent observations). Nevertheless, the recurrent and apparently spontaneous nature of (in particular) focal epileptic seizures is at least compatible with an underlying slow process guiding the transitions.^30–32 Here, we have used a two-state (sigmoidal) function as an underlying model of the slow process and then quantified to what degree this leads to a successful state separation within the period of interest. Conveniently, the same type of state separation analysis can be performed with both uni- and bivariate measures. The scoring of the state separation then offers a possibility of an informed ranking of the high-dimensional data (typically on the order of 100 recorded channels with around 5000 possible connections between them). As such, the workflow can be applied to any case where a state transition is seen in EEG data, e.g., the fragmentation at seizure offset,^33,34 the onset or offset of REM sleep,³⁵ and others.

Specifically, for the multivariate analysis, we have exploited a recurrence measure as a bivariate quantity that overcomes some of the shortcomings of visual analysis and thus offers additional insight. First, recurrence is not dependent on signal amplitude. It quantifies whether a motif recurs in state-space after normalization (similar to comparing two objects independent of their size). Second, the finding is independent of repetition (a certain dominant Fourier frequency) as it finds the best recurrence within a window independent of its respective distance in time. This can be of particular value when analyzing, e.g., an epileptic seizure that has large but irregular spiking at seizure onset.^36–38 The potential problem of finding an accidentally high recurrence similarity is partly accounted for by the fitting process which smoothens many individual values and leads to an increase in the level of the quantity only if there is a substantial change in recurrences over a larger period of time.

The analysis offers the possibility to detect the consistency of seizure-onset dynamics between any number of seizures in the data from a given patient. As the duration of the invasive recording is limited by clinical considerations, no systematic studies of the onset dynamics are possible. Nevertheless, typically a few seizures are being recorded and the analysis with the proposed across-seizure measure offers the calculation of an estimate of the consistency of the involvement of channels or channel pairs (and thus anatomical locations and connections between them) from the available number of seizures. Note that the suggested optimization methods help to increase the robustness of the qualitative results. Nevertheless, different data will lead to different results, e.g., applying the analysis to individual frequency bands will yield different results. However, this is not a shortcoming of the method but a reflection of the actual complexity of the distributed epileptic dynamics. One possibility to account for different frequency content and filter settings is to compare the results against an ensemble of independent stationary Fourier surrogates to assure that only findings are considered which score above chance level.

Finally, from the point of view of network theory, it has been argued that seizure onset is not necessarily the result of a strictly localized abnormality but that even a focally appearing onset may be the result of a network interaction happening at different (nearby or distant) locations.^39,40 We have, therefore, combined a univariate and bivariate analysis to automatically extract the most prominent network nodes and the strongest connections between them based on state separation scoring. The evolution of the resulting “ictogenic” network (cf. Figure 5), as well as the “swarming” of distances between potentials as in Fig. 8, then offers new ways to animate the complex dynamics accompanying focal seizure onset.

This research was supported by the National Key R&D Program of China (No. 2020YFC2004400), National Natural Science Foundation of China (NNSFC) (Grant Nos. 11932003, 11772019, 11832003, 11772016, and 11702008), and Key Project of Science and Technology of Beijing Municipal Education Commission (No. KZ201810005007). G.B. and L.Z. would like to thank Yujiang Wang and Peter Taylor for the discussion. G.B. thanks Otto Rössler for continued inspiration since 1982. All EEG recordings were performed with informed consent from the patients' legal guardians, and as necessitated by the patients’ clinical course. The use of anonymized EEGs from the clinical database for quantitative analysis was reviewed and approved by the UCL Great Ormond Street Institute of Child Health Joint Research and Development Office.

The data that support the findings of this study are available from the corresponding author upon reasonable request.