Detecting seizure onset and offset times using human intelligence: A critical-transitions-based approach Andrew Flynn1,2,* , Cian McCafferty3 , Klaus Lehnertz4,5,6 , François David7 , Vincenzo Crunelli8 , Gordon Lightbody2,9 , and Sebastian Wieczorek1 1 School of Mathematical Sciences, University College Cork, T12 XF62 Cork, Ireland. 2 INFANT Research Centre, University College Cork, T12 DC4A Cork, Ireland.
3 Department of Anatomy & Neuroscience, University College Cork, Cork, Ireland.
4 Department of Epileptology, University of Bonn Medical Centre, 53127 Bonn, Germany.
5 Helmholtz-Institute for Radiation and Nuclear Physics, University of Bonn, 53115 Bonn, Germany. 6 Interdisciplinary Center for Complex Systems, University of Bonn, 53175 Bonn, Germany.
7 Center for Interdisciplinary Research in Biology, Collège de France, 75005 Paris, France
arXiv:2607.27105v1 [math.DS] 29 Jul 2026
8 Department of Pharmacology and Neuroscience, Faculty of Medicine, University of Lisbon, Lisbon, Portugal 9 Department of Electrical and Electronic Engineering, University College Cork, T12 YF78 Cork, Ireland.
ABSTRACT Most existing seizure detection algorithms require extensive pre-processing of the data and rely on heuristic or currently unexplainable machine learning approaches. These approaches often struggle with balancing detection sensitivity and specificity in the presence of variable seizure morphologies, interictal epileptiform discharges, and artefacts. Here, we consider an alternative approach: our seizure detection algorithm, which is based on the concept of critical transitions and overcomes the aforementioned limitations. Specifically, we perform a receiver-operating-characteristic analysis to quantify the performance of our algorithm in terms of its agreement with expert annotations of seizure onset and offset times in the voltage recordings of seizure activity in epileptic rodents with different seizure morphologies. We demonstrate how performance depends on algorithm parameters and varies across different rodent recording sessions. We determine the optimal set of algorithm parameters for each recording session, with near expert-level performance achieved in most cases. Finally, we derive a single general set of algorithm parameters applicable across all recording sessions. The algorithm maintains its high performance in this general setting, demonstrating its versatility, robustness across varying seizure morphologies, and potential to complement machine learning algorithms.
Introduction Seizure detection algorithms have become a valuable tool in the diagnosis, monitoring, and management of epilepsy, a condition which effects 50 million people worldwide1 . However, despite decades of research and improvements in obtaining voltage recordings of seizure activity from the brain, significant challenges remain in developing robust and accurate seizure detection algorithms. Factors that add to this challenge include the susceptibility of voltage recordings to artefacts2 , the presence of interictal epileptiform discharges3, 4 , and the variability of seizure morphology across individuals, recording sites, and recording sessions. As a result, many algorithms struggle to balance sensitivity and specificity - correctly identifying seizure and nonseizure intervals - in clinically realistic settings, yet major progress has been made thanks to modern machine learning methods5, 6 . However, the price to pay for using certain machine learning techniques is generally their black-box nature and tendency to confabulate (hallucinate); they lack explainability and generate false, but sometimes plausible, predictions/information7 . Addressing the above challenges is essential for advancing personalised epilepsy care and realising the full potential of neurotechnological interventions, such as cortical stimulation8, 9 , which require real-time information on the state of the brain. In this paper, we conduct a more detailed analysis on the performance of a seizure detection algorithm introduced in our previous work, Flynn et al.10 . To provide a different viewpoint to and complement recent trends in using artificial intelligence to develop seizure detection algorithms, we use an ancient technology known as ‘human intelligence’. To develop this algorithm we first acknowledged that the epileptic seizures we analysed are states of high-amplitude, synchronous electrical activity in the brain with defined onset and offset11 . Further, while seizures have a wide range of associated symptoms, seizures constitute a sudden and large change in brain state and have a commensurate impact on the life of a person with epilepsy12 . From a mathematical perspective, a sudden and large change in the state of a complex system corresponds with Ashwin et al.’s13 definition of a critical transition (CT). The mathematical theory of CTs has been instrumental in describing and foreseeing sudden and large changes in climate and ecological systems14, 15 and more recently in shaping global policymaking16, 17 . Inspired by its impact on the environmental sciences, we consider the brain as a complex system and apply the theory of CTs to detect seizure activity in the brain. Specifically, we model seizure onset and offset as CTs between a non-seizure state (NS state) and seizure state (S state), and develop an algorithm to detect CTs between these states in voltage recordings of brain activity. Importantly, our algorithm is based entirely on mathematical definitions of these states and CTs between them. These definitions are informed from observations of seizure activity in local field potential (LFP) recordings but generalisable to electroencephalography (EEG) recordings. Our algorithm detects CTs through monitoring only a single feature of the data, the voltage recording at a given time, i.e., we simply use the data we are given as input to our algorithm. In our previous work, Flynn et al.10 , we introduced our algorithm as part of a framework for classifying seizure onset in terms of different types of CT. We applied this framework to voltage recordings of seizure activity in Genetic Absence Epilepsy Rats from Strasbourg (GAERS), specifically measurements of the LFP in mV. In the present paper, we conduct a more detailed analysis on the performance of our algorithm when applied to the same set of voltage recordings; see part M1 of the Methods section for details on how these voltage recordings are obtained and annotated by an expert. These voltage recordings provide a reasonable challenge for any seizure detection algorithm given there are numerous artefacts and interictal epileptiform discharges present and seizure morphology varies significantly between GAERS. In particular, the change in voltage amplitude from the NS to S state varies between GAERS, with some exhibiting relatively significantly larger changes than others. In practice, a well-designed seizure detection algorithm will identify seizure and non-seizure time intervals that agree with expert annotations. However, one must carefully choose what performance metrics are used to quantify this agreement. In this paper we take inspiration from Temko et al.5 and Mathieson et al.6 by using performance metrics based on receiveroperating-characteristic (ROC) curves. These metrics enable us to determine the optimal choice of algorithm parameters for each recording session we analyse, with very good performance achieved on average. Based on this we derive a general set of algorithm parameters applicable across all recording sessions. We show that the algorithm maintains much of its accuracy in this general setting, demonstrating its robustness across varying seizure morphologies. Overall, our results demonstrate the benefits of studying the epileptic brain from a CT perspective. We produce a new approach for accurate and automated seizure detection that is grounded in mathematics and can be adapted depending on definitions of NS and S states. Since our approach does not rely on heuristic or currently unexplainable machine learning approaches, it could complement these approaches, improving their overall performance beyond what can be achieved by adjusting the parameters of the machine learning algorithm.
Results Detecting seizure activity in voltage recordings In this subsection we briefly describe the main premise and technical aspects of our algorithm which are needed to discuss our results. 2/23
Main premise of our algorithm
We designed our algorithm through observation of generic characteristics and expert annotations of seizure activity. We use Fig. 8 from the Methods section as an example where, according to an expert, seizure onset happens at 𝑡 = 𝜏1 ≈ 8.5s and offset at 𝑡 = 𝜏2 ≈ 22s. From a CT perspective, we observe: (O1) a non-seizure state (NS), characterised by a small-amplitude and weakly-correlated type of fluctuation around 0 mV, (O2) a seizure state (S), characterised by a large-amplitude and strongly-correlated type of fluctuation around 0 mV, and (O3) two CTs between the NS and S states, characterised by a sudden change from one type of fluctuation to another. We consider these two CTs in Fig. 8 to be analogous to seizure onset and offset, respectively. This approach is motivated by similarities between the International League Against Epilepsy’s definition of a seizure11 and Ashwin et al.’s definition of a CT13 . Main technical aspects of our algorithm
Inspired by a ‘non-ideal relay’18 , we define CTs in terms of successive crossings of two voltage thresholds and use a ‘moving window analysis’ to assess how long the brain remains in the new state after successive crossings of the two voltage thresholds. By using two thresholds we avoid the shortcomings of a single threshold approach in the case of CTs to the S state; a single threshold can be crossed multiple times within a short time interval given the high frequency of voltage measurements or, naturally, when in the S state given its large-amplitude oscillatory-like nature, both of which could result in false detections. Briefly, our algorithm detects a CT from the NS to S state at time 𝑡 = 𝑡1 if the absolute value of the time series exceeds an upper voltage threshold 𝛼 and then continues to exceed a lower voltage threshold 𝛽 frequently enough for a period of at least 𝜏S . Similarly, our algorithm detects a CT from the S to NS state at time 𝑡 = 𝑡2 if the absolute value of the time series falls below 𝛽 and then does not exceed 𝛼 for a period of at least 𝜏NS . A window of length 𝜏𝑤 is moved along the time series in discrete steps of size Δ to check whether the thresholds are exceeded in each window throughout the specified period. Part M2 of the Methods section provides precise definitions of when the algorithm detects a CT, specifies the algorithm’s technical details, and how the algorithm is designed to mitigate the influence of artefacts and interictal epileptiform discharges. In Fig. 1 we use the actual voltage recording from Fig. 8 to illustrate in three simple steps how our algorithm detects a CT in voltage recordings. The remainder of the present paper is focused on tuning the parameters of our CT detection algorithm to maximise its agreement with the expert’s annotations in terms of the corresponding seizure and non-seizure time intervals. We denote seizure (𝐸) and non-seizure time intervals according to the expert by 𝐼S(𝐸) and 𝐼NS . The ends of these intervals are defined by the values of (𝐴) seizure onset times 𝜏1 and seizure offset times 𝜏2 . Similarly, we use 𝐼S(𝐴) and 𝐼NS to denote the seizure and non-seizure time intervals obtained from the algorithm, by detecting the times of different CTs, namely 𝑡1 and 𝑡2 ; see part M3 of the Methods section for precise definitions. In a given recording session, we denote the number of 𝐼S(𝐸) intervals by 𝑁S(𝐸) , and the number of (𝐸) (𝐸) (𝐸) (𝐴) 𝐼NS intervals by 𝑁NS , where 𝑁NS = 𝑁S(𝐸) − 1. We use a similar notation for the algorithm, namely 𝑁S(𝐴) and 𝑁NS .
Classification of time intervals, constructing ROC curves, and algorithm performance metrics In this paper we take inspiration from Temko et al.5 and Mathieson et al.6 and evaluate the performance of our algorithm using metrics based on receiver-operating-characteristic (ROC) curves, a long standing method from signal detection theory to evaluate sensitivity and specificity, introduced in the 1940s19 . The points on our ROC curves are derived from classifications of (𝐴) (𝐸) 𝐼S(𝐴) and 𝐼NS where the 𝐼S(𝐸) and 𝐼NS are taken as the ground truth. In this subsection, we briefly describe the time interval classification procedure, how we construct ROC curves using this procedure, and the performance metrics we obtain from ROC curves. Time interval classification procedure
In short, a given 𝐼S(𝐴) is classified as true positive (TP) if it overlaps with at least one 𝐼S(𝐸) . Otherwise 𝐼S(𝐴) is classified as false
(𝐴) (𝐸) (𝐴) positive (FP). Similarly, a given 𝐼NS is classified as true negative (TN) if it overlaps with at most one 𝐼NS . Otherwise 𝐼NS is classified as false negative (FN). We introduce two parameters, 𝜌S , 𝜌NS ∈ [0, 1], to quantify how much overlap is required to (𝐴) (𝐴) classify a given 𝐼S(𝐴) as TP and 𝐼NS as TN. When 𝜌S = 𝜌NS = 0 any overlap between intervals of 𝐼S(𝐴) and 𝐼S(𝐸) , and 𝐼NS and
(𝐸) (𝐴) 𝐼NS , will suffice to classify 𝐼S(𝐴) as a TP or 𝐼NS as a TN, whereas when 𝜌S = 𝜌NS = 1 the intervals need to completely overlap, (𝐴) (𝐸) i.e., an 𝐼S(𝐴) needs to contain an 𝐼S(𝐸) or vice-versa, and similarly for 𝐼NS and 𝐼NS . See part M3 of the Methods section for a more precise description of our classification procedure. Figure 2 illustrates how our classification procedure works. Specifically, we pick a portion of the voltage recordings from recording session T2M (shown in Fig. 2 (a)) that was annotated by the expert. These annotations are shown in Fig. 2 (b) via a
3/23
(d)
β crossed at t = t2. Potential CT at t = t2. Moving window already active.
LFP [mV]
0.1 (a) α exceeded at t = t1. Potential CT at t = t1. Moving window activated. 0.0
t1
τw
t 1 + τS
τw
t2
t2 + τNS
-0.1 0.1 (b) Window shifted by m∆.
(e)
LFP [mV]
Window shifted by m∆.
0.0
τw
t1 + m∆
t2 + m∆
τw
-0.1 (f)
CT from the S to NS state detected at t = t2. Moving window deactivated.
LFP [mV]
0.1 (c) CT from the NS to S state detected at t = t1. 0.0
-0.1 7
8
t1 + nS ∆
τw
9
10
t2 + nNS∆ t [s]
21
22
23
τw 24
t [s]
Figure 1. Illustrating how the CT detection algorithm works, using the voltage recordings from Fig. 8 as an example. (a)-(c) Detecting a CT from the NS to S state at time 𝑡 = 𝑡1 . (d)-(f) Detecting a CT from the S to NS state at time 𝑡 = 𝑡2 . The algorithm parameters are chosen as 𝛼 = 0.05, 𝛽 = 0.04, 𝜏S = 2, 𝜏NS = 3, 𝜏𝑤 = 1, and we use 𝑚 = 300 in (b) and (e). The red and green horizontal lines indicate the thresholds of 𝛼 and 𝛽, respectively. (𝐸) binary sequence representation of the 𝐼S(𝐸) and 𝐼NS with coloured boxes beneath indicating these intervals as the ground truth. We use our algorithm to detect CTs between NS and S states in Fig. 2 (a) with algorithm parameters chosen as follows: 𝛼 = 0.08, (𝐴) 𝛽 = 0.07, 𝜏NS = 3, 𝜏S = 2, and 𝜏𝑤 = 1. We show the resulting classification of 𝐼S(𝐴) and 𝐼NS via the coloured boxes beneath the binary sequence representation of these time intervals for 𝜌S = 𝜌NS = 0 in Fig. 2 (c) and 𝜌S = 𝜌NS = 1 in Fig. 2 (d). By comparing Fig. 2 (c) to (d) we see the main effect that 𝜌S and 𝜌NS have on the classifications - more time intervals are classified as TPs and TNs when 𝜌S = 𝜌NS = 0 as opposed to 𝜌S = 𝜌NS = 1. In later sections we make specific reference to the performance of our algorithm when 𝜌S = 𝜌NS = 0.75 as Mathieson et al.6 set their equivalent parameters to this value in accordance with standards set by clinicians.
Constructing ROC curves
The important quantities we obtain from our classification procedure are the total number of TPs, FPs, TNs, and FNs in a given recording session which we denote by 𝑁TP , 𝑁FP , 𝑁TN , and 𝑁FN , respectively. These are used to define the following quantities which are needed to construct ROC curves: True positive rate (TPR): 𝑁TP ∕(𝑁TP + 𝑁FN ) ∈ [0, 1], also referred to as the ‘sensitivity’, quantifies the agreement between the expert and the algorithm on detecting the presence of seizure time intervals. False positive rate (FPR): 𝑁FP ∕(𝑁FP + 𝑁TN ) ∈ [0, 1], also referred to as ‘fall-out’ or ‘1 - specificity’, quantifies the disagreement between the expert and the algorithm on detecting the presence of non-seizure time intervals.
4/23
LFP [mV]
0.1 (a) 0.0
-0.1 S (b) Expert
NS S (c) Algorithm ρS = ρNS = 0
NS S (d) Algorithm ρS = ρNS = 1
NS 1700
1750
1800 TP
1850 FP
1900 TN
1950 FN
2000
t [s]
Figure 2. Illustrating how the seizure and non-seizure time interval classification procedure specified in part M3 of the Methods section works and the effect that 𝜌S and 𝜌NS have on the classifications. (a) shows a portion of voltage recordings that were annotated by the expert. (b) shows a binary sequence representation of the non-seizure and seizure time intervals according to the expert. (c) and (d) show the corresponding sequence of the non-seizure and seizure time intervals according to the algorithm. The boxes beneath these sequences are coloured green for TN, red for TP, orange for FN, and black for FP. The algorithm’s parameters were chosen as 𝜏NS = 3, 𝜏S = 2, 𝜏𝑤 = 1, with 𝛼 = 0.08 and 𝛽 = 0.07 (indicated by the horizontal red and green lines in (a)), 𝜌S = 𝜌NS = 0 in (c), 𝜌S = 𝜌NS = 1 in (d).
Based on the above, the closer the TPR is to 1 and FPR is to 0 the better the performance of our algorithm. It is convenient to show this in a two-dimensional plane with FPR values assigned to the x-axis and TPR values to the y-axis. In our case, we construct ROC curves by connecting (FPR, TPR) points between (0, 0) and (1, 1) for increasing values of the FPR where each point corresponds to different choices of the upper threshold 𝛼 while the other algorithm parameters remain fixed. Thus, the closer a given point on the ROC curve is to (0,1), the better the algorithm performs. We refer to the (0,1) point as the ‘point of ideal performance’. Why vary 𝛼 to construct ROC curves: We vary the upper threshold 𝛼 because it is the parameter the algorithm is most sensitive to. Across different GAERS there is a much greater difference in terms of the amplitude of fluctuation than the temporal characteristics of the voltage recordings; seen in part M1 of the Methods section by comparing the significant differences in voltage amplitudes across GAERS in Fig. 9 to the marginal differences in the probability densities of seizure and non-seizure time interval durations across GAERS in Fig. 10. Removing misleading points from ROC curves: We set TPR = 0 if 𝑁TP + 𝑁FN = 0 and set FPR = 0 if 𝑁FP + 𝑁TN = 0 in order to avoid dividing by 0. This typically occurs when 𝛼 is not large enough. See part M4 of the Methods sections for details on additional steps taken to remove misleading points from ROC curves. Remark: The FPR is related to the ‘true negative rate’ (TNR) via FPR = 1 − TNR, where TNR = 𝑁TN ∕(𝑁TN + 𝑁FP ). Temko et al.5 and Mathieson et al.6 construct their ROC curves using the TNR as opposed to the FPR. Our preference is to use the 5/23
FPR, and the difference from using the TNR is a different orientation of the ROC curve. Algorithm performance metrics
We use the following two metrics to quantify the algorithm performance’s in terms of the values of the upper threshold 𝛼 that are used to construct ROC curves according to the method specified above. Distance from ideal (DFI): We quantify the algorithm’s performance at individual 𝛼 values by computing the distance between the corresponding point on the ROC curve and the point of ideal performance, i.e., the (0,1) point. We refer to this distance as the √ ‘distance from the ideal for a given 𝛼’ and denote it by DFI(𝛼) ∈√[0, 2]. When DFI(𝛼) = 0 this corresponds to ideal performance, meaning all intervals are correctly classified, and DFI(𝛼) = 2 corresponds to the worst possible performance, meaning all intervals are incorrectly classified. We say the algorithm achieves excellent performance when 0 < DFI(𝛼) ≤ 0.1, very good performance when 0.1 < DFI(𝛼) ≤ 0.2, good performance when 0.2 < DFI(𝛼) ≤ 0.3, fair performance when 0.3 < DFI(𝛼) ≤ 0.4, and poor performance when DFI(𝛼) > 0.4. We denote the 𝛼 that minimises DFI(𝛼) for a given ROC curve as 𝛼 ∗ . Area under the ROC curve (AUC): We quantify the algorithm’s performance across the range of 𝛼 values used to construct an ROC curve by computing the ‘area under the ROC curve’. We denote this quantity by AUC ∈ [0, 1]. When AUC = 1, this means the algorithm achieved ideal performance for at least one value of 𝛼 used to construct the ROC curve. We say the algorithm achieves excellent performance when 0.9 ≤ AUC < 1, very good performance when 0.8 ≤ AUC < 0.9, good performance when 0.7 ≤ AUC < 0.8, fair performance when 0.6 ≤ AUC < 0.7, and poor performance when DFI(𝛼) < 0.6. An AUC < 0.5 means the algorithm performed worse than flipping a coin each time a seizure time interval is detected as to whether it will be classified as TP or FP. Remark: We compute both metrics to provide a more comprehensive assessment of the algorithm’s performance. Further, relying solely on the widely used AUC metric may result in a misleading evaluation of the algorithm’s performance. For instance, when comparing two ROC curves, one may have a larger AUC but the other may have a smaller DFI(𝛼 ∗ ), meaning the former shows better performance across a range of 𝛼 values but the latter shows better performance at a specific 𝛼 value. In our experiments, we compare the DFI(𝛼 ∗ ) obtained from ROC curves that were constructed for the same time series of voltage recordings and range of 𝛼 values but using different values for the other algorithm parameters. We refer to the set of algorithm parameters corresponding to the smallest DFI(𝛼 ∗ ) as the ‘optimal set of parameters’ for detecting CTs between NS and S states in a given time series of voltage recordings. Evaluating the algorithm’s performance in a single recording session In this subsection we use the DFI(𝛼) and AUC metrics to evaluate the algorithm’s performance on a single time series of voltage recordings containing multiple seizure and non-seizure time intervals. Through a series of experiments, we highlight the insights these metrics provide and use them to identify the optimal set of algorithm parameters for this time series. Technical details of experiments
We conduct our experiments using recording session T2M (recordings from rat T during the M𝑡ℎ recording session on day 2 of recording). This time series consists of ≈ 7, 500s (≈ 2 hours) of continuous voltage recordings taken in steps of 0.001s, (𝐸) contains 154 seizure intervals according to the expert, i.e., 𝑁S(𝐸) = 154, and the length of the largest 𝐼S(𝐸) was ≈ 56s and 𝐼NS was ≈ 1200s. (𝐴) We apply our algorithm to the entire T2M time series and classify the resulting 𝐼S(𝐴) and 𝐼NS in each parameter setting specified in Table 3. More specifically, we consider three different settings of 𝜏NS , 𝜏S , and 𝜏𝑤 and refer to these parameter (𝐸) settings as 𝑃1 , 𝑃2 , and 𝑃3 . 𝑃1 is chosen to reflect the expert’s annotation criteria as the length of the shortest 𝐼NS was ≈ 0.7s
and shortest 𝐼S(𝐸) was ≈ 1s. 𝑃2 and 𝑃3 are chosen to be slightly larger to provide a comparison. Furthermore, for 𝑃1 , 𝑃2 , and 𝑃3 , we consider a fixed separation between 𝛼 and 𝛽, denoted by 𝛿𝛼𝛽 ≥ 0, such that 𝛽 = 𝛼 − 𝛿𝛼𝛽 and 𝛿𝛼𝛽 < 𝛽 ≤ 𝛼. We consider values of 𝛿𝛼𝛽 = [0.01, 0.06] and 𝛼 ∈ [𝛿𝛼𝛽 + 0.01, 0.1]. We also consider values of 𝜌S , 𝜌NS ∈ [0, 1]. (𝐴) After using our algorithm to detect CTs in the time series and classifying the resulting 𝐼S(𝐴) and 𝐼NS , we construct ROC curves based on TPR and FPR values corresponding to 𝛼 ∈ [𝛿𝛼𝛽 + 0.01, 0.1] in each of the parameter settings specified in Table 3.
Insight from ROC curves - 𝜌S = 𝜌NS = 0, 𝛿𝛼𝛽 = 0.01
In Fig. 3 we show the ROC curves obtained for 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c). In each case 𝜌S = 𝜌NS = 0 and 𝛿𝛼𝛽 = 0.01, meaning 𝛼 ∈ [0.02, 0.1]. The red data points correspond to (FPR, TPR) values obtained for corresponding values of 𝛼. The green data points correspond to the 𝛼 ∗ obtained in each case, whose value is specified in the lower right corner. The AUC is also specified in the lower right corner. Figure 3 shows the algorithm performs best in parameter setting 𝑃3 , achieving near ideal performance at several values of 𝛼. More specifically, while the algorithm performs reasonably well in parameter setting 𝑃1 with AUC ≈ 0.86, there is a significant 6/23
1.0 (a)
(b)
(c)
TPR
0.8 0.6 0.4 0.2 0.0 0.0
: α = 0.07 AUC = 0.8578 0.2
0.4 0.6 FPR
0.8
1.0 0.0
: α = 0.085 AUC = 0.9528 0.2
0.4 0.6 FPR
0.8
1.0 0.0
: α = 0.08 AUC = 0.9745 0.2
0.4 0.6 FPR
0.8
1.0
Figure 3. ROC curves obtained when using parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c) (see Table 3). In each panel, (FPR, TPR) points plotted in red correspond to different values of 𝛼, the green point corresponds to 𝛼 ∗ (the optimal 𝛼), the lower right corner specifies 𝛼 ∗ and the AUC. improvement in the algorithm’s performance for 𝑃2 with AUC ≈ 0.95 and further for 𝑃3 with AUC ≈ 0.975. The DFI(𝛼 ∗ ) for each parameter setting provides similar insight with DFI(𝛼 ∗ ) ≈ 0.29 for 𝑃1 , DFI(𝛼 ∗ ) ≈ 0.15 for 𝑃2 , and DFI(𝛼 ∗ ) ≈ 0.09 for 𝑃3 . We found that increasing the values of 𝜏NS and 𝜏S beyond 𝑃3 resulted in little to no improvement in performance but significantly less seizure time intervals were detected by the algorithm in comparison to the expert; see Fig. S-1 in [20] for further details. Note, the constraint to construct ROC curves for increasing FPR values ensures the AUC can be computed. However, this means that two successive points on a given ROC curve may not necessarily correspond to a successive change in 𝛼; see Figs. S-2-S-4 in [20] for further details. In Figs. S-5 and S-6 in [20] we present results from studies on the level of agreement achieved between the algorithm and the expert which are not accounted for by ROC curves. Specifically, the agreement between values of 𝑡1 and 𝜏1 . Further, Temko et al.5 and Mathieson et al.6 use additional metrics derived from their classification procedure to construct curves such as, ‘seizure detection rate vs. false detections per hour’. We construct similar curves and find they provide complementary insight to ROC curves; see Fig. S-7 in [20]. AUC for different choices of 𝜌S and 𝜌NS
In Fig. 4 we extend the results shown in Fig. 3 by plotting heatmaps of the AUC values obtained when using different choices of 𝜌S and 𝜌NS ranging from 0 to 1. More specifically, the top row of Fig. 4 shows how the AUC varies for different choices of 𝜌S and 𝜌NS when using parameter setting 𝑃1 in (a), 𝑃2 in (b), 𝑃3 in (c), where 𝛿𝛼𝛽 = 0.01 in each case. The same is shown in the lower row of Fig. 4 for 𝛿𝛼𝛽 = 0.06. The black box/point in each panel of Fig. 4 corresponds to the largest AUC value (up to a difference of 0.001%) obtained for a given choice of 𝜌S and 𝜌NS . There is a box as opposed to a point in some cases as a small number of seizure and non-seizure time intervals are incorrectly classified across a range of 𝜌S and 𝜌NS values. Similar to Fig. 3, Fig. 4 shows performance improves as we move the parameter settings away from 𝑃1 to 𝑃2 and 𝑃3 . Additionally, Fig. 4 shows performance worsens as 𝛿𝛼𝛽 increases. For instance, in Fig. 4 (c) (𝑃3 with 𝛿𝛼𝛽 = 0.01) the AUC is greater than 0.96 for a wider range of 𝜌S and 𝜌NS values than in Fig. 4 (f) (𝑃3 with 𝛿𝛼𝛽 = 0.06). This effect is more pronounced for 𝑃1 and 𝑃2 . Furthermore, in terms of the more clinically relevant choice of 𝜌S and 𝜌NS used by Mathieson et al.6 , specifically, 𝜌S = 𝜌NS = 0.75, the algorithm maintains excellent performance when using 𝑃3 , achieving an AUC ≈ 0.94. We deduce from Figs. 3 and 4 that the optimal parameter setting of the algorithm for this time series is parameter setting 𝑃3 with 𝛿𝛼𝛽 = 0.01 and 𝛼 = 0.08. Evaluating the algorithm’s performance across multiple different recording sessions In this subsection we use the AUC and DFI(𝛼 ∗ ) metrics to evaluate the performance of our algorithm across multiple recording sessions from different GAERS. Informed by the results of the previous subsection, we restrict our study to parameter setting 𝑃3 with 𝛿𝛼𝛽 = 0.01. Furthermore, we only consider recording sessions where 𝑁S(𝐸) > 60 in order to base our claims on a reasonably sized sample; see the left column of Table 1 for details. The bar chart in Fig. 5 shows the resulting AUC for each recording session when choosing 𝜌S = 𝜌NS = 0 (in black) and 𝜌S = 𝜌NS = 0.75 (in red). In general, Fig. 5 shows that performance varies depending on which recording session the algorithm 7/23
1
(a)
(b)
(c)
0.8 0.6 ρNS 0.4
0.9745
0.2
0.9528
0.8578
0 1
(d)
(e)
(f)
0.8 0.6 ρNS 0.4
0.9734
0.2 0
0.9454
0.8354 0
0.2
0.4 ρ 0.6 S
0.8
10
0.2
0.4 ρ 0.6 S
0.8
10
0.2
0.4 ρ 0.6 S
0.8
1
AUC 0.5
0.6
0.7
0.8
0.9
1
Figure 4. Heatmaps which show the AUC of ROC curves for a given choice of 𝜌S and 𝜌NS . Each column corresponds to parameter settings 𝑃1 , 𝑃2 , and 𝑃3 as specified in Table 3. Each row corresponds to a different choice of 𝛿𝛼𝛽 where 𝛿𝛼𝛽 = 0.01 in (a)-(c) and 𝛿𝛼𝛽 = 0.06 in (d)-(f). The maximum AUC value is highlighted in the lower left corner in (a)-(f). is applied to, and in some cases there can be a significant decrease in AUC when 𝜌S = 𝜌NS = 0.75 as opposed to 0. More specifically, when 𝜌S = 𝜌NS = 0 the algorithm performs reasonably well, achieving an average AUC ≈ 0.89 across the different recording sessions with the largest (best) AUC ≈ 0.97 for recording session T2M and smallest (worst) AUC ≈ 0.75 for recording session K2E. The algorithm still performs reasonably well when 𝜌S = 𝜌NS = 0.75, however, the average AUC decreases to ≈ 0.785 with the largest (best) AUC ≈ 0.95 for recording session T8C and smallest (worst) AUC ≈ 0.51 for recording session S4E. The bar chart in Fig. 6 shows the resulting DFI(𝛼 ∗ ) for each recording session when choosing 𝜌S = 𝜌NS = 0 (in green) and 𝜌S = 𝜌NS = 0.75 (in purple). Similar to Fig. 5, Fig. 6 shows that the DFI(𝛼 ∗ ) also varies depending on which recording session the algorithm is applied to. In contrast to the AUC results, there is generally a relatively smaller change in DFI(𝛼 ∗ ) when changing 𝜌S and 𝜌NS from 0 to 0.75. More specifically, the average DFI(𝛼 ∗ ) increases from ≈ 0.174 to ≈ 0.194, the smallest (best) DFI(𝛼 ∗ ) increases from ≈ 0.075 to ≈ 0.076 for recording session T8C (best in both cases), and the largest (worst) DFI(𝛼 ∗ ) increases from ≈ 0.32 to ≈ 0.367 for recording session K5M (worst in both cases). Figures 5 and 6 show the AUC and DFI(𝛼 ∗ ) metrics are generally in agreement, i.e., when AUC is relatively large DFI(𝛼 ∗ ) is relatively small. However, there are some exceptions. For instance, for recording sessions K2E and K2M, the DFI(𝛼 ∗ ) metric indicates the algorithm achieves very good performance when choosing 𝜌S = 𝜌NS = 0 or 0.75, but the AUC differs significantly for these recording sessions; the AUC for K2M indicates the algorithm achieves excellent performance, but the AUC for K2E indicates the algorithm achieves only good to fair performance depending on the values of 𝜌S and 𝜌NS . A more detailed breakdown of the results discussed above is provided in Table 1. Evaluating the algorithm’s performance using a general set of algorithm parameters across multiple recording sessions In this subsection we evaluate the performance of our algorithm across multiple recording sessions from different GAERS using a ‘general’ set of algorithm parameters in each case. This general set of parameters is based on the average of the 𝛼 ∗ values listed in Table 1. Specifically, the average of the 𝛼 ∗ values when choosing 𝜌S = 𝜌NS = 0.0 or 0.75 is ≈ 0.075, we denote this average 𝛼 by 𝛼. Thus, we consider parameter setting 𝑃3 with 𝛿𝛼𝛽 = 0.01 and 𝛼 as our general set of algorithm parameters. Figure. 7 shows the change in performance of our algorithm in terms of the change in the distance from the ideal (DFI) 8/23
1.0 ρS = ρNS = 0.00 ρS = ρNS = 0.75
ROC-AUC
0.9 0.8 0.7 0.6 0.5 T2M
T4M
T5M
T8C
K2E
K2M K3A recording session
K4A
K5M
S4E
S4L
S5M
Figure 5. AUC of ROC curves constructed from different recording sessions when (in black) 𝜌S = 𝜌NS = 0 and (in red) 𝜌S = 𝜌NS = 0.75. Remainder of algorithm parameters are chosen as 𝑃3 and 𝛿𝛼𝛽 = 0.01 (see Table 3).
ROC: DFI(α∗)
0.3
ρS = ρNS = 0.00 ρS = ρNS = 0.75
0.2
0.1
0.0
T2M
T4M
T5M
T8C
K2E
K2M K3A recording session
K4A
K5M
S4E
S4L
S5M
Figure 6. DFI(𝛼 ∗ ) of ROC curves constructed from different recording sessions when (in green) 𝜌S = 𝜌NS = 0 and (in purple) 𝜌S = 𝜌NS = 0.75. Remainder of algorithm parameters are chosen as 𝑃3 and 𝛿𝛼𝛽 = 0.01 (see Table 3). when using the general set of algorithm (parameters in comparison to the optimal set of algorithm parameters for each recording ) session. More specifically, we plot DFI 𝛼 − DFI(𝛼 ∗ ) for each recording session when choosing 𝜌S = 𝜌NS = 0.0 (in green) and 𝜌S = 𝜌NS = 0.75 (in purple). Figure. 7 shows that, for many recording sessions, our algorithm maintains much of its performance when using the general set of parameters, i.e., there is little change in the DFI. This is true when 𝜌S = 𝜌NS = 0 or the more clinically relevant setting of 𝜌S = 𝜌NS = 0.75. However, for rat S, there is a more noticeable decrease in performance (larger DFI(𝛼) than DFI(𝛼 ∗ )). This is mostly due to the differences in seizure morphology shown in Fig. 9; the amplitudes of fluctuation in the voltage recordings for rat S are significantly smaller than those for rats T and K. This is further evidenced by the values of 𝛼 ∗ listed in Table 1, i.e., the values of 𝛼 ∗ for rat S are significantly smaller than those for rats T and K .
Discussion The present paper provides a more detailed analysis on the performance of the seizure detection algorithm introduced in our recent work, Flynn et al.10 . Specifically, we applied our algorithm to the same set of voltage recordings of seizure activity studied above, local field potential (LFP) recordings of seizure activity in Genetic Absence Epilepsy Rats from Strasbourg (GAERS), and, taking inspiration from Temko et al.5 and Mathieson et al.6 , quantified performance in terms of metrics derived from receiver-operating-characteristic (ROC) curves. These metrics show that our algorithm (i) achieves very good performance on average across the voltage recordings we considered when using the optimal set of algorithm parameters in each recording session and (ii) maintains much of its accuracy when using a general set of algorithm parameters applicable across all recording sessions. 9/23
Rec. sess. (𝑁S(𝐸) )
AUC; 𝜌S = 𝜌NS = 0.0 (𝛼 ∗ : 𝑁S(𝐴) , 𝑁TP , 𝑁TN , DFI)
AUC; 𝜌S = 𝜌NS = 0.75 (𝛼 ∗ : 𝑁S(𝐴) , 𝑁TP , 𝑁TN , DFI)
T2M (154)
0.9745 (0.08: 100, 94, 92, 0.0925)
0.9440 (0.08: 100, 91, 90, 0.1279)
T4M (135)
0.9247 (0.07: 101, 99, 87, 0.1182)
0.8899 (0.07: 101, 96, 86, 0.1386)
T5M (82)
0.8793 (0.085: 52, 47, 40, 0.2198)
0.8421 (0.085: 52, 46, 40, 0.2329)
T8C (194)
0.9545 (0.085: 150, 149, 137, 0.0749)
0.9508 (0.085: 150, 148, 137, 0.0764)
K2E (95)
0.7435 (0.075: 81, 62, 76, 0.2090)
0.6303 (0.075: 81, 62, 75, 0.2155)
K2M (112)
0.9266 (0.08: 118, 106, 96, 0.1992)
0.9155 (0.08: 118, 105, 92, 0.2287)
K3A (83)
0.9541 (0.075: 75, 66, 66, 0.1615)
0.7301 (0.075: 75, 64, 65, 0.1901)
K4A (77)
0.8841 (0.07: 58, 47, 53, 0.1889)
0.7521 (0.07: 58, 47, 51, 0.2105)
K5M (244)
0.7642 (0.085: 124, 94, 97, 0.3205)
0.6888 (0.08: 125, 86, 98, 0.3673)
S4E (127)
0.8387 (0.07: 69, 62, 44, 0.3110)
0.5057 (0.06: 98, 71, 77, 0.3402)
S4L (64)
0.9406 (0.06: 47, 46, 41, 0.1009)
0.8008 (0.06: 47, 46, 41, 0.1009)
S5M (144)
0.9078 (0.06: 102, 102, 90, 0.0973)
0.7753 (0.06: 102, 101, 89, 0.1068)
Table 1. AUC values computed from ROC curves where 𝜌S = 𝜌NS = 0.0 (middle column) and 𝜌S = 𝜌NS = 0.75 (right column) for specific recording sessions (left column) using parameter setting 𝑃3 with 𝛿𝛼𝛽 = 0.01 (see Table 3). Information in parenthesis in left column: the no. of seizure time intervals according to the expert. Information in parenthesis in middle and right columns: the optimal 𝛼 value (denoted by 𝛼 ∗ ), and the following corresponding to 𝛼 ∗ : the no. of seizure intervals according to the algorithm, the no. of these time intervals that were classified as TP, and the DFI(𝛼 ∗ ) (denoted by DFI). While our algorithm only monitors a single feature of the voltage recordings, the voltage at a given time, we expect performance to further improve by monitoring additional features that are relevant to the different types of critical transition (CT) associated with seizure onset, such as the variance and autocorrelation of the time series10 . Furthermore, for seizures that involve a bifurcation, it may be possible to detect precursors of seizure onset before a CT by monitoring features associated with ‘early warning signals of critical slowing down’, such as increases in variance and autocorrelation in advance of a CT. However, one must conduct further analysis to determine the predictive power of such features21, 22 . To widen our algorithm’s applicability, in future work we may consider adaptive threshold techniques to account for seizures whose amplitude either increases or decreases during the course of a seizure23 . More generally, our algorithm can be used to detect similar CTs in complex systems beyond the brain, such as CTs found in thermoacoustic systems24 and aircraft wings25 . Overall, our results demonstrate the benefit of studying the epileptic brain through the lens of CTs. We show that seizure detection algorithms can be designed through a mathematics-based approach, without reliance on machine learning tools, and that robust and accurate detection can be achieved despite the presence of artefacts, interictal epileptiform discharges, and variations in seizure morphology.
Methods M1: Obtaining and annotating seizure data from GAERS The data discussed here was obtained by Cian McCafferty, François David, and Vincenzo Crunelli and was presented in 26. Obtaining the data: The electrophysiological data was acquired from male Genetic Absence Epilepsy Rats from Strasbourg (GAERS) aged between 4-7 months when in a state of ‘relaxed wakefulness’ where seizures were experienced more often. Silicon-site electrodes were used to sample voltages (20000/second) from the ventrobasal thalamus while the rats were able to move freely and alternate between waking, sleeping, and seizing states. This data was processed with a Plexon HST/32V-G20 VLSI-based preamplifier and associated digitization system and subsequently down-sampled to 1000 samples/second. Labelling the recording session: Each recording session was labelled using the following convention described through example, ‘S1K’: voltage recordings from rat S on day 1 of recording during the K𝑡ℎ recording session on that day. Annotating the data: We denote seizure onset times by 𝜏1 and offset times by 𝜏2 . The process used by the expert to obtain these is specified below. Spike-wave discharges (SWDs) are the electrical hallmark of absence seizures, they define electrical seizure onset and offset times. These were identified using Cambridge Electronic Design’s ‘Spike2’ software and the following 10/23
change in DFI: DFI(α) − DFI(α∗)
ρS = ρNS = 0.00 ρS = ρNS = 0.75
0.20 0.15 0.10 0.05 0.00
T2M
T4M
T5M
T8C
K2E
K2M K3A K4A recording session
K5M
S4E
S4L
S5M
Figure 7. The change in DFI when using an upper threshold of 𝛼 over 𝛼 ∗ for 𝜌S = 𝜌NS = 0 (in black) and 𝜌S = 𝜌NS = 0.75 (in red). Remainder of algorithm parameters are chosen as 𝑃3 and 𝛿𝛼𝛽 = 0.01 (see Table 3). LFP [mV]
0.1
0.0
-0.1
0
5
10
15
20
25
t [s]
Figure 8. Example of seizure activity in GAERS. Time series of voltage recordings plotted in terms of local field potential, denoted by LFP vs. time in seconds. Vertical lines correspond to (orange) seizure onset and (purple) seizure offset times according to expert annotations. Portion of the voltage recordings shown here is taken from session ‘S5A’. procedure. In all cases, EEG at 1000Hz was used for this step. EEG was acquired at 1000Hz for fMRI (see McCafferty et al.27 ) and behaviour was down-sampled by averaging for neuronal activity. Briefly, smoothening (voltage at time 𝑡 is set to the mean of voltages from time 𝑡 − 10ms to time 𝑡 + 10ms) and DC removal (voltage at time 𝑡 is set to the original voltage minus the mean of voltages from time 𝑡 − 0.1s to time 𝑡 + 0.11s) functions were used to reversibly visually clean the frontoparietal differential EEG. Then, a negative amplitude threshold (mean voltage minus 5-7 standard deviations of baseline non-SWD EEG) was used to detect putative spike-wave crossing points, defined as whenever the signal crossed this amplitude threshold. The crossing points were then grouped into events based on the intervals between them (maximum time between initial two crossings 0.2s, maximum time between any two crossings within an event 0.35s, minimum of 5 crossings per event) and the defined properties of SWDs (minimum duration 0.5s, minimum inter-SWD interval 0.5s merging any SWDs with shorter intervals), and subsequently, using a frequency threshold, these events were classified as SWDs (if > 75% of intercrossing intervals were within a 5–12Hz range) or other (e.g., noise, sleep). Labelled SWDs were then visually inspected for accuracy of SWD detection, as well as seizure onset and offset times. Periods of sleep were identified based on sharp increases in the 1–4Hz frequency band and were excluded from analysis. Periods of non-REM sleep were rare in the recordings and not informative due to their relatively short duration. Examples and statistics of seizure activity: See Fig. 8 for an example of expert annotations of voltage recordings. See Fig. 9 for representative examples of seizure activity from different GAERS, illustrating differences in seizure morphology. Specifically, the average amplitude of fluctuation in the S state is much larger in Fig. 9 (a) than in Figs. 9 (b) and (c), a difference that is consistently observed across recording sessions from each of these GAERS. While Fig. 9 shows differences in how the S state emerges, this is not unique to individual GAERS. We show in Flynn et al.10 that the S state can emerge through different types of CT, with similar frequencies of each transition type observed across GAERS. See Fig. 10 for a comparison of the probability density of residence times in the S and NS states for all seizure onset and offset times annotated by the expert in the voltage recordings of rats T, K, and S. Note, there are similar distributions for each 11/23
(a) T8C
(b) K4A
(c) S4E
LFP [mV]
0.1
0.0
-0.1 -5
T [s]
0
-5
0
T [s]
-5
0
T [s]
probability density [arb. units]
Figure 9. Examples of voltage recordings chosen to highlight some of the differences in seizure morphology between different GAERS (recording session is specified in top left corner). LFP is plotted here vs. 𝑇 = 𝑡 − 𝜏1 . Vertical orange lines correspond to 𝑇 = 0, i.e., the seizure onset time according to the expert.
-1
(a)
S state NS state
(b)
(c)
-3
-5
100 101 102 103 Rat T residence times [s]
100 101 102 103 Rat K residence times [s]
100 101 102 103 Rat S residence times [s]
Figure 10. Probability density of residence times in the S and NS states computed from expert annotations of seizure onset and offset times for all voltage recordings of rat T in (a), rat K in (b), and rat S in (c). The information presented here is based on 621, 1593, and 1136 transitions from each respective rat. rat. The average residence time in the S and NS state for rat T was ≈ 11 s and ≈ 86 s, for rat K was ≈ 13 s and ≈ 76 s, and for rat S was ≈ 10 s and ≈ 74 s. In summary, when comparing the insights on seizure activity in Fig. 9 to Fig. 10, we observe much greater variability in terms of the amplitude of fluctuation than the temporal characteristics of seizure activity across different GAERS. Artefacts and how they are accounted for: In the voltage recordings we analysed, artefacts are high-amplitude and highfrequency oscillatory-like events; see Fig. 12 for an example. Artefacts appear for different reasons, the most common being (i) movement of electrodes/wiring and (ii) movement of muscles close to the electrodes (generally for chewing). Some artefacts are generated by signal overload and noise from electrical devices around the recording setup, however, these are less common as the recording process was generally appropriately amplified and shielded. Artefacts are accounted for through the annotation method described above since these events do not align with the frequency profile of SWDs. Additionally, the visual inspection described above also includes the manual rejection of artefacts which were not excised from annotation. Entire recording sessions were excluded if more than 5% of the session consisted of artefact activity. Remark (comparison with human data): GAERS can express multiple absence seizures per minute28 , far exceeding the frequency of any absence seizure syndrome in humans. For instance, Gregorčič et al.29 found that, for children with treatmentresistant childhood absence epilepsy, the median number of seizures per day was three. For a review of the GAERS model, with attention to its similarities and differences to human absence seizures, see Depaulis et al.30 . M2: CT detection algorithm applied to voltage recordings of seizure activity - technical details For convenience, we describe the algorithm in terms of a variable 𝑥(𝑡), which in our case is the LFP at a given time 𝑡, i.e., LFP(𝑡). We use 𝛿 > 0 to denote the time interval between two consecutive data points in a given time series, in the case of voltage recordings in GAERS, we have 𝛿 = 0.001s. We introduce the following six parameters that our algorithm uses: the upper voltage threshold 𝛼 > 0, the lower voltage threshold 0 < 𝛽 < 𝛼, the size of the moving window 𝜏𝑤 > 0, the time step 12/23
LFP [mV]
0.1
0
-0.1 9589
9590
9591
e t1
9592
9593
9594
9595
t [s]
Figure 11. Example of an almost-occurring CT (interictal epileptiform discharges). Vertical dashed line indicates (in red) when a CT from the NS to S state almost occurs at time 𝑡 = 𝑡̃1 . Algorithm parameters chosen as 𝛼 = 0.07, 𝛽 = 0.06, 𝜏NS = 3, 𝜏S = 2, 𝜏𝑤 = 1, and Δ = 0.001. Horizontal lines indicate the thresholds of (in red) 𝛼 and (in green) 𝛽. size of the moving window 𝛿 ≤ Δ ≤ 𝜏𝑤 , the minimum time duration of larger-amplitude fluctuations, 𝜏S > 𝜏𝑤 , expressed as 𝜏S = 𝑛S Δ + 𝜏𝑤 , where 𝑛S is an integer, and the minimum time duration of smaller-amplitude fluctuations, 𝜏NS ≥ 𝜏S , expressed similarly as 𝜏NS = 𝑛NS Δ + 𝜏𝑤 , where 𝑛NS is an integer. In all our experiments we set Δ = 𝛿 = 0.001. The starting point: We start in the NS state, where |𝑥(𝑡)| < 𝛽 for the time duration of at least 𝜏NS . The moving window: When the brain is in the NS state and |𝑥(𝑡)| exceeds 𝛼 at time 𝑡 = 𝑡𝑗 , the moving window is activated and |𝑥(𝑡)| is examined within consecutive windows of duration 𝜏𝑤 that are shifted in time by Δ, starting with [𝑡𝑗 , 𝑡𝑗 + 𝜏𝑤 ], then [𝑡𝑗 + Δ, 𝑡𝑗 + 𝜏𝑤 + Δ], [𝑡𝑗 + 2 Δ, 𝑡𝑗 + 𝜏𝑤 + 2 Δ], and so on. The moving window is deactivated in two cases: (i) the brain is in the NS state, the window is activated, but the algorithm does not detect a CT to the S state, and (ii) the brain is in the S state and the algorithm detects a CT to the NS state. Critical transitions: The algorithm detects a CT from the NS to S state at time 𝑡 = 𝑡1 if: (a1) The brain is in the NS state just before 𝑡1 . (a2) |𝑥(𝑡)| exceeds 𝛼 at time 𝑡 = 𝑡1 , i.e., |𝑥(𝑡1 )| = 𝛼 and |𝑥(𝑡1 + 𝛿)| > 𝛼. (a3) Each of the 𝑛S consecutive positions of the moving window contains an |𝑥(𝑡)| > 𝛽. The algorithm detects a CT from the S to NS state at time 𝑡 = 𝑡2 if: (b1) The brain is in the S state just before 𝑡2 . (b2) |𝑥(𝑡)| falls below 𝛽 at time 𝑡 = 𝑡2 , i.e., |𝑥(𝑡2 )| = 𝛽 and |𝑥(𝑡2 + 𝛿)| < 𝛽. (b3) Each of the 𝑛NS consecutive positions of the moving window contains no |𝑥(𝑡)| ≥ 𝛼. In other words, the algorithm detects a CT from the NS to S state if |𝑥(𝑡)| exceeds the upper threshold 𝛼 and then continues to exceed the lower threshold 𝛽 frequently enough for a period of at least 𝜏S . Similarly, the algorithm detects a CT from the S to NS state if |𝑥(𝑡)| falls below the lower threshold 𝛽 and then does not exceed the upper threshold 𝛼 for a period of at least 𝜏NS . Almost-occurring critical transitions: The algorithm detects an almost-occurring CT from the NS to S state at time 𝑡 = 𝑡̃1 if (a1) and (a2) are satisfied but (a3) is not. Similarly, the algorithm detect an almost-occurring CT from the S to NS state at time 𝑡 = 𝑡̃2 if (b1) and (b2) are satisfied but (b3) is not. See Fig. 11 for an example of an almost-occurring CT from the NS to S state. Note, events like these are associated with interictal epileptiform discharges, short time intervals of seizure-like activity that do not meet the criteria to be considered as seizure activity. Accounting for artefacts
Voltage recordings of seizure activity are notoriously susceptible to artefacts. The data we analyse in this paper is no exception, see Fig. 12 for a typical example of artefact activity and part M1 of the Methods section for reasons why artefacts appear. Artefacts can be mistaken for seizure activity and seizure detection algorithms need to be designed accordingly. We now describe how we alter the above algorithm to account for artefacts. From Fig. 12 we observe that for 𝑡 ∈ [2, 5], the LFP jumps between −0.15 and 0.15 mV in short time intervals, often within consecutive measurements (0.001s). Based on this observation, we alter condition (a2) of our algorithm to account for artefacts as follows: 13/23
Artefact LFP [mV]
0.1
0
-0.1 t̂1 0
2
t̂2 4
6
8
t [s]
Figure 12. Example of artefact activity in voltage recordings. Vertical dashed lines indicate (in pink) the beginning and (in grey) end of artefact activity. Vertical solid lines indicate seizure onset time according to (in orange) the expert and (in red) the algorithm when applied to a portion of the voltage recordings from session ‘T8C’. Algorithm parameters chosen as 𝛼 = 0.055, 𝛽 = 0.04, 𝜏NS = 5, 𝜏S = 2, 𝜏𝑤 = 1, and Δ = 0.001. Horizontal lines indicate the thresholds of (in red) 𝛼 and (in green) 𝛽. [ ] (a2-1) |𝑥(𝑡)| exceeds 𝛼 at time 𝑡 = 𝑡′1 and |𝑥(𝑡 + 𝛿) − 𝑥(𝑡)| < 𝜉 for all 𝑡 ∈ 𝑡′1 , 𝑡′1 + 𝑡𝑤 , or, [ ] (a2-2) |𝑥(𝑡)| exceeds 𝛼 at time 𝑡 = 𝑡′1 and |𝑥(𝑡 + 𝛿) − 𝑥(𝑡)| ≥ 𝜉 for any 𝑡 ∈ 𝑡′1 , 𝑡′1 + 𝑡𝑤 . If (a2-1) is true then the algorithm continues as before to evaluate whether a CT from the NS to S state is detected at 𝑡′1 . On the other hand, if (a2-2) is true then we say artefact activity begins at time 𝑡̂1 = 𝑡′1 and we continue to monitor |𝑥(𝑡)| in moving [ ] windows. We say artefact activity ends at time 𝑡 = 𝑡̂2 if |𝑥(𝑡̂2 )| < 𝛼 and |𝑥(𝑡 + 𝛿) − 𝑥(𝑡)| < 𝜉 for all 𝑡 ∈ 𝑡̂2 , 𝑡̂2 + 𝑡𝑤 . We find most artefacts are accounted for by choosing 𝜉 = 0.2. Fig. 12 shows the result of applying our algorithm with (a2-1) and (a2-2) and parameters chosen as 𝛼 = 0.055, 𝛽 = 0.04, 𝜏NS = 5, 𝜏S = 2, 𝜏𝑤 = 1, and Δ = 0.001. The example in Fig. 12 is chosen to show that for the current choice of 𝜏S and 𝜏NS , the time when artefact activity begins would have been considered as the time when a CT from the NS to S state occurs unless artefacts are accounted for. Furthermore, Fig. 12 shows that by excluding these artefacts there is still a strong agreement between the expert and the algorithm on the seizure onset times. M3: Seizure and non-seizure time interval classification procedure Seizure and non-seizure time intervals
We introduce the following terminology to define seizure and non-seizure time intervals in time series of voltage recordings based on expert annotations of seizure onset and offset times and the times that our algorithm detects CTs. Expert: We denote the seizure onset and offset times according to the expert as 𝜏1(𝑖) and 𝜏2(𝑖) where 𝑖 = 1, 2, … , 𝑁S(𝐸) and
𝑁S(𝐸) is the number of seizure time intervals in the time series according to the expert. We define the 𝑖𝑡ℎ seizure time in( ) { } terval according to the expert as 𝐼S(𝐸) = 𝑡 ∈ ℝ ∶ 𝜏1(𝑖) ≤ 𝑡 < 𝜏2(𝑖) . We define the subsequent non-seizure time interval as ( ) { } 𝑖 (𝐸) (𝐸) 𝐼NS = 𝑡 ∈ ℝ ∶ 𝜏2(𝑖) ≤ 𝑡 < 𝜏1(𝑖+1) . Thus, for a given 𝑁S(𝐸) , there are 𝑁NS = 𝑁S(𝐸) − 1 non-seizure time intervals. 𝑖
Algorithm: We denote the times that our algorithm detects CTs from the NS to S state as 𝑡(𝑗) and CTs from the S to NS state as 1
𝑡(𝑗) where 𝑗 = 1, 2, … , 𝑁S(𝐴) and 𝑁S(𝐴) is the number of seizure time intervals in the time series according to the algorithm. We 2 ( ) { } ≤ 𝑡 < 𝑡(𝑗) use the same convention to define the 𝑗 𝑡ℎ seizure time interval according to the algorithm as 𝐼S(𝐴) = 𝑡 ∈ ℝ ∶ 𝑡(𝑗) . 1 2 𝑗 ( ) { } (𝐴) We define the subsequent non-seizure time interval as 𝐼NS = 𝑡 ∈ ℝ ∶ 𝑡(𝑗) ≤ 𝑡 < 𝑡(𝑗+1) . Similarly, for a given 𝑁S(𝐴) , there 2 1 𝑗
(𝐴) are 𝑁NS = 𝑁S(𝐴) − 1 non-seizure time intervals.
)| ( ) |( (𝑗) || (𝑗+1) (𝐴) || − 𝑡 𝐼 − 𝑡(𝑗) The length of the above time intervals is calculated as follows, || 𝐼S(𝐴) || = 𝑡(𝑗) , , and simi| NS 𝑗 | = 𝑡1 2 1 2 𝑗| | | | ( ) ( ) | | | (𝐸) | | using the corresponding 𝜏1 and 𝜏2 terms, where |𝐼| denotes the length of the interval 𝐼. We larly for || 𝐼S(𝐸) || and || 𝐼NS | 𝑖| 𝑖| | | )| |( (𝐴) ) | |( | as a residence time in the S state according to the algorithm, | 𝐼 (𝐴) | as a residence time in the NS refer to a given || 𝐼S | | | 𝑗| | | NS 𝑗 | 14/23
)| |( |( (𝐸) ) | | in terms of the expert’s annotations. state, and similarly for || 𝐼S(𝐸) || and || 𝐼NS | 𝑖| 𝑖| | | Classification of 𝐴S and 𝐴NS
The metrics used to evaluate the performance ( ) of our algorithm are based on the following time interval classification procedure: (𝐴) True positive: we classify a given 𝐼S as true positive (TP) if any of the following three conditions are true: C1: C2:
( (
𝐼S(𝐴) 𝐼S(𝐴)
)
(
𝑗
)
) ( ) ( ) is contained in some 𝐼S(𝐸) , that is 𝐼S(𝐴) ⊆ 𝐼S(𝐸) for some 𝑖. 𝑖
(
𝑗
𝑗
𝑗
)
(
contains at least one 𝐼S(𝐸) , that is 𝐼S(𝐴)
)
𝑖
𝑖
(
𝑗
)
⊃ 𝐼S(𝐸)
𝑖
for some 𝑖.
( ) ( ) C3: Neither C1 or C2 are true but there is sufficient overlap between 𝐼S(𝐴) and some 𝐼S(𝐸) , 𝑗 𝑖 ( ) ⋂( ) ( ) (𝐴) (𝐸) (𝐸) that is | 𝐼S 𝐼S |∕| 𝐼S | > 𝜌S for some 𝑖 and a given 0 ≤ 𝜌S ≤ 1. 𝑗
𝑖
𝑖
) False positive: if C1-C3 are not true then 𝐼S(𝐴) is classified as false positive (FP). 𝑗 ( ) (𝐴) True negative: we classify a given 𝐼NS as true negative (TN) if any of the following three conditions are true: D1: D2:
( (
(𝐴) 𝐼NS (𝐴) 𝐼NS
)
(
𝑗
𝑗
)
𝑖
(
)
𝑗
) ( ) ( (𝐸) (𝐴) (𝐸) is contained in some 𝐼NS , that is 𝐼NS ⊆ 𝐼NS for some 𝑖. (
(𝐸) contains one 𝐼NS
) 𝑖
𝑖
𝑗
(
and no 𝐼S(𝐸)
)
(
𝑘
(
(𝐴) , that is 𝐼NS
) 𝑗
( ) (𝐸) for some 𝑖. ⊃ 𝐼NS 𝑖
( ) (𝐸) (𝐴) D3: Neither D1 or D2 are true but there is an 𝐼NS with sufficient overlap with 𝐼NS , 𝑗 𝑖 ( ) ⋂( ) ( ) (𝐴) (𝐸) (𝐴) that is | 𝐼NS 𝐼NS |∕| 𝐼NS | > 𝜌NS for some 𝑖 and a given 0 ≤ 𝜌NS ≤ 1 and 𝑗 𝑖 ( ) ( )𝑗 ( ) ( ) (𝐴) (𝐴) 𝐼NS does not contain any 𝐼S(𝐸) , i.e., 𝐼NS ⊅ 𝐼S(𝐸) for all 𝑘. 𝑗
𝑘
(
(𝐴) False negative: if D1-D3 are not true then 𝐼NS
)
) 𝑗
𝑗
𝑖
is classified as false negative (FN).
We denote the total number of TPs with 𝑁TP , FPs with 𝑁FP , TNs with 𝑁TN , and FNs with 𝑁FN . Table 2 specifies how we check C1-C3 and D1-D3 in terms of the corresponding 𝑡1 , 𝑡2 , 𝜏1 , and 𝜏2 values. Note, in C3 and D3 we introduced the parameters, 𝜌S and 𝜌NS , which we use as thresholds to quantify the minimum amount (𝐴) of overlap between time intervals that is needed for a given 𝐼S(𝐴) or 𝐼NS to be classified as TP or TN. These parameters can be chosen in accordance with standards set by clinicians, for instance, in Mathieson et al.6 the authors set their parameters to 0.75. ( ) ( ) Remark: based on C2, we allow for a given 𝐼S(𝐴) to be classified as TP even if it contains more than one 𝐼S(𝐸) 𝑗 𝑖 ( ) (𝐸) and an 𝐼NS . However, according to certain performance metrics, this can result in misleadingly high performance for 𝑘 ( ) (𝐴) reasons that we outline in part M4 of the Methods Section. In contrast to C2, with D2 we only allow for a given 𝐼NS to be 𝑗 ( ) (𝐸) classified as TN if it contains no more than one 𝐼NS . 𝑖
M4: Mitigating the occurrence of misleading points on ROC curves In this subsection we (i) show that certain parameter settings of our algorithm can lead to misleading points on an ROC curve and (ii) present a method to prevent these misleading points from appearing. How misleading points occur: Misleading points appear when the algorithm detects very few seizure and non-seizure time intervals in comparison to the expert and these time intervals satisfy the TP (C1-C3) and TN (D1-D3) conditions. As a result, the corresponding point on the ROC curve incorrectly indicates that the algorithm performed well. This typically occurs when the upper voltage threshold, 𝛼, is set too small or too large, causing the algorithm to detect one or few seizure time intervals whose total duration is disproportionally long or short in comparison to the time series the algorithm is applied to. For example, 15/23
C1
𝜏1(𝑖) ≤ 𝑡(𝑗) < 𝑡(𝑗) ≤ 𝜏2(𝑖) . 1 2
C2
𝑡(𝑗) ≤ 𝜏1(𝑖) < 𝜏2(𝑖) ≤ 𝑡(𝑗) . 1 2
C3
𝑡(𝑗) ≤ 𝜏1(𝑖) < 𝑡(𝑗) ≤ 𝜏2(𝑖) and (𝑡(𝑗) − 𝜏1(𝑖) ) > 𝜌S (𝜏2(𝑖) − 𝜏1(𝑖) ) or, 1 2 2
D1
𝜏2(𝑖) ≤ 𝑡(𝑗) < 𝑡(𝑗+1) ≤ 𝜏1(𝑖+1) . 2 1
D2
𝜏1(𝑖) ≤ 𝑡(𝑗) ≤ 𝜏2(𝑖) and 𝜏1(𝑖+1) ≤ 𝑡(𝑗+1) ≤ 𝜏2(𝑖+1) . 2 1
D3
𝜏1(𝑖) ≤ 𝑡(𝑗) ≤ 𝜏2(𝑖) and 𝜏2(𝑖) ≤ 𝑡(𝑗+1) ≤ 𝜏1(𝑖+1) and (𝑡(𝑗+1) − 𝜏2(𝑖) ) > 𝜌NS (𝑡(𝑗+1) − 𝑡(𝑗) ) or, 2 1 1 1 2
𝜏1(𝑖) ≤ 𝑡(𝑗) < 𝜏2(𝑖) ≤ 𝑡(𝑗) and (𝜏2(𝑖) − 𝑡(𝑗) ) > 𝜌S (𝜏2(𝑖) − 𝜏1(𝑖) ). 1 2 1
𝜏2(𝑖) ≤ 𝑡(𝑗) ≤ 𝜏1(𝑖+1) and 𝜏1(𝑖+1) ≤ 𝑡(𝑗+1) ≤ 𝜏2(𝑗+1) and (𝜏1(𝑖+1) − 𝑡(𝑗) ) > 𝜌NS (𝑡(𝑗+1) − 𝑡(𝑗) ). 2 1 2 1 2
Table 2. Implementation of conditions used to classify seizure (C1-C3) and non-seizure (D1-D3) time intervals. Parameter setting
𝜏NS
𝜏S
𝜏𝑤
𝛿𝛼𝛽
𝛼
𝛽
𝜌S and 𝜌NS
𝑃1
0.7
1
0.5
[0.01, 0.06, 0.01]
[𝛿𝛼𝛽 + 0.01, 0.1, 0.005]
𝛽 = 𝛼 − 𝛿𝛼𝛽
[0, 1, 0.025]
𝑃2
2
1
1
"
"
"
"
𝑃3
3
2
1
"
"
"
"
Table 3. Parameter settings used when evaluating the performance of the CT detection algorithm described in part M2 of the Methods section. Values specified in square brackets correspond to [lower value, upper value, difference between each value]. it can happen that a given 𝐼S(𝐴) contains almost all the 𝐼S(𝐸) . We illustrate an example of such a scenario in Fig. 13 when the algorithm is used to detect CTs between NS and S states in the T2M voltage recordings, a portion of which is shown in (a). In this example the algorithm parameters are chosen as 𝛼 = 0.03, 𝛽 = 0.02, 𝜏NS = 2, 𝜏S = 1, 𝜏𝑤 = 1, and the time intervals are classified with overlap parameters set as 𝜌S = 𝜌NS = 0. In (b) and (c) we show a binary sequence representation of the seizure and non-seizure time intervals according to expert annotations and the algorithm. Beneath each sequence we use the same coloured box convention in Fig. 2 to indicate how the corresponding seizure and non-seizure time intervals are classified. In this example, the expert identified 154 seizure time intervals, i.e., 𝑁S(𝐸) = 154, while the algorithm identified only 4, i.e., 𝑁S(𝐴) = 4.
(𝐸) However, these 4 𝐼S(𝐴) contain most of the 𝐼S(𝐸) and 𝐼NS . Without adapting the classification procedure this results in (FPR, TPR) = (0,1), indicating that the algorithm achieves ideal performance. However, this is not the case in reality. Mitigating the occurrence of misleading points: To combat the issue described above, we construct our ROC curve using (FPR, TPR) points whose corresponding time intervals satisfy the following condition:
•
∑𝑁S(𝐴) ||( (𝐴) ) || ∑𝑁S(𝐸) ||( (𝐸) ) || 𝐼 ≤ ℎ 𝑢 𝑖=1 | 𝐼S | | for a suitably chosen ℎ𝑢 > 0. 𝑗=1 | S 𝑗| 𝑖| | |
In other words, the sum of residence times in the S state according to the algorithm must be less than or equal to a multiple, ℎ𝑢 , of the sum of residence times in the S state according to the expert. It was found empirically that by setting ℎ𝑢 = 2, this additional step to our classification procedure removes the misleading points on a given ROC curve. The result of applying this additional step for different choices of ℎ𝑢 is illustrated in Fig. 14. More specifically, we plot ROC curves based on applying our algorithm to recording session T2M with parameters setting 𝑃2 and 𝛿𝛼𝛽 = 0.02 and classifying the resulting seizure and non-seizure time intervals with 𝜌S = 𝜌NS = 0. Figure 14 (a) shows the ROC curve obtained when no steps are taken to mitigate the occurrence of misleading points, the (FPR, TPR) point obtained from using 𝛼 = 0.03 incorrectly indicates the algorithm achieves ideal performance. Figure 14 (b) shows the ROC curve obtained when setting ℎ𝑢 = 4, the misleading points corresponding to 𝛼 = 0.03 and 0.035 are removed. Figure 14 (c) shows that by setting ℎ𝑢 = 2, the misleading points corresponding to 𝛼 ∈ [0.03, 0.055] are removed. This also improves the overall smoothness of the ROC curve. What is 16/23
LFP [mV]
0.1 (a) 0.0
-0.1 S (b) Expert
NS S
(c)
Algorithm ρS = ρNS = 0
TP
TN
NS 4200
4210
4220
4230
4240
4250
t [s]
Figure 13. Example of misleading TP and TN classifications. (a) shows a portion of voltage recordings. (b) and (c) show binary sequence representations of the non-seizure and seizure time intervals in (a) according to expert annotations and our algorithm. The boxes beneath these sequences are coloured green for TN and red for TP. Algorithm parameters chosen as 𝛼 = 0.03, 𝛽 = 0.02, 𝜏NS = 2, 𝜏S = 1, 𝜏𝑤 = 1, with 𝜌S = 𝜌NS = 0 used in the classification of the seizure and non-seizure intervals. common across (a)-(c) is that the AUC is not highly sensitive to this method for ℎ ≥ 2. This is a desirable result as the sole motivation behind using this additional constraint was to mitigate the influence of these misleading points.
References 1. World Health Organisation. Online Article: Epilepsy. https://www.who.int/news-room/fact-sheets/detail/epilepsy, Last accessed on 01-06-2026. 2. Islam, M. K., Rastegarnia, A. & Yang, Z. Methods for artifact detection and removal from scalp EEG: A review. Neurophysiol. Clinique/Clinical Neurophysiol. 46, 287–305 (2016). 3. de Curtis, M., Jefferys, J. G. & Avoli, M. Interictal epileptiform discharges in partial epilepsy. Jasper’s Basic Mech. Epilepsies [Internet]. 4th edition (2012). 4. Niedermeyer, E. & da Silva, F. L. Electroencephalography: basic principles, clinical applications, and related fields (Lippincott Williams & Wilkins, 2005). 5. Temko, A., Thomas, E., Marnane, W., Lightbody, G. & Boylan, G. Performance assessment for EEG-based neonatal seizure detectors. Clin. Neurophysiol. 122, 474–482 (2011). 6. Mathieson, S. R. et al. Validation of an automated seizure detection algorithm for term neonates. Clin. Neurophysiol. 127, 156–168 (2016). 7. O’Hagan, J., Keane, A. & Flynn, A. Confabulation dynamics in a reservoir computer: Filling in the gaps with untrained attractors. Chaos 35, 093130 (2025). 8. Sun, F. T., Morrell, M. J. & Wharen, R. E. Responsive cortical stimulation for the treatment of epilepsy. Neurotherapeutics 5, 68–74 (2008). 9. Morrell, M. J. Responsive cortical stimulation for the treatment of medically intractable partial epilepsy. Neurology 77, 1295–1304 (2011). 17/23
1.0 (a)
(b)
(c)
TPR
0.8 0.6 0.4
0.0 0.0
hu = 4 : α∗ = 0.085 AUC = 0.9509
No adjustment : α∗ = 0.03 AUC = 0.9583
0.2 0.2
0.4 0.6 FPR
0.8
1.0 0.0
0.2
0.4 0.6 FPR
0.8
1.0 0.0
hu = 2 : α∗ = 0.085 AUC = 0.9528 0.2
0.4 0.6 FPR
0.8
1.0
Figure 14. Effect of ℎ𝑢 on removing misleading points from an ROC curve. (a) ROC curve without applying the ℎ𝑢 adjustment. (b) and (c) ROC curves after applying the ℎ𝑢 adjustment with ℎ𝑢 = 4 and ℎ𝑢 = 2, respectively. All ROC curves are generated using the same classifications of seizure and non-seizure time intervals in Fig. 3 for parameter setting 𝑃2 with 𝛿𝛼𝛽 = 0.02. In each panel, (FPR, TPR) points plotted in red correspond to different values of 𝛼, the green point corresponds to 𝛼 ∗ (the optimal 𝛼), the lower right corner specifies 𝛼 ∗ and the AUC. 10. Flynn, A. et al. Classifying seizure generation mechanisms: A critical transitions framework. arXiv:2511.20522 (2025).
arXiv preprint
11. Beniczky, S. et al. Updated classification of epileptic seizures: Position paper of the International League Against Epilepsy. Epilepsia 66, 1804–1823 (2025). 12. Kaye, D. et al. Impact of Prolonged Seizures on Patients’ and Caregivers’ Quality of Life (P1-9.014). In Neurology, vol. 104, 2708 (Lippincott Williams & Wilkins Hagerstown, MD, 2025). 13. Ashwin, P., Perryman, C. & Wieczorek, S. Parameter shifts for nonautonomous systems in low dimension: bifurcation-and rate-induced tipping. Nonlinearity 30, 2185 (2017). 14. Lenton, T. M. et al. Tipping elements in the Earth’s climate system. Proc. Natl. Acad. Sci. 105, 1786–1793 (2008). 15. Lenton, T. M. Tipping positive change. Philos. Trans. R. Soc. B 375, 20190123 (2020). 16. Armstrong McKay, D. I. et al. Exceeding 1.5 C global warming could trigger multiple climate tipping points. Science 377, eabn7950 (2022). 17. Lenton, T. M. et al. Global tipping points report 2025 (2025). University of Exeter. 18. Krasnosel’skii, M. A. & Pokrovskii, A. V. Systems with Hysteresis (Springer Science & Business Media, 2012). 19. Fawcett, T. An introduction to ROC analysis. Pattern recognition letters 27, 861–874 (2006). 20. Flynn, A. et al. Supplementary material from: “Detecting seizure onset and offset times using human intelligence: A critical-transitions-based approach”. (See attached document). 21. Lehnertz, K. Time-series-analysis-based detection of critical transitions in real-world non-autonomous systems. Chaos 34 (2024). 22. Ashwin, P., Bastiaansen, R., von der Heydt, A. S. & Ritchie, P. D. Early warning skill, extrapolation and tipping for accelerating cascades. Proc. Roy. Soc. A 481, 20250405 (2025). 23. Engel, J., Pedley, T. A. & Aicardi, J. Epilepsy: a comprehensive textbook, vol. 3 (Lippincott Williams & Wilkins, 2007). 24. Etikyala, S. & Sujith, R. Change of criticality in a prototypical thermoacoustic system. Chaos 27 (2017). 25. Ma, J. et al. Predicting tipping phenomenon in a conceptual airfoil structure under extreme flight environment. J. Sound Vib. 618, 119306 (2025). 26. McCafferty, C. et al. Cortical drive and thalamic feed-forward inhibition control thalamic output synchrony during absence seizures. Nat. neuroscience 21, 744–756 (2018). 27. McCafferty, C. et al. Decreased but diverse activity of cortical and thalamic neurons in consciousness-impairing rodent absence seizures. Nat. Commun. 14, 117 (2023). 18/23
28. Powell, K. L. et al. Seizure expression, behavior, and brain morphology differences in colonies of Genetic Absence Epilepsy Rats from Strasbourg. Epilepsia 55, 1959–1968 (2014). 29. Gregorčič, S. et al. Difficult to treat absence seizures in children: A single-center retrospective study. Front. Neurol. 13, 958369 (2022). 30. Depaulis, A., David, O. & Charpier, S. The genetic absence epilepsy rat from strasbourg as a model to decipher the neuronal and network mechanisms of generalized idiopathic epilepsies. J. Neurosci. Methods 260, 159–174 (2016).
Acknowledgements This publication has emanated from research conducted with the financial support of Taighde Éireann – Research Ireland under grant number [19/FFP/6782].
Author contributions statement A.F. and S.W. conceived the experiments, A.F. conducted the experiments, A.F. and S.W. analysed the results, C.Mc C., F.D., and V.C. provided the voltage recordings, C.Mc C., K.L., and G.L. provided valuable discussions. A.F., C.Mc C., K.L., G.L., and S.W. reviewed the manuscript.
Additional information The authors have no competing interests to disclose.
19/23
(a) τNS 16 12
AUC (b) τNS
DFI
0.96
0.3
0.92
8
16 12
0.2
8 0.88
4 2
4
6
0.1 4 2
τS
(c) τNS 16 12
6
τS
0
D̃|I| (d) τNS
D̃N
0.3
0.4
0
8
4
-0.3
4
16 12
0
8
-0.4
4 -0.6 2
4
6
τS
-0.8 2
4
6
τS
Figure S-1. Extension to Fig. 3 in [1]; illustrating how the values of 𝜏NS and 𝜏S influence the algorithm’s performance in ̃|𝐼| , and (d) 𝐷 ̃|𝑁| . terms of the following metrics (a) AUC, (b) DFI, (c) 𝐷
Supplementary information S1: Wider study on how 𝜏NS and 𝜏S influence the algorithm’s performance In Fig. S-1 we extend the results shown in Fig. 3 in the main text to show how different choices of 𝜏S ∈ [1, 8] and 𝜏NS ∈ [1, 23] influence the algorithm’s performance in terms of (a) the AUC, (b) the DFI(𝛼 ∗ ), (c) the relative difference between total time ̃|𝐼| (see Eq. (1)), and (d) the relative difference between total number of seizure time intervals, spent in the S state, denoted by 𝐷 ̃𝑁 (see Eq. (2)). Note, we compute 𝐷 ̃|𝐼| and 𝐷 ̃𝑁 based on ||𝐼 (𝐴) || and 𝑁 (𝐴) which are computed using the values of denoted by 𝐷 S | S | 𝑡1 and 𝑡2 obtained when using the optimal 𝛼 (i.e., 𝛼 ∗ ) for recording session T2M for a given choice of 𝜏S and 𝜏NS . To maintain ̃|𝐼| and 𝐷 ̃𝑁 is to indicate that the connection to Fig. 3 in [1], we keep 𝛿𝛼𝛽 = 0.01, 𝜏𝑤 = 1, and 𝜌S = 𝜌NS = 0. Note, the tilde in 𝐷 these are relative distances. Figures S-1 (a) and (b) show that while AUC > 0.97 and DFI < 0.1 for many choices of 𝜏NS when 𝜏S = 1, 2. In contrast Figs. S-1 (c) and (d) show that for larger 𝜏NS values, the total time spent in the S state according to the algorithm is significantly greater than according to the expert, and correspondingly, the algorithm detects significantly less seizure time intervals than the expert. Thus, the results presented in Fig. S-1 further inform our choice of 𝜏S = 2 and 𝜏NS = 3 to be optimal, reflecting a balance between performance and detecting a similar amount of seizure time intervals to the expert and that these intervals are of a similar lengths. ∑𝑁S(𝐴) ||( (𝐴) ) || ∑𝑁S(𝐸) ||( (𝐸) ) || 𝐼 | − 𝑖=1 | 𝐼S 𝑖 | 𝑗=1 | S 𝑗| | | |. ̃ 𝐷|𝐼| = ) (𝐸) |( | 𝑁 ∑ S | (𝐸) | 𝐼 | 𝑖=1 | S 𝑖| | 𝑁 (𝐴) − 𝑁S(𝐸) ̃𝑁 = S 𝐷 . 𝑁S(𝐸)
(1)
(2)
20/23
1
(a) P1
(b) P2
TPR FPR
0.8
(c) P3
0.6 0.4 0.2 0 0.02
0.04
0.06
0.08
α 0.02
0.04
0.06
0.08
α 0.02
0.04
0.06
0.08
α
Figure S-2. 𝛼 vs. TPR (in green) and FPR (in red) for parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c).
500
(a) P1
(b) P2
(E)
NS
NS , NTP, NFP
400
(c) P3
200
200
100
100
(A)
300 200 100 0 0.02
0.04
0.06
0.08
α
0 0.02
0.04
0.06
0.08
α
0 0.02
0.04
0.06
0.08
α
Figure S-3. 𝛼 vs. 𝑁S(𝐴) (in blue), 𝑁TP (in red), and 𝑁FP (in black) for parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c) and 𝛿𝛼𝛽 = 0.01. The green vertical lines corresponds to the optimal 𝛼. Dashed horizontal blue lines indicate 𝑁S(𝐸) for the T2M time series. Data points to the left of the red vertical line do not contribute to the corresponding ROC curves in Fig. 3 in [1].
S2: Analysis of quantities used to generate ROC curves in Fig. 3 in [1] Figures S-2-S-4 provide a more detailed breakdown on the different quantities that are used to construct the ROC curves in Fig. 3 in [1]. Specifically, we show how these quantities vary with respect to 𝛼 in parameter settings 𝑃1 , 𝑃2 , and 𝑃3 in the respective panels (a), (b), and (c) in Figs. S-2-S-4. Figure S-2 shows how the TPR and FPR vary with respect to 𝛼. Figure S-3 shows how 𝑁S(A) , 𝑁TP , and 𝑁FP vary with respect to 𝛼 and how these quantities compare to 𝑁S(E) . Similarly, Figure S-4 shows
(A) (E) how 𝑁NS , 𝑁TN , and 𝑁FN vary with respect to 𝛼 and how these quantities compare to 𝑁NS . Further, in Figs. S-3-S-4, the purpose of the red vertical line is to emphasise that all information to the left of the line is not used to construct ROC curves, and the purpose of the green vertical line is to indicate the quantities which correspond to the optimal choice of 𝛼, denoted by 𝛼 ∗ . Figures S-2-S-4 show there is a nonlinear relationship between 𝛼 and each quantity mentioned above. We find that while the algorithm achieves its best performance for 𝑃3 , the algorithm detects less seizure time intervals in comparison to using parameter sets 𝑃2 and 𝑃1 . This is to be expected since 𝜏S is smaller for 𝑃2 and 𝑃1 . Furthermore, the motivation to exclude points from the ROC curves (see part M4 of the Methods section in [1]) becomes clearer from Figs. S-3 and S-4 as these points correspond to cases where the algorithm detects a much larger number of seizure and non-seizure time intervals than the expert.
S3: Assessing agreement beyond ROC curves In this subsection, we examine the agreement between our algorithm and the expert’s annotations in ways beyond that captured by ROC curves, specifically the ROC curves shown in Fig. 3 in [1]. S3 (a): Comparing 𝑡1 and 𝜏1
In Fig. S-5 we compare the times that our algorithm detects CTs from the NS to S state, 𝑡1 , to the expert annotation of seizure onset times, 𝜏1 , for 𝐼S(A) that were classified as TPs when using parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c). For each TP, we compute the time difference 𝑇1 = 𝜏1 − 𝑡1 and show, for each value of 𝛼, the number of TPs associated with each 𝑇1 ∈ [−5, 5]. Darker shades of red indicate a larger number of TPs. 21/23
(a) P1
500
(b) P2
(E)
NNS
NNS , NTN, NFN
400
(c) P3
200
200
100
100
(A)
300 200 100 0 0.02
0.04
0.06
α
0.08
0 0.02
0.04
0.06
0.08
0 0.02
α
0.04
0.06
α
0.08
(𝐴) Figure S-4. 𝛼 vs. 𝑁NS (in blue), 𝑁TN (in red), and 𝑁FN (in black) for parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c) and
(𝐸) 𝛿𝛼𝛽 = 0.01. The green vertical lines corresponds to the optimal 𝛼. Dashed horizontal blue lines indicate 𝑁NS for the T2M time series. Data points to the left of the vertical red line do not contribute to the corresponding ROC curves in Fig. 3 in [1].
T1 (a) 2
# T1 (b)
# T1 (c)
#
25
25
25
20
2
20
2
20
0
15 0
15 0
15
-2
10 -2
10 -2
10
5
5
5
-4 0.04
0.06
0.08
α
0
-4 0.04
0.06
0.08
α
0
-4 0.04
0.06
0.08
α
0
Figure S-5. Heatmaps which show the number of TPs that take a particular 𝑇1 = 𝜏1 − 𝑡1 value for a given 𝛼 for parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c) and 𝛿𝛼𝛽 = 0.01. The green vertical line in each panel corresponds to the optimal 𝛼 for each parameter setting. Across Figs. S-5 (a)-(c) we observe the following common trends: 𝑇1 < 0 for most TPs when 𝛼 ≲ 0.055, meaning that our algorithm typically detects CTs from the NS to S state at times before the expert’s annotation of the corresponding seizure onset. 𝑇1 ≈ 0 for most TPs when 0.055 ≲ 𝛼 ≲ 0.085, thus showing close agreement between our algorithm and the expert annotations, with 𝜏1 and 𝑡1 differing by 0.5s in most cases. 𝑇1 > 0 for most TPs when 𝛼 ≳ 0.085, indicating that our algorithm typically detects most CTs from the NS to S state at times after the expert’s annotation of the corresponding seizure onset. S3 (b): Comparing residence times
Figure S-6 compares the probability densities of residence times in the (a) S and (b) NS states obtained ( ) from(our algorithm ) ( with ) (𝐸) those derived from the expert annotations. Specifically, we compare the lengths of the different 𝐼S(𝐸) and 𝐼NS to 𝐼S(𝐴) 𝑖 𝑖 𝑗 ( ) (𝐴) and 𝐼NS for the T2M time series studied in [1]. The residence times derived from the experts annotations are shown in 𝑗
black and the algorithm-derived residence times are computed from the 𝑡1 and 𝑡2 values obtained when using the optimal 𝛼 values, 𝛼 ∗ , for parameter settings 𝑃1 (in blue), 𝑃2 (in orange), and 𝑃3 (in green). Overall, the residence times obtained from our algorithm closely match those obtained from the expert annotations for both the S and NS states across all three parameter settings. S4: Analysis of PV curves (complimentary to analysis of ROC curves) Taking further inspiration from Temko et al. [2] and Mathieson et al. [3], we also consider the following metrics to examine the agreement between our algorithm and the experts annotations: Positive predictive value (PPV): 𝑁TP ∕𝑁S(𝐴) ∈ [0, 1], also referred to as the ‘seizure detection rate’, defines the proportion of 22/23
probability density [arb. units]
(a)
10
E
−1
P1
P2
P3
101 102 residence times in NS state
103
(b)
10−2 10−3 10−4 10−5
100
101 102 residence times in S state
103
100
Figure S-6. Probability density of residence times in (a) the S state and (b) the NS state based on seizure and non-seizure time intervals according to the expert (denoted by E) and the algorithm in the parameter settings 𝑃1 , 𝑃2 , and 𝑃3 for 𝛿𝛼𝛽 = 0.01, computed for the optimal 𝛼 in parameter setting.
1.0 (a)
(b)
(c)
PPV
0.8 0.6 0.4 : α∗ = 0.065 AUC = 0.8576
0.2 0.0 0.0
0.2
0.4 0.6 NPV
0.8
1.0 0.0
: α∗ = 0.085 AUC = 0.9528 0.2
0.4 0.6 NPV
0.8
1.0 0.0
: α∗ = 0.08 AUC = 0.9797 0.2
0.4 0.6 NPV
0.8
1.0
Figure S-7. PV curves obtained when using parameter setting 𝑃1 in (a), 𝑃2 in (b), and 𝑃3 in (c). In each panel, (NPV, PPV) points plotted in red correspond to different values of 𝛼, the green point corresponds to 𝛼 ∗ (the optimal 𝛼), the lower right corner specifies 𝛼 ∗ and the AUC. correctly classified seizure time intervals. (𝐴) Negative predictive value (NPV): 𝑁FP ∕𝑁NS ∈ [0, 1], the proportion of incorrectly classified non-seizure time intervals. Note, the NPV is a modified version of the ‘false detections per hour’ (FD/h) metric considered by Temko et al. and Mathieson et al. which is used to quantify the number of false positives that occur per hour in a given time series. The FD/h metric is an insightful clinical metric when the time series that are analysed are several hours long. Since the time series we analyse are much shorter, and are at most 1-2 hours long, we use the NPV to quantify similar behaviour. We use the PPV and NPV metrics to construct what we call ‘predictive value’ (PV) curves using the same procedure as ROC curves specified in [1] for increasing values of NPV. Similar to the TPR and FPR, the closer the PPV is to 1 and the NPV is to 0 the better the performance of our algorithm. In Fig. S-7 we show PV curves for the same parameter settings as in Fig. 3 in [1]. Each PV curve provides similar insight to its ROC counterpart.
References [1]. Flynn, A. et al. “Detecting seizure onset and offset times using human intelligence: A critical transitions based approach”. [2]. Temko, A. et al. “Performance assessment for EEG-based neonatal seizure detectors.” Clin. Neurophysiol. 122, 474–482 (2011). [3]. Mathieson, S. R. et al. “Validation of an automated seizure detection algorithm for term neonates.” Clin. Neurophysiol. 127, 156–168 (2016).
23/23