On the nature of seizure dynamics

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% O₂/5% CO₂) ice cold (4°C) artificial CSF containing: 126 mM NaCl; 3.5 mM KCl; 2 mM CaCl₂; 1.2 mM MgCl₂; 25 mM NaHCO₃; 1.2 mM NaHPO₄; 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 CO₂/O₂ gas mixture. The pH was controlled during each experiment. After a 30 min baseline recording in Mg²⁺-containing artificial CSF solution, the media was changed to one without added Mg²⁺ in the chamber containing the hippocampus. Under this condition, the extracellular concentration of Mg²⁺ does not necessarily decrease to zero because there is minimal Mg²⁺ 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 Mg²⁺ 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 MgCl₂; 0.1 mM CaCl₂; 1 mM EGTA; 2 mM MgATP; 0.5 mM Na₄GTP; 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 GABA_A 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

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 R², 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 Ca²⁺ recordings

Seizure-like events recorded when lowering extracellular Ca²⁺ were kindly provided by Dr. J. R. Jefferys and Dr. Jiruska Premysl. The experimental procedure is described in Jiruska et al. (2010).

Basic building blocks of seizure dynamics

As we searched for invariant features in seizure dynamics, any experimental model system may be used. Here, we analysed SLEs recorded in the intact immature hippocampus of mice in vitro, as this preparation is easily amenable to experimentation. When placed in continuous epileptogenic conditions, the preparation produces spontaneous recurrent SLEs (Fig. 1, unfiltered record). In Fig. 2A, we show one such typical spontaneous SLE (0.01 Hz high-pass filter). As observed in many types of seizure (Perucca et al., 2013), SLEs are characterized by a beginning (onset), various sequences of fast discharges and spike and wave events (SWEs), and an end (offset). Here, we consider fast discharges and SWEs as basic building blocks of SLEs, i.e. their temporal arrangement, embedding and amplitude can vary indefinitely without affecting the generality of the results. Similar dynamic patterns were found in a zebrafish model of hyperthermia-induced seizure (Fig. 2B) and in human patients with epilepsy (Fig. 2C and Supplementary Table 1).

Time-evolving phenomena can be formalized with differential equations and state variables, which describe the fundamental dynamics of the system. State variables are the smallest possible subset of system variables that can represent the entire state of the system at a given point in time (Supplementary material). They are not unique and often difficult to relate to directly measurable biological/biophysical variables. These variables can be plotted in the space spanned by the state variables (the ‘state space’) to demonstrate the characteristic ‘flow’, which determines how the trajectories evolve. In dynamic system theory, time scale separation allows the grouping of state variables into subsets acting on the same time scale (Haken, 1983) that can be considered as the building blocks of a dynamic system. For SLEs, the two building blocks, fast discharges and SWEs operate on different time scales (fast and slow, respectively). We will refer to each building block as an ‘ensemble’. Fast discharges necessitate a first ensemble with at least two state variables (x₁, y₁) 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 (x₂, y₂). 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 (x₁, y₁) and (x₂, y₂) 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 (x₁, y₁) and (x₂, y₂). 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, (x₁, y₁) and (x₂, y₂), 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 (x₁, y₁) and (x₂, y₂) that we capture via a bidirectional coupling f₁(x₁, x₂) and f₂(x₁, x₂). 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 (x₁, y₁) 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.

