A taxonomy of seizure dynamotypes

1Epilepsy is one of the most common neurological disorders with an estimated prevalence of 50 million worldwide. It is characterized by spontaneously recurring seizures, which are ‘a transient occurrence of signs and/or symptoms due to abnormal excessive or synchronous neuronal activity in the brain’. However, there are a vast array of signs, symptoms, and underlying causes of seizures. Thus, despite high prevalence and considerable morbidity and mortality, it has been challenging to characterize, treat, and understand seizures, which prevents the development of reasoned, mechanistic approaches to therapy and improved patient care. Seizure classifications to date have been purely descriptive of empirical data: clinical manifestations (e.g. focal vs. generalized) that are based upon the region of brain affected rather than the seizure itself, and visual descriptions of electroencephalogram (EEG) waveforms. These classifications have been subjected to numerous revisions. In its latest position paper, the International League Against Epilepsy states: ‘Because current knowledge is insufficient to form a scientifically based classification, the 2017 Classification is operational (practical)’. In effect, that classification is based upon the epilepsy phenotype—the clinical symptoms that arise during a seizure. Over the past several years, a separate approach has arisen: investigating the genotype of specific epilepsies, which may lead to more informed treatment decisions that match deficits with mechanisms. And for decades, intractable epilepsy has been treated with epilepsy surgery, which relies upon the high spatial resolution of imaging and implanted electrodes to find a seizure focus. These approaches are based upon our available tools: clinical expertise, genetics, imaging, pathology, and surgery. However, seizures are by definition dynamic phenomena, and none of these tools characterize the fundamental dynamics of seizures.

2In this work, we introduce an organizing principle of seizure dynamics based on nonlinear dynamics and bifurcation theory. Bifurcations are sudden qualitative changes in behavior, including onset and offset of oscillations. Here, we introduce the term ‘dynamotype’ to describe a seizure’s composite, observable, dynamic characteristics in electrophysiological recordings comprising seizure onset and offset. Together, dynamotype, phenotype and genotype provide a rich, multifaceted description of the dynamics, clinical manifestation, and underlying pathology of a seizure. The organization of seizures along dynamotypes leads naturally to a Taxonomy of Seizure Dynamics (TSD) providing practical, objective metrics for classification. As the periodic table of elements is a tabular display of chemical elements arranged according to proton number and electronic configuration, TSD is a tabular arrangement according to bifurcation type of seizure onset and offset. Furthermore, the organization of the periodic table can be used to derive relationships between the various element properties and predict chemical properties and behaviors of undiscovered or newly synthesized elements. Here, we explore the capacity of TSD to fulfill this second functional part of the analogy also and demonstrate the existence of all dynamotypes in human epilepsy. We then discuss TSD in the context of a canonical model in nonlinear dynamics and identify the relations amongst the seizure dynamotypes. TSD is available for immediate transfer to clinical practice, providing a rational method of characterizing seizures and subsequently a better understanding of the underlying principles governing seizure generation and termination.

3The basis of TSD is the observation that seizures are characterized by abrupt changes in the EEG waveform at seizure onset and offset, which we interpreted as bifurcations known from dynamic system theory. As a seizure evolves, the brain moves from a normal state into a seizure and back again. Recent work has focused on describing these transitions with empirically chosen visual patterns and has found interesting relationships with underlying pathology, surgical outcome, and sudden unexpected death in epilepsy. However, these transitions can also be described more rigorously and mathematically as bifurcations. Bifurcations represent qualitative changes that both define and constrain the system dynamics. The concept has been used to understand neuronal firing: when a neuron goes through a bifurcation, the emergent dynamics often comprise a new set of behaviors such as multiple stable fixed points (rest activity) or limit cycles (oscillatory activity). When a neuron oscillates quickly about a limit cycle, it produces fast repetitive activity known generically as bursting, which is a periodic wave form with periods of oscillatory spiking and periods of quiescence. These dynamics can be described by a set of differential equations and accompanying parameters. The concepts about periodic wave forms with alternating periods of oscillations and quiescence have been extended to EEG wave forms in seizures, where oscillatory (ictal) states and quiescent (non-ictal) states alternate. Two variables are the minimum necessary to generate oscillations, and in systems with two variables, there are only six types of bifurcations involved in bursting (Figure 1, for further details also see ). Four can be used to enter the bursting regime, and another four to exit, giving a total of 16 possible dynamotypes. A key benefit of this organizing principle is that it unambiguously identifies the invariant properties of individual events, which may provide mechanistic insight into the underlying causes and response to specific interventions. It also provides a model that not only accounts for the effects of noise on the system and multistability, but also generates a time series. Generalizing this to epilepsy, proposed the existence of 16 theoretically possible dynamotypes (i.e. seizure types), and found one seizure offset bifurcation that was present across multiple species, brain regions, and pathologies, including a small cohort of humans. Based on that initial work, we now expand and present a taxonomy of seizure dynamotypes. In the following, we begin with the definition of the different types of seizures based on dynamics at their onset and offset. Then, we show that seizures recorded from different centers in the world can be rigorously classified, and how classification can be performed in daily clinical practice. Next, we introduce a canonical model in nonlinear dynamics with two important properties: 1) the model is canonical, which means that under certain mild conditions (see Materials and methods) the behaviors of other models of arbitrary physiological detail can be represented and explained by the canonical model; 2) the canonical model captures all dynamotypes in a single mathematical representation. Transitions between types can be obtained through an ultra-slow modulation of the model’s parameters providing a map of the parameter space, which systematically predicts relations between dynamotypes, including a hierarchy across dynamotypes. We demonstrate from a large repertoire of empirical data that patients navigate this seizure map to express the different types of seizures. Finally, we discuss how TSD can be used in a wide range of novel applications in clinical care and research.

