## Abstract

The stochastic partial differential equations (SPDEs) stated by Steyn-Ross and co-workers constitute a model of mesoscopic electrical activity of the human cortex. A simplification in which spatial variation and stochastic input are neglected yields ordinary differential equations (ODEs), which are amenable to analysis by techniques of dynamical systems theory. Bifurcation diagrams are developed for the ODEs with increased subcortical excitation, showing that the model predicts oscillatory electrical activity in a large range of parameters. The full SPDEs with increased subcortical excitation produce travelling waves of electrical activity. These model results are compared with electrocortical data recorded at two subdural electrodes from a human subject undergoing a seizure. The model and observational results agree in two important respects during seizure: (i) the average frequency of maximum power, and (ii) the speed of spatial propagation of voltage peaks. This suggests that seizing activity on the human cortex may be understood as an example of pathological pattern formation. Included is a discussion of the applications and limitations of these results.

## 1. Introduction

Over 50 million people worldwide suffer from epilepsy, a debilitating condition of chronic unprovoked seizures. There exist many different classifications of seizures (e.g. simple partial seizures, complex partial seizures, generalized seizures), yet all seizures share the general characteristic of inducing a hyperexcitation of the cortex (the outer few millimetres of the brain). The microscopic behaviour of individual cortical neurons during this hyperexcitation is well documented and is the basis for some drug treatments. However, comprehensive microscopic recordings are infeasible, especially in human patients. Therefore, physicians record the macroscopic electrical behaviour of the seizing human cortex and collect the scalp electroencephalogram (EEG) and electrocorticogram (ECoG). At this macroscopic level, which involves interactions between millions of individual neurons, cortical seizures are less well understood (Quyen *et al*. 2003).

In this paper, we describe a model of the mesoscopic electrical behaviour of the seizing human cortex. We employ a recently developed continuum model of cortical electrical activity, which consists of a system of stochastic partial differential equations (SPDEs; Steyn-Ross *et al*. 2003). These SPDEs (and the related ordinary differential equations; ODEs) have been successfully applied to study the electrical properties of the anaesthetized human cortex (i.e. the *anaesthetodynamic model* of cortical function). In a series of papers, Steyn-Ross and colleagues showed that the model equations predict changes in the scalp EEG of anaesthetized patients consistent with experimental results (Steyn-Ross *et al*. 1999, 2001*a*,*b*, 2003).

Here, we consider whether these SPDEs can model the mesoscopic electrical behaviour recorded from the human cortex during a seizure. Our methods differ from previous discussions. We do not assume that some variables equilibrate much faster than others (the slow-membrane or adiabatic approximation; Steyn-Ross *et al*. 1999), nor do we apply techniques of stochastic calculus (such as the Ornstein–Uhlenbeck equation) to the linearized SPDEs (Steyn-Ross *et al*. 2003). Instead, we approach the SPDEs from a dynamical systems perspective and employ ideas and tools from bifurcation theory. We show that in some parameter regimes, the SPDEs predict pathological oscillatory behaviour—characterized by a bifurcation to a large-amplitude, coherent travelling wave against a background of spatio-temporal fluctuations—in some important ways consistent with ECoG recordings from the seizing cortex.

The organization of this paper is as follows. In §2, we discuss clinical ECoG data recorded at two subdural electrodes from a human subject during his typical complex partial seizures. In §3, we review the model equations, the dimensionless SPDEs we use here, and the variable of principal interest: the excitatory mean soma potential *h*_{e}. In §4, we consider a simplified version of the SPDEs that possesses neither spatial variance nor stochastic input; we call this system the dimensionless ODEs. In §§4.1 and 4.2, we compute bifurcation diagrams for the dimensionless ODEs at two different parameter values. We show that in both cases oscillatory activity (as expected during a seizure) occurs. For typical parameter values, the oscillatory activity is unstable and short-lived. For other parameter values (which may be consistent with the seizing cortex), the oscillatory activity is stable and long-lived. In §5, we analyse the complete dimensionless SPDEs. We show that results from this system agree qualitatively with results from the simpler dimensionless ODEs, and that travelling wave solutions exist. In §6, we compare the model solutions with the experimental recordings, and in §7, we discuss these results.

## 2. Observational data: ECoG seizure recordings

The cellular mechanisms of a seizure have been extensively studied and some general microscopic characteristics have been deduced (Dichter & Ayala 1987). Preceding a seizure, thousands of individual neurons in the seizure focus (the brain region where a seizure begins, also known as the epileptogenic zone) undergo depolarization shifts followed by an afterhyperpolarization. As long as this behaviour is confined to the seizure focus, there may be no clinical manifestation (although this synchronous activity can be detected as an interictal spike or a sharp wave in the EEG or ECoG). Gradually, as the seizure develops, the magnitude of the afterhyperpolarization decreases and individual neurons generate nearly continuous action potentials. The inhibition surrounding the seizure focus weakens, the seizure spreads to other cortical neurons, and a clinical seizure occurs. Here, we will not consider this microscopic behaviour of individual neurons, nor will we consider ictogenesis (the initiation of the seizure, which may occur in deeper brain regions). Instead, we investigate the mesoscopic characteristics of seizure activity on the cortical surface as made manifest in subdural ECoG recordings.

Large amounts of ECoG data are often recorded from the seizing cortex (i.e. ictal data) as part of clinical care of patients with intractable epilepsy. These patients whose seizures do not respond well to drug treatments may undergo surgery to remove a region of cortex that is believed to contain the seizure focus. Such resective surgery offers a chance to eliminate or ameliorate the seizures. In the process of planning such surgery, clinicians must accurately locate the seizure focus, and identify any surrounding areas of functional cortex that should be preserved to minimize post-operative deficits. If necessary, to define the seizure focus more precisely with respect to structural and functional anatomy, clinicians may implant subdural electrodes on the brain surface for extended time periods (typically one week) to obtain recordings during seizures and locate the seizure onset. The electrodes are placed over the region of cortex suspected to contain the seizure focus based on: (i) semiologic features of the patient's seizures (i.e. the observable features of the seizure, such as eye deviation or limb jerking), (ii) scalp EEG recordings made during typical seizures, (iii) abnormalities observed in the ECoG recorded intraoperatively at the time of electrode implantation (i.e. interictal recordings), or (iv) abnormalities seen on brain images. After the electrodes are implanted and the patient experiences his or her typical seizures, physicians locate the seizure focus through visual analysis of the ECoG recordings and formulate the surgical plan.

