## Abstract

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.

## 1. Introduction

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 GABA_{A} and GABA_{B} 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 pioneering work of Lopes da Silva

*et al*. (1974) and Freeman (1975), who developed the equations of motion, describing the behaviour of neural masses, based on detailed experimental studies.The differences in the time scales of activation/inactivation due to the GABA

_{A/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 *V*_{a}(** r**,

*t*), the average firing rate

*ς*

_{a}(

**,**

*r**t*) and the axonal field

*ϕ*

_{a}(

**,**

*r**t*) at position

**and time**

*r**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 *V*_{a} of a neural mass, under incoming postsynaptic potentials *P*_{a}. Equation (3.2) describes the average firing rate of the mass. If *V*_{a} 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 GABA_{A} and GABA_{B} synapses. Since the GABA_{B} functions via second messenger processes, its postsynaptic currents develop on a slower time scale than the currents of GABA_{A} (figs 5.8 and 5.9 in Destexhe & Sejnowski (2001)). Hence, we expect that, if thalamic neurons fire, the GABA_{B}-mediated connection to the specific relay nuclei will be delayed with respect to GABA_{A}. 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 *V*_{i}≈*V*_{e} and *ς*_{i}≈*ς*_{e} (Robinson *et al*. 2002). Furthermore, *a*_{a} 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.

## 4. Results

Our results are presented in §4*a*,*b*. In §4*a*, 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 §4*b*, 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 1*a*. These are the waveforms that are classically associated with the absence seizures. More general polyspike patterns, such as the ones displayed in figure 1*b*,*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 GABA_{B} coupling, we choose a delay of *τ*=100 ms. This represents the order of magnitude of rise and decay time of GABA_{B} currents (Destexhe & Sejnowski 2001, figs 5.8 and 5.9). The bifurcation diagram is shown in figure 3*a*.

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 3*a*,*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 3*c*. 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 3*b*). 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 GABA_{B} coupling? And what happens to *ϕ*_{e}(*t*) as *ν*_{se} is increased further?

### (b) Continuation in two parameters

The analysis presented in §4*a* 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 4*a* displays the various bifurcation curves, computed for our mean-field model (A 2). In figure 4*b*, 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 GABA_{B} synapse, the model in figure 4*b* assumed that the delay arises owing to synaptic transmission times for the signals travelling between the cortex and the thalamus.

In figure 4*a*,*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 4*a* or *τ*≃0.02 in figure 4*b*, we end up in a region where both models do not support the spike–wave solutions.

Other features of figure 4*a*,*b* are qualitatively different; in figure 4*a*, we find a clear relationship between the GABA_{B} 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 4*b* 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 4*b* 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 4*a*. 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 4*a* back to the time series of our mean-field model? Figure 5 provides a schematic overview. Figure 5*a*(ii)–*c*(ii) consists of two parts: figure 5*a*(i)–*c*(i) is a small section taken from figure 4*a*. Below each section, we display a bifurcation diagram in *ν*_{se}. As an example, in figure 5*a*(ii), we display a subsection of figure 4*a* at approximately *τ*=0.06 s. Hence, if the delay of the slow GABA_{B} 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 5*a*(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 5*a*(iii–v)–*c*(iii–v). For example, the time series in figure 5*a*(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 5*b*(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 5*a*(ii)–*c*(ii). We find that an increase in GABA_{B} 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 6*a* shows the EEG traces from our seizure database, taken from three different patients with CAE. Figure 6*b* 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.

## 5. Discussion

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 4*b*).

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

## Acknowledgments

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.

## Footnotes

One contribution of 10 to a Theme Issue ‘Delay effects in brain dynamics’.

- © 2009 The Royal Society