Automated fusion of multimodal imaging data for identifying epileptogenic lesions in patients with inconclusive magnetic resonance imaging.

Many methods applied to data acquired by various imaging modalities have been evaluated for their benefit in localizing lesions in magnetic resonance (MR) negative epilepsy patients. No approach has proven to be a stand-alone method with sufficiently high sensitivity and specificity. The presented study addresses the potential benefit of the automated fusion of results of individual methods in presurgical evaluation. We collected electrophysiological, MR, and nuclear imaging data from 137 patients with pharmacoresistant MR-negative/inconclusive focal epilepsy. A subgroup of 32 patients underwent surgical treatment with known postsurgical outcomes and histopathology. We employed a Gaussian mixture model to reveal several classes of gray matter tissue. Classes specific to epileptogenic tissue were identified and validated using the surgery subgroup divided into two disjoint sets. We evaluated the classification accuracy of the proposed method at a voxel-wise level and assessed the effect of individual methods. The training of the classifier resulted in six classes of gray matter tissue. We found a subset of two classes specific to tissue located in resected areas. The average classification accuracy (i.e., the probability of correct classification) was significantly higher than the level of chance in the training group (0.73) and even better in the validation surgery subgroup (0.82). Nuclear imaging, diffusion-weighted imaging, and source localization of interictal epileptic discharges were the strongest methods for classification accuracy. We showed that the automatic fusion of results can identify brain areas that show epileptogenic gray matter tissue features. The method might enhance the presurgical evaluations of MR-negative epilepsy patients.

INTRODUCTION

Epileptic seizures can be suppressed with medication in most patients. In 30% of cases, the seizures are resistant to antiepileptic drugs and surgery is considered as an optimal treatment option (Guerrini, Sicca, & Parmeggiani, 2003) in indicated patients, in particular in those with lesions that are clearly visible on magnetic resonance (MR) images. During the presurgical evaluation process, some candidates are revealed to be unsuitable for epileptic surgery, often due to the absence of any conclusive detectable lesion in standard presurgical MR. In these nonlesional focal epilepsy (NLE) patients, surgical resective treatment is possible, but invasive exploration is required and the postsurgical outcome is inferior to the results in patients with evident structural lesions (Smith et al., 1997).

Epilepsy is a network disease (Tellez‐Zenteno, Ronquillo, Moien‐Afshari, & Wiebe, 2010; Vytvarova, Marecek, Fousek, Strycek, & Rektor, 2017); nevertheless, a seizure onset zone (SOZ) based on a lesion is the primary cause of epileptic seizures in focal epilepsy and also in NLE cases (Dachet et al., 2015).

The precise identification of a lesion has long been a research goal in neuroscience. A series of imaging methods were evaluated for their potential benefits in localizing lesions or SOZ (Kini, Gee, & Litt, 2016). To the knowledge of the authors, no method has yet been published that has been widely accepted as a stand‐alone method for lesion detection in NLE. The usual conclusion is that some individual methods can significantly contribute to presurgical evaluation and that a fusion of results from a set of methods and modalities is desirable.

During the last decade, the fusion of results from a set of methods has moved from clinical consultations to automatic processing using statistical methods developed for solving general statistical problems of classification. The approach works well in patients who have lesions that are visible in standard clinical MR. Several studies that fused various metrics describing cortex geometry and local characteristics of gray matter showed up to 80% accuracy in classifying a subject as a healthy control (HC) or as a patient with or without mesial temporal sclerosis (Lai, Guo, Cheng, & Wang, 2017; Rudie, Colby, & Salamon, 2015). Even higher accuracy—up to 90%—was reported in studies that fuse various metrics of cortex morphology for automatic delineation of visible lesions (Ahmed et al., 2014; El Azami, Hammers, Costes, & Lartizien, 2013; Fellah et al., 2012).

The identification of SOZ has been less explored in NLE patients. A study with 13 NLE showed concordance with clinical hypothesis in 46% of cases by fusing T1 and FLAIR images (Kotikalapudi et al., 2018). Another study examined 15 pharmacoresistant patients by fusing quantitative T1, T2, and tissue fraction maps. The authors identified lesions in four NLE patients that were highly concordant with the patients' electroclinical presentations (Ma et al., 2019).

Remarkably, very few studies with NLE offer validation either by surgery treatment outcome, ensuing histopathology, or invasive electroencephalography (EEG). Moreover, almost all fusing studies employ rather small sets of methods, usually based on a single recording modality, and thus capture the tissue features with limited view angles.

We collected clinical, electrophysiological, MR, and nuclear imaging data from 150 patients with pharmacoresistant epilepsy who were diagnosed as NLE or had inconclusive MR findings. We fused a series of neuroimaging methods within the frame of a classification problem to reveal several classes of gray matter tissue. We used two disjoint subsets of a group of 32 patients who underwent surgical treatment with known postsurgical outcomes and histopathology to identify and validate classes that are specific for epileptogenic tissue. The main aim was to evaluate the classification accuracy of the proposed method at a voxel‐wise level, that is, at the scale of millimeters. Further, we evaluated the contribution of individual methods to the identification of epileptogenic tissue and recommended a battery of methods for identifying the SOZ in patients with focal MR‐negative epilepsy.

METHODS

Subjects

We recruited patients with refractory focal epilepsy who had been referred to the Brno Epilepsy Center as epilepsy surgery candidates. The whole cohort of patients comprised 150 subjects. A majority (127 patients) were ranked as NLE cases due to negative findings provided by an experienced neuroradiologist based on presurgical MR epilepsy imaging. In a minority subgroup (18 patients), the reading of MR imaging was either inconclusive, that is, the radiologist could neither fully confirm nor negate a lesion, or there was a suspected lesion (5 patients) that could not be identified as epileptogenic due to the divergence from other diagnostic methods, in particular ictal video‐EEG. Further, we recruited a group of 110 age‐ and gender‐matched HC subjects. For details on the HC group, see Supplementary Methods. All subjects signed an informed consent form before entering the study. The study was approved by a local ethic committee.

Out of the 150 patients, a subgroup of 32 underwent surgery treatment (surgery group). We evaluated the outcomes after at least 6 months of follow‐up care and collected the results of histopathological examination (Table S1). The other subgroup of 118 patients did not undergo surgery because they had inconclusive presurgical evaluation results, they were still on a waiting list, or the testing of new drugs was not yet closed (nonsurgery group).

We divided all patients to form reasonable groups for the classification task. We compiled a final training group (TRG) from the nonsurgery group, and patients from the surgery group who had either an outcome other than ENGEL I with positive histology findings or any ENGEL outcome with negative histology findings. The large cohort of TRG subjects enabled a robust estimation of the set of gray matter tissue classes. A subgroup from TRG—those patients with an outcome other than ENGEL I with positive histology findings (17), testing group A (TEGA)—was then used to identify classes that are specific to tissue with proven pathology. The other patients from the surgery group—those with an outcome of ENGEL I with positive histology findings (11), testing group B (TEGB)—was the crucial subgroup of patients in whom surgery successfully targeted the SOZ; this group was used to validate the identified epileptogenic gray matter tissue classes.

We excluded 13 patients due to insufficient data (technical problems during data acquisition or poor quality). The details on each group are presented in Table 1. The grouping is depicted in Figure 1.

TABLE 1

Details on study groups

	# Subjects	Age (years) median (Q2–Q3)	Gender M/F
All patients	137	33 (25–39)	78/59
TRG	126	33 (25–39)	74/52
TEGA	17	30 (21–39)	13/4
TEGB	11	33 (25–44)	4/7

Abbreviations: F, female; M, male; TEGA, testing Group A; TEGB, testing Group B; TRG, training group; Q, quartile.

FIGURE 1

The details on grouping of patients. NRES, nonresected; RES, resected; H+, positive histology; H−, negative histology; TRG, training group; TEGA, testing Group A; TEGB, testing Group B; the numbers denote number of patients in each group

Data

Methods of interest

The battery of methods used in the presented study was set according to several rules. All data were acquired using neuroimaging modalities widely available in epilepsy centers, that is, MR, high‐density scalp EEG (HDEEG), positron emission tomography (PET), and single‐proton emission computed tomography (SPECT). The methods capture tissue features that might be formed during the long‐term effects of epilepsy or by tissue activity specific for epilepsy in ictal or interictal states. We excluded methods that in principle reflect the instantaneous functional states of the brain during interictal periods with ambiguous relations to epilepsy, such as resting‐state functional MRI (fMRI) or any connectivity measures based on fMRI or HDEEG. Considering that the proposed method aims on nonlesional cases, there is no prior information on localization. We would have to deal with whole brain activity that might be an intractable problem. Finally, each method's outcome had to be in the form of a voxel‐wise parametric spatial map to facilitate method fusion. Our battery of methods comprised MR‐based methods, HDEEG‐based methods, and nuclear imaging methods.

MR‐based methods

We employed a voxel‐based morphometry approach (Ashburner & Friston, 2005; Kotikalapudi et al., 2018; Muhlau et al., 2009). It is based on the automatic segmentation of high‐resolution T1 into tissue compartments and can reveal gray matter changes in individual subjects (Kini et al., 2016). We included methods sensitive to local gray matter intensity changes (GMC) and local gray matter volume (GMV) changes. The segmentation outcome can be further used to assess gray‐white matter boundary blurring. The method called junction (JUN) can identify subtle blurring in focal cortical dysplasia (FCD) lesions (Bernasconi & Bernasconi, 2011; Colliot, Antel, Naessens, Bernasconi, & Bernasconi, 2006).

We used arterial spin labeling (ASL), which quantifies cerebral blood flow (Boscolo Galazzo et al., 2015) and has a potential to reveal perfusion abnormalities in patients with focal seizures (Lee, Kwon, & Lee, 2019).

Further, we acquired diffusion‐weighted imaging (DWI) data. DWI is sensitive to the extent and direction of water molecule diffusion and provides unique information on tissue microstructure measured in‐vivo (Colombo, Salamon, Raybaud, Ozkara, & Barkovich, 2009). We used fraction anisotropy (FA) and mean diffusivity (MD) methods that were shown to reveal diffusion alterations in gray matter in patients with temporal lobe epilepsy (Winston et al., 2020). Further, we used mean kurtosis (MK). The MK has been reported to identify gray matter microstructure alterations in temporal lobe epilepsy (Bonilha et al., 2015).

HDEEG‐based methods

HDEEG in combination with electrical source imaging algorithms (Michel & Brunet, 2019) mitigate the statement about low spatial resolution of scalp EEG. In the last decade, HDEEG has often been used to localize interictal epileptiform discharges (IED) generators (Megevand & Seeck, 2020).

Nuclear imaging methods

PET imaging is an established method widely used to localize SOZ. It shows high sensitivity to the hypometabolism specific to epileptogenic tissue (Juhasz & John, 2020). Interictal and ictal SPECT imaging can provide valuable information on SOZ localization (Krsek et al., 2013). We utilized a variant called STATISCOM that subtracts ictal and interictal SPECT and compares the resulting image with data from an HC cohort (Sulc et al., 2014).

All but one of these methods resulted in voxel‐wise T‐value statistical maps that show the normalized deviation of a single patient's data from a norm estimated based on cohorts of HC subjects. We based all ensuing processing on T‐values as they inherently take into account the size of control cohort and suppress the potential bias that could be introduced by control cohorts' sizes variability over methods. The IED is the only method that is based solely on the patient's interictal HDEEG data. Table 2 summarizes all 10 methods used in our study. For details on the preprocessing and statistical analyses for each method, see Supplementary Methods.

TABLE 2

The list of methods

Abbreviation	Modality	Method	Interpretation from the patient data perspective, relative to HC (excluding IED)
GMC	MR	Local gray matter concentration	Increase: Alterations in intensity of gray matter in structural images
GMV	MR	Local gray matter volume	Decrease: Gray matter atrophy
JUN	MR	Gray–white matter boundary blurring	Increase: Augmented intensity gradient between gray and white matter in structural images
ASL	MR	Quantified cerebral blood flow	Decrease: Decreased perfusion
IED	HDEEG	Source localization of interictal epileptic discharges
PET	PET	Metabolism alteration	Decrease: Hypometabolism
SPE	SPECT	STATISCOM—Evaluation of ictal vs. interictal perfusion	Increase: Source of early ictal activity
FA	MR	FA	Decrease: Decrease in directionality of molecule diffusion, unspecific microstructure alterations
MD	MR	MD	Increase: Enhancement of molecule diffusion, unspecific microstructure alterations
MK	MR	MK	Decrease: Enhancement of molecule diffusion due to loosing barriers, unspecific microstructure alterations

Abbreviations: ASL, arterial spin labeling; FA, fraction anisotropy; GMC, gray matter concentration; GMV, gray matter volume; HDEEG, high density electroencephalography; IED, interictal epileptic discharges; JUN, junction; MD, mean diffusivity; MK, mean kurtosis; MR, magnetic resonance; PET, positron emission tomography; SPE, single proton emission. SPECT, single‐proton emission computed tomography; STATISCOM, statistical ictal SPECT coregistered to MRI.

Data acquisition

The nuclear imaging data were acquired during the clinical evaluation of the patients (for details, please see Supplementary Methods).

All other imaging and electrophysiology data were acquired at CEITEC Masaryk University Neuroscience Center, which is equipped with a 3T Prisma Siemens MR machine and MR compatible 256 channel EGI EEG system. For details, please see Supplementary Methods.

Data preparation

The resulting statistical maps from all 10 methods differed in spatial resolution and the amount of smoothness due to the unequally set Gaussian spatial filters within the stage of individual method processing. To unify all data, we resampled each image into a space with an isotropic voxel‐size of 1.5 mm and smoothed it by a spatial filter with a Gaussian kernel (FWHM = 15 mm). The reason for smoothing was three‐fold: the smoothing helped to minimize the mutual intermethods registration inaccuracies (Mikl et al., 2008; Salmond et al., 2002), the Gaussian filter renders all data more normally distributed, and the equalization of smoothness ensures fairly the same impact on results smoothness. For details on smoothing, please see Supplementary Methods.

Finally, the data were restricted to voxels belonging to the cortex, hippocampi, and amygdale according to the gray matter template included in the SPM12 toolbox (www.fil.ion.ucl.ac.uk/spm/software/spm12/, probability threshold >0.1). We excluded the basal ganglia, thalamus, midbrain, and cerebellum as unlikely locations of SOZ.

The IED data were intensity normalized by the individual subject's maximal value.

Surgery mask

We constructed a surgery mask (SM) image for each subject from the testing groups. The SM had the character of a binary image that encloses the resected region. The SM was created by the manual delineation of postsurgery imaging.

Gray matter tissue classes estimation

Training process

We used a Gaussian mixture model to extract latent classes of gray matter tissue from the complete TRG. We employed an algorithm that can handle missing data (Delalleau, Courville, & Bengio, 2018) implemented in‐house in MATLAB R2017a (MathWorks, Inc.). We estimated the model for single to 10 classes and 10 times with randomly initialized parameters to preclude convergence into a local minimum. Each estimation was ended either when an average relative change in parameters dropped below 0.1% and a relative change in free energy dropped below 0.001% or when the algorithm failed to converge within 1,000 iterations. We watched Bayesian information criterion (BIC). An optimal number of classes (NC) was defined as the number where the BIC had the knee point, that is, where an additional class (and thus additional model parameters) does not explain notable part of previously unexplained variability; the solution with the optimal NC and the lowest BIC was used in ensuing analyses.

Voxel classification

We classified each voxel from each TEGA and TEGB subject using parameters of estimated classes (i.e., means and covariance matrices over methods). Each voxel was assigned probabilities of belonging to individual classes. As a result, we obtained a set of NC spatial probability maps for each subject.

Epileptogenic gray matter tissue classes

Selection of epileptogenic gray matter tissue classes

The selection and all ensuing testing and comparisons were based on the receiver operating characteristic (ROC) and area under ROC curve (AUC) approach: voxels inside SM were considered as actual positives; voxels outside SM were considered as actual negatives.

We computed ROC and AUC for each subject in TEGA based on its SM and class probability maps. This was repeated for each possible subset from the NC classes. When the subset comprised more than a single class, the voxel‐wise maximal probability value over classes was considered. We looked for an optimal subset of classes (OSC) with the highest average AUC over TEGA. We used a one‐sample t test to assess whether the AUC was significantly higher than 0.5, that is, better than random classification. The level of statistical significance was set to p < .05, Bonferroni corrected for the number of possible subsets (2^NC‐1).

Validation of epileptogenic gray matter tissue classes

The previously identified epileptogenic gray matter tissue classes, that is, the members of OSC, were validated using the TEGB group. We computed the AUC and statistics for the OSC in the same way as for TEGA group.

Methods importance evaluation

Using the OSC and all subjects from the testing groups (TEGA + TEGB), we computed AUC for all possible combinations of methods. We intended to examine the effect of missing data from a method and thus the importance of each method for the classification. We assessed the method's importance from two perspectives: What is the AUC change when data from a method is completely missing? and What is the systematic contribution of a method to the AUC level? The former question was addressed by a one‐sample t test applied to the differences between AUC computed with or without a method. To answer the latter question, we employed a linear mixed‐effects model that included interception, fixed effect factor for each method (with two levels for existing and missing data), and a subject factor as predictors and AUC values as response variable. The model was applied to the data with a subtracted value of 0.5. This transformation gave the intercept an interpretation of whether the overall classification was better than chance. The choice of the model was determined by the fact of missing data that precluded usage of simple N‐way repeated measures analysis of variance. The level of statistical significance was set to p < .05, corrected for the number of methods.

RESULTS

