BOLD mapping of human epileptic spikes recorded during simultaneous intracranial EEG-fMRI: The impact of automated spike classification

Objectives: Simultaneous intracranial EEG and functional MRI (icEEG-fMRI) can be used to map the haemodynamic (BOLD) changes associated with the generation of IEDs. Unlike scalp EEG-fMRI, in most patients who undergo icEEG-fMRI, IEDs recorded intracranially are numerous and show variability in terms of field amplitude and morphology. Therefore, visual marking can be highly subjective and time consuming. In this study, we applied an automated spike classification algorithm, Wave_clus (WC), to IEDs marked visually on icEEG data acquired during simultaneous fMRI acquisition. The motivation of this work is to determine whether using a potentially more consistent and unbiased automated approach can produce more biologically meaningful BOLD patterns compared to the BOLD patterns obtained based on the conventional, visual classification. Methods: We analysed simultaneous icEEG-fMRI data from eight patients with severe drug resistant epilepsy, and who subsequently underwent resective surgery that resulted in a good outcome: confirmed epileptogenic zone (EZ). For each patient two fMRI analyses were performed: one based on the conventional visual IED classification and the other based on the automated classification. We used the concordance of the IED-related BOLD maps with the confirmed EZ as an indication of their biological meaning, which we compared for the automated and visual classifications for all IED originating in the EZ. Results: Across the group, the visual and automated classifications resulted in 32 and 24 EZ IED classes respectively, for which 75% vs 83% of the corresponding BOLD maps were concordant. At the single-subject level, the BOLD maps for the automated approach had greater concordance in four patients, and less concordance in one patient, compared to those obtained using the conventional visual classification, and equal concordance for three remaining patients. These differences did not reach statistical significance. Conclusion: We found automated IED classification on icEEG data recorded during fMRI to be feasible and to result in IED-related BOLD maps that may contain similar or greater biological meaning compared to the conventional approach in the majority of the cases studied. We anticipate epileptiform detection and classification of IEDs on icEEG by clinical neurophysiologists can be consuming and is unreliable and Elger, 1999, Barkmeier et al., 2012, Gaspard et al., 2014, Sharma et al., 2017). Furthermore, the incorrect and inconsistent markings of IEDs have been shown to result in an excess of false positive and false negative BOLD clusters (Flanagan et al., 2009). correlates additional, independent IED classifications’ biological value. To investigate we analysed icEEG-fMRI data collected from patients who subsequently underwent resective surgery that resulted in a good postsurgical outcome, thereby providing some confirmation of the epileptogenic zone (EZ). We compared the BOLD maps obtained using two general linear models: one derived from the IEDs classified automatically and the other derived from the IEDs classified using the conventional visual method; this comparison was made in terms of localisation of the IED-related BOLD changesin relation tothe confirmed EZ. We hypothesised that an automated approach to IED classification can result in more or equally concordant, and therefore biologically meaningful, BOLD maps.


M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT

Introduction
Simultaneous EEG-fMRI has been shown to be a useful tool in mapping the regions associated with the generation of interictal epileptiform discharges (IEDs). For example, blood-oxygen level dependent (BOLD) mapping of IEDs detected on scalp EEG can provide added value to the localisation of the EZ (Zijlmans et al., 2007, Moeller et al., 2009, Thornton et al., 2010, Thornton et al., 2011, Pittau et al., 2012, An et al., 2013, Coan et al., 2016, Centeno et al., 2017. In one study, patients were reconsidered for surgery after identifying the BOLD correlates of IEDs (Zijlmans et al., 2007) and others have suggested its potential use in predicting postsurgical outcome (Thornton et al., 2010, An et al., 2013, Coan et al., 2016, Centeno et al., 2017. In these studies, IEDrelated BOLD changes have been found to be located in proximity, but also often remote from the EZ suggesting that their generators can involve a widespread network. This may partly reflect the area of synchronous neural activity required for their detection on scalp EEG (Tao et al., 2005). However, the sensitivity of detecting IEDs during an fMRI scan can be a limitation for scalp EEG-fMRI studies (Aghakhani et al., 2006, Salek-Haddadi et al., 2006. In contrast, the sensitivity of intracranial EEG (icEEG) to IEDs is much greater, due to the lack of attenuation and spatial integration from the scalp and the skull (Carreño and Lüders, 2001), which has led to the implementation of simultaneous icEEG and fMRI acquisitions (icEEG-fMRI) in an effort to better understand the haemodynamic changes associated with epileptic activity (Carmichael et al., 2012, Boucousis et al., 2012. To date, there have only been three studies investigating the whole-brain haemodynamic correlates of IEDs detected on icEEG (Vulliemoz et al., 2011, Aghakhani et al., 2015. In these studies, the IEDs were detected and classified by an EEG reviewer in the usual visual manner and a general linear model (GLM) was used to map the BOLD correlates of each IED class. From a theoretical viewpoint, in order to create the optimal model of the IED-related BOLD changes, it is important that once the IEDs have been detected, that they be classified so that each regressor represents the activity of a specific neuronal population to the exclusion of other generators (Liston et al., 2006). In the conventional approach the classification of IEDs generally requires the reviewer to distinguish different IEDs based on their field distribution: essentially the EEG channels in which they occur (Gotman, 1999, James et al., 1999. However as noted above, the high sensitivity of icEEG can result in much more abundant IEDs than in scalp recordings, and often more varied in their morphology and distribution than on scalp EEG (Spencer et al., 2015). Indeed, the M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT detection and classification of IEDs on icEEG by clinical neurophysiologists can be time consuming and is unreliable (Dümpelmann and Elger, 1999, Barkmeier et al., 2012, Gaspard et al., 2014, Sharma et al., 2017. Furthermore, the incorrect and inconsistent markings of IEDs have been shown to result in an excess of false positive and false negative BOLD clusters (Flanagan et al., 2009).
Automated algorithms for the analysis of IEDs on icEEG have been designed with the principal aim of reducing the subjectivity of EEG marking (Dümpelmann and Elger, 1999, Bourien et al., 2005, Valenti et al., 2006, Brown et al., 2007, Barkmeier et al., 2012, Gaspard et al., 2014. Such methods are able to detect IEDs but do not exploit the relationship between the activity across channels. Given that IED field distribution is an important consideration in the way IEDs are typically classified visually, we favour automated schemes that incorporate the waveforms across several channels (Hufnagel et al., 2000, Bourien et al., 2005, Pedreira et al., 2014, Sharma et al., 2017. To our knowledge, the only such classification scheme evaluated by comparison with multiple observers was Sharma et al. (2017). This study validated an automated neuronal spike classification algorithm, Wave_clus (WC) (Quian Quiroga et al., 2004, Pedreira et al., 2014 for the automated classification of icEEG IEDs (Sharma et al., 2017), by formally demonstrating that the WC results fall within the observed inter-rater variability for three expert EEG reviewers.
In the present work, our aim was to use this automated IED classification algorithm (WC) (Sharma et al., 2017) to determine whether the models of the IED-related BOLD changes obtained based on the resulting IED classification are comparable to those obtained using the traditional, visual approach. This follows our previous work on the evaluation of the value of the automated classification of IED on scalp EEG-fMRI based on their BOLD correlates (Pedreira et al, 2014);In essence we use the fMRI dataas an additional, independent means of assessing the IED classifications' biological value. To investigate this, we analysed icEEG-fMRI data collected from patients who subsequently underwent resective surgery that resulted in a good postsurgical outcome, thereby providing some confirmation of the epileptogenic zone (EZ). We compared the BOLD maps obtained using two general linear models: one derived from the IEDs classified automatically and the other derived from the IEDs classified using the conventional visual method; this comparison was made in terms of localisation of the IED-related BOLD changesin relation tothe confirmed EZ. We hypothesised that an automated approach to IED classification can result in more or equally concordant, and therefore biologically meaningful, BOLD maps.

Patients
Eight patients (6 males; age range: 32-42 years) who underwent simultaneous intracranial EEG-fMRI at the National Hospital for Neurology and Neurosurgery (UCLH NHS Foundation Trust, Queen Square, London, UK) were included in this study (Table   1). These patients were part of a group of 19 who underwent icEEG-fMRI; we retrospectively selected all those who subsequently underwent resective surgery with a good postsurgical outcome (ILAE 1;Wieser et al, 2001) (no seizure for at least 3 years post-resection). The patients gave informed written consent for participation in this study, which was approved by the joint research ethics committee of the National Hospital for Neurology and Neurosurgery and UCL Institute of Neurology, Queen Square, London, UK.
At the time of participation in this study all patients were undergoing invasive icEEG recordings for clinical purposes to assess their suitability for curative surgery. The need for invasive EEG recordings had been established in a multidisciplinary meeting to delineate the ictal onset zone and/or to perform direct electro-cortical stimulation for functional mapping. Prior to the invasive EEG recordings, the patients had undergone standard pre-surgical evaluation including long-term video-scalp EEG monitoring, structural MRI and other investigations such as positron emission tomography (PET), magnetoencephalography (MEG) or ictal single photon emission computed tomography (ictal SPECT) when available.

Simultaneous icEEG-fMRI acquisition
Between 31 and 91 electrode contacts were implanted in each patient with grid, strip and depth electrodes or a combination thereof (all Ad-Tech subdural epilepsy or Spencer Probe depth electrodes; Ad-Tech Medical Instrument Corporation, Racine, WI, USA); see Table 2. The electrodes were connected to an MR-compatible EEG amplifier system (BrainAmp MR+; Brain Products, Gilching, Germany). In accordance with our icEEG-fMRI protocol (Carmichael et al., 2012), echo planar images (EPI) were acquired using a 1.5T Siemens Avanto scanner (Erlangen, Germany) with a standard birdcage transmit/receive head-coil. Depending upon patient comfort inside the scanner and time constraints either one (for patients 1, 2, 3, 4) or two (patients 5, 6, 7, 8) useable 10minute resting-state EPI time series was obtained; a T1-weighted structural scan was also acquired. The icEEG signals were recorded at a sampling rate of 5 kHz. In the case of M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT patient 3 the scanning-related artefacts on icEEG during the first EPI sequence could not be corrected satisfactorily due to a technical problem at acquisition; therefore, only the second EPI sequence was included in this study. Patient 4 had a subclinical seizure during one of two EPI series (see Chaudhary et al., 2016 for a report on the analysis of this data). Therefore, 2/12 sessions were not included in this study.

icEEG pre-processing and analysis
Two IED classification strategies were employed for subsequent BOLD mapping using general linear modelling (GLM): firstly, conventional visual IED classification (Vulliemoz et al, 2011); and secondly, automated IED classification using a version of WC (Quian Quiroga et al., 2004) adapted specifically for epileptic discharges on icEEG (Sharma et al., 2017).
Pathological EEG patterns commonly found in invasive recordings include individual IEDs (spikes and sharp waves), repetitive IEDs (polyspikes) and paroxysmal fast activity (PFA) (Widdess-Walsh et al., 2007). Individual and repetitive IEDs can occur as single isolated epileptiform discharges (SED) or continuous epileptiform discharges (CED) (Palmini et al., 1995, Turkdogan et al., 2005. In two patients (# 3 and 8), there were CED with more than 1 IED/sec, which precluded the use of WC for their classification, due to the close proximity in time of the discharges, which causes difficulties in isolating them within the time window used in the WC process (see Sharma et al., 2017). In another patient (#4), polyspikes were observed that had the same field as individual IED; these were not classified using WC.
Prior to IED detection and classification, offline, average template subtraction-based, correction for MR scanning artefacts was applied to the EEG (Allen et al., 2000) and the resulting signals were down sampled to 250 Hz, and band-pass filtered (2-70Hz).

Visual IED marking and classification
EEG reviewer 'H1' (UJC) marked all clear or suspected individual SED and CED using the BrainVision Analyser software (Brain Products, Germany), based on their morphology and field using two types of event markers: point markers for individual IEDs, and onset and offset markers for repetitive IEDs (CED) and PFA.

Automated (Wave_clus) IED classification
Following our previous work, WC (Pedreira et al., 2014;Sharma et al., 2017) was used to M A N U S C R I P T

A C C E P T E D ACCEPTED MANUSCRIPT
automatically classify the IEDs identified by H1 (see above): SEDs and CEDs (< 1 IED/sec). In summary, Wave_clus is an automated neuronal spike classification algorithm that identifies and exploits small but consistent differences across multiple waveforms using wavelet decomposition and a superparamagnetic clustering algorithm (Quian Quiroga et al., 2004). For the purpose of automatically classifying icEEG IEDs, this process was carried out over the event waveforms captured on between 8 to 14 icEEG channels (the channels of interest) in which the IEDs were noted in the clinical EEG report as being most frequent and prominent (see Sharma et al., 2017 for a more detailed description). The result of this process was a set of IED classes (and corresponding labels) with every IED assigned to one of the classes. Polyspikes with the same field as a given Wave_clus class were subsequently added to this class.

Electro-clinical labelling of IED classes
In order to facilitate the interpretation and presentation of the findings, the IED classes were also labelled based on two criteria: 1) extent of their field distribution and 2) spatial relationship to the EZ.

Extent of field distribution
For field distribution, IEDs were labelled according to their spatiotemporal localization and distribution across the implanted electrodes (Lüders et al., 2000) as either: • Focal: if they involved up to 4 contiguous electrode contacts, • Regional: if they involved 5-10 contiguous electrode contacts, • Widespread: if they involved more than 10 contiguous electrode contacts, or: • Non-contiguous: if they had a focal or regional field but also extended to noncontiguous electrode contacts Relationship to the EZ IED classes were also labelled according to their spatial relationship to the EZ. All IED classes recorded in, and limited to, the brain area overlapping the EZ were labelled as IZ1 ('primary irritative zone'). All IED classes extending to brain areas outside the EZ were labelled as IZ2 ('secondary irritative zone') (Bettus et al., 2011). In this study we focused on the IZ1 IED, as they are the most clinically relevant (Cunningham et al., 2011, Vulliemoz et al., 2011, An et al., 2013.

fMRI data analysis
As mentioned previously, we created two GLMs to map the IED-related BOLD changes M A N U S C R I P T A C C E P T E D Each IED class was modelled as a separate effect of interest. For the CED with more than 1 IED/sec which were not classified using WC (one class for patients 3 and 8, as described above) the visual classification was used as an effect of interest in GLM1 and GLM2; furthermore, in patient 4 visually identified polyspikes (which could not be classified automatically) were added to a WC class that contained individual IED with the same field distribution, for inclusion in GLM2.
Individual IED markers were represented as a zero-duration stick function and polyspikes and paroxysmal fast activity were represented as variable duration blocks and convolved with the canonical hemodynamic response function, and its temporal and dispersion derivatives (canonical HRF+TD+DD).
The first two volumes of the fMRI time series data were discarded to account for the T1saturation effect; slice timing correction, scan realignment to the mean and spatially smoothed using an isotropic Gaussian kernel of 8mm Full Width Half Maximum (FWHM) were employed (Friston et al., 1995). For patients who underwent two EPI series (patients 5, 6, 7, 8), these were included in a single GLM as separate sessions.
Twenty-four inter-scan realignment parameters -6 realignment parameters from image pre-processing and a Volterra expansion of these (Friston et al., 1996) were included in the GLM as confounds to account for motion-related effects. We then applied a robust weighted least squares toolbox (Diedrichsen and Shadmehr, 2005) to reduce the influence of potential physiological and other sources of noise and artefacts.
For each IZ1 IED class, the presence of significant BOLD changes was assessed over the whole brain using a F contrast across the canonical HRF+TD+DD regressors, at a statistical threshold of p< 0.001 (uncorrected for multiple comparisons) and a cluster size threshold of 5 contiguous voxels (Grouiller et al., 2011, Centeno et al., 2017. The resulting SPMs were co-registered with pre-and post-surgical T1-weighted MRI scans using the rigid-body registration tool in SPM.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
For each IZ1 IED class, we evaluated the spatial concordance with the EZ, irrespective of the sign of BOLD change considering that both increases and decreases have been found in the epileptogenic zone (Thornton et al., 2011, Chaudhary et al., 2012Pittau et al., 2013). Similarly to our previous work (Coan et al., 2016) the concordance of each IEDrelated BOLD map with the EZ was assessed as either: • Concordant (C): when one or more BOLD clusters overlapped with the area of surgical resection or were within up to 2 cm (Cartesian distance, within the same lobe) of the resection margin. For each patient we calculated the proportion of C maps for GLM1 and GLM2: (number of C maps)/(number of maps); the Wilcoxon signed-rank test was used at the group level.

RESULTS
Six of the eight patients were diagnosed with frontal lobe epilepsy, one had temporal lobe epilepsy and one had parietal lobe epilepsy (see Table 1). The mean number of months they were seizure free post resective surgery was 64 months (SD: 12mo; range: 48 -87mo) (see Table 2). Across all subjects, the mean number of IED events detected by H1 was 1206 (SD: 989; range: 460-3567; Table 2) and the mean number of events classified by WC was 1105 (SD: 915; range: 277-3246; Table 3). In the following we summarise the IED classes obtained by the EEG reviewer and the automated algorithm.
In relation to IED field, for GLM1 there was a degree of concordance in 14/18 of Focal, 8/8 of Regional, 1/1 of Widespread, and 1/1 of Non-contiguous IED classes. For GLM2, there was a degree of concordance in 14/16 of Focal, 4/6 of Regional and 2/2 of Noncontiguous classes. See Supplementary Tables 2 and 3 for more complete descriptions of the GLM1 and GLM2 BOLD maps respectively.

Case Reports
To better illustrate the range of results obtained we present four case reports: two in which GLM2 showed greater BOLD concordance (patients 4 and 5), one in which concordance was the same for GLM1 and GLM2 (patient 6) and one in which the maps obtained from the automated classification were less concordant (patient 3). We report only on the results for the IZ1 IED.

Patient 4: Improved Concordance
This patient was diagnosed with temporal lobe epilepsy and the EZ was located in the left anterior temporal lobe.
Automated classification: The automated algorithm classified the IEDs into four classes, two of which were in IZ1 (classes 1 and 2) totalling 405 focal spikes (N=360 and 45, respectively; see Figure 2).

BOLD Results:
GLM1: the concordance for the three IZ1 classes was: C for classes 1 and 2, and D for classes 3 and 4. For classes 1 and 2 the global maximum was located in the left lateral occipital lobe; the GM was in the right superior frontal gyrus for class 3 and in the left orbital frontal cortex for class 4.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
GLM2: 67 visually classified polyspikes were found to match the field distribution and morphology of WC class 2 and therefore included in GLM2 in the same effect as WC class 2 (for a total of 112 IED). The concordance for the two IZ1 classes was: C. For class 1 the global maximum was located in the anterior cingulate cortex and for class 2 the global maximum was located within the EZ.

Patient 5: Improved Concordance
This patient was diagnosed with frontal lobe epilepsy and the EZ was located in the right anterior orbitofrontal region.
Visual classification: Reviewer H1 classified the IEDs into five classes, two of which were labelled IZ1 (classes 1 and 2) comprising a total of 223 IED; class 1 consisted of focal spikes N=142), while class 2 consisted of regional spikes; see Figure 3.
Automated classification: The automated algorithm classified the IEDs into three classes, two of which were labelled IZ1 (classes 1 and 2) comprising a total of 924 IED; class 1 consisted of focal spikes (N=333) and class 2 consisted of regional spikes (N=591).

BOLD Results
GLM1: the degree of concordance for the two IZ1 classes was: C for class 1 and D for class 2, and the global maximum was located in the right mesial occipital lobe for both.
GLM2: the degree of concordance for the two IZ1 classes was: C for classes 1 and 2. For IED class 1, the global maximum was located in the right occipital lobe, and for IED class 2 it was located in the left superior parietal lobe.

Patient 6: No Change in Concordance
This patient was diagnosed with frontal lobe epilepsy and the EZ was located in the right supplementary motor area and right superior frontal gyrus.
Visual classification: Reviewer H1 classified the IEDs into five classes, three of which were labelled IZ1 comprising a total of 733 IED (classes 1-3); classes 1 and 2 consisted of focal spikes (N = 211 and 46, respectively), while class 3 consisted of regional spikes (N=476); see Figure 4.
Automated classification: The automated algorithm classified the IEDs into three classes, two of which were labelled IZ1 (classes 1 and 2) and consisted of focal spikes (N=364 and 250, respectively for a total of 614).

BOLD results:
GLM1: the concordance for the three IZ1 classes was C. For class 1 the global maximum was located in the left orbital frontal cortex, for class 2 the global maximum was located in the right inferior frontal gyrus and for class 3 the global maximum was located in the right fronto-temporal lobe.
GLM2: the degree of concordance for the two IZ1 IED classes was C. For class #1 the global maximum was located in the left frontal pole, and in the right posterior temporal lobe for class 2.

Patient 3: Decreased Concordance
This patient was diagnosed with frontal lobe epilepsy and the EZ was in the left anterior inferior and middle frontal gyri.
Visual classification: Reviewer H1 classified the IEDs into 3 classes, two of which were labelled IZ1 (classes 1 and 2) comprising a total of 1035 IED; class 1 consisted of focal IEDs (N = 770) with (CED pattern of >1 per sec), and class #2 consisted of regional IEDs (N = 265); see Figure 5.
Automated classification: The CED pattern (visual class 1) was not classified using WC, but included in GLM2. WC classified the remaining 313 IED into three classes, one of which was labelled as IZ1 (class 1), consisting of focal 75 IEDs.

BOLD results:
GLM1: the concordances for the two IZ1 classes were: D for class 1 and C for class 2. For class 1 the global maximum located in the right lateral occipital lobe and in the EZ for class 2.
GLM2: the degree of concordance for the two IZ1 classes (Visual class 1 and WC class 1) was D. For Visual class 1 the global maximum was located in the posterior cingulate cortex and for WC class 1 the global maximum was located in the left posterior temporal lobe.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
Our aim was to determine whether classifying intracranially recorded IEDs with simultaneous fMRI using an automated spike classification algorithm could provide more biologically meaningful IED-related BOLD maps compared to those obtained following visual expert analysis in a group of patients with well characterised EZ.
To do this, IED-related BOLD maps obtained based on the automated algorithm were compared to the visual (classic / standard) classification (along the lines of our previous work on scalp EEG-fMRI;Pedreira et al., 2014) in terms of the presence or absence of BOLD changes in the EZ. To facilitate interpretation of the findings, we limited this analysis to the IEDs observed within the EZ (IZ1 IEDs) since we can be more confident of the location of their generators than the propagation IEDs.
One of the main challenges associated with the use of icEEG data for fMRI mapping is the difficulty of forming a parsimonious model of potential BOLD changes that reflects the complex spatio-temporal dynamics of IEDs. We found that the IED-related BOLD maps were concordant to the area of resection for 83% of the IZ1 classes for the automated algorithm compared to 75% for the visual classes (see table 3

Clinical and biological significance
To date, there have only been three studies focusing on the whole-brain mapping of IEDrelated BOLD changes using simultaneous icEEG-fMRI (Vulliemoz et al., 2011, Aghakhani et al., 2015. It is important to note that in our centre, the main motivation for the work based on icEEG-fMRI (Vulliemoz et al., 2011, Chaudhary et al, 2016, Murta et al, 2016, Ridley et al, 2017 is to gain a better understanding of the neurobiology of epilepsy, rather than to evaluate it as a potential clinical tool. Therefore, we propose that the present work's potential clinical significance lies mostly in its putative impact on the classification of IED in intracranial recordings; As explained previously and discussed further below, here the primary role M A N U S C R I P T

A C C E P T E D ACCEPTED MANUSCRIPT
of the simultaneously-acquired fMRI data is in providing a basis for the comparison of two types of IED classification.
In previous icEEG-fMRI studies of IED, these were detected and classified by EEG reviewers (Vulliemoz et al., 2011, Aghakhani et al., 2015. However, our previous study confirmed that the classification of icEEG IEDs can be inconsistent between EEG reviewers and that the automated classifications obtained using Wave_clus fall within inter-observer variability (Sharma et al., 2017).
Furthermore, we also previously demonstrated that the same algorithm, when applied to scalp EEG, showed an improved fMRI localisation compared to visual classification as indicate that the statistical analysis underlying WC classification may yet reveal a specific capability of classifying IEDs that is different from humans. Therefore, in the absence of a ground truth in terms of IED classification, and similar to Pedreira et al.
(2014), we sought to compare IED classification through independent data, namely the associated BOLD changes and in particular their relationship with the EZ.
Previous icEEG-fMRI studies involved in mapping IED-related BOLD across the whole brain assessed concordance based on location of a significant BOLD cluster to the most active spiking electrode contact (Vulliemoz et al., 2011, Aghakhani et al., 2015. The largest study carried out by Aghakhani et al. (2015) speculated that discordant (BOLD changes outside the lobe of the most active electrode) IED-related BOLD changes was due to the EZ being located elsewhere as this finding was found in patients with a poor postsurgical outcome (Aghakhani et al., 2015). This is consistent with previous scalp EEG-fMRI studies that indicate the absence of statistically significant BOLD changes in the area of resection can be predictive of a poor postsurgical outcome (Thornton et al., 2010, An et al., 2013, Coan et al., 2016, Centeno et al., 2017. Therefore, we assessed the presence of IZ1 IED-related BOLD changes using the best available gold standard defined as the area of resection in good postsurgical outcome patients (Grouiller et al., 2011).
In this study, any significant BOLD cluster was labelled concordant if it overlapped or was within 2cm of (within the same lobe as) the area of surgical resection (Coan et al., 2016). Some previous scalp EEG-fMRI studies on the clinical relevance of IED mapping  An et al., 2013). This is mainly based on two considerations: the cluster containing the statistical maximum is unique for any given map, which is a convenient simplification, and it represents the most likely location of the generator by some, albeit indirect, measure of the statistically strongest haemodynamic change. However, the sensitivity of the GM in localising the EZ is still inconclusive, as some studies have shown that it localises the EZ (defined as the area of resection in good postsurgical outcome patients) in less than half of the patients (Grouiller et al., 2011, An et al., 2013, Centeno et al., 2017. As mentioned previously, we defined concordance similar to Coan et al. (2016), as the presence of any significant BOLD cluster within or close to the EZ, and who showed that the presence of a significant BOLD cluster within or in close proximity to the EZ is strongly associated with good outcome (and absence of any cluster in the EZ, with a poor outcome). Indeed, in the present study only a small percentage of IED datasets showed the GM to be within or in the immediate proximity of the EZ for BOLD maps obtained using the visual and automated approach (16% and 17% respectively) (see Table 3) raising further questions on the notion of the global maximum as a marker of EZ. Furthermore, as discussed under Methodological considerations, there are physical reasons why the choice of the GM as part of the concordance assessment may be suboptimal for icEEG-fMRI data. The most important property of any BOLD-based assessment scheme is that provides an even playing field based on the spatial relationship between regions of BOLD changes and the epileptogenic generators.
In this study, four patients showed an improved concordance following automated IED classification (patients 4, 5, 7 and 8). As mentioned previously, our main focus was comparing concordance levels for the classes obtained by WC and visual classification and not comparing the spatial distribution of the EEG classes themselves. However, the advantage of using our automated approach can be seen for patient 5 for example (see Figure 5). In this patient, the EEG reviewer and WC identified two IZ1 IED classes involving similar electrodes with same distribution: one with a focal distribution involving the FP electrodes and the second with a regional distribution involving FP and AM electrodes (see Table 3 and Figure 3). Close examination of the two classifications shows that WC identified more IEDs for both classes (FP-focal: N = 333, FP-AM regional: N = 591) compared to the visual classification (FP-focal: N = 142, FP-regional: N = 81).

Methodological considerations and future work
The WC algorithm is based on the statistical characterisation of the entire IED population (in any given recording), and therefore it differs quite fundamentally from human expert classification, which is mostly done in sequence and focused on individual events, with a certain degree of consideration of the rest of the IED population. This is a way of integrating inherent patient-specific characteristics into the expert assessment, along with more general knowledge; in the statistical approach, the statistical weight of any feature (number of events, size of the difference with the rest of the sample) is factored in in a more global and objective/categorical way. This means that it is possible that certain statistically minor (sparse) events, which may be deemed clinically important by any given expert, may be overlooked by the WC feature extraction process.
Our automated classification pipeline was designed for the automated classification of individual IEDs as reflected in its use of a fixed-length time window centred around the peak of the IED to incorporate the maximum information of the IEDs with smallest possible time span to minimise the impact of noise (Sharma et al., 2017). Therefore it is not suitable in its current form to be applied to the EEG patterns with a fast repetitive nature such as polyspikes, PFA and very continuous CEDs; commonly observed EEG patterns in FCD patients (Palmini et al., 1995, Turkdogan et al., 2005, Widdess-Walsh et al., 2007. Therefore, to ensure that our BOLD models were as statistically efficient as possible, we included a manual step in our classification using Wave_clus. For example, regarding CEDs the EEG reviewer made a note of the frequency of the CEDs and if they occurred >1 per second, they were not automatically classified but modelled in the GLM; this was carried out in two patients (patients 3 and 8) (see Figure 5 for an example for patient 3). The automated classification of IEDs for patient 3 resulted in more discordant maps than the visual classification (see Table 3 and Figure 5). We note that uniquely for this patient the channels of interest, which were selected entirely based on information gained from the clinical report, included those showing frequent CEDs (channels DA3, DA4; see Figure 5 and Table 2). Therefore, the continuous nature of this M A N U S C R I P T

A C C E P T E D ACCEPTED MANUSCRIPT
EEG pattern, which the algorithm is not designed to handle, may have resulted in suboptimal classification of the individual IEDs.
The loss of MR signal in the vicinity of the electrodes in our data dataset due to the susceptibility artefact (Carmichael et al., 2012) is likely to result in reduced statistical power to detect local BOLD changes, which may partly explain why the GM was found to be distant to the EZ in some cases, providing further justification for the concordance scheme used here.
However, one of the main potential advantages of icEEG-fMRI compared to scalp EEG-fMRI is the opportunity to map the haemodynamic response associated with very focal discharges, as observed on depth EEG for example. In this study we defined focal discharges as IEDs occurring in up to 4 contiguous electrode contacts. The IED-related BOLD maps acquired using WC show that 87% of the WC IZ1 focal IEDs (13/15) show significant BOLD changes in the EZ as well as distant from the EZ (see Supplementary   Table 3 and Figure 3 for an example). These findings show that even for IEDs that are classified as being focal, distant regions maybe involved in their generation. Therefore, our results further reinforce the notion that even for very focal discharges, the epileptic network can be widespread, involving regions responsible for generating the IED and seizures, and regions remote from this area.
The practical benefits of our approach include time saving and dispensation of the need for reproducibility studies (Flanagan et al., 2009) and therefore, the uncertainty associated with changes in human EEG rater marking. Similar to previous studies using WC (Pedreira et al., 2014, Sharma et al., 2017 this study focused on the classification of IEDs detected by an EEG reviewer. In this regard it would be useful to also automate the detection of IEDs and then classify the markings using WC. In this work, we chose to map the BOLD changes associated with each IED class separately, even though it may have been possible to combine some classes as a single effect of interest based on their field distributions and morphologies: essentially, it is conceivable that some classes represent variations in the activity of a single generator (for example, all the IED recorded in the left temporal lobe of patient 4). However, this is a complex issue that belongs more specifically to the field of fMRI modelling and which we have previously studied for IED recorded on the scalp during fMRI (Salek-Haddadi et al, 2006).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
The application of our technique to a group of patients with varied surgical outcome would allow us to better study the biological of the observed BOLD patterns. Such studies could also include a consideration of their relationship with the seizure onset zone, in addition to that with the EZ, which is generally more extensive, and perhaps shaed additional light on their associated brain networks.

Conclusion
We found automated IED classification on icEEG data recorded during fMRI to be feasible and to result in IED-related BOLD maps that may contain similar or greater biological meaning compared to the conventional approach in the majority of the cases studied. Furthermore, we have shown that the BOLD changes associated with even very

M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT