The ictal wavefront is the spatiotemporal source of discharges during spontaneous human seizures

Abstract

The extensive distribution and simultaneous termination of seizures across cortical areas has led to the hypothesis that seizures are caused by large-scale coordinated networks spanning these areas. This view, however, is difficult to reconcile with most proposed mechanisms of seizure spread and termination, which operate on a cellular scale. We hypothesize that seizures evolve into self-organized structures wherein a small seizing territory projects high-intensity electrical signals over a broad cortical area. Here we investigate human seizures on both small and large electrophysiological scales. We show that the migrating edge of the seizing territory is the source of travelling waves of synaptic activity into adjacent cortical areas. As the seizure progresses, slow dynamics in induced activity from these waves indicate a weakening and eventual failure of their source. These observations support a parsimonious theory for how large-scale evolution and termination of seizures are driven from a small, migrating cortical area.

ARTICLE

Received 16 Sep 2015 | Accepted 19 Feb 2016 | Published 29 Mar 2016

Elliot H. Smith1, Jyun-you Liou2, Tyler S. Davis3, Edward M. Merricks4, Spencer S. Kellis5, Shennan A. Weiss6, Bradley Greger7, Paul A. House3, Guy M. McKhann II1, Robert R. Goodman8, Ronald G. Emerson9,Lisa M. Bateman10, Andrew J. Trevelyan4,* & Catherine A. Schevon10,*

The extensive distribution and simultaneous termination of seizures across cortical areas has led to the hypothesis that seizures are caused by large-scale coordinated networks spanning these areas. This view, however, is difcult to reconcile with most proposed mechanisms of seizure spread and termination, which operate on a cellular scale. We hypothesize that seizures evolve into self-organized structures wherein a small seizing territory projects high-intensity electrical signals over a broad cortical area. Here we investigate human seizures on both small and large electrophysiological scales. We show that the migrating edge of the seizing territory is the source of travelling waves of synaptic activity into adjacent cortical areas. As the seizure progresses, slow dynamics in induced activity from these waves indicate a weakening and eventual failure of their source. These observations support a parsimonious theory for how large-scale evolution and termination of seizures are driven from a small, migrating cortical area.

1 Department of Neurological Surgery, Columbia University Medical Center, New York, New York 10032, USA. 2 Department of Physiology and Cellular Biophysics, Columbia University, New York, New York 10032, USA. 3 Department of Neurosurgery, University of Utah, Salt Lake City, Utah 84132, USA.

4 Institute of Neuroscience, Newcastle University, Newcastle upon Tyne NE1 7RU, UK. 5 California Institute of Technology, Division of Biology and Biological Engineering, Pasadena, California 91125, USA. 6 Department of Neurology, UCLA David Geffen School of Medicine, Los Angeles, California 90095, USA.

7 School of Biological and Health Systems Engineering, Arizona State University, Tempe, Arizona 85287, USA. 8 Department of Neurosurgery, Icahn School of Medicine at Mount Sinai, New York, New York 10029, USA. 9 Department of Neurology, Weill Cornell Medical College, New York, New York 10021, USA.

10 Department of Neurology, Columbia University Medical Center, New York, New York 10032, USA. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to C.A.S. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

DOI: 10.1038/ncomms11098 OPEN

The ictal wavefront is the spatiotemporal source of discharges during spontaneous human seizures

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

Despite over 50 years of research into seizure electro-physiology, there remains a limited understanding of the evolving spatial and temporal structure of seizures. This

knowledge gap undoubtedly contributes to the enduring mystery surrounding many open questions in epilepsy research pertaining to seizure origin, evolution and termination. Clinical electroencephalography (EEG) recordings of focal seizures consistently show a progressive expansion of synchronized pathological activity to large cortical regions, followed by abrupt self-termination that occurs simultaneously across all regions1,2. The question of how the synchronized pathological activity is coordinated over expansive cortical areas remains unanswered, yet is a fundamental issue for both basic research and clinical studies that rely on the use of EEG to identify sites of seizure origin.

It is currently presumed that seizure activity in humans is generated from a cohesive cortical network spanning potentially several centimetres26; however, the linking mechanism for nodes in these networks is yet undetermined. Animal studies, on the other hand, have catalogued a range of spatially focused cellular or ionic changes during epileptic events711. The challenge is to provide a coherent explanation that encompasses the discrepancy in scale between these animal models of seizures and the clinical observations12.

In our previous work, we described how extensive cortical local eld potentials (LFPs) can be driven from a small cortical area exhibiting strong neuronal population ring13,14. We therefore proposed a spatiotemporal structure for seizures that corresponded with observations from both animal models1517 and spontaneous human seizures13,18. In this structure, the seizing territory is led by a slowly advancing, sharply demarcated, narrow (o2 mm) band of continuous (tonic) multiunit ring, termed the ictal wavefront13,19.

Behind this wavefront, in the seizing territory, there are synchronized rhythmic discharges that give rise to the classic EEG signature of seizures over broad areas of cortex. While the ictal wavefront is difcult to identify in standard or intracranial EEG recordings, the synchronized activity in the seizing territory can be detected from high-frequency LFP14.

Here we conduct a detailed investigation of the spatiotemporal dynamics of the seizing territory in the human brain, to characterize mechanisms of seizure evolution and termination. We hypothesize that, as previously reported in a murine slice model15, fast-moving travelling waves originating from the ictal wavefront trigger intense depolarizations and associated multiunit ring in the seizing territory. That these ictal discharges arise from the ictal wavefront indicates that the ictal wavefront is the focal point for the mechanisms of seizure origin, propagation and termination. Specically, we hypothesize that gradual reduction in the excitatory currents generated from a slowly weakening ictal wavefront is sufcient to explain the evolution of both eld potentials and multiunit activity (MUA) in the seizing territory, including the seizures eventual spontaneous termination.

ResultsNatural epochs in spontaneous human seizures. Recordings of clinical seizures were analysed in ve patients implanted with microelectrode arrays (MEAs), as previously described13,20,21. These were 4 4 mm Utah-style arrays, with 1-mm long

electrodes designed to record unit activity and LFPs from adjacent tissue in neocortical layers 4 and 5, and 400 mM spacing between the electrodes13,22. In three of these patients (patients 3, 4 and 5; see Supplementary Table 2; six seizures), the MEA was located within the cortical region invaded by the seizure. Following terminology we introduced previously13, we

term this region the ictal core. In these three patients, the MEA captured the expansion of the ictal core after seizure initiation, which was not captured directly. While all of the MEAs were implanted into the clinically dened seizure onset zone, in two patients, the seizures did not invade the sites of MEA implantation (patients 1 and 2, 10 seizures); that is, the MEAs were in penumbral sites, which remained outside the ictal core for the entire duration of each seizure (Supplementary Fig. 1). As previously reported, these penumbral areas were characterized by an absence of intense, synchronized MUA despite exhibiting low-frequency LFP that was otherwise typical of seizure activity13.

To address our hypotheses regarding the seizure generating territory, this article focuses on the three patients in which the MEAs were recruited into the ictal core. For these patients, seizure invasion at the MEA site was identied by the abrupt appearance of transient, intense and asynchronous multiunit ring referred to as the ictal wavefront, marking the moment of failure of strong inhibitory restraint13,19,23. This period of intense ring was followed by a rhythmic train of ictal discharges, that is, high-amplitude eld potential deections associated with transient bursts of multiunit ring and high-g activity, alternating with periods of relative silence between discharges (Fig. 1b,c). The dominant rhythm of ictal discharges gradually decreased in frequency towards the end of the seizure (Fig. 1d). On the basis of these observations, the MEA data from these seizures were subdivided into the following readily identiable phases (Fig. 1a): (1) a pre-recruitment epoch, in which ictal discharges may begin to appear in the low-frequency LFP;(2) the steady multiunit ring characterizing the ictal wavefront and recruitment into the seizure; (3) a post-recruitment epoch, following passage of the ictal wavefront, in which the inter-discharge intervals (IDIs) remain stable; and (4) a pre-termination epoch (means.d. duration 16.808.35 s) in

which the pace of the ictal discharge rhythm gradually slows and becomes irregular, until the seizure abruptly terminates (Fig. 1; Supplementary Fig. 2a).

Such a transition from regular ictal discharging to slowing, irregular ictal discharges has been observed as a signature of approaching seizure termination for decades24. Here we quantied when this signature appeared (Methods), and operationally dened said time as the division between post-recruitment and pre-termination epochs. The range of division times was consistently small across MEA channels (s.d. 3.72 s) and electrocorticography (ECoG) channels for each

seizure (s.d. 6.07 s). Furthermore, IDI distributions were not

signicantly different between seizures recorded on MEA and ECoG (two-sample KolmogorovSmirnov tests, all PZ0.52;

Supplementary Fig. 2b). There was no signicant difference in the post-recruitment/pre-termination division times determined from MEA recordings compared with those determined from ECoG (MannWhitney U-tests, both PZ0.3; Supplementary

Fig. 2c), indicating that this phase transition occurs simultaneously in all tested channels.

Ictal discharges form travelling waves. Previous work has shown that travelling waves of synaptic potentials propagate across the brain under normal sensory processing25,26 (reviewed in ref. 27). Travelling waves have also been demonstrated by assessing site-to-site delays in high-amplitude peaks of epileptiform discharges2830. Travelling waves appeared in this study during ictal discharges: high-amplitude LFP events with duration up to 200 ms, moving rapidly across the cortical surface. In contrast, the slowly propagating ictal wavefront is distinguished from travelling waves by its composition of steady ring lasting several seconds with minimal effect on low-frequency activity13.

2 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098 ARTICLE

a

MEA

LFP

MUA

Prerecruitment

Ictal wavefront

Postrecruitment

Pre-termination

b

c

Ictalwavefront Termination

Prerecruitment Pre-termination

Postrecruitment

500 V

20 APs/s

d

Electrodes

e

IDI (s)

0.5

0.25

0 0 10 20 30 40 50 60 70

Time relative to seizure onset (s)

Figure 1 | Progression of seizure activity recorded from a microelectrode array in the ictal core. Colours indicate the different epochs in the seizure: light blue is the pre-recruitment epoch, orange is the ictal wavefront epoch, purple is the post-recruitment epoch and pink is the pre-termination epoch. Seizure termination is labelled with a dark blue dotted line. (a) Cartoon of hypothesized spatial organization of seizure activity relative to the microelectrode array (MEA) location. (b) Raw LFP traces recorded from a single microelectrode during a seizure. (c) Averaged ring rate over electrodes on the array during the seizure in b. (d) Multiunit raster plot over MEA channels during the seizure in b. (e) Inter-discharge intervals (IDIs) for each discharge through the duration of the same seizure. The abrupt transition from regular to irregular IDIs at 36 s after seizure onset marks the transition from post-recruitment (period following passage of the ictal wavefront) to the pre-termination epoch. In both b and c, the abrupt cessation of dischargesin both LFP and MUA, which is the dening event marking seizure termination, is evident.

We examined ictal discharges individually in the MEA recordings, to assess their spatiotemporal evolution on the ne (400 mm) scale afforded by the MEA conguration. Analysis of low-frequency (o50 Hz) LFP during ictal discharges consistently revealed variable latencies across the MEA (Fig. 2a). The direction of travel of each discharge was determined by ranking the minima of the rst derivative (that is, the maximally decreasing slope) of the low-frequency LFP in time. Using these ranked minima, a velocity vector was dened for each discharge, indicating the direction taken by the travelling wave (Fig. 2b). The extent of the measured delays across the MEA, which along its maximum diagonal spans 5.67 mm, was 25.9617.03 ms (means.d.; range 4135 ms; Fig. 2c). Travelling wave speeds estimated

from these delays were 0.210.31 m s 1 (means.d.; range 0.421.42 m s 1). Corresponding latencies were noted

in the associated MUA bursts (Supplementary Fig. 6a,b). These latencies determined from MUA were not signicantly different from those determined from low-frequency LFP (Supplementary Fig. 6c,d; repeated measures analysis of variance, P 1).

Examining travelling wave speeds in each epoch revealed that speeds increased signicantly from the pre-recruitment to post-recruitment and pre-termination epochs (means.d. speed for pre-recruitment: 0.0860.07 m s 1, post-recruitment:0.260.24 m s 1, pre-termination: 0.310.17 m s 1; analysis of variance with Tukey post hoc tests, Po0.01; Fig. 2c), as measured by both the magnitude of the velocity vector (Tukey post hoc test, Po0.01; Supplementary Fig. 3) and travelling wave delays across the array (Tukey post hoc test, Po0.01; Supplementary Fig. 4). This speed change is consistent with the collapse of inhibition

at the ictal wavefront. After recruitment, travelling wave speed remained fairly constant, with no signicant difference between the post-recruitment and pre-termination epochs (Tukey post hoc test, P 0.31).

Travelling waves arise from the ictal wavefront. Previous studies from mouse neocortical slices in vitro showed that ictal discharges propagate away from the ictal wavefront along the cortical space in a brain slice15. The hypothesis follows that an analogous propagation pattern of fast-moving synaptic potentials arises from the ictal wavefront, as it slowly and radially expands across the two-dimensional cortical sheet (see Supplementary Movie 2 in Schevon et al.13 for a clear example of ictal wavefront propagation across the MEA). To test this hypothesis, we rst determined the direction of propagation of the ictal wavefront from the propagation of asynchronous multiunit ring across the MEA. We then examined the directions of rapid travelling waves, relative to the direction of the slow ictal wavefront, for each discharge during the pre- and post-recruitment epochs and the pre-termination epoch.

Travelling waves in each epoch showed statistically signicantly preferred directions across the spatially xed layouts of each MEA and ECoG grid (Rayleighs test of non-uniformity, all Po0.01). Two predominant, opposing directions relative to the propagation direction of the ictal wavefront were observed across the MEA (Fig. 3a). Travelling wave directions before and after passage of the ictal wavefront through the MEA were signicantly different both within (Kuipers test, Po0.01) and across (Kuipers test, Po0.01) patients. Pre-recruitment travelling waves were

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

a b

Pre-recruitment

Post-recruitment

4 mm

Electrodes

4 mm

0.5 mV

10 ms

Electrodes

4 mm

Electrodes

4 mm

1 mV

10 ms

Electrodes

4 mm

Pre-termination

Electrodes

4 mm

1 mV

10 ms

Electrodes

c

80 d

Wave speed (m s1 )

*n.s.

80

80

0.3

*

Discharge count

60

60

60

40

40

40

0.2

20

20

20

0.1

0 0 50 100 150

0 0 50 100 150

0 0 50 100 150

Delay (ms)

Delay (ms)

Delay (ms)

0

Figure 2 | Ictal discharges form travelling waves across the microelectrode array. (a) Example, low-frequency LFP recorded from each microelectrode during three ictal (EEG) discharges, one from each seizure epoch as indicated by the coloured asterisks in the seizure traces, colour coded by when the maximally negative slope of the travelling wave occurs on each microelectrode. Seizures from two patients (patient 3 and 5) are shown. (b) The footprint of the microelectrode array corresponding to the three ictal discharges in a, with electrode positions colour coded the same way as the LFP in a. Vectors indicating travelling wave direction are superimposed on the microelectrode array footprint in white. (c) Histograms of delays between the rst and last discharges on array microelectrodes during each epoch (pre-recruitment in light blue, post-recruitment in purple and pre-termination in pink).(d) Measures of discharge speed across all patients, during the three epochs shown in a. Post hoc Tukeys range test determined that pre-recruitment bursts (N 273) were signicantly slower than both post-recruitment (N 315) and pre-termination (N 282) bursts (P 9.5 10 10 and P 9.5 10 10).

Post-recruitment and pre-termination bursts did not exhibit signicantly different speeds (P 0.82). Error bars indicate s.e.

aligned with the direction of the slow-moving ictal wavefront (one-sample test for mean angle, all P40.05; Fig. 3b).

Post-recruitment and pre-termination travelling wave directions were oriented in the opposite direction of the slow movement of the ictal wavefront (one-sample test for mean angle, Po0.01 for travelling waves in both post-recruitment and pre-termination epochs; Fig. 3c,d).

This perceived ip in travelling wave direction from the vantage point of the xed cortical recording site of the MEA, coinciding with ictal wavefront passage, indicates that synaptic travelling waves during ictal discharges originate from the ictal wavefront in space and time. Specically, before the arrival of the ictal wavefront at the MEA, discharges travelled in the same direction as the ictal wavefront. After ictal wavefront passage, the discharges abruptly changed direction, that is, were directed back into the ictal core (Supplementary Movie 2).

Travelling waves during ictal discharges were also observed extending across the ECoG (subdural) grid. Figure 4a shows examples of ictal discharges recorded across ECoG electrodes in two patients, colour coded by when the maximally decreasing slope occurs, and the ranked delays across ECoG electrodes. These travelling waves across ECoG electrodes also showed consistent directions within each seizure (Rayleighs test of nonuniformity, all Po0.01), and matched the direction of travelling waves seen in the simultaneous MEA recording. Figure 4b shows the median ranks of travelling wave delays superimposed on brain reconstructions for patients 3 and 5. Travelling waves on ECoG appeared to extend from the ictal wavefront as well, having overlapping relative direction distributions with pre-recruitment waves on the MEAs (Fig. 4c; Kuipers test, all PZ0.1). Low-frequency LFP was also increasingly synchronized across the ECoG grid as seizures approached termination (means.d. slope

4 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098 ARTICLE

a

90

270

b

90

270

Vector length

Traveling wave direction (degrees)

2

1.5

180

0

180

Ictal wavefront propagation direction

Ictal wavefront propagation direction

Prerecruitment

1

Postrecruitment

Pre-recruitment

Traveling wave direction (degrees)

c

90

270

d

90

270

180

180

Post-recruitment

Normalized discharge count

Pre-termination

Figure 3 | Travelling wave directions during seizure epochs are dependent on the location of the ictal wavefront. (a) Example velocity vectors for one seizure before and after recruitment. The black line indicates the direction of ictal wavefront passage across the MEA. The shift in direction between pre-recruitment and post-recruitment discharges that occurs precisely at the time of ictal wavefront passage is apparent. In each case, the ictal wavefronts location relative to the MEA site is opposite to the predominant travelling wave direction. (bd) Distributions of travelling wave directions across all seizures, relative to the direction of the ictal wavefront (black line), coloured by epoch (b: pre-recruitment, c: post-recruitment, d: pre-termination).

of linear regression: 0.0110.0056; Fig. 4d). This pre-termination synchrony increase was apparent in 15 out of 17 seizures examined (Supplementary Fig. 4). These results show that rapidly propagating travelling waves from the ictal wavefront are also apparent across ECoG electrodes, where they give rise to the commonly observed hypersynchrony during seizures.

Simulated ictal wavefront reproduces termination seizure dynamics. That the speed and direction of travelling waves during ictal discharges are dependent on the relative location of the ictal wavefront suggests the ictal wavefront is the agent of both seizure propagation and termination. The hypothesis follows that a dissipation, weakening or cessation of the ictal wavefront should result in seizure termination. To acquire a theoretical understanding of how input from the ictal wavefront affects the neural dynamics of the ictal core, we examined a simple mean-eld computational model of the ictal core. Specically, we modelled the contribution of the ictal wavefront as input to a small recurrently connected network of simulated neurons. The model used experimentally determined parameters of layer 5 pyramidal neurons. The specics of the model, including these parameters, are described in detail in Supplementary Information. Similar approaches have been used as simple models of spinal cord neuron dynamics and working memory traces31,32.

The qualities of the MUA dynamics we observed in the ictal core (Fig. 5a,b) were reproduced well in our model by simply decreasing the external input to the simulated network, as would happen with a gradual dissipation of the ictal wavefronts excitatory drive. Specically, with high external input, a seizure-

like pattern was initiated and the model was locked into a tonically ring state. As input decreased, the model activity rst underwent a transition from tonic ring to a pattern of rhythmic discharging. As input continued to decrease, peak ring rate, IDI and discharge width all increased steadily until the seizure terminated with an abrupt cessation of discharging. These trends were present both in the model and in the recorded seizure data (Fig. 5ce).

Signatures of a weakening ictal wavefront. We next sought to quantify the predicted effects of our hypothesized seizure spread and termination model in multi-scale electrophysiological recordings of human seizures. Due to the effect of seizure recruitment on action potential waveforms33, multiunit events were detected by thresholding the ltered (5003000 Hz) MEA data as in prior publications13,14,33. To assess whether the predicted increase in multiunit ring rates towards seizure termination were present in the ictal core, we quantied trends in peak MUA and high-g amplitude during the pre-termination epoch using linear regression. Linear regression slopes were tested against a zero-slope null hypothesis to determine signicance of trends in activity leading up to seizure termination. In the ictal core, peak MUA averaged across all recorded microelectrode channels increased signicantly as the seizure progressed from recruitment to termination (Wilcoxon signed-rank test, Poo0.01;

that is, Fig. 6a). In contrast, no signicant trend in peak MUA was observed across MEA channels in the penumbra before seizure termination (Wilcoxon signed-rank test, P 0.41;

Supplementary Fig. 1c).

Our hypothesis predicts that synchronization of unit ring gradually decreases in the ictal core as a seizure approaches termination. This trend occurs despite the increase in multiunit ring rate. In the MEA data, there was a signicant decrease in the high-g envelope across MEA channels in the ictal core (Wilcoxon signed-rank test, Poo0.01; Fig. 6a). In contrast, high-g recorded at penumbral MEA sites exhibited a signicant increase in pre-termination high-g power (Wilcoxon signed-rank test, Poo0.01; Supplementary Fig. 1b). Peak high-g frequency also decreased in seizures recorded from the ictal core from 183.925.4 Hz at the beginning of the post-recruitment period to 141.516.1 Hz at seizure termination (Wilcoxon signed-rank test, Po0.01; Fig. 6b).

Although the MEA data are highly specic, they can provide information only from the 4 4-mm coverage area. The clinical

ECoG grid, on the other hand, samples a region spanning several centimetres, but with lower (one cm) spatial resolution. While MUA is not accessible from these large subdural electrodes, the location of the ictal core can be inferred from high-g activity recorded on ECoG, which acts as a surrogate for intense neuronal bursting14. Paralleling the MEA recording, we expected to nd decreases in high-g frequency and amplitude across clinical ECoG electrodes that exhibited high-g activity indicative of recruitment into the ictal core. We therefore tested all ECoG recording sites at which the high-g amplitude exceeded ve times the s.d. of pre-ictal recordings (15.915.0 electrodes per seizure). High-g amplitude decreases were observed in these ECoG electrodes in all ve patients, including the two patients with the MEA positioned outside the ictal core (Wilcoxon signed-rank test, Po0.01; Fig. 6c), with decreases noted individually in 15 of 17 seizures. High-g peak frequency also decreased in recordings from the same ECoG electrodes during the pre-termination period, from 138.529.9 to 109.232.1 Hz (Fig. 6d).

Pre-termination multiunit and high-c activity correlation. The progressive dissociation in the trends between high-g and ring

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

c

a b

Early

Late

1 mV

ECoG MEA

Direction of ictal wavefront propagation

40 ms

MEA

1 mV

MEA

d

1 mV

0.5

MEA

1

Variance explained

0.9

0.8

0.7

0.6

40 30 20 10 0 10Time relative to seizure termination (s)

Figure 4 | Travelling waves propagating across the microelectrode array correspond to those propagating across the ECoG grid. (a) Examples of ictal discharges recorded across ECoG electrodes in two patients colour coded by when the maximally decreasing slope occurs, with cyan indicating earlier electrodes and magenta later electrodes. Colour-coded circles above the discharges indicate a spatial representation of delays over the ECoG grid, with the same colour code as the voltage traces. (b) Reconstructions of ECoG grid and MEA locations for the same patients in a. Average delays across ECoG electrodes are colour coded as in a. Black arrows indicate MEA locations. (c) Distributions of MEA pre-recruitment travelling waves (black) and travelling waves ECoG (orange). (d) Three channels of LFP recorded from ECoG electrodes during one seizure (patient 3) and a representative increase in variance explained by the rst-principal component, across all ECoG electrodes on the grid shown in b (lower), during the pre-termination epoch. The black dotted line shows linear regression used to determine the slope of the variance explained during the pre-termination epoch for the example seizure.

rate, as a seizure approaches termination, indicates that the previously observed correlation between the two signals34,35 may have multiple contributing factors. It has been proposed that action potential synchrony provides a substantial contribution to the amplitude of the high-g signal3638. To address the relative contribution of action potential synchrony to high-g amplitude, correlations among peak multiunit ring rate, peak high-g amplitude and multiunit discharge width were examined among individual discharges during the pre-termination epoch. Examining these metrics for each discharge controls for the variable intervals between successive discharges and allows for measurement of synchrony over discrete time periods, dened by the natural progression of the seizure. Again, there was a signicant anti-correlation between per-discharge ring rate and high-g amplitude as would be expected from the aforementioned trend (Spearmans r 0.21, Po0.01; Fig. 7b). However, a

greater anti-correlation was observed between per-discharge high-g amplitude and the multiunit discharge width (Spearmans r 0.60, Poo0.01; Fig. 7c), which was dened

as the s.d. of a Gaussian t to the MUA during each discharge (Fig. 7a). This result suggests that the degree of temporal dispersion of MUA is correlated with the amplitude of the high-g signal recorded across the MEA.

To examine the relationship between multiunit synchrony and high-g amplitude explicitly, mutual information for multiunit event coincidence was measured among all pairs of micro-electrodes during the pre-termination epoch. Mutual information was signicantly correlated with high-g amplitude over discharges (Spearmans r 0.38, Poo0.01; Fig. 7d). Together, these

results show that action potential ring in the ictal core desynchronizes as a seizure progresses towards termination. Desynchronization therefore manifests as decreasing high-g amplitude and peak frequency in both MEA and ECoG recordings.

DiscussionThe analyses of spontaneous human seizures and simulations of ictal neurons presented here support the hypothesis that the migrating ictal wavefront is the primary source of ictal activity, and consequently also the agent of seizure termination (Fig. 8). Thus, the seizure is sustained by activity in newly recruited territories, much like a forest re is sustained by fresh wood. We describe a straightforward electrophysiological mechanism that is sufcient to explain the major electrophysiological hallmarks of populations of neurons in the ictal core leading up to seizure termination: increased multiunit ring rate, desynchronization of

6 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098 ARTICLE

a

a

Microelectrode

Pre-

recruitment

Ictal

wavefront

Post-

recruitment

Pre-

termination

100

High-gamma

amplitude (V)

(spikes per s)

1,600

1,400

80

60

Mean firing rate (spikes per s)

1,200

100

3

1,000

800

Firing rate

600

40

Input strength (nS)

400

20

200

Model Data

50

0

0 5 10 15 20 25 30 35 40 45

b

300

15,000

100

Frequency (Hz)

0

2

200

0

5

10

Time (s)

Power (V2 )

15

20

30

10,000

b

100

5,000

100

Firing rate (spikes per s)

80

0 0 5 10 15 20 25 30 35 40 45

60

Time (s)

40

c

150

0 0 5 10 15 20 25 30 35 40 45

ECoG

20

100

0 40

30

20

10

0

High-gamma

amplitude (V)

Time from seizure termination (s)

25

50

c d

Input strength (nS)

2.82.72.62.52.42.32.2

Input strength (nS) e

2.82.72.62.52.42.32.2

0

0

2.82.72.62.52.42.32.2

0

0 5 10 15 20 25 30 35 40 45

1

0.2

Input strength (nS)

d

200

400

Normalized firing rate

Inter-discharge

interval (s)

Discharge width

(seconds)

0.8

0.8

1 Model Data

150

300

0.4

0.4

Frequency (Hz)

0.1

100

Power (V2 )

0.6

0.6

200

0.2

Model

Data

0.2

Model

Data

50

0 10

5

0

0 10

5

0 10

5

Time relative to seizure termination (s)

Time (s)

Figure 5 | Decreasing input to a network of simulated neurons recreates mean-eld dynamics observed in the data. (a) Mean ring rate of the simulated ictal network progresses from tonic ring, to rhythmic discharging, to self-termination, with decreasing extra-network input. The colours correspond to the epochs labelled in Fig. 1 as apparent in the models ring rate. (b) Firing rate from all microelectrode channels during one spontaneous human seizure, showing the progression from tonic ring during ictal wavefront passage, to rhythmic discharging, to seizure termination. The coloured bars above the seizures correspond to the epochs labelled in Fig. 1. (c) Peak mean ring rate during each discharge in the pre-termination period for the simulated neurons (green) and the seizure shown in b (purple). (d) Inter-discharge intervals during the pre-termination period, for the simulated neurons (green) and the seizure shown in b (purple). (e) Discharge width during the pre-termination period, for the simulated neurons (green) and the seizure shown in b (purple).

Figure 6 | High-c activity in the seizure core decreases in amplitude and frequency before seizure termination. (a) High-g instantaneous amplitude and multiunit ring recorded from an example microelectrode through the duration of the same seizure shown in Fig. 5b. The blue dotted line indicates linear regression of the peak high-g amplitude for each discharge during the pre-termination epoch, and the orange dotted line similarly indicates peak ring rate. (b) Spectrogram of the same seizure in a. The white dotted line indicates peak high-g frequency through the duration of the seizure, and the blue dotted line indicates linear regression of the peak high-g frequency during the pre-termination epoch. c and d show the same seizure as in a and b recorded on the ECoG electrode adjacent to the microelectrode array. Dotted lines indicate the same linear regression measurementsas in a and b.

MUA, decreases in both high-g amplitude and frequency, changes in the directions of travelling waves of ictal synaptic activity, and a slowing rate and increasing width of ictal discharges. An important conclusion is that the wide-ranging effects of seizures appear to be triggered from a small cortical region, which migrates during the seizure. That is, the sustenance and eventual termination of the seizure may depend on the intense activity of small, spatially restricted populations of neurons. The simultaneous, wide-area seizure termination can therefore be parsimoniously explained by sufcient dissipation of input from the ictal wavefront.

Previous studies have shown increasing low-frequency synchronization towards the end of seizures3,5 and point to hypersynchrony as a primary electrophysiological hallmark of seizures7,39. On the basis of these observations, it has been proposed that seizures are mediated by a large-scale network structure with coordinated long-range neuronal interactions2,5,40. This model would suggest that the MUA associated with ictal discharges should also demonstrate increasingly tight synchronization over extended cortical territories. We found

that the opposite is the case. Further, the large-scale eld potential synchronization in ECoG recordings can be explained as a by-product of travelling waves moving rapidly across large areas of cortex (Fig. 4). We theorize that these travelling waves resemble those previously studied under normal conditions: both eld potential and multiunit ring components are generated locally in response to transmitted synaptic barrages25,27. The rapid travelling wave speeds, up to 1.42 m s 1, together with the relatively long (up to 500 ms) duration of ictal discharges29,41, result in substantial temporal overlap and thus high levels of measured synchronization. Furthermore, as IDIs increase leading up to seizure termination, there are gaps between discharges that are markedly more synchronous than the discharges themselves (Supplementary Fig. 5).

The discrepancy between the apparent increases in LFP synchrony and the decreases in MUA synchrony proved to be a consistent feature of seizures regardless of recording modality or scale. While multiunit ring bursts were aligned with ictal discharges, these bursts became gradually desynchronized as the seizure approached terminationdespite a corresponding

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

increase in ring rate. Although multiunit desynchronization could only be directly observed in the area sampled by the MEA, indirect evidence of the same process was noted across a much wider area: correlated with the decreasing multiunit synchrony was a decline in high-g amplitude and frequency that was present both in the microelectrode recordings and in the ECoG recorded from nearby clinical subdural electrodes. Accordingly, we observed that all clinical ECoG recordings in our study sample demonstrating signicant ictal high-g activity showed similar decreases in high-g amplitude and frequency, indicating that the desynchronization process observed in the MEA recordings occurs throughout the recruited cortex.

The local neuronal population dynamics of the ictal core, including the increasing IDIs leading up to seizure self-termination, can be explained as the effect of gradually decreasing excitatory input from the ictal wavefront. The computational model presented here indicates that a dissipating or weakening ictal wavefront is sufcient to produce all the observed temporal dynamics of the ictal core, including self-termination, across a wide cortical region. The ictal wavefront is o2 mm wide13,23 because its intense tonic ring can only be sustained by very-high-intensity excitatory input that results from the failure of strong feedforward inhibition. Furthermore, the model predicts that only a slight drop in excitatory drive is sufcient to end the tonic ring phase and for the local neuronal population to convert to a repetitive bursting pattern42,43. The model predicts that further decreases in input actually increase ring rate up until the point where bursting stops abruptly. All of these dynamics are present in MUA recorded from the ictal core in spontaneous human seizures.

Future experimental and modelling work will be needed to characterize the ictal wavefronts reduced input mechanism fully, as these passive recordings in epilepsy patients undergoing clinically necessitated monitoring do not permit destructive experimental manipulations. For example, while the current study provides evidence that paroxysmal discharges originate from the ictal wavefront, there is little data apart from computational studies42,43 addressing how neurons undergo a transition from intense tonic ring to patterns of rhythmic discharging. Although the current results support our hypothesis that input from the ictal wavefront weakens in the pre-termination epoch, the mechanisms of this weakening remain to be determined. For instance, the wavefront may run into a blind alley, surrounded by tissue that is already recruited and becoming increasingly refractory, or it may encroach onto other areas that have a superior inhibitory restraint. The detailed characterization of pre-termination seizure electrophysiology we have provided, however, can be used to direct experimental manipulations capable of testing these hypotheses.

Many proposed mechanisms for seizure termination in acute animal studies operate over a limited cortical region44, leaving open the question of how they could subserve large-scale, simultaneous seizure termination. These mechanisms include acidosis11, changes in intracellular and extracellular ion concentrations7,45, increased membrane potassium conductance46,47, suppression from subcortical structures48,49, and neurotransmitter50 or ATP depletion51,52. The thesis presented herein implies that much of observed ictal electrophysiological activity, both core and penumbral, occurs as a response to electrical signals broadcast from the ictal wavefront. While these results do not directly implicate a specic biophysical mechanism for human focal seizure termination, they provide fresh insights into the processes leading up to seizure termination in humans and explain how cellular-scale processes can be compatible with clinical EEG observations.

A salient nding presented here is the progressive dissociation between high-g and MUA leading up to seizure termination.

High-g activity, that is, high-frequency, broadband LFP ranging from B70 to 300 Hz, depending on the electrodes and ampliers used to record the electrical signals53,54, is functionally distinct from band-limited gamma (3050 Hz) oscillations and is thought to reect summated postsynaptic output from the local neuronal population36,55,56 in a submillimetre domain20,57. High-g power has been shown to be highly correlated with ring rate in studies of evoked visual stimuli in monkey visual cortex38, and in microwire recordings from the human medial temporal lobe35 or neocortically implanted microelectrodes and standard clinical electrodes14,5860. In contrast, our results indicate that changes in action potential synchrony, independent of ring rate, can affect the high-g signal.

However, further investigation of this effect in normal burst ring versus epileptic bursting is warranted.

It has been proposed that population action potential synchrony is an important contributing factor to high-g power36,37. Furthermore, simulations show that increased action potential synchrony increases LFP power, especially in the high-g range42,61. The current results support these mechanistic understandings of high-g by showing, in spontaneous human recordings, a decrease in multiunit synchrony that is highly correlated with a decrease in high-g amplitude and peak frequencydespite marked ring rate increases towards the end of a seizure. Both the increased ring rate and the decreased multiunit synchronization may be a consequence of the increased IDIs in the pre-termination phase. As discharges in a seizure separate temporally, there is a larger window for action potential ring to occur before the next discharge, or more time for fast metabolic and ionic recovery. In this scenario, more action potentials could re with greater temporal dispersion because there is a larger time window in which to do so. The strong relationship between MUA synchrony and the high-g amplitude suggests that action potential synchrony provides a substantial contribution to the high-g signal.

Finally, our results have implications for treatment and prevention of epilepsy. Therapies targeting the source of seizures have always been hampered by the lack of knowledge of what this is precisely. Our results suggest novel avenues for improving the detection and localization of sites that drive seizure activity using the speed and direction of travelling waves to localize the ictal wavefront. Recent developments in high-density recording technology have provided the means to record many single units across the surface of the human brain from a thin-lm substrate62. Such recordings could detect travelling wave directions, as well as the ictal wavefront and core, over broad cortical areas with high spatial resolution.

Future investigations into the mechanisms of ictal wavefront dissipation and seizure termination could yield new approaches for therapies to terminate seizures earlier, by hastening the dissipation of the ictal wavefront, for example. Stopping seizures earlier could greatly reduce the burden of epilepsy, as longer seizures have been associated with more severe postictal paralysis63, peri-ictal respiratory dysfunction64 and are likely associated with more severe postictal maladies in general.

The current results support an analogous model of human seizure propagation to that observed in brain slice experiments and examine the models implications for the mechanisms of seizure termination. We have shown that local desynchronization leading up to seizure termination is accompanied by a breakdown in the commonly observed correlation between action potential ring and high-g activity. Local desynchronization in the ictal core occurs despite increasing large-scale synchronization of

8 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098 ARTICLE

a

b

Early

Late

800

Max MI Max MI

Filtered LFP (300 3,000 Hz)

Early

Late

40

Filtered LFP (0.3 50 Hz)

400

1,000

0.03

High-gamma

amplitude (V) Multiunit

activity

(events per s)

Gaussian Fit

Filtered LFP (0.3 50 Hz)

0.01

1,000

100

20

200

600

500

50

0

0

0

0

20

200

400

500

40

400

50 60

1,000

100

600

80

800

200

20

1,500

150 100

10

68.62 68.64 68.66 68.68

96.4 96.45 96.5 96.55 96.6

0

0

Time from seizure onset (s)

0 50 100Time (ms)

0 50 100 Time (ms)

c d e

45

Early

30

Early

0.14

Early

Peak multiunit firing rate (APs/s)

Multiunit discharge width (ms)

Mutual Information (z-score)

40

25

0.12

35

0.1

30

20

25

0.08

15

20

0.06

15

10

0.04

10

5

5

0.02

Late

Late

Late

0 0 500 1,000 1,500

0 0 500 1,000 1,500

0 0 500 1,000 1,500

High-gammma amplitude (V)

High-gammma amplitude (V)

High-gammma amplitude (V)

Figure 7 | Local desynchronization in the ictal core precedes seizure termination. (a) Examples of ltered data for two sample discharges during the pre-termination epoch. Low-frequency ltered data (0.350 Hz) are shown in grey and ltered data for multiunit activity detection (5003000 Hz) are shown in teal. The threshold for multiunit event detection is indicated with a blue line. Bursts are indicated on the seizure in grey using black arrows.(b) Examples of high-g (orange), ring rate (blue) and Gaussian t (black dotted lines) for burst width superimposed on the multiunit activity raster plot (teal) for two bursts in one seizure: one at the beginning of the pre-termination epoch (early) and one at the end of the pre-termination epoch (late).

Maximum mutual information values for each discharge are shown in boxes. (c) Relationship between peak multiunit ring rate and peak high-g amplitude for each discharge. (d) Relationship between high-g amplitude and multiunit discharge width for each discharge. (e) Relationship between peak high-g amplitude and multiunit synchrony in each discharge. Colours in c, d, and e represent when each discharge occurred during the pre-termination epoch, with black dots occurring earlier and copper dots occurring later.

Ictal wavefront

Pre-recruitment Pre-termination

Post-recruitment

Unrecruited

a

b

Ictal wavefront

Ictal wavefront Ictal wavefront Dissipating

