BOLD changes occur prior to epileptic spikes seen on scalp EEG

www.elsevier.com/locate/ynimg NeuroImage 35 (2007) 1450 – 1458

BOLD changes occur prior to epileptic spikes seen on scalp EEG Colin S. Hawco, a Andrew P. Bagshaw, a,b Yingli Lu, a François Dubeau, a and Jean Gotman a,⁎ a

Montreal Neurological Institute and Hospital, McGill University, 3801 University Street, Montreal, Quebec, Canada H3A 2B4 School of Psychology, University of Birmingham, Birmingham, B15 2TT, UK

b

Received 11 August 2006; revised 30 November 2006; accepted 5 December 2006 Available online 25 January 2007 This study examined BOLD changes prior to interictal discharges in the EEG of patients with epilepsy. From a database of 143 EEG– fMRI studies, we selected the 16 data sets that showed both strong fMRI activation in the original analysis and only a single spike type in the EEG. Scans were then analyzed using seven model HRFs, peaking 3 or 1 s before the event, or 1, 3, 5, 7, or 9 s after it. An HRF was calculated using a deconvolution method for all activations seen in each analysis. The results showed that seven data sets had HRFs that peaked 1 s after the event or earlier, indicating a BOLD change starting prior to the spike seen on the scalp EEG. This is surprising since the BOLD change is expected to result from the spike. For most of the data sets with early peaking HRFs, the maximum activation in all of the statistical maps was when the model HRF peaked 1 s after the event, suggesting that the early activation was at least as important as any later activation. We suggest that this early activity is the result of neuronal changes occurring several seconds prior to a surface EEG event, but that these changes are not visible on the scalp. This is the first report of a BOLD response occurring several seconds prior to an interictal event seen on the scalp and could have important implications for our understanding of the generation of epileptic discharges. © 2007 Elsevier Inc. All rights reserved.

Introduction Epilepsy is a pervasive neurological disorder that affects up to 3% of people by age 75 (Hauser et al., 1993). It is characterized by observable clinical seizures, often involving uncontrolled motor activity and/or loss of consciousness, and abnormalities in the electroencephalogram (EEG). These abnormalities take the form of continuous rhythmic activity during the ictal period (when a seizure is actually occurring) and are often present between clinical seizures (the interictal period) in the form of EEG discharges, isolated spikes, or bursts of spiking activity. These interictal events (or spikes) are not accompanied by clinical manifestations, loss of consciousness, or seizure activity and can be useful in the ⁎ Corresponding author. Fax: +1 514 398 8106. E-mail address: [email protected] (J. Gotman). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter © 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2006.12.042

diagnosis of the disorder and the localization of the seizure focus (the area of the brain from which the seizures originate). Interictal events can be divided into two basic categories: focal spikes, thought to originate from a specific brain region, and generalized spikes, typically thought of as originating from a more widespread cortical network (Blumenfeld, 2005). Functional magnetic resonance imaging (fMRI) is a noninvasive method that allows indirect examination of neuronal activity. In many cases, fMRI studies utilize exogenous experimental stimuli to examine a neurological process. However, an endogenous event, such as an interictal epileptic discharge, can also be studied using fMRI (Gotman et al., 2004; Hamandi et al., 2004). Recording the EEG simultaneously with fMRI allows the timing of the spikes to be locked to the blood oxygen level dependant (BOLD) signal, allowing the examination of BOLD changes corresponding to the epileptic activity. Typically, the hemodynamic response to an event is modeled based on an expected response. One of the most commonly used hemodynamic response functions (HRF) is that from Glover (1999). This assumes a rise after the event with a peak at 5.4 s and a smaller undershoot peaking after approximately 12 s. Several studies examining the variability of the BOLD response have suggested that a single HRF model may not be optimal. Aguirre et al. (1998) examined the variability of the HRF across different subjects, and across different scanning sessions within the same subject. They found substantial variability between subjects in the shape of the HRF, but little variability within the same subject across different scanning sessions. Smaller, but not negligible, differences have also been observed in the HRF from different brain regions of the same subject (Miezin et al., 2000; Neumann et al., 2003; Handwerker et al., 2004). Others have noted that the amplitude of the HRF can vary not only across subjects, but also across time within a single subject during one scanning session (Menz et al., 2006), though again the variability was found to be less than the variability across subjects. Most of the work on HRF variability has involved cognitive tasks with healthy volunteers, and the effects of brain abnormalities on the BOLD signal are poorly understood. However, studies suggest that HRFs associated with epileptic events do not differ substantially from the response to exogenous stimuli (Lemieux et al., 2001; Bénar et al., 2002; Lu et al., 2006).

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

Bénar et al. (2002) examined the hemodynamic response in epileptic patients in voxels found to be active in an analysis using the standard Glover HRF. Despite the fact that the use of a single linear model causes a bias for the detected voxels to have an HRF similar to the model, substantial variability was noted across patients. Bénar et al. (2002) also noted that there were differences in the amplitude of the BOLD response between patients. Lu et al. (2006) used a deconvolution method which makes no assumptions about the shape of the hemodynamic response and also found considerable variation. In order to take the variability of the HRF into consideration, Bagshaw et al. (2004) conducted a study in which four model HRFs were used, with peak times of 3, 5, 7, or 9 s after the event. It was found that the use of four models significantly improved the detection of activation. This range of model HRFs helps compensate for variation across individuals, and within individuals across different brain regions. Recent evidence has shown that BOLD changes may occur prior to the onset of ictal activity as assessed by visual inspection of the scalp EEG. Federico et al. (2005) scanned three patients who had a typical partial epileptic seizure while in the scanner, and changes in the BOLD signal were seen in the period leading up to the seizures. The localization of these changes was consistent with the clinical findings for these patients. The observed changes in BOLD signal occurred as early as 20 min prior to the EEG seizure onset. A pre-ictal change in BOLD signal is consistent with other findings that have suggested the existence of a pre-ictal period using SPECT (Baumgartner et al., 1998) and EEG (Elger and Lehnertz, 1998; Navarro et al., 2005). Changes in BOLD responses have also been seen to precede spiking in an animal model of penicillin-induced epilepsy (Makiranta et al., 2005). Penicillin was injected into the somatosensory cortex while simultaneous EEG–fMRI was recorded in pigs under deep anesthesia. The timing of changes following the injection was noted in both the EEG and the BOLD signals. An increase in BOLD signal was observed almost immediately after the injection of penicillin, prior to the emergence of spikes on the EEG (which were observed 49 s to 2 min 44 s after injection). Spikes were very frequent and continuous after onset, and the BOLD signal increase did not plateau until some time after spike onset. The signal increases of 9.4% to 21.3% were larger than those typically seen in combined EEG–fMRI studies of interictal events in humans, which tend to be on the order of 1–2% (Lemieux et al., 2001; Bénar et al., 2002). The continuous nature of this spiking was also different from the isolated and sporadic events typically examined in EEG–fMRI studies in humans. While several studies have examined HRFs in response to interictal spikes (Lemieux et al., 2001; Bénar et al., 2002; Bagshaw et al., 2005; Lu et al., 2006), none has reported BOLD changes occurring before the observation of the interictal event on scalp EEG in humans. There is no evidence of a neuronal event preceding an interictal discharge. During examination of our data we observed that some patients seem to have a BOLD signal change prior to an interictal EEG discharge. Following this observation, we decided to study a group of patients who were known to have had activation in a previous analysis (using the methods of Bagshaw et al., 2004) to see if these early changes could be observed, and to determine how they compared with the activations seen in standard EEG–fMRI studies.

1451

Materials and methods Subjects Data sets were selected from a database of 143 scanning sessions from 130 patients according to the following selection criteria: (1) Consistent event localization in the EEG during the EEG– fMRI scanning session (i.e. all discharges during the session were localized to the same electrodes). (2) Strong fMRI activation on a previous analysis (using the methods of Bagshaw et al., 2004), defined as a cluster of 5 contiguous voxels above a t value of 3.1 (p < 0.01 for cluster extent, corrected for multiple comparisons), with at least one voxel in the cluster with a t value above 4.8 (at least one voxel in each cluster having a p < 0.05 at the level of a single voxel, corrected for multiple comparisons). This ensured that all data sets had a very significant activation and minimized the risk of the response being artifactual. (3) No movement in excess of 5 mm (1 voxel). Sixteen data sets from 15 patients met the above criteria and were included in the analysis, with one patient having been scanned twice. The patients studied had a variety of clinical syndromes and spike morphologies and included 7 patients with idiopathic generalized epilepsy (IGE), 1 with secondary generalized epilepsy (who was scanned twice), and 7 with focal epilepsy. Clinical details and spike distribution are presented in Table 1. Data sets were divided into two groups: those in which all events were single isolated spikes (7 data sets), and those that had bursts of spikes (i.e. some or all events consisted of a run of spikes; 9 data sets). Data sets in the second group could have a combination of bursts and single spikes. This separation of data sets was done to distinguish patients with prolonged events, who are known to have a BOLD response that varies as a function of the length of the event (Bagshaw et al., 2005), from patients with isolated events, who are expected to have a more consistent BOLD response. Patients gave their informed consent for this study, in agreement with the policies of the Research Ethics Board of the Montreal Neurological Institute and Hospital. Data acquisition Continuous EEG–fMRI data were acquired for each patient. All patients underwent one study except one who had a second study approximately 3 months after the first (data sets G and G2). Functional MRI images were acquired on one of two 1.5 T MRI scanners (Vision or Sonata, Siemens, Erlangen, Germany) using echo-planar imaging (voxel dimensions 5 × 5 × 5 mm, 25 slices, 64 × 64 matrix, TE of 50 ms, TR of 3 s and a flip angle of 90°). Images were acquired in runs of 120 frames of 25 slices lasting approximately 6 min with a short pause between runs when no fMRI data were collected. Between 4 and 13 runs were collected for each patient. A T1-weighted scan was also obtained for anatomical localization of the functional data (1-mm slice thickness, 256 × 256 matrix, TE = 9.2 ms, TR = 22 ms, flip angle 30°). The total scan time for each patient was up to 2 h. Images were motion corrected, with all images realigned to the third frame

1452

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

Table 1 Data set

Clinical syndrome

Spike distribution

MRI

Age at onset

Age at testing

Events

Runs

Duration of bursts (s)

A

TLE

LT

35

42

44

13

0

B C D E F

TLE? PLE? IGE IGE TLE FCE

LP Generalized LT RT LF

5 19 ? 8 3

50 66 ? 45 44

25 12 4 16 5

10 10 9 10 8

0 0 0 0 0

G G2 H I J K L M N O

SGE SGE IGE TLE IGE TLE IGE IGE CPE IGE

