Analysis of the EEG–fMRI response to prolonged bursts of interictal epileptiform activity

www.elsevier.com/locate/ynimg NeuroImage 24 (2005) 1099 – 1112 Analysis of the EEG–fMRI response to prolonged bursts of interictal epileptiform activity Andrew P. Bagshaw,* Colin Hawco, Christian-G. Benar, Eliane Kobayashi, Yahya Aghakhani, ´ Francois Dubeau, G. Bruce Pike, and Jean Gotman ¸ Department of Neurology and Neurosurgery, Montreal Neurological Institute, McGill University, Room 786, 3801 University Street, Montreal, ´ Quebec, Canada H3A 2B4 ´ Received 9 May 2004; revised 22 July 2004; accepted 7 October 2004 Available online 1 December 2004 The use of combined EEG–fMRI to study interictal epileptiform activity is increasing and has great potential as a clinical tool, but the haemodynamic response to epileptiform activity remains incompletely characterised. To this end, 19 data sets from 14 patients with prolonged bursts of focal or generalised interictal epileptiform activity lasting up to 15 s were analysed. To determine whether the inclusion of the durations of the epileptic events in the general linear model resulted in increased statistical significance of activated regions, statistical maps were generated with and without the event durations. The mean differences when including the durations were a 14.5% increase in peak t value and a 29.5% increase in volume of activation. This suggests that when analysing EEG–fMRI data from patients with prolonged bursts of interictal epileptiform activity, it is better to include the event durations. To determine whether the amplitudes and latencies of the measured responses were consistent with the general linear model, the haemodynamic response functions for bursts of different durations were calculated and compared with the model predictions. The measured amplitude of the response to the shortest duration events was consistently larger than predicted, which is consistent with studies in normal subjects. For the two data sets with the widest range of event durations, the measured amplitudes increased with the durations of the events without evidence of the plateau that was expected from the general linear model. There were no consistent differences between the measured and modelled latencies. D 2004 Elsevier Inc. All rights reserved. Keywords: EEG–fMRI; Interictal epileptiform activity; Haemodynamic response Introduction In the last few years, the simultaneous measurement of electroencephalography (EEG) and functional magnetic resonance * Corresponding author. Fax: +1 514 398 8106. E-mail address: andrew.bagshaw@mcgill.ca (A.P. Bagshaw). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2004.10.010 imaging (fMRI) has been used by several groups to study interictal activity in patients with epilepsy (Aghakhani et al., 2004; Al Asmi et al., 2003; Archer et al., 2003a; Boor et al., 2003; Jager et al., ¨ 2002; Krakow et al., 2001; Lazeyras et al., 2000; Lemieux et al., 2001). This work differs from most applications of fMRI in that the stimuli are internally generated, random, usually not very frequent, and of limited duration (from a fraction of a second to rarely more than 10 s). This makes the study of interictal epileptiform events extremely challenging in terms of the number and type of events that can be investigated in a particular patient. The long-term goal of such work is to use continuous EEG– fMRI recordings as one of a battery of clinical tools to locate the irritative zone, the area of tissue responsible for the generation of the interictal epileptiform activity (Rosenow and Luders, 2001), ¨ particularly in those patients who are candidates for surgical intervention. In order to achieve this goal, it will be necessary to understand the nature of the haemodynamic response to interictal epileptiform discharges in order to optimise the modelling of the fMRI data and hence increase the sensitivity of the method. To date, although the haemodynamic response to various types of stimuli has been reasonably well characterised in normal subjects, it has been much less studied in patient populations (Bagshaw et al., 2004; Barch et al., 2003; Benar et al., 2002; Buckner et al., ´ 2000; Renshaw et al., 1994). In addition, the extent to which the assumptions employed when analysing data from normal subjects are true in a patient population has not been assessed. In a standard experimental activation paradigm, the haemodynamic response to the stimulus is modelled by linear convolution of a time-invariant impulse response function (IRF) with a boxcar function of the same duration as the stimulus (Dale and Buckner, 1997; Miezin et al., 2000). A linear system is defined as one in which the scaling and superposition properties are satisfied, such that the response to long stimulus bursts can be predicted by shifting and adding appropriate bursts of shorter duration, or more precisely by convolving with a suitably scaled boxcar function. Departures from linearity in the blood oxygenation level dependent (BOLD) signal have been noted in normal subjects during a variety of tasks, and generally suggest that the response to stimuli of short 1100 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Table 1 Patients’ clinical details Patient A Age (years)/sex 22/F Diagnosis and clinical detailsa SGE Bi perisylvian polymicrogyria Seizures: head drops with CPS, L hand dystonic posturing, and GTCS R FLE Seizures: absences during infancy, but after age 10 seizures were characterised by head deviation to L, L arm elevation R TLE or FLE Seizures: staring, lip smacking, wringing of hands, GTCS L H focus with max T Nonlesional Seizures: sensation of shock, headache, then loses consciousness and falls IGE (JAE with EM) L F (SMA) L MCA congenital ischemia Seizures: blank episodes, tonic elevation of R arm and sometimes leg IGE (CAE) R T–O Nodular heterotopia Seizures: sensation of head spinning, then head turns to right, automatisms, moaning IGE (JAE) L TLE Multiple cavernous angiomas Seizures: CPS with automatisms, hears music, GTCS Partial Epilepsy (P T F involvement) Seizures: CPS leading to GTCS IGE (JAE) R FLE Seizures: focal seizures with secondary generalization, respiratory difficulty followed by head/eye deviation to the R IGE (GTCS) duration is larger than would be expected from analysing the response to stimuli of long duration (Birn et al., 2001; Boynton et al., 1996; Liu and Gao, 2000; Vasquez and Noll, 1998). It has also recently been suggested that, in generalised epilepsy with bursts of 30-Hz spike and wave activity, higher statistical values can be obtained by not including the duration of the discharge in the fMRI analysis, i.e., by taking the onset of the burst alone and modelling it as an instantaneous event (Ho et al., 2003). Given the limitations of working with stimuli that cannot be controlled experimentally, it is important that the analysis and modelling are as accurate as possible in order to maximise the sensitivity of the technique. To this end, the current study explored the issue of the optimum way to analyse prolonged bursts of epileptic activity. The first part of the study examined the effect on the statistical maps of including or ignoring the durations of the bursts of epileptic activity during the analysis. The second part compared the amplitude and latency (time to peak) of the measured responses to bursts of different duration with the expected responses from within the framework of the general linear model (GLM, Friston et al., 1995a,b; Worsley et al., 1996). The goals of the study were to characterise more thoroughly the haemodynamic response and to test the assumption of linearity in patients with epilepsy. B 24/M C 23/F D 35/M E F 24/F 19/F G H 18/F 25/F Materials and methods Patients Patients in whom prolonged bursts of interictal epileptiform activity were observed on the EEG recorded during fMRI scanning were selected from a database of 91 patients who had undergone continuous EEG–fMRI monitoring at the Montreal Neurological Institute and Hospital. As described below, the epileptiform activity observed in each patient was placed into temporal bins dependent upon the duration of the burst (i.e., 0–1 s, 2–3 s, etc). Since the purpose of the study was to compare the fMRI response to bursts of different duration, patients were excluded who did not have at least two bins with at least three events in each bin. In addition, one patient was excluded because the bursts were very long, often over 30 s in length, and may have represented electrographic seizures rather than interictal activity. Twenty-nine patients were initially selected due to the presence of bursts of epileptiform activity, and 14 were included in the study after applying the above inclusion criteria. Of these 14 patients, five were clinically diagnosed with idiopathic generalised epilepsy (IGE) and the remainder with partial epilepsy. Table 1 contains clinical details of the patients. Written informed consent was obtained from all patients in accordance with the regulations of the Research Ethics Board of the Montreal Neurological Institute and Hospital. Data acquisition Continuous EEG–fMRI data were acquired in all patients. Functional MRI images were acquired in one of two 1.5-T MR scanners (Vision and Sonata, Siemens, Erlangen, Germany) using an echo-planar imaging (EPI) sequence (voxel dimensions 5 Â 5 Â 5 mm, 25 slices, 64 Â 64 matrix, TE = 50 ms, TR = 3 s, flip angle 908). The fMRI data were collected in runs of 120 frames (one frame being a complete measurement of 25 slices) lasting I J 33/F 20/M K L M 25/M 22/M 29/F N a 37/F Abbreviations: R—right, L—left, F—frontal, T—temporal, P—parietal, H—hemisphere, Bi—bilateral, SGE—secondary generalised epilepsy, IGE—idiopathic generalised epilepsy, TLE—temporal lobe epilepsy, FLE—frontal lobe epilepsy, CPS—complex partial seizures, GTCS— generalised tonic-clonic seizures, JAE—juvenile absence epilepsy, CAE—childhood absence epilepsy, EM—eyelid myoclonus, SMA— supplementary motor area, MCA—middle cerebral artery. approximately 6 min. Between 7 and 12 runs were collected for each patient (see Table 3). A T1-weighted scan was also collected for anatomical localisation of the functional data sets (1-mm slice thickness, 256 Â 256 matrix, TE = 9.2 ms, TR = 22 ms, flip angle 308), leading to a total scanning time of approximately 2 h for each patient. EEG data were recorded via 21 MR compatible Ag/AgCl electrodes placed according to the International 10–20 system, using an EMR32 amplifier (Schwarzer, Munich, Germany). In order to minimise movement during the course of the scanning session, the patient’s head was immobilised with 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 artefact generated by the MR scanner A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 1101 (FEMR software, Schwarzer; Hoffmann et al., 2000). The spatial distribution, morphology, time of onset and duration of epileptic activity were subsequently marked by an experienced electroencephalographer (YA or EK). In all cases, the spatial distribution of the epileptic activity was stable without propagation during the course of a prolonged burst. Image preprocessing and statistical analysis The fMRI data were motion corrected and smoothed [Gaussian kernel 6-mm full width at half maximum (FWHM)] using in-house software. The first three frames of each run were not included in the analysis to ensure that the magnetisation was in a steady state, and differences in the slice acquisition time were corrected for. Models and signals were prewhitened with an autoregressive filter of order 1, and low frequency drifts in the signal were modelled with a third-order polynomial fitted to each run. Statistical analysis was performed using the methods and software of Worsley et al. (1996, 2002). Epileptic events were considered as different types based on the spatial distribution and morphology observed on the scalp EEG, resulting in 19 data sets from the 14 patients (i.e., five patients had more than one event type). In the following, the two data sets from, for example, patient C are designated as C1 and C2. Each data set was analysed with four monophasic, single gamma function IRFs peaking after 3, 5, 7 and 9 s (P3, P5, P7 and P9) to allow for some variation in the latency of the response while retaining information about its expected shape (Bagshaw et al., 2004; Buckner et al., 1998). The full width at half maximum of the gamma functions was 5.2 s. Composite statistical maps were created by taking the maximum value from the four analyses at each voxel, similar to the method employed by Saad et al. (2001) when looking at the delay of the response. Only positive responses were considered for further analysis. Comparison of composite maps with and without durations In order to determine whether the inclusion of the event durations improved the statistical analysis, two composite maps were created, one with EEG burst durations included in the model and one without. For each data set, depending on the extent of the activation, up to three separate clusters were selected and compared in the two composite maps. The clusters were selected in order of decreasing peak t value, starting with the cluster with the highest overall t value. For each cluster, the maximum t value and the volume of activation were calculated in the two t maps. The mean changes in t value and volume between the two maps were then calculated for all of the clusters in each patient and subsequently for all patients. In all cases, changes in t value and volume are expressed relative to the statistical maps generated without including the durations of the events. ROI selection and fitting of HRFs To compare the amplitude and latency of the measured data with the expected responses from the GLM, the measured haemodynamic response functions (HRFs) were estimated by fitting the fMRI signal within a region of interest (ROI) using a Fourier basis set (Josephs et al., 1997; Kang et al., 2003). The relative amplitudes and latencies of the fitted responses to events of different durations were then calculated and compared with the modelled responses. For each data set, an ROI was defined by selecting the voxel with the highest t value from the two composite maps. Four additional contiguous voxels above a t value of three were included to form a ROI of five voxels (625 mm3). Subsequent analysis of each data set thus concentrated on the data from these five voxels. For each patient, events with a common spatial distribution on the EEG were separated into time bins according to the duration of the burst. The bins were from 0 to 1 s, 2–3 s, 4–5 s, 6–7 s, 8–9 s, and greater than or equal to 10 s but less than 15 s. Only one patient had bursts that lasted for more than 15 s (patient H had eight events with durations from 19 to 50 s). These events were excluded from the analysis as there were few of them and they were spread over a wide range. Table 3 shows the number of events in each bin for each data set. In order to avoid a single abnormal event skewing the results, only bins with three or more events were further analysed. The fMRI response to each bin was assessed by fitting the fMRI data in a 64-s window surrounding each event (15 s prior to the event onset to 48 s after). For each data set, this resulted in a mean HRF for each bin. To avoid including data that did not have a robustly fitted response, amplitude and latency values were only calculated for HRFs with a signal to noise ratio (SNR) of greater than three. The SNR was calculated by taking the amplitude of the peak (in percentage signal change) and dividing it by the standard deviation of the background, defined as 10–5 s prior to and 30–48 s after the event. Comparison of fitted data and GLM A model of expected amplitude and latency for each bin was created by taking a standard IRF (Glover, 1999) and convolving it with a boxcar function of length equal to the median duration of each bin, using the software of Worsley et al. (2002). The model for each bin is shown in Fig. 1. The modelled amplitudes in Fig. 1 have been scaled so that the amplitude in the second bin is equal to 1. For comparison with the model, the amplitudes of the fitted HRFs were also scaled so that the fitted HRFs matched the modelled amplitude in the second temporal bin (i.e., the fitted amplitudes in the bin corresponding to bursts of 2–3 s were scaled to have a value of 1). A scaling factor was thus calculated based on the fitted data in the second bin, and subsequently applied to the amplitude values of Fig. 1. The haemodynamic responses predicted by the GLM for bursts of activity of different durations. The amplitudes have been scaled so that the amplitude of the response to bin 2–3 has a value of 1. 1102 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 the fitted responses for the other bins. This allows for a direct comparison between data sets from different patients of the relative changes in amplitude from one bin to the next. In a similar way, the latencies (times to peak) were scaled to match the time to peak of the second modelled bin. The absolute values of the latencies of the modelled data are dependent on the choice of IRF. In the present case, the time to peak of the modelled data in the second bin is close to 6.5 s. The second bin was selected as the reference bin as it contained a large number of events in most patients and lead to a robust fitted HRF. In one patient (data set H), the second bin was excluded due to a fitted HRF with a low SNR, so the third bin (event durations from 4 to 5 s) was used as the reference. Effect of choice of ROI on comparison of fitted data and GLM To assess the consistency of the results as a function of the choice of ROI, in eight data sets up to two additional ROIs of five voxels were selected. The eight data sets were chosen on the basis that the fitted HRFs from the original ROI were robust in several different bins. The additional ROIs were centred around the voxels with the next highest peak t values after the original ROI, provided that the peak voxels were separated by at least five intervening voxels. This was to ensure that the ROIs did not overlap. The amplitudes and latencies of the fitted HRFs from the additional ROIs were compared with the expected responses as for the original ROI. durations of 62.5% (data set C1) and a minimum of 0.4% (data set B1). The two data sets (C2 and E) which had higher t values when durations were ignored had mean differences of 7.4% and 7.5%, respectively. The volumes of activation did not necessarily increase when the peak t value increased. For example, for the data sets that had higher t values when durations were included, the mean change in cluster volume was an increase of 23.4%, but the range of values went from an increase in volume of 328.5% to a decrease in volume of 74.5%. Considering the generalised and focal epilepsy patients separately, the mean increase in t value for the five generalised epilepsy patients was 4.9% (range À7.5% to 18.1%) and 19.3% for the focal epilepsy patients (range À7.4% to 62.5%). The mean change in volume of activation was 16.9% (range À11.0% to 58.2%) for the generalised epilepsy patients and 35.7% (range À74.3% to 328.5%) for the focal epilepsy patients. These differences were not significant (two-tailed t test, P N 0.1 for both the change in t value and volume). Comparison of fitted data and GLM From the 19 data sets that had a sufficiently large number of epileptic events of different durations to be included in the study, five were subsequently excluded as a result of the SNR criterion applied to the fitted HRFs (data sets B1, H, J, M and N2). In addition, most other patients had at least one bin removed for the same reason (see Table 3). The bins that were kept for analysis and comparison with the model are marked in bold text in Table 3. Six data sets had only two bins remaining in the analysis following the application of the inclusion criteria. Fig. 3 has examples of the fitted HRFs that passed and failed the SNR criterion. The results of the comparison between the measured and modelled amplitudes of the BOLD response are shown in Fig. 4, with the analogous data for peak latencies in Fig. 5. For each patient, the scaled peak amplitude and latency are plotted with the corresponding data from the GLM, for each temporal bin. The method of fitting the data with a Fourier basis set allows an estimation of the error in the peak amplitude but does not lend itself to an estimation of the error in the calculation of peak latency. The error bars in Fig. 4 indicate plus and minus one standard deviation of the fitted amplitude (Worsley et al., 2002). Some general observations can be made concerning the agreement between the fitted and modelled amplitudes. The two data sets with the most bins (A with six and C1 with five) demonstrate similar behaviour in that the amplitude of the fitted HRFs continues to increase with the duration of the epileptic activity, and does not show the plateau that is evident for the modelled data (see Fig. 4). For eight data sets, the fitted HRFs from the shortest duration bin were of sufficient quality to be included in the study and for all of them the fitted amplitude of the response in the 0–1 second bin is larger than the expected modelled amplitude. In general, the comparison between the measured and modelled latencies of the responses demonstrates reasonably good agreement between the two (see Fig. 5). There are some exceptions that generally show that the measured latency did not increase as much as expected (see the results for data sets C1, I and L). In Fig. 6, the mean amplitudes and latencies measured for each bin and averaged across all data sets are plotted. In addition, the average values measured from patients with generalised (IGE) and partial epilepsy are plotted separately. When data were available from several ROIs within a particular data set (see the next section), these were also included in the calculation of the mean values. The Results Comparison of composite maps with and without durations Table 2 contains details of the comparison of the composite statistical maps generated with and without the durations of the events. For each patient, the change in peak t value and volume of activation are given for up to three separate clusters of activation. In four data sets (B2, H, M and N2), no clusters were detected that were significant and common to both composite statistical maps. This does not allow a meaningful numerical comparison of the effect of including the durations of the events. The results for these data sets are mixed, with B2 and H slightly better without durations, and M and N2 slightly better with. In the remaining 15 data sets with at least one cluster that was significant in both composite maps, 13 showed a mean increase in t value when durations were included and two showed a mean decrease. The mean change across all data sets was a 14.5% increase in t value and a 29.5% increase in volume when the durations of the events were included. Fig. 2 shows examples of the effects on the statistical maps of including or ignoring the durations of the events. In general, the areas of activation are similar in the two t maps. There is some suggestion of regional variability in the effect of including the durations of the events, with not all significant clusters affected in the same way (see, e.g., data set G in Fig. 2). This is also evident in the changes in t values and volumes that are given for each cluster in Table 2. However, overall the changes are relatively consistent, with either increasing or decreasing values across the clusters within a particular data set. For those data sets that had higher t values when durations were included, the mean increase in t value was 17.9% with a maximum difference between the composite maps including and ignoring A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Table 2 A direct comparison of up to three clusters per data set in activation maps created with and without the durations of the epileptic events Data set Cluster t value with durations 19.3 5.7 9.8 7.6 6.2 6.2 9 na na 17.9 8.2 6.3 5.2 18.9 9 5.9 4.6 4.3 6.1 5.7 5.6 10.8 6.5 5.6 11.1 na 4 na na 9.5 6.2 5.7 5.9 5.8 5.1 9.9 7.7 6.9 8.3 7.7 6.9 12.1 9.6 7.3 22.6 na 6.9 6 6.4 3.3 t value without durations 10.8 5.1 7 na 7.4 5.3 na 6.6 6.9 9.8 5.8 6.6 5.8 17.4 8.5 6.4 5.4 4.3 5.6 5.4 4.8 10.9 6 5.1 9.4 5.1 na 4.7 4.6 8.8 6.1 6.1 5.5 4.5 3.4 9.6 7.2 6.5 4.7 7.2 5.5 11.1 8 6.8 na 5.3 6.6 5.9 6.6 na Volume with durations (voxels) 397 9 85 36 62 16 13 na na 4584 39 3640 14 1134 7 2140 13 7 190 25 31 454 19 68 1963 na 5 na na 52 8 16 171 78 33 820 20 30 36 1485 402 149 200 12 32 na 83 29 17 5 Volume without durations (voxels) 2604 8 66 na 286 54 na 120 25 4207 40 3352 4 1350 13 1837 14 6 219 25 46 540 55 51 2206 13 na 11 13 44 8 13 124 16 5 786 18 28 22 986 373 393 52 23 na 8 72 25 15 na t value change per cluster (%) 78.7 11.8 40.0 À16.2 17.0 – – – 83.7 41.4 À4.5 À10.3 8.6 5.9 À7.8 À14.8 0.0 8.9 5.6 16.7 À0.9 8.3 9.8 18.1 – – – – 8.0 1.6 À6.6 7.3 28.9 50.0 3.1 6.9 6.2 76.6 6.9 25.5 9.0 20.0 7.4 – – 4.5 1.7 À3.0 – Volume change per cluster (%) À84.8 12.5 28.8 À78.3 À70.4 – – – 9.0 À2.5 8.6 250.0 À16.0 À46.2 16.5 À7.1 16.7 À13.2 0.0 À32.6 À15.9 À65.5 33.3 À11.0 – – – – 18.2 0.0 23.1 37.9 387.5 560.0 4.3 11.1 7.1 63.6 50.6 7.8 À62.1 284.6 À47.8 – – 15.3 16.0 13.3 – Mean t value change per data set (%) 1103 Mean volume change per data set (%) A B1 B2 C1 C2 D E F1 F2 G H I J K1 K2 L M N1 N2 a b c a b c a b c a b a b a b a b c a b c a b c a b a b c a b c a b c a b c a b c a b c a b a b c a 43.5 À14.5 0.4 À74.3 – 62.5 À7.4 7.3 – 3.2 129.3 À31.1 À7.5 8.7 10.4 À15.3 5.7 18.1 À16.0 À11.0 – – 1.0 13.8 28.7 328.5 5.4 7.5 36.3 40.7 12.1 – 58.2 – 1.1 – 14.9 – na = no activation. error bars indicate plus and minus 1 standard deviation of the fitted mean values. Effect of choice of ROI on comparison of fitted data and GLM In six of the eight data sets for which HRFs were fitted to data from three ROIs, all three resulted in fits that were of sufficient quality to be included. In one data set (A), only one additional ROI was included, and in another (G) no additional ROIs resulted in data of sufficient quality to be included. Figs. 7a and b contain the comparisons between the model and fitted data. Note that seven data sets are shown due to the absence of data set G. In general, the results of the comparison of the measured and modelled amplitudes demonstrate that within the uncer- 1104 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Fig. 2. Examples of statistical maps generated with and without the durations of the events included in the model. The maps for data sets A and D are clearly better when the durations are included, both in terms of the increased t values and volumes of activation. For data sets F2 and G, the differences are more subtle but there are slight increases in t value and volume when the event durations are ignored. tainty associated with the fitting process, the estimated values are consistent within a given data set. This is particularly evident in those data sets for which bin 0–1 is included, which show a consistent trend for the amplitude of bin 0–1 to be larger than expected. In addition, the long duration bins from data sets A and C1 demonstrate amplitudes that are consistently above the predictions of the model. The latencies of the fitted data from multiple ROIs are also consistent within a particular data set. With the single ROI, data set C1 had latency values that increased more slowly than the model, and this behaviour is replicated by the data from the two additional ROIs. Discussion Summary The results suggest that in general it is better to include the duration of the epileptic events in the statistical analysis. In the majority of the data sets analysed, the statistical significance of areas of activation increased when durations were included. In 15 data sets, several clusters of activation could be directly compared in statistical maps calculated with and without durations, and in 13 of them an increase in t value was noted when the durations were included. The changes in volumes of Table 3 Details of the scalp EEG inside the scanner and the number of events in each temporal bin Data set A B1 B2 C1 C2 D E F1 F2 G H I J K1 K2 L M N1 N2 a b Number of fMRI runs 10 10 10 11 11 7 10 12 12 10 9 10 11 7 7 10 8 10 10 Spike morphology and spatial distributiona Gen sharp and slow wave with L F–T predominance Bi F spikes R F–C sharp and slow wave Bi F–C–T sharp and slow wave R N L R F–T sharp and slow wave or slow L H spikes with F–T predominance Gen spike and slow wave with F predominance Bi F spikes L F spikes Gen spike and slow wave R T sharp and slow wave Gen spike and slow wave L T spikes L F–C–T spikes or slow wave Bi T sharp and slow wave with L predominance Gen spike and wave Bi F sharp and slow wave with shifting side predominance Gen spike and slow wave at times R N L R F–T sharp or slow sharp Bin (no. events)b 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 0–1 (100), 2–3 (61), 4–5 (21), 6–7 (5), 8–9 (4), 10+ (10) (43), 2–3 (8), 4–5 (2), 6–7 (0), 8–9 (1), 10+ (0) (5), 2–3 (3), 4–5 (6), 6–7 (0), 8–9 (1), 10+ (0) (55), 2–3 (44), 4–5 (24), 6–7 (15), 8–9 (6), 10+ (1) (22), 2–3 (28), 4–5 (8), 6–7 (1), 8–9 (0), 10+ (0) (7), 2–3 (7), 4–5 (10), 6–7 (6), 8–9 (2), 10+ (2) (14), 2–3 (20), 4–5 (2), 6–7 (0), 8–9 (0), 10+ (0) (0), 2–3 (5), 4–5 (4), 6–7 (0), 8–9 (0), 10+ (0) (3), 2–3 (22), 4–5 (6), 6–7 (0), 8–9 (0), 10+ (0) (4), 2–3 (19), 4–5 (15), 6–7 (4), 8–9 (1), 10+ (0) (16), 2–3 (17), 4–5 (13), 6–7 (3), 8–9 (3), 10+ (3) (11), 2–3 (9), 4–5 (4), 6–7 (1), 8–9 (5), 10+ (2) (476), 2–3 (76), 4–5 (28), 6–7 (3), 8–9 (1), 10+ (0) (166), 2–3 (6), 4–5 (0), 6–7 (0), 8–9 (0), 10+ (0) (66), 2–3 (49), 4–5 (13), 6–7 (2), 8–9 (3), 10+ (1) (15), 2–3 (14), 4–5 (25), 6–7 (2), 8–9 (0), 10+ (0) (35), 2–3 (41), 4–5 (18), 6–7 (10), 8–9 (0), 10+ (3) (20), 2–3 (24), 4–5 (1), 6–7 (0), 8–9 (0), 10+ (0) (5), 2–3 (4), 4–5 (0), 6–7 (0), 8–9 (0), 10+ (0) Abbreviations: R—right, L—left, F—frontal, C—Central, T—temporal, H—hemisphere, Bi—bilateral, Gen—generalised. Bins in bold text are those that were included in the analysis of the amplitudes and latencies of the HRFs. A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 1105 Fig. 3. Examples of the fitted HRFs. For data sets B2 and K2, bin 0–1 and bin 8–9, respectively, failed the inclusion criterion specifying that the SNR should be greater than three (see red arrows). The ROI from which the data were fitted is also given. activation were less consistent, but this is likely to be a result of activation regions sometimes appearing as large clusters and sometimes splitting into several smaller volumes. The mean magnitude of the increase in t value was 17.9%, compared with a decrease of 7.5% for the two data sets that were better without durations. This suggests that on an individual basis, it is likely to be more profitable to include the durations of the events in the analysis, since the possible gain in terms of statistical significance of the activations is considerably larger than the possible loss. This observation is true both for the patients with focal epilepsy and those with generalised epilepsy, with no statistically significant difference between the results for the two types of patients. This is in contrast to the abstract presented by Ho et al. (2003) in which a reduction in the statistical values was noted in analyses of data from five patients with generalised epilepsy when the event durations were included. The reason for this discrepancy is unclear, but may be due to differences in the methods of analysis, details of which were not clear from the abstract. Analysing the BOLD signal directly and comparing the fitted HRFs with those expected within the GLM leads to some evidence of nonlinearity, but in general is quite consistent with the observation that including the duration of the bursts is likely to lead to a better model of the BOLD signal change, and hence increased t values. There are no consistent differences between the measured and modelled latencies, although the measured values for some data sets do not increase as much as expected. The significance of this result is not clear. Some evidence of a disagreement between the data and the model can be observed in the amplitudes of the fits to data sets A and C1, which have six and five temporal bins, respectively. The measured amplitude increases monotonically with the duration of the epileptic events, and does not show evidence of the plateau predicted by the GLM. However, as discussed below, there are issues related to the form of the impulse response function used in the model, which mean that apparent discrepancies between the data and the model at long durations should be approached with caution. Another interesting observation is that, when the shortest duration bin is present, the mean of the fitted amplitude is always larger than predicted by the GLM. The modelled amplitude of the shortest duration bin is much less dependent on the choice of IRF than the amplitudes of the long duration bins (see below). In addition, this observation is consistent with previous studies of nonlinearity in various situations, which tend to show that the response to short duration stimuli is larger than would be predicted based on the response to long duration stimuli (Birn et al., 2001; Boynton et al., 1996; Liu and Gao, 2000; Vasquez and Noll, 1998). Methodological considerations There is an inherent problem in attempting to assess the linearity of the BOLD signal in a population of patients with 1106 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Fig. 4. Comparison of the amplitudes of the fitted and the predicted HRFs for each data set. The two data sets with more than four included bins have a linear trend line fitted to the amplitude values. The error bars indicate the uncertainty in the fitting procedure. epilepsy owing to the random and uncontrollable nature of the stimuli, and the fact that pathological events are generated by mechanisms that are not as well defined as those responsible for the functioning of the normal brain. One of the limitations of the current study is the number of events in each temporal bin, and the number of temporal bins that are available in each patient. As can be seen from Table 3, all of the patients had many more events in the short than in the long duration bins, and in six patients it was only possible to compare the haemodynamic responses in two of the six bins. There is no easy solution to this problem. Scanning the patients for longer would increase the number of events but is unlikely to be feasible, given that the scanning sessions already last for almost 2 h. Alternatively, it is possible that specific patients have unusually frequent bursts of several lengths, but this would A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 1107 call into question the extent to which the results could be applied to a general patient population. From an analytical point of view, using a more sophisticated method to estimate the form of the HRF, such as a Bayesian approach (Marrelec et al., 2003), might lead to more robustly fitted responses, which could result in fewer bins being rejected and hence allow the issue of linearity to be assessed more conclusively. The ROIs that were used for comparison with the predictions of the GLM were selected by performing an initial analysis with a constrained form for the IRF within a linear framework, which limits the type of responses that can be detected. Using several IRFs gives some flexibility in terms of the timing of the response and selecting the ROIs based on the voxels with the highest t values in either of the two composite maps will minimise the extent 1108 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Fig. 6. Mean amplitudes and latencies for each bin averaged across all data sets, and separately for generalised and focal epilepsy patients. The error bars indicate the spread in the mean fitted values across all data sets. (Only one data set had bin 10+ included, and so bin 10+ is not plotted here.) to which the results are biased towards voxels that clearly show linear behaviour. A more rigorous way of avoiding bias in the voxel selection would be to use an initial analysis that does not assume linearity or a specific form for the IRF, such as independent component analysis. However, this method remains under development and has not been shown to have the sensitivity necessary for the detection of brief and sporadic stimuli. Several methods have been used to compare measured data with the predictions of a linear convolution model, including the shifted summation of the response to short duration bursts to predict the response to longer bursts and the use of a short duration response as the IRF for subsequent linear convolution (Birn et al., 2001; Glover, 1999; Liu and Gao, 2000; Miller et al., 2001). More sophisticated approaches have also been used to estimate directly the nonlinear relationship between the stimulus presentation and the consequent haemodynamic response (Friston et al., 1998, 2000; Kershaw et al., 2001). The approach taken in the current study to model the expected responses was to use a standard IRF (Glover, 1999) and to convolve it with boxcar functions of appropriate lengths. The relative changes in amplitude and latency were then compared with the measured data. An alternative strategy would have been to determine a patient-specific IRF, either by using the measured response for the shortest duration bin or deconvolving the measured response to bursts of longer duration. The latter method has mainly been employed using the measured response to a block design, which leads to a robust estimate of the IRF. Fig. 7. Comparison of data and model from multiple ROIs within seven data sets, for (a) the amplitude and (b) the latency values. A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 1109 Fig. 7 (continued). However, in the present case, the fitted data were rather noisy, particularly for long-duration events, and would have resulted in a very noisy IRF. For this reason, it was felt that a standard response function would be more appropriate. The negative aspect of this choice is that the modelled values of amplitude and latency are dependent on the shape of the IRF. If the choice of IRF is inappropriate, the discrepancy apparent between the measured and modelled data could erroneously be attributed to nonlinearity. However, in the absence of a patient-specific IRF due to the insufficient amount of data, the standard response measured by Glover (1999) is a reasonable compromise, as it was efficient at detecting positive activation in patients with epilepsy (Bagshaw et al., 2004). Fig. 8 demonstrates the effect on the modelled amplitudes of changes in the IRF. The mean measured data across all data sets is also plotted. It can be seen that both removing the undershoot, and subsequently increasing the width of the first peak, delay the appearance of the plateau region, although the latter has more impact than the former and results in modelled data that more closely match the measured data. This suggests that caution should be exercised when interpreting the results of the comparison bet- Fig. 8. The effect of changing the shape of the modelled impulse response function on the predicted amplitudes and latencies as a function of duration. Results are given for the standard HRF (Model Glover), and for the standard HRF minus the undershoot (Model Glover No Undershoot) and with the FWHM subsequently doubled (Model Glover FWHM*2). Note that doubling the FWHM after having removed the undershoot does not change the latency values compared with only removing the undershoot, and so these data have not been plotted. The mean data across all data sets are also plotted for comparison. 1110 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 ween the measured data and the model for events with long durations. However, Fig. 8 shows that changes in the IRF do not have a significant effect on the amplitude or latency of the shortest duration bin, and hence the observation of the discrepancy between model and data remains. It should also be noted that an IRF with a width of 10.4 s (i.e., twice the FWHM of the standard IRF) is extremely unlikely physiologically (see, e.g., Aguirre et al., 1998), suggesting that not all of the discrepancy between the modelled and fitted data can be attributed to the inadequacies of the IRF model. Physiological considerations Data from patients with focal and generalised epilepsy were analysed, provided that prolonged bursts of epileptiform activity were observed on the scalp EEG. Although this raises the question of whether the mechanisms producing the epileptic activity are different in focal and generalised epilepsy, and hence whether different types of experimental stimuli are being considered, it does not alter the conclusions that can be drawn concerning the BOLD response for a particular patient. It would seem to be a reasonable assumption that for a given patient the epileptic activity is produced by the same mechanism, independently of the duration of the burst. Under this assumption, therefore, the bstimulusQ producing the BOLD response for a given event type is constant, and the only variable parameter is the length of the stimulus. However, the mechanisms producing runs of spikes and bursts of spike and wave are not completely understood and it is possible that the amplitude of the neural event (stimulus) producing a burst of prolonged activity is greater than that responsible for an analogous but shorter event, in which case two parameters vary simultaneously. Most neuroimaging techniques have been used to investigate the haemodynamic/metabolic consequences of prolonged epileptic activity in humans, including PET, SPECT and 133Xe inhalation (Abou-Khalil et al., 1987; Bajc et al., 1992; Bittar et al., 1999; Chugani et al., 1993; Engel et al., 1985; Ochs et al., 1987; Prevett et al., 1995; Sikai et al., 1978; Sperling and Skolnick, 1995; Theodore et al., 1985). However, all have a low temporal resolution, with an acquisition time of at least several minutes, meaning that interictal states usually comprise both normal and abnormal EEG. Some recent studies using EEG–fMRI have investigated patients with spike and wave discharges, without specifically examining the effect of events of different durations (Aghakhani et al., 2004; Archer et al., 2003b; Hill et al., 1999; Salek-Haddadi et al., 2003). Transcranial Doppler ultrasound can measure the changes in blood flow velocities as a result of epileptic activity, as shown by Diehl et al. (1998) who observed changes in the middle cerebral arteries in 13 patients. They demonstrated that bursts of generalised spike and wave with a mean duration of 8.4 F 7.2 s (range 1.5–30 s) lead to cortical perfusion changes, but did not investigate the relative changes caused by bursts of different durations, other than to note that short duration bursts (mean duration 0.92 F 0.2 s, range 0.5–1.2 s) did not cause significant flow changes. This may explain the fact that bin 0–1 had to be excluded from the analysis on many occasions due to a fitted HRF with a low SNR, despite the fact that there were often many events. The mechanisms responsible for the generation of generalised spike and wave discharges and runs of focal spikes are different. The most commonly accepted theory of spike and wave generation involves an interplay between the thalamus and the cortex by which normal thalamic discharges are sent to a slightly hyperexcitable cortex that subsequently responds with spike and wave activity (Gloor, 1969; Jasper and Droogleever-Fortuyn, 1947). A recent EEG–fMRI study has shown that generalised spike and wave discharges in humans are usually associated with thalamic activation and both cortical activation and deactivation (Aghakhani et al., 2004), a pattern very similar to that observed in an experimental model of absence seizures (Tenney et al., 2003). Focal epileptic discharges, on the other hand, are usually considered local phenomena, although the current EEG–fMRI investigations tend also to implicate regions remote from the apparent focus. Combined EEG–fMRI is the only currently available neuroimaging approach with a sufficiently high temporal resolution to clearly differentiate the haemodynamic response to different kinds of epileptic spikes in human subjects, and hence is potentially useful not only as a clinical tool but also as a method of understanding the basic mechanisms responsible for epileptic discharges. The evidence available from the current group of patients suggests that the BOLD response to both focal and generalised epileptic activity behaves in much the same way (see Fig. 6), despite the different mechanisms that are undoubtedly at work. This is an interesting observation, and one that may warrant further investigation. Future work Many questions remain concerning the nature of the haemodynamic response to prolonged bursts of epileptiform activity. In the current work, increases in the peak t value and extent of activation were considered as indicative of a more accurate description of the region of cortex activated by the epileptic activity. However, it is known that gradient echo EPI sequences are more sensitive to BOLD signal changes in large veins than the capillary bed (Lai et al., 1993; Lee et al., 1995), and hence as with most fMRI studies the activation being detected may include a substantial contribution from veins at some distance from the site of neuronal activity (Turner, 2002). This highlights the need for fMRI acquisition protocols that are less sensitive to venous signal changes and detailed comparison with invasive or noninvasive neurophysiological recordings. This study concentrated on positive BOLD responses, but negative responses are also widely seen in patients with epilepsy (Archer et al., 2003b; Bagshaw et al., 2004). The interpretation of negative responses is less clear than for positive, and it would be interesting to examine the effect of the duration of the bursts on deactivations, although this may require better fitting techniques and more data, since negative responses tend to be of lower amplitude (Shmuel et al., 2002). The selection of voxels to form the ROIs was based purely on their statistical significance, and no effects associated with their location were examined. Regional differences in the properties of the haemodynamic response have been noted, for example between motor and visual cortex (Birn et al., 2001; Miller et al., 2001; Waldvogel et al., 2000) and calcarine and fusiform cortex (Huettel and McCarthy, 2001). A similar type of study could be undertaken in patients with epilepsy, particularly since widespread activation was often seen in the present group of patients associated with the often widespread EEG abnormalities. One problem is that the interpretation of the different regions of activation is much less clear than in, for example, visual and motor tasks. As more patients are scanned, it will be possible to determine whether there are consistent differences between the haemodynamic responses in A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 1111 patients with different types of focal and generalised epilepsy, including the effect of cortical lesions and modifications to the vasculature as a result of surgery. Conclusion This study characterises the haemodynamic response to prolonged bursts of interictal epileptiform activity, and represents a more thorough investigation of this issue than has been performed previously. The results suggest that when analysing EEG–fMRI data from patients with prolonged bursts of interictal epileptiform activity, it is better to include the durations of the events in the statistical analysis. In general, the comparison between the measured amplitudes and latencies and those expected from the GLM shows reasonable agreement, although there is evidence of disagreement for the shortest and longest duration bursts. Acknowledgments This work was supported by grant MOP 38079 of the Canadian Institutes of Health Research. The authors would like to thank B. Stefanovic for helpful comments. CGB was funded by a CIHR doctoral research award. EK was funded by a Preston Robb Fellowship from the Montreal Neurological Institute. References Abou-Khalil, W.B., Siegel, G.J., Sackellares, J.C., Gilman, S., Hichwa, R., Marshall, R., 1987. Positron emission tomography studies of cerebral glucose metabolism in chronic partial epilepsy. Ann. Neurol. 22, 480 – 486. Aghakhani, Y., Bagshaw, A.P., Benar, 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. Aguirre, G.K., Zarahn, E., D’Esposito, M., 1998. The variability of human, BOLD hemodynamic responses. NeuroImage 8, 360 – 369. Al Asmi, A., Benar, C.-G., Gross, D.W., Agha Khani, Y., Andermann, F., ´ Pike, B., Dubeau, F., Gotman, J., 2003. fMRI activation in continuous and spike-triggered EEG–fMRI studies of epileptic spikes. Epilepsia 44 (10), 1328 – 1339. Archer, J.S., Briellmann, R.S., Syngeniotis, A., Abbott, D.F., Jackson, G.D., 2003a. Spike-triggered fMRI in reading epilepsy. Neurology 60, 415 – 421. Archer, J.S., Abbott, D.F., Waites, A.B., Jackson, G.D., 2003b. fMRI ddeactivationT of the posterior cingulate during generalized spike and wave. NeuroImage 20, 1915 – 1922. 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. ˇ´ ˇ Bajc, M., Basia, M., Hajnsek, S., Topuzovia, N., 1992. Regional ´ brina perfusion changes in patients with primary generalized epilepsy assessed by Tc99m HM-PAO and SPECT. Neurol. Croat. 41, 13 – 20. Barch, D.M., Mathews, J.R., Buckner, R.L., Maccotta, L., Csernansky, J.G., Snyder, A.Z., 2003. Hemodynamic responses in visual, motor and somatosensory cortices in schizophrenia. NeuroImage 20, 1884 – 1893. Benar, 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. Birn, R.M., Saad, Z.S., Bandettini, P.A., 2001. Spatial heterogeneity of the nonlinear dynamics in the fMRI BOLD response. NeuroImage 14, 817 – 826. Bittar, R.G., Andermann, F., Olivier, A., Dubeau, F., Dumoulin, S.O., Pike, G.B., Reutens, D.C., 1999. Interictal spikes increase cerebral glucose metabolism and blood flow: a PET study. Epilepsia 40 (2), 170 – 178. Boor, S., Vucurevic, G., Pleiderer, C., Stoeter, P., Kutschke, G., Boor, R., 2003. EEG-related functional MRI in benign childhood epilepsy with centrotemporal spikes. Epilepsia 44 (5), 688 – 692. Boynton, G.M., Engel, S.A., Glover, G.H., Heeger, D.J., 1996. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16 (13), 4207 – 4221. Buckner, R.L., Koutstaal, W., Schacter, D.L., Dale, A.M., Rotte, M., Rosen, B.R., 1998. Functional–anatomic study of episodic retrieval II Selective averaging of event-related fMRI trials to test the retrieval success hypothesis. NeuroImage 7, 163 – 175. Buckner, R.L., Snyder, A.Z., Sanders, A.L., Raichle, M.E., Morris, J.C., 2000. Functional brain imaging of young, nondemented and demented older adults. J. Cogn. Neurosci. 12 (Suppl 2), 24 – 34. Chugani, H.T., Shewmon, A., Khanna, S., Phelps, M.E., 1993. Interictal and postictal focal hypermetabolism on positron emission tomography. Pediatr. Neurol. 9, 10 – 15. Dale, A.M., Buckner, R.L., 1997. Selective averaging of rapidly presented individual trials using fMRI. Hum. Brain Mapp. 5, 329 – 340. Diehl, B., Knecht, S., Deppe, M., Young, C., Stodieck, S.R.G., 1998. Cerebral hemodynamic response to generalized spike-wave discharges. Epilepsia 39 (12), 1284 – 1289. Engel, J., Lubens, P., Kuhl, D.E., Phelps, M.E., 1985. Local cerebral metabolic rate for glucose during petit mal absences. Ann. Neurol. 17, 121 – 128. Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J.B., Frith, C.D., Frackowiak, R.S.J., 1995. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapp. 2, 189 – 210. Friston, K.J., Holmes, A.P., Poline, J.B., Grasbky, P.J., Williams, S.C., Frackowiak, R.S., Turner, R., 1995. Analysis of fMRI time-series revisited. NeuroImage 2, 45 – 53. Friston, K.J., Josephs, O., Rees, G., Turner, R., 1998. Nonlinear eventrelated responses in fMRI. Magn. Reson. Med. 39, 41 – 52. Friston, K.J., Mechelli, A., Turner, R., Price, C.J., 2000. Nonlinear responses in fMRI: the balloon model, Volterra kernels and other hemodynamics. NeuroImage 12, 466 – 477. Gloor, P., 1969. Neurophysiological basis of generalized seizures termed centrocephalic. In: Gastaut, H., Jasper, H., Bancaud, J., Waltregny, A. (Eds.), The Physiopathogenesis of the Epilepsies. Thomas, Springfield, pp. 209 – 236. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. NeuroImage 9, 416 – 429. Hill, R.A., Chiappa, K.H., Huang-Hellinger, F., Jenkins, B.G., 1999. Hemodynamic and metabolic aspects of photosensitive epilepsy revealed by functional magnetic resonance imaging and magnetic resonance spectroscopy. Epilepsia 4 (7), 912 – 920. Ho, L., Lemieux, L., Salek-Haddadi, A., Diehl, B., 2003. The relationship of the BOLD response to the temporal distribution of epileptiform discharges in EEG–fMRI [abstract]. Presented at the 9th International Conference on Functional Mapping of the Human Brain, June 19–22, 2003, New York, NY. Available on CD-Rom in NeuroImage 19(2). Hoffmann, A., J7ger, L., Werhahn, K.J., Jaschke, M., Noachter, S., Reiser, M., 2000. Electroencephalography during functional echo-planar imaging: detection of epileptic spikes using post-processing methods. Magn. Reson. Med. 44, 791 – 798. Huettel, S.A., McCarthy, G., 2001. Regional differences in the refractory period of the hemodynamic response: an event-related fMRI study. NeuroImage 14, 967 – 976. J7ger, L., Werhahn, K.J., Hoffmann, A., Berthold, S., Scholz, V., Weber, J., Noachtar, S., Reiser, M., 2002. Focal epileptiform activity in the brain: 1112 A.P. Bagshaw et al. / NeuroImage 24 (2005) 1099–1112 Frank, L.R., Buxton, R.B., 2001. Nonlinear temporal distribution of the cerebral blood flow response. Hum. Brain Mapp. 13, 1 – 12. Ochs, R.F., Gloor, P., Tyler, J.L., Wolfson, T., Worsley, K., Andermann, F., Diksic, M., Meyer, E., Evans, A., 1987. Effect of generalized spike-andwave discharge on glucose metabolism measured by positron emission tomography. Ann. Neurol. 21, 458 – 464. Prevett, M.C., Duncan, J.S., Jones, T., Fish, D.R., Brooks, D.J., 1995. Demonstration of thalamic activation during typical absence seizures using H215O and PET. Neurology 45, 1396 – 1402. Renshaw, P.F., Yurgelun-Todd, D.A., Cohen, B.M., 1994. Greater hemodynamic response to photic stimulation in schizophrenic patients: an echo planar MRI study. Am. J. Psychiatry 151, 1493 – 1495. Rosenow, F., Lqders, H., 2001. Presurgical evaluation of epilepsy. Brain 124, 1683 – 1700. Saad, Z.S., Ropella, K.M., Cox, R.W., DeYoe, E.A., 2001. Analysis and use of fMRI response delays. Hum. Brain Mapp. 13, 74 – 93. Salek-Haddadi, A., Lemieux, L., Merschhemke, M., Friston, K.J., Duncan, J.S., Fish, D.R., 2003. Functional magnetic resonance imaging of human absence seizures. Ann. Neurol. 53, 663 – 667. 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. Sikai, F., Meyer, J.S., Narritomi, H., Hsu, M.-C., 1978. Regional cerebral blood flow in patients with epilepsy. Arch. Neurol. 35, 648 – 657. Sperling, M.R., Skolnick, B.E., 1995. Cerebral blood flow during spikewave discharges. Epilepsia 36 (2), 156 – 163. Tenney, J.R., Duong, T.Q., King, J.A., Ludwig, R., Ferris, C.F., 2003. Corticothalamic modulation during absence seizures in rats: a functional MRI assessment. Epilepsia 44 (9), 1133 – 1140. Theodore, W.H., Brooks, R., Margolin, R., Patronas, N., Sato, S., Porter, R.J., Mansi, L., Bairamian, D., DiChiro, G., 1985. Positron emission tomography in generalized seizures. Neurology 35, 684 – 690. Turner, R., 2002. How much cortex can a vein drain? Downstream dilution of activation-related cerebral blood oxygenation changes. NeuroImage 16, 1062 – 1067. Vasquez, A.L., Noll, D.C., 1998. Nonlinear aspects of the BOLD response in functional MRI. NeuroImage 7, 108 – 118. Waldvogel, D., van Gelderen, P., Immisch, I., Pfeiffer, C., Hallett, M., 2000. The variability of serial fMRI data: correlation between a visual and motor task. NeuroReport 11, 3843 – 3847. 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. detection with spike-related functional MR imaging–preliminary results. Radiology 223, 860 – 869. Jasper, H.H., Droogleever-Fortuyn, J., 1947. Experimental studies on the functional anatomy of petit mal epilepsy. Res. Publ.-Assoc. Res. Nerv. Dis. 26, 272 – 298. Josephs, O., Turner, R., Friston, K., 1997. Event-related fMRI. Hum. Brain Mapp. 5, 243 – 248. Kang, J.K., Benar, C.-G., Al-Asmi, A., Agha Khani, Y., Pike, G.B., ´ Dubeau, F., Gotman, J., 2003. Using patient-specific hemodynamic response functions in combined EEG–fMRI studies in epilepsy. NeuroImage 20, 1162 – 1170. Kershaw, J., Kashikura, K., Zhang, X., Abe, S., Kanno, I., 2001. Bayesian technique for investigating linearity in event-related BOLD fMRI. Magn. Reson. Med. 45, 1081 – 1094. Krakow, K., Lemieux, L., Messina, D., Scott, C.A., Symms, M.R., Duncan, J.S., Fish, D.R., 2001. Spatio-temporal imaging of focal interictal epileptiform activity using EEG-triggered functional MRI. Epileptic Disord. 3 (2), 67 – 74. Lai, S., Hopkins, A.L., Haacke, E.M., Li, D., Wasserman, B.A., Buckley, P., Friedman, L., Meltzer, H., Hedera, P., Friedland, R., 1993. Identification of vascular structures as a major source of signal contrast in high resolution 2D and 3D functional activation imaging of the motor cortex at 1.5T: preliminary results. Magn. Reson. Med. 30, 387 – 392. Lazeyras, F., Blanke, O., Perrig, S., Zimine, I., Golay, X., Delavelle, J., Michel, C.M., de Tribolet, N., Villemure, J.-G., Seeck, M., 2000. EEGtriggered functional MRI in patients with pharmacoresistant epilepsy. J. Magn. Reson. Imaging 12, 177 – 185. Lee, A.T., Glover, G.H., Meyer, C.H., 1995. Discrimination of large venous vessels in time course spiral blood-oxygen-level-dependent magnetic resonance neuroimaging. Magn. Reson. Med. 33, 745 – 754. 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. Liu, H.-L., Gao, J.-H., 2000. An investigation of the impulse response functions for the nonlinear BOLD response in functional MRI. Magn. Reson. Imaging 18, 931 – 938. Marrelec, G., Benali, H., Ciuciu, P., Pelegrini-Isaac, M., Poline, J.-B., 2003. ´´ Robust Bayesian estimation of the hemodynamic response function in event-related fMRI using basic physiological information. Hum. Brain Mapp. 19, 1 – 17. Miezin, F.M., Maccotta, L., Ollinger, J.M., Petersen, S.E., Buckner, R.L., 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. Miller, K.L., Luh, W.-M., Liu, T.T., Martinez, A., Obata, T., Wong, E.C.,