Scaling effects and spatio-temporal multilevel dynamics in epileptic seizures.

Epileptic seizures are one of the most well-known dysfunctions of the nervous system. During a seizure, a highly synchronized behavior of neural activity is observed that can cause symptoms ranging from mild sensual malfunctions to the complete loss of body control. In this paper, we aim to contribute towards a better understanding of the dynamical systems phenomena that cause seizures. Based on data analysis and modelling, seizure dynamics can be identified to possess multiple spatial scales and on each spatial scale also multiple time scales. At each scale, we reach several novel insights. On the smallest spatial scale we consider single model neurons and investigate early-warning signs of spiking. This introduces the theory of critical transitions to excitable systems. For clusters of neurons (or neuronal regions) we use patient data and find oscillatory behavior and new scaling laws near the seizure onset. These scalings lead to substantiate the conjecture obtained from mean-field models that a Hopf bifurcation could be involved near seizure onset. On the largest spatial scale we introduce a measure based on phase-locking intervals and wavelets into seizure modelling. It is used to resolve synchronization between different regions in the brain and identifies time-shifted scaling laws at different wavelet scales. We also compare our wavelet-based multiscale approach with maximum linear cross-correlation and mean-phase coherence measures.

Introduction

Trying to predict epileptic seizures using time series analysis has been an important research topic for decades. In particular, the now wide-spread use of EEG (electroencephalography) techniques to acquire data has been a major driving force of the subject. The review article [1] and the recent book [2] provide perspectives what has been achieved in seizure prediction. The main goal was to identify and characterize a pre-ictal phase occurring before the onset and to design measures that approximately predict the critical starting time of the seizure [3]. Since research has focused in this direction there are still gaps [4] in our understanding of seizures from a dynamical systems perspective [5]–[7]. In this paper, we are going to address this issue and focus on dynamical mechanisms as e.g. in [8] instead of aiming at a predictive technique for seizures.

The main themes of our results are the deep links to mathematical multiscale techniques [9], [10] and the observation of scaling laws at different spatio-temporal levels. From models based on biophysical principles of brain dynamics it is expected that multiple spatial [11] and multiple time scales [12] play an important role for epileptic seizures [13]. Based on a combination of analyzing epileptic seizure patient data and neuron modelling we split the problem into three spatial scales and show that at each individual spatial level the problem exhibits multiple time scale behaviour. We point out that our approach to verify the existence of multiscale phenomena is primarily data-driven and complements modelling approachs (see e.g. [14]).

On the smallest spatial scale, we employ model-based analysis of single neurons [15], [16] using a multiple time scale stochastic FitzHugh-Nagumo model [17]–[19] with a focus on early-warning signs [20] of spiking and scaling laws. In particular, we investigate three different cases of spiking and provide the first results of scaling laws in critical transition theory [21] for neurons in an excitable state. Scaling law results for systems without equlibria near bifurcations have recently been applied successfully in climate modeling [22], [23] and in ecological systems [24], [25]. Apparently these techniques have not been applied to neuroscience problems yet although the phenomenon of slowing down has been found in neuronal systems [26]. We analyze three different regimes for the relationship between noise and time scale separation and show that the variance can be a precursor of spiking in some parameter regimes while it fails in the low noise case. In this context, we point out that the distributions of interspike intervals [27] has been studied extensively in single neuron models but that our work only studies the time series locally near a bifurcation and does not require multiple events.

The second spatial scale which we consider are clusters/regions of neurons [28], [29]. Here we use electrocorticogram (ECoG) data; see Materials section. We examine the onset of the epileptic seizure using the variance as a simple univariate measure. We observe that during a certain period before the seizure the variance shows oscillations. Furthermore, very close to the transition to a seizure the inverse of the variance displays a linear scaling law. Based on critical transition theory, these observations are generically characteristics for Hopf bifurcation [30]. It is very important to note that many seizure models [31]–[35] suggest a Hopf bifurcation as a main mechanism as the transition point. Therefore, our results not only provide a first application of local scaling laws near bifurcations to data but also validate the proposed bifurcation mechanism arising from biophysical principals. Similar to the individual neurons scale, we point out that distributions of interseizure intervals have been studied [36] but that we do not require multiple events.

On the largest spatial scale we analyze the synchronization and correlation between different brain regions [37]. Several bivariate measures have been proposed [1] to study epileptic seizures but the underlying complex network structure makes the problem difficult [38]. Our approach utilizes a recent technique calculating phase-locking intervals (PLIs) [39] based on wavelet transforms [10]. Wavelet-based methods have been applied previously in the context of epileptic seizures [40] but our approach is the first to investigate PLIs and associated phase-locking. We show that our wavelet-based method [10], [39] measures increasing phase-locking and resolves a multiple time scale structure near the seizure onset. Furthermore, we observe a linear scaling law of average phase-locking and that phase-locking at different scales often starts at different times. These results apply near the seizure onset and could potentially relate to recently observed rapid discharges [41], [42]. We also compare our results to other bivariate measures such as maximum linear cross-correlation [43], [44] and mean phase coherence [45].

In summary, our study introduces two recently developed methods (critical transitions, PLIs) into the analysis of epileptic seizures. Using critical transitions theory we give the first analysis of early-warning signs for excitable neurons, identify a potential Hopf bifurcation as the seizure onset mechanism from data and find a new scaling law of single-event time series data at the cluster level. For the wavelet-based phase-locking technique, we provide a comparative study to other bivariate measures and discover a scaling law occurring at time-shifted onset times. On each of the three spatial levels we also identified a multiple time scales structure, based on a data-driven time series approach.

Results

Single Neurons

We start on the level of single neurons. Clearly it is very problematic to get data in this case before epileptic seizures so that we resort to model neurons. The main question will be whether we can predict a spike in the voltage time trace of the model neuron before it occurs. The FitzHugh-Nagumo (FHN) model [17], [18], [46] is a simplification of the Hodgkin-Huxley equations [47] which model the action potential in a neuron. We point out that the methods we are going to present here are going to apply to a much wider class of excitable neuronal models than the FHN equation such as the original Hodgkin-Huxley model [48] or the Morris-Lecar system [49] since these models have similar bifurcation structure and multiple time scale properties [15].

There are several forms of the FHN-equation [50]. One possible version suggested by FitzHugh is the Van der Pol-type [51] modelwhere represents voltage, is the recovery variable and , , are parameters. We think of as an external signal or applied current [52] and assume that the time scale separation satisfies so that is the fast variable and the slow variable. The dynamics of (1) can be understood using a fast-slow decomposition [53]–[55]. Setting in (1) yields a differential equation on the slow time scale defined on the algebraic constraintWe call the critical manifold; see Figure 1. Differentiating implicitly with respect to we find so that the differential equation on can be written aswhich we refer to as slow flow. Observe that the slow flow is not well-defined at the two points . Applying a time re-scaling to the fast time to (1) gives

Figure 1

Simulation of (3) with and using an Euler-Maruyama numerical SDE solver [?]; red curves are deterministic trajectories with and blue curves are sample paths with .

Systems have always been started at . The critical manifold is shown in grey and the -nullcline as a dashed black curve. (a) , the equilibrium for the full system lies on . (b) , the equilibrium lies on near the fold point . The deterministic trajectory has only one spike while noise-induced escapes produce repeated spiking for the stochastic system.

Setting in 2 gives the fast flow where implies that is viewed as a parameter in this context. Observe that consists of equilibrium points for the fast flow and that the points are fold (or saddle-node) bifurcation points [56] in this context. The critical manifold naturally splits into three partswhere are attracting equilibria and are repelling equilibria for the fast flow. We view as the refractory state and as the excited state for the neuron. For trajectories are concatenations of the fast and slow flows. We will consider two different situations for the parameters . In the first situation we chose the parameters so that (1) has a single equilibrium point on where is the -nullcline of the FHN-equation; see Figure 1(a1)–(a2). For suppose that ; then the slow flow moves the system to , a jump via the fast subsystem to occurs, the slow flow on brings the system to and another jump returns it to . This is the classical relaxation oscillation [55], [57]. However, in neuroscience one often also considers the excitable regime [15] where the global equilibrium for the system is stable and lies on close to ; see Figure 1(b1)–(b2). In this case, a trajectory of (1) can generate, depending on , at most one excursion/spike to the excitable state before returning to . Repeated spiking in the excitable regime can be obtained using the more general stochastic FHN-equationwhere is delta-correlated white noise and is a parameter representing the noise level. We can now ask whether individual neuron spiking activity already has precursors. This viewpoint should provide new insights how neurons are able to control synchronization and how control failure occurs. Recent results on predicting critical transitions [20] suggest that statistical precursors can be used to predict events similar to spiking in neurons from a time series without knowing their exact location. The detailed mathematical theory can be found in [21], [30].

Here we present the first application of this theory in the context of single neurons. We want to predict a spiking transition from a neighborhood of to and consider the variance as an early-warning sign

Observe that we can view also as a function of , and write , since the mapping between and is bijective when restricting to . In the relaxation oscillation regime (see Figure 1(a1)–(a2)) and if are sufficiently small it can be shown [30], [58] thatfor some constant and where is small. Therefore an increase in fast voltage-variable variance can potentially be used to predict and to control spiking if no equilibrium exists near . Here we extend the results of [30] by investigating the excitable regime. Figure 2 shows an average of the variance computed over 100 sample paths using a sliding window technique [21]. Figure 2(a) shows the relaxation oscillation regime where we can confirm the theoretical prediction (4).

Figure 2

Average of the variance (black curves) over 100 sample paths starting for at up to a final time .

The green curves are fits of using (4) with fitting parameters and . Fixed parameter values are . (a) Relaxation oscillation regime with . (b) Excitable regime with ; sample paths can exhibit oscillations around the stable focus equilibrium which are visible in the variance. (c) Excitable regime with where larger noise regularizes the variance similar to (a). (d) Excitable regime where smaller noise does not allow fast escapes from and yields decreasing variance.

The excitable regime is much more interesting since the equilibrium point can lead to a variety of distinct regimes depending on the noise level. In Figure 2(b) the noise is at an intermediate level so that deterministic oscillations around the equilibrium are visible in the variance before an escape; hence the prediction (4) is not a good prediction of a spike but one should rely on the oscillatory mechanism before escapes. In Figure 2(c) the noise is larger which provides a regularizing effect for the variance via noise-induced escapes. This relates to the well-known mechanism of coherence resonance [19]. In Figure 2(d) the noise is very small so that sample paths need exponentially long times to escape and are metastable near . This causes a decrease in variance and will make predictions very difficult. The different scaling regimes for noise level and time scale separation are discussed in more detail in [30], [59]–[61].

Based on our results we can conclude predictability of a spiking event and hence also its external control by input currents depend crucially on noise level and statistical properties of the state of a neuron. In particular, in the excitable state already a small change in the noise level or system parameters can result in a substantial loss of control due to unpredictable spiking. This could cause undesirable synchronization and continuous spiking. Let us point out that this is just one possible explanation for a potential prediction/control failure during epileptic seizures but our results show that prediction at neuronal level can already be extremely complicated. We proceed to look at the next scale in our analysis and move from single neurons to clusters/regions of neurons.

Local Data and Clusters

On the level of regions, we can start to analyze data obtained before epileptic seizures. The eight time series we use are described in detail in the Materials section. A natural extension of our previous strategy is to compute the variance for each time series using a sliding window technique and to understand the scaling laws associated with the variance on the cluster level.

Figure 3 shows the results of this computation. We plot the inverse of the variance for since this makes it easier to understand the scaling of near the seizure point at . Vertical lines are drawn for orientation purposes in Figure 3 separating a region of low variance from a high-variance regime, giving an indication where the seizure roughly occured; see also Materials. Furthermore, for we have marked several local maxima which have been found by subdividing each time series into 20 equal time intervals and checking whether the local maximum in is also a maximum for . All four plots have several important features in common:

Figure 3

The eight plots show the average channel activity (top, blue) and the average of the inverse variance (bottom, black) for the eight time series ; the horizontal axis is the time axis where the labels correspond to the sample point number.

The sliding window length corresponds to the length of the initial gap in (5000 points). The green dots mark some local maxima of which correspond to local minima of . The fitted red curves are linear and demonstrate that the variance increases near the epileptic seizure. The black dashed vertical lines are inserted for orientation purposes, separating the two regions of low and high variance.

decreases near the seizure. The scaling law seems to be given by

There are multiple local maxima and minima for before approaching the seizure point. This indicates that we should expect oscillations in statistical indicators near epileptic seizures. Remarkably, also the number of local maxima varies only slightly between to .

The last local maximum before shows that there is a period of low variance close to a seizure.

The last local maximum before is already very close to the seizure. This means that predictions could be very difficult just based on a calculation of the variance.

The next problem to consider is what types of dynamical models can reproduce the behavior we have observed from the data analysis i.e. we look for a model for the variance in clusters/regions of neurons that displays the observed oscillatory behavior and scaling law. At first glance, the dynamics in Figures 3(a)–(h) could be interpreted as a summation of voltage traces from Figure 2(b) i.e. of neurons that are (almost) in synchrony where the coherent spiking originates from the noise-induced escape of a spiral sink. However, the real problem in understanding the dynamical mechanism of epileptic seizures is shown in Figure 4 where we also plot the inverse of the variance near a critical transition. The similarities to the data in Figure 3 are clear; all four observations (A)–(D) also apply in Figure 4. The data in Figure 4 have been generated using a simple model for a Hopf critical transition [21], [30]:where are independent white noise processes that satisfy for . The model (6) was first analyzed in the context of delayed Hopf bifurcation [62], [63]. Observe that the deterministic part of the fast variables is the normal form of a (subcritical) Hopf bifurcation [64]. The slow variable can also be viewed as time since . For the simulation in Figure 3 we have chosenwith a deterministic initial condition . It is known that near a sub- or supercritical Hopf bifurcation a scaling law of the form (5) holds [30]. Obviously the scales differ between Figure 3 and Figure 4 but those can be re-scaled to match. Therefore we have found a dynamical model that could potentially explain the qualitative features of a single variance time series for a cluster of neurons.

Figure 4

Time series of the fast variable (top, blue) and the associated inverse of the variance (bottom, black) for a Hopf critical transition model (6) with parameter values 7; cf. also Figure 3.

It is very important to note that we have obtained the conjecture that a Hopf bifurcation is involved in the transition to a seizure without a detailed biophysical model. In fact, several mean-field models for various types of epileptic seizures do exhibit Hopf bifurcations [31]–[35] that form a boundary between a regular equilibrium (non-seizure) an oscillatory (seizure) regime. However, there are several mean-field models available [14] and also other bifurcation mechanisms have been identified to play a role near seizure onset [65].

Our methods also have another important implication regarding the distinction between a preictal and a proictal state [36]. From a dynamical perspective, it was suggested that one can differentiate between models that show a distinct preictal state with a parameter driving the system to a bifurcation or systems showing a proictal state where noise-induced escapes play a dominant role [6]. A subcritical Hopf bifurcation is a model that can interpolate between the two cases. Consider (6) in the following two cases:

, and : the equilibrium is a stable focus for the deterministic dynamics but it is well-known [66] that a finite-time noise-induced escape always occurs. This can be viewed as the transition beyond a basin boundary given by the unstable limit cycles [36]. If we include another (seizure-state) attractor beyond this basin boundary we can view the situation near as a “purely proictal” state. It is well-known how to calculate the probabilistic likelihood of this Hopf transition and also for many other bifurcations involving metastability [67], [68].

, : If the noise is sufficiently small then we will reach the Hopf bifurcation point with high probability [69] and our prediction method via scaling of the variance and critical transitions applies. We are in a “purely preictal” situation.

Obviously there is a continuum of possibilities in between these two situations [60], [69] depending on the scaling of noise and time scale separation. In fact, the results shown in Figure 2 illustrate the variation in such a continuum situation for the saddle-node bifurcation. A study of intermediate regimes for all bifurcations, including the Hopf bifurcations on a mean-field level, could certainly be carried out similar to the strategy employed in [30]. Let us also point out that several models have been proposed to account for this problem in the context of epileptic seizures [70]. However, these models are usually based on introducing global dynamics as well as using global measures, such as interseizure intervals, for validation. Historically similar dynamical systems attempts have been made in other disciplines, for example for multi-mode oscillations [71] in chemistry. Later on, it turned out [53] that the local mechanisms and scaling laws are much more important as they often form the truly mathematically generic [72] building blocks of the dynamics. The Hopf bifurcation normal form (6) as well as the local dynamics near the fast subsystem saddle-node bifurcation in (1) are the most generic - i.e. codimension 1 [72], [73] - phenomena available. Therefore it is absolutely necessary to investigate the link between these phenomena and epileptic seizures first as demonstrated by our scaling law results.