4The goal of this work is to characterize seizures by their underlying onset and offset dynamics, which depends upon identifying reliable, canonical dynamic features. While the dynamics of a single neuron have already been described, linking that behavior to a seizure generated by millions of neurons is complex. We chose to analyze the EEG signal from standard intracranial electrodes, as it is the most clinically relevant and widely studied method to measure brain dynamics. These patterns, which are visible within the EEG waveform, identify the bifurcations that define the invariant properties, and thus the first rigorous classification of seizure dynamics. Figure 1 demonstrates these different bifurcations, showing how the signal changes in terms of amplitude and frequency of successive spikes and may contain a shift in the baseline as the seizure starts or stops. Of note, in dynamical terms, a ‘spike’ is defined as any prominent sharp transient associated with the dynamical process. For human EEG, we assume this includes all fast transients < 200 ms with amplitude that is distinguishable from the background. Of note this dynamical definition also includes the fast, low amplitude spiking seen at the beginning of some seizures.

5In this work, we present algorithms to measure the invariant properties, which can then be used to classify the seizure types. While the theory behind this classification has been proven mathematically, measurement of these values under real conditions is challenging. This is because 1) EEG recordings of the brain are much more complex than single bursting cells, 2) EEG is notoriously noisy, and 3) there is limited understanding of the underlying physiology that produces the EEG waveforms. Despite these limitations, we previously found strong evidence that at least one dynamotype exists across multiple species. Herein, we present both an automated algorithm and a visual method to analyze these noisy data. We find that visual analysis is quite reliable and often preferable under clinical conditions, as demonstrated by recent work in other noisy neural signals. Therefore, while we do present the algorithm results as validation, the final clinical analysis is based upon the visual classifications.

6As described in Figure 1, there are four onset bifurcations. Two of these progress gradually from resting state into seizure via specific scaling laws for the amplitude or frequency: in the supercritical Hopf bifurcation (SupH), the amplitude of the oscillations starts at zero and increases proportionally to the square root of the distance from the bifurcation point; the Saddle-Node on an Invariant Circle (SNIC) bifurcation has oscillations that increase in frequency as the square-root of the same distance. The other two, Saddle-Node (SN) and subcritical Hopf (SubH) bifurcations, have abrupt amplitude and frequency changes that do not follow specific scaling laws. The SN can contain a jump in the signal baseline (i.e. direct current (DC) shift), but in the absence of a detectable DC shift these two types can be difficult to distinguish even theoretically. SN without DC shift (‘SN (-DC)”) and SubH are thus grouped together in this work. Further demonstration of the generation of these time series can be found in;.

7Two of the bifurcations are characterized by decreasing frequency, following square-root scaling for the SNIC and logarithmic scaling for the Saddle-Homoclinic (SH) bifurcation. Both these dynamics manifest as slowing of the seizure down to zero near its end, which is the well-known clinical hallmark of a seizure. This ‘slowing’ at seizure termination has been identified as an inherent characteristic of seizures across all spatial scales, and in multiple species and brain regions. Due to the small number of spikes near the end of seizures, it is difficult to distinguish between these two scaling laws when using only frequency of spikes; throughout this work, when we refer to ‘slowing-down’ it refers to logarithmic or square root scaling, implying scaling down to zero. The other two bifurcations do not require slowing at termination: the SupH bifurcation with square-root-scaled decreasing amplitude and no scaling law for the frequency, and the Fold Limit Cycle (FLC) that has no specific scaling law for either frequency or amplitude. Note that these two offset types may appear somewhat atypical to a clinician, as ‘slowing down’ is expected at the end of most seizures. The only dependence on a DC shift is that only the SH can have a shift, but the absence of a shift the slowing down pattern can be SH (-DC) or SNIC.

8We analyzed seizures from 120 patients recorded on intracranial EEG in seven centers worldwide (Appendix I.1) to identify the bifurcations at onset and offset. All patients had focal onset seizures. Results are shown in Figure 2.

9The canonical features necessary to distinguish the bifurcations are the dynamics of the amplitude and interspike intervals (ISI) and the presence/absence of a DC shift (Figure 2A). We identified the spike timing and amplitude to allow for both visual and mathematical analysis (Appendix I.2). In order to validate the analysis, we first compared the results of three human reviewers (Appendix I.3) and an automated algorithm (Appendix I.4) in a gold-standard computational model that generated 60 seizures of each type. These methods were very successful in distinguishing the different bifurcations in the model data, although we had to consider SN (-DC) and SubH as a single group for seizure onset, and SH (-DC) and SNIC as a single group for offsets because there is no method to distinguish them. We then compared these same methods on 120 human seizures. We found that concordance was also reliable in human data (Appendix I.5). These results show that the chosen features are capable of distinguishing the different bifurcations reliably for both human visualization and algorithms, that human seizure dynamics are consistent with the modeled bifurcations, and that human reviewers can use the methods described in the next paragraphs to classify onset and offset bifurcations reliably. In addition, we found that human reviewers were more reliable than the automated algorithm in noisy clinical data. We then used the human markings on the clinical data for the taxonomy in Figure 2. Note, however, that the identification of bifurcations from empirical data is notoriously difficult and generally cannot unambiguously prove that a given bifurcation is present, although it allows assessment of self-consistency. Further investigation may use additional tools, such as perturbations of the system, to corroborate these results further.

10Onset dynamics - To analyze onset dynamics, we first investigated all 51 patients that had been recorded using equipment capable of visualizing DC shifts. Of note, DC shift alone has previously been shown to be highly correlated to the seizure onset zone. Many seizures (41%) started with constant amplitude spikes and a DC shift, signifying SN bifurcation (Figure 2E). The second most frequent was similar amplitude/frequency dynamics without a DC shift (either SubH or SN, 37%), followed by SupH (14%) and SNIC (6%). Ambiguities in the classification were treated systematically as detailed in the Appendix (I.6), and only 1/51 seizures could not be agreed upon by the reviewers. Despite the inability to distinguish SN (-DC) and SubH, we demonstrate that at least three of the predicted bifurcation types are present at onset in human seizures. Importantly, some seizures display complex dynamics because they go through more than one bifurcation as the seizure begins. For example, nearly half of the seizures labeled as SN onset progressed into square-root amplitude scaling after 2–5 s, consistent with a switch to SupH dynamics; we only labeled the initial bifurcation herein. For this reason, we have slightly modified the definition of onset bifurcation as compared to the literature. Here the onset bifurcation is the one causing a departure from the resting state, even if it is not directly causing the onset of oscillations as there can be intermediate states (Appendix II.1).

11We next looked at the onset bifurcations in the other 69 seizures that were recorded with non-DC coupled hardware. Of the four onset bifurcations, SN and SubH become indistinguishable without a DC shift available. Combining these two into a single group, we found that 67% were SN-SubH, 26% were SupH, 6% were SNIC, and only one seizure did not have reviewer consensus.

12Offset dynamics – We first examined offsets in those patients with DC recordings for the most robust classification. The analysis of the interspike intervals (ISI) and spike amplitudes revealed a logarithmic/square-root slowing-down with constant amplitude in 20/51 patients, and of those 10 had DC shifts at offset. Thus, 10/51 (20%) were SH (+DC), while the remaining 10 (20%) were potentially SH(-DC) or SNIC (Figure 2F). The remaining 31 patients did not have slowing at the end of their seizure (53% FLC, 6% SupH, 2% no reviewer consensus).

13For non-DC coupled data, we grouped SH and SNIC in the remaining 69 patients. The majority of seizures had arbitrary ISI and/or amplitude (54% FLC), while 41% had slowing down characteristic of SH(-DC)-SNIC. The remaining seizures either had constant ISI with amplitude that decreased as a square root (SupH, 3%) or had no consensus among reviewers (3%). To supplement the examples in Figure 2, as well as to clarify how to approach several challenging scenarios, there is a primer with examples of all the different bifurcations in Appendix I.7.

14One significant challenge in the above analyses was the presence of noise. Within the theoretical model, ‘arbitrary’ dynamics refer to abrupt changes at seizure onset/offset without clear scaling laws to or from zero. We would still expect a smooth behavior for amplitude and frequency close to (not at) the onset/offset point, and specific trends (increasing, constant, decreasing) are possible even though not prescribed (Appendix V). However, in human data analysis ‘arbitrary’ includes a wide range of other behaviors, especially noise. It is important to note that the taxonomy above includes the first seizure that could be analyzed from every patient—we did not restrict the analysis to ‘clean’ seizures. Some of the seizures were noisy, either from technical concerns or physiological effects of the seizure. We chose this method in order to provide a robust, real-world demonstration of this analysis. Because noise can be classified as ‘arbitrary,’ this analysis may overestimate the numbers of SubH-SN (-DC) onsets and FLC offsets. However, it was clear that this limitation was not merely technical: several patients had complex physiological dynamics. Four of the FLC seizures were highly unusual from a clinical perspective: one had increasing amplitude for the last 10 s, one consisted of low-voltage fast activity that ended abruptly without any other change, one ended with irregular spike waves, and one had accelerating frequency at the end (Appendix I.8 and V). These examples highlight the vast heterogeneity of seizures, as well as the need for a valid taxonomy that allows for scientific discussion of the critical dynamical properties.