BiF BiF Generalized LF Generalized LT Generalized Generalized BiF Generalized

Multiple cavernomas, the 3 more prominent are on L frontopolar region, L anterior T region, and splenium Subcortical LF cyst Hydrocephaly Normal RTL resection Porencephalic cyst lateral to lateral ventricle surrounded by gliosis at R RT arachnoid cyst RT arachnoid cyst Normal Gliosis on LT pole Normal BPPMG Normal Normal R perisylvian PMG, RT resection Normal

5 5 10 17 ? 0.5 12 14 0.5 12

30 30 42 25 ? 26 22 26 37 36

14 12 76 24 4 21 56 70 8 8

6 8 10 10 4 10 10 10 9 8

0 0–6 0–2 0–10 1–4 0–5 1–6 0–4 1–4 4–7

Abbreviations: L = left, R = right, F = frontal, T = temporal, P = parietal, BiF = bifrontal, PMG = polymicrogyri. IGE = idiopathic generalized epilepsy, SGE = secondary generalized epilepsy, CPE = central parietal epilepsy, BPPMG = Biparietal polymicrogyri.

of the first run. A 6 mm FWHM smoothing kernel was applied to the data. EEG data were collected using 21 MRI compatible Ag/AgCl electrodes placed according to the international 10–20 system using either an EMR32 (Schwarzer, Munich, Germany) or a BrainAmp (BrainProducts, Munich, Germany) amplifier. The patient's head was immobilized during the scan using a plastic bag filled with small polystyrene spheres, in which a vacuum was obtained by air suction (S&S X-ray products, Brooklyn, NY). The EEG was filtered offline to remove the artifact generated by the MRI scans, allowing the entire EEG trace to be examined. FEMR software (Schwarzer, Munich, Germany; Hoffman et al., 2000) was used for the EEG collected using the EMR32 data, while Vision Analyzer version 1.04 (BrainProducts, Munich, Germany) was used for the BrainAmp EEG data (Allen et al., 2000). An experienced electroencephalographer marked the spatial distribution, morphology, duration and timing of the epileptiform discharges. In the cases with bursts of epileptiform discharges, the spatial distribution was stable, with no propagation during the course of the burst. No patients had ictal activity while being scanned.

previous analysis (peaks 3, 5, 7, and 9) as well as early peaking HRFs (peak− 3, peak− 1, and peak1) for the investigation of both “standard” and early responses. The full width at half maximum for all model HRFs was 5.2 s. Due to limitations inherent in the statistical package used, it was necessary to remove events occurring in the first 6 s of each run to accommodate the transformations necessary to create the early peaking model HRFs. In the case of an event that consists of a burst of spikes, the linear model predicts that the HRF can be calculated by convolving the HRF for a single spike with a boxcar function of length equal to the duration of the event. This has been shown to be valid for epileptic discharges (Bagshaw et al., 2005). The duration of the event was convolved with the model for bursts of spikes. The correction for the multiple comparisons due to the large number of voxels (8000 voxels in the brain, assuming a brain size of 1000 cc and using voxel size 5 mm × 5 mm × 5 mm) was done using random field theory (Worsley et al., 1996), while a Bonferroni

Statistical analysis Statistical analysis of the fMRI data was performed using the methods of Worsley et al. (1996, 2002). Models were pre-whitened with an AR model of order 1, and slow fluctuations in the data were taken into consideration by fitting a third order polynomial to each run. The first 3 frames were removed from each run to ensure the magnetization was in a steady state. Timing of the spikes, in seconds, was used for an event related analysis, with effects of the order of slice acquisition taken into account. For each patient, seven statistical maps were created using seven model HRFs consisting of a single positive gamma function peaking 3 s before spike onset (peak− 3), 1 s before spike onset (peak− 1), or 1, 3, 5, 7, or 9 s after spike onset (peak1, peak3, peak5, peak7, and peak9 respectively; see Fig. 1). We included the HRFs used in the

Fig. 1. The seven model hemodynamic response functions (HRFs) used in the analysis. Each HRF was used in a separate, distinct analysis of the data.

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

correction was used to account for the 7 analyses. For each of the seven statistical maps generated for each data set, all voxels showing a t value higher than 3.1 (p < 0.001 non corrected) were selected to test the spatial extent of the clusters using random field theory (Worsley et al., 1996). The correction for multiple comparisons stemming from the presence of 8000 voxels in the brain resulted in the requirement of a cluster of 3 contiguous voxels above 3.1 in order to obtain a corrected p < 0.01. In order to take into account the 7 analyses, we performed a Bonferonni correction. This results in requiring a significance level of 0.01/7 = 0.0014 for each cluster. This is fulfilled by requiring 5 voxels above 3.1, which therefore results in a corrected significance level of 0.01. To exclude weak responses and help ensure that all considered activations are robust, we applied a second criterion requiring at least one voxel within the cluster to have a t value over 4.8, which corresponds to p < 0.05 for a single voxel, also corrected for 8000 voxels. It also ensured that the voxel with the highest t value in any cluster of activation, which was used to calculate an HRF (see below), would be significant on its own, independently of its neighbors in the cluster. Maps were masked, such that all voxels outside the brain were excluded. Statistical maps were obtained independently for each model HRF. Clusters with a positive t value, corresponding to an increase in BOLD signal, were considered activations, while clusters with a negative t value were considered deactivations. Deactivations were not considered further in this study for a number of reasons. The HRF associated with deactivations tends to be noisier than that for activations, making it more difficult to quantify the results. Although recent work has confirmed the link between negative BOLD responses and decreases in neuronal activity in monkey V1 (Shmuel et al., 2006), the neurophysiological origin of deactivations is generally less well understood than that of activations and they may be the result of a number of different physiological phenomena (Shmuel et al., 2002). In addition, the localization of negative BOLD responses to interictal discharges tends to be less well correlated with the expected cortical origin of the epileptic discharge, making their relationship to the epilepsy unclear (Kobayashi et al., 2006a). Given the difficulties in quantifying and interpreting the HRF in response to deactivations, it was decided to focus on activations. Activations were examined to determine concordance with the expected location of the epileptic focus (as determined by an experienced epileptologist). An activation was considered concordant if it was located in the same lobe as the expected focus. In the case of patients with generalized bursts, activity was considered to be concordant if it was non-focal (i.e. it did not involve only a specific cortical area, but rather had a more widespread and symmetrical cortical distribution).