In this paper, we consider data collected from a 49-year-old male who had electrode implantation as part of his care at the University of California, San Francisco (UCSF) Epilepsy Center. He had medically refractory complex partial seizures that began with staring and manual automatisms and frequently generalized into convulsions. Prior scalp EEG recordings had suggested that his seizures originated in the right temporal lobe, but MRI showed left mesial temporal sclerosis, raising the possibility that his seizures arose from the left temporal lobe instead. To better lateralize seizure onset, subdural and depth electrodes were implanted, with subdural electrode strips placed over the left and right lateral frontal regions and the left and right subtemporal regions. In addition, bilateral hippocampal depth electrodes were placed, although we will not discuss these hippocampal recordings here. Each subdural electrode strip contained six 4 mm diameter platinum–iridium discs arrayed in a single row, embedded in a 1.5 mm thick silastic sheet with 2.3 mm diameter exposed surfaces and 10 mm spacing between the discs. After electrode implantation, the patient was brought to the video–telemetry unit, where his antiseizure medications were slowly withdrawn. ECoG data were recorded continuously at 400 Hz from all electrodes for 5.5 days. During this time, observations of eight typical seizures were captured. These recordings were of good quality and were determined, based on clinical review by board-certified clinical neurophysiologists, to show seizure onset from right medial temporal regions, with a pattern characteristic of seizures arising from this region. The patient went on to have a right anterior temporal lobectomy based on these clinical data.

ECoG epochs containing six of the patient's seizures were extracted from the clinical record and reviewed for research purposes in accordance with UCSF and University of California, Berkeley human subjects guidelines. (Two of the seizure data files were corrupted and no longer available for extraction.) In figure 1, we show ECoG data leading up to and during a seizure recorded at two neighbouring subdural electrodes above the right lateral frontal region. We refer to the time-series recorded at these two electrodes as *X* (lower curve) and *Y* (upper curve). The ECoG data shown in this figure possess three notable features. First, we observe that large-amplitude voltage oscillations occur during the seizure (*t*>17.5 s). We determine the frequency of these oscillations by computing the windowed power spectrum (WPS) of *X* and *Y*. To compute the WPS, we partition the 50 s of data plotted in figure 1 into 100 overlapping windows of 1.0 s duration and 0.5 s overlap. For example, the first window includes data for 0.0 s≤*t*≤1.0 s, the second 0.5 s≤*t*≤1.5 s, the third 1.0 s≤*t*≤2.0 s and so on. We then multiply the data in each window by the Hanning function and compute the power spectrum. We plot the WPS results for *X* and *Y* in figure 2*a*,*b*, respectively. The horizontal axis in each figure corresponds to the centre times of the windows, and the solid vertical line to the time of seizure onset. We note, in the seconds preceding the seizure, the suppression of power at both electrodes for all frequencies. After seizure onset, the regions of largest power occur at frequencies between 1 and 10 Hz. The second notable feature in the observed data is the abrupt transition from normal ECoG activity (for *t*<14 s) to seizing activity (for *t*>17.5 s). This can be observed both in the original time-series data (figure 1*a*) and in the WPS (figure 2). Lastly, we find that the oscillations recorded at the two electrodes are slightly out of phase. One may deduce this conclusion directly from figure 1*b*; we note that peaks in the oscillations in *X* (lower trace) appear to precede those in *Y* (upper trace) until *t*≈29 s. At this point, the relationship between the two time-series becomes more complicated. To determine how the phase relationship between these two electrodes changes over time, we compute the windowed cross-correlation (WCC). In a manner similar to the WPS calculation, we partition the data into overlapping windows and compute the cross-correlation between *X* and *Y* in each window. We show the results in figure 3. Following seizure onset (denoted by the vertical line in the figure), two regions of strong correlation occur: the first from 17.5 s<*t*<27 s, and the second from 28 s<*t*<50 s. In the first region, *X* leads *Y* (i.e. the lag time is positive) while in the second region, *X* lags *Y* (i.e. the lag time is negative but near zero). In §6, we quantify the WPS and WCC results for six seizures recorded from the subject. We compare these results with the solutions to the SPDEs, discussed in §3.

## 3. Model: dimensionless SPDEs

An ideal model of the human cortex would describe the electrical behaviour of each individual neuron and its surrounding extracellular environment. This discrete model—applied over the entire three-dimensional cortex—would contain over 10^{10} dynamical variables and would be intractable. Fortunately, the physiology of the cortex (e.g. dense local connections) and numerous experimental results suggest that neurons act in populations or assemblies (Singer 1993; Singer & Gray 1995). Moreover, most methods for observing the human cortex (e.g. the ECoG time-series of interest in this work) record the summed electrical activity from populations of order 10^{5} cortical neurons (Nunez 1995). For these reasons, researchers have developed continuum models of cortical electrical activity.

Such continuum models are not new; some of the earliest were developed by Freeman (1964) and Wilson & Cowan (1972, 1973). Here, we utilize the system of SPDEs stated by Steyn-Ross *et al*. (2003). These SPDEs were introduced in Steyn-Ross *et al*. (1999), where the authors included stochastic input to the system of PDEs first stated by Liley *et al*. (1999). To develop this model, Liley *et al*. (2002) considered a spatially averaged approximation to the human cortex. Microscopic elements (like individual neurons) were averaged over columnar volumes (perpendicular to the cortical surface) whose diameter was chosen to lie below the spatial resolution of EEG or ECoG recordings. The resulting mesoscopic model (with the addition of stochastic input) consists of 14 coupled, nonlinear SPDEs (Steyn-Ross *et al*. 2003).

By solving the SPDEs in Steyn-Ross *et al*. (2003) numerically, one computes solutions for all 14 variables as functions of space and time. One of these variables, *h*_{e}, is the spatially averaged excitatory soma membrane potential. Researchers have demonstrated that the deviation of *h*_{e} from rest is proportional to the sign-reversed value of the extracellular local field potential (LFP). Because the ECoG represents the spatially averaged LFP, we assume that *h*_{e} is linearly related to the ECoG (Liley *et al*. 2002). In this way, the model variable *h*_{e} is related to the observational ECoG data. In what follows, we compare *h*_{e} calculated in numerical solutions to the SPDEs with ECoG data recorded during seizure. We show that increasing the subcortical excitatory input to the model cortex produces behaviour in *h*_{e} that mimics ECoG data recorded from the seizing cortex.

In appendix A, we discuss a procedure to non-dimensionalize the SPDEs and associated functions. The main advantage of recasting the equations in dimensionless form is a reduction in the number of parameters. There are 20 parameters in the dimensionless SPDEs whereas in the original SPDEs there are 29. We define each dimensionless variable and parameter in terms of its dimensional counterparts from Steyn-Ross *et al*. (2003) in tables 1 and 2, respectively. In what follows, we study a simplified ODE version of the dimensionless model in §4, and compute numerical solutions to the complete dimensionless SPDE model in §5. Readers less interested in the mathematical details may like to skip to §6 where we compare the model solutions and observational results.

## 4. Simulations: dimensionless ODEs

Our goal is to determine whether the dimensionless SPDEs in equation (A 1*a*–*h*) can be used to model the mesoscopic electrical activity observed on the seizing cortical surface, as discussed in §2. However, the size of the system (14 first-order differential equations and 20 parameters), the stochastic input in equation (A 1*c*–*f*) and the spatial dependence in equation (A 1*g*,*h*) make this a challenging model. Therefore, to gain what insight we can from a simpler model, we ignore for the moment the spatial dependence and stochastic input in equation (A 1*a*–*h*). The resulting equations form a system of ODEs, which we call the dimensionless ODEs. Although there are many tools available for the study of ODEs, we use Auto (continuation and bifurcation software for ODEs) to determine fixed points, their stability type, limit cycles and bifurcations in the phase portraits (Doedel *et al*. 2000).

During a seizure, the cortex typically enters a state of hyperexcitation, manifest in ECoG data through large-amplitude oscillations (Niedermeyer & DaSilva 1999). This activity corresponds to spatially coherent oscillations in the variable *h*_{e}. In this direction, we investigate whether oscillations in *h*_{e} occur in the dimensionless ODEs owing to changes in the excitatory parameters. We note that different types of seizures produce different oscillatory patterns in EEG and ECoG recordings (Spencer *et al*. 1992; Ebersole & Pacia 1996). Here, we seek to model the seizing cortex in a qualitative way suggestive of the typical ECoG data shown in figure 1. Specifically, we require that: (i) the model produces stable oscillations in *h*_{e}, (ii) the frequency of the oscillations agree (roughly) with clinical observations, and (iii) the transition to oscillatory behaviour occurs abruptly.

The general statement that increased excitation (or decreased inhibition) incites cortical seizure activity masks the numerous associated physiological changes that occur in the seizing neuronal assemblies (Dichter & Ayala 1987; Morimoto *et al*. 2004). At the cellular level, these physiological changes are observed in experiment and can be compared to results computed from detailed computational models of a single neuron (Traub *et al*. 2001). It is not clear how the cellular mechanisms or single neuron models relate to mesoscopic seizure recordings or continuum models. Seizures induced in animal models allow mesoscopic ECoG recordings and some control over parameters related to continuum models. For example, Freeman (1992) discusses electrocortical data recorded from a seizing animal's olfactory bulb. He compares these recordings with his KIII model of the olfactory system and finds that a parameter connecting excitatory subsets of the model must increase to induce seizures. Analogies between such animal models and human seizure activity can be made but, as for the cellular models, these relationships are not clear.

In our analysis, we vary only two parameters, *P*_{ee} and *Γ*_{e}, both related to the excitation of the model. We choose these parameters for two reasons. First, as mentioned previously, the general claim that increased excitation incites seizures is well known. Thus, we select two parameters related to the excitation of the model. Second, an increase in the level of the membrane potential of a neuronal population is thought to be an important control factor in inducing seizures (Dichter & Ayala 1987; da Silva *et al*. 2003). We show that increases in *Γ*_{e} raise the excitatory mean soma potential *h*_{e} of the stable fixed points. A similar result holds for *P*_{ee} but is not shown. We fixed the remaining 18 dimensionless parameters at the typical values shown in table 2. In terms of dimensional variables, the dimensionless parameter *P*_{ee} is: (i) directly proportional to *p*_{ee} (the subcortical excitatory spike input to the excitatory neurons of the cortex), and (ii) inversely proportional to *S*^{max} (the maximum firing rate induced by the soma voltage). Thus, an increase in *P*_{ee} represents either an increase in the subcortical excitation of the cortex or a decrease in the maximum firing rate. To model an increase in excitatory input (from a subcortical region, such as the thalamus) to the cortex, the parameter *P*_{ee} is increased. The dimensionless parameter *Γ*_{e} is: (i) directly proportional to *S*^{max} and *G*_{e} (the peak excitatory postsynaptic potential; EPSP), and (ii) inversely proportional to *γ*_{e} (the EPSP neurotransmitter rate constant) and (the magnitude of the difference between the excitatory reversal and resting potentials). Thus, an increase in *Γ*_{e} represents either an increase in peak EPSP amplitude or maximum firing rate, or a decrease in the difference between the excitatory reversal and resting potentials or the EPSP neurotransmitter rate constant. The latter corresponds to an increase in the EPSP duration. We note that the typical values of *P*_{ee} and *Γ*_{e} are 11.0 and 1.42×10^{−3}, respectively (see table 2).