Correlations between clusters

In the preceding two sections we investigated neuronal dynamics at different spatial scales, from single model neurons to neurophysiological data from clusters of neurons, using the variance as a univariate measure. In the following section, we will focus on the dynamics from many clusters of neurons encompassing a larger spatial scale. In systems with spatial degrees of freedom an increase in the noise level can produce spatiotemporal order characterized by more regular activity patterns [74], [75].

In contrast to the previous, bivariate measures for the activity between different clusters will be used. Bivariate measures can take into account the correlation of two signals. Information about the correlation of neuronal activity between different anatomical regions can give insights into the state of the network as a whole. With regard to epilepsy, correlation based measures such as mean phase coherence (MPC) and maximum linear cross correlation (MLCC) have yielded promising results in identifying pre-ictal states [29], [44], [76], [77].

In this section, we will start by considering the maximum linear cross correlation for the ECoG data used in the preceding parts, reviewing and confirming some recent observations. We will then continue to extend the bivariate analysis to wavelet-based synchronization measures able to resolve pairwise correlations at different frequency bands. We will compare these results to those obtained using MLCC and MPC. Our focus is again on multiscale character of the system with the goal of identifying scaling relationships at each level of observation.

Maximum linear cross-correlation

The maximum linear cross-correlation (MLCC) quantifies the similarity between two time series and . MLCC is a linear measure of lag-synchronization which captures the normalized product of two time series dependent on a lag [43]:whereis the linear cross-correlation function. As a measure of synchronization between activity in different anatomical areas, MLCC has been proposed and successfully applied as a precursor for pre-ictal brain activity [77], [78]. We computed the MLCC of 5 randomly chosen signal pairs for each time window (5000 sampling steps, a consecutive time window being shifted 50 sampling steps forward). Figure 4 shows the average over the 5 pairs for each of the 8 time series considered in the preceding sections.

In most of the depicted time-series (patients 1, 2, 3, 4, 6, 7) an increase in the MLCC can be observed with the seizure onset (Fig. 5) which is in agreement with the general observation of increased synchronization during a seizure [37]. Prior to epileptic seizures a decrease in MLCC values has been reported and used to identify a preseizure state [77], [78]. To relate to these reports and later also compare MLCC to the wavelet-based synchronization measure (see following section), we calculated MLCC for a pre-ictal and an inter-ictal time interval. Figure 6 (left column) depicts the time series of MLCC values of patient 4 during a pre-ictal (top) and an exemplary inter-ictal interval (middle), an interval being at least 6 hours apart from the next seizure attack. Average values of MLCC are plotted left in the bottom row illustrating the comparably lower values during pre-ictal intervals. MLCC levels are lower during the pre-ictal compared to the inter-ictal interval confirming recent reports of decreased synchronization as one characteristic precursor for a seizure.

Figure 5

Maximum linear cross-correlation for eight pre-ictal time series .

Vertical lines indicate the approximate onset of the seizure attack.

Figure 6

Decrease of synchronization measures during a pre-ictal interval.

Left column: time series of maximum linear cross correlation during a pre-ictal (top) and an inter-ictal (middle) interval. Right column: time series of for three scales during a pre-ictal (top) and an inter-ictal (middle) period are depicted. Vertical dashed lines indicate the onset of the seizure attack. Averages over the first 150000 sample points of each time series indicate a distinct decrease of each synchronization measure during the pre-ictal interval (bottom row). Error bars show standard deviations.

A lot of effort has been put forward to utilize the observed synchronization drop in predicting seizure attacks, most of these works addressing the question whether it could be used to identify a preseizure state [29], [44], [45], [76]–[78]. In this work we are not addressing this issue but focus on the dynamics and scaling relations of correlation measures near the seizure onset. For this purpose we extend the analysis to wavelets able to resolve correlations between clusters for different frequency bands.

Wavelets

Wavelet analysis has been applied in neuroscience research for some time [79], [80]. Wavelet coefficients provide a frequency-dependent moving average over a time series which can be used to derive a time-resolved frequency-profile for the data given. This capacity has also been made use of in the detection of seizures [81] and the investigationen of frequency profiles of epileptic seizures in humans and animals [40], [82]. Wavelet analysis [10] can also be used as an elegant tool to identify intervals of phase synchronization (or phase-locking) between neurophysiological time series. The phase definition can thereby be used for broad-band synchronization analysis or analysis of a specific frequency of interest.