1453

the event. This removes the confound of events with varying durations within a single data set. The HRF was calculated for a time window of 64 s, from 15 s prior to the event to 48 s post event, for the single maximal voxel at each cluster of activation. Calculated HRFs were interpolated to one data point per second. To avoid including data that did not have a robustly fitted response, only HRFs with a signal to noise ratio (SNR) greater than 3 were analyzed. The SNR was calculated by taking the amplitude of the peak and dividing it by the standard deviation of the background, defined as 12 to 6 s prior to and 30 to 48 s after the event. The peak was defined as the point of highest amplitude in the time range between 4 s before the event and 15 s after the event. The time to peak for all calculated HRFs was measured with respect to the onset of the scalp EEG event. Calculated HRFs were considered to have an early peak if the time to peak was 1 s after the EEG event, at the same time as the EEG event, or prior to the EEG event. With a FWHM of 5.2 s, the model HRF has its onset 4 or 5 s before its peak. This indicates that an HRF peaking 1 s after the event has an onset 2 or 3 s before the event, or in other words that the change in BOLD response precedes the electrical activity observed on the scalp EEG. Results Statistical maps Table 2 summarizes the results of the statistical maps. All data sets had significant activation, as this was a criterion for inclusion in this study. In all cases, the model HRF that resulted in the statistical map with the largest t statistic was the same as the model HRF that resulted in the largest volume of activations, or the two model HRFs differed by at most 2 s (for instance, largest t statistic in peak3 map and largest volume in peak5 map). The maximum t statistic was found in peak1 in four data sets (25%), in peak3 for five (31%), in peak5 for five (31%), and in peak7 for two (13%). The largest total volume of activation (which included all clusters for each statistical map) was found in peak1 for four (25%), in peak3 for four (25%), in peak5 for three (19%), and in peak7 for four (25%), and it was equal in peak5 and peak7 in one (6%). No data sets had a maximum t statistic or a maximum volume of activations using peak− 3, peak− 1, or peak9. Two data sets without bursts (A and E) and four data sets with bursts (I, K, M, and O) had no activation in the statistical maps generated with the early model HRFs (peak−3, peak−1, peak1). The other five data sets without bursts and five data sets with bursts showed at least one activation cluster in the maps generated with the early model HRFs.

Calculated HRFs Calculated HRFs In order to exclude the possibility that activations in the early statistical maps were related to later peaking BOLD changes (i.e. regions that showed an activation with peak t values in the later maps but with significant voxels in the early maps) and to conduct a closer examination of the timing of the BOLD response, an HRF was calculated for every significant activation in each statistical map for all data sets. The HRF was obtained using a deconvolution method (Cox, 1996; Ward, 2000). This method uses the event timing to estimate the HRF from the BOLD data, with no constraints or assumptions on its shape. The deconvolution method allows us to calculate the instantaneous HRF given the duration of

Of a total of 437 clusters of activation found in all data sets, 330 (75.5%) resulted in a calculated HRF that passed the SNR criterion. In one data set (M), none of the calculated HRFs passed the SNR criterion. Five of the seven data sets that did not include bursts of interictal discharges (B, C, D, F, and G1) and two of the nine data sets with bursts (H and J) had at least one early peaking calculated HRF (with a peak time of 1 s after the event or earlier, as defined above). Of these seven data sets with early peaking calculated HRFs, four had idiopathic generalized epilepsy, one had secondary generalized epilepsy, one had fronto-central epilepsy,

1454

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

Table 2 Results of the fMRI analysis Data set

Total volume

Maximum t statistics

Early HRF?

Expected focus location

Main fMRI activations

5.6 5.6 13.3 5.4 5.2 10.3 7.7

n y y y n y y

LT LTP* n/a n/a* RCT* R hemisphere Bilateral, maybe left

Peak5

8.1

n

Bilateral, maybe left

38.25 16.88 90.50 68.25 46.13 14.00 33.13

Peak3 Peak3 Peak1 Peak7 Peak5 Peak3 Peak3

8.2 5.9 8.8 9.3 12.3 7.1 10.2

y n y n n n/a n

n/a LT* n/a LTP n/a n/a* RFC

24.75

Peak7

6.8

n

n/a*

