Abstract

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.

Introduction

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

Experimental procedures

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%].

Data analysis

Noise

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.

Raw data

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.

Statistics

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 recordings

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.

Zebrafish recordings

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).

Results

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).

Figure 1

Experimental model of spontaneous SLEs. (A) When placed in persistent epileptogenic conditions, SLEs are generated at regular intervals in the isolated mouse whole hippocampus at postnatal Day 6 in vitro. Direct current recordings show a DC shift at SLE onset, which reverses after SLE offset. (B–D) Correspond to insets b–d in A.

Figure 1

Experimental model of spontaneous SLEs. (A) When placed in persistent epileptogenic conditions, SLEs are generated at regular intervals in the isolated mouse whole hippocampus at postnatal Day 6 in vitro. Direct current recordings show a DC shift at SLE onset, which reverses after SLE offset. (B–D) Correspond to insets b–d in A.

Figure 2

Seizure patterns conserved across species and brain regions. All displayed recordings were obtained in alternating current mode, thus filtering very slow variations of the field potential. (A) SLE recorded in mouse whole hippocampus displaying a typical sequence of tonic and tonic-clonic patterns. Two generic patterns can be distinguished: fast discharges (#) and large spike and wave events (*). Note that fast discharges can be embedded in the wave (*#). (B) Hyperthermia-induced SLE recorded in vivo in zebrafish display similar patterns, with fast discharges shown in panel I and SWE in panel II. (C) Spontaneous seizure recorded in an epileptic patient (Supplementary Table 1) displaying a fast discharge followed by the occurrence of spike and wave events. Note that in all three species, there is a slowing down of the activity when reaching seizure offset.

Figure 2

Seizure patterns conserved across species and brain regions. All displayed recordings were obtained in alternating current mode, thus filtering very slow variations of the field potential. (A) SLE recorded in mouse whole hippocampus displaying a typical sequence of tonic and tonic-clonic patterns. Two generic patterns can be distinguished: fast discharges (#) and large spike and wave events (*). Note that fast discharges can be embedded in the wave (*#). (B) Hyperthermia-induced SLE recorded in vivo in zebrafish display similar patterns, with fast discharges shown in panel I and SWE in panel II. (C) Spontaneous seizure recorded in an epileptic patient (Supplementary Table 1) displaying a fast discharge followed by the occurrence of spike and wave events. Note that in all three species, there is a slowing down of the activity when reaching seizure offset.

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.

Table 1

A taxonomy of seizure-like events

Bifurcations SNIC Saddle homoclinic Supercritical Hopf Fold limit cycle 
Saddle-node Fold/circle Fold/homoclinic Fold/Hopf Fold/fold cycle 
SNIC Circle/circle Circle/homoclinic Circle/Hopf Circle/fold 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 
Saddle-node Fold/circle Fold/homoclinic Fold/Hopf Fold/fold cycle 
SNIC Circle/circle Circle/homoclinic Circle/Hopf Circle/fold 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 forumla denotes the distance to the bifurcation point forumla, then the changes in the discharges will scale with forumla in ways that are characteristic for each bifurcation type. The frequency and amplitude of discharges may be either constant and independent of forumla (denoted by ‘Fixed’ in Tables 2 and 3), completely arbitrary (‘Arbitrary’) or scale from zero following a square root or logarithmic behaviour [‘Zero (forumla)’, ‘Zero (forumla)’]. 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.

Table 2

Properties of oscillatory onset bifurcations

Bifurcation of equilibrium  Behaviour  Frequency  Amplitude 
Saddle-node  Bistable  Fixed  Fixed 
SNIC  Monostable  Zero (forumla Fixed 
Supercritical Hopf  Monostable  Fixed  Zero (forumla
Subcritical Hopf  Bistable  Fixed  Arbitrary 
Bifurcation of equilibrium  Behaviour  Frequency  Amplitude 
Saddle-node  Bistable  Fixed  Fixed 
SNIC  Monostable  Zero (forumla Fixed 
Supercritical Hopf  Monostable  Fixed  Zero (forumla
Subcritical Hopf  Bistable  Fixed  Arbitrary 

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 forumla represents the distance to the bifurcation point. SNIC = saddle-node on invariant cycle.

Table 3

Properties of oscillatory offset bifurcations

Bifurcation of oscillations  Behaviour  Frequency  Amplitude 
SNIC  Monostable  Zero (forumla Fixed 
Supercritical Hopf  Monostable  Fixed  Zero (forumla
Fold Limit cycle  Bistable  Fixed  Arbitrary 
Saddle homoclinic  Bistable  Zero (forumla Fixed 
Bifurcation of oscillations  Behaviour  Frequency  Amplitude 
SNIC  Monostable  Zero (forumla Fixed 
Supercritical Hopf  Monostable  Fixed  Zero (forumla
Fold Limit cycle  Bistable  Fixed  Arbitrary 
Saddle homoclinic  Bistable  Zero (forumla 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 forumla 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.

The Epileptor

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:  

formula
where  
formula
and x0 = −1.6; y0 = 1; τ0 = 2857; τ1 = 1; τ2 = 10; Irest1 = 3.1; Irest2 = 0.45; γ = 0.01. Here the state variables x1 and y1 comprise the first subsystem responsible for fast oscillations and x2 and y2 the second subsystem involved in spike wave events. The slow permittivity variable is z. The characteristic time scales are τ0 of the permittivity variable, τ1 of ensemble 1 and τ2 of ensemble 2, where the time scale hierarchy is τ0 ≫ τ2 ≫ τ1. The time constant of τ1 does not appear in the equations explicitly as it is equal to 1, hence we omitted it for simplicity. Note that the integral coupling function g(x1) can be rewritten as an ordinary differential equation, which then technically introduces a sixth state variable (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).

Figure 3

Caricatures of the flows in state space of ensemble 1 as the slow permittivity variable z changes. Rows Ia and IIa indicate metaphorically the bistability of the ‘normal’ (left minimum, Ia) and seizure (right minimum, IIa) state, as well as its loss. Note that ‘normal’ brain trajectories are displayed as a fixed point for the sake of illustration (it does not reflect the diversity of possible trajectories). Rows Ib and IIb show the corresponding flows in state space. As the permittivity z decreases (from left to right), rows Ia and Ib show how the interictal state loses its stability and the transition occurs towards the ictal state (seizure onset) via a saddle-node bifurcation. Rows IIa and IIb show the equivalent situation for increasing values of z and the homoclinic bifurcation leading to seizure offset.

Figure 3

Caricatures of the flows in state space of ensemble 1 as the slow permittivity variable z changes. Rows Ia and IIa indicate metaphorically the bistability of the ‘normal’ (left minimum, Ia) and seizure (right minimum, IIa) state, as well as its loss. Note that ‘normal’ brain trajectories are displayed as a fixed point for the sake of illustration (it does not reflect the diversity of possible trajectories). Rows Ib and IIb show the corresponding flows in state space. As the permittivity z decreases (from left to right), rows Ia and Ib show how the interictal state loses its stability and the transition occurs towards the ictal state (seizure onset) via a saddle-node bifurcation. Rows IIa and IIb show the equivalent situation for increasing values of z and the homoclinic bifurcation leading to seizure offset.

Figure 4

Bifurcation diagram of the Epileptor. (A) The set of fixed points form curves, where the solid line indicates the stable fixed point. A branch of limit cycles terminates at the homoclinic bifurcation point (HB), whereas the fixed points lose stability via saddle-node bifurcations (SN). The system displays bistability between the left saddle-node bifurcation point (SN) and the homoclinic bifurcation point (HB). (B) The projection of the Epileptor trajectory is plotted onto the bifurcation diagram.

Figure 4

Bifurcation diagram of the Epileptor. (A) The set of fixed points form curves, where the solid line indicates the stable fixed point. A branch of limit cycles terminates at the homoclinic bifurcation point (HB), whereas the fixed points lose stability via saddle-node bifurcations (SN). The system displays bistability between the left saddle-node bifurcation point (SN) and the homoclinic bifurcation point (HB). (B) The projection of the Epileptor trajectory is plotted onto the bifurcation diagram.

Figure 5

Slow permittivity state variable and seizure topology. (A, left) Seizure generated by the Epileptor with five state variables. Seizure onset, time course and offset are controlled by the permittivity state variable evolving slowly in time (red). Note that the SLE occurs with a rapid and large shift of the potential. Right, the seizure trajectory (expressed in terms of −x1 + x2) is approximated in a 3D space defined by the first state variables (X = −x1 and Y = x2) and by the slow permittivity variable (Z = z). Note that the values of the z-variable have been shifted upwards for plotting purposes. (B, left), simultaneous recording in the hippocampus of a SLE in low Mg2+ conditions in DC mode, O2 levels in the preparation (yellow) and NADH levels (red), which indirectly reflect ATP use. Note the large DC shift during the SLE, as predicted by Epileptor. The time course of oxygen and NADH is similar to that of the slow permittivity variable. Right, the 3D representation of a seizure in a delayed space (X and Y), with Z the extracellular potassium concentration measured simultaneously (Supplementary Fig. 8) is very similar to that obtained by the Epileptor in A.

Figure 5

Slow permittivity state variable and seizure topology. (A, left) Seizure generated by the Epileptor with five state variables. Seizure onset, time course and offset are controlled by the permittivity state variable evolving slowly in time (red). Note that the SLE occurs with a rapid and large shift of the potential. Right, the seizure trajectory (expressed in terms of −x1 + x2) is approximated in a 3D space defined by the first state variables (X = −x1 and Y = x2) and by the slow permittivity variable (Z = z). Note that the values of the z-variable have been shifted upwards for plotting purposes. (B, left), simultaneous recording in the hippocampus of a SLE in low Mg2+ conditions in DC mode, O2 levels in the preparation (yellow) and NADH levels (red), which indirectly reflect ATP use. Note the large DC shift during the SLE, as predicted by Epileptor. The time course of oxygen and NADH is similar to that of the slow permittivity variable. Right, the 3D representation of a seizure in a delayed space (X and Y), with Z the extracellular potassium concentration measured simultaneously (Supplementary Fig. 8) is very similar to that obtained by the Epileptor in A.

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.

Figure 6

Homoclinic bifurcation at seizure offset in various species. (A) In a mouse hippocampus, interspike intervals display logarithmic scaling typical of a homoclinic bifurcation. The last spike of the seizure is used as our reference time point (red squares mark seizure durations). After accounting for uncertainty in seizure onset time and clonic firing at seizure termination, log scaling fit the data better than other potential models (Supplementary Fig. 9A). The interspike interval from this reference displays a logarithmic scaling, which characterizes a homoclinic bifurcation in the three species. Red line: a log equation fit to the last seven datapoints (t < 10) and extrapolated. (B) Summary of all measures performed in mice hippocampi (n = 16), zebrafish (n = 2) and human (n = 24). Logarithmic scaling is preserved in all species. Note the similarity between the different human subjects, who had a wide array of epilepsy pathologies (Supplementary Table 1). The slope difference relates to differences in seizure duration in the various conditions. (C) Logarithmic scaling in zebrafish. (D) In human, we show independent seizures simultaneously recorded in the right and left hemisphere with different dynamics, both showing log scaling. Inset: corrected equation fit (red) ignores fast spiking between clonic bursts in the last 30 s of seizure. LSS = Left somatosensory cortex; RIT = Right inferior temporal lobe.

Figure 6

Homoclinic bifurcation at seizure offset in various species. (A) In a mouse hippocampus, interspike intervals display logarithmic scaling typical of a homoclinic bifurcation. The last spike of the seizure is used as our reference time point (red squares mark seizure durations). After accounting for uncertainty in seizure onset time and clonic firing at seizure termination, log scaling fit the data better than other potential models (Supplementary Fig. 9A). The interspike interval from this reference displays a logarithmic scaling, which characterizes a homoclinic bifurcation in the three species. Red line: a log equation fit to the last seven datapoints (t < 10) and extrapolated. (B) Summary of all measures performed in mice hippocampi (n = 16), zebrafish (n = 2) and human (n = 24). Logarithmic scaling is preserved in all species. Note the similarity between the different human subjects, who had a wide array of epilepsy pathologies (Supplementary Table 1). The slope difference relates to differences in seizure duration in the various conditions. (C) Logarithmic scaling in zebrafish. (D) In human, we show independent seizures simultaneously recorded in the right and left hemisphere with different dynamics, both showing log scaling. Inset: corrected equation fit (red) ignores fast spiking between clonic bursts in the last 30 s of seizure. LSS = Left somatosensory cortex; RIT = Right inferior temporal lobe.

Seizure topology

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.

Figure 7

Different paths to seizure onset. (A) Stimulations (red stars, traces b, c, d and e) given to Epileptor between two SLEs (trace a shows Epileptor regular behaviour) either failed to trigger a SLE (traces b and c) or generated one before expected (traces d and e). This predicts the presence of a refractory period occurring right after SLE offset and that the system can be pushed toward SLE onset. (B) Experimental verification. Top: The whole hippocampus and the septum connected to each other were placed in two different chambers (inset). After generating a series of SLEs in the hippocampus in low Mg2+ conditions, the extracellular concentration of Mg2+ was raised to 0.4 mM, which maintained the hippocampus below the SLE threshold. The septum was bathed with normal artificial CSF. Bottom: A stimulating electrode was placed in the septum to stimulate axons projecting to the hippocampus. The stimulation generated a small DC shift followed by a SLE. The same train applied after seizure offset failed to evoke a SLE. After waiting >10 min, stimuli of equal magnitude produced a SLE (not shown). (C) Top: noise was progressively increased in Epileptor until seizure onset was reached leading to the prediction that synaptic noise is sufficient to drive the system to the bifurcation. Note that the number of spikes before seizure onset scales with the distance to the bifurcation. (D) Experimental validation. The hippocampus (H) was placed in subthreshold conditions as in B. Top: a hippocampal neuron was recorded in voltage clamp mode at +10 mV to measure GABAergic currents. The extracellular concentration of K+ was raised by 5 mM in the septum (S), leading to increased cell firing there, which correlated with an increase in synaptic activity received by the neuron. This led to the occurrence of a SLE. The septum was then returned to normal artificial CSF conditions. Changing osmolarity in the hippocampus with 50 mM mannitol was sufficient to induce a SLE without increasing synaptic noise (Supplementary Table 2). Bottom: Both procedures were synergistic. Raising [K+] by 2 mM in the septum or adding 10 mM mannitol were not sufficient to trigger a SLE by themselves. When they were both combined, a SLE could be evoked, demonstrating that multiple different trajectories can lead to seizure onset. aCSF = artificial CSF; LFP = local field potential; stim = stimulation; TU = arbitrary time units; VC = voltage clamp.

Figure 7

Different paths to seizure onset. (A) Stimulations (red stars, traces b, c, d and e) given to Epileptor between two SLEs (trace a shows Epileptor regular behaviour) either failed to trigger a SLE (traces b and c) or generated one before expected (traces d and e). This predicts the presence of a refractory period occurring right after SLE offset and that the system can be pushed toward SLE onset. (B) Experimental verification. Top: The whole hippocampus and the septum connected to each other were placed in two different chambers (inset). After generating a series of SLEs in the hippocampus in low Mg2+ conditions, the extracellular concentration of Mg2+ was raised to 0.4 mM, which maintained the hippocampus below the SLE threshold. The septum was bathed with normal artificial CSF. Bottom: A stimulating electrode was placed in the septum to stimulate axons projecting to the hippocampus. The stimulation generated a small DC shift followed by a SLE. The same train applied after seizure offset failed to evoke a SLE. After waiting >10 min, stimuli of equal magnitude produced a SLE (not shown). (C) Top: noise was progressively increased in Epileptor until seizure onset was reached leading to the prediction that synaptic noise is sufficient to drive the system to the bifurcation. Note that the number of spikes before seizure onset scales with the distance to the bifurcation. (D) Experimental validation. The hippocampus (H) was placed in subthreshold conditions as in B. Top: a hippocampal neuron was recorded in voltage clamp mode at +10 mV to measure GABAergic currents. The extracellular concentration of K+ was raised by 5 mM in the septum (S), leading to increased cell firing there, which correlated with an increase in synaptic activity received by the neuron. This led to the occurrence of a SLE. The septum was then returned to normal artificial CSF conditions. Changing osmolarity in the hippocampus with 50 mM mannitol was sufficient to induce a SLE without increasing synaptic noise (Supplementary Table 2). Bottom: Both procedures were synergistic. Raising [K+] by 2 mM in the septum or adding 10 mM mannitol were not sufficient to trigger a SLE by themselves. When they were both combined, a SLE could be evoked, demonstrating that multiple different trajectories can lead to seizure onset. aCSF = artificial CSF; LFP = local field potential; stim = stimulation; TU = arbitrary time units; VC = voltage clamp.

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.

Figure 8

Possible physiological correlates of ensembles 1 and 2 of Epileptor. Cell attached recording A and whole cell current clamp recording B of two stratum oriens interneurons in the CA1 region during a SLE. (A) The GABA neuron fired at SLE onset during the large spike and wave (1), stopped during the fast oscillation (2) and resumed firing when spike and wave reappeared (3). Note that during the late spike and wave event (3), the GABA neuron stops firing during the wave when the fast oscillation occurs. (B) The current clamp recording shows that the cell stops firing as it enters into depolarization block. (C) Whole cell recording in voltage clamp mode of GABAergic currents. Note the presence of large GABAergic inputs during the large spike and wave before SLE onset (1), their loss during the fast oscillation (2), and their re-occurrence during the late part of the SLE (3). Note that GABAergic currents are absent during the fast oscillation of the late spike and wave complex (3). (D) Whole cell voltage clamp recording of synaptic glutamatergic currents received by a GABA neuron during a SLE. Note that the cell receives strong glutamatergic inputs during all phases of the SLE, including spikes (1 and 3) and fast oscillations (2). Note the remarkable synchrony between the glutamatergic currents and the field during the fast oscillation (2).

Figure 8

Possible physiological correlates of ensembles 1 and 2 of Epileptor. Cell attached recording A and whole cell current clamp recording B of two stratum oriens interneurons in the CA1 region during a SLE. (A) The GABA neuron fired at SLE onset during the large spike and wave (1), stopped during the fast oscillation (2) and resumed firing when spike and wave reappeared (3). Note that during the late spike and wave event (3), the GABA neuron stops firing during the wave when the fast oscillation occurs. (B) The current clamp recording shows that the cell stops firing as it enters into depolarization block. (C) Whole cell recording in voltage clamp mode of GABAergic currents. Note the presence of large GABAergic inputs during the large spike and wave before SLE onset (1), their loss during the fast oscillation (2), and their re-occurrence during the late part of the SLE (3). Note that GABAergic currents are absent during the fast oscillation of the late spike and wave complex (3). (D) Whole cell voltage clamp recording of synaptic glutamatergic currents received by a GABA neuron during a SLE. Note that the cell receives strong glutamatergic inputs during all phases of the SLE, including spikes (1 and 3) and fast oscillations (2). Note the remarkable synchrony between the glutamatergic currents and the field during the fast oscillation (2).

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.

Discussion

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.

Acknowledgements

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.

Funding

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

Supplementary material is available at Brain online.

Abbreviations

    Abbreviations
  • DC

    direct current

  • SLE

    seizure-like event

  • SWE

    spike and wave event

References

Andrew
RD
Seizure and acute osmotic change: clinical and neurophysiological aspects
J Neurol Sci
 , 
1991
, vol. 
101
 (pg. 
7
-
18
)
Arya
R
Kabra
M
Gulati
S
Epilepsy in children with down syndrome
Epileptic Disord
 , 
2011
, vol. 
13
 (pg. 
1
-
7
)
Bak
P
Paczuski
M
Complexity, contingency, and criticality
Proc Natl Acad Sci USA
 , 
1995
, vol. 
92
 (pg. 
6689
-
96
)
Bartolomei
F
Chauvel
P
Wendling
F
Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG
Brain
 , 
2008
, vol. 
131
 
Pt 7
(pg. 
1818
-
30
)
Bazhenov
M
Timofeev
I
Frohlich
F
Sejnowski
TJ
Cellular and network mechanisms of electrographic seizures
Drug Discov Today Dis Models
 , 
2008
, vol. 
5
 (pg. 
45
-
57
)
Beall
CM
Two routes to functional adaptation: Tibetan and Andean high-altitude natives
Proc Natl Acad Sci USA
 , 
2007
, vol. 
104
 
Suppl 1
(pg. 
8655
-
60
)
Blanco
JA
Stead
M
Krieger
A
Stacey
W
Maus
D
Marsh
E
, et al.  . 
Data mining neocortical high-frequency oscillations in epilepsy and controls
Brain
 , 
2011
, vol. 
134
 
Pt 10
(pg. 
2948
-
59
)
Bragin
A
Wilson
CL
Staba
RJ
Reddick
M
Fried
I
Engel
J
Jr
Interictal high-frequency oscillations (80-500 Hz) in the human epileptic brain: entorhinal cortex
Ann Neurol
 , 
2002
, vol. 
52
 (pg. 
407
-
15
)
Breakspear
M
Roberts
JA
Terry
JR
Rodrigues
S
Mahant
N
Robinson
PA
A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis
CerebCortex
 , 
2006
, vol. 
16
 (pg. 
1296
-
313
)
Brodie
MJ
Barry
SJ
Bamagous
GA
Norrie
JD
Kwan
P
Patterns of treatment response in newly diagnosed epilepsy
Neurology
 , 
2012
, vol. 
78
 (pg. 
1548
-
54
)
Brunel
N
Wang
XJ
What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance
J Neurophysiol
 , 
2003
, vol. 
90
 (pg. 
415
-
30
)
Chauviere
L
Doublet
T
Ghestem
A
Siyoucef
SS
Wendling
F
Huys
R
, et al.  . 
Changes in interictal spike features precede the onset of temporal lobe epilepsy
Ann Neurol
 , 
2012
, vol. 
71
 (pg. 
805
-
14
)
Cook
MJ
O'Brien
TJ
Berkovic
SF
Murphy
M
Morokoff
A
Fabinyi
G
, et al.  . 
Prediction of seizure likelihood with a long-term, implanted seizure advisory system in patients with drug-resistant epilepsy: a first-in-man study
Lancet Neurol
 , 
2013
, vol. 
12
 (pg. 
563
-
71
)
Cross
MC
Hohenberg
PC
Pattern Formation outside of Equilibrium
Rev Mod Phys
 , 
1993
, vol. 
65
 pg. 
851
 
Cressman
JR
Ullah
G
Ziburkus
J
Schiff
SJ
Barreto
E
The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics
J Comput Neurosci
 , 
2009
, vol. 
26
 (pg. 
159
-
70
)
Deco
G
Jirsa
VK
McIntosh
AR
Emerging concepts for the dynamical organization of resting-state activity in the brain
Nat RevNeurosci
 , 
2011
, vol. 
12
 (pg. 
43
-
56
)
Deco
G
Jirsa
VK
Robinson
PA
Breakspear
M
Friston
K
The dynamic brain: from spiking neurons to neural masses and cortical fields
PLoS Comput Biol
 , 
2008
, vol. 
4
 pg. 
e1000092
 
Destexhe
A
Rudolph
M
Pare
D
The high-conductance state of neocortical neurons in vivo
Nat Rev Neurosci
 , 
2003
, vol. 
4
 (pg. 
739
-
51
)
El-Hassar
L
Milh
M
Wendling
F
Ferrand
N
Esclapez
M
Bernard
C
Cell domain-dependent changes in the glutamatergic and GABAergic drives during epileptogenesis in the rat CA1 region
J Physiol
 , 
2007
, vol. 
578
 
Pt 1
(pg. 
193
-
211
)
Ermentrout
GB
Terman
DH
Mathematical Foundations of Neuroscience
 , 
2010
New York, Dordrecht, Heidelberg, London
Springer
Foss
J
Longtin
A
Mensour
B
Milton
J
Multistability and delayed recurrent loops
Phys Rev Lett
 , 
1996
, vol. 
76
 (pg. 
708
-
11
)
Foss
J
Milton
J
Multistability in recurrent neural loops arising from delay
J Neurophysiol
 , 
2000
, vol. 
84
 (pg. 
975
-
85
)
Freeman
WJ
Mass action in the nervous system
 , 
1975
New York, San Francisco, London
Academic Press
Friedman
D
Honig
LS
Scarmeas
N
Seizures and epilepsy in Alzheimer's disease
CNS Neurosci Ther
 , 
2012
, vol. 
18
 (pg. 
285
-
94
)
Frohlich
F
Bazhenov
M
Iragui-Madoz
V
Sejnowski
TJ
Potassium dynamics in the epileptic cortex: new insights on an old topic
Neuroscientist
 , 
2008
, vol. 
14
 (pg. 
422
-
33
)
Frohlich
F
Sejnowski
TJ
Bazhenov
M
Network bistability mediates spontaneous transitions between normal and pathological brain states
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
10734
-
43
)
Gerstner
W
Kistler
WM
Spiking neuron models. Single neurons, populations, plasticity
 , 
2002
Cambridge
Cambridge University Press
Gnatkovsky
V
Librizzi
L
Trombin
F
de Curtis
M
Fast activity at seizure onset is mediated by inhibitory circuits in the entorhinal cortex in vitro
Ann Neurol
 , 
2008
, vol. 
64
 (pg. 
674
-
86
)
Haken
H
Synergetics, an Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology
 , 
1983
New York
Springer-Verlag
Hale
J
Koçak
H
Dynamics and bifurcations
 , 
1991
New York
Springer
Heinemann
U
Albrecht
D
Beck
H
Ficker
E
von
HD
Stabel
J
Delayed K+ regulation and K+ current maturation as factors of enhanced epileptogenicity during ontogenesis of the hippocampus of rats
Epilepsy Res Suppl
 , 
1992
, vol. 
9
 (pg. 
107
-
14
)
Heinemann
U
Konnerth
A
Pumain
R
Wadman
WJ
Extracellular calcium and potassium concentration changes in chronic epileptic brain tissue
Adv Neurol
 , 
1986
, vol. 
44
 (pg. 
641
-
61
)
Hindmarsh
JL
Rose
RM
A model of neuronal bursting using three coupled first order differential equations
Proc R Soc Lond B Biol Sci
 , 
1984
, vol. 
221
 (pg. 
87
-
102
)
Huberfeld
G
Menendez de la
PL
Pallud
J
Cohen
I
Le Van
QM
Adam
C
, et al.  . 
Glutamatergic pre-ictal discharges emerge at the transition to seizure in human epilepsy
Nat Neurosci
 , 
2011
, vol. 
14
 (pg. 
627
-
34
)
Hunt
RF
Hortopan
GA
Gillespie
A
Baraban
SC
A novel zebrafish model of hyperthermia-induced seizures reveals a role for TRPV4 channels and NMDA-type glutamate receptors
Exp Neurol
 , 
2012
, vol. 
237
 (pg. 
199
-
206
)
Ikeda
A
Taki
W
Kunieda
T
Terada
K
Mikuni
N
Nagamine
T
, et al.  . 
Focal ictal direct current shifts in human epilepsy as studied by subdural and scalp recording
Brain
 , 
1999
, vol. 
122
 
Pt 5
(pg. 
827
-
38
)
Isomura
Y
Fujiwara-Tsukamoto
Y
Takada
M
A network mechanism underlying hippocampal seizure-like synchronous oscillations
Neurosci Res
 , 
2008
, vol. 
61
 (pg. 
227
-
33
)
Izhikevich
EM
Neural excitability, spiking, and bursting
Int J Bifurcation Chaos
 , 
2000
, vol. 
10
 (pg. 
1171
-
266
)
Jacobs
J
Levan
P
Chatillon
CE
Olivier
A
Dubeau
F
Gotman
J
High frequency oscillations in intracranial EEGs mark epileptogenicity rather than lesion type
Brain
 , 
2009
, vol. 
132
 
Pt 4
(pg. 
1022
-
37
)
Jansen
BH
Rit
VG
Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns
Biol Cybern
 , 
1995
, vol. 
73
 (pg. 
357
-
66
)
Jett
DA
Chemical toxins that cause seizures
Neurotoxicology
 , 
2012
, vol. 
33
 (pg. 
1473
-
5
)
Jiruska
P
Csicsvari
J
Powell
AD
Fox
JE
Chang
WC
Vreugdenhil
M
, et al.  . 
High-frequency network activity, global increase in neuronal activity, and synchrony expansion precede epileptic seizures in vitro
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
5690
-
701
)
Jiruska
P
de Curtis
M
Jefferys
JG
Schevon
CA
Schiff
SJ
Schindler
K
Synchronization and desynchronization in epilepsy: controversies and hypotheses
J Physiol
 , 
2013
, vol. 
591
 
Pt 4
(pg. 
787
-
97
)
Karafin
M
St Louis
EK
Zimmerman
MB
Sparks
JD
Granner
MA
Bimodal ultradian seizure periodicity in human mesial temporal lobe epilepsy
Seizure
 , 
2010
, vol. 
19
 (pg. 
347
-
51
)
Khalilov
I
Esclapez
M
Medina
I
Aggoun
D
Lamsa
K
Leinekugel
X
, et al.  . 
A novel in vitro preparation: the intact hippocampal formation
Neuron
 , 
1997
, vol. 
19
 (pg. 
743
-
9
)
Khalilov
I
Holmes
GL
Ben Ari
Y
In vitro formation of a secondary epileptogenic mirror focus by interhippocampal propagation of seizures
Nat Neurosci
 , 
2003
, vol. 
6
 (pg. 
1079
-
85
)
Kramer
MA
Truccolo
W
Eden
UT
Lepage
KQ
Hochberg
LR
Eskandar
EN
, et al.  . 
Human seizures self-terminate across spatial scales via a critical transition
Proc Natl Acad Sci USA
 , 
2012
, vol. 
109
 (pg. 
21116
-
21
)
Krook-Magnuson
E
Armstrong
C
Oijala
M
Soltesz
I
On-demand optogenetic control of spontaneous seizures in temporal lobe epilepsy
Nat Commun
 , 
2013
, vol. 
4
 pg. 
1376
 
Kuznetsov
YA
Elements of applied bifurcation theory
 , 
1998
Springer
New York
Lopes da Silva
F
Blanes
W
Kalitzin
SN
Parra
J
Suffczynski
P
Velis
DN
Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity
Epilepsia
 , 
2003a
, vol. 
44
 
Suppl 12
(pg. 
72
-
83
)
Lopes da Silva
FH
Blanes
W
Kalitzin
SN
Parra
J
Suffczynski
P
Velis
DN
Dynamical diseases of brain systems: different routes to epileptic seizures
IEEE Trans Biomed Eng
 , 
2003b
, vol. 
50
 (pg. 
540
-
8
)
Lopes da Silva
FH
Pijn
JP
Wadman
WJ
Dynamics of local neuronal networks: control parameters and state bifurcations in epileptogenesis
Prog Brain Res
 , 
1994
, vol. 
102
 (pg. 
359
-
70
)
Luttges
MW
McGaugh
JL
Permanence of retrograde amnesia produced by electroconvulsive shock
Science
 , 
1967
, vol. 
156
 (pg. 
408
-
10
)
Marder
E
Taylor
AL
Multiple models to capture the variability in biological neurons and networks
Nat Neurosci
 , 
2011
, vol. 
14
 (pg. 
133
-
8
)
Marten
F
Rodrigues
S
Benjamin
O
Richardson
MP
Terry
JR
Onset of polyspike complexes in a mean-field model of human electroencephalography and its application to absence epilepsy
Philos Trans A Math Phys Eng Sci
 , 
2009
, vol. 
367
 (pg. 
1145
-
61
)
Milton
JG
Epilepsy as a dynamic disease: a tutorial of the past with an eye to the future
Epilepsy Behav
 , 
2010
, vol. 
18
 (pg. 
33
-
44
)
Mody
I
Lambert
JD
Heinemann
U
Low extracellular magnesium induces epileptiform activity and spreading depression in rat hippocampal slices
J Neurophysiol
 , 
1987
, vol. 
57
 (pg. 
869
-
88
)
Moody
WJ
Futamachi
KJ
Prince
DA
Extracellular potassium activity during epileptogenesis
Exp Neurol
 , 
1974
, vol. 
42
 (pg. 
248
-
63
)
Nakken
KO
Solaas
MH
Kjeldsen
MJ
Friis
ML
Pellock
JM
Corey
LA
Which seizure-precipitating factors do patients with epilepsy most frequently report?
Epilepsy Behav
 , 
2005
, vol. 
6
 (pg. 
85
-
9
)
Nunez
PL
The brain wave equation: a model for EEG
Math Biosci
 , 
1974
, vol. 
21
 (pg. 
279
-
97
)
Osorio
I
Manly
B
Sunderam
S
Toward a quantitative multivariate analysis of the efficacy of antiseizure therapies
Epilepsy Behav
 , 
2010
, vol. 
18
 (pg. 
335
-
43
)
Paz
JT
Davidson
TJ
Frechette
ES
Delord
B
Parada
I
Peng
K
, et al.  . 
Closed-loop optogenetic control of thalamus as a tool for interrupting seizures after cortical injury
Nat Neurosci
 , 
2013
, vol. 
16
 (pg. 
64
-
70
)
Penovich
PE
Helmers
S
Catamenial epilepsy
Int Rev Neurobiol
 , 
2008
, vol. 
83
 (pg. 
79
-
90
)
Perucca
P
Dubeau
F
Gotman
J
Intracranial electroencephalographic seizure-onset patterns: effect of underlying pathology
Brain
 , 
2013
, vol. 
137
 
Pt 1
(pg. 
183
-
96
)
Pitkanen
A
Sutula
TP
Is epilepsy a progressive disorder? Prospects for new therapeutic approaches in temporal-lobe epilepsy
Lancet Neurol
 , 
2002
, vol. 
1
 (pg. 
173
-
81
)
Pumain
R
Menini
C
Heinemann
U
Louvel
J
Silva-Barrat
C
Chemical synaptic transmission is not necessary for epileptic seizures to persist in the baboon Papio papio
Exp Neurol
 , 
1985
, vol. 
89
 (pg. 
250
-
8
)
Quilichini
PP
Le Van
QM
Ivanov
A
Turner
DA
Carabalona
A
Gozlan
H
, et al.  . 
Hub GABA neurons mediate gamma-frequency oscillations at ictal-like event onset in the immature hippocampus
Neuron
 , 
2012
, vol. 
74
 (pg. 
57
-
64
)
Raol
YH
Brooks-Kayal
AR
Experimental models of seizures and epilepsies
Prog Mol Biol Transl Sci
 , 
2012
, vol. 
105
 (pg. 
57
-
82
)
Rinzel
J
Teramoto
E
Yamaguti
M
A formal classification of bursting mechanisms in excitable systems
Mathematical topics in population biology, morphogenesis and neurosciences, vol. 71 of lecture notes in biomathematics
 , 
1987
New York
Springer
Robinson
SJ
Childhood epilepsy and autism spectrum disorders: psychiatric problems, phenotypic expression, and anticonvulsants
Neuropsychol Rev
 , 
2012
, vol. 
22
 (pg. 
271
-
9
)
Rodrigues
S
Barton
D
Szalai
R
Benjamin
O
Richardson
MP
Terry
JR
Transitions to spike-wave oscillations and epileptic dynamics in a human cortico-thalamic mean-field model
J Comput Neurosci
 , 
2009
, vol. 
27
 (pg. 
507
-
26
)
Roy
D
Ghosh
A
Jirsa
VK
Phase description of spiking neuron networks with global electric and synaptic coupling
Phys Rev E Stat Nonlin Soft Matter Phys
 , 
2011
, vol. 
83
 
5 Pt 1
pg. 
051909
 
Shu
Y
Hasenstaub
A
Badoual
M
Bal
T
McCormick
DA
Barrages of synaptic activity control the gain and sensitivity of cortical neurons
J Neurosci
 , 
2003
, vol. 
23
 (pg. 
10388
-
401
)
Somjen
GG
Mechanisms of spreading depression and hypoxic spreading depression-like depolarization
Physiol Rev
 , 
2001
, vol. 
81
 (pg. 
1065
-
96
)
Stacey
W
Le Van Quyen
M
Mormann
F
Schulze-Bonhage
A
What is the present-day EEG evidence for a preictal state?
Epilepsy Res
 , 
2011a
, vol. 
97
 (pg. 
243
-
51
)
Stacey
WC
Kellis
S
Patel
PR
Greger
B
Butson
CR
Signal distortion from microelectrodes in clinical EEG acquisition systems
J Neural Eng
 , 
2012
, vol. 
9
 pg. 
056007
 
Stacey
WC
Krieger
A
Litt
B
Network recruitment to coherent oscillations in a hippocampal computer model
J Neurophysiol
 , 
2011b
, vol. 
105
 (pg. 
1464
-
81
)
Stacey
WC
Lazarewicz
MT
Litt
B
Synaptic noise and physiological coupling generate high-frequency oscillations in a hippocampal computational model
J Neurophysiol
 , 
2009
, vol. 
102
 (pg. 
2342
-
57
)
Stefanescu
RA
Jirsa
VK
A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons
PLoS Comput Biol
 , 
2008
, vol. 
4
 pg. 
e1000219
 
Suffczynski
P
Lopes da Silva
FH
Parra
J
Velis
DN
Bouwman
BM
van Rijn
CM
, et al.  . 
Dynamics of epileptic phenomena determined from statistics of ictal transitions
IEEE Trans Biomed Eng
 , 
2006
, vol. 
53
 (pg. 
524
-
32
)
Suh
M
Ma
H
Zhao
M
Sharif
S
Schwartz
TH
Neurovascular coupling and oximetry during epileptic events
Mol Neurobiol
 , 
2006
, vol. 
33
 (pg. 
181
-
97
)
Takens
F
Rand
DA
Young
LS
Detecting strange attractors in turbulence
Dynamical systems and turbulence, lecture notes in mathematics
 , 
1981
, vol. 
898
 
Springer-Verlag
(pg. 
366
-
81
)
Tallgren
P
Vanhatalo
S
Kaila
K
Voipio
J
Evaluation of commercially available electrodes and gels for recording of slow EEG potentials
Clin Neurophysiol
 , 
2005
, vol. 
116
 (pg. 
799
-
806
)
Teramae
JN
Fukai
T
Local cortical circuit model inferred from power-law distributed neuronal avalanches
J Comput Neurosci
 , 
2007
, vol. 
22
 (pg. 
301
-
12
)
Trevelyan
AJ
Sussillo
D
Yuste
R
Feedforward inhibition contributes to the control of epileptiform propagation speed
J Neurosci
 , 
2007
, vol. 
27
 (pg. 
3383
-
7
)
Truccolo
W
Donoghue
JA
Hochberg
LR
Eskandar
EN
Madsen
JR
Anderson
WS
, et al.  . 
Single neuron dynamics in human focal epilepsy
Nat Neurosci
 , 
2011
, vol. 
14
 (pg. 
635
-
43
)
Valentin
A
Alarcon
G
Honavar
M
Garcia Seoane
JJ
Selway
RP
Polkey
CE
, et al.  . 
Single pulse electrical stimulation for identification of structural abnormalities and prediction of seizure outcome after epilepsy surgery: a prospective study
Lancet Neurol
 , 
2005
, vol. 
4
 (pg. 
718
-
26
)
Vanhatalo
S
Holmes
MD
Tallgren
P
Voipio
J
Kaila
K
Miller
JW
Very slow EEG responses lateralize temporal lobe seizures: an evaluation of non-invasive DC-EEG
Neurology
 , 
2003
, vol. 
60
 (pg. 
1098
-
104
)
Walker
MC
The potential of brain stimulation in status epilepticus
Epilepsia
 , 
2011
, vol. 
52
 
Suppl 8
(pg. 
61
-
3
)
Wang
Y
Goodfellow
M
Taylor
PN
Baier
G
Phase space approach for modeling of epileptic dynamics
Phys Rev E Stat Nonlin Soft Matter Phys
 , 
2012
, vol. 
85
 
6 Pt 1
pg. 
061918
 
Wendling
F
Bartolomei
F
Bellanger
JJ
Chauvel
P
Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition
Eur J Neurosci
 , 
2002
, vol. 
15
 (pg. 
1499
-
508
)
White
H
A Heteroskedasticity-consistent covariance matrix estimator and a direct test for Heteroskedasticity
Econometrica
 , 
1980
, vol. 
48
 (pg. 
817
-
38
)
Wilson
HR
Cowan
JD
Excitatory and inhibitory interactions in localized populations of model neurons
Biophys J
 , 
1972
, vol. 
12
 (pg. 
1
-
24
)
Worrell
GA
Gardner
AB
Stead
SM
Hu
S
Goerss
S
Cascino
GJ
, et al.  . 
High-frequency oscillations in human temporal lobe: simultaneous microwire and clinical macroelectrode recordings
Brain
 , 
2008
, vol. 
131
 
Pt 4
(pg. 
928
-
37
)
Zakaria
T
Noe
K
So
E
Cascino
GD
Wetjen
N
Van Gompel
JJ
, et al.  . 
Scalp and intracranial EEG in medically intractable extratemporal epilepsy with normal MRI
ISRN Neurol
 , 
2012
, vol. 
2012
 pg. 
942849
 
Zhang
H
Tan
J
Reynolds
E
Kuebler
D
Faulhaber
S
Tanouye
M
The Drosophila slamdance gene: a mutation in an aminopeptidase can cause seizure, paralysis and neuronal failure
Genetics
 , 
2002
, vol. 
162
 (pg. 
1283
-
99
)
Zhao
M
Nguyen
J
Ma
H
Nishimura
N
Schaffer
CB
Schwartz
TH
Preictal and ictal neurovascular and metabolic coupling surrounding a seizure focus
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
13292
-
300
)

Author notes

*These authors contributed equally to this work.
See doi:10.1093/brain/awu147 for the scientific commentary on this article.