Detailed descriptions of the various bifurcation types and their normal forms (canonical models) can be found in Kuznetsov (1998). Note that normal forms classify local bifurcations around equilibria in state space, whereas global bifurcations involving ‘larger’ invariant sets such as periodic orbits are significantly more complicated. To aid in the classification of seizure types in the framework of a taxonomy of SLEs, the scaling properties of the seizure onset/offset bifurcations are listed as bifurcation of equilibrium in Table 2 and as bifurcation of oscillations in Table 3. As the bifurcation point of seizure onset/offset is approached through the slow change of the permittivity z, the discharges change their behaviour. If denotes the distance to the bifurcation point , then the changes in the discharges will scale with in ways that are characteristic for each bifurcation type. The frequency and amplitude of discharges may be either constant and independent of (denoted by ‘Fixed’ in Tables 2 and 3), completely arbitrary (‘Arbitrary’) or scale from zero following a square root or logarithmic behaviour [‘Zero ()’, ‘Zero ()’]. Each bifurcation may also be distinguished through the number of states that can coexist for the same permittivity value (‘Bistable’ and ‘Monostable’). In the following, we use the taxonomy of SLEs in Table 1 to characterize the time course and recurrence of the SLEs in our experimental conditions.

Table 2.

Properties of oscillatory onset bifurcations

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

Open in a new tab

The four types of bifurcation of equilibrium identify the rows. The column Behaviour indicates if bistability between the oscillatory and equilibrium state is possible, or if the system remains fixed in one state (monostable). The columns Frequency and Amplitude indicate how the onset oscillation scales after the bifurcation of the equilibrium. Here the parameter represents the distance to the bifurcation point. SNIC = saddle-node on invariant cycle.

Table 3.

Properties of oscillatory offset bifurcations

Bifurcation of oscillations	Behaviour	Frequency	Amplitude
SNIC	Monostable	Zero ()	Fixed
Supercritical Hopf	Monostable	Fixed	Zero ()
Fold Limit cycle	Bistable	Fixed	Arbitrary
Saddle homoclinic	Bistable	Zero ()	Fixed

Open in a new tab

The four types of bifurcation of oscillations identify the rows. The column Behaviour indicates if bistability between the oscillatory and equilibrium state is possible. The columns Frequency and Amplitude indicate how the offset oscillation scales before the bifurcation to the equilibrium. Here the parameter represents the distance to the bifurcation point. SNIC = saddle-node on invariant cycle.

Saddle-node bifurcation at seizure onset

In the experiment, the onset of SLEs is characterized by the abrupt appearance of fast discharges (Fig. 2A). As these do not scale up from zero amplitude or frequency, this transition is limited to either a subcritical Hopf or a saddle-node bifurcation (Table 2). A saddle-node bifurcation requires a baseline shift of the measured time-dependent variable, whereas a subcritical Hopf bifurcation results in a transition centered on the baseline and does not have a baseline shift. DC shifts from the baseline field potential occurred systematically in our experimental data (Figs 1, 5B and 7B), ruling out the subcritical Hopf bifurcation.

Homoclinic bifurcation at seizure offset

As SLEs progress towards offset, there is a DC shift back to baseline (Figs 1, 5B and 7B). Only the fold limit cycle and the homoclinic bifurcations show a DC shift at the bifurcation point (indicated by ‘Bistability’ in Table 3), although secondary bifurcations can be introduced to mimick this behaviour (Izhikevich, 2000). The fold limit cycle bifurcation maintains a constant frequency towards the seizure offset, which was not found in the majority of our experiments. The predominant candidate at seizure offset is thus the homoclinic bifurcation, which, amongst others, predicts that the interspike intervals scale logarithmically as seizure offset approaches. The prediction of a logarithmic scaling will be validated experimentally in the following section.

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, (x₁, y₁) and (x₂, y₂), 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 (x_1, y₁) and the mathematical form of an excitable system with a saddle node on invariant circle bifurcation for ensemble 2 with (x_2, y₂) 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, x₁ + x₂, 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 x₁ + x₂ as a function of time bore striking resemblance with the field potential. The complete system of the Epileptor equations reads then as follows:

where