Left superior temporal, biparietal Bilateral superior central–parietal Bilateral frontal and thalamus Anterior commissure Left frontal Widespread cerebral (more left than right) and thalamus Bilateral fronto-temporal, occipital midline, left fronto-central Bilateral fronto-temporal (maximum left), occipital, bilateral globus pallidus Bilateral thalamus, midline occipital, and centro-pareital Left frontal, right temporal Anterior, widespread occipital and parietal Superior left temporal, right temporal, bilateral parietal Bilateral occipital and midline frontal Periventricular, thalamus Right parietal (near lesion), left temporal–parietal, left anterior temporal, midline parietal Periventricular, thalamus

Model

Volume (cm3)

Model

t statistics

A B C D E F G

Peak7 Peak5 Peak1 Peak3 Peak7 Peak1 Peak3

30.13 8.00 157.63 5.50 2.25 123.75 39.13

Peak5 or peak7 Peak7 Peak1 Peak1 Peak7 Peak1 Peak5

G2

Peak5

247.13

H I J K L M N

Peak1 Peak3 Peak1 Peak5 Peak3 Peak5 Peak3

O

Peak5

Total volumes include the total volume (calculated by counting the total number of active voxels in each statistical map and converting to cm3) and model HRF that was used for that statistical map; the t statistic is shown for the highest t statistic from each data set, along with the model HRF in which it is found. A “y" for Early HRF means that the patient has at least one HRF that peaks at or before 1 s after the spike. Data set N has “n/a" as no HRFs passed the SNR criteria. The expected focus and actual areas of activation are shown. A “*" indicates that the data set's fMRI results were not concordant with the EEG localization. Abbreviations: L = left, R = right, C = central, F = frontal, T = temporal, P = parietal.

and in one data set it was uncertain if the patient had temporal lobe or parietal lobe epilepsy. Though the small sample size prevented any meaningful statistical analysis, no other factor, such as number of events, number of runs, age at seizure onset, age at time of scanning, or maximum t statistic or volume seem to account for or correlate with the presence or absence of early peaking calculated HRFs. Fig. 2 shows some early peaking activations and their HRFs. Comparison of HRFs and statistical maps Overall, the peak time of the calculated HRFs was similar to the peak time of the model HRF used to generate the statistical map from which the voxels were selected for HRF calculation. For example, the HRFs calculated from clusters seen on the statistical map created using the peak5 model HRF tended to have a peak time of around 5 s. This is to be expected, but lends support to the validity of the calculated HRFs. Activations were only found in a subset of statistical maps (i.e. for a given data set, no cluster of activation was seen in all of the statistical maps created with the seven model HRFs). Early activations were sometimes seen only in early statistical maps (Fig. 3). Some activations seen in early statistical maps did not have an early peaking calculated HRF. These activations generally also had a higher t statistic in maps generated using a later model HRF. An important question is whether the patients with at least one early peaking calculated HRF show their main BOLD response prior to the spike, or whether this early HRF represents a minor event and the main response occurs after the spike. Of the data sets with at least one early peaking calculated HRF, three have maximum t statistic and maximum volume in the statistical map generated using peak1 (C, F, and J), one had a maximum t statistic in peak1 and a maximum volume in peak3 (D), and one had a

maximum t statistic in peak3 and a maximum volume in peak1 (H). Two data sets with an early peaking calculated HRF (B and G) did not have a maximum t statistic nor a maximum volume in peak1: data set B had a maximum t statistic in peak7 and a maximum volume in peak5, and G had a maximum t statistic in peak5 and maximum volume in peak3. Fig. 4 shows the maximum t statistic and maximum volume in each statistical map for all data sets, divided into data sets with at least one early peaking calculated HRF and those with only later peaking calculated HRFs. It appears that, when an early activation is present, it most often represents the main activation. Concordance of EEG and fMRI data Three data sets without bursts (B, D, and E) and three with bursts (I, M, and O) had poor concordance between the expected focus and the overall pattern of fMRI activations. Visual inspection of the activations suggests the concordance between the area of activation and the expected focus is not dependant on the model HRF used. The activations on the early peaking statistical maps appear to be as concordant, qualitatively, as activations on the later peaking statistical maps. The activation associated with early peaking calculated HRFs was concordant with the EEG in 5 of 7 data sets with an early peaking calculated HRF: In the three data sets with generalized spikes (C, H, and J), the early maps showed widespread activity with a posterior distribution; data set G had left frontal activation in the peak1 statistical map, and bilateral frontal activation in peak3, which was concordant with bifrontal spikes; data set B had bilateral parietal activation, with spikes in the left parietal region. The remaining two data sets with early calculated HRFs (F and D) had a poor concordance with the EEG. Data set F had

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

1455

Fig. 2. Activations with early HRFs, showing the area of activation and the calculated HRF. The red arrow indicates the voxel associated with the calculated HRFs. As the amplitude of the calculated HRFs is arbitrary using this method, calculated HRFs are baseline corrected, with the first value set to zero, and scaled such that the maximum amplitude using each method is equal to one. These calculated HRFs show a clear change in the BOLD signal prior to the epileptic scalp discharge.

widespread, scattered activation with left sided predominance. Some activation in peak1 is concordant with the EEG (left frontal), but there are several activations that are not. Data set F also showed thalamic activation in the peak− 3, peak− 1, and peak1 statistical maps. The early activation in data set D is in the corpus callosum, whereas spikes were left temporal. Some left temporal activation was seen on the early model statistical maps, but the calculated HRF for this region did not pass the SNR criterion. Discussion The results demonstrate that, in some patients, BOLD activity is present prior to the observation of an interictal epileptic event on scalp EEG. Five of the seven data sets with early peaking calculated HRFs showed a maximum t statistic or volume of activation in the statistical map generated using the model HRF peaking at 1 s, while the other two showed a maximum t statistic and volume using later model HRFs. This suggests that, in the majority of patients with early peaking calculated HRFs, the early BOLD response is more significant than any later activity that may be present. The implication is that the early BOLD responses seen in these patients

represent the main BOLD manifestation of the spike. Even though this is seen in a minority of patients, it is very intriguing. Neurophysiological implications Why would some patients show a BOLD response that starts several seconds before the scalp EEG event? It has always been assumed that the BOLD signal changes as a consequence of the excessive neuronal discharge constituting the epileptic spike. This is obviously the case in many patients, but in an important subpopulation the change in BOLD cannot be the result of the discharge, although it is tightly coupled with it. Two hypotheses appear plausible: either the BOLD signal change that precedes the neuronal discharge is its cause, or there is a change in neuronal activity that occurs prior to the discharges seen on the scalp, and which results in the BOLD signal change. The first hypothesis would imply that the primary, causal, event in an epileptic discharge is an event that results in a BOLD change but does not result from a change in neuronal activity. This initiating event would therefore be a change in blood flow, blood volume, or oxygen content. There is no evidence in the literature implicating

1456

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

Fig. 3. Two examples of activation that are seen only using the early peaking model HRFs (peaks − 3, − 1, 1). The yellow arrow shows the location of the activation, absent for peak 3 and later peaks.

Fig. 4. Maximum t statistic and total volume of activation for each statistical map for each patient. The graphs are divided into data sets from patients with at least one early peaking calculated HRFs (left), and data sets from patients with no early peaking calculated HRF (right). The black vertical line separates early model HRFs from later model HRFs. The volumes have been normalized to allow visualization of the wide range of values. The largest volume in each data set has been scaled to 1, with all other volumes for that data set scaled as a percentage of the largest volume. Data sets with at least one early peaking calculated HRF predominantly show activation in the early model HRFs (the left of the black line), while data sets with no early peaking HRFs show activation predominantly in later peaking model HRFs. This suggests the early activation represents the main activation in patients with early peaking HRFs.

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

such a mechanism, and the available evidence points towards preserved neurovascular coupling as a result of epileptic discharges (Stefanovic et al., 2005). Examinations of the BOLD response in epilepsy have generally found that the HRF associated with spikes does not differ significantly in quality or overall shape from those of other phenomena, further suggesting an intact neurovascular coupling. Vascular lesions may certainly be the cause of chronic epilepsy (Mottolese et al., 2001; Ettinger, 1994) but abrupt vascular events have not been implicated as a cause for epileptic discharges. We therefore consider this hypothesis highly unlikely. The second hypothesis invokes a neuronal event that is not seen in the scalp EEG but occurs a few seconds before a visible epileptic discharge and results in the “early” BOLD signal change we have illustrated. One can hypothesize two reasons for the “invisibility” of this neuronal event on the scalp EEG. Either the event has spatial or amplitude characteristics such that it is not visible on the scalp but might be visible with intracranial electrodes, or the event does not consist of a synchronized discharge and would not be visible with any type of EEG electrode. Our study does not allow the separation of these two possibilities since no intracranial electrodes were present. Studies comparing scalp and intracranial EEG have most often emphasized the spatial characteristics of the intracranial EEG at the time of a scalp event (Pacia and Ebersole, 1997; Tao et al., 2005), but have not specifically addressed the question raised here: is it frequent that a scalp discharge is preceded by a discharge seen intracranially? It is worth noting, however, that when epileptic spikes propagate, it is usually over a few tens of milliseconds (Emerson et al., 1995), a duration that does not explain the BOLD response starting a few seconds before the scalp spike. We have not seen any report describing spikes propagating over seconds. The second possibility, implying a nonsynchronized neuronal or glial event, is impossible to prove with EEG recordings. This possibility finds support, however, in the study of BOLD changes resulting from an experimental penicillininduced epileptic focus, in which Makiranta et al. (2005) found BOLD increases before the first spike appeared, thus implying that a change in the level of neuronal activity was taking place before it resulted in a synchronized discharge. It was also invoked in our study of a patient with multiple short electrographic epileptic seizures, in whom the most significant BOLD increase took place in a region remote from the EEG discharge, a region that was close to the scalp but showed no EEG change (Kobayashi et al., 2006b), and the results of Federico et al. (2005), which showed BOLD changes long before the onset of a seizure as measured on the EEG. In general, given the different neurophysiological and cellular origins of scalp EEG and the BOLD fMRI signal, it is conceivable that events could occur which are observed with one modality but not the other (Nunez and Silberstein, 2000). It is possible that any early neuronal activity that was not seen on the scalp EEG would effectively result in a prolonged burst which would be reflected in the shape of the overall BOLD response (i.e. unseen plus detected activity separated by some seconds would be merged into one BOLD response due the sluggishness of the HRF). It is difficult to assess if this is the case using the current methods and with the current amount of data. Given the variability of the size and shape of the HRFs calculated in this and other studies (Lu et al., 2006; Bagshaw et al., 2005), it is difficult to draw a distinction between an HRF that could be the response to a burst of pre-spike activity, and one that happens to have a particular shape, or to draw conclusions as to the specific effect of any early, unseen electrical activity on the BOLD signal,

1457

other than to say that early BOLD activity is present. Future studies would need to examine this question with more detailed analysis of the BOLD signal, possibly on an event by event basis to determine the consistency of the early responses and their relationship to the observed scalp spikes. The level of concordance between activations associated with early peaking HRFs and EEG localization is similar to that of later peaking activations and EEG localization. This suggests that the early peaking activations are not better than the later ones for localization purposes, which might have been suspected if it was considered that the early responses represented the initial epileptiform activity, while the later responses included propagated activity (in any event this is unlikely due to the fast propagation of epileptic activity, as discussed above). Patients with early peaking HRFs seem to have their maximal response to early model HRFs in the same way that other patients have a maximal response to later model HRFs, and the level of concordance between the EEG and fMRI data in patients with at least one early peaking calculated HRF and patients with no early peaking calculated HRF is relatively similar. This suggests that the early activations do not represent the origin of the epileptiform discharges any better than the later activations. In fact, it seems that these early activations are more or less equivalent to later activations in their quality and appearance. In general, concordance in our EEG–fMRI studies has been around 60% (Gotman et al., 2006), which is consistent with what is seen in this study. Analysis issues The good concordance between the model used in the statistical map and the peak timing of the calculated HRFs further demonstrates that a single model HRF is inadequate to detect all fMRI activations related to interictal spikes. This is consistent with the results of Bagshaw et al. (2004) and Lu et al. (2006). The patterns of activation seen in this study suggest that a range of HRFs is necessary. The optimal model HRFs should probably peak between − 1 and 7 s as little activity was seen with HRFs peaking at − 3 and 9 s. However, if the inclusion of HRFs peaking at these extremes serves to detect even a small number of activations, it may be worthwhile to consider their use. A model HRF peaking at 9 s has been demonstrated to be useful in detecting deactivations (Bagshaw et al., 2004, Aghakhani et al., 2004). The issue of HRF shape, in particular the FWHM, has not been evaluated systematically, but we saw in this study that some data sets (such as B, G, and G2) have responses that may have a larger FWHM than the model used to detect them. It is clear that the inclusion of the early peaking models (peak1, peak− 1, and peak− 3) was useful in detecting some activations that were not present in later peaks (see Fig. 3). In conclusion, we have established that in some cases there is a clear BOLD response starting several seconds before the short scalp EEG epileptic event used as a trigger to measure that change. This is a surprising finding, which may imply that some epileptic discharges are preceded by an energy-consuming phenomenon that lasts several seconds and may not consist of a synchronized neuronal discharge. References Aghakhani, Y., Bagshaw, A.P., Bénar, C.G., Hawco, C., Andermann, F., Dubeau, F., Gotman, J., 2004. fMRI activation during spike and wave discharges in idiopathic generalized epilepsy. Brain 127, 1127–1144.

1458

C.S. Hawco et al. / NeuroImage 35 (2007) 1450–1458

Aguirre, G.K., Zarahn, E., D'Esposito, M., 1998. The variability of human BOLD hemodynamic responses. NeuroImage 8, 360–369. Allen, P.J., Josephs, O., Turner, R., 2000. A method for removing imaging artifact from continuous EEG recorded during functional MRI. NeuroImage 12, 230–239. Bagshaw, A.P., Aghakhani, Y., Benar, C.G., Kobayashi, E., Hawco, C., Dubeau, F., Pike, G.B., Gotman, J., 2004. EEG–fMRI of focal epileptic spikes, analysis with multiple haemodynamic functions and comparison with gadolinium-enhanced MR angiograms. Hum. Brain Mapp. 22, 179–192. Bagshaw, A.P., Hawco, C., Benar, C.G., Kobayashi, E., Aghakhani, Y., Dubeau, F., Pike, G.B., Gotman, J., 2005. Analysis of the EEG–fMRI response to prolonged bursts of interictal epileptiform activity. NeuroImage 24, 1099–1112. Baumgartner, C., Serles, W., Leutmezer, F., Pataraia, E., Aull, S., Czeck, T., et al., 1998. Preictal SPECT in temporal lobe epilepsy, regional cerebral blood flow is increased prior to electroencephalography-seizure onset. J. Nucl. Med. 39, 978–982. Bénar, C.G., Gross, D.W., Wang, Y., Petre, V., Pike, B., Dubeau, F., Gotman, J., 2002. The BOLD response to interictal epileptiform discharges. NeuroImage 17, 1182–1192. Blumenfeld, H., 2005. Cellular and network mechanisms of spike-wave seizures. Epilepsia 46 (Suppl. 9), 21–33. Cox, R.W., 1996. AFNI software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res. 29, 162–173. Elger, C.E., Lehnertz, K., 1998. Seizure prediction by non-linear time series analysis of brain electrical activity. Eur. J. Neurosci. 10, 786–789. Emerson, R.G., Turner, C.A., Pedley, T.A., Walczak, T.S., Forgione, M., 1995. Propagation patterns of temporal spikes. Electroencephalogr. Clin. Neurophysiol. 94, 338–348. Ettinger, A.B., 1994. Structural causes of epilepsy. Tumors, cysts, stroke, and vascular malformations. Neurol. Clin. 121, 41–56. Federico, P., Abbott, D.F., Briellmann, R.S., Harvey, A.S., Jackson, G.D., 2005. Functional MRI of the pre-ictal state. Brain 128, 1811–1817. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. NeuroImage 9, 416–429. Gotman, J., Benar, C.G., Dubeau, F., 2004. Combining EEG and FMRI in epilepsy, methodological challenges and clinical results. J. Clin. Neurophysiol. 21, 229–240. Gotman, J., Kobayashi, E., Bagshaw, A.P., Benar, C.G., Dubeau, F., 2006. Combining EEG and fMRI, a multimodal tool for epilepsy research. J. Magn. Reson. Imaging 23, 906–920. Hamandi, K., Salek-Haddadi, A., Fish, D.R., Lemieux, L., 2004. EEG/ functional MRI in epilepsy, the Queen Square Experience. J. Clin. Neurophysiol. 21, 241–248. Handwerker, D., Ollinger, J., D'Esposito, M., 2004. Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses. NeuroImage 21, 1639–1651. Hauser, W.A., Annegers, J.F., Kurland, L.T., 1993. Incidence of epilepsy and unprovoked seizures in Rochester, Minnesota: 1935–1984. Epilepsia 34, 453–468. Kobayashi, E., Bagshaw, A.P., Grova, C., Dubeau, F., Gotman, J., 2006a. Negative BOLD responses to epileptic spikes. Hum. Brain Mapp. 27, 488–497.

Kobayashi, E., Hawco, C.S., Grova, C., Dubeau, F., Gotman, J., 2006b. Widespread and intense BOLD changes during brief focal electrographic seizures. Neurology 66, 1049–1055. Lemieux, L., Salek-Haddadi, A., Josephs, O., Allen, P., Toms, N., Scott, C., Krakow, K., Turner, R., Fish, D.R., 2001. Event-related fMRI with simultaneous and continuous EEG, description of the method and initial case report. NeuroImage 14, 780–787. Lu, Y., Bagshaw, A.P., Grova, C., Kobayashi, E., Dubeau, F., Gotman, J., 2006. Using voxel-specific hemodynamic response function in EEG– fMRI data analysis. NeuroImage 32, 238–247. Makiranta, M., Ruohonen, J., Suominen, K., Niinimaki, J., Sonkajarvi, E., Kiviniemi, V., Seppanen, T., Alahuhta, S., Jantti, V., Tervonen, O., 2005. BOLD signal increase precedes EEG spike activity, a dynamic penicillin induced focal epilepsy in deep anesthesia. NeuroImage 27, 715–724. Menz, M.M., Neumann, J., Muller, K., Zysset, S., 2006. Variability of the BOLD response over time: an examination of within-session differences. NeuroImage 32, 1185–1194. Miezin, F., Maccotta, L., Ollinger, J., Petersen, S., Buckner, R., 2000. Characterizing the hemodynamic response: effects of presentation rate, sampling procedure, and the possibility of ordering brain activity based on relative timing. NeuroImage 11, 735–759. Mottolese, C., Hermier, M., Stan, H., Jouvet, A., Saint-Pierre, G., Froment, J.C., Bret, P., Lapras, C., 2001. Central nervous system cavernomas in the pediatric age group. Neurosurg. Rev. 24, 55–73. Navarro, V., Martinerie, J., Le Van Quyen, M., Baulac, M., Dubeau, F., Gotman, J., 2005. Seizure anticipation, do mathematical measures correlate with video-EEG evaluation? Epilepsia 46, 385–396. Neumann, J., Lohmann, G., Zysset, S., von Cramon, D., 2003. Within-subject variability of BOLD response dynamics. NeuroImage 19, 784–796. Nunez, P.L., Silberstein, R.B., 2000. On the relationship of synaptic activity to macroscopic measurements, does co-registration of EEG with fMRI make sense? Brain Topogr. 13, 79–96. Pacia, S.V., Ebersole, J.S., 1997. Intracranial EEG substrates of scalp ictal patterns from temporal lobe foci. Epilepsia 38, 642–654. Shmuel, A., Yacoub, E., Pfeuffer, J., Van de Moortele, P.F., Adriany, G., Hu, X., Ugurbil, K., 2002. Sustained negative BOLD., blood flow and oxygen consumption response and its coupling to the positive response in the human brain. Neuron 36, 1195–1210. Shmuel, A., Augath, M., Oeltermann, A., Logothetis, N.K., 2006. Negative functional MRI response correlates with decreases in neuronal activity in monkey visual area V1. Nat. Neurosci. 9, 569–577. Stefanovic, B., Warnking, J.M., Kobayashi, E., Bagshaw, A.P., Hawco, C., Dubeau, F., Gotman, J., Pike, G.B., 2005. Hemodynamic and metabolic responses to activation, deactivation and epileptic discharges. NeuroImage 28, 205–215. Tao, J.X., Ray, A., Hawes-Ebersole, S., Ebersole, J.S., 2005. Intracranial EEG substrates of scalp EEG interictal spikes. Epilepsia 46, 669–676. Ward, B.D., 2000. Deconvolution analysis of FMRI time series data, AFNI, 3d deconvolve documentation. Medical College of Wisconsin. Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C., Friston, K.J., Evans, A.C., 1996. A unified statistical approach for determining significant signals in images of cerebral activation. Hum. Brain Mapp. 4, 58–73. Worsley, K.J., Liao, C.H., Aston, J., Petre, V., Duncan, G.H., Morales, F., Evans, A.C., 2002. A general statistical analysis for fMRI data. NeuroImage 15, 1–15.