Epidynamics characterize and navigate the map of seizure dynamics

30 Seizures are a disruption of normal brain activity present across a vast range of species, diseases, 31 and conditions. Here we introduce a mathematical theory, Epidynamics, which provides a 32 conceptual framework to characterize how seizures start, evolve, and terminate. We provide the 33 first objective taxonomy of seizures based on a straightforward analysis of EEG data. Analyzing 34 over 2000 focal-onset seizures recorded from 7 epilepsy centers on five continents, we find 35 evidence of the predicted 16 Dynamic Classes of seizures. The theory also enables drawing a 36 map of brain dynamics that includes most seizure classes and status epilepticus. We demonstrate 37 that patients navigate the map during their lifetime and verify critical predictions of the theory. 38 Epidynamics not only provides a way to stratify patients in complement to present practical 39 classifications but also guides biophysically based mechanistic approaches and provides a 40 language to describe the most critical features of seizure dynamics. 41 42 Impact statement: 43 Epidynamics, a mathematically-derived method of classifying seizure dynamics, provides a 44 rigorous method for classifying and quantifying seizures and a theoretical framework for 45 understanding seizure initiation and propagation. 46


INTRODUCTION
single group for seizure onset, and SH without DC shift and SNIC as a single group for offsets 198 because there is no method to distinguish them. We then compared these same methods on 120 199 human seizures. We found that concordance was also reliable in human data (Appendix I.5). 200 These results show that the chosen features are capable of distinguishing the different 201 bifurcations reliably for both human visualization and algorithms, that human seizure dynamics 202 are consistent with the modeled bifurcations, and that human reviewers can use the methods 203 described in the next paragraphs to classify onset and offset bifurcations reliably. In addition, we 204 found that human reviewers were more reliable than the automated algorithm in noisy clinical 205 data. We then used the human markings on the clinical data for the taxonomy in Fig. 2. 206 207 Onset dynamics -To analyze onset dynamics, we first investigated all 51 patients that had been 208 recorded using equipment capable of visualizing DC shifts. Of note, DC shift alone has 209 previously been shown to be highly correlated to the seizure onset zone (Ikeda et al., 1996;Ikeda 210 et al., 1999;Kanazawa et al., 2015). Many seizures (41%) started with constant amplitude spikes 211 and a DC shift, signifying SN bifurcation (Fig. 2E). The second most frequent was similar 212 amplitude/frequency dynamics without a DC shift (either SubH or SN, 37%), followed by SupH 213 (14%) and SNIC (6%). Ambiguities in the classification were treated systematically as detailed in 214 the Appendix (I.6) , and only 1/51 seizures could not be agreed upon by the reviewers. Despite 215 the inability to distinguish SN (-DC) and SubH, we demonstrate that at least three of the 216 predicted bifurcation types are present at onset in human seizures. Importantly, some seizures 217 display complex dynamics because they go through more than one bifurcation as the seizure 218 begins. For example, nearly half of the seizures labeled as SN onset progressed into square-root 219 amplitude scaling (likely due to spread of the seizure to neighboring tissue) after 2-5 seconds, 220 consistent with a switch to SupH dynamics; we only labeled the initial bifurcation herein. For this 221 reason, we have slightly modified the definition of onset bifurcation as compared to the literature. 222 Here the onset bifurcation is the one causing a departure from the resting state, even if it is not 223 directly causing the onset of oscillations as there can be intermediate states (Appendix II.1). 224 225 We next looked at the onset bifurcations in the other 69 seizures that were recorded with non-DC 226 coupled hardware. Of the four onset bifurcations, SN and SubH become indistinguishable 227 without a DC shift available. Combining these two into a single group, we found that 67% were 228 SN-SubH, 26% were SupH, 6% were SNIC, and only one seizure did not have reviewer 229 consensus. 230 231 Offset dynamics -We first examined offsets in those patients with DC recordings for the most 232 robust classification. The analysis of the interspike intervals (ISI) and spike amplitudes revealed 233 a logarithmic/square-root slowing-down with constant amplitude in 20/51 patients, and of those 234 10 had DC shifts at offset. Thus, 10/51 (20%) were SH, while the remaining 10 (20%) were 235 potentially SH or SNIC (Fig. 2F). The remaining 31 patients did not have slowing at the end of 236 their seizure (53% FLC, 6% SupH, 2% no reviewer consensus). Analyzing arbitrary dynamics. One significant challenge in the above analyses was the 245 presence of noise. Within the theoretical model, "arbitrary" dynamics refer to abrupt changes 246 without clear scaling laws to or from zero. However, in human data analysis "arbitrary" includes 247 a wide range of other behaviors, especially noise. It is important to note that the taxonomy above 248 includes the first seizure that could be analyzed from every patient-we did not restrict the 249 analysis to "clean" seizures. Some of the seizures were noisy, either from technical concerns or 250 physiological effects of the seizure. We chose this method in order to provide a robust, real-world 251 demonstration of this analysis. Because noise can be classified as "arbitrary," this analysis may 252 overestimate the numbers of SubH-SN onsets (no DC shift) and FLC offsets. However, it was 253 clear that this limitation was not merely technical: several patients had complex physiological 254 dynamics. Four of the FLC seizures were highly unusual from a clinical perspective: one had 255 increasing amplitude for the last 10 seconds, one consisted of low-voltage fast activity that ended 256 abruptly without any other change, one ended with irregular spike waves, and one had 257 accelerating frequency at the end (Appendix I.8 and V). Correlation between clinical data and seizure type 283 We compared all available clinical metadata from patients with their Dynamic Class and found 284 no correlation between seizure class and patient gender, pathology, or localization. There was a 285 correlation with age, as older patients tended to have more SupH onsets (Appendix I.9). We also 286 compared these results with a prior visual classification that identifies 7 basic seizure onset 287 patterns (Perucca et al., 2014), and found 6/7 patterns without any apparent relationship to 288 clinical data or pathology (Appendix I.10). There were no significant similarities between the 289 Dynamic Class and the visual classification. 290 291 Seizure dynamics vary in human epilepsy 292 While analyzing this dataset, we noted that one patient had two consecutive seizures belonging to 293 different classes: one supH/supH and one supH/SH, raising the possibility that an individual may 294 express different types of seizures. This finding was surprising, as many clinicians assume that a 295 person's seizure should be "stereotyped", i.e. consistent over time. In  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint 25 seconds) to measure the offset ISI. To be conservative, we only classified seizures as slowing-305 down or constant if such a determination was unequivocal and labeled the rest as "not assessed," 306 meaning that the determination was not readily evident on brief visual inspection. As seen in Fig.  307 3, all 13 patients expressed at least two offset patterns. Note that this is likely an underestimation 308 of the heterogeneity in seizure classes: these recordings did not contain DC coupling so onset 309 bifurcations were not assessed, and we did not distinguish the offsets into all four types. recordings. This greatly reduces the degrees of freedom necessary to describe the observed 320 activity, i.e. a small number of differential equations are sufficient to describe the collective 321 behavior (Fig. 4A). We here consider a system with the minimum number of variables necessary 322 to produce oscillatory activity, two. Based on the parameter values, two states can be 323 distinguished: resting (fixed point) or oscillatory (limit cycle). When these two states coexist for ictal duration. We previously validated this approach with the "Epileptor," a set of five 330 differential equations able to account for the most dominant seizure class (Jirsa et al., 2014) 331 (SN/SH, also known as 'square-wave' bursting (Rinzel, 1987)) . In the Epileptor, the transition 332 from "normal" to seizure state and back again, as well as the seizure dynamics, are controlled by 333 a collective permittivity variable that evolves on a slow time scale. However, the Epileptor only 334 accounted for a single Dynamic Class and is not sufficient to explain the data presented above. 335 The fact that individuals can express different classes of seizures over time leads to two 336 predictions that must be included within a model of human seizure dynamics: different classes of 337 seizures must coexist in the same model, and there must be an endogenous mechanism by which 338 the brain can transition slowly between the different classes. can be placed (Dumortier et al., 1991). In effect, the map is a representation of the range of states 345 in which a brain region can exist, oscillatory (ictal) and non-oscillatory (interictal), and the 346 transitions between them. The oscillatory state produces spiking activity that is described by the 347 fast variables (millisecond scale activity). However, on a much slower time scale of the order of 348 seizure length, the brain can move towards a transition to an interictal state, as described by a 349 slow variable. Migrations to different locations on the map can occur on a usually even slower 350 timescale (10's-1000's of seconds), which we call here "ultraslow". The use of an ultraslow 351 variable allows full exploration of the map (Saggio et al., 2017). Applying these general 352 mathematical principles to epilepsy implies that any brain region able to generate SN/SH seizures 353 can potentially generate other classes by navigating to different dynamical regimes (i.e. changing 354 the parameters of differential equations) (Saggio et al., 2017). The map related to the 355 organization of seizure classes is shown in Figure 4B-C. The state of a brain region at any 356 moment can be represented as a location on the map, which will define its dynamical properties. 357 Regions in the map that correspond to different regimes, including resting and ictal states, are 358 separated by bifurcation curves. Seizures are represented as black arrows, each arrow 359 corresponding to one specific seizure class. To produce a seizure, the system, which is initially in 360 the resting state within the bistability region, heads towards the onset bifurcation curve. When 361 this curve is reached, the resting state disappears and the system is forced to go into the 362 oscillatory seizure regime within the bistability region. This transition in state causes an inversion 363 in the trajectory of brain state, with the system now heading towards the offset bifurcation curve. 364 When the offset is reached, the system goes back to rest and inverts direction again. Movement 365 along the black arrow is produced by slow (of the order of the ictal length) mechanisms leading 366 to seizure offset. Note that the movement towards the onset and offset bifurcations occurs in both 367 cases from within the bistability region. Ultraslow movements on the map are responsible for 368 changing the location of the brain state while at rest and enable the expression of different classes 369 of seizures as observed clinically. This theory thus provides a potential explanation for the 370 clinical observation of multiple classes of seizures in a single patient, and the map provides a 371 hypothesis to describe how a patient's current state (i.e. location on the map) can affect seizure 372 dynamics (whether a seizure is likely to occur, and what type is most likely). Figure 4C  It is important to note that the topology of the map in Fig. 4 was initially proven mathematically 383 to be generic and rigorously valid for bursting (Dumortier et al., 1991;Baer et al., 2006). This 384 invariance establishes the ground truth to define the relationships between the different 385 bifurcations in the proximity of the SN/SH class, which leads to a key prediction: transitions 386 between certain classes may be more common due to their proximity on the map. . We analyzed the corresponding trajectories on the map to 418 determine how this had occurred. Prolonged seizures occurred when the brain state crossed the 419 SN onset curve but was unable to return to rest through the offset bifurcation and remained 420 mainly in the violet "seizure only" region in Fig. 5G. The slow variable naturally drives the state 421 towards offset, but in these cases was continually overridden by noise, causing the state to 422 'escape' from the bistability region. We then analyzed this effect by simulating various levels of 423 noise and showed that there was a clear correlation between the noise variance and the likelihood 424 of entering status epilepticus (Appendix IV). We compared these results with our clinical data, 425 which had two examples of non-convulsive status epilepticus. In both cases, the seizures began in 426 typical fashion, but instead of terminating began having long periods of constant ISI with varying 427 periods of amplitude fluctuations. There were many abrupt transition periods during which the 428 ISI and amplitudes became arbitrary. After these brief periods of disorganization, the dynamics 429 returned to constant ISI. We compared the dynamics of our model results with these human 430 seizures and found that the transition between different dynamics is quite similar. In Fig. 5D, we 431 show a portion of two human seizures and one example of the simulation (Fig. 5E-F)  Accelerating seizure: We then analyzed the seizure offset that increased in frequency described in 440 Appendix I.8. In this case, we explored conditions on the map that could produce "speeding up" 441 at the end of the seizure. We identified multiple trajectories in the map of brain states capable of 442 producing these unusual seizure dynamics. As in the case of status epilepticus, these unusual 443 patterns are dependent upon the relative position within the brain map, in this case occurring 444 when brain states along a trajectory are affected by multiple bifurcations that are in close 445 proximity (Appendix V). These results further validate our theory of Dynamic Class and show 446 how it provides a rational explanation for a wide range of physiological dynamics. 447