It is perhaps easier to interpret the parameters *P*_{ee} and *Γ*_{e} directly from equation (A 1*a*–*h*). The parameter *P*_{ee} appears only in equation (A 1*c*). This parameter controls the strength of the (non-stochastic) excitatory subcortical input to excitatory neurons in the cortex. The parameter *Γ*_{e} appears in equation (A 1*a*,*b*), and controls the influence of excitatory input on the mean soma membrane potentials. For example, when *Γ*_{e} is large, the effect of the excitatory inputs ( and ) on the potentials ( and ) is increased. We will show that solutions of the model equations agree qualitatively with ECoG data recorded from the seizing cortex when *P*_{ee} is dramatically increased, and *Γ*_{e} is slightly decreased.

In what follows, we compute bifurcation diagrams and numerical solutions of the dimensionless ODEs at two values of *P*_{ee}. We first consider the results for *P*_{ee}=11.0 (the typical value) and show that oscillations in *h*_{e} occur but are unstable and short lived. We then show that for *P*_{ee}=548.066 (nearly 50 times the typical subcortical excitation), *h*_{e} undergoes large-amplitude, stable oscillations, and that similar oscillations occur over a wide range of parameter values. The results in each case are compared with the clinical data discussed in §2.

### 4.1 Example: dimensionless ODEs at *P*_{ee}=11.0

In the first example, we fix *P*_{ee} at its typical (dimensionless) value of 11.0 and vary the parameter *Γ*_{e}. In figure 4*a*, we plot *h*_{e} for the fixed points of the dimensionless ODEs versus the parameter *Γ*_{e}. Note that in this figure, we plot the dimensional variable *h*_{e} (with units mV), which is related to by the scale factor . The solid lines and dashed line in figure 4*a* correspond to the stable and unstable fixed points of the dimensionless ODEs, respectively. We note that an increase in *Γ*_{e} produces an increase in both curves of the stable fixed point value of *h*_{e} (i.e. increased steady-state values of the spatially averaged excitatory soma membrane potential).

