Seizures can occur spontaneously and in a recurrent manner, which defines epilepsy; or they can be induced in a normal brain under a variety of conditions in most neuronal networks and species from flies to humans. Such universality raises the possibility that invariant properties exist that characterize seizures under different physiological and pathological conditions. Here, we analysed seizure dynamics mathematically and established a taxonomy of seizures based on first principles. For the predominant seizure class we developed a generic model called Epileptor. As an experimental model system, we used ictal-like discharges induced in vitro in mouse hippocampi. We show that only five state variables linked by integral-differential equations are sufficient to describe the onset, time course and offset of ictal-like discharges as well as their recurrence. Two state variables are responsible for generating rapid discharges (fast time scale), two for spike and wave events (intermediate time scale) and one for the control of time course, including the alternation between ‘normal’ and ictal periods (slow time scale). We propose that normal and ictal activities coexist: a separatrix acts as a barrier (or seizure threshold) between these states. Seizure onset is reached upon the collision of normal brain trajectories with the separatrix. We show theoretically and experimentally how a system can be pushed toward seizure under a wide variety of conditions. Within our experimental model, the onset and offset of ictal-like discharges are well-defined mathematical events: a saddle-node and homoclinic bifurcation, respectively. These bifurcations necessitate a baseline shift at onset and a logarithmic scaling of interspike intervals at offset. These predictions were not only confirmed in our in vitro experiments, but also for focal seizures recorded in different syndromes, brain regions and species (humans and zebrafish). Finally, we identified several possible biophysical parameters contributing to the five state variables in our model system. We show that these parameters apply to specific experimental conditions and propose that there exists a wide array of possible biophysical mechanisms for seizure genesis, while preserving central invariant properties. Epileptor and the seizure taxonomy will guide future modeling and translational research by identifying universal rules governing the initiation and termination of seizures and predicting the conditions necessary for those transitions.
Epilepsy is characterized by the occurrence of spontaneous seizures, which arise when large regions of the brain produce uncontrolled, synchronous neural activity. Partial-onset seizures form the most common form of epilepsy, which is most likely to be drug-resistant (Brodie et al., 2012). Epilepsy can exist by itself or be associated to other neurological disorders including Alzheimer’s disease (Friedman et al., 2012), autism (Robinson, 2012) and Down’s syndrome (Arya et al., 2011). Despite the fact that so many different pathologies, conditions and network reorganizations can result in partial epilepsy (Pitkanen and Sutula, 2002), one observation is particularly striking: the electrophysiological signature of different seizures is remarkably similar from case to case, even among primitive laboratory models. For example, different seizure-onset patterns are common to different epileptogenic lesions (Perucca et al., 2013). Although seven seizure-onset patterns can be distinguished, two major categories of activities are consistently found: fast oscillations and spikes with or without waves (Perucca et al., 2013).
Another striking property of seizures across species (from flies to humans) is the possibility of triggering them in any ‘normal’ brain using an array of inducing conditions. In humans, sleep deprivation, stress, electroshock treatment or toxins can evoke seizures (Luttges and McGaugh, 1967; Nakken et al., 2005; Jett, 2012). In animals, electrical stimulation (kindling) or administration of various chemical compounds is commonly used to trigger seizures in vivo (Raol and Brooks-Kayal, 2012). A large variety of protocols can be used in vitro to produce ictal-like events, including in human slices (Huberfeld et al., 2011), demonstrating that even small neuronal networks can be forced into a ‘seizure’ state. Although the conditions needed to induce them may be very different, the electrophysiological signature of such seizures is remarkably similar to those recorded in vivo, including the presence of fast oscillations and spikes.
As seizures can occur under such diverse conditions, including in ‘normal’ networks, they belong to the dynamic repertoire of brain activities, as do other types of oscillations (e.g. theta, gamma etc.). Based on their apparent stereotypy, we hypothesize the existence of dynamical properties that would be invariant in most spontaneous and evoked seizures across brain regions and species. The notion of invariance is key to our approach. In non-linear dynamics, one form of invariance produces bifurcations, which are transitions from one type of behaviour to another. The most generic type of bifurcations are local, which are invariant under transformation and can often be described by canonical models (typically differential equations). Canonical models define the minimal requirements necessary for a generic behaviour to arise (Hale and Koçak, 1991; Kuznetsov, 1998).
In the first part of this article, we test our hypothesis of invariance in seizures. We use a simple in vitro model system and systematically characterize transitions between normal and epileptic states. From this characterization, we develop a taxonomy of epileptic seizures. For one particularly prominent class of seizures, we derive the bifurcations at seizure onset and offset and construct a set of differential equations for the complete seizure like event (SLE), which defines our model for seizure evolution, the ‘Epileptor’. This process of developing canonical models based upon generic properties of a system has been successful in describing dynamical phenomena such as ferromagnetism, lasers and superconductivity in physics (Cross and Hohenberg, 1993) and neuronal discharge patterns (spiking and bursting) in biology (Ermentrout and Terman, 2010). Their interest lies in their abstract nature, as they do not depend upon a detailed knowledge/identification of biophysical properties and still enable the identification of general rules. In the second part of this article, we test the rules derived from the Epileptor to account for seizure dynamics in other species and diverse brain regions for different types of epilepsies in patients. Finally, predictions from Epileptor will be tested experimentally, particularly exploring different pathways to seizure onset. Our results reveal that seizures are a simple, conserved behaviour of brain networks that can be systematically classified through their onset and offset bifurcations, suggesting that anti-seizure strategies could involve altering dynamical properties rather than specific pathways.
Materials and methods
All protocols were designed and approved according to INSERM and international guidelines for experimental animal care and use. Experiments were performed on intact hippocampal-septum preparations taken from FVB NG mice between postnatal Day 5 and 7 (day of birth = Day 0). Animals were sacrificed by rapid decapitation and brains were extracted and transferred to oxygenated (95% O2/5% CO2) ice cold (4°C) artificial CSF containing: 126 mM NaCl; 3.5 mM KCl; 2 mM CaCl2; 1.2 mM MgCl2; 25 mM NaHCO3; 1.2 mM NaHPO4; and 10 mM D-glucose (pH 7.3). The two hemispheres were separated and dissected to obtain hippocampal-septum preparations (Khalilov et al., 1997). They contained the hippocampal formation, septum and parts of adjacent neocortical areas. Preparations were transferred to a chamber with artificial CSF at room temperature. After at least 2 h incubation, preparations were transferred to the recording chamber. The hippocampus and the still-connected septum were placed into two chambers and perfused with different media (Khalilov et al., 2003). To ensure a high level of oxygenation, the preparation was perfused at 15 ml/min with artificial CSF warmed to 33°C bubbled with CO2/O2 gas mixture. The pH was controlled during each experiment. After a 30 min baseline recording in Mg2+-containing artificial CSF solution, the media was changed to one without added Mg2+ in the chamber containing the hippocampus. Under this condition, the extracellular concentration of Mg2+ does not necessarily decrease to zero because there is minimal Mg2+ contamination from other constituents of the artificial CSF, reaching up to 0.08 mM (Mody et al., 1987). This very low residual level does not prevent the genesis of epileptiform events. In other experiments, external [K+] was increased by adding increments of 1 mM KCl to artificial CSF. Mannitol was added in the same way in increments of 10 mM. Extracellular recordings were performed using glass extracellular electrodes filled with low Mg2+ artificial CSF placed in the CA1 stratum oriens region. Field potentials were amplified with isoDAM-8A differential amplifiers (World Precision Instruments), which allows direct current (DC) recordings (low pass filter set to 3 kHz), digitalized with DigiData1200 converter (Molecular Devices) stored on the hard drive of the personal computer using PClamp 8.2 software (Molecular Devices).
Interneurons and pyramidal cells were blindly recorded or identified using infrared-differential interference contrast microscopy through a ×60 water immersion objective. Microelectrodes had a resistance of 4–8 MΩ, and an internal solution of the following composition was used to record excitatory and inhibitory postsynaptic currents: 135 mM Cs-gluconate; 2 mM MgCl2; 0.1 mM CaCl2; 1 mM EGTA; 2 mM MgATP; 0.5 mM Na4GTP; 10 mM HEPES; and 0.5% biocytin (pH 7.3; 270–280 mOsm). Access resistance was monitored throughout the experiments (range 12–30 MΩ). Experiments were discarded if series resistance increased by more than 20%. Cell attached recordings were first performed to record the firing activity of the recorded cells. For voltage clamp experiments, cells were kept at −60 mV or +10 mV for the analysis of glutamatergic or gamma aminobutyric acid (GABA)-ergic spontaneous postsynaptic currents, respectively. These currents were sensitive to D-APV/NBQX [D-2-amino-5-phosphonovalerate/2,3-dihydroxy-6-nitro-7-sulfamoyl-benzo(F)quinoxaline] and bicuculline, antagonists of NMDA/AMPA (N-methyl-D-aspartate/α-amino-3-hydroxy-5-methyl-4-isoxazole propionic acid) and GABAA receptors respectively. Some cells were recorded in current clamp to monitor membrane voltage. In these experiments, the Cs-gluconate-based solution was replaced by a K-gluconate-based solution. All data were acquired using an analog-digital converter (Digidata 1322B, Molecular Devices) and analysed using Clampfit (Molecular Devices) or Matlab (Mathworks)-based software. During recordings, all neurons were passively filled with biocytin for post hoc morphological identification (Quilichini et al., 2012).
Recordings of changes in the extracellular [K+] were performed with double-barreled K+-sensitive and reference microelectrodes manufactured and calibrated as described previously (Heinemann et al., 1992). In brief, electrodes were pulled from double-barrelled theta glass. The reference barrel was filled with 154 mM NaCl solution, the ion sensitive barrel with potassium ionophore I cocktail A60031 (Fluka) and 100 mM KCl. Ion-sensitive microelectrodes with a sensitivity of ∼8 mV/mM of [K+] were used for experiments.
Glass oxygen microelectrodes (OX-10, Unisense) were calibrated and then placed at a 75–100-µm depth into the CA1 stratum oriens region of the hippocampus close to the field electrode. The data were then aquired synchronously with the local field potential.
Changes in NADH (or FAD) fluorescence in hippocampal preparations were monitored using a 330±40 nm (450±40 nm) band pass excitation filter and 420 nm (520 nm) long pass filter for emission (OmegaOptical). The light source was the Intensilight C-HGFI illuminator (Nikon Instruments Europe B.V.) equipped with a mercury arc lamp. Hippocampi were epi-illuminated and imaged through a Nikon upright microscope (FN1, Eclipse) with ×4/0.10 Nikon Plan objective (Nikon Instruments Europe B.V.). Images were acquired using a linear, cooled 12-bit charge-coupled device camera (Sensicam, PCOAG) with a 640 × 480 digital spatial resolution. Because of a low level of fluorescence emission for the fluorophores, NADH and FAD images were acquired every 600 ms as 8 × 8 binned images. The exposure time was adjusted to obtain a fluorescence intensity between 2000 and 3000 optical intensity levels. Images were stored in a computer as 12 bit files (effective spatial resolution of 80×). The recording sites were extracted in three to five regions of interest using ImageJ software (developed by Wayne Rasband, National Institutes of Health). Data were expressed as the percentage changes in fluorescence over a baseline [(ΔF / F) × 100%].
Physiological recordings contain many types of noise, but for the purpose of the paper we concentrate on the random activity produced by uncorrelated synaptic events, or ‘synaptic noise’, which is the activity we modulate by adding KCl to the septum above. We measured several characteristics of this noise, comparing the levels just before an event (pre-ictal, 1–60 s) with the baseline levels (interictal, >60 s). Both current clamp (mV) and voltage clamp (pA) recordings provide similar data for noise analysis, though only in the former can cellular depolarization be assessed.
Amplitude distributions were fit to a Gaussian curve. Power spectral density was obtained using the ‘pwelch’ method in Matlab with 10 000 samples and 100 overlap. Noise colour was tested by fitting the Power spectral density > 10 Hz to the equation: Power spectral density = A/frequency^k (Destexhe et al., 2003), with typical values of k ∼ 1–2.5 for neural noise activity. The White test was used to determine whether noise increased as seizures approached, i.e. testing for heteroskedasticity (White, 1980). Cross correlations of the signals were used to assess whether there was any non-random autocorrelation. Noise intensity was measured as the second moment of the patch clamp data about the median amplitude (Stacey et al., 2009), summed over 0.1 s. Total variance in the pre/interictal periods was the average variance of local field potential of all 0.1 s bins.
Spike time data
In this analysis and elsewhere, ‘spikes’ refer to population transients. For noise analysis, the detected spikes referred to detected afferent signals that produced post-synaptic potentials. Raw data were bandpass filtered (0.5–20 Hz, bidirectional elliptic), passed through a peak detector and spikes identified by manually assigning a threshold. These values were chosen manually to maximize automated detection of post-synaptic potentials based upon visual inspection of 2-s segments. The time between all successive spikes (interspike interval) was fitted to the lognormal distribution to assess whether it the spikes occurred as a Poisson process.
Mann Whitney U and t-tests were performed to compare pre- versus interictal data: raw voltage, spike time, noise intensity and local field potential voltage.
Interspike interval scaling
All species tested (human, mouse, rat and zebrafish) were analysed in identical fashion. The spikes detected in this analysis refer to the epileptic ‘spike and slow wave’ activity seen on traditional EEG or local field potential recordings. The recorded voltage was first bandpassed (4–100 Hz, bidirectional elliptic) and peaks/troughs of epileptic spikes identified using the ‘findpeaks’ function (Matlab, Mathworks) and manual thresholding. The parameters were chosen manually to maximize automated detection of epileptic spikes (filtering out slow wave) based upon visual inspection of 10 s segments of recording. A trained clinical epileptologist manually determined those parameters, as well as seizure start and end times by visual analysis of every seizure. The relationship between the interspike intervals and time until the end of the seizure was evaluated by equation fitting. For this relationship, the interspike interval was the latency (s) between consecutive spikes, and the time was the duration from the first spike in the pair to the end of the seizure (s). All times were positive numbers. Data were fit to several equations using a least-squares fit algorithm in the ‘cftool’ (Matlab, Mathworks). Suitability of logarithmic scaling was assessed by comparing summed squared error, adjusted R2, and qualitative features among the different models (Supplementary material). For visualization purposes, the plots of these data were oriented the same way as the EEG data (Fig. 6). Thus, the x-axis for the interspike interval plots is reversed with the lowest numbers (end of the seizure) at the far right. The first interspike interval in a seizure occurs at the far left, with an x-value equal to the total duration of the seizure.
Human seizures were recorded using standard clinical procedures for intracranial monitoring, and the patients consented to a data sharing agreement as approved by the local ethical committee. Data were recorded using a standard 128-channel clinical acquisition system: XLTek, Inc.: 0.1 Hz high pass filter, 100 Hz low pass filter, 512 Hz sampling rate (Worrell et al., 2008). As is typical for such recordings, the data are recorded with alternating current coupling, which removes any DC component. The data were de-identified, posted and downloaded from the International EEG Database (http://www.ieeg.org). All clinical data for each patient shown is available from the IEEG database. Studies used and clinical data are listed in Supplementary Table 1.
Hyperthemia-induced seizure-like events were kindly provided by Dr S. Baraban. The experimental procedure is described in (Hunt et al., 2012).
Low Ca2+ recordings
Seizure-like events recorded when lowering extracellular Ca2+ were kindly provided by Dr. J. R. Jefferys and Dr. Jiruska Premysl. The experimental procedure is described in Jiruska et al. (2010).
Basic building blocks of seizure dynamics
As we searched for invariant features in seizure dynamics, any experimental model system may be used. Here, we analysed SLEs recorded in the intact immature hippocampus of mice in vitro, as this preparation is easily amenable to experimentation. When placed in continuous epileptogenic conditions, the preparation produces spontaneous recurrent SLEs (Fig. 1, unfiltered record). In Fig. 2A, we show one such typical spontaneous SLE (0.01 Hz high-pass filter). As observed in many types of seizure (Perucca et al., 2013), SLEs are characterized by a beginning (onset), various sequences of fast discharges and spike and wave events (SWEs), and an end (offset). Here, we consider fast discharges and SWEs as basic building blocks of SLEs, i.e. their temporal arrangement, embedding and amplitude can vary indefinitely without affecting the generality of the results. Similar dynamic patterns were found in a zebrafish model of hyperthermia-induced seizure (Fig. 2B) and in human patients with epilepsy (Fig. 2C and Supplementary Table 1).
Time-evolving phenomena can be formalized with differential equations and state variables, which describe the fundamental dynamics of the system. State variables are the smallest possible subset of system variables that can represent the entire state of the system at a given point in time (Supplementary material). They are not unique and often difficult to relate to directly measurable biological/biophysical variables. These variables can be plotted in the space spanned by the state variables (the ‘state space’) to demonstrate the characteristic ‘flow’, which determines how the trajectories evolve. In dynamic system theory, time scale separation allows the grouping of state variables into subsets acting on the same time scale (Haken, 1983) that can be considered as the building blocks of a dynamic system. For SLEs, the two building blocks, fast discharges and SWEs operate on different time scales (fast and slow, respectively). We will refer to each building block as an ‘ensemble’. Fast discharges necessitate a first ensemble with at least two state variables (x1, y1) due to their oscillatory nature (Hale and Koçak, 1991). SWEs comprise large amplitude spikes followed by long lasting wave components reminiscent of a large class of systems collectively termed excitable systems (Gerstner and Kistler, 2002), which necessitate a second ensemble with two state variables (x2, y2). Our goal is now to write the differential equations of ensembles 1 and 2 that account for SLE genesis, time course and recurrence in our experimental conditions.
Epileptor: a dynamic multiscale model of seizures with five state variables
By themselves, the two sets of differential equations for (x1, y1) and (x2, y2) can neither generate the time-evolving event that is a SLE nor its recurrence. At least one additional state variable, z, acting on a very slow time scale (slower than that of SWEs), is necessary to capture the time course of the alternating sequence of SLEs. This slow state variable z guides the entire system, not only between SLEs, but also throughout the SLE time course. At the same time the dynamics of z depends on the states of the other state variables (x1, y1) and (x2, y2). We call z the slow permittivity variable, as it describes the systemic effects that dictate how close the system is to the seizure threshold. As we will show later, z likely includes a large number of extracellular processes that occur on an ultraslow time scale and presumably influence the likelihood of seizure occurrence. In order to build the model, we used the special case where the system is placed in continuous epileptogenic conditions. In subsequent sections, we shall demonstrate that z can in fact cover general cases, including clinically relevant situations.
All three sets of state variables, z, (x1, y1) and (x2, y2), depend upon each other, indicating that the differential equations describing their time evolution are coupled. The analysis of our data set of SLEs revealed that fast discharges are either not present during the SWE (Fig. 2A) or only present during the wave and not the spike (Fig. 2A). This property defines a directional link between (x1, y1) and (x2, y2) that we capture via a bidirectional coupling f1(x1, x2) and f2(x1, x2). It is important to note that fast discharges must not be mistaken for high frequency oscillations, which often override spikes at seizure onset (Bragin et al., 2002). High frequency oscillations are related to epileptogenic networks (Jacobs et al., 2009), but are not completely specific to epilepsy (Blanco et al., 2011) and are not necessarily causally linked to seizure genesis (Quilichini et al., 2012). Rather, the fast discharges described by (x1, y1) are analogous to the low voltage fast activity that is still an active area of basic research (Jiruska et al., 2013) and has great value in localizing focal seizures (Zakaria et al., 2012).
At some distance before and after seizures, the brain appears to operate ‘normally’ and expresses its rich dynamic repertoire of diverse brain states, which may vary greatly in different models and species. The ictal state, however, represents a clear departure from the normal baseline behaviour. The transitions from this normal state to a seizure and back constitute two bifurcations. The invariance of the bifurcations at seizure onset and offset is a manifestation of this nonlinear coupling. By identifying all possible bifurcation combinations we will now develop a taxonomy of SLEs and then discuss the most prominent class.
Taxonomy of seizure-like events
The changes of dynamics at SLE onset and offset have certain characteristic features such as variations of amplitude and frequency of the discharges, which produce distinct flow changes in state space and are the bifurcations we seek to identify. In the mathematical literature, the mode of operation characterized by the alternation of a silent phase of near-equilibrium point activity and an active phase of rapid discharges is called bursting (Ermentrout and Terman, 2010). The first classification of bursters was proposed by Rinzel (1987) and systematically extended, among others, by Izhikevich (2000). Different classes of bursters correspond to the different transitions between the silent and active phase of the burst cycle. It turns out that there are only four types of bifurcations of equilibria (fixed points at seizure onset) and four types of bifurcations of oscillations (planar limit cycles at seizure offset) resulting in 16 classes in total (Izhikevich, 2000; Ermentrout and Terman, 2010). This mathematical classification is the basis of our proposed taxonomy of SLEs. Note that theoretically more complicated behaviours may exist, but are significantly less probable to be encountered in nature (Izhikevich, 2000). The 16 different classes are summarized in Table 1, where we follow the naming scheme of Izhikevich (2000), in which each class is labelled according to the type of bifurcation of equilibria/bifurcation of limit cycle. The four bifurcations of equilibria are saddle-node (fold), saddle-node on invariant circle, supercritical Hopf and subcritical Hopf bifurcation. The four bifurcations of limit cycles are saddle-node on invariant circle, saddle-homoclinic, supercritical Hopf and fold cycle.
|Bifurcations||SNIC||Saddle homoclinic||Supercritical Hopf||Fold limit cycle|
|Supercritical Hopf||Hopf/circle||Hopf/homoclinic||Hopf/Hopf||Hopf/fold cycle|
|Subcritical Hopf||SubHopf/circle||SubHopf/homoclinic||SubHopf/Hopf||SubHopf/fold cycle|
|Bifurcations||SNIC||Saddle homoclinic||Supercritical Hopf||Fold limit cycle|
|Supercritical Hopf||Hopf/circle||Hopf/homoclinic||Hopf/Hopf||Hopf/fold cycle|
|Subcritical Hopf||SubHopf/circle||SubHopf/homoclinic||SubHopf/Hopf||SubHopf/fold cycle|
Based on the classification scheme of codimension-1 planar bursters (Izhikevich, 2000). The rows identify bifurcations of equilibria and the columns bifurcations of oscillations (more specifically, planar limit cycles). SNIC = saddle-node on invariant cycle.
Detailed descriptions of the various bifurcation types and their normal forms (canonical models) can be found in Kuznetsov (1998). Note that normal forms classify local bifurcations around equilibria in state space, whereas global bifurcations involving ‘larger’ invariant sets such as periodic orbits are significantly more complicated. To aid in the classification of seizure types in the framework of a taxonomy of SLEs, the scaling properties of the seizure onset/offset bifurcations are listed as bifurcation of equilibrium in Table 2 and as bifurcation of oscillations in Table 3. As the bifurcation point of seizure onset/offset is approached through the slow change of the permittivity z, the discharges change their behaviour. If denotes the distance to the bifurcation point , then the changes in the discharges will scale with in ways that are characteristic for each bifurcation type. The frequency and amplitude of discharges may be either constant and independent of (denoted by ‘Fixed’ in Tables 2 and 3), completely arbitrary (‘Arbitrary’) or scale from zero following a square root or logarithmic behaviour [‘Zero ()’, ‘Zero ()’]. Each bifurcation may also be distinguished through the number of states that can coexist for the same permittivity value (‘Bistable’ and ‘Monostable’). In the following, we use the taxonomy of SLEs in Table 1 to characterize the time course and recurrence of the SLEs in our experimental conditions.
|Bifurcation of equilibrium||Behaviour||Frequency||Amplitude|
|Supercritical Hopf||Monostable||Fixed||Zero ()|
|Bifurcation of equilibrium||Behaviour||Frequency||Amplitude|
|Supercritical Hopf||Monostable||Fixed||Zero ()|
The four types of bifurcation of equilibrium identify the rows. The column Behaviour indicates if bistability between the oscillatory and equilibrium state is possible, or if the system remains fixed in one state (monostable). The columns Frequency and Amplitude indicate how the onset oscillation scales after the bifurcation of the equilibrium. Here the parameter represents the distance to the bifurcation point. SNIC = saddle-node on invariant cycle.
|Bifurcation of oscillations||Behaviour||Frequency||Amplitude|
|Supercritical Hopf||Monostable||Fixed||Zero ()|
|Fold Limit cycle||Bistable||Fixed||Arbitrary|
|Saddle homoclinic||Bistable||Zero ()||Fixed|
|Bifurcation of oscillations||Behaviour||Frequency||Amplitude|
|Supercritical Hopf||Monostable||Fixed||Zero ()|
|Fold Limit cycle||Bistable||Fixed||Arbitrary|
|Saddle homoclinic||Bistable||Zero ()||Fixed|
The four types of bifurcation of oscillations identify the rows. The column Behaviour indicates if bistability between the oscillatory and equilibrium state is possible. The columns Frequency and Amplitude indicate how the offset oscillation scales before the bifurcation to the equilibrium. Here the parameter represents the distance to the bifurcation point. SNIC = saddle-node on invariant cycle.
Saddle-node bifurcation at seizure onset
In the experiment, the onset of SLEs is characterized by the abrupt appearance of fast discharges (Fig. 2A). As these do not scale up from zero amplitude or frequency, this transition is limited to either a subcritical Hopf or a saddle-node bifurcation (Table 2). A saddle-node bifurcation requires a baseline shift of the measured time-dependent variable, whereas a subcritical Hopf bifurcation results in a transition centered on the baseline and does not have a baseline shift. DC shifts from the baseline field potential occurred systematically in our experimental data (Figs 1, 5B and 7B), ruling out the subcritical Hopf bifurcation.
Homoclinic bifurcation at seizure offset
As SLEs progress towards offset, there is a DC shift back to baseline (Figs 1, 5B and 7B). Only the fold limit cycle and the homoclinic bifurcations show a DC shift at the bifurcation point (indicated by ‘Bistability’ in Table 3), although secondary bifurcations can be introduced to mimick this behaviour (Izhikevich, 2000). The fold limit cycle bifurcation maintains a constant frequency towards the seizure offset, which was not found in the majority of our experiments. The predominant candidate at seizure offset is thus the homoclinic bifurcation, which, amongst others, predicts that the interspike intervals scale logarithmically as seizure offset approaches. The prediction of a logarithmic scaling will be validated experimentally in the following section.
In conjunction, the saddle-node (or fold) bifurcation at seizure onset and the homoclinic bifurcation at seizure offset define the predominant class of an SLE, the fold/homoclinic class, also called square wave burster (Ermentrout and Terman, 2010). Given the constraints imposed by the nature of the bifurcations and the links between the five state variables z, (x1, y1) and (x2, y2), we can now develop the full mathematical model of the SLE, which we call the Epileptor. There are standard models of square wave bursters, which we consider here as a starting point for model development. Typically these models relate to neural discharges on time scales from 10 ms to seconds and are based on some ionic current mechanisms producing slow negative feedback in models of electrical activity in pancreatic beta-cells and respiratory rhythms within the pre-Bötzinger complex (Ermentrout and Terman, 2010). The slow variable of these standard models, however, acts on a different time scale than considered for SLEs, which suggests that the involved biophysical mechanisms may not be the same (see our later discussion on the likely candidate biophysical mechanisms of the permittivity variable). Still, the standard models of square wave bursters provide mathematical guidance in the development of the Epileptor model. For zero coupling between the two ensembles, we adapt the mathematical form of the standard model in Hindmarsh and Rose (1984) for ensemble 1 with (x1,y1) and the mathematical form of an excitable system with a saddle node on invariant circle bifurcation for ensemble 2 with (x2,y2) as shown by Roy et al. (2011). Then the following modifications are introduced: a linear inhibitory coupling from ensemble 2 to 1 and a low-pass filtered excitatory coupling from ensemble 1 to 2 to generate the SWE and interictal spikes; the negative feedback coupling of the permittivity z to the ensemble 2 to bias the preictal spikes towards SLE onset; and changes of the non-linearities of the original standard models to guarantee the structural stability of the bifurcations. Structural stability was tested computationally for all used parameters. The final parameter values were chosen to fit the Epileptor against the experimental data, where the sum of the two ensemble x-variables, x1 + x2, was matched visually against the electrographic signatures of a SLE. Although each of the state variables may reflect a diversity of biophysical variables, we found that plotting x1 + x2 as a function of time bore striking resemblance with the field potential. The complete system of the Epileptor equations reads then as follows:Supplementary material). The initial conditions for the numerical simulation are x1 = 0; y1 = −5; z = 3; x2 = 0; y2 = 0. Noise is introduced into each equation as linear additive Gaussian white noise with zero mean and a variance of 0.025 for the first subsystem and 0.25 for the second subsystem. For the solution of the stochastic equations we employ the Euler-Maruyama method.
Dynamics of the Epileptor
In Epileptor, ‘normal’ brain function and seizures occupy two different regions of state space, being isolated from one another by a divergent flow acting as a barrier (Fig. 3). This barrier between the two states establishes bistability and is known as the separatrix, which in this case essentially describes the seizure threshold (Frohlich et al., 2010; Kramer et al., 2012). We show the bifurcations of seizure onset (Fig. 3, rows Ia and Ib) and seizure offset (Fig. 3, rows IIa and IIb) as cartoons of six potential landscapes (rows Ia and IIa) and six vector fields (rows Ib and IIb). The potential landscapes metaphorically illustrate the transitions between seizure and non-seizure states for changing values of permittivity z and the vector fields show the corresponding flows in the two-dimensional state space of the Epileptor’s ensemble 1. In rows Ia and IIa, two minima are separated by a potential well. The left minimum is the normal state, the right minimum is the seizure state. As the permittivity variable changes, the potential landscape changes and one minimum becomes a maximum resulting in a transition from one state to the other. The flows in state space show the same behaviour (Ib and IIb). Here we plot trajectories for various initial conditions, where the arrows indicate the direction of the flow and circles indicate the equilibrium points (fixed points). A stable fixed point attracts trajectories in its neighborhood (full circle), whereas unstable fixed points deflect the trajectories (empty circle). A so-called saddle point is a fixed point with a stable and unstable direction (empty circle). In the first figure of Ib (from left to right) a saddle point separates the state space in two regions with a stable fixed point to its left and a stable limit cycle to its right. As the permittivity is changed towards seizure onset the separatrix moves closer to the stable fixed point (middle figure, row Ib) until the separatrix collides with the stable fixed point and the two disappear (right figure, row Ib) via a saddle-node bifurcation. Row IIb shows the corresponding scenario for seizure offset: as the permittivity changes, the separatrix moves towards the limit cycle (middle figure) until separatrix and limit cycle collide and annihilate each other (right figure). This event is called the homoclinic bifurcation and leaves the system with only the stable fixed point. A bifurcation diagram captures the hitherto described dynamics more quantitatively in Fig. 4, which has been generated using analytical calculation of fixed points and numerical continuation. Fig. 5A traces out time series from simulations of the Epileptor model. Note the presence of SWEs before seizure onset, even far from it. The presence of more interictal spikes (and even fast oscillations) reflects that the Epileptor is getting close to the separatrix. The experimental/clinical analogy would be that network excitability increases and generates spikes and/or oscillations. These activities have their own dynamics; they can take various shapes (reflecting the diversity of interictal/preictal states) as the system moves close or away from the separatrix. This concept of getting close to the separatrix and away from it without reaching the bifurcation point is analogous to the ‘preictal state’ that has been theorized for many years (Stacey et al., 2011a) and may have been identified in a recent clinical trial (Cook et al., 2013).
Because the Epileptor was constructed on the basis of an experimental model of SLEs (an immature hippocampus placed in continuous epileptogenic conditions in vitro), we first tested the generality of its predictions in other brain regions and species, and then explored the concept of the separatrix to understand how seizure onset can be reached in more realistic conditions.
Validating model predictions: bifurcations and invariance
The Epileptor model makes two important predictions about seizure onset and offset in the seizure class fold/homoclinic (Table 1) that we can validate directly: (i) seizure onset can only occur in the presence of a DC shift of the field potential; and (ii) the interspike intervals show a logarithmic scaling approaching seizure offset.
Presence of a direct current shift at seizure onset
We used the presence of a DC shift at seizure onset and its reversal after seizure offset in our experimental model to identify the onset and offset bifurcations (Fig. 1). This DC shift is smaller (−0.26 ± 0.08 mV, n = 15 SLEs from five hippocampi, P < 0.01) and shorter lasting (1–3 min) than that found in spreading depression (Somjen, 2001). Like most clinical and experimental EEG, seizures are recorded in alternating current mode, thus filtering any slow variations of the field potential (Fig. 2). The DC shift at onset of a partial spontaneous seizure, is actually well known in several species including baboons (Pumain et al., 1985) and humans (Ikeda et al., 1999; Vanhatalo et al., 2003), though it is under-recognized clinically because clinical electrodes are poor at recording DC (Tallgren et al., 2005; Stacey et al., 2012) and seizures are rarely recorded with DC coupling. To the best of our knowledge, we are not aware of published spontaneous focal seizures recorded in DC mode without a DC shift.
Logarithmic scaling of interspike intervals approaching seizure offset
The presence of a homoclinic bifurcation at seizure offset imposes a logarithmic scaling of interspike intervals. Introducing noise in Epileptor enabled us to generate numerous SLEs and verify that SLE offsets in the Epileptor show the logarithmic scaling of homoclinic bifurcations (Supplementary Fig. 1). In all SLEs recorded in our experimental conditions in vitro (n = 16 hippocampi), the interspike intervals exhibited a logarithmic scaling at seizure offset (Fig. 6A and B), verifying the prediction that seizure offset corresponds to a homoclinic bifurcation. We then tested the interspike intervals in each of the first 24 patients with identifiable seizures in the International EEG database (www.ieeg.org). This arbitrary sample includes several different seizure types and brain regions (Supplementary Table 1). We verified logarithmic scaling in 20 (83%) of the patients (Supplementary material), the majority of which had remarkably similar dynamics despite their clinical disparities (Fig. 6B, Supplementary Fig. 2 and Supplementary Table 1). Similar findings, which appear visually as spikes slowing down near the end of seizures, were recently identified in live human recordings on multiple spatial scales from EEG down to multiunits (Kramer et al., 2012). The same property was found in other species, such as zebrafish (Fig. 6B and C, Supplementary Fig. 9), and is apparent also in flies [see Fig. 2C in Zhang et al. (2002)]. These invariant scaling laws argue strongly in favour that a homoclinic bifurcation at offset is an invariant property of focal seizures.
The Epileptor model allows us to unravel the topology and trajectories of seizures in the state space defined by the five state variables. Since five-dimensional representations are not practicable, we used a tri-dimensional representation for illustration, plotting x1, x2 and z as a function of time. The resulting topology corresponds to spirals on a cone (Fig. 5A). Although −x1+x2 bears analogy with the field potential, the precise biophysical equivalent of the z variable is unknown and will be likely complex. One of its distinguishing features is its slow time-dependent evolution. Several biophysical parameters are known to evolve slowly in time during seizures, including extracellular ions (Heinemann et al., 1986) and oxygen (Suh et al., 2006). We measured extracellular [K+] during SLEs and plotted it as the z variable; the topology of experimental SLEs were remarkably similar to the theoretical one (Fig. 5B). As the measurement of a slow variable was not available in the humans and zebrafish, we constructed the trajectories in a space of delayed state variables, which approximates the original state space (Takens, 1981) and allows unfolding the resulting trajectories for illustration of the flow topology. Even though the trajectory lines are distorted and difficult to read, they still unveil the basic phases of the seizures, showing clear similarities between the Epileptor and SLEs (Supplementary Fig. 3).
Validating model predictions: how seizures start
Seizures can occur with regularity in some patients, as for catamenial epilepsy (Penovich and Helmers, 2008) or during the sleep-wake cycle (Karafin et al., 2010), in keeping with the possibility that ultra-slow systemic effects, analogous to the slow permittivity variable, can push the system periodically to seizure onset. However, seizures arise, most of the time, without obvious causes or predictive factors. In such instances, we propose that brain trajectories are close to the separatrix, and the closer the ‘normal’ brain state is to the separatrix, the easier it will be to have a seizure. Linked to the distance to the separatrix (i.e. proximity to seizure threshold) is the build-up of preictal spiking in the Epileptor (Fig. 7C), which is analogous to physiological fluctuations around baseline such as increased preictal spiking (Huberfeld et al., 2011), other preictal waveforms (Stacey et al., 2011a) or predictive features in the EEG (Cook et al., 2013). However, these fluctuations are quite variable, do not inevitably progress to a seizure (Cook et al., 2013), and are still dominated by the normal baseline activity. Thus, these excursions have not yet crossed the separatrix and provide evidence that (i) the system is close to threshold; and (ii) there are many potential trajectories from the ‘normal’ state to the ‘seizure’ state. The stereotypical appearance of a seizure despite the numerous trajectories leading into it provides further evidence that it is a distinct state of the system. In the Epileptor, it is possible to know the distance from the separatrix (bifurcation diagram in Fig. 4), and force the system pass its threshold. We shall explore two possibilities for provoking seizures and validate them experimentally: large external stimulation, as well as timely internal noise.
Triggering seizures with electrical stimulation
The most straightforward method to produce an excursion of brain trajectories toward the separatrix is to manipulate the system directly. In Epileptor, this can be simply done by adding current, which moves the system towards the separatrix and triggers a SLE before its expected time of occurrence (Fig. 7A, traces d and e). Electrical stimulations consistently trigger seizures in humans: following electroconvulsive shocks in any ‘normal’ brain (Walker, 2011) and following cortical stimulation during presurgical evaluation of epileptic patients (Valentin et al., 2005). Another prediction, which is not as well described clinically, is that there is a refractory period after an SLE offset in which the same ‘kick’ cannot produce a seizure (Fig. 7A, traces b and c). The system will be refractory outside of the interval of bistability (Fig. 4). To test both predictions experimentally, we isolated the whole mouse hippocampus and the septum connected to it (Fig. 7B). We placed the hippocampus in one chamber and the septum in another, but maintained the physical connection between them. Each chamber could be perfused with different solutions without exchange between both compartments (Supplementary Fig. 4). Bathing the hippocampus in low Mg2+ artificial CSF resulted in recurrent SLEs with a typical DC shift (Fig. 7B). The concentration of Mg2+ was then raised (0.4–0.6 mM) to prevent the occurrence of spontaneous SLEs. Hippocampal circuits were thus maintained close to the separatrix, i.e. the z variable does not drive the system to the separatrix but maintains it in its vicinity. A train of electrical stimuli (10 s, 10 Hz, 170–230 µA) applied to the septum triggered a SLE in the hippocampus (Fig. 7B). After SLE offset, the same trains failed to elicit SLEs until at least 10 min had passed (n = 3 independent experiments), demonstrating a refractory period.
Synaptic noise as a physiological seizure trigger
In non-linear dynamics, noise enables a system to explore its state space (Deco et al., 2011) and in computational models has been shown to initiate and spread epileptiform activity (Suffczynski et al., 2006; Stacey et al., 2011b). Introducing noise in the Epileptor gave rise to different SLE patterns, still keeping the generic features of fast discharges and SWEs as basic building blocks of activity, but organized differently (Supplementary Fig. 5). This shows that seizures can be organized differently in terms of discharge patterns, while keeping the same invariant features. Next, we varied the distance to the saddle-node bifurcation systematically in the Epileptor and investigated its sensitivity to noise (Fig. 7C). Introducing noise made the Epileptor generate interictal spikes and reach seizure onset before expected from the bifurcation diagram, thus generating a SLE. In other words, although the system was at some distance from the bifurcation point, while both normal and seizure states coexisted, noise made the system explore its surroundings in state space until it passed the separatrix. The number of interictal spikes was directly related to the distance to the bifurcation (Fig. 7C) reflecting the slow transformations occurring in network dynamics as the system approaches the seizure (Trevelyan et al., 2007).
Noise was also able to evoke SLEs in the experimental preparation. We found a systematic increase of spontaneous synaptic noise received by neurons before SLE onset when hippocampi were placed in continuous epileptogenic conditions (Supplementary Fig. 6). To demonstrate causality, we used the dual-chamber as above, with the hippocampus maintained in conditions that did not enable the genesis of spontaneous SLEs. Raising external [K+] in the septum chamber increased the firing activity of septal neurons, thus sending more synaptic activity in their target hippocampal cells (Supplementary Fig. 7). The septum could thus be used as a generator sending synaptic noise to the hippocampus. When reaching a critical value, synaptic noise was sufficient to trigger SLEs in the hippocampus (Fig. 7D and Supplementary Table 2).
Increasing synaptic noise is not the only way to pass the threshold. The conceptual framework provided by Epileptor enables exploring other putative mechanisms. For example, change in osmolarity can occur in several clinical situations and is associated with increased seizure susceptibility, as it alters synaptic transmission, cell volume, and ephaptic communication (Andrew, 1991). Using hippocampi maintained in subthreshold conditions, we found that changing the osmolarity in the hippocampal chamber also generated SLEs (Fig. 7D and Supplementary Table 2). Importantly, there was no noise increase before SLEs, demonstrating that the seizure onset was approached via a different route. However, changes in noise and osmolarity were synergistic, producing SLEs when subthreshold levels of each were applied simultaneously (Fig. 7D). In such a two-state non-epileptic/epileptic system, there are thus multiple ways or combinations to reach seizure onset, such as, but not limited to, stimulation, noise increase and change in osmolarity. Clinically, this prediction means that there are many different paths to reach seizure onset.
Biophysical correlates of Epileptor state variables
The state variables in the Epileptor, although abstract, are the key to quantifying seizure dynamics. Recognizing biological processes that have similar behaviours is technically very challenging, but helps understand the underlying processes that comprise a seizure. We now present a strategy to identify these processes experimentally, acknowledging that their identification may only apply to the intact immature mouse hippocampus in low Mg2+.
Slow permittivity variable
In Epileptor, the key property of the permittivity variable z is its evolution on a slow time scale. Several biophysical parameters are known to change slowly during or preceding seizures, including extracellular levels of ions (Heinemann et al., 1986), oxygen (Suh et al., 2006) and metabolism (Zhao et al., 2011). We thus measured the levels of extracellular [K+] (Bazhenov et al., 2008; Frohlich et al., 2008), oxygen and intracellular NADH/FAD+, which reflects the activity of energy metabolism in hippocampi placed in epileptogenic conditions. The return to baseline of [K+], pO2, and NADH/FAD+ corresponded to the initiation of the next bifurcation, as predicted by the Epileptor (Fig. 5B and Supplementary Fig. 8). We did not identify a biophysical variable changing during the interictal period, i.e. a variable that would drive the system to seizure threshold. However, [K+], pO2, and NADH/FAD+ appear to contribute to SLE time course and offset (but not its initiation), with a time evolution compatible with z during SLEs. It is interesting to note that seizures recorded in vivo in cats (Moody et al., 1974) and awake baboons (Pumain et al., 1985) are characterized by the same time-dependent evolution of extracellular [K+].
Cells recapitulate state variable behaviour
The Epileptor predicts that the fast subsystem composed of ensemble 1 with variables (x1, y1) should be active during SWEs with the fast oscillations occurring only during the wave part of the event. As experimental data suggests that glutamatergic and GABAergic cells are important for fast discharges and SWEs, respectively (Isomura et al., 2008), we predicted that GABAergic cells would be more active during SWEs than during fast discharges, whilst glutamatergic cells would display a reverse pattern of activity. Hence we associate the involvement of glutamatergic and GABAergic activity with the faster variable x1 and slower variable x2, respectively. These hypotheses were tested experimentally in hippocampi placed in continuous epileptogenic conditions. Consistent with the prediction, GABAergic neurons recorded in the cell attached configuration fired action potentials during SWEs, stopped firing during the fast discharge, and resumed firing when SWEs re-occurred during SLEs (Fig. 8A). Current clamp recordings revealed that they stopped firing during the fast discharge because they entered into depolarization block (Fig. 8B). Single cell recordings provide only a partial picture of the activity of GABA neurons. If GABA neurons are more active during SWEs than fast discharges, we predicted that the GABAergic synaptic drive (produced by GABA neuron firing) received by hippocampal neurons should be predominant during SWEs. In keeping with this prediction, whole cell recordings revealed that hippocampal neurons received a strong barrage of synaptic GABAergic currents before SLE onset, which was largely reduced during the fast discharge (Fig. 8C). The GABAergic drive recovered upon recurrence of SWEs, only during the spike component, but not during the fast oscillation embedded in the wave as predicted. Conversely, hippocampal neurons received a barrage of high frequency glutamatergic inputs synchronized with the fast discharge in the field (Fig. 8D), in keeping with the predicted behaviour of pyramidal cells. It is important to note that the presence of strong glutamatergic currents during the spike component of SWEs indicates that pyramidal cells also contribute to (x2, y2), which is why the biophysical correlate of the Epileptor is more complex than the simple interpretation of ensembles 1 and 2 as excitatory and inhibitory cells.
The previous results show that it is possible to construct hypothesis-driven search of the biophysical variables contributing to the state variable, which may be specific to experimental conditions (as demonstrated hereafter), and in the clinic, to any specific patient.
Seizures: conserved dynamics from many trajectories
It is important to stress here again that the Epileptor model does not impose constraints upon the nature of the biological variables; only upon their dynamics and respective relationships. This means that multiple parameter configurations can theoretically give rise to the same system behaviours. For instance, work in the guinea pig brain suggests that the fast discharges therein are produced by inhibitory interneurons, rather than by pyramidal cells (Gnatkovsky et al., 2008). This concept, that there are many functional biological pathways that produce the same outcome, was developed and verified experimentally in simple neuronal networks (Marder and Taylor, 2011) and is present in human biology as well (Beall, 2007). This issue is particularly relevant both clinically and experimentally, as numerous anatomical, molecular, electrophysiological and functional modifications have been described in different types of epilepsies, experimental models and species; without clear consensus emerging about major contributing parameters. Hence the possible correlates of the five state variables in the low Mg2+ model of SLEs in the intact hippocampus may be only valid for these very specific experimental conditions. To illustrate this key concept, we used a drastic situation, analysing SLEs recorded in the absence of Ca2+, a condition that abolishes neurotransmitter release (Jiruska et al., 2010). In these conditions, the scheme proposed above with pyramidal cells/GABA neurons and intact neurotransmission cannot apply. Yet, SLEs, in low Ca2+, were also characterized by a strong DC shift, and an offset still characterized by a homoclinic bifurcation with logarithmic scaling (Supplementary Fig. 9). Interestingly, seizures recorded in awake baboons are similarly characterized by a large decrease in extracellular Ca2+ compatible with a lack of neurotransmission (Pumain et al., 1985). Hence, even in the absence of neurotransmission, SLEs still maintain their invariant features, further supporting the generic nature of our findings.
We have identified several invariant features of seizures that are preserved across different species, models, and brain regions. Invariant features serve as constraints to define landmarks in high-dimensional parameter spaces, allowing the researcher to describe inherent features that are independent of specific parameters. Here we used the scaling behaviour of frequency and amplitude during seizure onset and offset to identify the specific bifurcations underlying these transitions between brain states. This approach allows a systematic characterization of SLEs into a taxonomy of 16 classes. Based on the analyses of our experimental data, we identified the predominant class of SLEs as the seizure onset/offset pair ‘fold/homoclinic bifurcation’ and modeled its dynamics through the Epileptor model.
Where does the Epileptor stand in the model literature? Models can be broadly separated in ‘analogies’ and ‘biophysical’. Scaling laws are often used to characterize system properties as a whole by making analogies to models from physics, such as coupled mechanical pendula (Osorio et al., 2010) or sand piles (Bak and Paczuski, 1995). Such analogies allow for some intuition of a functional mechanism and demonstrate the existence of certain general system properties, including multistability (Lopes da Silva et al., 1994, 2003a; Shu et al., 2003; Teramae and Fukai, 2007; Milton, 2010), bifurcations (Lopes da Silva et al., 2003b), and delayed recurrent loops (Foss et al., 1996; Foss and Milton, 2000). But analogies do not predict which multistable states or which bifurcations are present in a given neuroscience system. More importantly, they do not generate a time series that specifically models the data. On the other extreme, biophysical models are derived from commonly accepted microscopic neuron models, such as the Hodgkin-Huxley equations and/or biophysical principles and laws, e.g. Nernst equation of electrochemical equilibrium, conservation of mass and/or energy, or mean fields (Deco et al., 2008), leading to neural population or ‘mass’ models (Wilson and Cowan, 1972; Nunez, 1974; Freeman, 1975; Jansen and Rit, 1995; Brunel and Wang, 2003; Stefanescu and Jirsa, 2008). Many of these models find applications in epilepsy (Wendling et al., 2002; Breakspear et al., 2006; Kramer et al., 2012; Wang et al., 2012). In these models, parameters have a biophysical meaning and can often be independently measured, allowing for experimental validation. But the biophysical parameter space is vast and many parameter configurations can give rise to the same system behaviours (Marder and Taylor, 2011). Our approach has been to identify invariant landmarks in the parameter space obtained from the experimental data. These landmarks are the bifurcations that unambiguously define how a system qualitatively changes its behaviour. The taxonomy of SLEs assembles these invariant dynamic elements of seizure evolution in one framework that can be used to guide interpretation of experimental results or as a map for choosing parameters in biophysical models. It is remarkable that the vast majority of seizures in patients corresponded to the class of the bifurcation pair saddle-node/homoclinic. Though our findings are consistent across all data sets, species and brain regions we studied, other seizure types will likely exist. Indeed, seizures in four patients did not display homoclinic bifurcations at seizure offset because their interspike intervals remained essentially constant throughout the seizure (Supplementary Table 1). This behaviour is unusual for seizures, but is potentially consistent with the fold/Hopf and fold/fold classes (assuming the saddle-node (fold) bifurcation holds for seizure onset) in the SLE taxonomy (Table 1). Another distinction is made regarding absence seizures. Previous mathematical models of absence epilepsy have predicted seizure onsets in which there is no warning or ‘slow process’ (Lopes da Silva et al., 2003a; Breakspear et al., 2006) and the ictal state is entered through a Hopf bifurcation with a continuously evolving poly-spike pattern (Marten et al., 2009; Rodrigues et al., 2009). Absence seizures are characterized by spike and wave discharges resembling the time scales observed for SWEs of the Epileptor, which are slow compared to the fast discharges in non-absence seizures.
One of the most intriguing predictions is that seizure onset, time course and offset are controlled by a slow permittivity variable. This concept, that focal seizures inherently are influenced by slow processes that govern when they are likely to occur, emphasizes the important role of extracellular effects for discharges of neuronal populations. In our experimental model, the levels of extracellular [K+], oxygen and ATP consumption show the same time-dependent changes as the permittivity variable during seizures. Such correlation does not imply causality, and many other parameters may also display slowly evolving modifications, e.g. release of molecules, synaptic vesicle depletion, phosphorylation processes etc. The varying release of neuromodulators-neurotransmitters during biological cycles (circadian, sleep/wake), alterations in glucose/O2 supply etc. could fit the dynamic constraints of this slow permittivity variable. They may constitute key biomarkers for seizure prediction and may be investigated as such. As some of these parameters may not give rise to an electrophysiological signature, specific sensors (e.g. glucose, ATP and adenosine) would be needed to measure them. It is important to note that if some variables may push the system toward the bifurcation, many others oppose this movement (collectively referred to as protective/anti-seizure mechanisms), e.g. peptides (such as vasoactive intestinal polypeptide), activation of adenosine A1 receptors etc. From the Epileptor standpoint, the balance between slowly acting pro- and anti-seizure mechanisms are exactly what constitutes the slow permittivity variable. A challenge for future development will be the bridging of different spatiotemporal scales, such as the levels of single neuron and neural population dynamics using mean field techniques (Deco et al., 2008). Given the heterogeneity of neuronal behaviour during seizure initiation (Truccolo et al., 2011), we anticipate a complex interplay among groups of neurons that present different types of spiking patterns, probably with fast and slower time scales as encountered in the two Epileptor ensembles. Furthermore, environmental effects and electrotonic couplings will likely play a prominent role in light of our results in absence of synaptic transmission (see Cressman et al., 2009 for biophysical candidate mechanisms). Understanding which sets of parameters control seizure time course (for example with optogenetics to control the activity of specific cell types) will be crucial to design the best strategies to stop focal-onset seizures as soon as they start with closed-loop systems (Krook-Magnuson et al., 2013; Paz et al., 2013).
The concept of a slow permittivity variable may provide alternate (but non-exclusive) explanations to long-standing debates. For example, whether or not interictal spikes are causally related to seizure genesis remains unclear. A study using tissue slices obtained from epileptic patients showed that a build-up of large pre-ictal spikes preceded seizure-like events when slices were placed in continuous epileptogenic conditions (Huberfeld et al., 2011). The transition was slow (30 min) and required the activation of NMDA receptors. Rather than driving seizures, we propose that the occurrence of pre-ictal spikes may just reflect the modifications occurring within the network on a slow time scale. In other words, as the slow permittivity drives the system close to the bifurcation, the conditions for generating such pre-ictal spikes are met. Hence, interictal spikes may just signify a specific position of the system in its state space. In keeping with this view, interictal spikes show complex dynamics in the days preceding the first spontaneous seizure when networks undergo complex reorganizations in vivo (El-Hassar et al., 2007; Chauviere et al., 2012), and they tend to disappear over time when slices are bathed in continuous epileptogenic conditions (Trevelyan et al., 2007). This proposal does not rule out the possibility for interictal spikes to act as a driving force toward the bifurcation, but is not part of the Epileptor mechanisms.
Nearly every brain region can be driven out of the ‘healthy’ subspace to produce seizures, depending upon the severity and the widespread diffusion of the process leading to the seizure. In a ‘healthy’ brain, the trajectories of brain activities are far from the seizure threshold, and need strong external interventions (like an electroconvulsive shock) to reach seizure state. In pathological conditions, we propose that the reorganization of the underlying circuits move normal brain trajectories closer to the separatrix (which corresponds to how much the seizure threshold has been lowered), increasing the likelihood for seizure occurrence. Many factors can potentially move the system in such a way, which would explain why network reorganizations are often brain region-, model-, time- and epilepsy type-dependent (Pitkanen and Sutula, 2002). Hence, there are a large number of possible combinations, all leading to the same functional outcome: bringing normal brain trajectories in the vicinity of the separatrix. A recent clinical trial, in which patients with epilepsy were recorded continuously over several months, found potential evidence of brain approaching the separatrix (Cook et al., 2013). In this study, several patients had unique, characteristic changes in their EEG before seizures that could be identified by an automated algorithm. In these patients, the algorithm was quite sensitive in predicting oncoming seizures; however, there were many ‘false alarms’ in which a seizure did not occur. The Epileptor predicts that these changes may not have been failures of the algorithm, but rather a reflection of the system approaching the bifurcation.
We performed the experimental analysis of seizure dynamics on data from single regions, mostly the hippocampus. However, seizures often involve large networks of networks, with initiation and propagation zones (Bartolomei et al., 2008). Epileptor may serve as a building block of coupled dynamical models to study the network mechanisms of seizure propagation over larger brain regions. Because seizure propagation is the most detrimental factor to the patient’s quality of life, the identification of potential invariances in seizure propagation would be particularly beneficial in the clinic.
We conclude that seizures belong to the possible repertoire of brain activities. They can occur under stringent conditions in the ‘normal’ brain, but their probability of occurrence is increased as the underlying reorganizations bring the system close to the bifurcation. This may explain why so many different pathological conditions (e.g. Alzheimer’s disease, stroke, autism, brain trauma etc.) are also associated with seizures, yet seizures generally have conserved dynamics that are recognizable to clinicians regardless of pathology. Their invariance from flies to rodents to humans clearly argues that seizures are a universal behaviour of neural systems.
We thank P. Jiruska, J.R. Jefferys and S. Baraban for providing the data on zero Ca2+ and zebrafish seizures and K. El Houssaini for help in preparing Fig. 4. We thank for carefully reading the manuscript: A. Gulyas, G. Kalamangalam, W. Loescher, R. McIntosh, P. Ritter, J. Terry and D. Turner.
Supported by INSERM, ANR MOTION, Fondation pour la Recherche sur l’Epilepsie (FFRE), Fondation pour la Recherche sur le Cerveau (FRC), CURE, the Brain Network Recovery Group through the James S. McDonnell Foundation, the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement #602102, and the FP7-ICT BrainScales, and NIH K08NS069783.
Supplementary material is available at Brain online.
spike and wave event