Cerebral Cortex Advance Access originally published online on November 9, 2005
Cerebral Cortex 2006 16(9):1296-1313; doi:10.1093/cercor/bhj072
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
A Unifying Explanation of Primary Generalized Seizures Through Nonlinear Brain Modeling and Bifurcation Analysis
1 School of Physics, University of Sydney, NSW 2006, Australia, 2 Brain Dynamics Centre, Westmead Hospital and University of Sydney, Westmead, NSW 2145, Australia, 3 School of Psychiatry, University of New South Wales, Randwick, NSW 2031, Australia, 4 Department of Mathematical Sciences, Loughborough University, Loughborough, LE 113TU UK, 5 Department of Neurology, Westmead Hospital, Westmead, NSW 2145, Australia, and 6 The Black Dog Institute, Hospital Road, Randwick, NSW 2031, Australia
Address correspondence to email: mbreak{at}unsw.edu.au.
 |
Abstract |
|---|
 Top Abstract IntroductionMethodsResultsConclusionAppendixReferences |
|---|
The aim of this paper is to explain critical features of the human primary generalized epilepsies by investigating the dynamical bifurcations of a nonlinear model of the brain's mean field dynamics. The model treats the cortex as a medium for the propagation of waves of electrical activity, incorporating key physiological processes such as propagation delays, membrane physiology, and corticothalamic feedback. Previous analyses have demonstrated its descriptive validity in a wide range of healthy states and yielded specific predictions with regards to seizure phenomena. We show that mapping the structure of the nonlinear bifurcation set predicts a number of crucial dynamic processes, including the onset of periodic and chaotic dynamics as well as multistability. Quantitative study of electrophysiological data supports the validity of these predictions. Hence, we argue that the core electrophysiological and cognitive differences between tonic–clonic and absence seizures are predicted and interrelated by the global bifurcation diagram of the model's dynamics. The present study is the first to present a unifying explanation of these generalized seizures using the bifurcation analysis of a dynamical model of the brain.
Key Words: bifurcation • neural modeling • nonlinear dynamics • primary generalized epilepsy • time series analysis
|
Introduction |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
Primary generalized seizures are pathological brain rhythms that, by definition, involve all cortical regions and are associated with a gross disruption of cognitive activity. Absence (Petit Mal) and tonic–clonic (Grand Mal) seizures are the 2 main generalized seizures in humans. Several features critically distinguish these two types of seizures. Tonic-clonic seizures are associated with markedly different pre- and postictal electroencephalographs (EEG), lengthy duration, and dynamically evolving waveforms that occur within each seizure. Although subjects are typically awake and conscious prior to the seizure, they are invariably unconscious postictally. In contrast, absence seizures have a brief on–off quality, similar pre- and postictal EEG, and a well-structured periodic spike and slow-wave shape that slows only slightly during the seizure but does not significantly alter in its morphology. Remarkably, cognitive function is only minimally disrupted after the seizure. The aim of this paper is to study the mechanisms underlying the onset, evolution, and offset of these seizures using a physiologically motivated model of the brain's dynamics. More specifically, we test and extend a number of specific predictions (Robinson and others 2002
) based on the premise that generalized seizures represent a transition from stable linear dynamics to unstable nonlinear behavior by studying the nature of the model's bifurcations. In doing so, we seek a formal and unified explanation of the generalized seizures.
A bifurcation is a sudden change in a dynamical system's activity, such as from steady-state to periodic behavior. The transition from laminar to turbulent fluid flow is a well-known physical example. Nonlinear instabilities and bifurcations in large-scale neural activity may be of special significance to brain dynamics. Depending on the context, timing, and extent, such phenomenon may be either adaptive—allowing flexible switches in cognitive set (Wright and others 1985; Friston 2000; Breakspear 2002; Freeman and Rogers 2002; Breakspear and others 2003; Breakspear and others 2004) and behavior (Kelso and others 1992; Daffertshofer and others 2000; Fuchs and others 2000)—or disruptive, such as that at the onset of a generalized seizure (Arnhold and others 1999; Stam and van Dijk 2002). Understanding such nonlinear instabilities may provide a unique window into the nature of neurophysiological processes occurring in neural systems because they denote particular types of dynamical processes (Crevier and Meister 1998; Izhikevich 2001), such as those with strong feedback and nonlinearity. Finally, a bifurcation from linear to nonlinear brain dynamics would render any data analysis method grounded in stochastic linear theory problematic. Hence, the elucidation of bifurcations in large-scale neuronal systems has cognitive, physiological, and methodological significance for neuroscience.
Bifurcations occur in a variety of forms, such as from linear (steady state) to nonlinear dynamics or from one type of nonlinear oscillation to another (Abraham and Shaw 1992). The former can occur if the strength of a feedback process rises above a critical value. Several empirical studies have concluded that large-scale electrical activity in the healthy human brain is, indeed, predominantly a linear/stochastic phenomenon with intermittent instances of weakly nonlinear fluctuations in the alpha frequency (8–13 Hz) range (e.g., Stam and others 1999; Breakspear and Terry 2002). The nonlinear neural model employed in the present study is able to predict and explain a variety of resting- and sleeping-state EEG when evolving in a stable, weakly damped linear regime (Robinson and others 2001, 2004). This view (of only occasional and weak nonlinearities) also finds strong support in the success of classic functional neuroscience algorithms that are rooted in a stochastic/linear framework (Friston and others 1994). In contrast, clinical research suggests that several pathological processes, such as seizures (Andrezejak and others 2001) and abnormal rhythms in pathological states (Stam and Pritchard 1997; Stam and others 1997) have a strong nonlinear component. In Robinson and others (2002) it was proposed that the transition from resting-state EEG to seizure activity may be viewed as a bifurcation from linear to nonlinear oscillations in the brain's electric activity, whereas different types of seizures may be viewed as bifurcations between distinct types of nonlinear dynamics (see also Wendling and others 2002; Lopes da Silva and others 2003; Perez Valezquez and others 2003).
The present paper undertakes to formally study this hypothesis in a physiologically based model of brain activity (Robinson and others 2001). Two instabilities are studied—one at approximately 3 Hz and the other within the alpha (10 Hz) rhythm. These represent the respective frequencies at which absence and tonic–clonic seizures are initiated. The behavior of the model at and beyond such instabilities is compared with 3 EEG data sets—1 of young subjects with absence seizures, 1 with tonic–clonic seizures, and 1 of resting-state healthy adult subjects. The first 2 data sets allow comparison between modeled and experimental seizure data. The latter data have been previously shown to be associated with a weak nonlinear structure (Stam and others 1999; Breakspear and Terry 2002) and hence permit analysis of the model in a stable, resting state—although in the vicinity of an instability. The local bifurcation diagrams are studied for each instability and then compared with the observed phenomena. The global bifurcation diagram—obtained by examining the structure of the bifurcation set under changes in multiple parameters—permits local features of each instability to be understood from a global and unifying perspective. We hence investigate whether the topology of the global bifurcation diagram explains the essential difference between absence and tonic–clonic seizures by interrelating them within this broader framework.
|
Methods |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
There are 3 components to the methodology: 1) The nonlinear corticothalamic model that forms the basis for the bifurcation and time series analysis. A description of this model is given below and follows Robinson and others (2001
, 2002). 2) Scalp EEG data—which permits empirical testing of the predictions generated by the model—is then described. 3) Nonlinear techniques are employed in order to study numerical data generated from the model and the observed scalp EEG data. These are described last. Corticothalamic Brain Model
Large-scale neural activity arises from interactions between several neural populations, notably excitatory and inhibitory cortical neurons and specific subcortical nuclei such as the thalamus. The corticothalamic model studied here is based on the evolution of several dynamical variables within each of these populations. The variables represent the local mean value of a physiological process at position r in these neural systems, averaged over a small patch (
0.3 mm) of surrounding neuropil. Hence, this model belongs to the class of "lumped" or "mean field" neural models (e.g., Nunez 1974; Freeman 1975; Jirsa and Haken 1996; Robinson and others 1997).
General (Spatially Continuous) Model
We first describe the model in its general form. The dynamical variables within each neural population are the local mean cell-body potentials Va, the mean rate of firing at the cell-body Qa, and the propagating axonal fields 
a. The subscript a refers to the neural population (e = excitatory cortical; i = inhibitory cortical; s = specific thalamic nucleus; r = thalamic reticular nucleus; n = nonspecific subcortical noise). The firing rates Qa are related to the potentials Va according to the sigmoid activation function Qa(r, t) = S[Va(r, t)], where S is a smooth sigmoidal function that increases from 0 to Qmax as Va increases from –
to . We model S as
 | (1) |