The S-shape and stability of the fixed points is similar to that discussed in Steyn-Ross *et al*. (1999, 2003). In those works, the authors varied the (dimensional) inhibitory neurotransmitter rate constant *γ*_{i}. Here, the parameter *Γ*_{e} is inversely proportional to the (dimensional) excitatory neurotransmitter rate constant *γ*_{e}. We have found, but do not show, similar bifurcation diagrams for the *dimensional* ODEs with *γ*_{i}→*γ*_{i}/*λ* or *γ*_{e}→*γ*_{e}*λ* and parameter *λ* varied between 0.1 and 1.5.

In addition to the saddle–node bifurcations of fixed points, the dimensionless ODEs at *P*_{ee}=11.0 also undergo a Hopf bifurcation at *Γ*_{e}=1.20×10^{−3}. (We note that this represents a 15% decrease in the typical value of *Γ*_{e}=1.42×10^{−3}.) We mark this Hopf bifurcation in figure 4*a* with an asterisk. We are particularly interested in Hopf bifurcations because, at a Hopf bifurcation, the dynamics of *h*_{e} can change from stationary behaviour to oscillatory behaviour. To illustrate the oscillatory behaviour of *h*_{e} near the Hopf bifurcation, we compute a numerical solution to the dimensionless ODEs near the Hopf bifurcation at *Γ*_{e}=1.21×10^{−3} and *P*_{ee}=11.0. We choose the initial conditions so that the dynamics begin just outside the basin of attraction of the stable fixed point and compute the trajectory using a fourth-order Runge–Kutta method with a time-step of 0.4 ms. We plot in figure 4*b* (dimensional) *h*_{e} as a function of dimensional time *t*. For 0 s<*t*<9.5 s, *h*_{e} oscillates at approximately 8 Hz. The amplitude of the oscillations steadily increases until *t*=9.5 s, at which point *h*_{e} abruptly moves to the stable fixed point near −84 mV. We show this transition with finer resolution in figure 4*c*.

We have shown that for *P*_{ee}=11.0 the dynamics of *h*_{e} undergo a Hopf bifurcation at *Γ*_{e}=1.20×10^{−3}. Near the Hopf bifurcation there exist transient oscillations in *h*_{e} that increase in amplitude until the dynamics undergo a transition to a stable fixed point near −84 mV. We have not shown but note that adding noise to the system decreases the time of the transient oscillations. These unstable oscillations in *h*_{e} do not mimic the electrical activity of the seizing cortex where the oscillations maintain a large amplitude and are necessarily stable to incessant perturbations from other cortical, as well as deeper, brain regions.

### 4.2 Example: dimensionless ODEs at *P*_{ee}=548.066

Here, we consider an example more closely related to the seizing cortex. To increase the excitation of the model cortex, we fix *P*_{ee}=548.066 (nearly 50 times the typical value). This can be interpreted as increased excitatory input from deeper brain regions to the cortex, say. As in §4.1, we vary the parameter *Γ*_{e} and plot in figure 5*a* the bifurcation diagram in *h*_{e} for the dimensionless ODEs. The solid lines and dashed line in figure 5*a* correspond to stable and unstable fixed points of the dimensionless ODEs, respectively, and the asterisks to Hopf bifurcations. We note that an increase in *Γ*_{e} produces an increase in both curves of stable fixed points of *h*_{e}.

There are additional curves in figure 5*a* absent from figure 4*a*. These curves indicate the extremal values and stability type of the limit cycles born in the Hopf bifurcations. The dash–dot lines in figure 5*a* represent the maxima and minima of *h*_{e} achieved during a *stable* limit cycle. The dotted lines represent the maxima and minima of *h*_{e} achieved during an *unstable* limit cycle. Thus, for 0.67×10^{−3}<*Γ*_{e}<0.96×10^{−3}, the dynamics of *h*_{e} feature a stable limit cycle of large amplitude. We note that a 32% decrease in the typical value of *Γ*_{e}=1.42×10^{−3} is required to enter this range. The stable, oscillatory behaviour of *h*_{e} satisfies one of our criteria for modelling the seizing cortex.

To illustrate the stable oscillations in *h*_{e} for 0.67×10^{−3}<*Γ*_{e}<0.96×10^{−3}, we compute a numerical solution to the dimensionless ODEs at *Γ*_{e}=0.96×10^{−3} and *P*_{ee}=548.066 using the fourth-order Runge–Kutta method with a time-step of 0.4 ms. In figure 5*b*, we plot *h*_{e} (with dimensions mV) as a function of dimensional time *t*. After an initial transient, *h*_{e} is entrained by a large-amplitude limit cycle with a dominant frequency near 7.5 Hz. We note that these oscillations are not sinusoidal. The frequency of the *h*_{e} oscillations is qualitatively consistent with the clinical observations from a seizing patient shown in figure 1.

The true environment of the human cortex is continually changing, for example, in response to sensory stimuli. Therefore, any model of the seizing cortex must behave properly over a broad range of parameter values. Here, we require that oscillatory activity in *h*_{e} occur over extended regions of the parameters *P*_{ee} and *Γ*_{e}. To confirm this, we compute numerical solutions to the dimensionless ODEs for 11.0<*P*_{ee}<1000.0 and 0.5×10^{−3}<*Γ*_{e}<1.4×10^{−3} using the fourth-order Runge–Kutta method with a time-step of 0.4 ms. We then determine the difference between the maximum and minimum achieved by the solution *h*_{e} after transient behaviour has decayed. If *h*_{e} approaches a fixed point, then the maximum and minimum are nearly equal and their difference approaches zero. However, if *h*_{e} is entrained by a limit cycle (e.g. figure 5*b*), then the difference between the maximum and minimum achieved by *h*_{e} is non-zero. In figure 6, we plot the difference between the maximum and minimum achieved by *h*_{e} as a function of the parameters *P*_{ee} and *Γ*_{e}. The difference is plotted in greyscale with white representing a 0 mV difference and black representing a 50 mV difference. We find that oscillations in *h*_{e} (represented by the dark regions in figure 6) extend over a broad range of parameter values beginning near *P*_{ee}=250.0 and *Γ*_{e}=1.3×10^{−3}. These regions of oscillatory activity in *h*_{e} illustrate the parameter values at which the dimensionless ODEs ‘seize’. We have found that the Hopf bifurcation (and thus the oscillatory activity) is born in a codimension two bifurcation of the Takens–Bogdanov type at *P*_{ee}=53.2 and *Γ*_{e}=2.98×10^{−3} (Kramer 2005). We note that increasing *P*_{ee}, and thus raising the subcortical excitatory input to the model cortex, enlarges the region of *Γ*_{e} over which oscillations in *h*_{e} occur.

## 5. Simulation: dimensionless SPDEs

We have considered in some detail the dimensionless ODEs. There is no direct conclusion that the analysis of this simplified system allows one to draw for the full dimensionless SPDEs. The stochastic inputs or the spatially distributed dynamics may destroy the interesting features we observed in the dimensionless ODEs. We have shown in §4.2 that the dimensionless ODEs undergo a Hopf bifurcation near *P*_{ee}=548.066 and *Γ*_{e}=0.96×10^{−3}, and that oscillations in *h*_{e} persist over a broad region of parameter values. Here, we shall examine whether the full dimensionless SPDEs exhibit similar oscillatory behaviour. We will consider three examples. For each example, we compute numerical solutions to the dimensionless SPDEs using the Euler–Maruyama algorithm with fixed steps in space and time, 14 mm and 0.1 ms, respectively (Higham 2001). We consider only one spatial dimension ; the system can now be visualized as describing a line of closely spaced electrodes, such as the strip of subdural electrodes described in §2. In each example, we enforce periodic boundary conditions in space, fix *Γ*_{e} and *α* (representing the noisy input from subcortical sources to the cortex) to be uniform in and , and vary the spatial and temporal dependence of *P*_{ee}. In the first example, we set *P*_{ee} to be uniform in space and time. In the second, we set *P*_{ee} to be Gaussian in space and constant in time, while in the third, we set *P*_{ee} to be Gaussian in space and non-uniform in time. In §6, we compare the results of these simulations and the clinical seizure data.

We start by computing numerical solutions to the dimensionless SPDEs at *Γ*_{e}=1.04×10^{−3} and *P*_{ee}=548.066, near the parameter values used in §4.2. In this example, both *P*_{ee} and *Γ*_{e} are uniform in space and time. We note that, at these parameter values, the noise-free dimensionless ODEs are near the Hopf bifurcation. In figure 5*a*, the simplified system would be just to the right of the rightmost asterisk. In figure 7*a*, we show the solution to the dimensionless SPDEs with *α*=0.001. In this figure, we plot (dimensional) space on the horizontal axis, (dimensional) time on the vertical axis, and the value of *h*_{e} in linear greyscale; white corresponds to *h*_{e}=−65 mV, and black to *h*_{e}=−45 mV. From this figure it is clear that the value of *h*_{e} is approximately constant (near −53 mV, the fixed point of the dimensionless ODEs) in both space and time.

In figure 7*b*, we show the solution to the dimensionless SPDEs with *α*=0.002. Here, we have increased the strength of the spatially uniform stochastic input, and plotted *h*_{e} in linear greyscale with white corresponding to *h*_{e}=−100 mV, and black to *h*_{e}=0.0 mV. Otherwise, the parameters, initial conditions and plotting scheme are identical to those used to compute and display the previous numerical solution. There are two interesting features to note in figure 7*b*. First, travelling waves are engendered in *h*_{e}. In the figure, these are represented by repeated white and black ridges inclined with respect to the horizontal. These waves travel rapidly to the right with an approximate speed of 2.2 m s^{−1} and oscillate in time with an approximate frequency of 10 Hz. We compare this wave speed and frequency with that derived for the clinical data in §6. Second, these waves appear quite abruptly (near *t*=2000 ms) and show a high degree of spatial organization, similar to the clinical seizure data discussed in §2. Careful examination reveals that the oscillations appear to be triggered near 200 mm<*x*<300 mm, and subsequently spread across the whole domain.

Waves in continuum models of the human cortex are not new. Steyn-Ross *et al*. (2003) investigate the *dimensional* SPDEs and show that spatial modes develop when the corticortical e→i diffusivity dominates the corticortical e→e diffusivity. Jirsa & Kelso (2000) show that spatio-temporal patterns (reminiscent of the seizing cortex) occur in a continuum model of cortical electrical activity with a single inhomogeneous connection. By adjusting the extent of this inhomogeneous connection, the stability of the spatio-temporal patterns changes. In addition, Pinto & Ermentrout (2001) show how travelling waves may develop in a model derived from the physiology of a neocortical slice.

In the second example, we consider numerical solutions to the dimensionless SPDEs for *P*_{ee} non-uniform in space and constant in time. Here, we set *P*_{ee} to be a Gaussian function in *x* with maximum at *x*=350 mm, full width at half maximum of 46 mm, and minimum *P*_{ee}=11.0. The localized region of hyperexcitation near *x*=350 mm may be thought of as the seizure focus in our simple model. We solve the dimensionless SPDEs with *Γ*_{e}=0.87×10^{−3}, *α*=0.001 (both uniform in space and time) and show the solution in figure 8. We find that localized oscillations (the dark ridges in figure 8) emerge from the region of hyperexcitation near *x*=350 mm, travel outward with approximate speed 1.2 m s^{−1} and approximate frequency 7.5 Hz, and then dissolve in the regions of lower excitation, where *P*_{ee}→11.0. This localized oscillatory activity is more representative of a localized seizure than the global oscillations illustrated in figure 7*b*.

In the final example, we compute numerical solutions to the dimensionless SPDEs for *P*_{ee} non-uniform in space and time. As in the previous example, we fix *Γ*_{e}=0.87×10^{−3} and *α*=0.001 (both uniform in space and time). However, unlike the previous example, we vary *P*_{ee} in time. At *t*={0, 5000} ms, we fix *P*_{ee}=11.0 uniform in space. For 500 ms<*t*<5000 ms, we set *P*_{ee} to be a Gaussian function in *x* with maximum at *x*=350 mm, full width at half maximum of 46 mm, and minimum *P*_{ee}→11.0. At *t*={500, 1000, 1500, 2000, 2500} ms, we increase in steps of 100 so that , respectively. These increases are marked by the solid vertical lines in figure 9*a*. We then decrease in steps of 100 so that for *t*={3000, 3500, 4000, 4500} ms, , respectively. These decreases are marked by the dashed lines in figure 9*b*. We show in figure 9 that localized oscillations in *h*_{e} (the dark ridges) begin at *t*=1900 ms, where . The frequency of these oscillations increases when we increase to 510 at *t*=2500 ms. Then, as we decrease between 3000 ms<*t*<5000 ms, the oscillations in *h*_{e} decrease in frequency and eventually disappear.

This simple example provides a crude model for the evolution of the seizing cortex. At normal levels of subcortical excitation (*P*_{ee}=11.0) there exist only disorganized spatio-temporal fluctuations in *h*_{e} and no pathological oscillatory activity. Only when exceeds a threshold level do localized oscillations in *h*_{e} appear. As the subcortical excitation is increased further, the frequency of the oscillations increases. Then, as is decreased, the oscillations slow and eventually disappear.

## 6. Results

In §2, we showed ECoG data recorded at two neighbouring electrodes during seizure, and in §§3–5, we computed numerical solutions to a continuum model of cortical electrical activity. Here, we compare the observed ECoG data with the model solutions for *h*_{e}. Specifically, we compute two quantities: (i) the average frequency at which the power spectrum achieves its maximum during seizure, and (ii) the average speed at which the voltage peaks propagate between two subdural electrodes during seizure. We denote these quantities *f*_{0} (the peak frequency) and *v* (the propagation speed), respectively. We show that *f*_{0} and *v*, computed from the model solutions and the ECoG seizure data, are roughly in agreement.

In our determination of *f*_{0}, we consider data collected at the same two subdural electrodes (e.g. *X* and *Y*) for each of six seizures from the same patient. Part of one such seizure is shown in figure 1. To compute *f*_{0}, we first low pass filter the ECoG data *X* and *Y* below 55 Hz, and then compute the WPS. For example, in figure 2*a* we show the WPS for the time-series *X* recorded during a typical seizure. In this case, the maximum power occurs between 1 and 10 Hz following seizure onset (at *t*≈17.5 s). To determine *f*_{0}, we calculate the average frequency of maximum power in the WPS for the duration of the seizure (here, for 17.5 s<*t*<50 s). We list *f*_{0} and its standard deviation over the duration of the seizure for both electrodes (i.e. *X* and *Y*) and each seizure in table 3. For reference, we label the seizures *S*_{1} to *S*_{6}, where *S*_{1} corresponds to the data discussed in §2. We also compute the average of *f*_{0} over the six seizures for *X* and *Y* and find 4.1±0.8 and 5±1 Hz, respectively.

Propagating waves of electrical activity have been observed in numerous mammalian systems and during seizures in rats and in cats (Chervin *et al*. 1988; Connors & Amitai 1993; Pinto & Ermentrout 2001). Here, we assume that, during each seizure, voltage peaks propagate between the two neighbouring subdural electrodes (with time-series *X* and *Y*). That this assumption is valid can be inferred from the data shown in figure 1*b*. We note for 17.5 s<*t*<27 s the similarity of the wave forms in *X* and *Y* and that peaks in *X* precede peaks in *Y*. To determine the average propagation speed *v* of these voltage peaks, we first low pass filter the time-series below 55 Hz. We then compute the WCC between *X* and *Y*, as shown in figure 3 for *S*_{1}. At each time in the WCC, a maximum correlation occurs for some lag time. For example, at *t*=20 s, the maximum correlation occurs at a lag of approximately 25 ms. We compute this lag at which the maximum correlation occurs for the duration of the seizure. We find that, in general, the lag values are consistent over two temporal intervals. The first interval (*I*1) includes the time of seizure onset and the 10 s that follow. The second interval (*I*2) includes all later times (i.e. all times 10 s after the seizure onset until seizure termination). For example, in *S*_{1} the two time-intervals over which the lag of maximum correlation is consistent are: *I*1=17.5 s<*t*<27 s and *I*2=28 s<*t*<50 s. Such intervals can be chosen for each seizure except *S*_{6}. For this seizure, the lag values of maximum correlation are near zero for all times following seizure onset. Therefore, we do not include *S*_{6} in our analysis of *v*. We determine *v* in each interval by dividing the known electrode separation (10 mm) by the average lag in each interval. We list *v* and its standard deviation in each interval and for each seizure in table 4. We note that for *S*_{3} the average lag is 0 ms in *I*2 and therefore *v* is infinite. We compute the average and standard deviations of *v* in *I*1 (over *S*_{1} to *S*_{5}) and *I*2 (over *S*_{1}, *S*_{2}, *S*_{4} and *S*_{5}) to find 0.5±0.1 and −3±1 m s^{−1}, respectively. We note that travelling waves of electrical activity on animal cortices have been reported with speeds of 0.06–0.09 m s^{−1} (Chervin *et al*. 1988).

To compare these observational results with the model results, we compute *f*_{0} and *v* from the numerical solutions to the model cortex for the variable *h*_{e}. We tabulate the results in table 5 for the dimensionless ODEs simulation shown in figure 5*a*, and the dimensionless SPDEs simulations shown in figures 7*b* and 8. We note that the results of the simulations agree with the observed data for *f*_{0} (within a factor of 2) and *v* in *I*1 (within a factor of 5).

## 7. Discussion

Our aim was to consider the potential of the SPDEs to provide a model of the mesoscopic electrical activity on the seizing cortex. In our analysis of the dimensionless ODEs, we have shown that there do exist limit cycles in *h*_{e} associated with Hopf bifurcations in the dynamics. We have also shown that the simpler dimensionless ODEs suggest behaviour in the complicated dimensionless SPDEs that is consistent with the clinical seizure data. Specifically, we have shown that at typical parameter values *Γ*_{e}=1.42×10^{−3} and *P*_{ee}=11.0, a reduction in *Γ*_{e} (to *Γ*_{e}=1.20×10^{−3}) produced oscillatory behaviour in *h*_{e}. However, this oscillatory behaviour was transient and did not produce the large-amplitude, stable oscillations characteristic of the seizing cortex. Moreover, this oscillatory behaviour occurred at a single value of *Γ*_{e}. Thus, the parameter values of the model cortex would have to be carefully tuned to reach this point.

To increase the excitability of the model cortex and make it ‘seize’, we increased *P*_{ee}. We show in figures 5*a* and 6 that for *P*_{ee}>250.0 there exist stable, large-amplitude limit cycles in *h*_{e} over a wide range of *Γ*_{e}. We suggest in §5 that these oscillations are associated with travelling wave dynamics in the dimensionless SPDEs. Thus to produce behaviour in *h*_{e} consistent with that observed in the seizing cortex, we must: (i) increase *P*_{ee}, and (ii) decrease *Γ*_{e} in the dimensionless SPDEs. In dimensional terms, we must increase the subcortical excitation *p*_{ee} or decrease the maximum firing rate *S*^{max}, and we must increase the EPSP neurotransmitter rate constant *γ*_{e} or the difference between the reversal and resting potentials or decrease the peak excitatory postsynaptic potential *G*_{e} or *S*^{max}.

That we must decrease *Γ*_{e} to induce a seizure seems counterintuitive. We expect that an increase in the strength of excitatory inputs would promote seizing activity. However, the decrease in *Γ*_{e}, from 1.42×10^{−3} to 0.87×10^{−3} used to create figure 9, represents only a 39% change. The change in *P*_{ee} required to induce a seizure is much larger. We show in figure 9 that this parameter must increase by at least 2700%. We may interpret the changes in these two parameters required for seizure induction following the models in da Silva *et al*. (2003). Alterations of the epileptic brain, owing to genetic or environmental (e.g. injury) effects, for example, predispose it to seizures. In modelling the epileptic cortex, we may interpret this predisposition as a permanent (model I) or a gradual (model II or model III) change of a model parameter. The results shown in figure 9 illustrate the latter case; with *Γ*_{e} fixed, a gradual increase in *P*_{ee} induces seizing activity in *h*_{e}. Because the change in *P*_{ee} is gradual and affects the dynamics of *h*_{e} we could try to anticipate these types of seizures. We may also interpret these results following the model I framework. For example, consider a permanent increase in *P*_{ee} from 11.0 to 310.0. When modelling the cortical activity of any individual (either healthy or epileptic), we may assume that the model parameters undergo routine fluctuations in their normal values. In healthy individuals (with the typical parameter value *P*_{ee}=11.0), a small negative fluctuation in the model parameter *Γ*_{e} may induce a harmless, transient oscillation (like that shown in figure 4) but not a seizure. This same fluctuation in *Γ*_{e} will cause the predisposed epileptic cortex to seize. In this case, the increased value of *P*_{ee} allows the dynamics to cross the separatrix (or boundary) between the normal and seizing states (Robinson *et al*. 2002). Thus, only a small change in *Γ*_{e} is necessary to produce pathological oscillations in *h*_{e}. Because a random fluctuation in *Γ*_{e} induces the seizure, these types of seizures cannot be anticipated. From the results presented here, we suggest a large change in *P*_{ee} and a small change in *Γ*_{e} may provide a model of the seizing cortex.

In tables 3–5 we compute two results, the average peak frequency *f*_{0} and the average propagation speed *v*, derived from clinical ECoG recordings and solutions to the SPDEs. We find that the values derived from the observational and model data approximately agree (*f*_{0} within a factor of 2 and *v* in *I*1 within a factor of 5), but note the following important issues. First, the subdural electrode strips used in the clinical recordings consist of six electrodes each separated by 10 mm. Thus, the spatial sampling of the cortex is poor and we cannot observe detailed wave behaviour in the voltage. It may be true that electrodes *X* and *Y* lie on either side of a sulcus (a groove or furrow in the cortex); the precise locations are not known. It may also be true that a wave in electrical potential is propagating obliquely with respect to the shortest path between electrodes *X* and *Y*. Future recordings that utilize high-density electrode grids with small interelectrode spacing would provide better results to compare with the theory. Second, the lead/lag relationships determined from the WCC are ambiguous. For example, we showed in figure 3 that for 17.5 s<*t*<27 s, the oscillations at electrode *X* lead the oscillations at electrode *Y* by approximately 25 ms. During this interval, these time-series are nearly sinusoidal with a well-defined frequency near 4 Hz; see the data in figure 1*b* or the WPS in figure 2. The cross-correlation of two sinusoids is itself a sinusoid, and therefore ambiguous; although we suggest *X* leads *Y* by approximately 25 ms, it is also true that *Y* leads *X* by approximately 225 ms (the period of the 4 Hz cycle minus 25 ms). We chose the case *X* leads *Y* for two reasons: (i) the amplitude of oscillations in *X* is bigger than in *Y*, and (ii) a correlation at 25 ms between two closely spaced electrodes is more reasonable than a correlation of 225 ms. It does not seem possible to resolve this ambiguity of the cross-correlation. Third, we have assumed that the same voltage wavefront passes through both electrodes. If the wave source lies between the two electrodes, then this assumption is incorrect. Smaller electrode spacing may help locate the wave source more accurately. Fourth, given such a complicated model (14 nonlinear first-order stochastic differential equations and 20 parameters), we could potentially produce any desired behaviour in *h*_{e} by properly tuning the parameters. Here, we chose to adjust two parameters, both related to the excitation of the model cortex. In this way, we used the physiology of the seizing cortex to constrain the adjustable parameters in the model. Finally, we do not mean to suggest that the SPDEs provide an accurate model of all human cortical seizure activity. Here, we have only compared the model results with data collected from a single subject, and shown that the model and observational data agree in an approximate way.

Although confounding factors exist both in the model and observations, these results still illustrate an important use of the SPDEs and other continuum models: establishing links between the known ECoG data and unknown cortical physiology. Traditionally, neurophysiologists visually analyse ECoG data recorded from a seizing patient. No attempt is made to link the changes observed in these time-series with corresponding changes in cortical physiology. Continuum models of cortical electrical activity allow one to consider such links in a quantitative way. This is especially important for human cortical data, where invasive procedures are inappropriate. For example, Robinson *et al*. (2002) consider an ODE model of the corticothalamic network. Unlike the SPDEs model presented here, these models explicitly include the interactions of cortical neurons with the thalamus. The authors show that the model can produce the 3 Hz spike wave characteristic of absence seizures and how the model parameters relate to the period of the seizures. A corticothalamic model is used to study other types of epilepsy in da Silva *et al*. (2003). In this work, we relate the emergence of pathological oscillations in the measured cortical voltage to changes in two parameters: a small decrease in *Γ*_{e} and a large increase in *P*_{ee}. In this way, we suggest a link between the ECoG data and the physiology of the cortex.

Measures of dynamical quantities may also be computed from observed and simulated time-series and compared. For example, the correlation dimension, largest Lyapunov exponent (LLE), and synchronization have been shown in observations to decrease preceding a seizure (Quyen *et al*. 2003; Iasemidis 2003). Such quantities are calculable from the model equations. Dafilis *et al*. (2001) vary two *dimensional* parameters (*p*_{ee} and *p*_{ei}) and compute the LLE and Kaplan–Yorke dimension for the *dimensional* ODEs model from Liley *et al*. (1999). A similar analysis may be performed in which one varies the parameters *P*_{ee} and *Γ*_{e} of the dimensionless ODEs model presented here. Although the dimension and LLE determined from noisy ECoG data are of limited accuracy, a comparison of the theoretical and observational results may provide verification or suggest improvements to the model (Freeman 1992).

Understanding both the ECoG patterns of cortical ictal activity and the related changes in cortical physiology may help shape new strategies for seizure detection and focus location. This, in turn, may allow clinicians to deliver abortive seizure therapy that is targeted both in time and in space. For example, experimental methods such as focal cooling, brief pulse cortical electrical stimulation and local drug delivery all depend on reliable algorithms for the detection and prediction of seizures (Eder *et al*. 1997; Lesser *et al*. 1999; Yang *et al*. 2002). The combination of parameters one can change in a well-motivated, physiologically relevant model to simulate these dynamics may suggest strategies for therapeutic intervention.

## Acknowledgments

M.A.K. would like to acknowledge the support of a Graduate Fellowship from the U.S. National Science Foundation, and a Berkeley Fellowship from the University of California. The authors are grateful to Nicholas M. Barbaro, MD, for electrode placement and surgical management.

## Footnotes

- Received November 23, 2004.
- Accepted January 29, 2005.

- © 2005 The Royal Society