15Sixteen seizure dynamotypes– The preceding data validate that at least three types of onset and offset are systematically present in human focal epilepsy. As detailed above, real clinical data are challenging: the lack of DC shift makes it difficult to distinguish some bifurcations, and the noisiness of EEG is hard to distinguish from arbitrary dynamics. Nevertheless, these results show robust evidence that human seizures conform to both the onset and offset bifurcations predicted by our framework. These combinations lead to a taxonomy containing 16 dynamotypes of electrographic seizures. We identified the dynamotypes in patients with DC recordings. Two patients did not achieve reviewer consensus (one onset, one offset), leaving 49 patients (Figure 2F). We identified 12 different dynamotypes, with the limitation that several of the dynamotypes cannot be fully distinguished in the absence of a DC shift (e.g. SH(-DC) – SNIC offsets). The taxonomy was dominated by seizures with either SN or SubH onsets and slowing (SH-SNIC) or arbitrary dynamics (FLC) at the end. In this cohort of focal onset seizures, the SupH and SNIC onsets were less common, accounting for all four dynamotypes that were absent.

16We compared all available clinical metadata from patients with their dynamotype and found no correlation between seizure type and patient gender, pathology, or localization. There was a correlation with age, as older patients tended to have more SupH onsets (Appendix I.9). We also compared these results with a prior visual classification that identifies seven basic seizure onset patterns, and found 6/7 patterns without any apparent relationship to clinical data or pathology (Appendix I.10). There were no significant similarities between the dynamotype and the visual classification.

17While analyzing this dataset, we noted that one patient had two consecutive seizures belonging to different types: one supH/supH and one supH/SH, raising the possibility that an individual may express different types of seizures. This finding was surprising, as many clinicians assume that a person’s seizure should be ‘stereotyped’, that is consistent over time. In fact, multiple medical devices have been designed under the presumption that a patient’s seizures would have similar appearance over time. To test whether individuals display different types of seizures over time, we used a unique dataset from Melbourne in which patients had intracranial EEG recorded continuously for many months. We analyzed over 2000 seizures from 13 patients. Given the size of the sample, we limited the analysis to the most straightforward metric: the ISI at seizure offset to determine whether there was slowing at seizure termination. This allowed us to differentiate the SH/SNIC from the supH/FLC bifurcations (i.e. slowing-down or constant ISI). There were 658 seizures of sufficient length (>25 s) to measure the offset ISI. To be conservative, we only classified seizures as slowing-down or constant if such a determination was unequivocal and labeled the rest as ‘not assessed,’ meaning that the determination was not readily evident on brief visual inspection. As seen in Figure 3, all 13 patients expressed at least two offset patterns. Note that this is likely an underestimation of the heterogeneity in seizure types: these recordings did not contain DC coupling so onset bifurcations were not assessed, and we did not distinguish the offsets into all four types. Nevertheless, we can unambiguously conclude that individuals have seizures from different dynamotypes over time.

18The clinical data above show that seizures can be classified based on their onset/offset bifurcations and that a patient’s seizures may display multiple dynamotypes. To gain a deeper understanding of the relations between dynamotypes, we formalize these findings within a single unifying mathematical framework, which can account for all these behaviors. During seizures, the firing activity of neurons becomes organized, enabling the emergence of oscillatory activity that can be observed in electrographic recordings. This greatly reduces the degrees of freedom necessary to describe the observed activity, that is a small number of differential equations are sufficient to describe the collective behavior (Figure 4A). We here consider a system with the minimum number of variables necessary to produce oscillatory activity, two. Based on the parameter values, two states can be distinguished: resting (fixed point) or oscillatory (limit cycle). When these two states coexist for the same range of parameter values (bistability), transitions between them can be promoted by noise if the system is sufficiently close to a bifurcation. However, the statistics of ictal durations points to the existence of a deterministic process governing this transition for seizure offset and possibly for the onset. This can be achieved with the addition of a third variable acting on the timescale of ictal duration. We previously validated this approach with the ‘Epileptor,’ a set of five differential equations able to account for the most dominant dynamotype (SN/SH, also known as ‘square-wave’ bursting ). In the Epileptor, the transition from ‘normal’ to seizure state and back again, as well as the seizure dynamics, are controlled by a collective permittivity variable that evolves on a slow time scale. However, the Epileptor accounts in principle for only a single dynamotype and is not sufficient to explain the data presented above (note that systematic parameter variations also show a range of other bifurcations, although these variations are model-specific and not canonical). The fact that individuals can express different types of seizures over time leads to two predictions that must be included within a model of human seizure dynamics: different dynamotypes must coexist in the same model, and there must be an endogenous mechanism by which the brain can transition slowly between dynamotypes.