ictal wavefront

Apparent directions of rapidly propagating travelling waves, relative to the slow-moving ictal wavefront, coloured by seizure epoch

Figure 8 | Summary model of spatial dynamics of seizure evolution and termination. (a) Schematic showing snapshots of seizure progression through the four seizure epochs. Colours represent the spatial organization of the epochs shown in time series data of Fig. 1. Dashed-line squares represent the MEA footprint, and dashed-line circles represent ECoG electrodes. Yellow ECoG electrodes represent those exhibiting hypersynchronous ictal activity in the low-frequency (0.350 Hz) LFP. As the recruited region (seizure core) gradually expands to envelope the MEA site, the ictal wavefront approaches the site, passes through it and then continues to move away from the site. (b) Schematic depictions of travelling wave activity across the MEA microelectrodes during each of the four seizure epochs. Rapidly moving travelling waves across the MEA change direction based on the MEAs location relative to the slowly moving ictal wavefront.

low-frequency LFP across the brain, which is a by-product of rapidly travelling waves originating at the ictal wavefront. These results describe the spatiotemporal structure of seizure generation and termination in human chronic epilepsy and have implications for targeted treatments aimed at stopping seizures.

Methods

Patients and ethics statement. Study participants consisted of adults with pharmacoresistant focal epilepsy who underwent chronic invasive EEG studies to help identify the epileptogenic zone for subsequent removal. During this process, these patients had arrays of ECoG electrodes implanted across large sections of their brains to localize their seizure focus. These patients also had a 4 4 mm MEA

with 96 (10 10 electrode grid) 1-mm long penetrating microelectrodes (MEA,

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

Blackrock Microsystems, Inc., Salt Lake City, UT) implanted into the clinically dened seizure onset zone, along with the subdural clinical ECoG electrodes. Histology work has previously conrmed that these electrodes were positioned in layers 4 and 5 in the human-association cortex13. However, it is likely that there is some variability due to individual variations in cortical topology, as well as to tissue changes stemming from the implant itself. House et al.65 provide a detailed description of the surgical procedures for implanting this particular type of MEA into humans. The study was approved by the Institutional Review Boards of the Columbia University Medical Center and University of Utah, and informed consent was obtained by each patient before implantation.

Signals from the MEA were acquired continuously at 30 kHz per channel(0.3 Hz7.5 kHz bandpass, 16-bit precision, range 8 mV). The reference was either subdural or epidural, chosen dynamically based on recording quality. Subdural EEG signals were acquired using a standard clinical video EEG system (XLTek, Natus Medical Inc., Oakville, ON, Canada) at 5002000 Hz per channel(0.5 Hz high pass, low-pass set to 14 sampling rate, 24-bit precision), and then subsequently downsampled to 500 Hz. The two data sets were aligned using a pulse-coded signal delivered simultaneously to digital input ports of both recording systems.

Although all of the MEAs examined in this study were implanted into the clinically dened seizure onset zone, none recorded seizure initiation, and only three MEAs were recruited into the locally seizing territory, termed the ictal core. The ictal core, dened as the locus of intensely discharging neuronal ring during seizures, was identied as described in ref. 13. Six seizures from three patients in which the MEA was positioned within the ictal core were examined. In one of these seizures, the MEA recording (but not the ECoG recording) was interrupted after recruitment into the seizure and before seizure termination. This seizure is therefore excluded from direct analyses of the pre-termination epoch. Ten seizures from two patients in which the MEA was not recruited into the ictal core were also examined.

Data collection and preprocessing. As sorting single units is not possible in the ictal core due to the loss of recognizable waveforms after ictal wavefront passage33, timestamps for MUA were instead derived using previously established detection methods13. MUA was segregated from the raw LFP by band-pass ltering the raw LFP between 500 and 3,000 Hz using a 150th order nite impulse response (FIR) lter. Negative peaks in this signal were detected and those peaks that exceeded four times the s.d. of the signal in the negative direction, and that occurred 41 ms after immediately preceding peaks, were retained as multiunit timestamps and waveforms. Any peaks that exceeded eight times the s.d. of the ltered MUA signal were rejected as artefact. Firing rate was calculated from these multiunit timestamps in 100-ms windows every 25 ms through the duration of the seizure.

High-g was segregated from the raw LFP by band-pass ltering the raw LFP between 50 and 300 Hz for MEA data and 50 and 150 Hz for ECoG data, based on acquisition sampling rates, using a 150th order FIR lter. The lower range of this interval was chosen to include the low high-g frequencies observed towards the end of seizures. High-g instantaneous amplitude was dened as the absolute valueof the Hilbert transform of this signal. To examine LFP spectral content, instantaneous amplitude and phase were rst extracted by Morlet wavelet decomposition on 84 scales from 1 to 300 Hz for MEA data and 1 to 150 Hzfor ECoG data. The absolute values of these wavelet decompositions were then normalized by the inverse of the frequency spectrum analysed.

Natural epochs in human seizures. To understand the dynamics of the ictal core leading up to seizure termination, it was necessary to identify transitions local ictal activity. Accordingly, the MEA seizure recordings were segmented into four distinct epochs between the xed times of seizure onset and termination. Seizure onset and termination were determined clinically from the low-frequency ECoG, by consensus of two epileptologists (C.A.S. and L.M.B.).

The ictal wavefront epoch was determined by smoothing the ring rate of each channel with a Gaussian kernel with a 200-ms s.d. The width of this kernel emphasized the tonic ring of the ictal wavefront and de-emphasized the faster bursting activity of the ictal core. The ictal wavefront epoch was dened as the time of the rst channels mean minus its s.d. until the last channels mean plus its s.d. The pre-recruitment period was operationally dened as the time of seizure onset until the start of the ictal wavefront epoch. The post-recruitment and pre-termination periods comprised the interval between the end of the ictal wavefront epoch and seizure termination. The division point between the post-recruitment and pre-termination epochs was determined by calculating the coefcient of variation over the IDIs after end of the ictal wavefront epoch. The coefcient of variation was calculated as the s.d. of groups of 30 IDIs divided by the mean of those IDIs, sliding 3 IDIs at a time through the duration of the seizure. If the coefcient of variation exceeded 5% for any group of 30 IDIs, then the time of the rst IDI in that group dened the division between post-recruitment and pre-termination epochs.

Travelling waves and large-scale LFP synchrony. The direction of the travelling wave in each discharge was determined by ranking MEA or ECoG channels according to the time of occurrence of each ictal discharge, using the broadband LFP signal. Time of occurrence was measured using the minimum of the rst derivative of the low-frequency LFP (o50 Hz); that is, the point of maximally

decreasing slope. A similar method, which emphasizes fast components ofthe discharge, has been previously employed to detect interictal discharge propagation15,28. Positive and negative weighted centroids were then determined for the resulting channel rankings across the MEA or ECoG grid footprint space. The resultant of these two points dened the velocity vector for the travelling wave. Rayleighs test for non-uniformity was used to determine whether the vectors were normally distributed on the unit circle. Kuipers test was used to compare distributions of wave direction angles. Descriptive circular statistics were calculated using the circular statistics toolbox for Matlab66.

The direction of the ictal wavefront and directions of travelling waves were determined in similar ways, but with different temporal granularity. The peak ring rate (calculated using the 200 ms kernel mentioned above) was ranked for each channel during the ictal wavefront epoch and the centroid method mentioned above was carried out on these peaks to calculate propagation direction for the ictal wavefront. Directions of individual travelling waves could also be determined from MUA in a similar manner, but using a kernel with a 10-ms s.d. to emphasize the fast ictal discharges (Supplementary Fig. 6).

Broadband LFP synchrony was evaluated in a similar manner to previous studies67. Principal component analysis was performed on 1 s of data, from all ECoG electrodes on the subdural grids overlying the MEAs, every 10 ms through the duration of the seizure. Synchrony was measured as the amount of variance explained by the rst-principal component (that is, how much variance in the ECoG signal can be explained by a single covariance dimension). Trends in the slopes and mean values of this measure of synchrony were quantied during the pre-termination epoch using linear regression.