and x₀ = −1.6; y₀ = 1; τ₀ = 2857; τ₁ = 1; τ₂ = 10; I_rest1 = 3.1; I_rest2 = 0.45; γ = 0.01. Here the state variables x₁ and y₁ comprise the first subsystem responsible for fast oscillations and x₂ and y₂ the second subsystem involved in spike wave events. The slow permittivity variable is z. The characteristic time scales are τ₀ of the permittivity variable, τ₁ of ensemble 1 and τ₂ of ensemble 2, where the time scale hierarchy is τ₀ ≫ τ₂ ≫ τ₁. The time constant of τ₁ 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(x₁) 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 x₁ = 0; y₁ = −5; z = 3; x₂ = 0; y₂ = 0. Noise is introduced into each equation as linear additive Gaussian white noise with zero mean and a variance of 0.025 for the first subsystem and 0.25 for the second subsystem. For the solution of the stochastic equations we employ the Euler-Maruyama method.

Dynamics of the Epileptor

In Epileptor, ‘normal’ brain function and seizures occupy two different regions of state space, being isolated from one another by a divergent flow acting as a barrier (Fig. 3). This barrier between the two states establishes bistability and is known as the separatrix, which in this case essentially describes the seizure threshold (Frohlich et al., 2010; Kramer et al., 2012). We show the bifurcations of seizure onset (Fig. 3, rows Ia and Ib) and seizure offset (Fig. 3, rows IIa and IIb) as cartoons of six potential landscapes (rows Ia and IIa) and six vector fields (rows Ib and IIb). The potential landscapes metaphorically illustrate the transitions between seizure and non-seizure states for changing values of permittivity z and the vector fields show the corresponding flows in the two-dimensional state space of the Epileptor’s ensemble 1. In rows Ia and IIa, two minima are separated by a potential well. The left minimum is the normal state, the right minimum is the seizure state. As the permittivity variable changes, the potential landscape changes and one minimum becomes a maximum resulting in a transition from one state to the other. The flows in state space show the same behaviour (Ib and IIb). Here we plot trajectories for various initial conditions, where the arrows indicate the direction of the flow and circles indicate the equilibrium points (fixed points). A stable fixed point attracts trajectories in its neighborhood (full circle), whereas unstable fixed points deflect the trajectories (empty circle). A so-called saddle point is a fixed point with a stable and unstable direction (empty circle). In the first figure of Ib (from left to right) a saddle point separates the state space in two regions with a stable fixed point to its left and a stable limit cycle to its right. As the permittivity is changed towards seizure onset the separatrix moves closer to the stable fixed point (middle figure, row Ib) until the separatrix collides with the stable fixed point and the two disappear (right figure, row Ib) via a saddle-node bifurcation. Row IIb shows the corresponding scenario for seizure offset: as the permittivity changes, the separatrix moves towards the limit cycle (middle figure) until separatrix and limit cycle collide and annihilate each other (right figure). This event is called the homoclinic bifurcation and leaves the system with only the stable fixed point. A bifurcation diagram captures the hitherto described dynamics more quantitatively in Fig. 4, which has been generated using analytical calculation of fixed points and numerical continuation. Fig. 5A traces out time series from simulations of the Epileptor model. Note the presence of SWEs before seizure onset, even far from it. The presence of more interictal spikes (and even fast oscillations) reflects that the Epileptor is getting close to the separatrix. The experimental/clinical analogy would be that network excitability increases and generates spikes and/or oscillations. These activities have their own dynamics; they can take various shapes (reflecting the diversity of interictal/preictal states) as the system moves close or away from the separatrix. This concept of getting close to the separatrix and away from it without reaching the bifurcation point is analogous to the ‘preictal state’ that has been theorized for many years (Stacey et al., 2011a) and may have been identified in a recent clinical trial (Cook et al., 2013).

Because the Epileptor was constructed on the basis of an experimental model of SLEs (an immature hippocampus placed in continuous epileptogenic conditions in vitro), we first tested the generality of its predictions in other brain regions and species, and then explored the concept of the separatrix to understand how seizure onset can be reached in more realistic conditions.

Validating model predictions: bifurcations and invariance

The Epileptor model makes two important predictions about seizure onset and offset in the seizure class fold/homoclinic (Table 1) that we can validate directly: (i) seizure onset can only occur in the presence of a DC shift of the field potential; and (ii) the interspike intervals show a logarithmic scaling approaching seizure offset.

Presence of a direct current shift at seizure onset

We used the presence of a DC shift at seizure onset and its reversal after seizure offset in our experimental model to identify the onset and offset bifurcations (Fig. 1). This DC shift is smaller (−0.26 ± 0.08 mV, n = 15 SLEs from five hippocampi, P < 0.01) and shorter lasting (1–3 min) than that found in spreading depression (Somjen, 2001). Like most clinical and experimental EEG, seizures are recorded in alternating current mode, thus filtering any slow variations of the field potential (Fig. 2). The DC shift at onset of a partial spontaneous seizure, is actually well known in several species including baboons (Pumain et al., 1985) and humans (Ikeda et al., 1999; Vanhatalo et al., 2003), though it is under-recognized clinically because clinical electrodes are poor at recording DC (Tallgren et al., 2005; Stacey et al., 2012) and seizures are rarely recorded with DC coupling. To the best of our knowledge, we are not aware of published spontaneous focal seizures recorded in DC mode without a DC shift.

Logarithmic scaling of interspike intervals approaching seizure offset

The presence of a homoclinic bifurcation at seizure offset imposes a logarithmic scaling of interspike intervals. Introducing noise in Epileptor enabled us to generate numerous SLEs and verify that SLE offsets in the Epileptor show the logarithmic scaling of homoclinic bifurcations (Supplementary Fig. 1). In all SLEs recorded in our experimental conditions in vitro (n = 16 hippocampi), the interspike intervals exhibited a logarithmic scaling at seizure offset (Fig. 6A and B), verifying the prediction that seizure offset corresponds to a homoclinic bifurcation. We then tested the interspike intervals in each of the first 24 patients with identifiable seizures in the International EEG database (www.ieeg.org). This arbitrary sample includes several different seizure types and brain regions (Supplementary Table 1). We verified logarithmic scaling in 20 (83%) of the patients (Supplementary material), the majority of which had remarkably similar dynamics despite their clinical disparities (Fig. 6B, Supplementary Fig. 2 and Supplementary Table 1). Similar findings, which appear visually as spikes slowing down near the end of seizures, were recently identified in live human recordings on multiple spatial scales from EEG down to multiunits (Kramer et al., 2012). The same property was found in other species, such as zebrafish (Fig. 6B and C, Supplementary Fig. 9), and is apparent also in flies [see Fig. 2C in Zhang et al. (2002)]. These invariant scaling laws argue strongly in favour that a homoclinic bifurcation at offset is an invariant property of focal seizures.

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 x₁, x₂ and z as a function of time. The resulting topology corresponds to spirals on a cone (Fig. 5A). Although −x₁+x₂ bears analogy with the field potential, the precise biophysical equivalent of the z variable is unknown and will be likely complex. One of its distinguishing features is its slow time-dependent evolution. Several biophysical parameters are known to evolve slowly in time during seizures, including extracellular ions (Heinemann et al., 1986) and oxygen (Suh et al., 2006). We measured extracellular [K⁺] during SLEs and plotted it as the z variable; the topology of experimental SLEs were remarkably similar to the theoretical one (Fig. 5B). As the measurement of a slow variable was not available in the humans and zebrafish, we constructed the trajectories in a space of delayed state variables, which approximates the original state space (Takens, 1981) and allows unfolding the resulting trajectories for illustration of the flow topology. Even though the trajectory lines are distorted and difficult to read, they still unveil the basic phases of the seizures, showing clear similarities between the Epileptor and SLEs (Supplementary Fig. 3).

Validating model predictions: how seizures start