The TRG and TEGB did not differ in age (Mann–Whitney Z = 0.483, p = .63) or in gender balance (Fisher's exact, Chi2(1) = 2.06, p = .15).

All patients had available data from at least five methods (TRG 7, 7–8, median, quartiles; TEGA 8, 7–9; TEGB 8, 7–9). The amount of available data in TRG did not differ from TEGB (Mann–Whitney Z = 1.32, p = .21). No method was unequally frequent in TEGA versus TEGB (GMC, GMV, JUN: Fisher's exact Chi2(1) = 0, p = 1; ASL: Chi2(1) = 1.01, p = .31; IED: Chi2(1) = 1.62, p = .20; PET: Chi2(1) = 0.67, p = .41; SPE: Chi2(1) = 0.05, p = .82; FA, MD, MK: Chi2(1) = 0.001, p = .97), Figure S1.

The training algorithm converged before the maximal allowed iterations in all solutions. According the BIC, the optimal NC was six (Figure S2). Figure S3 shows resulting classes C1 to C6 from the best solution and their frequency in TRG data. The C2 class had a notably lower frequency of occurrence than the others.

The OSC with respect to voxel classification in TEGA was the combination of C2 and C4 with the AUC = 0.73 ± 0.18. The AUC was significantly higher than 0.5 (one‐sample t test: T(16) = 5.40, one‐sided p < .001). The C2 class comprises voxels with high deviation from HC in GMV (gray matter loss in patients), ASL (decreased perfusion), PET (strong hypometabolism), SPE (ictal perfusion higher than interictal), MK (decrease of barriers for water diffusion), and MD (enhancement of water diffusion). The C4 class comprises voxels with a high deviation from HC in SPE, MK, and MD. Moreover, C4 voxels had high values in IED, that is, likely source of IEDs (Figure 2).

FIGURE 2

Epileptogenic gray matter tissue classes C2 (red dots) and C4 (violet dots). The upper plot depicts each method's deviation from the healthy control norm in terms of T values. The lower plot depicts normalized current density estimated by source localization of interictal epilepsy discharges. The dots and lines depict mean and SD of the distribution over training voxels belonging to each respective class. For abbreviations, please see Table 2

The OSC showed even better classification in the TEGB group results than in TEGA: AUC = 0.82 ± 0.14. It was significantly higher than 0.5 (one‐sample t test: T(10) = 7.66, one‐sided p < .001). An outcome of voxel classification is depicted in Figure 3 (whole brain situation in Figure S4, results overlaid on standard clinical MR images in Figure S5).

FIGURE 3

Results of voxel classification in patients with ENGEL I surgery outcome and proven pathology (11 subjects from the testing group). The images show voxels belonging to epileptogenic gray matter tissue classes (C2 and C4; voxels where the highest probability had C2 are in orange and C4 are in violet; threshold set to show 1% voxels with highest values). The surgery extent is delineated by a yellow line. The top row shows cases with proven hippocampal sclerosis, the bottom row at left are cases with focal cortical dysplasia, and the rest are cases with gliosis

We estimated AUC for TEGA and TEGB subjects for all possible combinations of methods to simulate the missing data from each method. Considering just relevant cases, the total amount of AUC values through all subjects was 8,804 out of 28,644 (210 = 1,023 combinations of methods × TEGA(17) + TEGB(11) = 28 subjects). PET was the only method that caused a significant decrease in AUC when missing (Table 3a).

TABLE 3

(a) Decrease of classification accuracy in simulated situations with missing methods. The simulated missing PET significantly decreased classification accuracy (one‐sample T test; p < .05, Bonf. corrected, df =27). (b) Effect of each method on classification accuracy. (c) Effect of each method on classification accuracy with discarded data from nuclear imaging methods. For (b) and (c), the denoted methods significantly increase the accuracy (linear mixed effect model, p < .05, Bonf. corrected, df = 8,793). On average, the classification accuracy is better than chance as represented by significant effect of intercept

Method		Intercept	GMC	GMV	JUN	ASL	IED	PET	SPE	FA	MD	MK
(a)	EE	—	0.006	−0.008	0.003	0.034	−0.053	−0.056	−0.015	0.002	−0.013	0.002
	p‐Value	—	.34	.52	.58	.02	.34	.004	.20	.82	.25	.77
(b)	EE	0.140	−0.001	0.007	<0.001	−0.062	0.048	0.076	0.043	−0.006	0.047	−0.003
	p‐Value	<.001	.92	.55	.99	.01	.22	<.001	<.001	.68	.001	.69
(c)	EE	0.110	0.003	0.011	0.006	−0.022	0.154	—	—	0.012	0.069	0.001
	p‐Value	<.001	.76	.45	.57	.11	<.001	—	—	.51	<.001	—

Note: All significant p‐values are in bold.

Abbreviations: ASL, arterial spin labeling; EE, effect's estimate; FA, fraction anisotropy; GMC, gray matter concentration; GMV, gray matter volume; IED, interictal epileptic discharges; JUN, junction; MD, mean diffusivity; MK, mean kurtosis; PET, positron emission tomography; SPE, single proton emission.

The results of a linear mixed effect model showed significant effects of intercept, PET, SPE, and MD on a level of AUC. For all three methods, the effect was positive, that is, the presence of the method's data significantly increased the accuracy of voxel‐wise classification. The significant positive intercept proves that overall, the classification is higher than 0.5, that is, better than chance (Table 3b).

We also computed the linear mixed effect model in a variant without PET and SPE data. The results showed a significant and positive effect of intercept, IED, and MD (Table 3c). Please note that the overall accuracy of classification decreased only slightly, from 0.64 to 0.61 (intercept +0.5).

DISCUSSION

The search for appropriate surgical targets in patients with pharmacoresistant epilepsy has a long history. It started with the separate evaluations of various methods. Developments in this area have since shifted to the fusion of results derived from a set of methods (Kini et al., 2016). Interestingly, the studies thus far have comprised limited selections of methods restricted mostly to methods from single modalities. The studies have been rather vague in terms of the localization of the focus. Mostly they aim at differential diagnoses (Lai et al., 2017; Rudie et al., 2015) or the localization is assessed at the lobar level (Lascano et al., 2016).

We present results on automated fused data from several different modalities and derived methods. The proposed approach fuses the individual methods directly at the voxel level and assigns each voxel a membership in either a healthy or an epileptogenic gray matter tissue class. We analyzed data from 137 patients with pharmacoresistant epilepsy who were ranked as nonlesional, that is, MR negative or with inconclusive MR finding. In five patients, there was a suspected lesion that could not be identified as being epileptogenic due to the divergence with other diagnostic methods, in particular ictal video‐EEG. The analysis comprised postprocessing of structural MR images, MR perfusion images, MR DWI, source localization of IED, and nuclear imaging. We employed a Gaussian mixture model to reveal six latent classes of gray matter tissue in a large training subset of 126 subjects. The method was completely blind to the input data, that is, no clinical information regarding the considered localization of focus was supplied. Using a subgroup of patients with surgery outcome and known histopathology, we identified two classes that are specific for epileptogenic tissue. The voxel classification into epileptogenic or other class was 82% accurate in the patients with favorable outcome.

Epileptogenic gray matter classes

We identified two classes. The C2 class comprises voxels where patients had gray matter loss, decreased perfusion, hypometabolism, higher perfusion in ictal state compared to interictal state, decrease of barriers for water diffusion, and enhanced water diffusion. This is in accordance with previous findings: the GMV loss was shown especially in temporal‐lobe patients with HS (Bernhardt, Hong, Bernasconi, & Bernasconi, 2015); decreased perfusion was reported in children (Lee et al., 2019); the hypometabolism seen in PET and interictal/ictal‐specific perfusion assessed by SPECT are well‐known markers of epileptogenic tissue (Krsek et al., 2013; Rathore, Dickson, Teotonio, Ell, & Duncan, 2014); and microstructure abnormalities in epileptogenic tissue identified by diffusivity metrics were reported to be sensitive to FCD (Feindel, 2013; Winston et al., 2020). The C4 class is similar to that of C2. The difference is that it demonstrates voxels with high values in IED, which can be interpreted as an irritative zone, that is, a source of IEDs. Those two classes merged had the highest accuracy of classification. In a majority of the successfully operated (TEGB) patients, the C4 localized the surgery region. In two patients, it was C2. In all but one patient, the method identified a region inside the resected area. In subject S22, the identified region was located within proximity to the resected area.

Implications for clinical evaluation of patients

Our study shows that the automated fusing of methods can result in highly accurate SOZ identification, up to 82% in our data. Importantly, this result was achieved with no prior information on SOZ localization and even at the voxel level.

Further, we showed that in our data that nuclear and DWI methods, specifically, PET, STATISCOM, and MD, make the highest contribution to classification accuracy. This is in line with published high positive predictive value for accurate lobe identification using PET and SPECT data (Lascano et al., 2016). Other publications found increased intra cortical MD in patients with epilepsy (Winston et al., 2020). Our results suggest that even if nuclear imaging is achievable, it is worth supplying DWI. In our dataset, the MD marker significantly improves the classification accuracy beyond the level of PET and SPECT alone. On the other hand, nuclear imaging is invasive method that exposes the patient to radiation. Ictal SPECT imaging requires nonstop patient monitoring and precise timing of tracer administration and ensuing data acquisition. The analysis of the importance of methods using the discarded data from nuclear imaging methods showed only a slight decrease in classification accuracy; MD and source localization of IED emerged as a method that significantly increases accuracy. This is in accordance with Lascano et al., who showed that the source localization of IED with MR imaging was the best combination of methods (Lascano et al., 2016). Overall, our results suggest favoring nuclear imaging, high‐quality DWI, and source localization of IED when planning advanced imaging in the stage of searching for epileptogenic focus in nonlesional, MR‐negative epilepsy. We recommend to use also the other methods, including morphometry metrics, ASL, FA, and MK because they may become important for classification when the data from nuclear imaging, DWI and HDEEG are not available.

We consider our results as a proof of concept. Although the two identified classes are specific to our facility, our environment, and even our patient cohort, and are not directly transferable to other centers, we suppose that very similar epileptogenic classes would result from another epilepsy center's data. We recommend establishing site‐specific norms using HC groups and to train the GMM classifier with site‐specific epilepsy cohorts, ensuring all data are measured at the same machines.

Technical considerations

Missing data

The missing data situation was very frequent in our dataset. Overall, there were an average of eight accessible methods out of the 10 for TEGA and TEGB subjects. Therefore, we chose an algorithm that can handle missing data (Delalleau et al., 2018). As a consequence, the algorithm effectively widens the less frequent methods' SDs and thus decreases their impact in deciding which class a voxel belongs to (see Figure 2, SPE method). On the other hand, it allowed for the inclusion of almost all patients in the training and testing process.

Number of classes

The optimal number of six classes was chosen according to BIC criteria. We identified a number where the change in BIC reached its plateau. This was to some extent a subjective selection. We have therefore examined the results of classification with five and seven classes. We were able to identify classes very similar to C2 and C4 with similar classification results (Supplementary Results).

Smoothing

The presmoothing of the input data had several impacts on the results. While it could suppress some interesting localized effects, we believe the benefits it brought were more important (e.g., the reduction of intersubject and intermethod spatial registration inaccuracies). It is even recommended to apply rather strong spatial filtering when producing single subject versus group comparisons (Salmond et al., 2002), which is a predominant strategy in our data.

Methods selection

The selection of methods is a limitation of presented study. We restricted ourselves to methods that are commonly used in epilepsy, have an established methodology, and are easily implemented at the voxel‐wise representation level. We are aware that other methods have been published with promising impacts on lesion detection. Most importantly, we have not included surface‐based morphology metrics (Ahmed et al., 2014; Whelan et al., 2018) or T1/T2 intensity alterations metrics (Noth et al., 2020; Shultz, O'Brien, Stefanidou, & Kuzniecky, 2014). We also did not include simultaneous EEG–fMRI methods due to a lack of IED in almost all EEG data acquired during fMRI sessions (Pardoe & Kuzniecky, 2014; Yamazoe et al., 2019). All these methods might improve classification accuracy or could render some classes to be specific to any pathology.

CONCLUSION

We showed that the automatic fusion of results from a battery of individual methods results in highly accurate epileptogenic tissue localization. Importantly, the proposed method works at the voxel level, that is, at the scale of millimeters. The data and MATLAB code are available upon request from authors.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

ETHICS STATEMENT

We confirm that we have read the journal's position on issues involved in ethical publication and affirm that this report is consistent with those guidelines.

Table S1 Details on surgery group

Table S2: Methods' smoothness

Figure S1. Data availability. The numbers at diagonal show number of subjects with available data from respective modality. The off‐diagonal numbers show number of subjects who have data available from both respective modalities. TRG .. training group; TEGA .. testing group A; TEGB .. testing group B

Figure S2. The BIC of all Gaussian Mixed model realizations.

Figure S3. Tissue fingerprints of estimated classes (left) and their frequency of occurrence in training data (right).

Figure S4: Results of voxel classification in patients with ENGEL I surgery outcome and proved pathology. The images show voxels belonging to epilepsy‐specific gray matter tissue classes (C2 and C4; voxels where the highest probability showed either C2 or C4 are in orange or violet respectively; threshold set to show 1% voxels with highest values). The surgery extent is delineated by yellow line. The first to third rows (top to bottom) show cases with proved HS, the fourth and fifth are cases with proved FCD, the sixth are cases with proved gliosis.

Figure S5. Results of voxel classification in patients with ENGEL I surgery outcome and proved pathology overlaid on standard clinical MR images. The images show voxels belonging to epilepsy‐specific gray matter tissue classes (C2 and C4; voxels where the highest probability showed either C2 or C4 are in orange or violet respectively; threshold set to show 1% voxels with highest values). The surgery extent is delineated by yellow line. For subjects' details, please see Table S1.

Figure S6. The epilepsy‐specific gray matter tissue classes C52 (red dots) and C54 (dark red dots). The upper plot depicts each method's deviation from the healthy control norm in terms of T values. The lower plot depicts normalized current density estimated by source localization of interictal epilepsy discharges. The dots and lines depict mean and standard deviation of the distribution over training voxels belonging to each respective class. For abbreviations, please see Table 2.

Figure S7. The epilepsy‐specific gray matter tissue classes C72 (red dots) and C77 (dark red dots). The upper plot depicts each method's deviation from the healthy control norm in terms of T values. The lower plot depicts normalized current density estimated by source localization of interictal epilepsy discharges. The dots and lines depict mean and standard deviation of the distribution over training voxels belonging to each respective class. For abbreviations, please see Table 2.

Click here for additional data file.