Ictal core and wavefront simulation. To gain a theoretical understanding of the ictal wavefront on seizure electrophysiology, the relatively small network of reactive excitatory cells in the ictal core was simulated using a mean-eld approach based on previous computational studies31,32. The simulated ictal core had recurrent, intra-network and extra-network inputs. Intra-network input was subject to ring rate-dependent synaptic depression, as in previous studies31,32,61. The extra-network input was used to model the contribution of the ictal wavefront to the network of cells in the ictal core. Extra-network input was linearly decreased in a biophysically plausible range for seizure activity, and the temporal dynamics of the simulated ictal core were observed. Details of the model, including cellular parameters (Supplementary Table 1), can be found in Supplementary Information.

Per-discharge analyses. Due to the typical ictal pattern featuring trains of transient high-amplitude discharges with relative attenuation of signal in between, all analyses were conducted on a per-discharge basis. Ictal discharges were rst isolated by identifying all local maxima in both mean high-g amplitude and mean multiunit ring rate across MEA channels that exceeded twice the s.d. above the mean computed from the entire seizure. Some peaks in multiunit or high-g activity appeared during the seizure without a peak in the other signal. The operational denition of discharges in the ictal core therefore required nearly simultaneous peaks in both the high-g activity and MUA; specically within a window dened by the midpoint of the IDI preceding and following each discharge. Over all seizures in which the MEA was recruited into the ictal core, 1.7% of discharges were rejected due to this criterion (12 peak detections out of 708 peak detections, recorded from six seizures). Per-discharge correlations among ring rate and high-g measures were evaluated with Spearmans r. IDIs were calculated from the mean ring rate and high-g amplitude over all recorded MEA channels by subtracting each successive time at which peak ring rate or high-g occurred.

We used the discharge width to quantify the temporal dispersion of multiunit ring. Discharge width was measured by rst converting multiunit timestamps into N M binary matrices, where N represented the number of channels on the array

and M represented the number of milliseconds from seizure start through termination. Each element of the matrix was either one or zero, representing the presence or absence of an action potential on that specic channel during that specic millisecond. The sum of this matrix over the channel dimension for each discharge, as dened by the aforementioned criteria, was t to a Gaussian and the s.d. of the t was operationally dened as the discharge width (for example, black dotted lines in Fig. 7a). Discharges with R2 values o0.9 were excluded from this analysis (9.2% of discharges, or 64 of 696 discharges).

Mutual information was used as a measure of multiunit timing coincidence68. Mutual information for action potential timing was measured by taking 100-bit binary words from each channel in the aforementioned action potential matrix. Mutual information was then calculated for these 100-bit words every 25 ms through the duration of the seizure for each pairwise combination of channels across the MEA. The discharge amplitude was dened as the peak median mutual information across combinations of channels normalized by the ring rate atthat point.

References

1. Peneld, W. & Erickson, T. C Epilepsy and Cerebral Localization (ed. Thomas, C. C.) 623 (Oxford, UK, 1941).

2. Kramer, M. A. & Cash, S. S. Epilepsy as a disorder of cortical network organization. Neuroscientist 18, 360372 (2012).

10 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098 ARTICLE

3. Kramer, M. A. et al. Coalescence and fragmentation of cortical networks during focal seizures. J. Neurosci. 30, 1007610085 (2010).

4. Kramer, M. A. et al. Human seizures self-terminate across spatial scales via a critical transition. Proc. Natl Acad. Sci. USA 109, 2111621121 (2012).5. Schindler, K., Elger, C. E. & Lehnertz, K. Increasing synchronization may promote seizure termination: evidence from status epilepticus. Clin. Neurophysiol. 118, 19551968 (2007).

6. Jiruska, P. et al. Synchronization and desynchronization in epilepsy: controversies and hypotheses. J. Physiol. 591, 787797 (2013).

7. Timofeev, I. & Steriade, M. Neocortical seizures: initiation, development and cessation. Neuroscience 123, 299336 (2004).

8. Krishnan, G. P. & Bazhenov, M. Ionic dynamics mediate spontaneous termination of seizures and postictal depression state. J. Neurosci. 31, 88708882 (2011).

9. Lado, F. A. & Mosh, S. L. How do seizures stop? Epilepsia 49, 16511664 (2008).

10. Pinto, D. J., Patrick, S. L., Huang, W. C. & Connors, B. W. Initiation, propagation, and termination of epileptiform activity in rodent neocortex in vitro involve distinct mechanisms. J. Neurosci. 25, 81318140 (2005).

11. Ziemann, A. E. et al. Seizure termination by acidosis depends on ASIC1a. Nat. Neurosci. 11, 816822 (2008).

12. DAmbrosio, R. & Miller, J. W. What is an epileptic seizure? Unifying denitions in clinical practice and animal research to develop novel treatments. Epilepsy Curr. 10, 6166 (2010).

13. Schevon, C. A. et al. Evidence of an inhibitory restraint of seizure activity in humans. Nat. Commun. 3, 1060 (2012).

14. Weiss, S. A. et al. Ictal high frequency oscillations distinguish two types of seizure territories in humans. Brain 136, 37963808 (2013).

15. Trevelyan, A. J., Baldeweg, T., van Drongelen, W., Yuste, R. & Whittington, M. The source of afterdischarge activity in neocortical tonic-clonic epilepsy.J. Neurosci. 27, 1351313519 (2007).16. Ma, H., Zhao, M. & Schwartz, T. H. Dynamic neurovascular coupling and uncoupling during ictal onset, propagation, and termination revealed by simultaneous in vivo optical imaging of neural activity and local blood volume. Cereb. Cortex 23, 885899 (2013).

17. Ellender, T. J., Raimondo, J. V., Irkle, A., Lamsa, K. P. & Akerman, C. J. Excitatory effects of parvalbumin-expressing interneurons maintain hippocampal epileptiform activity via synchronous afterdischarges. J. Neurosci. 34, 1520815222 (2014).

18. Truccolo, W. et al. Neuronal ensemble synchrony during human focal seizures.J. Neurosci. 34, 99279944 (2014).19. Trevelyan, A. J., Sussillo, D., Watson, B. O. & Yuste, R. Modular propagation of epileptiform activity: evidence for an inhibitory veto in neocortex. J. Neurosci. 26, 1244712455 (2006).

20. Schevon, C. A. et al. Spatial characterization of interictal high frequency oscillations in epileptic neocortex. Brain 132, 30473059 (2009).

21. Kellis, S. et al. Multi-scale analysis of neural activity in humans: implications for micro-scale electrocorticography. Clin. Neurophysiol. 127, 591601 (2015).

22. Normann, R. A. Technology insight: future neuroprosthetic therapiesfor disorders of the nervous system. Nat. Clin. Pract. Neurol. 3, 444452 (2007).

23. Trevelyan, A. J., Sussillo, D. & Yuste, R. Feedforward inhibition contributes to the control of epileptiform propagation speed. J. Neurosci. 27, 33833387 (2007).

24. Peneld, W. & Jasper, H. Epilepsy and the functional anatomy of the human brain (Little, Brown & Co, 1954).

25. Benucci, A., Frazor, R. A. & Carandini, M. Standing waves and traveling waves distinguish two circuits in visual cortex. Neuron 55, 103117 (2007).

26. Zanos, T. P., Mineault, P. J., Nasiotis, K. T., Guitton, D. & Pack, C. C.

A sensorimotor role for traveling waves in primate visual cortex. Neuron 85, 615627 (2015).

27. Sato, T. K., Nauhaus, I. & Carandini, M. Traveling waves in visual cortex. Neuron 75, 218229 (2012).

28. Alarcon, G. et al. Origin and propagation of interictal discharges in the acute electrocorticogram. Implications for pathophysiology and surgical treatment of temporal lobe epilepsy. Brain 120, 22592282 (1997).

29. Emerson, R. G., Turner, C. A., Pedley, T. A., Walczak, T. S. & Forgione, M. Propagation patterns of temporal spikes. Electroencephalogr. Clin. Neurophysiol. 94, 338348 (1995).