Seizures can occur with regularity in some patients, as for catamenial epilepsy (Penovich and Helmers, 2008) or during the sleep-wake cycle (Karafin et al., 2010), in keeping with the possibility that ultra-slow systemic effects, analogous to the slow permittivity variable, can push the system periodically to seizure onset. However, seizures arise, most of the time, without obvious causes or predictive factors. In such instances, we propose that brain trajectories are close to the separatrix, and the closer the ‘normal’ brain state is to the separatrix, the easier it will be to have a seizure. Linked to the distance to the separatrix (i.e. proximity to seizure threshold) is the build-up of preictal spiking in the Epileptor (Fig. 7C), which is analogous to physiological fluctuations around baseline such as increased preictal spiking (Huberfeld et al., 2011), other preictal waveforms (Stacey et al., 2011a) or predictive features in the EEG (Cook et al., 2013). However, these fluctuations are quite variable, do not inevitably progress to a seizure (Cook et al., 2013), and are still dominated by the normal baseline activity. Thus, these excursions have not yet crossed the separatrix and provide evidence that (i) the system is close to threshold; and (ii) there are many potential trajectories from the ‘normal’ state to the ‘seizure’ state. The stereotypical appearance of a seizure despite the numerous trajectories leading into it provides further evidence that it is a distinct state of the system. In the Epileptor, it is possible to know the distance from the separatrix (bifurcation diagram in Fig. 4), and force the system pass its threshold. We shall explore two possibilities for provoking seizures and validate them experimentally: large external stimulation, as well as timely internal noise.

Triggering seizures with electrical stimulation

The most straightforward method to produce an excursion of brain trajectories toward the separatrix is to manipulate the system directly. In Epileptor, this can be simply done by adding current, which moves the system towards the separatrix and triggers a SLE before its expected time of occurrence (Fig. 7A, traces d and e). Electrical stimulations consistently trigger seizures in humans: following electroconvulsive shocks in any ‘normal’ brain (Walker, 2011) and following cortical stimulation during presurgical evaluation of epileptic patients (Valentin et al., 2005). Another prediction, which is not as well described clinically, is that there is a refractory period after an SLE offset in which the same ‘kick’ cannot produce a seizure (Fig. 7A, traces b and c). The system will be refractory outside of the interval of bistability (Fig. 4). To test both predictions experimentally, we isolated the whole mouse hippocampus and the septum connected to it (Fig. 7B). We placed the hippocampus in one chamber and the septum in another, but maintained the physical connection between them. Each chamber could be perfused with different solutions without exchange between both compartments (Supplementary Fig. 4). Bathing the hippocampus in low Mg²⁺ artificial CSF resulted in recurrent SLEs with a typical DC shift (Fig. 7B). The concentration of Mg²⁺ 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 Mg²⁺.

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⁺], pO₂, 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⁺], pO₂, 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 (x₁, y₁) 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 x₁ and slower variable x₂, 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 (x₂, y₂), which is why the biophysical correlate of the Epileptor is more complex than the simple interpretation of ensembles 1 and 2 as excitatory and inhibitory cells.

The previous results show that it is possible to construct hypothesis-driven search of the biophysical variables contributing to the state variable, which may be specific to experimental conditions (as demonstrated hereafter), and in the clinic, to any specific patient.

Seizures: conserved dynamics from many trajectories

It is important to stress here again that the Epileptor model does not impose constraints upon the nature of the biological variables; only upon their dynamics and respective relationships. This means that multiple parameter configurations can theoretically give rise to the same system behaviours. For instance, work in the guinea pig brain suggests that the fast discharges therein are produced by inhibitory interneurons, rather than by pyramidal cells (Gnatkovsky et al., 2008). This concept, that there are many functional biological pathways that produce the same outcome, was developed and verified experimentally in simple neuronal networks (Marder and Taylor, 2011) and is present in human biology as well (Beall, 2007). This issue is particularly relevant both clinically and experimentally, as numerous anatomical, molecular, electrophysiological and functional modifications have been described in different types of epilepsies, experimental models and species; without clear consensus emerging about major contributing parameters. Hence the possible correlates of the five state variables in the low Mg²⁺ 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 Ca²⁺, 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 Ca²⁺, 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 Ca²⁺ 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.