In this paper, we introduce a modification of a mean-field model used to describe the brain's electrical activity as recorded via electroencephalography (EEG). The focus of the present study is to understand the mechanisms giving rise to the dynamics observed during absence epilepsy, one of the classical generalized syndromes. A systematic study of the data from a number of different subjects with absence epilepsy demonstrates a wide variety of dynamical phenomena in the recorded EEG. In addition to the classical spike and wave activity, there may be polyspike and wave, wave spike or even no discernible spike–wave onset during seizure events. The model we introduce is able to capture all of these different phenomena and we describe the bifurcations giving rise to these different types of seizure activity. We argue that such a model may provide a useful clinical tool for classifying different subclasses of absence epilepsy.
Epilepsy is the most common serious primary brain disease affecting 380 000 people in the UK, with 30 000 new cases per year. Epilepsy carries significant mortality (Cockerell et al. 1994) and morbidity (Buck et al. 1997) as well as reduced quality of life (Devinsky et al. 1995). A recent National Sentinel Audit revealed an annual mortality in the UK of more than 2000 cases per year. Epilepsy has high costs, estimated at £1.93 billion per year in the UK in 1994 (Cockerell et al. 1994). The focus of the present paper is on a class of primary generalized seizures, called absence seizures, which typically affect children and young adolescents. The term ‘absence’ arises from the clinical presentation of these seizures whereby cognitive abilities are lost, often resulting in brief periods of behavioural arrest (or the absences). When monitoring subjects with epilepsy using electroencephalography (EEG), dynamically evolving patterns may be observed in cortical electrical activity; the classical pattern associated with absence seizures is an approximately 3 Hz spike and wave discharge that is synchronous across all channels of activity. However, closer inspection of EEG recorded from different subjects (e.g. Hrachovy & Frost 2006) demonstrates a number of different dynamical phenomena, e.g. polyspike and wave discharges and wave–spike discharges, in addition to the classical spike–wave profile (figure 1). Previous theoretical work by Breakspear et al. (2006) and Rodrigues et al. (2006, 2007) focused on understanding the transitions from healthy states to the classical spike and wave discharge, demonstrating that the oscillatory behaviour (the wave part of the complex) arises as a result of a Hopf bifurcation (HB) upon increasing excitatory connections between the cortex and the thalamus. Further increasing this parameter gives rise to the spike due to the appearance of an inflection point in the vector field. In addition, these previous studies have demonstrated the predictive and descriptive validity of a mean-field approach for modelling a wide range of healthy states as well as generalized seizures in humans (Breakspear et al. 2006). The studies of EEG recordings collected during the absence seizures lend themselves naturally to this approach, since during an absence seizure there exists a hypersynchronous state that entrains the populations of neurons to the same dynamic activity.
In the present work, we introduce a mean-field model of cortico-thalamic interactions and examine the mechanisms by which transitions from healthy states to those observed from a number of different absence seizures can occur. Two model parameters are indicated to be crucial in this respect. One is the parameter that governs the strength of excitatory cortical–thalamic interactions, which has been demonstrated to be crucial in our previous studies. The other parameter is a time delay that we introduce to describe the different time scales of the interactions for GABAA and GABAB receptors in the specific relay nuclei. We demonstrate that this delay is crucial in determining the number of spikes per cycle during seizure activity in the model and that the excitatory coupling plays a role in the transition from healthy to seizure dynamics. It is natural to consider such interactions as several human and animal studies of absence seizures have implicated abnormalities in cortico-thalamic loops in the generation of seizure activity (Crunelli & Leresche 2002).
Our work is motivated by a desire to elucidate potential mechanisms for the onset of seizure activity, which could in turn inform novel experimental and clinical approaches to studying epilepsy, e.g. the potential to develop seizure prediction algorithms based on the bifurcation tracking of parameters of interest that could provide an alternative approach to the standard one based on signal processing (see Mormann et al. (2005) for a review).
2. Clinical data
The data presented in this paper came from a database of 48 seizures from 20 subjects, who were part of a consecutive series having absence seizures during routine outpatient EEG investigation. The EEG data were obtained using conventional clinical equipment. Standard silver–silver chloride disc electrodes attached to the scalp were used, and electrode placement followed the Maudsley variant of the international 10–20 system. The EEG recordings typically lasted for approximately 1 hour, with periods of wakefulness with eyes open and closed, as well as drowsiness and sleep in some patients. Attempts were made in all cases to provoke the seizures using hyperventilation and photic stimulation (stroboscopic stimulation at a range of frequencies). All the EEG data were reviewed by a trained clinical encephalography technician, who identified spike–wave discharges for further analysis here.
In addition to the EEG data, medical records were reviewed to establish the epilepsy syndrome according to the classification of the Commission on Classification and Terminology of the International League Against Epilepsy (1989). In each case, the patients were assigned to one of the four syndromes: childhood absence epilepsy (CAE); juvenile absence epilepsy (JAE); juvenile myoclonic epilepsy; and generalized tonic–clonic seizures on awakening. While these syndromes are clearly defined in terms of the age of onset, other seizure types experienced by the patient and other clinical criteria, there is a considerable overlap between these syndromes. Although absence seizures are seen more often in children than in adults, persistence of absence epilepsy from childhood into adulthood is well recognized, as is de novo adult onset (Trinka 2005). The patients studied here had an age range of 3–65 years.
A typical absence is a non-convulsive epileptic seizure. It consists of a sudden, brief impairment of consciousness accompanied by a generalized, synchronous, bilateral, 2.54 Hz spike and slow wave discharge seen on the EEG, and followed by abrupt recovery (Crunelli & Leresche 2002). Typical absence seizures occur in several different idiopathic generalized epilepsy syndromes, and detailed phenomena during the absence seizure may vary between syndromes and between patients, for example in duration and accompanying motor features such as myoclonic jerks of limbs or eyelids. The ictal EEG usually shows a single spike followed by a wave, but may show polyspikes. The discharges are traditionally regarded as synchronous across all the channels and usually bilaterally symmetrical. Typically, the spike–wave discharge frequency is highest at the start of the burst and falls through the seizure. As well as inter-individual variability in the EEG and in the clinical presentation, there is considerable variability in response to the drug treatment (Duron et al. 2005).
3. Mathematical description
(a) The mean-field model
In order to model the EEG activity in patients with absence seizures, we investigate the neural dynamics in two distinct brain regions, namely the thalamus and the cerebral cortex. The motivation for this approach is that the thalamo-cortical interactions are found to play an essential role in the generation of this (poly)spike–wave EEG activity. A variety of experimental evidence is summarized in §8.1 of Destexhe & Sejnowski (2001). The model that we present may be thought of as a mean-field, neural mass model. It describes the dynamics of large interacting groups of neurons, which are referred to as the neural masses. The model arises as a result of merging a number of theoretical viewpoints.
The differences in the time scales of activation/inactivation due to the GABAA/B receptors in the thalamus as described in figs 5.8 and 5.9 in the book by Destexhe & Sejnowski (2001).
Incorporation of the cortico-thalamic loop, shown in studies (Destexhe & Sejnowski 2001) to be important in the generation of sleep spindles and abnormal seizure activity.
The wave-like equation for the propagation of cortical activity, as described by Jirsa & Haken (1996).
Merging these approaches together, we consider an extension of an existing mean-field approach, employed in previous studies (Breakspear et al. 2006; Rodrigues et al. 2006), consisting of both cortical and thalamic components. The thalamic component is assumed to consist of two neural masses; an excitatory mass of specific relay nuclei (s) and an inhibitory mass of reticular nuclei (r). Similarly, the cortical component incorporates both a mass of excitatory pyramidal cells (e) and inhibitory interneurons (i).
Each of the neural masses (e, i, s, r) is described by three dynamical variables: its average membrane potential Va(r, t), the average firing rate ςa(r, t) and the axonal field ϕa(r, t) at position r and time t. The subscript (a) is used to refer to each of the different masses. These variables obey the following dynamical rules:(3.1)(3.2)and(3.3)
Equation (3.1) describes the changes in the average membrane potential Va of a neural mass, under incoming postsynaptic potentials Pa. Equation (3.2) describes the average firing rate of the mass. If Va exceeds a threshold θa, the neural mass fires action potentials with an average rate ςa. The sigmoid shape of ςa ensures that the firing rate never exceeds a maximum value . This is crucial, because the real neurons cannot fire infinitely quickly. Equation (3.3) describes the propagation of action potential fields ϕa via axons of the neural mass. The axonal fields influence other neural masses via the synaptic connections, hence they can lead to new postsynaptic potentials P(r, t).
In figure 2, we present a schematic of all the connections within the mean-field model. We view the cortex as one compartment, with axonal connections to and from the two thalamic subcompartments. The excitatory connections via the fields ϕe, ϕs are displayed as arrows, while the inhibitory connections ϕr are given by lines with dots. The specific relay nuclei also receive the sensory input ϕn, which is modelled as a constant field in our research. Incoming synaptic fields ϕa, received by each compartment, are assumed to be linearly summed in the dendritic tree (Robinson et al. 2002), hence(3.4)for neural mass (a), where νab are coupling constants for the average strength of synaptic connections, which closes the set of rules of the mean-field model.
A further change from the previous studies is the nature of the delay. In the present work, we are motivated by the studies of Destexhe & Sejnowski (2001) and introduce two connections ϕr(t) and ϕr(t−τ) from the reticular nuclei to the specific relay nuclei. These connections are designed to represent the inhibitory GABAA and GABAB synapses. Since the GABAB functions via second messenger processes, its postsynaptic currents develop on a slower time scale than the currents of GABAA (figs 5.8 and 5.9 in Destexhe & Sejnowski (2001)). Hence, we expect that, if thalamic neurons fire, the GABAB-mediated connection to the specific relay nuclei will be delayed with respect to GABAA. In research by Coombes et al. (2002), a similar delay was obtained as follows: the wave-like equation (3.3) is derived from an integral between the firing rate ςa(r, t) and a synaptic kernel. Depending on the particular synapse modelled, this integral can lead to a delay in the right-hand side of equation (3.3). Combining this with equation (3.4), and the computations in appendix A, gives a delayed connection ϕr(t−τ).
In summary, equations (3.1)–(3.4) are the theoretical framework of our research. Since such a large system of DDEs is computationally expensive, it is desirable to make some simplifications. First of all, the cortical inhibitory interneurons are relatively sparse compared with pyramidal cells (a ratio of 4 : 1 in favour of the pyramidal cells). The excitatory (e) neurons are believed to be the main source of the EEG, hence we approximate Vi≈Ve and ςi≈ςe (Robinson et al. 2002). Furthermore, aa represents a ratio between an axonal velocity, and a typical length scale for axonal propagations. Since the inhibitory neurons only have short-range projections, compared with the excitatory neurons (Robinson et al. 1997, 2002), we effectively have γi,r→∞. A similar argument, explained in Rodrigues et al. (2008), can be used to approximate γs→∞.
In Robinson et al. (2002), it is shown that removing the time derivative in the right-hand side of (3.3) does not influence the numerical results; hence, it can be dropped for all the neural masses. If we combine this with γi,r,s→∞, equation (3.3) reduces to ϕi,r,s=ςi,r,s. The excitatory cortical neurons (e) cannot be approximated this way. However, because we model a generalized seizure, with global cortical activity, we neglect the spatial dependence of fields and voltages in our model. The above approximations reduce our model to a system of eight first-order delay differential equations, given in appendix A. By making this reduction, our model becomes computationally more tractable.
(b) Numerical methods
Large systems of the delay differential equations, similar to our model (A 2), can produce a variety of complex dynamics. In our case, we are especially interested in 2–4 Hz periodic solutions, which resemble the EEG activity during the absence seizures. Firstly, numerical integrators can be used to obtain the time series of the model (A 2). The results can be qualitatively compared with EEG traces of the patients. The numerical integration is performed through a fourth-order Adams–Bashforth scheme, adjusted to model the delay differential equations (Baker et al. 1994).
Secondly, since we wish to understand the onset, change and termination of the absence seizure dynamics, and in particular the onset and termination of the (poly)spike patterns, we employ the Matlab package DDE-BIFTOOL (Engelborghs et al. 2001) to find and continue the branches of these solutions, as the parameters of the model are varied. A further package PDDE-CONT (Szalai 2005) is also used in instances where the continuing branches of periodic solutions becomes intractable using Matlab.
Our results are presented in §4a,b. In §4a, we demonstrate that the mean-field model can produce periodic patterns that closely resemble the EEG data recorded during human absence seizures. We considered a database of 48 seizures and observed that, as well as the classical single spike and wave per cycle, the dynamics in many seizures consisted of multiple (poly)spikes per cycle. We illustrate that, in the model, the formation of spikes occurs as a result of inflection points in the solutions of our mean-field equations. Moreover, the existence of these inflection points can be studied using one-parameter bifurcation diagrams.
We extend these results in §4b, using the numerical continuation packages described previously to track the bifurcations and the birth of inflection points in our model. This generalizes our previous analysis by providing a detailed understanding of the mechanisms leading to the onset and deformation of the (poly)spike patterns. We proceed to link these results back to the time-series data simulated from the model and draw a comparison between the model and the clinical data from our absence seizure database.
(a) EEG and model: time series and bifurcations
Previous research relating the mean-field models and the absence seizures has focused on reproducing single spike and wave patterns (Breakspear et al. 2006; Rodrigues et al. 2006), an example of which is given in figure 1a. These are the waveforms that are classically associated with the absence seizures. More general polyspike patterns, such as the ones displayed in figure 1b,c, have been studied in only limited cases, e.g. in Wendling et al. (2002). Our aim is to build on this research, and explain the mechanisms that cause the onset of such polyspike activity, which we observe frequently in our absence seizure database.
The approach we follow extends the earlier research of Rodrigues et al. (2006), where it was shown that the coupling strength νse, between pyramidal cells in the cortex and specific cells in the thalamus, is a crucial parameter for observing the onset of spike–wave activity. Essentially, this parameter is related to the weight of the synaptic interactions between these two neuronal populations, which could conceivably vary with the changes in underlying neurophysiology and as such is a sensible parameter to vary. The main variable of interest is ϕe(t). We recall that this variable describes the axonal fields of the cortical pyramidal cells, which give the main contribution to the EEG signals. Hence, we assume a functional relationship between ϕe and the scalp EEG voltage, which to the first approximation may be considered linear. However, we stress that, even in this simple case, the constant of proportionality is patient specific.
The bifurcation diagrams provide a systematic way to study the dynamics of ϕe(t). Here, the aim is to find qualitative changes in ϕe under variations of the parameter νse, e.g. the transitions from ‘healthy’ steady-state neuronal activity to ‘seizure’ (poly)spike–wave patterns. We construct the diagrams by integrating the system of delay equations (A 2) as described in §3. During the integration, all model parameters are kept fixed to the values in table 1. In Robinson et al. (2004), the estimates of the model parameters were obtained from a variety of sources, including the physiological estimates where available. We used these estimates as a guideline for our parameters, combined with numerical simulation to find the relevant regions of parameter space containing seizure-like activity. Once finished, we add a small increment δ to the value of νse. Equations (A 2) are then integrated again, with the previous solution as the initial condition. This procedure is repeated hundreds of times. During each step, we record the local maxima and minima of ϕe(t) and plot them against νse. For the GABAB coupling, we choose a delay of τ=100 ms. This represents the order of magnitude of rise and decay time of GABAB currents (Destexhe & Sejnowski 2001, figs 5.8 and 5.9). The bifurcation diagram is shown in figure 3a.
The dynamics we find are similar to past investigations (Rodrigues et al. 2006). For νse<1.48×10−3 V s, we observe that the system (A 2) resides in a steady state. Thus, all of its eight variables are constant in time, including ϕe(t). This can be compared with the clinical EEG data before the seizure occurs (see arrows from figure 3a,b). Under normal conditions, the EEG signals are noisy patterns around a steady-state average voltage. At νse≈1.48×10−3 V s, we find the onset of an approximate 3 Hz oscillation, through a HB. The emanating two curves in the diagram represent its local maxima and minima. Some patients develop similar oscillatory EEG signals at the onset of absence seizures (Rodrigues et al. 2006). Their patterns do not contain any spikes and just consist of approximately 3 Hz sine-like oscillations.
An interesting change is observed at νse≈1.66×10−3 V s. It seems like ϕe(t) undergoes a second bifurcation at this point. However, a closer look at the time series reveals that the solution developed a small local maximum, through an inflection point. An example is sketched in figure 3c. This maximum develops into the familiar spike, which characterizes the spike–wave patterns. This is by no means a period doubling (PD) bifurcation; the stability and the period of the solution do not change along νse≈1.66×10−3 V s. Rather, this point represents a local deformation of the periodic solution profile (Rodrigues et al. 2008). We term this as a ‘ghost bifurcation’, since it mimics the shape of a bifurcation in our diagram.
Ghost bifurcations have been observed in past research (Rodrigues et al. 2006, fig.1), where the classification of the solution was left undetermined. The diagrams in these investigations always displayed at most one single ghost bifurcation, which led to the classic spike–wave profiles. An interesting new feature of our present model is the occurrence of a second ghost bifurcation at νse≈1.8×10−3 V s. This corresponds to the appearance of a double spike–wave solution. In the real patient EEG, such patterns can also develop during an absence seizure (arrow to figure 3b). Based on these observations, we wish to address the following questions: how many spikes can appear in the solutions of our model? Is their existence related to the delayed GABAB coupling? And what happens to ϕe(t) as νse is increased further?
(b) Continuation in two parameters
The analysis presented in §4a has demonstrated that the coupling between the cortex and specific cells, νse, is a crucial parameter in this model. Increases in νse lead to the onset of seizure dynamics, as well as transitions from single spike and wave to the polyspike solutions. Recent numerical studies of the mean-field models of the general form presented in this paper have demonstrated that the delay τ also plays a crucial role in the observance of spikes. For example, in the absence of the delay, only a transition to oscillatory behaviour via a HB was observed, without the spikes arising (Rodrigues et al. 2008). Similarly, we do not observe the spikes arising in the time series produced from the model for values of τ<0.04 suggesting that the delay is again significant in the generation of spike and wave solutions. To investigate this further, we employ two-parameter continuation, where we study νse and the delay τ simultaneously. This method allows us to draw two-dimensional curves of the (ghost) bifurcations in a (νse, τ) parameter plane.
These diagrams, also known as the activity maps, are extremely important. The bifurcation curves are boundaries between the different dynamical regions of the model, such as the healthy regions, where the model exhibits the steady-state solutions, and the seizure regions, where the model produces the polyspike patterns. In previous research (Wendling et al. 2002), similar activity maps were obtained by simulating the EEG activity on a two-parameter grid. We extend this method using the numerical continuation tools discussed in §2. This technique is more useful for creating the activity maps, since finding the boundaries of a parameter region is more efficient than finding the entire region step by step and calculates the boundary to a higher resolution. We also wish to find curves of ghost bifurcations in the (νse, τ)-plane. The standard continuation methods are not appropriate in this case as the branch of solutions here does not correspond to a standard bifurcation. To overcome this, we employ a special detection method to find the inflection points (this approach was first used in Rodrigues et al. 2008).
Figure 4 shows the activity maps obtained. Figure 4a displays the various bifurcation curves, computed for our mean-field model (A 2). In figure 4b, we performed the same computation for the original cortico-thalamic model, which has been discussed in the works of Breakspear et al. (2006) and Rodrigues et al. (2006, 2008). Exemplars of the solution profiles produced by each model demonstrate the need to consider modifications of the model. In the original model, although multiple inflection points exist, they do not develop into spikes as is the case for the model studied in this paper (which we return to in due course). The main difference between the two approaches is the incorporation of the delay τ. Whereas we propose in the present study a delay to describe the slow GABAB synapse, the model in figure 4b assumed that the delay arises owing to synaptic transmission times for the signals travelling between the cortex and the thalamus.
In figure 4a,b, we find a single curve of the HB and various ghost bifurcations, given by solid black curves. We also find the PD and the saddle-node bifurcations of limit cycles (SL), given by dashed and solid grey curves, respectively. Their significance will be discussed below. The different regions in figure 4, bounded by the (ghost) bifurcations, are marked by various shades of grey. The Hopf curves define the transition from the model's steady-state (lightest grey) into the periodic solutions (darker grey). In both parts of the figure, we find that an increase in delay τ leads to a decrease in the value of νse at the HB. Thus, if the delay in both models is increased, a weaker coupling strength νse is required to initiate the oscillatory behaviour. If τ is decreased below τ≃0.04 in figure 4a or τ≃0.02 in figure 4b, we end up in a region where both models do not support the spike–wave solutions.
Other features of figure 4a,b are qualitatively different; in figure 4a, we find a clear relationship between the GABAB delay τ and the number of spikes observed. An increase in τ leads to a stepwise increase of the number of spikes. The spiking regions are separated by black curves of ghost bifurcations. Crossing such a curve means that a spike is added or removed from the model solution. Conversely, figure 4b displays a different scenario. There exists a region of spike waves between νse≈3.8×10−3 and 6.2×10−3 V s. At approximately τ≈0.08 s, another region labelled ‘two-spike wave’ exists, bounded by a curve of the ghost bifurcations. The solution ϕe(t) contains a wave, followed by a small spike, a plateau of constant ϕe, and a larger spike. A similar region can be seen at (νse=3×10−3, τ=0.16). The original cortico-thalamic model is mainly capable of qualitatively reproducing single spike–wave solutions.
If τ is fixed and νse is increased in figure 4, eventually one of the grey branches will be crossed. This implies that the (poly)spike solutions either lose stability via a PD bifurcation or vanish completely in a saddle-node bifurcation (SL). A detailed analysis of the PD bifurcations in figure 4b is given in the recent work of Rodrigues et al. (2008). It was found that small isles of the stable periodic solutions and small chaotic regions exist, on the right of the (PD)-curve. We have found similar dynamics in figure 4a. These small regions of more complex dynamics are beyond the scope of the present study and remain under investigation.
How can we relate the activity map in figure 4a back to the time series of our mean-field model? Figure 5 provides a schematic overview. Figure 5a(ii)–c(ii) consists of two parts: figure 5a(i)–c(i) is a small section taken from figure 4a. Below each section, we display a bifurcation diagram in νse. As an example, in figure 5a(ii), we display a subsection of figure 4a at approximately τ=0.06 s. Hence, if the delay of the slow GABAB coupling is fixed at 0.06 s, and we increase νse stepwise, the diagram in the lower part of (a(ii)) is obtained. In figure 5a(ii)–c(ii), ϕe(t) develops into periodic 2–4 Hz (poly)spike patterns, which eventually destabilize via the saddle-node (SL) or the PD bifurcations. The latter bifurcation occurs in a tiny window of the parameter space, and is followed by an even tinier chaotic attractor. These features could only be captured by DDE-BIFTOOL and not our numerical integrator. Past the (SL) and (PD) points, ϕe(t) then converges to a steady state ϕe=250 s−1. This steady state, which coexists with all the periodic solutions of our model, corresponds to the maximal neuronal firing rate at .
Some of the various periodic (poly)spike patterns of ϕe(t) are given in figure 5a(iii–v)–c(iii–v). For example, the time series in figure 5a(iii–v) correspond to a fixed connection delay τ=0.06 s, and, from (iii–v), νse=0.0017, 0.002 and 0.0024 V s. Similarly, figure 5b(iii–v),c(iii–v) corresponds to τ=0.10 and 0.16 s, where νse=0.0017, 0.002, 0.0026 V s and νse=0.0016, 0.00165, 0.002 V s from (iii–v), respectively. These time series correspond to sections of the bifurcation diagrams in figure 5a(ii)–c(ii). We find that an increase in GABAB delay not only results in more spikes, but also a decrease in frequency. At τ=0.06, we find patterns with a frequency of approximately 3 Hz. For τ=0.10 and 0.16 the frequency decreases to 2.3 and 1.8 Hz, respectively. The first two frequencies lie in the range of the absence seizure EEG, of which the main frequency component is situated between 2 and 4 Hz, respectively.
Our final step is to make a link between the above-mentioned time series, and the real patient EEG data. Some examples are shown in figure 6. In §3, we pointed out that ϕe can be approximated to be linearly related to the EEG voltages. An often used convention is that the constant of proportionality is taken to be negative (Robinson et al. 2002). Figure 6a shows the EEG traces from our seizure database, taken from three different patients with CAE. Figure 6b displays various time traces of our model, using the activity map in figure 4 to determine the parameter regimes. In each case, we find an interesting qualitative agreement between the model and the data, suggesting that there is a significant variation in the value of the delay τ between the patients. It remains an important question of how to interpret the biophysical meaning of such a variation.
The work presented in this paper was motivated by a desire to understand the mechanisms responsible for transitions between the different types of seizure activity observed in the EEG recordings from human subjects with absence seizures. Studying a database of 48 seizures from 20 subjects, we observed not only the classical single spike and wave activity, but a wide variety of polyspike and wave solutions (up to three spikes per cycle). Building on our earlier work in this area (Breakspear et al. 2006; Rodrigues et al. 2006), we isolated the transition to single spike and wave activity, identifying an inflection point in the system dynamics as causing the spike in the model (Rodrigues et al. 2008). Subsequently, we investigated the addition of further spikes per cycle, but observed that the successive inflection points occurring in the original model of Robinson et al. (2002) did not accurately reproduce the data recordings from human subjects (figure 4b).
Consequently, we developed a modified cortico-thalamic model, which incorporated the differences in GABAergic channels in the connections between the inhibitory and the excitatory thalamic neuronal populations. Studying the transitions in this model, using a combination of bifurcation theory and numerical continuation methods, we were able to map out in (νse, τ)-space branches of the inflection points giving rise to different numbers of spikes per cycle. These activity maps demonstrate that the delay plays a crucial role in these transitions. A HB was observed to give a transition between the healthy activity and the oscillatory behaviour preceding the onset of spike and wave activity. From the analysis of the model, we were able to obtain time series from our model that mimic those recorded from human subjects (figure 5). In this regard, a further important question is to understand the mechanisms that control the frequency of oscillations within each window corresponding to the solutions with different numbers of spikes. We are currently investigating this phenomenon and the results will be presented elsewhere.
These findings also suggest a number of hypotheses of potential clinical relevance that should be investigated. Firstly, the onset of the seizure in our model corresponds to the transition between the steady-state and the oscillatory activity, rather than the onset of spike–wave behaviour. This suggests that clinical symptoms of the absence seizures may occur before the observation of spike–wave activity in the EEG recordings. This could potentially be investigated using a continuous performance task in human subjects with absence seizures, such as following a moving trace on a screen. Secondly, the variation in delay between the patients as observed when fitting model time series to the data suggests that it is important to establish a mapping between the macroscopic models and the microscopic parameters. Finally, the techniques presented here provide a means for predicting when a seizure might occur, by estimating and tracking parameters of the model implicated in seizure onset from the clinical data. For such an approach to be successful, the model must be capable of accurately reproducing a wide variety of clinically observed seizure types (exemplars of this can be observed in figure 6). Techniques for the tracking of bifurcations from the experimental data have recently been proposed and validated in an electronic circuit by Sieber & Krauskopf (2008) and this approach is a strong candidate to be used in this case. This would provide an alternative method for seizure prediction than more traditional methods based on signal processing (Mormann et al. 2005).
All authors acknowledge financial support from the EPSRC via grant no. EP/D068436/01 and the Leverhulme Trust via the Theoretical Neuroscience Network. Useful conversations with Gabor Stepan are also gratefully acknowledged.
One contribution of 10 to a Theme Issue ‘Delay effects in brain dynamics’.
- © 2009 The Royal Society