In this study, we investigated broad-band phase-locking between pairs of signals as introduced in [83]. There, the original signal is decomposed with respect to multiple scales related to frequency bands of decreasing size. To derive a scale-dependent estimate of the phase difference between two time series, we follow the approach described in [39] using Hilbert transform derived pairs of wavelet coefficients [83]. The instantaneous complex phase vector for two signals and is defined as:where denotes the -th scale of a Hilbert wavelet transform and its complex conjugate. A local mean phase difference in the frequency interval defined by the -th wavelet scale is then given bywithbeing a less noisy estimate of averaged over a brief period of time [39]. One can then identify intervals of phase-locking (PLI) as periods when is smaller than some arbitrary threshold which we set to here. Furthermore, we require the modulus squared of the complex time average, , to be greater than 0.5, limiting the analysis to phase difference estimates above this level of significance. We denote phase-locking intervals between two signals and as . To obtain a measure of frequency-specific phase-locking in a defined time window, we calculate the sum of for all pairs of signals and normalize this expression to confine the measure to the interval :where is the number of signals and the number of time steps in the time window under consideration.

We analyzed data for each patient for 3 different scales, referring to frequency bands 12–25, 6–12 and 3–6 Hz for patients 1–3, 5–8 and 16–32, 8–16 and 4–8 Hz for patient 4, respectively. The computation of was done for time windows of 5000 sampling steps, consecutive time windows were shifted forward by 50 sampling steps. Figure 7 shows the results of this computation.

Figure 7

Phase-locking measure for the eight time series .

Colors correspond to different scales. The vertical dashed lines indicate the approximate onset of the epileptic seizure attack.

In all 8 patients, comparably low values of () are observed for all scales. Similarly to MLCC synchronization as measured by the phase-locking intervals for different scales is decreased during an pre-ictal interval compared to an inter-ictal one (Fig. 6, right column). The observed decrease in synchronization measure suggests that application of could also prove useful in preseizure state detection algorithms, similar to the MLCC.

As mentioned earlier, our focus is on the dynamical behavior near the seizure onset. Aside from the aforementioned low values of , some characteristic features can be observed:

Phase-locking measured by increases around seizure onset times. (Similar to MLCC, this is seen less clearly in patients 5 and 8.)

The increase of for different scales often starts at different times.

The increase of appears to be linear.

Point A reflects the fact of increased synchronization between cortical regions observed during seizures. We observed that starts to incease at different times for different scales. Figure 8 depicts the exemplary behavior of of patient 1 near seizure onset time. Furthermore, appeared to increase linearly. Fitting a linear function close to seizure onset times provided the better fit compared to power-law or exponential relationships (Fig. 8).

Figure 8

Comparison between and mean phase coherence for patient 1.

Both measures based on phase-synchronization show a similar behavior with an increase around seizure onset time. Colored vertical lines indicate the beginning of the increase in synchronization near seizure onset (black dashed vertical line). The increase appears to be linear (grey dotted lines) and starts at different times for different scales.

Another nonlinear measure based on phase synchronization is the mean phase coherence [45], [76]. For two pairs of neurophysiological time series and it is given bywith being the phase difference between the two signals at time and denoting the average over time. We calculated for all pairs of signals using the wavelet-derived, scale-dependent phase differences for each patient. The average over all showed a similar time course as (Fig. 8). Near seizure onset, the same temporal order of the increase in synchronization was observed indicating independence from the specific measure of phase-synchronization. Direct comparison of both nonlinear synchronization measures and to MLCC suggests that the frequency resolved measures add new information at the onset of the seizure. Therefore such multiscale measures may potentially be better suited to explain the dynamical process that causes a seizure attack.

Discussion

In the present paper we aimed for a better understanding of the dynamical processes involved in seizure generation. Our approach extended over three spatial scales involving two recently developed methods (critical transitions and wavelet derived phase-lock intervals). We showed for the first time that the theory of critical transitions [21], [30] can be applied in the context of excitable neurons operating near the spiking threshold. On the level of clusters of neurons we identified a potential Hopf bifurcation as the seizure onset mechanism from data based on this theory and found a new scaling law of single-event time series data. On the largest spatial scale we observed a scaling law occurring at time-shifted onset times and compared our wavelet-based phase-locking measure to other bivariate measures.

One of our main results is the observation of scaling laws on different spatial scales – for individual neurons (4), for activity of clusters of neurons (5) and for the increase of phase-locking near the seizure onset. A recent publication highlighted five power-law scaling laws related to epileptic seizures and their analogy to earthquakes (the Gutenberg-Richter distribution of event sizes, the distribution of interevent intervals, the Omori and inverse Omori laws and the conditional waiting time until next event) [84]. Other works investigating scaling laws of ictal and interictal epochs reported similar inter-seizure-interval statistics in genetically altered rats while in human data no power-law distribution was observed [34], [85].

The observation of such scaling laws is important because it may guide new models of seizure dynamics by allowing insights into the dynamical processes that may have generated the underlying data. Many of the scaling laws reported here and elsewhere [84] exhibit power laws. The observation of similar scaling laws on different spatial scales, from single neurons to the size distribution of different seizures, strongly emphasizes the multi-level character of epileptic seizure generation. More importantly, it yields insights into the dynamical properties of the underlying system [86]. The strong analogies between seismic shocks and brain seizures have previously been pointed out and hypothesized to emerge from the structural commonality of the two systems: both are composed of interacting nonlinear threshold oscillators and are far from equilibrium [87]. Critical dynamics is believed to be a consequence of these structural properties in both these systems. Recent findings in preparations of rat cortex [88] and primate brain in vivo [89] exhibiting power-law statistics of activity, a hallmark of phase transitions [90]–[92], have led to the hypothesis that also human brain dynamics is poised at a phase transition [39], [93]. Although such statistics can result from different processes, the self-similar behavior captured by the diverse scaling laws on different levels might potentially be related to the notion of criticality in brain dynamics. Models describing epilepsy should also resemble these multi-level scaling laws and take into account critical brain dynamics.

Decomposition into different spatial scales showed oscillations in a pre-seizure state at all levels. Observation of such oscillations in real world data offers characteristics to be useful when testing future models. As we showed here, based on critical transition theory, the variance's oscillations along with its scaling law are generically characteristics for Hopf bifurcation. These results therefore validate previous seizure models assuming a Hopf bifurcation as a main mechanism as the transition point [31], [33], [34]. While the goal in seizure prediction is to predict large events, there is growing consensus about the key role played by small events, from precursor oscillations to subclinical seizures [1], [84]. Future models and predictor systems should encompass those as prediction algorithms unable to account for such small oscillations would be ill-adapted and likely provide incorrect seizure forecasts.

Near seizure onset we observed a time shifted increase in phase-locking. In a recent study, wavelet analysis of spike-wave discharges, a different form seizure activiy, revealed changes in the time-frequency dynamics during discharges. While initially a short period with the highest frequency value was observed, the frequency later decreased [40], [94]. Other studies showed high frequency oscillations specifically at seizure onset [41], [42], see [13] for a comprehensive overview. Together these studies demonstrate dynamic changes in the time-frequency domain of seizures with higher dominating frequencies at seizure onset. One could speculate that the time shifts in phase locking reported here are related to these observations suggesting a frequency-dependent, shifted start of synchronization near seizure onset.

Materials and Methods

Eight patients undergoing surgical treatment for intractable epilepsy participated in the study. Patients underwent a craniotomy for subdural placement of electrode grids and strips followed by continuous video and electrocorticogram (ECoG) monitoring to localize epileptogenic zones. Solely clinical considerations determined the placement of electrodes and the duration of monitoring. All patients provided informed written consent. The study protocols were approved by the Ethics Committee of the Technical University Dresden. ECoG signals were recorded by the clinical EEG system (epas 128, Natus Medical Incorporated) and bandpass filtered between Hz and Hz. Data were continuously sampled at a frequency of Hz (patients and ) and Hz (patient , [95]) with two electrodes used as reference. We always indicate the sampling point number on the time axis if we use the data. No claims regarding a large-scale statistical validity of the data set is made since the total patient sample size is rather small. Although this is an important issue [29] we focus here on identifying the dynamical mechanisms and new time series analysis techniques in the context of epileptic seizures. Furthermore, we also do not claim that the vertical lines we use in the plots of the data indicate exact seizure onset as determined by neurophysiologists or direct monitoring of patient symptoms.