19Addressing these predictions within the framework of bifurcation analysis provides the entry point to propose a general taxonomy of seizure dynamics and postulate the existence of an ultraslow modulation. Previous mathematical work demonstrated that the procedure to build a minimal model for the SN/SH type provides a two-dimensional map for the parameter space of the fast variables on which all the six bifurcations can be placed. In effect, the map is a representation of the range of states in which a brain region can exist, oscillatory (ictal) and non-oscillatory (interictal), and the transitions between them. The oscillatory state produces spiking activity that is described by the fast variables (millisecond scale activity). However, on a slower time scale of the order of seizure length, the brain can move toward a transition to an interictal state, as described by a slow variable. Migrations to different locations on the map can occur on a usually even slower timescale (10’s-1000’s of seconds), which we call here ‘ultraslow’. Saggio et al. showed that the use of an ultraslow variable allows full exploration of the map. Applying these general mathematical principles to epilepsy implies that any brain region able to generate SN/SH seizures can potentially generate other types by navigating to different dynamical regimes (i.e. changing the parameters of differential equations). That work also showed that a large number of physiological neuron and neural population models can be mapped upon a canonical dynamic model under certain mild conditions (existence of a Bogdanov-Takens point). All physiological parameters are then absorbed in only three generic parameters, which span a three-dimensional parameter space, in which all bifurcations are represented. Detailed bifurcation analysis reveals that all neighborhood relations between bifurcations can be displayed (without loss of generality) as projections onto the spherical surface within the parameter space, yielding a canonical 2D map, shown in Figure 4B–C. This map displays the basic topology of all possible relations between bifurcation lines, including identity of the bifurcation and the organization of its nearest neighborhood including proximity and intersection of bifurcations. As one (or multiple) of the parameters is continuously varied, trajectories are traced out in this map, eventually connecting two bifurcation lines and thus establishing a seizure’s dynamotype. The navigation of the canonical 2D map (that is a flat projection of a spherical surface) through continuous parameter changes selectively generates dynamotypes and justifies that we call it a seizure map. The state of a brain region at any moment can be represented as a location on the map, which defines its dynamical properties. Regions in the map that correspond to different regimes, including quiescent and ictal states, are separated by bifurcation curves. Seizures are represented as black arrows, each arrow corresponding to one dynamotype. To produce a seizure, the system, which is initially in the quiescent state within the bistability region, heads toward the onset bifurcation curve. When this curve is reached, the quiescent state disappears and the system is forced to go into the oscillatory seizure regime within the bistability region. This transition in state causes an inversion in the trajectory of brain state, with the system now heading toward the offset bifurcation curve. When the offset is reached, the system goes back to rest and inverts direction again. Movement along the black arrow is produced by slow (of the order of the ictal length) mechanisms leading to seizure offset. Note that the movement towards the onset and offset bifurcations at this timescale occurs in both cases from within the bistability region. Ultraslow movements on the seizure map are responsible for changing the location of the brain state while at rest (as may happen during the night and day cycle) and enable the expression of different types of seizures as observed clinically. This framework thus provides a potential explanation for the clinical observation of multiple types of seizures in a single patient, and the seizure map provides a hypothesis to describe how a patient’s current state (i.e. location on the map) can affect seizure dynamics (whether a seizure is likely to occur, and what type is most likely). Figure 4C depicts paths (black arrows) for seven of the 16 dynamotypes placed on a two-parameter map. Adding one additional parameter allows this map to be extended and create seven other types in three dimensions, while the final two types require even higher dimensions to create. These higher dimensional types require very fine parameter tuning, and thus are less likely to occur. TSD does not predict the likelihood of dynamotypes, but in conjunction with the seizure map and choice of slow dynamics, a hierarchy of seizures can be established, which is supported by our clinical data (Appendix II. 5), for instance the dynamotypes that occurred the most (e.g. SN/FLC and SN/SH) were predicted to be among the most likely to occur.

20It is important to distinguish the two types of fluctuations within the brain map. The slow permittivity variable affects the general brain state, or position on the map, on the scale of minutes to hours. It represents underlying, and sometimes varying, conditions of the system than determine the position in the map, which has broad physiological implications. There are also fast fluctuations on the scale of ms to s, better described as perturbations of the brain state from its current location on the map. These perturbations are modeled as ‘noise’ within the model, but in reality they also include many physiological phenomena such as afferent signals and neural potentials—in effect anything that perturbs the system. Thus, although the model refers to the addition of ‘noise’ to the system, these effects can readily be attributed to physiological neural activity. Both types of fluctuations can lead to seizures by pushing the system across the bifurcation.

21It is important to note that the topology of the map in Figure 4 was initially proven mathematically to be generic and rigorously valid for bursting. This invariance establishes the ground truth to define the relationships between the different bifurcations in the proximity of the SN/SH type, which leads to a key prediction: transitions between certain types may be more common due to their proximity on the map. For example, considering the bistability region in the upper part of the map, we note that the offset curves of SH and SupH approach until they meet. When the curves are very close, even small fluctuations in the parameters can cause a transition between types. If fluctuating internal conditions allow individuals to move around these regions of convergence, patients may have seizures belonging to different types over time, as observed in our longitudinal analysis. The model predicts that transitions between specific types are more likely to occur if these types are close in the map, in the sense just shown of bifurcation curves approaching each other within the same bistability region. On the contrary, transitions between types belonging to distant bistability regions require stronger changes in the ultraslow permittivity variable(s) and are thus less likely to occur.

22Two dynamotypes are paired when they share the same seizure onset or the same offset, and a continuous change of parameters in the seizure map can lead from one type to the other, sometimes even during a seizure. The seizure map predicts possible pairings as motions toward different bifurcations while a seizure is ongoing. As proof of concept, we found several examples of such fluctuations in our cohort. In Figure 5A, one patient’s seizure had constant ISI and square-root amplitude scaling for approximately 70 s, properties exhibited when approaching the SupH bifurcation. The seizure appeared to be terminating, but then it abruptly restarted to terminate with slowing-down ISI and constant amplitude (SH/SNIC bifurcation). We found five examples of this behavior in our data (out of >2000 seizures) and reproduced it with our model. By definition, the dynamotype includes only the onset and final offset bifurcation, but the behavior during this seizure is intriguing and can be explained by the model. We considered a path for the SN/SupH type with an ultraslow drift of the offset point that changed the path to SN/SH and added noise to all variables to simulate fluctuations. With these settings, we ran 100 simulations, which generated several different dynamotypes, predominantly SN/SH, SupH/SH, and SN/SupH (Appendix III). Several had transitions during the seizure from one bifurcation to another, 10 of which in the same manner as the data in Figure 5A, switching from SupH to SH offset. Thus, the clinical example is one of the most favorable combinations within the model, as the SupH and SH bifurcations are so close that small fluctuations can cause the switch.

23Status epilepticus: Explorations of the seizure map also demonstrated another effect sometimes seen clinically: status epilepticus. Simulations with the previous settings in some cases produced continuous seizures that did not resolve by the end of the simulation, equivalent to status epilepticus. We analyzed the corresponding trajectories on the map to determine how this had occurred. Prolonged seizures occurred when the brain state crossed the SN onset curve but was unable to return to rest through the offset bifurcation and remained mainly in the violet ‘seizure only’ region in Figure 5G. The slow variable naturally drives the state toward offset, but in these cases was continually overridden by noise, causing the state to ‘escape’ from the bistability region. We then analyzed this effect by simulating various levels of noise and showed that there was a clear correlation between the noise variance and the likelihood of entering status epilepticus (Appendix IV). We compared these results with our clinical data, which had two examples of non-convulsive status epilepticus. In both cases, the seizures began in typical fashion, but instead of terminating began having long periods of constant ISI with varying periods of amplitude fluctuations. There were many abrupt transition periods during which the ISI and amplitudes became arbitrary. After these brief periods of disorganization, the dynamics returned to constant ISI. We compared the dynamics of our model results with these human seizures and found that the transition between different dynamics is quite similar. In Figure 5D, we show a portion of two human seizures and one example of the simulation (Figure 5E–F) and movement on the map (Figure 5H–I). Further demonstration of the patients’ status epilepticus is provided in Appendix IV. The patterns of organized- alternating with disorganized- firing, often known clinically as ‘waxing and waning seizures,’ are entirely consistent with the model: the seizure undergoes periods in which it progresses toward termination, then because of noise it reverts to a point farther away from the offset bifurcation, as described previously.

24Accelerating seizure: We then analyzed the seizure offset that increased in frequency described in Appendix V. In this case, we explored conditions on the map that could produce ‘speeding up’ at the end of the seizure. We identified multiple trajectories in the map of brain states capable of producing these unusual seizure dynamics. As in the case of status epilepticus, these unusual patterns are dependent upon the relative position within the brain map, in this case occurring when brain states along a trajectory are affected by multiple bifurcations that are in close proximity (Appendix 1—figures 17–19). These results demonstrate the explanatory value of the seizure map and show how it provides a rational explanation for a wide range of physiological dynamics.

25Seizures have been recognized clinically for millennia, but after nearly a century of electrographic recordings we still do not have a translatable method of characterizing their dynamics. We here address this issue and provide the first principled approach toward the organization of seizures in a Taxonomy of Seizure Dynamics (TSD). TSD establishes 16 dynamotypes of seizures, which could be extended to more exotic dynamotypes when considering non-planar bifurcations (see Appendix II.6). As TSD provides the classification of seizures, it remains completely unbiased to each seizure type. This invariance is broken by the seizure map, which establishes relations between dynamotypes and introduces a bias in the taxonomy, laying the grounds for a hierarchy of dynamotypes. The hierarchy is based on the mathematical consequences of how bifurcations are related to each other. The relations can be considered as structural in the sense that they rely on the static properties of location, shape, branching and topology of bifurcation curves in the seizure map. The implications are functional in the sense that they determine the non-static properties of a seizure’s discharge patterns including frequency, acceleration/deceleration, amplitude and amplitude growth. As such TSD and the seizure map provide another example of the ubiquitous link of structure and function in biology.

26Our classification aims at precisely identifying the seizure type in terms of dynamics, without any dependence upon specific symptoms, pathology, or localization. Thus, it is highly complementary to the classical operational classifications used by clinicians to diagnose and treat patients, which are based upon those factors without addressing dynamics. The dynamotype describes the behavior of the seizure itself, while the clinical classification describes the patient’s symptoms: together, both classifications are synergistic and can be used to improve patient stratification, providing more insight into diagnosis and treatment. TSD is based on simple, invariant, objective metrics that have compelling scientific rationale. With DC-coupled recordings, it is possible to distinguish the types with high fidelity, even with visual inspection. This method is thus readily available to clinicians, as many standard EEG acquisition devices now have excellent resolution near DC (<0.1 Hz).

27Our interpretation of the results relies on some key assumptions. First, we are assuming that the onset and offset of seizures are brought about by bifurcations. Another common mechanism in the literature is noise-induced transitions, which our model can reproduce (see Appendix II.3 for a discussion of how this would affect our classification). Second, we rely on the assumption of timescale separation between the dynamics of the spikes within seizures and the slower dynamics controlling seizure threshold and termination (both in theory and in the simulated sample of data). If this assumption does not hold, different phenomena could occur and the scaling laws could be impossible to identify. Third, we only considered planar bifurcations for simplification (Appendix II.6). These assumptions were the axioms for developing the theory and data analysis. Future work will address the validity and consequences of these simplifying assumptions.