449
Seizures have been recognized clinically for millennia, but after nearly a century of 450 electrographic recordings we still do not have a scientifically-based method of characterizing 451 them. We here address this issue and provide the first general mathematical theory that explains 452 the basic features of seizure onset and offset dynamics, which we term collectively the 453 'Epidynamics' of seizures. In canonical terms, there are 16 Dynamic Classes of seizures, and our 454 theory provides a general dynamical map that describes how a brain region can move between 455 these conditions. Our classification aims at precisely identifying the seizure type in terms of 456 dynamics, without any dependence upon specific symptoms, pathology, or localization. Thus which clinicians report the frequency and morphology of spikes. This method is helpful to 488 identify primary generalized epilepsies, but within focal seizures has limited clinical use. Two 489 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; II.5). In this work with focal epilepsies, the most common Classes were the SN (with DC 521 shift)/FLC and SN/SH (which are also likely Classes). 522 523 When comparing our results with a past visual classification system of spike frequency,(Perucca 524 et al., 2014) we found no correlation with pathology in our cohort of 120 patients. However, our 525 cohort did not have any patients with tuberous sclerosis, which was the only pathology associated 526 with burst suppression in that prior work. Combining the data from both studies, the different 527 patterns appear to be either evenly distributed or too rare to find robust correlations with 528 pathology. Similarly, Epidynamics is not strongly correlated with pathology. In terms of the map 529 of brain activity, we hypothesize that what determines the seizure dynamics is not the pathology 530 per se, but the location of the brain on the map. Specific pathologies may predispose to certain 531 regions, but there are many complex dynamics affecting brain state and many conditions that can 532 produce similar dynamics. This coincides with the idea of this brain map showing the full range 533 of potential seizure onset and offset activity. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Those results and ours strongly support the proposal that patients move closer and farther away 580 from seizure threshold (i.e. "travel the map") during their lifetime. This theory may also be 581 helpful in assessing a brain's current proximity to seizure bifurcations, i.e. predict the risk of 582 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. Our theory describes a map of brain states that accounts for the potential Dynamic Classes of 594 seizures, and how the resting brain can navigate between them. This is not a unifying theory to 595 describe seizures completely, as we defend a complementary approach, combining it with 596 traditional operational classifications. However, it provides a unique avenue to classify seizures 597 based upon their key dynamical features, while providing insight into how seizures become more 598 or less likely to occur at a given time. A necessary corollary is that, although our theory has been 599 developed in the context of seizure dynamics, it likely extends to physiological function of the 600 healthy brain (e.g. alterations between REM and slow wave sleep, and the appearance of gamma 601 frequencies or ripples during slow wave sleep) and stipulates the existence of at least two time 602 scales in any theory of the brain. Slow time scales are present in many theories of brain function, 603 but typically have been limited to the domain of learning and adaptation, thus functionally 604 separated from fast processes. Here the functional integration and co-evolution of the fast 605 neuroelectric and slow permittive time scales suggests emergent and inherent properties of brain 606 processes. 607 608 609 610 611 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Acknowledgements: 612 We thank Masao Matsuhashi for efforts in data acquisition. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint 651 652 653 654 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint between healthy and ictal states (yellow), coexistence between healthy and 'active' rest (a non-695 oscillatory state with a different baseline than the healthy state (gray)). When a seizure starts 696 because the system crosses an onset bifurcation, the slow variable enables movement along the 697 black arrow (a path in the map) to bring the system back to rest. Note that in this model, the 698 system alternates between the resting and seizing states within the bistability regions. The shape 699 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint of the arrow is meant to better show the trajectory followed by the system, however movement in 700 the model occurs back and forth along the same curve. Insets show example of timeseries for 701 different paths. SubSH: Subcritical Saddle Homoclinic, an unstable limit cycle that occupies a 702 small portion of the map but is incapable of starting/stopping seizures, and thus is not included in 703 of the modelled seizures in C,F. The seizure from C begins (red triangle) and naturally moves 723 upwards towards the supH bifurcation, but a change in the ultraslow drift prior to termination 724 pushes the system downward, changing the offset from SupH to SH. Conversely, the seizure 725 from F begins (black triangle) and repeatedly moves towards the SH termination due to the 726 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. is defined as the time between sequential spikes. Amplitude for a given spike is defined as the 981 absolute maximum peak-to-peak difference in a window of time ranging from the halfway point 982 of the first spike and the halfway point of the latter spike (see Fig. 2A in main body). 983 Beginning/ending spikes utilize only the halfway point of the nearest spike. All analyzed data 984 were decimated to ~200 Hz for efficient analysis. In order to remove slow transients and identify 985 the local spike amplitudes, raw EEG data were first high-pass filtered (MATLAB's 'highpass', 1 986 Hz). Peaks were found in Matlab using the 'findpeaks' function, after manually optimizing 987 amplitude, time, and prominence for each patient. This process required iterating the features 988 until visually confirming that spikes were correctly detected, using plots in Matlab with the spike 989 detections superimposed upon the EEG signal (see Appendix - Fig. 1). Due to the wide 990 variability between patients, we were unable to develop a reliable automated method to determine 991 the spike locations in every patient. 992 Sentinel spikes. Some seizures began with a single high amplitude 'sentinel spike,' which 993 typically was present in many electrodes and was not directly related to the seizure dynamics 994 thereafter. Additionally, it was sometimes difficult to know whether these events were actually 995 true sentinel spikes versus a transient caused by a filtered DC shift (e.g. the first spikes seen in 996 Fig. 2B&C, Appendix -Fig. 1A&C). This spike was used as the initial seizure onset time but was 997 not included as a detected spike for seizure dynamics, as this phenomenon is not a canonical 998 feature of any of the onset bifurcations. 999 Spike-wave complexes. Some seizures had spike-wave complexes in which the aftergoing 1000 slow wave was very prominent, similar to the waveforms seen in absence epilepsy. In these 1001 cases, the highest amplitude of the whole complex was used, even if it was the slow wave, and 1002 only a single event was counted from each complex (see Appendix - Fig. 1C). 1003 Clonic spiking. As described previously (Jirsa et al., 2014), when there was clonic 1004 spiking at the end seizures, the interclonic interval was used to evaluate the presence of slowing 1005 down, and the ISI between successive fast spikes within each clonic burst were ignored (see Fig.  1006 2B & Appendix - Fig. 1F). 1007 It is important to note that this analysis does not address every feature of the seizures. 1008 There are many phenomena, such as the three listed above, the shape and frequency of the spikes, 1009 and other complex patterns, that are not invariant features of the bifurcations. While these are 1010 important to the clinical description and can be relevant to understand the underlying dynamics, 1011 in the framework here proposed they do not contribute to defining the Epidynamics of seizure 1012 onset and/or offset. 1013 1014

I.3 Visual classification of seizure dynamics 1015
The key to differentiating the bifurcation type is to identify the invariant dynamical 1016 features, which can be summarized as the presence of a DC shift and the behavior of the ISI and 1017 amplitude (see Fig. 1). These features are typically quite easy to distinguish. The only prominent 1018 ambiguity is that it is not feasible under clinical conditions to distinguish between the logarithmic 1019 and square root functions at offset, as previously described (Jirsa et al., 2014). Thus, our first test 1020 was to determine if human reviewers can classify the different bifurcations visually using simple 1021 rules. 1022 The presence of a DC shift and the general trends of ISI and amplitude can readily be 1023 determined upon visual analysis. The basis for this analysis is to determine whether the amplitude 1024 and ISI scale to zero. For ISI, this appears as a decreasing frequency as T approaches 0, i.e. 1025 slowing down at seizure offset, or speeding up after seizure onset. For amplitude, it appears as a 1026 gradual change in the spike amplitudes, with the spike at T=0 being very small compared to the 1027 baseline, then increasing further away from T=0. This description is qualitative but is readily 1028 applied to typical EEG data. Since the onset/offset dynamics are typically defined by only 5-10 1029 spikes and there is considerable noise and variability in real EEG signals, rigorous curve fitting is 1030 rarely possible (though we included it in the examples in Appendix - Fig. 1). Just as with clinical 1031 EEG reading, we found that a much simpler and more reliable classification system was to 1032 visualize the ISI and amplitude plots of the first and last 10 spikes and determine whether the 1033 trends were scaling to zero, constant, or arbitrary, and if there was a DC shift. For amplitude, we 1034 defined 'scaling to zero' as steadily diminishing to less than 3 times the background level near 1035 T=0. For ISI, we define it as steadily larger ISI near T=0, with the last two ISI > 50% larger than 1036 the mean ISI 10 seconds prior. All analysis for onset and offset concentrated on the first/last 5 1037 seconds of data, but occasionally used up to 15 seconds to observe the full patterns. Using these 1038 definitions, we developed a visual classification system (Appendix -Table 1). 1039 1040 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020.

I.4 Features for automated classification of seizure dynamics 1051
We also developed an analytical tool to use quantitative features and machine learning to 1052 identify the Dynamic Class. The goal of this analysis was to determine if the features used in the 1053 qualitative study were robust. We designed features based upon the visual classification system in 1054 Appendix - Table 1, focused on quantifying DC shift, amplitude trends, and interspike interval  1055 trends. Definitions, feature computation, and feature descriptions are as follows: 1056 1057 Baseline Definition: Several features require definition of the baseline, i.e. seizure activity vs. 1058 non-seizure activity. All analyzed waveforms include a period of baseline, followed by seizure 1059 activity, followed by more baseline. Computationally, the term 'baseline' below is defined here 1060 as the segment of the waveform that started before/after a seizure. For an onset baseline, this 1061 segment was taken as the start of the waveform up until seizure onset. For an offset baseline, this 1062 segment was taken as the point after seizure offset to the end of the waveform. 1063 DC vs. Non-DC: Data acquired at Kyoto (Nihon Kohden EEG 1100 amplifier) or Michigan 1064 (Natus Quantum amplifier), both of which record down to 0.016 Hz, were included in the DC 1065 cohort, and all other data was considered "non-DC". Non-DC data were first filtered with 1 Hz 1066 highpass filter (Matlab 'highpass'), then all features extracted. On the DC data, features were 1067 computed on the raw, unprocessed EEG. 1068 Preprocessing steps: 1) Onset/offset times, as well as the relative bifurcation window lengths, 1069 were determined by a trained epilepsy specialist. 2) Spikes were identified using findpeaks.m 1070 (Matlab) to locate upper (maxima) and lower (minima) spikes. 3) Seizure polarity was 1071 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint determined by determining whether the median of the upper or lower spikes had a greater 1072 absolute difference from the baseline median. The value with larger difference was chosen as the 1073 true "spikes" and were used for all further analysis. All amplitudes were taken as absolute 1074 values. This step accounts for the fact that spikes can be either positive or negative in intracranial 1075 electrodes. 1076 Feature 1 -ISI trend: The ISI of all spikes were computed and plotted consecutively. A simple 1077 line was fit to the data with a least squares algorithm and the slope of the line fit was extracted as 1078 the overall ISI trend. For onset (offset), the first (last) 5-15 spikes were used, using as many as 1079 possible until a clear inflection point in the line. The order was reversed for offset. 1080 Feature 2 -Amplitude trend: This feature was computed exactly as in ISI trend, except the peak-1081 to-peak spike amplitudes were used in place of the ISI. 1082 Features 3&4normalized upper and lower peak median: The signed distance between the 1083 median of the upper peaks and the baseline median was computed (upm). That same was done for 1084 the lower peaks (lpm). If |upm| > |lpm|, then normalized values were nupm = 1, nlpm = upm/lpm. 1085 If |lpm| > |upm|, then nlpm = -1, and nupm = -lpm / upm (note that lpm is negative in all cases). 1086 These features identified DC shifts. 1087 1088

I.5 Validation of classification methods 1089
In our prior work with the Epileptor model (the SN/SH Class), we proved the goodness-1090 of-fit to the logarithmic equation for terminal ISI (Jirsa et al., 2014). However, in the current 1091 work we found that GoF methods were not robust when discerning between multiple seizure 1092 types in noisy data. This is because the number of samples (spikes) is small, while the combined 1093 uncertainty in the noise as well as the classification of the different Classes can be significant. 1094 The result was that the variability due to noise was often more prominent than the differences 1095 between bifurcation types. This difficulty is not surprising, as seizures are notoriously difficult 1096 for automated algorithms to identify under clinical (noisy) conditions, and classification will be 1097 even more difficult. More importantly, even under ideal conditions any such analysis would be 1098 limited due to the lack of a gold standard. We therefore decided on an alternate and more 1099 rigorous approach, appropriate for the data quality of clinical seizures. 1100 In order to validate our classification system, we used a multi-step approach. The key to 1101 this validation was the generation of a gold-standard with simulated data, in which we know 1102 which bifurcation is present. Using the gold standard, we first used machine learning tools to 1103 assess how accurately our chosen features are capable of identifying the differences between the 1104 bifurcations. We then tested how accurately human reviewers could identify each bifurcation 1105 using visual review. After proving that the human reviewers were accurate and reliable in the 1106 gold standard, we tested their reliability in the true clinical dataset. These steps are detailed 1107 below. 1108 1109 Step 1: Feature validation and bootstrapping: 1110 OBJECTIVE: Determine whether the chosen features capable of discerning between the 1111 four bifurcations in the gold standard. 1112 METHODS: Using the simulated seizures (240 onsets, 240 offsets), supervised learning 1113 was chosen to quantify how well the simple features in I.4 captured the differences between the 1114 different Dynamic Classes. This was performed by fitting labeled features to a statistical model to 1115 generate a goodness-of-fit measurement. Specifically, we used a multinomial logistic regression 1116 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; model ('mnrfit', Matlab), which outputs a model parameter 'deviance' that estimates the 1117 Goodness-of-Fit. The GoF was compared with the gold standard in the simulation data. Note that 1118 this method is not used for classification nor to determine the "best" potential fits. Rather, we 1119 used the GoF statistic to assess whether our chosen features were capable of discerning between 1120 the four bifurcations. While there certainly could be better features, we chose these because they 1121 recapitulate what experts use for visual analysis, which is what this test is trying to validate. 1122 To provide reference for the GoF computation, we performed a permutation test by 1123 randomly scrambling the bifurcation labels and re-fitted the data to the model. This was repeated 1124 10,000 times to provide the distribution of random GoF, which can then be compared with the 1125 index case and provide a p-value, a process known as bootstrapping. 1126 RESULTS: For both onset and offset bifurcations, the GoF with true labels was better 1127 than all 10,000 permutations (p < 1e-4). This result indicates that the 4 features are very effective 1128 in capturing the differences between the 4 onset and offset bifurcations and extremely unlikely to 1129 be due to chance. 1130 CONCLUSION: The chosen features are capable of discerning between the bifurcations. 1131 1132 Step (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ;

DC (N=51)
Non-DC (N=69) Onset P = 1.78e-15 P = 4.51e-12 Offset P = 4.51e-11 P = 9.72e-4 DC Non-DC Onset P < 1e-4 P = 0.0969 Offset P = 7e-4 P < 1e-4 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint confused the reviewers. However, this situation most likely is not pertinent to 1242 physiological seizures, as it is an artifact of the model construction. 1243 The following phenomena are not described by the canonical features of the bifurcations, and 1244 were difficult for the reviewers to classify. The majority of disagreements between reviewers 1245 contained at least one of the following: 1246 c. Polyspikes or spike waves at seizure onset, especially prior to low voltage fast activity 1247 d. Preictal spikes 1248 e. Noisy data, especially 'extraneous' spikes that disrupt scaling laws and "Unusual 1249 dynamics" (#5 in preceding section) 1250 f. Ambiguous seizure onset/offset times 1251 g. Low voltage fast activity that speeds up at seizure onset 1252 These patterns (a-g) are inherent challenges of seizure data, and such conditions have the 1253 potential for ambiguous classification. However, we do not consider this to be a failure of our 1254 method-rather, Epidynamics provides a method of identifying and quantifying these conditions. 1255 While these phenomena are not uncommon, they are not invariant features that define the 1256 bifurcation. In effect, the fact that they deviate from the basic Epidynamics not only allows us to 1257 identify them as atypical, but also to quantify how and why they are atypical. 1258 We include these graphical fits and RMSE for constant and square root/logarithmic equations to 1259 demonstrate how this tool may be used to determine parameters for these equations, which can be 1260 used for quantitative analysis, comparison, and robust discussion of seizure dynamics in future 1261 work. 1262 1263 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Appendix - Fig. 1A (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Appendix - Fig. 1B (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ;

1284
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Appendix - Fig. 1C (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Appendix - Fig. 1D (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint Appendix - Fig. 1E: SN

I.10 Comparison with previous classification method 1398
We also compared our results with a previous classification method (Appendix -Table 5). 1399 The low-frequency high-amplitude periodic spikes, (iii) sharp activity, (iv) spike-and-wave activity, 1401 (v) burst of high-amplitude polyspikes, (vi) burst suppression, and (vii) delta brush. Additionally, 1402 we created the category (U) to indicate seizure onsets that were not explained by the Perucca 1403 classification. In our analysis, we did not find any seizure onsets that corresponded with (vi). 1404 There was no significant correlation between the onset types and any of the patient metadata (chi-1405 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; square). In addition, we compared these classifications with the onset and offset bifurcations. A 1406 chi-squared test was performed on all combinations of Jirsa onset (i.e. based upon bifurcations in 1407 (Jirsa et al., 2014) ), Jirsa offset, and Perucca onset. There was no significance between any of the 1408 pairings. Note that two patients had type (vii) (delta brush, see Appendix - Fig. 1E In this work we used the model in (Saggio et al., 2017) and we refer to that paper for 1427 details. It consists of three ordinary differential equations: 1428 1429 Equation 1 Takens-Bogdanov singularity (Dumortier et al., 1991). They depend upon three parameters 1434 ( 2 , 1 , ). We can consider a sphere, with radius R, centered at the origin in the three-1435 dimensional parameter space. On the spherical surface there are curves of bifurcations that divide 1436 the surface into regions with different sets of attractors. It can be shown that, up to a certain value 1437 of the radius, the bifurcation diagrams on the spherical surface do not change (Dumortier et al., 1438(Dumortier et al., 1991. For this reason, the number of parameters relevant to describe the bifurcations in the 1439 system can be reduced, from the three Cartesian parameters ( 2 , 1 , ) to the two spherical 1440 parameters (, ) keeping the radius fixed. In the present work we use a radius R=0.4. This 1441 reduction allows for easier analysis and visualization of the parameter space. A flat sketch of the 1442 spherical surface and its bifurcations can be found in the map in Fig. 4C. A flat projection of the 1443 real map, as shown in Fig. 5, is obtained with Lambert equal area azimuthal projection. Details 1444 on how to reconstruct these maps can be found in (Saggio et al., 2017) . Movement on the sphere 1445 is promoted by making the three Cartesian parameters (or the two spherical ones) depending on a 1446 third variable, z, acting on a slower timescale described by the parameter 0 < ≪ 1. Movement 1447 is implemented so that, when the fast subsystem ( , ) is in the resting state ( , 0), it moves 1448 towards the onset bifurcation, and when the distance from the resting state is bigger than * it 1449 moves towards the offset bifurcation. When * >0 this gives periodic bursting. 1450 1451 The shape of the path along which the fast subsystem moves is taken to be the arc of the 1452 great circle linking the offset point to the onset point on the sphere. It is described by the 1453 following parameterization: 1454 1455 Equation 2: 1456 = ( 2 − 1 )

1457
( ) = ( cos + sin ) 1458 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint where = / and = (( × ) × )/‖( × ) × ‖. An example of path for each class is ability to identify a seizure trajectory prior to crossing the oscillatory bifurcation has intriguing 1500 implications for future work. 1501 1502 In Appendix - Fig. 4 we show bifurcation diagrams for the classes in the bistability region 1503 in the lower portion of the map. Here classes occur without baseline shifts, even when starting 1504 with a SN bifurcation. This makes it impossible to discern SN from SubH onsets based on Fig. 1.  1505  1506  1507 1508 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. ; https://doi.org/10.1101/2020.02.08.940072 doi: bioRxiv preprint III. Switching between Dynamic Classes due to fluctuations and ultra-slow 1675 modulations 1676 In the model above, the path along which the fast subsystem slowly moves is a simple arc 1677 linking the offset and onset points. However, movements promoted by real changes of parameters 1678 can be more complex. As a proof of concept of the effects that this can have on the system, we 1679 considered modifications of the path due to (i) ultra-slow drifting of the offset and onset points 1680 and and (ii) fluctuations produced by noise. 1681 The ultra-slow drift was obtained as in ( (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020. (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted February 9, 2020.  CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made