is the mean neural firing threshold, 
is the standard deviation of this (approximately normally distributed) threshold, and Qmax is the maximum firing rate.
In each neural population, firing rates Qa are propagated outward as fields a according to the damped wave equation
 | (2) |
 | (3) |
The parameter 
a = va/ra—where ra and va are the characteristic range and conduction velocity of axons of type a—governs the dispersion of propagating waves, and 
2 is the Laplacian operator (the second spatial derivative). The system of equations is closed by introducing the effect of incoming axonal inputs to neurons at r from other neural populations. The cell-body potential Va results after postsynaptic potentials have been filtered in the dendritic tree and then summed. For excitatory and inhibitory neurons within the cortex, this is modeled using a second-order delay-differential equation (Robinson and others 2001)
 | (4) |
 | (5) |
and ß are the inverse decay and rise times, respectively, of the cell-body potential produced by an impulse at a dendritic synapse. Note that input from the thalamus to the cortex is delayed in equation (4) by a propagation time t0/2. For neurons within the specific and reticular nuclei of the thalamus, it is the input from the cortex that is time delayed and hence
 | (6) |
The default values of all parameters are given in Table 1. All parameter values chosen in this study have been used in previously published studies (Robinson and others 2002), emphasizing that the results of the numerical analysis are used in a predictive rather than exploratory vein in relationship to the EEG data. A schematic representation of the main populations and loops in the model is presented in Figure 1.
|
|
Setting all spatial and temporal derivatives in equations (1)–(6) to 0 determines global (spatially invariant) corticothalamic steady states. Small perturbations around these states (representing noisy influx from the brainstem) obey a linear wave equation, the study of which has been employed to explain resting-state EEG temporal (Robinson and others 1997
, 2001, 2004) and spatial (O'Connor and others 2002) spectra. Global (Spatially Invariant) Model
A full nonlinear analysis of this system is a highly nontrivial task. However, in certain circumstances, particularly in the study of generalized seizures, brain activity may be dominated by very large scale—or even whole-brain—processes, in which case the dynamical variables may not depend greatly on spatial position r. This "global model" can be studied by setting the spatial gradient term in equation (3) to 0, yielding
 | (7) |
The variables Va, Qa, and a now depend solely on t and not r. The smallness of ri also lets us set i 
, yielding a set of 8 first-order delay-differential equations. These equations permit a computationally parsimonious method of studying large-scale brain dynamics and are employed for the majority of the present study. They are given in the Appendix (First-Order Delay-Differential Equations for the Spatially Uniform Model).
Studying the linear stability criteria for this system permits mapping of the boundary that marks the transition between steady-state behavior and nonlinear oscillations. Previous exploration of the model for realistic parameter ranges revealed a small number of key instabilities that constrain the way nonlinear oscillations may arise (Robinson and others 2002). As well as the 3-Hz and alpha instability studied in the present study, slow-wave (<1 Hz) and spindle instabilities (approximately 12 Hz) may also arise. The occurrence of only a small number of instabilities suggests that it may also be possible to study the dynamics and stability of the brain in a phase space of low dimensionality. Indeed, formal analysis of low-frequency instabilities suggests that 3 variables x, y, and z—parameterizing corticocortical, corticothalamic, and intrathalamic instability—capture the parameter combinations at which the brain model loses stability (Robinson and others 2002). These combine dynamic variables and state parameters and are given in the Appendix (Reduced Parameter Space). They enable visualization of where in parameter space the model becomes unstable. In Figure 2 the instability boundary within this truncated space is presented. The zone within the tent-shaped region is associated with stability of the model, whereas points outside are associated with nonlinear oscillations. The boundaries of the tent-shaped region correspond to the points at which the fixed points lose stability in a bifurcation diagram. The latter also contain information concerning transitions between different types of nonlinear oscillations.
|
In the present study, the general model is only required in the section dealing with weak nonlinear interdependence in the resting-state data. The spatially invariant model is employed to study the onset of the generalized seizures. Numerical integration was performed using a fourth-order Runge–Kutta integrator. A cubic-spline interpolator was employed in order to estimate the time-delayed values of the midpoints required for the Runge–Kutta algorithm.
Scalp EEG Data
Three EEG data sets are studied. All are scalp EEG data with electrode placement following the 10–20 international system and linked earlobe reference. Ethics approval was obtained prior to data collection for each data set, according to the Ethics Committee at each institution.
- "Healthy human EEG alpha data" collected from 40 adults (age 20–54 years) who disavowed psychiatric or neurological illness at Westmead Hospital, Sydney. Skin resistance at each site was <5 k
. Data were collected at a rate of 250 Hz through a SynAmps amplifier and filtered with a 50-Hz, low-pass, third-order Butterworth filter. Artifacts caused by eye movement were corrected offline according to the method of Gratton and others (1983). Data were collected from each subject during 130 s of a resting eyes-open paradigm and 130 s during a resting eyes-closed paradigm.
- "Absence seizure data" drawn from a database of 13 adolescent epileptic patients from the Department of Neurology, Westmead Hospital, Sydney. Data were collected at a rate of 200 Hz and filtered with a 70-Hz low-pass filter.
- "Generalized tonic–clonic seizure data" drawn from a database of 7 epileptic patients (Quian Quiroga and others 2002). These patients were admitted to an epilepsy-monitoring unit with diagnosis of pharmacoresistant epilepsy and no other accompanying disorders. Antiepileptic drugs were gradually tapered after admission. Each signal was digitized at 409.6 Hz through a 12-bit A/D converter and filtered with an antialiasing 8-pole low-pass Bessel filter with a cutoff frequency of 50 Hz. The signal was digitally filtered with a 1- to 50-Hz bandwidth filter and stored at 102.4 Hz. The scalp recording was measured bipolarly from the T4–T6 (right temporal) locations.
Nonlinear Data Analysis
EEG Data
In order to test whether seizure activity is associated with nonlinearity—implicit within the hypothesis that a seizure occurs after a nonlinear bifurcation—a nonlinear prediction algorithm (Casdagli 1989) was employed. Briefly, the data are divided into discrete time windows. In each window a "nonlinear prediction error" is calculated. This error reflects the ability of a local nonlinear model to predict the amount of uncertainty in the data. Low errors indicate a good fit and hence a nonlinear structure. Using a bootstrap or resampling scheme (applied to the original data), an ensemble of prediction errors is then calculated to represent the null hypothesis that the values of the errors are due to purely linear correlations within the data (Theiler and others 1992). The data are said to contain a nonlinear structure if the observed (experimentally derived) prediction error lies outside the distribution of these "surrogate" errors (for further details, see Terry and Breakspear 2003). Nineteen surrogates were constructed to allow for nonparametric statistical inference at 95% confidence within each window. For graphical clarity, the inverse of the prediction errors—which we term the "nonlinear prediction indices"—is plotted.
A modification to this algorithm allows multichannel EEG data to be tested for evidence of nonlinear interdependence—a nonlinear equivalent to the coherence function (Schiff and others 1997; Terry and Breakspear 2003). A multivariate surrogate algorithm (Prichard and Theiler 1994; Rombouts and others 1995) is then employed to exclude the contributions of purely linear correlations. This can be used in conjunction with the general version of the corticothalamic model to estimate the predicted nonlinear interdependences.
Numerical Data
Bifurcation diagrams illustrate the asymptotic behavior of dynamical systems across a range of parameter values and, hence, enable one to visualize sudden changes in a dynamical system subsequent to incremental changes in a state parameter. In order to construct these, numerical integration was performed across a varying range of vse—the excitatory influence of cortical pyramidal cells on the specific thalamic nuclei. This parameter was chosen because of its simple physical meaning and the prior implication of excitatory corticothalamic feedback in the pathophysiology of generalized tonic–clonic (Zifkin and Dravet 1997; McCormick and Contreras 2001) and absence (Destexhe and Sejnowski 2001; McCormick and Contreras 2001; Meeren and others 2002) seizures. Local maximums and minimums of the resulting numerical time series are plotted for each parameter value. For the time series comparison of the model and EEG data, we plot the macroscopic excitatory fields e as these best represent the cortical correlate of scalp potentials up to a linear transformation of the amplitudes (Nunez 1995). Specifically, scalp potential is proportional to the cortical potential, which is proportional to the mean cellular membrane currents, which are in turn proportional to the firing rates. Hence, apart from a (dimensional) constant of proportionality and the effects of volume condition, scalp EEG signals correspond closely to e (Robinson and others 2004). Parameter values are given in Table 1. To better simulate a real physiological system, small amplitude (signal–noise ratio = 0.90) autocorrelated stochastic terms were added to the parameters (system noise) for the numerical time series plots. Such noise was generated according to
 | (8) |
by
 | (9) |
Such a model of parameter noise represents the simplest method of simulating an autocorrelated stochastic process. It should be noted that such noise is only used in the time series plots and has no impact on the calculation of the model's bifurcation diagram.
|
Results |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
The analysis of the results is presented in 3 sections: 1) The bifurcation occurring at the 3-Hz instability is first studied as a model for absence seizures. The analysis of the model is given first, permitting prediction of the electroencephalographic data that follow. 2) The alpha instability is then presented. In order to underline that the model predicts the emergence of seizure activity from resting-state data and because previously published data are available, we first present an analysis of the model very close to the alpha instability—a "weak" instability. Subsequently, we present the full bifurcation of the global model to predict generalized tonic–clonic seizures. As with the 3-Hz instability, the model and then the EEG data are studied. 3) The global bifurcation is then presented in order to permit potential unification of the 3-Hz and alpha instabilities.
Absence Seizure: Corticothalamic Model
We first study the nonlinear instability occurring at approximately 3 Hz. This is achieved by choosing physiologically plausible parameters that place the system in the vicinity of a weak 3-Hz instability (Robinson and others 2002) as given in Table 1. As stated above, vse is then varied in order to study the geometry of the bifurcation set.
The results are presented in Figure 3a–j. The bifurcation diagram—panel (a)—exhibits a Hopf bifurcation to periodic dynamics with an initial (supercritical) instability at vse 1.8 x 10–3 V s and further period-doubling instabilities at vse 3.4 x 10–3 V s and vse 4.2 x 10–3 V s. In this noise-free plot, it can be seen that only periodic oscillations occur. An exemplar time series, with added system and measurement noise is given in panel (b). This was created by dynamically ramping vse from the linearly stable (weakly damped region) upward into the region of linear instability (panel c). The dashed lines in panel (a) show the extreme values of the ramp function. Close-up images of the onset (panel d) and offset (panel e) of the seizure exhibit a number of key phenomena: 1) Shortly after the onset of ramp-up of vse at t = 5 s periodic oscillations of growing amplitude appear in the field potentials. These occur as the system passes through the periodic regime in the bifurcation plot. 2) After 2–3 cycles, spike and slow-wave oscillations appear, heralding the onset of period-doubling limit-cycle oscillations. These continue throughout the seizure, although their amplitudes are modulated by the combined effects of system and measurement noise. 3) During the ramping down of vse at t = 19.5 s, the amplitude of the spike and wave oscillations diminish. The spikes disappear at approximately t = 20.5 s as the system passes through the simple period 1 regime. 4) Finally, the remaining oscillations are damped away, and the system returns directly to the same pre-ictal EEG state, governed by stable damped stochastic fluctuations. This is reflected in the spectrum of the seizure (panel f), showing similar pre- and postictal spectra. The seizure spectrum is dominated by the 3-Hz spike and wave oscillation and its harmonics.
|
A (3-dimensional) time-delay–embedded phase portrait of a simulated (noise-free) seizure is shown in panel (g). Gray arrows indicate the flow of orbits away from the unstable fixed point subsequent to the bifurcation and onto the limit-cycle attractor. The growth in amplitude of the spike (at the far side of the attractor) can be seen as the orbits spiral outward. The orbits follow the same unstable manifold (but spiraling downward) at the conclusion of the seizure (not shown). Panel (h) depicts the morphology of the seizure in the space spanned by the corticocortical, corticothalamic, and intrathalamic "stability" variables, x, y, and z (as given in the Appendix, Reduced Parameter Space). As expected, the seizure is located outside the tent-shaped stability zone. Because the seizure corresponds to limit-cycle dynamics, it can be embedded within this 3 dimensional space without crossing of the orbits (i.e., with uniqueness of the solution curves). This indicates that, during such a seizure, a dynamical system of reduced dimensionality should be able to sufficiently describe the macroscopic neural dynamics, or equivalently, a relatively small number of physiological processes may be responsible for the onset and maintenance of the seizure activity.
Panel (i) shows the 3 principal fields a (a = e, r, and s) during the model seizure, normalized by their own means, to illustrate the underlying mechanism from which the waveform is generated. We choose the peak of s as our starting point (t 10.0 s) for convenience. It must be stressed that this is an arbitrary choice because the periodic nature of the dynamics imply that the waveforms have no particular beginning or end but are an emergent feature of the entire corticothalamic system. When s (dashed line) is at its highest value, e (solid) is at an intermediate value and r (dotted) is low. The thalamus is in its positive-feedback, excitatory state. The reticular nucleus responds immediately to the high activity in the specific nuclei and in addition to the incoming signal leaving the cortex time t0/2 earlier. This sudden increase in r acts to suppress s rapidly, and the thalamus switches to its negative-feedback, inhibitory state (t 10.05 s). At time t0/2 later, the peak in activity in the specific nuclei reaches the cortex and e rises accordingly, with the peak slightly broadened by synaptic rise and decay times. Another t0/2 later, the spike in activity in the cortex reaches the thalamus, exciting the reticular nucleus and hence further suppressing the specific nuclei (t 10.10 s). This represents negative corticothalamic feedback, and the response of r results in a further broadening of the signal. With near-silence in the specific nuclei, the cortical neurons relax. At t0/2 later, there are no inputs to the reticular nucleus and hence it too relaxes, leading to a period of near-silence in all 3 fields (t 10.20 s). In the meantime, the specific nuclei have been receiving the external input stimulus n, which is sufficient to cause the specific nuclei to reactivate once the reticular nucleus has been suppressed, and so s rises. This corresponds to the thalamus switching back to its positive-feedback state. Immediately, r responds slightly, enough to subdue the growth of s, but not before the increased activity is substantial enough to excite e a time t0/2 later to the original intermediate value. Meanwhile, the specific nuclei have been receiving external stimuli continuously, and so while r is still low, they are able to fire, completing a full cycle (t 10.35 s). The delicate balance of this process highlights the sensitivity of the shape of the waveform to changes in the various parameters. Indeed, a variety of spike and wave and polyspike morphologies are possible under modest changes in parameters. This may explain the varied morphologies observed in clinical absence seizure EEG recordings. It is also consistent with a previous study of absence seizures employing autoregressive nonlinear signal analysis methods that showed that similar combinations of time-delayed nonlinear (quadratic) terms could explain a variety of spike and wave morphologies (Schiff, Victor, Canel, and Labar 1995).
For completeness, the cortical inhibitory field i (dot-dashed) is plotted with the excitatory field e (solid) in panel (j). Note that the dynamics of i closely match those of e, reflecting the net effect of cortical excitation and inhibition. Recall that i does not project onto the thalamus, and hence, it only affects s and r indirectly through its effect on e.
Absence Seizure: Scalp EEG Data
The first step of the experimental EEG data analysis was to test the key premise of the paper, namely, that seizures correspond to bifurcations from linearly stable to nonlinear oscillations. As discussed above, this was achieved by using a measure of nonlinear predictability and comparing EEG with phase-randomized surrogate data. An exemplar absence seizure (Fz electrode) is shown in Figure 4a. In Figure 4b, the normalized predictability index is given as the solid line. Plots obtained from 19 surrogate sets are dashed lines. It can be seen that the occurrence of the seizure is coincident with a sudden and large increase in this index of nonlinear structure, consistent with the appearance of nonlinear oscillations. Also noteworthy is that the pre- and postictal EEGs are associated with intermittent and weak nonlinearity, evident as occasional increases in the nonlinear predictability of the real compared with surrogate data (arrows). This is consistent with noisy perturbations of a weakly damped nonlinear system. In other words, pre- and postictal states represent a system close to a bifurcation. Comparable results were observed for all absence seizures studied.
|
Figure 4c,d illustrates an example of a numerically simulated absence seizure, integrated over a comparable time frame to the real seizure. Because the seizure was generated by pushing the system through nonlinear period-doubling bifurcations, it is not surprising that a coincident increase in nonlinear structure is seen. Prior to and following the seizure, the system has been set just below the Hopf bifurcation (see Fig. 3a). Hence, the occasional slight increases in nonlinear structure in the model over the surrogate data time series (arrows) during these times are to be expected, when fluctuations toward the bifurcation occur.
The nature of the nonlinear oscillations in the absence seizure is studied in Figure 5a–f. Panel (a) shows a complete seizure, revealing an approximately symmetrical appearance. The spectrum of this seizure (panel b) is dominated by the 3-Hz spike and wave morphology. Notably, the postictal spectrum returns rapidly to its relatively featureless preictal form. In panels (c) and (d), the onset and offset of simple periodic, then spike and wave oscillations are clearly visible. A time-delay–embedded reconstruction of this seizure is given in panel (e). The seizure onset (t = 3–5.8 s) is given in blue and the remainder of the seizure (t = 5.8–19.5 s) in black. This plot directly illustrates the outward periodic spiral (blue) of the system onto a large-amplitude oscillatory "attractor," comparable with Figure 3g.
|
It finally remains to determine the nature of the oscillations during the clinical seizure. We refrain from calculating dynamic "invariants" such as Lyapunov spectra as these are notoriously unreliable in short, noisy time series data (Dammig and Mitschke 1993
). Rather, we study the first-return map, which is an assumption-free method of visualizing the nature of oscillatory dynamics (Perez Valezquez and others 2003). Briefly, level crossings of a dynamical variable are noted. Denoting the temporal period between successive crossings as Tn, the graph of Tn against Tn–1 gives a truncated representation of the full dynamics. A plot of the seizure in panel (a) is given in panel (f). For ease of interpretation, we plot fn = 1/Tn—the pointwise frequency of the system. The first feature to be noted is that the points fall close to the line fn = fn–1. That is, the system is nearly periodic (the level crossing was chosen below the spike so that in effect we produce a second-return plot—i.e., a true full period of the oscillation). As indicated, the frequency starts above 3 Hz and quickly falls to approximately 2.6 Hz. It then varies around 2.6 Hz in a manner without any obvious geometrical structure—that is, not confined to any obvious low-dimensional invariant. We interpret this as the noisy modulation of a fixed point in this plot—that is, a noisy limit cycle of the full dynamics. This is consistent with the findings of a purely empirical study using a different methodology (Feucht and others 1998). In summary, the evolution of this seizure quantitatively matches that of the numerical seizure generated by the corticothalamic model in a number of crucial aspects, namely, the period-doubling manner of onset and offset, the similar pre- and postictal spectra, and the periodic (2) nature of the full seizure. The bifurcation plot of Figure 3a explains these phenomena.
Weak Nonlinear Alpha Instability
Previous analysis of scalp EEG data revealed an increase in amplitude and sharpness of the alpha peak during the instances of nonlinear interdependence (Breakspear and Terry 2002). Empirical studies alone, however, are unable to test competing explanations for this phenomenon, whereas the current model permits its explicit investigation. The system parameters are first chosen to yield resting-state eyes-closed dynamics, using previously published parameters (Robinson and others 2002), to yield a strongly damped, stable state. Increasing vse toward (but not beyond) the linear instability boundary yields a very weakly damped, marginally stable system. Comparing the spectral plots of the model in these 2 states with previously published empirical data investigates whether the resting brain functions in the vicinity of a 10-Hz nonlinear "corticothalamic" instability. The results are presented in Figure 6a–h.
|
The empirical results from the database of healthy subjects' EEG recordings are presented in panels (a)–(d). Panel (a) shows a pair of bipolar derivations (O1–P3 and O2–P4) that exhibit strong nonlinear interdependence according to the nonlinear prediction algorithm discussed above. Panel (b) shows data that do not exhibit such a property. Note the "cleaner", higher amplitude alpha oscillations of (a) compared with (b). Panels (c) and (d) show the linear and nonlinear properties of these epochs: in (c) are shown the cross-spectral density plots of all epochs showing nonlinear interdependence (solid) compared with those that do not (dashed). The increased amplitude and sharpness of the alpha peak of the "nonlinear" epochs are clearly visible. The mean cross-spectral power values (±standard error) for the nonlinear epochs at frequencies 7, 8, 9, and 10 Hz are 1210 (±200), 6520 (±1500), 17 600 (±4200), and 6010 (±1400) µV2 Hz–1, respectively. The mean power values of the "linear" epochs at the corresponding frequencies are 840, 1960, 4500, and 3920 µV2 Hz–1. Hence, the power of the nonlinear epochs is most strongly and significantly increased at 8 and 9 Hz. There is also a secondary peak at the first harmonic of this frequency, 18 Hz, evident principally in the nonlinear epochs. It should be recalled that although the cross-spectra is a purely linear measure, a system that is only very weakly damped (and driven by noise) typically exhibits an increase in power and sharpness at the least damped frequency and its harmonics, consistent with the results shown in this panel.
Panel (d) shows the "nonlinear interdependence" prediction errors H plotted against the length of the forward prediction iteration H, for the EEG data (solid lines with crosses) compared with an ensemble of 19 surrogate data (dashed line). Lower prediction errors correspond to stronger nonlinear interdependence. The plain solid line denotes the boundary of the null (purely linear) distribution, defined as the minimum of the surrogate data. The data pair in panel (a) yielded the lower prediction curve, consistent with strong nonlinear interdependence. The pair in (b) crosses into the null distribution after 6 time steps and is hence classified as containing only linear cross-dependence (as determined by comparing the root mean square of the 20 prediction errors of the measured and surrogate data).
The corresponding plots for the corticothalamic model are shown in panels (e)–(h). Setting vse = 10.4 x 10–4 V s places the system close to the alpha instability (see below) and yielded data as shown in panel (e), which exhibits the same waveform as the "nonlinear" scalp data of panel (a). In comparison, vse = 10.1 x 10–4 V s places the system further from the instability, hence ensuring strong damping of any nonlinear perturbations. This yields the noisier data of panel (f). Comparing the linear cross-spectra of these data (panel g) reveals similar increased amplitude and sharpening of the alpha peak derived from the weakly nonlinear data (solid) compared with the strongly damped data (dashed). The peaks at higher harmonics in the numeric data are even more obvious than those observed in the experimental data. The nonlinear interdependence prediction plots are given in panel (h). The "weakly nonlinear" time series in (e) clearly exhibits nonlinear interdependence (solid line with crosses) as the prediction errors are consistently smaller than their surrogate counterparts. The strongly damped epoch (solid line with stars) yields prediction errors well within the null distribution.
Tonic–Clonic Seizure: Corticothalamic Model
We now study the nonlinear dynamics that follow the instability occurring at approximately 10 Hz. This is achieved by increasing vse from that studied in the last section. The results are presented in Figure 7a–j. Bifurcation diagrams are given for both the macroscopic fields e (panel a) and the pyramidal cell potential Ve (panel b) as the latter reveals additional dynamical structure. Both these show a region of bistability from vse 
9.5 x 10–4 V s to vse 10.35 x 10–4 V s. To illustrate this, we first plotted the bifurcation from left to right in black and then from right to left in red. A "continuation technique" is employed whereby the final system values from a run with the previous parameter value are used as the initial conditions when that parameter is changed (either up or down). Hence, the region of bistability is revealed when both red and black solutions can be viewed separately. Otherwise black solutions are overlaid by red.
|
What are the implications of this bistable region for the dynamics? An exemplar numerical time series is given in (c). The corresponding values of vse are given in (d). The system is initialized on the linearly stable (black) arm within the bistable zone, close to the linear instability vse
10.2 x 10–4 V s. Next, vse is increased to vse = 10.45 x 10–4 V s. The system bifurcates from the fixed point at vse 10.3 x 10–4 V s (arrow 1). At this point, the fixed point becomes an unstable spiral and the orbits hence grow exponentially in amplitude (panel f), reaching the large-amplitude attractor seen to the right of the instability. The structure of Ve reveals that this is a period 6 attractor. If vse is then decreased (arrow 2), the system does not become immediately unstable, but rather tracks back through parameter space on this large-amplitude attractor. This explains the ongoing existence of large-amplitude oscillations in the time series. When vse < 9.8 x 10–4 V s, the attractor does again become unstable and the orbits collapse back onto the fixed point (arrow 3). Although the system is back on the fixed-point attractor, it has reached this state in quite a different (more strongly damped) region than its initial configuration: in order to return to the original state, it is necessary for vse to be increased again to vse = 10.25 x 10–4 V s. Thus, the system exhibits hysteresis. The above "loop" through phase space creates a distinctive fingerprint in the spectral plot (panel e). Prior to the onset of the seizure, the spectrum displays the characteristic alpha peak of resting-state EEG. The power within this peak increases after the bifurcation has been passed, and as the orbits grow exponentially toward the large-amplitude attractor. This power is particularly strongly expressed during the seizure. Due to the nonlinear nature of the oscillations, higher order harmonics (at 20, 30, 40 Hz, etc.) are evident. Due to the bistability, this pattern remains evident even when vse drops below its initial value (panel d). Once vse < 9.8 x 10–4 V s, the attractor loses stability and the system returns to the stable linear regime. However, because of the strong nonlinear damping, the overall spectral power is much diminished. Note that power across all frequencies is lower and the alpha peak is absent. The alpha peak returns when vse is restored to its initial value of 10.25 x 10–4 V s.
Figure 7g shows the large-amplitude attractor (in black) that occurs when vse is increased past the bifurcation point. The red orbits show the exponential growth in amplitude toward this attractor once the fixed point has become unstable. In (h) the attractor is depicted again in relationship to the tent-shaped stability zone in the truncated space spanned by the stability variables, x, y, and z. For this panel, all (system and parameter) noise has been withheld, revealing the underlying periodic attractor. An additional feature of the spectral plot (e) during the seizure, most notable at higher frequencies (e.g., between 40 and 50 Hz), is the existence of subharmonics at approximately 1/3 and 2/3 of the fundamental frequency. These reflect the period 6 character of the seizure attractor. Interestingly, similar subharmonics are also visible in the spectral plots of an experimental seizure, shown in Figure 1 of Schiff and others (2000).
Figure 7i shows the dynamics of the 3 principle fields a (a = e, r, and s) during the modeled seizure. As in the analysis of the absence seizure, we arbitrarily choose the peak s as a starting point (e.g., t 9.1 s). Whereas s is at its maximum, e is descending. The high activity in the specific nuclei—in addition to the high e that had occurred time t0/2 earlier (t 9.05 s)—induces a rise in the reticular nucleus r. The subsequent peak in r—which is broadened by synaptic rise and decay times—in turn suppresses the specific nuclei. The spike in the relay nuclei reaches the cortex t0/2 later (t 9.15 s), eliciting a peak in e. Now, from equations (2) and (4) (for a = e), it may be shown that in the linear regime (a reasonable approximation for this brief portion of the dynamics) Ve undergoes a damped oscillation with approximate frequency 
where is the linearized sigmoid slope. Because the lagged value of s is negligible for the duration of the oscillation, this is a purely cortical effect, as expected, because the frequency is too high to involve the corticothalamic loop. Analysis of equation (2) (for a = e) in the linear regime suggests that for 
>> e, long-range cortical axonal damping filters out frequencies as –4, and thus the fast oscillation manifests itself in e as only a slight bump. This bias toward lower frequencies also explains the smearing of e compared with the sharp peak in Ve that drives it. Another t0/2 after the cortex is excited, the e activity reactivates the thalamus and the cycle continues.
Panel (j) shows the excitatory (solid) and inhibitory (dot-dashed) cortical activities for this same seizure. In comparison with the absence case, it can be seen that the cortical inhibitory field i is not closely tied to the excitatory field e. Also apparent is that the high-frequency (75 Hz) damped oscillation in Ve discussed above is more strongly evident in i than e. Such observations permit explicit comparison with experimental data in future studies and show the predictive power of the modeling approach.
Tonic–Clonic Seizure: Scalp EEG Data
As with the absence seizure, we first investigate whether tonic–clonic seizures correspond to the expression of nonlinear structure in scalp EEG. Figure 8 shows a seizure (panel a) together with its nonlinear time series analysis (b). Prior to the seizure (t < 80 s) there is evidence of a nonlinear structure in the EEG only weakly at t 53 s. Immediately after the seizure onset (80–110 s), there exists evidence of nonlinearity in 3 of 4 consecutive epochs. There then follows a period (110–140 s) during which time, according to our nonlinear prediction technique, there is no evidence of nonlinear structure, despite appreciably large oscillations. Notably, the variance of the surrogate (null) distribution is particularly narrow during this period. The significance of this feature and the limitations of the nonlinear algorithm are discussed below. Strong nonlinear structure is then evident until the seizure terminates (155 s). There is no nonlinearity evident for the remainder of the recording.
|
A corresponding model simulation is given in panels (d)–(f): nonlinear structure is weakly evident in the first half of the time series as vse is held just below the bifurcation point (panel f). Strong nonlinear structure appears as soon as vse crosses the bifurcation point and remains evident until the seizure offset. Subsequently, there is a period of time during which nonlinear structure is strongly suppressed. It weakly reemerges once the alpha power is restored toward the end of the simulation (t > 45 s).
Hence, there is evidence for an increase in nonlinear structure in the tonic–clonic data, although it is much less striking than the absence seizure. Figure 8c shows the spectral density plot for the experimental data. As with the numerical seizure, increased power in the 10-Hz frequency peak coincides with the seizure onset (approximately 90 s). A peak persists throughout the seizure, although the peak frequency slows considerably. At one stage (100 s), 2 such peaks are visible. In comparison, the single peak frequency of the model seizure falls only marginally prior to the seizure offset (Fig. 7e). Of particular note is that, as with the model seizure, the experimental spectrum shows a marked postictal suppression of power across all frequencies (at 155 s). As discussed above, this is the critical fingerprint of traversing the bistable 10-Hz bifurcation diagram. Such a feature is also reflected in the highly asymmetrical character of both the modeled and experimental seizures.
Another critical feature of the epilepsy spectrum is the extremely strong increase in broad-spectrum (10–50 Hz) power from t 120 s until seizure termination. Given that this is a scalp EEG recording, such power putatively contains electromyelographic (EMG) contamination corresponding to the motor output of the seizure. It can be seen that this corresponds to the period of time during which the nonlinear prediction algorithm failed to find nonlinear structure in the data (panel b). It is probable that any putative nonlinear structure would be obscured by the high-amplitude noise, which also explains the small variance of the nonlinear index in the surrogate data. By comparison, no such EMG effect occurs in absence activity, so that the signal–noise ratio is more favorable.
The properties of the seizure are explored in further detail in Figure 9. Panel (a) shows the onset of the seizure. A rapid nonlinear growth of orbit amplitude is clearly visible and this proves to be close to exponential, consistent with dynamical evolution away from an unstable spiral outset as observed in the model. Comparing the seizure onset with the seizure termination in panel (b) reveals a number of features. Firstly, in contrast to the absence seizure, the waveform has changed markedly—suggesting nonstationarity of the dynamics. The frequency has fallen, and the waveform now has spike and slow waves. Secondly, the abrupt seizure offset is visible, and thirdly, the postictal EEG suppression can be seen. The strong nonstationarity of the dynamics represents an additional obstacle for the nonlinear prediction algorithms employed in Figure 8. However, the marked temporal asymmetry of each oscillation visible just prior to seizure termination is a classic feature of nonlinear dynamics (Stam and others 1998).
|
In the model seizure presented above, only the single parameter vse was varied. It is inevitable that, during a lengthy seizure, other physiological parameters may vary as a consequence of the abnormally high activity and/or hypoxia resulting from disruptions to normal respiratory effort. Other parameters may additionally be varied by regulatory mechanisms in order to terminate the seizure (Engel and others 1997
; McCormick and Contreras 2001). One parameter that would be likely to vary is the corticothalamic delay term t0 as this depends upon axonal and synaptic transit times that may be hindered by both sustained high activity (Poolos and others 1987) and hypoxia—leading, for example, to changes in ion concentration and diminished presynaptic neurotransmitter supply. The total transit time for a complete seizure also depends on the membrane time constants and ß—which would also be affected by these factors. For example, increased inhibitory influence of gamma-aminobutyric acid type B receptors (GABAB) would cause a decrease in and ß, increasing the round-trip time. However, we do not explicitly model these effects here. A numerical simulation is presented in Figure 10a–c. This was obtained by increasing t0 linearly from 80 ms to a "grossly pathological" 250 ms during the seizure, and then allowing it to return to its normal value after the seizure has terminated. The time series (a) and spectrum (b) reveal that the peak seizure frequency falls appreciably, as observed experimentally. The bistability fingerprint remains evident in the postictal spectral suppression. Introducing a second-parameter nonstationarity has the additional effect of increasing the nonstationarity of the seizure waveform, so that a spike and slow-wave morphology is now also observed in the model seizure (panel c).
|
It must be emphasized that this is an exploratory illustration of a varying one additional parameter after seizure initiation to simulate the secondary effect of a seizure (sustained high-frequency neuronal activity in a hypoxic environment) on neuronal physiology. A similar approach may also permit modeling of the decreasing frequency observed in the empirical absence seizure. In fact, it has previously been shown that drifting frequencies and cycle-to-cycle variations seen in absence seizures may contribute markedly to their nonlinear signal characteristics (Schiff, Victor, and Canel 1995
). Due to the intimate involvement of the intrathalamic loop in the absence waveform, changing membrane constants in this loop in addition to the corticothalamic loop could additionally contribute to a slowing peak frequency. An attempt to capture all the variance in the seizure dynamics through manipulation of a large number of system parameters (within their pathophysiological ranges) will be the subject of future research. Global Bifurcation Diagram
Figures 3a and 7a illustrate local bifurcation diagrams obtained by varying the single parameter vse. Each can be thought of as a 1-dimensional cross-section through the larger "global" bifurcation set spanned by all of the model's free parameters. It is natural to enquire as to the nature of such a set, which may shed light on a range of possible transitions to seizure activity in addition to being of mathematical interest in its own right.
The contrasting shape of the 2 local bifurcations suggests that they may be related through a simple "unfolding" of the bistability region of Figure 7a into the Hopf bifurcation of Figure 3a. That is, the region of bistability would be expected to grow smaller until the fixed point became unstable to the left (rather than the right) of the onset of the large-amplitude oscillations of Figure 7a. This hypothesis is motivated by the observation that such an unfolding is one of the basic geometric forms of all global bifurcation sets (Thom 1975). In order to test this hypothesis we generated a series of bifurcation diagrams stretching between those of Figures 3a and 7a in parameter space. This was achieved by linear interpolation of all parameters that vary between these 2 figures, as given in columns 2 and 3 of Table 1. Pertinent intermediate plots are given in Figure 11. In panel (a) is the tonic–clonic bifurcation of Figure 7 and in panel (e) is the "absence" bifurcation of Figure 3. The critical transition between the 2 types of bifurcations occurs at approximately 60–65% of the distance (in parameter space) between these 2 panels. Panel (b) shows the interpolation at 61% of the distance. Surprisingly, it can be seen that the fixed point has undergone a Hopf-type bifurcation within the zone of bistability. To the left of this instability, there hence exists both fixed-point and high-amplitude 10-Hz (aperiodic) oscillations. To the right of this instability, the same 10-Hz oscillations coexist with smaller amplitude 3-Hz periodic oscillations. Panel (c) shows interpolation at 63% of the distance. It can be seen that the Hopf bifurcation to 3-Hz activity occurs to the left of the region of bistability. Hence, 10-Hz periodic and 3-Hz periodic oscillations coexist within this region. By 64% of the distance (panel d), the Hopf and subsequent period-doubling bifurcations are now the only apparent nonlinear instabilities. The 10-Hz periodic branch of the tonic–clonic seizure cannot be accessed in this region of parameter space. Such a period-doubling–type bifurcation set continues until the absence case (panel e) is reached.
|
Hence, rather than there being a smooth transition between the 2 bifurcations, we instead see that they intersect orthogonally. The tent diagram (Fig. 2) provides an alternative means of visualizing this observation, in that the surfaces corresponding to the alpha and theta instabilities collide in the xyz parameter space outside the tent. This implies that both seizure types may coexist in the same region of parameter space and reflects upon the complexity of the high-dimensional dynamics of brain systems. This is consistent with clinical findings that patients with absence seizures are more likely than average to additionally have tonic–clonic seizures and vice versa (Stefan and Snead 1997
; Zifkin and Dravet 1997).
The global bifurcation diagram also provides some insight into the stability of the various epilepsy waveforms. The cross-sections given in Figure 11 and additional ones (not presented) reveal that there are large, open region sets of parameter space devoid of bifurcations. Changing parameters in these regions will merely result in smooth deformations of the amplitudes and slopes of the seizure oscillations. Only changes in the vicinity of a bifurcation—which are exceptional—will lead to a new waveform. This observation argues that the present approach is sufficiently robust to account for individual variability and the limitations of estimating exact underlying parameter values. However, in some regions of parameter space, such as the right half of Figure 11d, the system undergoes numerous bifurcations with only modest—or at times infinitesimal—alterations in one or more parameters. In these regions, we observe multiple periodic and chaotic spike and polyspike waveforms. The existence of these nonlinearly unstable regions is consistent with studies of "generic" nonlinear systems (Barreto and others 1997) and as noted above, could explain the variety of spike and wave morphologies observed clinically.
|
Conclusion |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
By performing the bifurcation analysis of a model of large-scale brain activity, this paper presents a parsimonious and unifying explanation of the defining features of the 2 major human generalized seizures—and their main differences. In summary, we used a Hopf bifurcation commencing the transition from healthy resting EEG to absence (3Hz) seizures. This yielded time series realizations with periodic spike and wave morphology with close similarity to scalp EEG data taken from an absence seizure data-base. Moreover, the nature of the bifurcation set yields a symmetrical on–off character that is also found in the EEG data. In contrast, the bifurcation diagram for tonic–clonic (10 Hz) seizures shows a sudden, discontinuous transition from damped (fixed-point) dynamics to large-amplitude nonlinear oscillations. The bifurcation set was of a "subcritical" nature—that is, associated with bistable behavior. This yielded time series realizations compatible with tonic–clonic seizures. Critically, once the seizure commenced, the system displays high-amplitude periodic oscillations and is required to traverse a considerable distance through parameter space before a stable (damped) linear regime reemerges. This asymmetry provides the explanation for the differences between pre- and postictal EEG spectra and cognition that are the defining features of tonic–clonic seizures. This is the central finding of this study. It extends upon and tests the predictions concerning the onset of epileptic activity in Robinson and others (2002)
that seizure phenomena arise when corticothalamic dynamics lose linear stability in specific regions of parameter space. We are not aware of any prior unifying explanations for these key features of the generalized seizures.
A common criticism of many physiologically based models of neural activity is that their numerous parameters undermine their explanatory power. Although this is an important consideration, a number of features of the current study are relevant in this regard. First, the model we employed was not formulated specifically to generate seizure waveforms but rather to incorporate the critical features of corticothalamic dynamics in order to explain the EEG temporal and spatial spectra in healthy resting and sleep states—when the activity is weakly damped and the nonlinearities can be neglected (Robinson and others 1997, 2001, 2002; O'Connor and others 2002). In the model, as in the clinical setting, seizures arise from "background" resting EEG states when—due to a change in an underlying physiological property—the nonlinearity surpasses a critical threshold. The behavior of the model's dynamical variables during nonlinear dynamics permits predictions to be made regarding real physiological processes. In the present study, such prior predictions (Robinson and others 2002) were further elaborated and compared to physiological data. Second, we observed that setting the model within the "resting state" regime but in the vicinity of a bifurcation point yielded weakly nonlinear alpha oscillations consistent with previous nonlinear analyses of resting EEG data (Stam and others 1999; Breakspear and Terry 2002). Third, there are no free parameters in the sense that they all correspond to discrete physiological processes that have been constrained by independent empirical estimates and matched against a large database of EEG spectra (Robinson and others 2004). We used previously published parameters and did not perform any a priori exploration of parameter space. Fourth, our numeric simulations yielded a number of phenomena that lent themselves to potential refutation, such as the periodic nature of the absence seizures and the spectral properties of pre- and postictal EEG. Finally, we employed a nonlinear prediction algorithm in order to verify that the onset of generalized seizures corresponds to a bifurcation from damped to strongly nonlinear behavior. The present study thus represents a predictive application of an existing model in a novel direction, rather than an exploratory or confirmatory study.
Several interesting questions remain to be answered. First, does the onset of the spike and slow-wave morphology represent the smooth deformation of a limit-cycle attractor or the superposition of 2 cycles? The present study suggests the latter possibility, although proof will require explicit analysis of the nonlinear stability properties of the limit cycle. Second, the corticothalamic model employed here does not explicitly include T-channels within the thalamus. Although the present model is not inconsistent with such channels, their explicit incorporation would be required to model the sensitivity of absence seizures to T-channel blockade with ethosuximide. T-channels, which show a refractory period following sustained bursting, have been incorporated into related neural population models of 3-Hz absence seizures (Destexhe and Sejnowski 2001). Destexhe and Sejnowski treated axonal propagation times as negligible, and hence, there was no time delay. In comparison, the presence of time-delayed corticothalamic feedback is a crucial ingredient in EEG spectra and the absence waveform we observed. In fact, qualitatively similar spike–wave oscillations can be generated by our model with a variety of time-delayed feedback loops that differ in some detail from those presented above. This is interesting because it is known experimentally that absence seizures can arise as a result of changes in a number of different neuronal pathways. However, the present model incorporates the important components of the corticothalamic system together with the time delays that—due to finite axonal conduction speeds—are present in the brain (Meeren and others 2002). Interestingly, although there are some differences in comparison to the Destexhe and Sejnowski model, both models emphasize the increased excitatory loops between the cortex and the specific and reticular nuclei of the thalamus underlying the generalized seizures. Third, recent analysis reveals that spikes can be generated without any time-delayed loops at all. These may provide an explanation for purely cortically generated spikes (McCormick and Contreras 2001).
For this study, we investigated a larger number of clinical seizures than were presented in the paper. However, the properties of these seizures that are pertinent to the validity of the model—such as the periodic waveform of the absence seizure—are the common and defining features of these seizure types and not at all limited to the examples presented. Furthermore, some features of the corticothalamic model seizures, such as the higher harmonics and subharmonics of the tonic–clonic seizure (Figs 7e and 10b), were not apparent in our EEG examples. However, these are quite obvious in intracranial seziure recordings (Schiff and others 2000), suggesting that they are simply obscured by the EMG activity present in scalp recordings. The second low-frequency peak seen in the empirical data (Fig. 8c) at t 100 s does not permit such a straightforward explanation. Possible explanations include 2 sequential crossings of the instability boundary at slightly different regions due to rapidly changing parameter values or the secondary excitation of a mode not captured by the present model. Further modeling and empirical study is required to establish the significance of this second peak and, if robust, seek a valid explanation for its origin.
It also important to note that the present paper employed primarily numerical analysis of the proposed corticothalamic model. Numerical analysis of delay-differential equations has its own caveats, such as the choice of the time-delayed values of the numerical interpolations. We chose a third-degree (cubic) interpolate and observed almost identical results for a variety of time steps, implying convergence. Increasing the time step or decreasing the accuracy of the time-delayed interpolates has the effect, as expected, of decreasing the values of vse at which bifurcations occur. This is consistent with the presence of multiplicative numerical noise. We are currently undertaking a formal stability (Lyapunov) analysis in order to complement the results presented here.
To illustrate the effect of hypoxia and neural fatigue on the modeled tonic–clonic seizure, we linearly increased the corticothalamic conduction delay (Poolos and others 1987). This had the effect of slowing the peak seizure frequency, imparting the "chirp" property observed in the clinical recordings (Schiff and others 2000). Although it is unlikely that the corticothalamic delay would increase as markedly (from 80 to 250 ms) as in the exemplar seizure we present, it is probable that other physiological processes would also be affected by the impact of a tonic–clonic seizure. For example, active inhibitory mechanisms are likely to involve an increase in the number of activated GABA receptors, which in the context of our model corresponds to an increase in synaptodendritic rise and decay times. Rather than attempting to study the influence of many changing parameters simultaneously, we choose to vary only a single likely target of hypoxia and neural fatigue—corticothalamic conduction time. Future work could study the manipulation of other model parameters. Such work, together with the incorporation of other physiological mechanisms into the model (such as T-channels) may have the effect of improving the match between the model seizure waveform and the properties of the clinical EEG. Mismatches between a model and observed phenomena are, in fact, to be expected whenever relevant physiological mechanisms are omitted, as they indicate precisely where the model requires refinement. Simply adding more and more detail to a model, however, without testing simpler realizations, adds to the validation problem discussed above. Critically, the present paper illustrates that a relatively simplified corticothalamic model is able to capture and explain the key features of the generalized seizures in a novel, unified manner.
|
Appendix |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
In this appendix we first give the 8 first-order delay-differential equations describing the brain in the spatially uniform case, followed by a brief description of the stability variables x, y, and z.
First-Order Delay-Differential Equations for the Spatially Uniform Model
We assume that intracortical connectivities are proportional to the numbers of synapses involved, implying Vi = Ve and Qi = Qe (Robinson and others 2002), so that the cortical inhibitory and excitatory processes are mutually enslaved. The smallness of ri and the thalamic nuclei allow us to set c , c = i, r, s, yielding the local approximation from equations (2) and (7), c(t) = S[Vc(t)].
Thus, from equations (2) and (4–7), we obtain
 | (A1) |
 | (A2) |
 | (A3) |
 | (A4) |
 | (A5) |
 | (A6) |
 | (A7) |
 | (A8) |
Reduced Parameter Space
Linear stability analysis suggests a 3-dimensional parameterization of low-frequency instabilities, defining an xyz coordinate space by (Robinson and others 2002)
 | (A9) |
 | (A10) |
 | (A11) |
where gain Gab = avab is the response in neurons a to unit input from neurons b, sigmoid slope a = dS(Va)/dVa, with Gese = GesGse, Gesre = GesGsrGre, and Gsrs = GsrGrs for convenience. In previous work, these parameters have only been used in the steady state where they were derived. Here we take them to define a coordinate transformation of the dynamical variables Va(t). The parameters x, y, and z give a measure of the level of activity in the brain structures comprising the model. Specifically, x describes purely cortical activity, y describes activity in the 2 corticothalamic loops (i.e., via the specific nuclei and via both the reticular and specific nuclei), and z describes activity in the intrathalamic loop (i.e., between the reticular and specific nuclei). Thus, x, y, and z relate to cortical, corticothalamic, and intrathalamic stabilities, respectively.
|
Acknowledgments |
|---|
The authors would like to thank A. Bleasel for assisting with absence data, R. Quian Quiroga for providing the tonic–clonic data, and 2 anonymous reviewers for helpful comments on an earlier version of the manuscript. MB was a recipient of an American Psychiatric Association/AstraZeneca Research Award. JRT would like to acknowledge funding from the Nuffield Foundation via a newly appointed lecturer grant and the Leverhulme Trust for funding the Theoretical Neuroscience Network. SR would like to acknowledge the financial support of the Loughborough University Sleep Research Centre and to the Bristol Laboratory for Applied Dynamical Engineering. The authors also acknowledge the support of an Australian Research Council grant and a University of Sydney SESQUI Research and Development grant.
|
References |
|---|
|
TopAbstractIntroductionMethodsResultsConclusionAppendixReferences |
|---|
Abraham RH, Shaw CD. 1992. Dynamics—the geometry of behavior. Redwood City, CA: Addison-Wesley.
Andrezejak RG, Laehnertz K, Mormann F, Rieke C, David P, Elger CE. 2001. Indications of nonlinear deterministic and finite-dimensional structure in time series of brain electrical activity: dependence on recording region and brain state. Phys Rev E 64:061907.
Arnhold J, Grassberger P, Lehnertz K, Elger C. 1999. A robust method for detecting interdependencies: application to intracranially recorded EEG. Physica D 134:419–430.[CrossRef]
Barreto E, Hunt B, Greborgi C, Yorke J. 1997. From high dimensional chaos to stable periodic orbits: the structure of parameter space. Phys Rev Lett 78:4561–4564.[CrossRef]
Breakspear M. 2002. Nonlinear phase desynchronization in human electroencephalographic data. Hum Brain Mapp 15:175–198.[CrossRef][Web of Science][Medline]
Breakspear M, Terry JR. 2002. Detection and description of nonlinear interdependence in normal multichannel human EEG. Clin Neurophysiol 113:735–753.[CrossRef][Web of Science][Medline]
Breakspear M, Terry JR, Friston K. 2003. Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a nonlinear model of neuronal dynamics. Network Comput Neural Syst 14:703–732.[CrossRef]
Breakspear M, Williams L, Stam CJ. 2004. Topographic analysis of phase dynamics in neural systems reveals formation and dissolution of ‘dynamic cell assemblies’. J Comput Neurosci 16:49–68.[CrossRef][Web of Science][Medline]
Casdagli M. 1989. Nonlinear prediction of chaotic time series. Physica D 35:335–357.[CrossRef][Web of Science]
Crevier DW, Meister M. 1998. Synchronous period-doubling in flicker vision of salamander and man. J Neurophysiol 79:1869–1878.
Daffertshofer A, Peper CE, Frank TD, Beek PJ. 2000. Spectral analysis of event-related encephalographic signals. Hum Mov Sci 19:475–498.[CrossRef][Web of Science]
Dammig M, Mitschke F. 1993. Estimation of Lyapunov exponents from time series: the stochastic case. Phys Lett A 178:385–394.[CrossRef]
Destexhe A, Sejnowski TJ. 2001. Thalamocortical assemblies. Oxford: Oxford University Press.
Engel J, Dichter MA, Schwartzkroin PA. 1997. Basic mechanisms of human epilepsy. In: Engel J, Pedley TA, editors. Epilepsy: a comprehensive textbook. Philadelphia, PA: Lippincott-Raven.
Feucht M, Moller U, Witte H, Schmidt K, Arnold M, Benninger F, Steinberger K, Friedrich MH. 1998. Nonlinear dynamics of 3Hz spike-and-wave discharges recorded during typical absence seizures in children. Cereb Cortex 8:524–533.
Freeman W. 1975. Mass action in the nervous system: examination of the neurophysiological basis of adaptive behaviour through the EEG. New York: Academic Press.
Freeman WJ, Rogers LJ. 2002. Fine temporal resolution of analytic phase reveals episodic synchronization by state transitions in gamma EEGs. J Neurophysiol 87:937–945.
Friston KJ. 2000. The labile brain. I. Neuronal transients and nonlinear coupling. Philos Trans R Soc Lond B 355:215–236.
Friston KJ, Jezzard PJ, Turner R. 1994. Analysis of functional MRI time-series. Hum Brain Mapp 1:153–171.[CrossRef]
Fuchs A, Mayville JM, Cheyne D, Weinberg H, Deecke L, Kelso JAS. 2000. Spatiotemporal analysis of neuromagnetic events underlying the emeregence of coordinative instabilities. Neuroimage 12:71–84.[CrossRef][Web of Science][Medline]
Gratton G, Coles M, Donchin E. 1983. A new method for off-line removal of ocular artifact. Electroencephalogr Clin Neurophysiol 55:468–484.[CrossRef][Web of Science][Medline]
Izhikevich E. 2001. Synchronization of elliptic bursters. SIAM Rev 43:315–344.[CrossRef]
Jirsa VK, Haken H. 1996. Field theory of electromagnetic brain activity. Phys Rev Lett 77:960–963.
Kelso J, Bressler S, Buchanan S, DeGuzman G, Ding M, Fuchs A, Holroyd T. 1992. A phase transition in human brain and behaviour. Phys Lett A 169:134–144.[CrossRef]
Lopes da Silva F, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. 2003. Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia 44(Suppl 12):72–83.
McCormick DA, Contreras D. 2001. On the cellular and network bases of epileptic seizures. Annu Rev Physiol 63:815–846.[CrossRef][Web of Science][Medline]
Meeren HK, Pijn JP, Van Luijtelaar EL, Coenen AM, Lopes da Silva FH. 2002. Cortical focus drives widespread corticothalamic networks during spontaneous absence seizures in rats. J Neurosci 22:1480–1495.
Nunez PL. 1974. The brain wave equation: a model for the EEG. Math Biosci 21:279–297.[CrossRef]
Nunez PL. 1995. Neocortical dynamics and human EEG rhythms. Oxford: Oxford University Press.
O'Connor SC, Robinson PA, Chiang AKI. 2002. Wavenumber spectra of EEG signals. Phys Rev E 66:061905.[CrossRef]
Perez Valezquez JL, Cortez MA, Carter Snead O III, Wennberg R. 2003. Dynamical regimes underlying epileptiform events: role of instabilities and bifurcations in brain activity. Physica D 186:205–220.[CrossRef]
Poolos NP, Mauk MD, Kocsis JD. 1987. Activity-evoked increases in extracellular potassium modulates presynaptic excitability in the CA1 region of the hippocampus. J Neurophysiol 58:404–416.
Prichard D, Theiler J. 1994. Generating surrogate data for time series with several simultaneously measured variables. Phys Rev Lett 73:951–954.[CrossRef][Web of Science][Medline]
Quian Quiroga R, Garcia H, Rabinowich A. 2002. Frequency evolution during tonic-clonic seizures. EMG Clin Neurophysiol 42:323–331.
Robinson PA, Rennie CJ, Rowe DL. 2002. Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys Rev E 65:041924.[CrossRef]
Robinson PA, Rennie CJ, Rowe DL, O'Connor SC. 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum Brain Mapp 23:53–72.[CrossRef][Web of Science][Medline]
Robinson PA, Rennie CJ, Wright JJ. 1997. Propagation and stability of waves of electrical activity in the cerebral cortex. Phys Rev E 56:826–840.[CrossRef]
Robinson PA, Rennie CJ, Wright JJ, Bahramali H, Gordon E, Rowe DL. 2001. Prediction of electroencephalographic spectra from neurophysiology. Phys Rev E 63:021903.[CrossRef]
Rombouts S, Keunen R, Stam C. 1995. Investigation of nonlinear structure in multichannel EEG. Phys Lett A 202:352–358.[CrossRef]
Schiff ND, Victor JD, Canel A. 1995. Nonlinear autoregressive analysis of the 3/s ictal electroencephalogram: implications for underlying dynamics. Biol Cybern 72:527–532.[Web of Science][Medline]
Schiff ND, Victor JD, Canel A, Labar DR. 1995. Characteristic nonlinearities of the 3/s ictal electroencephalogram identified by nonlinear autoregressive analysis. Biol Cybern 72:519–526.[Web of Science][Medline]
Schiff S, Colella D, Jacynac GM, Hughes E, Creekmore JW, Marshall A, Bozek-Kuzmicki M, Benkec G, Gaillard WD, Conry J, Weinstein SR. 2000. Brain chirps: spectrographic signatures of epileptic seizures. Clin Neurophysiol 111:953–958.[CrossRef][Web of Science][Medline]
Schiff S, So P, Chang T, Burke R, Sauer T. 1997. Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble. Phys Rev E 54:6708–6724.[CrossRef]
Stam CJ, Pijn J, Pritchard W. 1998. Reliable detection of nonlinearity in experimental time series with strong periodic components. Physica D 112:361–380.[CrossRef][Web of Science]
Stam CJ, Pijn J, Suffczynski P, Lopes da Silva FH. 1999. Dynamics of the alpha rhythm: evidence for non-linearity? Clin Neurophysiol 110:1801–1813.[CrossRef][Web of Science][Medline]
Stam CJ, Pritchard W. 1997. Dynamics underlying rhythmic and non-rhythmic variants of abnormal, waking delta activity. Int J Psychophysiol 34:5–20.
Stam CJ, van Dijk BW. 2002. Synchronization likelihood: an unbiased estimate of generalized synchronization in multivariate data sets. Physica D 163:236–251.[CrossRef]
Stam CJ, van Woerkom T, Keunen R. 1997. Nonlinear analysis of the EEG in Creutzfeldt-Jakob disease. Biol Cybern 77:247–256.[CrossRef][Web of Science][Medline]
Stefan H, Snead OC. 1997. Absence seizures. In: Engel J, Pedley TA, editors. Epilepsy: a comprehensive textbook. Philadelphia, PA: Lippincott-Raven.
Terry JR, Breakspear M. 2003. An improved algorithm for the detection of nonlinear interdependencies Biol Cybern 88:129–136.[CrossRef][Web of Science][Medline]
Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer J. 1992. Testing for nonlinearity: the method of surrogate data. Physica D 58:77–94.[CrossRef][Web of Science]
Thom R. 1975. Structural stability and morphogensis. Reading, MA: W.A. Benjamin.
Wendling F, Bartolomei F, Bellanger JJ, Chauvel P. 2002. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur J Neurosci 15:1499–1508.[CrossRef][Web of Science][Medline]
Wright JJ, Kydd R, Lees G. 1985. State changes in the brain viewed as linear steady-states and nonlinear transitions between steady-states. Biol Cybern 53:11–17.[CrossRef][Web of Science][Medline]
Zifkin B, Dravet C. 1997. Generalized convulsive seizures. In: Engel J, Pedley TA, editors. Epilepsy: a comprehensive textbook. Philadelphia: Lippincott-Raven.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  F. Freyer, K. Aquino, P. A. Robinson, P. Ritter, and M. Breakspear Bistability and Non-Gaussian Fluctuations in Spontaneous Cortical Activity J. Neurosci., July 1, 2009; 29(26): 8512 - 8524. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 V. K Jirsa Neural field dynamics with local and global connectivity and time delay Phil Trans R Soc A, March 28, 2009; 367(1891): 1131 - 1143. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
F. Marten, S. Rodrigues, O. Benjamin, M. P Richardson, and J. R Terry Onset of polyspike complexes in a mean-field model of human electroencephalography and its application to absence epilepsy Phil Trans R Soc A, March 28, 2009; 367(1891): 1145 - 1161. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 L. J. Voss, J. W. Sleigh, J. P. M. Barnard, and H. E. Kirsch The Howling Cortex: Seizures and General Anesthetic Drugs Anesth. Analg., November 1, 2008; 107(5): 1689 - 1703. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||