28Similar to sleep, seizures are universal from insects to humans, leading to the proposal that seizures are an inherent property of a brain, that is they are endogenous to the brain, perhaps as an emergent property of complex neuronal networks. This may explain why, despite the vast range of genetic, structural, chemical, and developmental conditions that cause epilepsy, seizures have a remarkably limited set of dynamical behaviors. It is therefore not surprising that elementary mathematical laws can describe their electrophysiological signature. However, clinical interpretation of seizure dynamics has been almost universally based on simple observation, in which clinicians report the frequency and morphology of spikes. This method is helpful to identify primary generalized epilepsies, but within focal seizures has limited clinical use. TSD by itself does not have any bias either regarding the dynamotype, but when linked to the seizure map shows certain dynamotypes to be more prevalent than others. Empirically, low voltage fast activity and focal DC shifts have been found to be highly predictive of the true seizure focus, which are both patterns corresponding to common dynamotypes, suggesting clinical utility.

29One limitation of previous clinical descriptions of seizures is that it has been unclear which dynamical features are relevant. There is high variability in the frequency and morphology of spikes due to individual fluctuations and noise, as seen in spontaneous seizures recorded in humans and experimental models. Our analysis identifies the invariant dynamics that impose important constraints on the system. A crucial aspect of our approach is that it allows us to disentangle characteristics that are necessary to describe the dynamics from other seizure related phenomena that are not fundamental to our simplified model (e.g. spike and wave complexes, preictal spikes, sentinel spikes (see Appendix I.2)). This is not to suggest that such biomarkers are not relevant to epilepsy: they have well-known correlation with the epileptogenic zone and can predict the occurrence of the first spontaneous seizure during experimental epileptogenesis. The Epileptor model, which comprises the SN-SH dynamotype, generates both spikes and seizures. The spikes are important indicators of the organization of the networks, but not part of the generic features of onset and offset dynamics and thus not appear in the generic model. If additional mechanisms are introduced based on other forms of reasoning, as in via a second population, spikes can be included in the dynamics. It is important to note that we only analyzed the particular case of drug-resistant epilepsies investigated with invasive intracranial recordings; however, it is unlikely that this theory is specific to such epilepsies as seizures forced in non-epileptic networks follow the same universal rules. We note that prior work focusing on individual bifurcations are all also entirely consistent with TSD, which encompasses all these onset and offset possibilities and further shows how they interact (Appendix VII).

30This work does not include data from generalized onset epilepsies, as these are not typically recorded with intracranial EEG. Our taxonomy is, however, fully consistent with past work on generalized epilepsy dynamics. Absence seizures, for instance, begin with sudden onset and offset of ~3 Hz large amplitude oscillations without a DC shift and terminate abruptly without slowing down to zero, which would point to a SubH/FLC type (which is the most likely dynamotype, see Appendix II.5). In this work with focal epilepsies, the most common dynamotypes were the SN (+DC)/FLC and SN/SH (which are also likely types).

31When comparing our results with a past visual classification system of spike frequency, we found no correlation with pathology in our cohort of 120 patients. However, our cohort did not have any patients with tuberous sclerosis, which was the only pathology associated with burst suppression in that prior work. Combining the data from both studies, the different patterns appear to be either evenly distributed or too rare to find robust correlations with pathology. Similarly, dynamotypes are not strongly correlated with pathology. In terms of the seizure map, we hypothesize that what determines the seizure dynamics is not the pathology per se, but the location of the brain on the map. Specific pathologies may predispose to certain regions, but there are many complex dynamics affecting brain state and many conditions that can produce similar dynamics. This coincides with the idea of the seizure map showing the full range of potential seizure onset and offset activity.

32There is great clinical and research potential in characterizing a seizure’s dynamotype, as it provides a unique perspective on brain networks. The current standards of epilepsy care focus on phenotype, genotype, and the time/location of seizure onset. While those methods have obvious utility, they do not address the underlying dynamics and thus have left several questions unanswered for decades. How do seizures start, stop and spread? How do we tell the difference between inter- or pre-ictal spiking and seizure initiation? Is it possible to measure the distance to seizure threshold, that is determine seizure risk at a given moment? How do we compare two different seizures? Is it possible to measure if a treatment is working by moving the brain ‘farther away’ from seizure onset, rather than waiting to see if seizures recur? These questions all require an understanding of the dynamics—an understanding that is not addressed by the current clinical tools. This is where the utility of the dynamotype is manifest. At its most basic, the dynamotype is a quick description of the key dynamics of a particular seizure, a clinical language that focuses on the aspects that are most important. This would supplement current visual descriptions, which typically are limited to amplitude and frequency. But there are many deeper applications of this tool as well. We have previously demonstrated that very different biophysical mechanisms can produce the same dynamotype. Here, we show that seizures from 120 patients contain almost the entire taxonomy of dynamotypes, and that a wide array of focal pathologies can be grouped into similar dynamotypes. Our interpretation is that this is because the seizures depend heavily on ‘local’ dynamics, that is the current brain state (location on the map) and acute perturbations (noise in the system), more than that a single pathology would necessarily predispose to a specific location on the seizure map.

33There are many other potential applications for TSD in basic research as well. For instance, our group recently published an analysis quantifying how epileptogenesis progresses in the tetanus toxin model in rats. That work showed that the dynamotype evolved over time, starting with SN and moving into SNIC (Appendix 1—figure 4) and sometimes SupH onsets over the course of weeks. A clinical trial is currently underway in France using the SN-SH dynamotype to model seizure foci and spread. Future versions of such tools could utilize the whole taxonomy to be much more comprehensive, tailoring models to the key underlying dynamics of specific patients. These models would greatly enhance modern network analytic tools, which would be greatly enhanced with a rational model to describe the underlying dynamics.

34One novel aspect of dynamotype is that understanding the underlying dynamics can help in the design of strategies to control seizures, such as with electrical stimulation. Studies on neuronal bursters, which are organized in similar dynamic types, demonstrate that types have different sensitivity to stimulation. For example, SubH onset acts as resonators, which require a resonant frequency in the stimulus to trigger oscillations, while SN onset behaves as an integrator in which the nature of the stimulus (excitatory or inhibitory) rather than the frequency plays a key role. There is a long history of using perturbation to probe the proximity of a nearby bifurcation in disciplines such as electrical power and reservoirs. Past work on stimulation to assess epileptogenicity is similar to such work and would be greatly enhanced with the insight gained from this model to understand the nearby bifurcations. There is a long history of bifurcation research in other fields that may also be applicable to seizures, such as using perturbations to assess proximity to a SubH or SupH bifurcation. Further theoretical and clinical work is necessary to assess whether knowledge of the dynamotype could also improve the ability to abort seizures with tailored stimulation. The second important prediction is that the synchronization properties of coupled bursters are bifurcation-dependent. This is a key issue for seizure propagation, as it predicts that the ability of a seizure to spread is dependent upon its type. Since the spatiotemporal organization of the seizure is part of the data features used to personalize brain network models (Virtual Epileptic Patient; ) and functional connectivity based approaches, the choice of the right dynamotype is critical for successful patient modeling and clinical translation.

35Another significant contribution of this work concerns the dynamics of the slow permittivity variable to explain how slow changes in the behavior/state of a brain region can bring it closer or farther away from different bifurcations, that is seizure threshold. The fact that all 13 patients had seizures belonging to at least two types implies that the permittivity variable moves on a dynamical map in which different types of bifurcations can occur. Each parameter of the map should be considered as a representation of a manifold of physiological variables that cooperate to produce a particular change in the system. Given the slow timescale at which these changes occur, neurochemical substances (e.g. hormones, neuromodulators etc.) are the best candidates. Within the permittivity variable we can here distinguish two timescales: a slow timescale of the order of the ictal length, and an ultra-slow timescale of the order of the interictal length (hours, days, months, years). Typical examples include the circadian regulation of seizures and catamenial epilepsy. Interestingly, both males and females display ultraslow (weeks) modulation of seizure probability both in rats and humans, further suggesting that these results are species- and sex-independent. Those results and ours strongly support the proposal that patients move closer and farther away from seizure threshold (i.e. ‘travel the map’) during their lifetime. This interpretation may also be helpful in assessing a brain’s current proximity to seizure bifurcations, that is predict the risk of seizures occurring. Several features are altered when nearing the onset bifurcations, such as preictal spikes, variance of the signal, and reaction to electrical probing of cortical excitability. Recent work has shown that interictal discharges act like system perturbations that behave like the slow approach to bifurcations, just as predicted by our model. Seizure forecasting, based on electrographic recordings, is already enhanced when circadian rhythms are used to inform the model. If the ultraslow physiological correlates of the map’s parameter could be identified, measured, and manipulated, this would open new possibilities to assess when the patient is moving toward unsafe regions of the map and to alter their trajectory, that is control seizures before they occur.

36Our proposed framework provides an organizational principle of seizure dynamics, which, when linked to canonical dynamic systems, identifies a generic seizure map charting out characteristics of dynamotypes including prevalence of a dynamotype and possible pairings of dynamotypes that can occur during a seizure, as well as others that are prohibited. TSD does not describe all possible seizure features, but relies on seizure onset and offset classification, which is why we defend a complementary approach, combining it with traditional operational classifications. However, it provides a unique avenue to classify seizures based upon their key dynamical features, while providing insight into how seizures become more or less likely to occur at a given time. A corollary is that, although TSD has been developed in the context of seizure dynamics, it likely extends to physiological function of the healthy brain (e.g. alterations between REM and slow wave sleep, and the appearance of gamma frequencies or ripples during slow wave sleep) and stipulates the existence of at least two time scales in any theory of the brain. Slow time scales are present in many theories of brain function, but typically have been limited to the domain of learning and adaptation, thus functionally separated from fast processes. Here, the functional integration and co-evolution of the fast neuroelectric and slow permittive time scales suggests emergent and inherent properties of brain processes.