30. Gonzlez-Ramrez, L. R., Ahmed, O. J., Cash, S. S., Wayne, C. E. & Kramer, M.A. A biologically constrained, mathematical model of cortical wave propagation preceding seizure termination. PLoS Comput. Biol. 11, e1004065 (2015).31. Mongillo, G., Barak, O. & Tsodyks, M. Synaptic theory of working memory. Science 319, 15431546 (2008).

32. Vladimirski, B. B., Tabak, J., ODonovan, M. J. & Rinzel, J. Episodic activity in a heterogeneous excitatory network, from spiking neurons to mean eld.J. Comput. Neurosci. 25, 3963 (2008).

33. Merricks, E. M. et al. Single unit action potentials in humans and the effect of seizure activity. Brain 138, 28912906 (2015).

34. Ray, S. & Maunsell, J. H. R. Differences in gamma frequencies across visual cortex restrict their possible use in computation. Neuron 67, 885896 (2010).

35. Manning, J. R., Jacobs, J., Fried, I. & Kahana, M. J. Broadband shifts in local eld potential power spectra are correlated with single-neuron spiking in humans. J. Neurosci. 29, 1361313620 (2009).

36. Miller, K. J., Sorensen, L. B., Ojemann, J. G., Nijs & den, M. Power-law scaling in the brain surface electric potential. PLoS Comput .Biol. 5, e1000609 (2009).

37. Zanos, T. P., Mineault, P. J. & Pack, C. C. Removal of spurious correlations between spikes and local eld potentials. J. Neurophysiol. 105, 474486 (2011).

38. Ray, S. & Maunsell, J. H. R. Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. PLoS Biol. 9, e1000610 (2011).

39. Cymerblit-Sabba, A. & Schiller, Y. Development of hypersynchrony in the cortical network during chemoconvulsant-induced epileptic seizures in vivo.J. Neurophysiol. 107, 17181730 (2012).40. Bettus, G. et al. Interictal functional connectivity of human epileptic networks assessed by intracerebral EEG and BOLD signal uctuations. PLoS ONE 6, e20071 (2011).

41. Engel, J. & Pedley, T. A. Epilepsy (Wolters Kluwer Health/Lippincott Williams and Wilkins, Philadephia, 2008).

42. Beverlin, B., Kakalios, J., Nykamp, D. & Netoff, T. I. Dynamical changes in neurons during seizures determine tonic to clonic shift. J. Comput. Neurosci. 33, 4151 (2011).

43. Abbott, L. & Van Vreeswijk, C. Asynchronous states in networks of pulse-coupled oscillators. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 48, 14831490 (1993).

44. Engel, J. Models of focal epilepsy. Suppl. Clin. Neurophysiol. 57, 392399 (2004).

45. Bragin, A., Csicsvri, J., Penttonen, M. & Buzski, G. Epileptic afterdischarge in the hippocampal-entorhinal system: current source density and unit studies. Neuroscience 76, 11871203 (1997).

46. Alger, B. & Nicoll, R. Epileptiform burst afterhyperolarization: calcium-dependent potassium potential in hippocampal CA1 pyramidal cells. Science 210, 11221124 (1980).

47. Timofeev, I., Grenier, F. & Steriade, M. Contribution of intrinsic neuronal factors in the generation of cortically driven electrographic seizures.J. Neurophysiol. 92, 11331143 (2004).48. Bal, T. & McCormick, D. A. What stops synchronized thalamocortical oscillations? Neuron 17, 297308 (1996).

49. Norden, A. D. & Blumenfeld, H. The role of subcortical structures in human epilepsy. Epilepsy Behav. 3, 219231 (2002).

50. Mitsuya, K., Nitta, N. & Suzuki, F. Persistent zinc depletion in the mossy ber terminals in the intrahippocampal kainate mouse model of mesial temporal lobe epilepsy. Epilepsia 50, 19791990 (2009).

51. Hamil, N. E., Cock, H. R. & Walker, M. C. Acute down-regulation of adenosine A(1) receptor activity in status epilepticus. Epilepsia 53, 177188 (2012).

52. Van Gompel, J. J. et al. Increased cortical extracellular adenosine correlates with seizure termination. Epilepsia 55, 233244 (2014).

53. Crone, N. E., Miglioretti, D. L., Gordon, B. & Lesser, R. P. Functional mapping of human sensorimotor cortex with electrocorticographic spectral analysis. II. Event-related synchronization in the gamma band. Brain 121, 23012315 (1998).

54. Buzski, G., Horvth, Z., Urioste, R., Hetke, J. & Wise, K. High-frequency network oscillation in the hippocampus. Science 256, 10251027 (1992).

55. Einevoll, G. T., Kayser, C., Logothetis, N. K. & Panzeri, S. Modelling and analysis of local eldpotentials for studying the function ofcortical circuits. Nat. Rev. Neurosci. 14, 770785 (2013).

56. Miller, K. J. Broadband spectral change: evidence for a macroscale correlate of population ring rate? J. Neurosci. 30, 64776479 (2010).

57. Bragin, A., Benassi, S. K., Kheiri, F. & Engel, J. Further evidence that pathologic high-frequency oscillations are bursts of population spikes derived from recordings of identied cells in dentate gyrus. Epilepsia 52, 4552 (2011).

58. Modur, P. N., Vitaz, T. W. & Zhang, S. Seizure localization using broadband EEG: comparison of conventional frequency activity, high-frequency oscillations, and infraslow activity. J. Clin. Neurophysiol. 29, 309319 (2012).

59. Zelmann, R., Lina, J. M., Schulze-Bonhage, A., Gotman, J. & Jacobs, J. Scalp EEG is not a blur: it can see high frequency oscillations although their generators are small. Brain Topogr. 27, 683704 (2014).

60. Bragin, A., Engel, J. & Staba, R. J. High-frequency oscillations in epileptic brain. Curr. Opin. Neurol. 23, 151156 (2010).

61. Wilson, H. R. & Cowan, J. D. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 124 (1972).

62. Khodagholy, D. et al. NeuroGrid: recording action potentials from the surface of the brain. Nat. Neurosci. 18, 310315 (2015).

63. Gallmetzer, P. et al. Postictal paresis in focal epilepsies--incidence, duration, and causes: a video-EEG monitoring study. Neurology 62, 21602164 (2004).

NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11098

64. Bateman, L. M., Li, C.-S. & Seyal, M. Ictal hypoxemia in localization-related epilepsy: analysis of incidence, severity and risk factors. Brain 131, 32393245 (2008).

65. House, P. A., MacDonald, J. D., Tresco, P. A. & Normann, R. A. Acute microelectrode array implantation into human neocortex: preliminary technique and histological considerations. Neurosurg. Focus 20, E4 (2006).

66. Berens, P. CircStat: a MATLAB toolbox for circular statistics. J. Stat. Softw. 31, 120 (2009).

67. Schindler, K., Leung, H., Elger, C. E. & Lehnertz, K. Assessing seizure dynamics by analysing the correlation structure of multichannel intracranial EEG. Brain 130, 6577 (2007).

68. Shannon, C. A mathematical theory of communications. Bell Syst. Tech. J. 27, 379423 (1948).

Acknowledgements

We thank Lawrence F. Abbott and Charles E. Schroeder for helpful advice on the development of this manuscript. This work funded by NIH/NINDS R01 NS084142 (C.A.S.).

Author contributions

E.H.S., A.T. and C.A.S. wrote the manuscript. E.H.S. analysed the seizure data. J.-y.L. did the computational modelling. E.H.S., C.A.S., R.G.E., T.S.D. and B.G. collected the data.

C.A.S., E.M.M., B.G., L.M.B. and A.J.T. provided developmental ideas. P.A.H., R.R.G. and G.M.M. implanted the MEAs. All authors contributed in editing the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions

Web End =reprintsandpermissions/

How to cite this article: Smith, E. H. et al. The ictal wavefront is the spatiotemporal source of discharges during spontaneous human seizures. Nat. Commun. 7:11098 doi: 10.1038/ncomms11098 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

12 NATURE COMMUNICATIONS | 7:11098 | DOI: 10.1038/ncomms11098 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications