Clinical utility of distributed source modelling of interictal scalp EEG in focal epilepsy

Clinical utility of distributed source modelling of interictal scalp EEG in focal epilepsy

Clinical Neurophysiology 121 (2010) 1726–1739 Contents lists available at ScienceDirect Clinical Neurophysiology journal homepage: www.elsevier.com/...

3MB Sizes 1 Downloads 89 Views

Clinical Neurophysiology 121 (2010) 1726–1739

Contents lists available at ScienceDirect

Clinical Neurophysiology journal homepage: www.elsevier.com/locate/clinph

Clinical utility of distributed source modelling of interictal scalp EEG in focal epilepsy C. Plummer a,b,*, M. Wagner c, M. Fuchs c, S. Vogrin a, L. Litewka a, S. Farish b, C. Bailey d, A.S. Harvey d,e, M.J. Cook a,b a

Centre for Clinical Neurosciences and Neurological Research, St. Vincent’s Hospital, 5th Floor Daly Wing, 35 Victoria Parade, Fitzroy, Victoria 3065, Australia Department of Medicine, University of Melbourne, Grattan Street, Parkville, Victoria 3052, Australia Compumedics Neuroscan, Heussweg 25, 20255 Hamburg, Germany d Department of Neurology, Royal Children’s Hospital, Flemington Road, Parkville, Victoria 3052, Australia e Department of Paediatrics, University of Melbourne, Grattan Street, Parkville, Victoria 3052, Australia b c

a r t i c l e

i n f o

Article history: Accepted 3 April 2010 Available online 8 May 2010 Keywords: EEG source localization Distributed modelling Dipole BFEC MTLE Focal epilepsy

a b s t r a c t Objective: Assess the clinical utility of non-invasive distributed EEG source modelling in focal epilepsy. Methods: Interictal epileptiform discharges were recorded from eight patients – benign focal epilepsy of childhood (BFEC), four; mesial temporal lobe epilepsy (MTLE), four. EEG source localization (ESL) applied 48 forward–inverse–subspace set-ups: forward – standardized, leadfield-interpolated boundary element methods (BEMs, BEMi), finite element method (FEMi); inverse – minimum norm (MNLS), L1 norm (L1), low resolution electromagnetic tomography (LORETA), standardized LORETA (sLORETA); subspace – whole volume (3D), cortex with rotating sources (CxR), cortex with fixed sources (CxN), cortex with fixed extended sources (patch). Current density reconstruction (CDR) maxima defined ‘best-fit’. Results: From 19,200 CDR parameter results and 2304 CDR maps, the dominant variables on best-fit were inverse model and subspace constraint. The most clinically meaningful and statistically robust results came with sLORETA–CxR/patch (lower Rolandic in BFEC, basal temporal lobe in MTLE). Computation time was inverse model dependent: sub-second (MNLS, sLORETA), seconds (L1), minutes (LORETA). Conclusions: From the largest number of distributed ESL approaches compared in a clinical setting, an optimum modelling set-up for BFEC and MTLE incorporated sLORETA (inverse), CxR or patch (subspace), and either BEM or FEMi (forward). Computation is efficient and CDR results are reproducible. Significance: Distributed source modelling demonstrates clinical utility for the routine work-up of unilateral BFEC of the typical Rolandic variety, and unilateral MTLE secondary to hippocampal sclerosis. Ó 2010 International Federation of Clinical Neurophysiology. Published by Elsevier Ireland Ltd. All rights reserved.

1. Introduction The increasing sophistication of the distributed modelling method for EEG source localization (ESL) has heightened debate on its clinical translatability to routine epilepsy practice. Even dipole modelling (an allied field boasting more clinical studies to date) remains largely tied to the research domain (Plummer et al., 2008). Indeed, the distributed model can be regarded as a three-dimensional latticework of multiple, point-like dipole models; each dipole carries its own orientation, location, and strength per time instant (also termed CDR, ‘current density reconstruction’) (Wagner et al., 2001).

* Corresponding author at: Centre for Clinical Neurosciences and Neurological Research, St. Vincent’s Hospital, 5th Floor Daly Wing, 35 Victoria Parade, Fitzroy, Victoria 3065, Australia. Tel.: +61 392883045; fax: +61 392883350. E-mail address: [email protected] (C. Plummer).

One theoretical advantage held by distributed (over dipole) modelling is that its algorithms address the inverse problem with fewer lead-in assumptions, or ‘priors’, mainly in respect of the number of sources that explain the measured EEG signal. Due to the ill-posed nature of the inverse problem however, (Helmholtz, 1853), solutions are highly under-determined (Fuchs et al., 1999). Post-process constraints are needed to arrive at a ‘best-fit’ solution. The mathematical tools employed to do this is essentially what distinguishes one distributed modelling algorithm from another. Constraints may also be anatomically informed (subspace constraint) whereby the three-dimensional solution space is specified (cortical versus head space). While there is a body of work in the biophysics literature on the application of distributed modelling to experimental simulations (Baillet and Garnero, 1997; Grova et al., 2006; Kim et al., 2002; Nummemnaa et al., 2007), comparatively little work has assessed the usefulness of this technology in a systematic way in the routine

1388-2457/$36.00 Ó 2010 International Federation of Clinical Neurophysiology. Published by Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.clinph.2010.04.002

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

clinical epilepsy setting. There are several immediate issues. First, much like functional Magnetic Resonance Imaging (f-MRI), CDR maps are typically threshold-set in order to curb the display of those elements of the map that are more likely to be representing noise. However, no standardized or agreed-upon threshold value has been determined to date. The second practical issue is that distributed modelling is generally more computationally demanding than is the case for dipole modelling. This is effectively why, at least at present, CDR solutions can only be calculated for the net distributed current per time instant, making it difficult for the investigator to tease apart activity that might be coming from multiple sources that overlap in space and time during the course of a spike or seizure discharge (Scherg et al., 1999). Third, while constraining the solution to a particular subspace within the volume conductor (or forward model) might minimize the final fit error in respect of the individualized cortical anatomy, it should be emphasized that however spatially constrained the final fit becomes, little will be gained if the choice of the forward or inverse model is inappropriate in the first place. And fourth, most clinical ESL studies select a particular source modelling set-up (a forward–inverse modelling combination with or without a subspace constraint) and test its performance across a particular patient cohort (a series of lesion positive epilepsy surgery cases for instance). ESL studies which compare the performance of multiple different types of distributed modelling approaches across the same patient cohort are genuinely lacking in the clinical literature. To our knowledge, the present study incorporates the largest number of distributed modelling set-ups to have been cross-examined in the clinical setting. A chief criticism of ESL generally, and of distributed modelling in particular, is that the mathematical constraints that may be applied to solve the inverse problem are not known to genuinely mirror actual electrophysiological behaviour. The algorithms that underpin the various distributed modelling approaches are technically loaded and, while these methods have produced encouraging results in the relatively high signal-to-noise environment of the simulated experiment when the true source configuration is known, their translatability to the routine clinical setting remains largely undefined. This is the primary motivation for the present study. More directly, what difference is made to the final distributed fit if the BEMi is used instead of the FEMi forward model; if LORETA is used instead of the sLORETA inverse model; or if CxR is used instead of the patch subspace constraint? How reproducible are the results? Are the results clinically meaningful when the number of electrodes applied is only in the order of 19–21, as is often the case for routine interictal EEG recordings? Is the modelling too cumbersome and time costly to carry out for prospective use in clinical practice? We have selected benign focal epilepsy of childhood (BFEC), the prototypic idiopathic focal epilepsy, and mesial temporal lobe epilepsy (MTLE), the prototypic symptomatic focal epilepsy, in the present study for two reasons. They represent two of the most common types of epilepsy encountered in routine clinical practice and they are arguably the most well characterized epilepsies in the field of non-invasive ESL at this point in time (Wong, 1998; Ebersole, 2000; Huiskamp et al., 2004; Pataraia et al., 2005; Ebersole and Hawes-Ebersole, 2007; Plummer et al., 2007).

1727

temporal lobe epilepsy (MTLE), 12–16 years (mean 14.3 years) secondary to hippocampal sclerosis on Magnetic Resonance Imaging (MRI) and histopathology (Engel class 1, 5–7 year follow-up). Stereotypic interictal epileptiform discharges (IEDs) of a single morphology and topography were present on visual assessment of the EEG in each patient (30-min 19-electrode recordings BFEC; prolonged 21-electrode pre-operative scalp recordings MTLE). International 10–20 electrode positions (Jasper, 1958) were used for both groups with accurate scalp measurement. In MTLE, two extra electrodes were placed at points one-third the distance from the outer canthus to the external auditory meatus (T1, T2) (Binnie et al., 1982). Standard 10 mm gold plated disc electrodes (pure silver cast cup, two millimetre central port) were attached to the scalp with SLE Collodion AdhesiveÒ and injected with high conductivity electrode gel (Parker Signa GelÒ). Recordings were at 256 Hz digital acquisition (16 bit ADC resolution) were performed with Compumedics ProFusionÒ software; 0.15 Hz (high pass) and 105 Hz (low pass) hardware filters. Study approval by the Ethics in Human Research Committees of the Royal Children’s Hospital and St. Vincent’s Hospital, Melbourne. Informed consent was obtained in all cases. 2.2. EEG source localization (ESL) EEG recordings were uploaded to Scan 4.3Ò (NeuroscanÒ, El Paso, Texas, USA) for analysis (common average reference at vertex). Single IEDs were epoched from 200 ms to +500 ms relative to the spike peak and uploaded to CURRY 5.0Ò (CompumedicsÒ, Melbourne, Australia). Notch 50 Hz and bandpass 0.5–70 Hz filters were applied. Mean and incremental (every four milliseconds) SNR calculations for the spike interval were based on the pre-spike interval (identical number of noise sampling points per ESL operation). Electrode positions were label match co-registered using predefined electrode locations for the three Montreal Neurological Institute (MNI)-based realistic models. The interval for ESL analysis was marked from spike onset-to-peak latency using a butterfly plot (see Supplementary Figures S1 and S2). Spike onset was defined by the first significant deflection (surface positive or negative) from baseline and spike peak was defined by the peak (surface negative) amplitude of the last overlapping spike waveform on the butterfly plot. Mean global field orthogonal signal components of the spike with an SNR >1.0 identified by principal component analysis (PCA) underwent independent component analysis (ICA). The five single spikes for each patient were later electrographically averaged based on automated detection and overlay of surface-negative maxima. Averaged IEDs were subjected to the same ESL operation as described for single IEDs. A total of 48 forward–inverse–subspace modelling set-ups (Fig. 1) were derived from three forward models (Fuchs et al., 2002) – standardized boundary element method (BEMs), leadfield-interpolated boundary element method (BEMi), leadfield-interpolated finite element method (FEMi); four inverse models – minimum norm least squares (MNLS) (Fuchs et al., 1999), minimum L1 norm (L1) (Wagner et al., 1998), low resolution electromagnetic tomography (LORETA) (PascualMarqui et al., 1994), standardized LORETA (sLORETA) (PascualMarqui, 2002); four subspace constraints – whole volume (3D), cortex with rotating sources (CxR), cortex with fixed sources normal to cortex (CxN), cortex with fixed extended sources using a 20 millimetre patch (patch) (Wagner et al., 2001).

2. Methods 2.3. Results analysis (quantitative) 2.1. Patients and EEG recordings Eight patients underwent interictal scalp EEG – four with benign focal epilepsy of childhood with centrotemporal spikes (BFEC), 6–10 years (mean 8.5 years); four with unilateral mesial

Each CDR map was quantified by seven output parameters describing the spatiotemporal behaviour of that mini-dipole (as one of several thousand mini-dipoles per CDR map) that reached the highest current density (cdmax) at some time-point during

1728

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Fig. 1. Outline showing how the 48 distributed modelling approaches were derived from the three separate modelling variables (at left) and how the quantitative and qualitative criteria relate to each other (at right). For illustrative purposes, the cdmax unitary ‘mini-dipole’ (amplitude or strength in micro-Amperes per millimetre squared), upon which the location and orientation rely, is marked in blue. The cdnet (current density in micro-Amperes per millimetre cubed) relates to the entire ‘mini-dipole’ distributed solution (or network) for a given forward–inverse–subspace ESL operation. The fwhm and rd are based on the cdnet, and the snr is a reflection of the latency at which the best-fit CDR solution is calculated. Corresponding qualitative analyses are based on the colour-scaled cdnet distribution against the MNI cortical surface or the 3D grid. Abbreviations: cd, current density; rd, residual deviation; fwhm, full width half maximum; snr, signal to noise ratio; MNI, Montreal Neurological Institute averaged MRI brain; 3D, whole grid space.

the IED epoch (best-fit): current density (cd, as micro-Amperes per millimetre squared for MNLS, L1, LORETA; F-distribution probability function for sLORETA), three location parameters (x, y, and z axes, millimetres), three orientation parameters (nx, ny, and nz vectors, fractions of ±1.0). Three additional parameters quantified net CDR map behaviour at the same cdmax latency (cdnet, as microAmperes per millimetre cubed for MNLS, L1, and LORETA; F-distribution for sLORETA): SNR, full width half maximum (fwhm, millilitres), residual deviation (rd) (Fig. 1). The cdnet was therefore the quantitative expression of the CDR map in that it was the calculated sum total current density of the solution for the pre-defined subspace, giving it cubic dimensions. (The cdmax carried squared dimensions as it represented the single mini-dipole fit that held the highest strength within the total cdnet subspace). The rd constituted the percentage of the electrode-recorded signal that was unaccounted for by the cdnet (viz. CDR map). That is, an rd of 15% would indicate that the CDR map does not explain 15% of the measured data seen by the electrodes at the given time-point. Two-way analysis of variance (ANOVA) and general linear method (GLM) sub-analysis stratified the relative impact of each modelling variable on CDR outcome.

2.4. Qualitative Comparisons between modelling approaches were based on the CDR map corresponding to each best-fit distributed modelling solution. Maps were colour-scaled to reflect the current density gradient. Lowest (black) and lower (red) scale signal represented activity removed from the CDR maximum. Higher (yellow) and highest (white) scale signal represented activity close to or at the CDR maximum. To avoid bias, map displays were not thresholdset. All upper scale (yellow–white) signal, irrespective of its final distribution (hemispheric, lobar, sublobar) was considered proximate to the source so modelled by a given forward–inverse–subspace set-up (Fig. 1). It follows then that final results analysis for every forward–inverse–subspace modelling operation was based on the quantitative data in the context of the qualitative data. More specifically, statistical significance for any forward, inverse, or subspace modelling effect on a given ESL attribute (location, orientation, strength, width/volume, variance) was always interpreted relative to the nature of that difference given by the corresponding CDR map. This was to avoid the problems met by previous studies where excessive emphasis was placed on millimetre scale differences that might meet statistical significance but that, in practice, are unlikely

to meet the same degree of clinical significance (Plummer et al., 2008) – the localization afforded by scalp recorded EEG is on a centimetre scale at best, corresponding to the sublobar scale anatomically (Ebersole, 2003). CDR solution ‘focality’, which was characterized quantitatively by the fwhm in respect of the cdnet, was similarly analysed. Highly ‘focal’ quantitative solutions (high cdnet but narrow/low volume fwhm), while reflecting higher ‘localizability’ of the source, may prove difficult to reconcile with the currently understood spatiotemporal behaviour of cortical sources in MTLE and BFEC – if well removed anatomically from the anticipated source configuration (e.g., to the contralateral occipital lobe in typical BFEC or to the lower brainstem in MTLE). 2.5. Reliability Intra-operator and inter-operator levels of concordance were assessed by the following method. Operator A selected 10 spikes (one spike from three patients and two spikes from the fourth patient from each group) for modelling. Each spike was modelled by Operator A twice (20 operations). Operator B independently performed a single trial on each spike (10 operations). For each ESL operation, operator A pre-selected the following criteria: patient IED, forward model, inverse model, and subspace constraint. The following ESL steps were therefore freely chosen by both operators: pre-spike interval (effectively the number of noise sampling points) for snr calculation, spike segment for analysis (with the aim of selecting spike onset-to-peak interval), and the number of allowable PCA and subsequent ICA components used to guide ESL. Measures of intra- and inter-operator concordance were based on the Fleiss scale as follows, kappa > +0.75: excellent, +0.4 to +0.75: fair, <0.4: poor (Fleiss, 1981). 2.6. Software factors Computation times and interface problems were logged. All ESL operations were carried out on a standard personal computer with an Intel Pentium 4Ò processor (2 GB RAM) running under a Windows XP Professional operating system. 3. Results Best-fit CDR parameters (from 19,200 results) are shown as Supplementary data in Supplementary Table S1 (BFEC) and Supplementary Table S2 (MTLE). Table 1 shows the forward model impact for a given inverse–subspace set-up. From 160 ANOVA comparisons

1729

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Table 1 Effect of forward model on best-fit CDR output values for a given inverse–subspace modelling set-up in BFEC (upper panel) and MTLE (lower panel). Inter-patient variance has been accounted for in each case by two-way ANOVA. Calculations are based on the datasets displayed in the Supplementary Table S1 (BFEC) and the Supplementary Table S2 (MTLE). Probability values <0.05 are given in bold. p Values could not be reliably calculated for the fwhm parameter in MTLE for L1–CxR and L1–CxN due to the minute fwhm means (order of 10 20) for these L1–subspace set-ups. Abbreviations: L1, minimum L1 norm; LORETA, low resolution electromagnetic tomography; MNLS, minimum norm least squares; sLORETA, standardized LORETA; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; MTLE, mesial temporal lobe epilepsy; 3D, whole grid space; CxR, cortex with rotating sources; CxN, cortex with fixed sources; patch, cortex with fixed extended sources; cd, current density maximum; fwhm, full width half maximum; x, y, and z, location axes; nx, ny, and nz, orientation vectors; snr, signal to noise ratio; rd, residual deviation; CDR, current density reconstruction; ANOVA, analysis of variance. L1

LORETA

3D BFEC cd fwhm X y z nx ny nz snr rd

0.24 0.92 0.59 0.19 0.47 0.027 0.011 0.35 0.99 0.99

MTLE cd fwhm X y z nx ny nz snr rd

0.19 0.15 0.51 0.56 0.53 0.59 0.77 0.65 0.89 0.66

CxR

CxN

patch

MNLS

sLORETA

3D

CxR

CxN

patch

3D

CxR

CxN

patch

3D

CxR

CxN

patch

0.08 0.17 0.68 0.27 0.80 0.91 0.50 0.75 0.98 0.86

0.008 0.33 0.66 0.23 0.84 0.17 0.94 0.72 0.99 0.88

0.002 0.84 0.88 0.06 0.94 0.71 0.93 0.56 0.99 0.96

0.09 0.94 0.99 0.91 0.91 0.98 0.98 0.92 0.99 0.99

<0.001 <0.001 0.40 <0.001 <0.001 0.65 0.84 0.43 0.99 0.69

<0.001 0.10 0.99 0.87 0.94 0.92 0.88 0.95 0.99 0.97

<0.001 0.78 0.45 0.49 0.90 0.41 0.47 0.89 1.00 0.96

<0.001 0.66 0.63 0.16 0.64 0.48 0.37 0.09 0.99 0.99

0.002 0.011 0.79 0.43 0.50 0.30 0.91 0.93 0.99 0.99

<0.001 0.044 0.76 0.13 0.97 0.78 0.06 0.91 0.99 0.99

<0.001 0.92 0.99 0.93 0.55 0.87 0.63 0.71 0.99 0.99

0.80 0.001 0.90 0.98 0.09 0.59 0.74 0.041 0.99 0.99

0.69 0.46 0.98 0.72 0.86 0.33 0.48 0.78 0.99 0.99

0.35 0.38 0.73 0.96 0.87 0.67 0.87 0.24 0.33 0.89

0.69 0.46 0.16 0.11 0.13 0.09 0.97 0.047 0.99 0.99

<0.001 – 0.37 0.76 0.53 0.68 0.92 0.43 0.97 0.91

<0.001 – 0.12 0.66 0.61 0.38 0.87 0.60 0.95 0.65

<0.001 0.23 0.09 0.61 0.50 0.55 0.69 0.89 0.93 0.49

0.17 0.80 0.90 0.85 0.82 0.97 0.90 0.88 0.99 0.99

0.008 0.45 0.74 0.63 0.016 0.13 0.97 0.77 0.99 0.54

0.10 <0.001 0.69 0.98 0.81 0.92 0.99 0.92 0.99 0.95

0.025 0.870 0.98 0.99 0.98 0.91 0.12 0.98 0.97 0.63

0.09 0.36 0.89 0.92 0.85 0.56 0.24 0.45 0.99 0.99

0.54 0.010 0.78 0.98 0.52 0.82 0.98 0.35 0.99 0.98

0.06 <0.0001 0.16 0.98 0.43 0.90 0.49 0.78 0.99 0.93

0.012 0.99 0.79 0.91 0.99 0.98 0.07 0.98 0.99 0.95

0.83 0.06 0.63 0.92 0.72 0.039 0.84 0.79 0.99 0.99

0.97 0.62 0.99 0.78 0.76 0.33 0.82 0.70 1.00 0.99

0.99 0.17 0.96 0.94 0.95 0.89 0.67 0.21 1.00 0.99

0.99 0.75 0.72 0.99 0.94 0.53 0.92 0.99 1.00 0.99

per group, 19 (12%) in BFEC and 11 (7%) in MTLE reached p < 0.05. For the inverse model impact for a given forward–subspace set-up (Table 2), 42 (35%) in BFEC and 49 (41%) in MTLE met p < 0.05 (120 ANOVAs). For the subspace constraint impact for a given forward–inverse set-up (Table 3), 66 (55%) in BFEC and 25 (21%) in MTLE met p < 0.05 (120 ANOVAs).

3.1. Forward model impact (BEMi, BEMs, FEMi) 3.1.1. cdmax, cdnet, fwhm In BFEC, cdmax means were 1.5–3.0-fold lower for FEMi versus BEMi, BEMs (BEM) (Supplementary Table S1, Table 1) as also reflected in the corresponding CDR map current density (cdnet)

Table 2 Effect of inverse model on best-fit CDR output parameters for a given forward–subspace modelling set-up in BFEC (upper panel) and MTLE (lower panel). Probability values <0.05 are given in bold. Abbreviations: BEMi, leadfield-interpolated boundary element method; BEMs, standardized boundary element method; FEMi, leadfield-interpolated finite element method; 3D, whole grid space; CxR, cortex with rotating sources; CxN, cortex with fixed sources; patch, cortex with fixed extended sources; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; MTLE, mesial temporal lobe epilepsy; cd, current density; fwhm, full width half maximum; x, y, and z, location axes; nx, ny, and nz, orientation vectors; snr, signal to noise ratio; rd, residual deviation; CDR, current density reconstruction. BEMi

BEMs

FEMi

3D

CxR

CxN

patch

3D

CxR

CxN

patch

3D

CxR

CxN

patch

BFEC cd fwhm x y z nx ny nz snr rd

<0.001 0.08 0.08 0.66 <0.001 0.64 0.71 0.52 0.99 0.026

<0.001 0.018 0.09 <0.001 <0.001 0.69 0.002 0.013 0.99 0.70

<0.001 0.042 0.56 0.25 0.51 0.69 0.28 0.16 0.99 0.83

<0.001 0.038 0.18 0.80 0.61 0.50 0.69 0.54 0.99 0.76

<0.001 0.16 0.11 0.66 <0.001 0.27 0.80 0.79 0.99 0.020

<0.001 0.031 0.10 0.025 <0.001 0.99 <0.001 0.56 0.99 0.91

<0.001 0.06 0.67 0.22 0.71 0.70 0.30 0.35 0.99 0.89

<0.001 0.034 0.19 0.89 0.67 0.45 0.81 0.72 0.99 0.93

<0.001 <0.001 0.12 0.013 <0.001 0.06 0.014 0.001 0.99 0.023

<0.001 <0.001 0.026 <0.001 <0.001 0.84 <0.001 0.39 0.99 0.19

<0.001 0.60 0.19 0.22 0.48 0.06 0.27 <0.001 0.41 0.51

<0.001 0.041 0.10 0.41 0.10 0.07 0.78 0.044 0.99 0.75

MTLE cd fwhm x y z nx ny nz snr rd

0.007 <0.001 0.91 0.67 <0.001 0.69 0.25 <0.001 0.98 0.007

0.008 <0.001 0.96 0.82 0.001 0.09 0.24 0.08 0.99 0.001

0.011 <0.001 0.86 0.45 0.040 0.11 0.61 0.022 0.99 0.09

0.009 <0.001 0.52 0.70 0.009 0.71 0.58 0.09 0.99 0.18

0.007 <0.001 0.79 0.99 <0.001 0.99 0.59 0.002 0.98 0.010

0.007 <0.001 0.94 0.93 <0.001 0.10 0.20 0.24 0.99 0.002

0.011 <0.001 0.89 0.27 0.013 0.61 0.60 0.006 0.99 0.044

0.008 <0.001 0.61 0.63 0.015 0.32 0.63 0.29 0.99 0.12

0.012 <0.001 0.47 0.92 <0.001 0.07 0.08 <0.001 0.99 0.07

0.005 <0.001 0.76 0.56 0.040 0.90 0.65 0.029 0.97 <0.001

0.006 <0.001 0.62 0.69 0.003 0.67 0.08 0.67 0.99 0.007

0.004 <0.001 0.28 0.87 0.008 0.27 0.08 0.20 0.99 0.004

1730

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Table 3 Effect of subspace constraint on best-fit CDR output parameters for a given forward–inverse modelling set-up in BFEC (upper panel) and MTLE (lower panel). Probability values <0.05 are given in bold. Abbreviations: L1, minimum L1 norm; LORETA, low resolution electromagnetic tomography; MNLS, minimum norm least squares; sLORETA, standardized LORETA; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; MTLE, mesial temporal lobe epilepsy; BEMi, leadfield-interpolated boundary element method; BEMs, standardized boundary element method; FEMi, leadfield-interpolated finite element method; cd, current density; fwhm, full width half maximum; x, y, and z, location axes; nx, ny, and nz, orientation vectors; snr, signal to noise ratio; rd, residual deviation; CDR, current density reconstruction. L1

LORETA

MNLS

sLORETA

BEMi

BEMs

FEMi

BEMi

BEMs

FEMi

BEMi

BEMs

FEMi

BEMi

BEMs

FEMi

BFEC cd fwhm x y z nx ny nz snr rd

0.003 <0.001 0.001 0.13 0.048 0.64 0.014 0.040 1.00 0.99

0.005 <0.001 0.002 0.06 0.54 0.61 0.020 0.003 0.99 0.99

0.002 <0.001 0.027 0.64 0.39 0.16 <0.001 <0.001 0.99 0.93

0.029 <0.001 0.35 0.024 <0.001 0.78 0.036 0.012 1.00 0.34

0.033 <0.001 0.74 0.06 <0.001 0.95 0.006 0.16 0.99 0.19

<0.001 <0.001 0.56 0.001 0.035 0.52 0.040 0.022 0.99 0.23

<0.001 <0.001 <0.001 0.025 0.39 0.74 0.030 0.010 0.99 0.99

<0.001 <0.001 <0.001 0.020 0.45 0.76 0.021 0.006 0.99 0.99

<0.001 <0.001 <0.001 0.74 0.82 0.19 <0.001 <0.001 0.99 0.99

0.08 <0.001 <0.001 0.047 <0.001 0.41 0.019 0.023 1.00 0.99

0.05 <0.001 <0.001 0.034 <0.001 0.48 0.014 0.016 1.00 0.99

0.22 <0.001 <0.001 0.019 <0.001 0.61 0.15 0.018 0.43 0.96

MTLE cd fwhm x y z nx ny nz snr rd

<0.001 <0.001 0.81 0.57 0.67 0.52 0.31 0.77 0.97 0.58

<0.001 <0.001 0.20 0.98 0.56 0.93 0.64 0.77 0.99 0.82

<0.001 <0.001 0.24 0.65 0.07 0.87 0.79 0.34 0.97 0.74

0.93 <0.001 0.74 0.23 <0.001 0.54 0.76 0.07 0.99 0.40

0.89 <0.001 0.31 0.16 <0.001 0.31 0.61 0.14 0.99 0.42

0.12 <0.001 0.54 0.10 <0.001 0.62 0.78 0.17 0.99 0.30

0.033 <0.001 0.51 0.88 0.47 0.73 0.65 0.07 0.99 0.99

0.029 <0.001 0.26 0.82 0.70 0.92 0.49 0.05 0.99 0.99

0.07 <0.001 0.007 0.36 0.60 0.47 0.42 0.56 1.00 0.99

0.81 <0.001 0.005 0.41 0.66 0.006 0.35 0.37 1.00 0.99

0.73 <0.001 0.010 0.48 0.84 0.028 0.47 0.47 1.00 0.99

0.97 <0.001 0.06 0.82 0.70 0.26 0.98 0.56 0.99 1.00

Fig. 2. Series of CDR maps taken from BFEC patient 2 (Supplementary Figure S1 IED, averaged) for subspace 3D. Note that CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 8.598 ms). Abbreviations: IED, interictal epileptiform discharge; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; CDR, current density reconstruction; 3D, whole grid space.

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

1731

Fig. 3. Series of CDR maps taken from BFEC patient 2 (Supplementary Figure S1 IED, averaged) for subspace CxR. Note that CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 8.598 ms). Abbreviations: IED, interictal epileptiform discharge; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; CDR, current density reconstruction; CxR, cortex with rotating sources.

scales (Figs. 2–5). FEMi fwhm means were 1.5–2.0-fold narrower for LORETA–CxR, MNLS–CxR, MNLS–CxN, yet 1.5-fold wider for sLORETA–3D (Figs. 3 and 4). In MTLE, cdmax means were 2.0–3.0-fold lower for the FEMi versus the BEM set-ups (Supplementary Table S2, Table 1). Differences were striking for the L1–subspace (Figs. 6–9). No differences were seen for the sLORETA–subspace. The fwhm means were 1.5– 2.0-fold narrower for FEMi with LORETA–CxN, MNLS–CxR, MNLS–CxN (Figs. 8 and 9) but corresponding CDR maps were highly variable (cdnet shifted to lateral left temporal lobe in MTLE patient 2 (MNLS–CxR), towards a different lobe altogether (ipsilateral frontal) in MTLE patient 4 (LORETA–CxN), and towards a complete loss of signal definition in MTLE patient 1 (MNLS–CxN). 3.1.2. Location (x, y, z) and orientation (nx, ny, nz) In BFEC, with one exception (LORETA–CxR), there was no forward modelling effect on location (Table 1). This was evident in the consistent location-wise CDR mapping by forward model. For LORETA–CxR, FEMi solutions were in the order of 35 mm posterior and 40 mm rostral to corresponding BEM locations (Supplementary Table S1). However, LORETA–CxR maps gave ill defined, scattered cdnet locations (Figs. 2–5). Orientation differences were limited to L1–3D (nx, ny) and sLORETA–3D, –patch (nz) where FEMi and BEM means carried opposing polarities; however, p values just met significance and standard deviations were wide (Supplementary Table S1, Table 1). In MTLE, a forward modelling influence on location was again restricted to LORETA–CxR (z parameter) (Table 1). The FEMi z mean was 8.8 mm and 9.6 mm rostral to the corresponding BEMi and

BEMs results. Qualitatively, the FEMi was associated with a signal shift from mesial temporal to frontobasal cortex (MTLE patients 2 and 3) and from anterior temporal tip to lateral frontal convexity (MTLE patient 4) (Fig. 7). For orientation, only sLORETA–3D saw a forward effect (nx) but mean polarities remained concordant across the three forward models (Supplementary Table S2, Table 1). 3.1.3. snr, rd For both groups, snr and rd were not influenced by forward model selection (Supplementary Tables S1 and S2, Table 1). 3.2. Inverse model impact (L1, LORETA, MNLS, sLORETA) 3.2.1. cdmax, cdnet, fwhm In BFEC, cdmax was influenced by inverse model for 12 forward–subspace set-ups (Table 2). LORETA and MNLS cdmax means were of the same magnitude, while L1 cdmax means were several hundred fold greater for BEM-based set-ups (Supplementary Table S1). Values for fwhm were generally widest for LORETA but narrowest for L1 and sLORETA. Qualitatively, LORETA and MNLS maps for CxN (Fig. 4) and patch (Fig. 5) were similar. For 3D, MNLS maps were less well defined (than LORETA) at the outer grid surface (Fig. 2). For CxR, LORETA and MNLS maps were heterogeneous and tended to cross lobes, while sLORETA maps were more contiguous, localizing to the lateral frontal lobe and antero-basal aspects of the parietal lobe rostral to Sylvian fissure (Fig. 3). L1 gave pinpoint localization (–3D, –CxR, –CxN, less so for –patch) (Figs. 2–5). In MTLE, cdmax and fwhm were heavily influenced by inverse model (principally L1) for all forward–subspace set-ups (Table 2).

1732

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Fig. 4. Series of CDR maps taken from BFEC patient 2 (Supplementary Figure S1 IED, averaged) for subspace CxN. Note that CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 8.598 ms). Abbreviations: IED, interictal epileptiform discharge; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; CDR, current density reconstruction; CxN, cortex with fixed sources.

No consistent differences between LORETA and MNLS were seen. L1 activation was discrete and difficult to perceive against the MNI cortical surface for CxR (Fig. 7) and CxN (Fig. 8). L1 fwhm means were narrowest, LORETA fwhm means were widest for BEM–3D, but sLORETA means were widest for FEMi–3D and for all forward–CxR set-ups (Supplementary Table S2). 3.2.2. Location (x, y, z) and orientation (nx, ny, nz) In BFEC, y locations were influenced by inverse model for all forward–CxR set-ups (Table 2); differences were attributable to LORETA where y locations were 15–20 mm anterior-most for BEMi/s– CxR in conjunction with a positive ny mean polarity (posteriorward vector) (Supplementary Table S1). For forward–3D and – CxR, z locations were most rostral for sLORETA (nz vector inferior-ward) and most caudal for LORETA (nz vector superior-ward) (Supplementary Table S1, Table 2). Qualitatively, LORETA–CxR locations were discordant, taking in brainstem, cerebellum and occipital cortices; the sLORETA cdnet distribution straddled mid to lower Rolandic sulcus, abutting the Sylvian fissure (Fig. 3). In MTLE, z location was influenced by inverse model for all forward–subspace set-ups (Table 2). The most caudal positions were consistently given by sLORETA (Supplementary Table S2) with differences greatest against LORETA for 3D and against L1 for CxR. Vector nz was consistently negative-most (inferior-ward) for sLORETA. Qualitatively, sLORETA–3D maps delineated grid points at the basal temporal surface (LORETA delineated superior temporal and basal frontal surfaces) (Fig. 6). For CxR, sLORETA favoured a source encompassing the basolateral temporal surface and either the anterior temporal tip (Patients 1 and 3) or the infero-basal temporal surface (Patients 2 and 4). MNLS favoured frontal lobe con-

vexity in each case (Fig. 7). For patch, the L1 cdnet (often scattered and multifocal) favoured the supero-lateral temporal and/or lateral margin of the ipsilateral frontal lobe (Fig. 9). 3.2.3. snr, rd The snr results in BFEC and MTLE were unaffected by inverse model. The rd was only affected at the 3D level in BFEC (Table 2) for which LORETA returned the lowest percentages (around 18%, other models around 24%) (Supplementary Table S1). In MTLE, rd values were around 10% higher for L1 against the other inverse algorithms (Supplementary Table S2 and Table 2). 3.3. Subspace constraint impact (3D, CxR, CxN, patch) 3.3.1. cdmax, cdnet, fwhm In BFEC, patch gave the lowest cdmax means for all forward–L1 and forward–LORETA set-ups (Supplementary Table S1, Table 3). CxR gave the lowest cdmax when MNLS was applied. There was little variation in the sLORETA F-distribution map across the four subspaces (Figs. 2–5). By comparison, there was a 50–100 fold reduction in the L1 cdnet scale for patch (and 3D) against CxR and CxN. The fwhm parameter was influenced by subspace across all forward–inverse set-ups (Table 3). CxN fwhm means were narrowest and 3D widest. In MTLE, patch reduced L1 cdmax relative to the other subspaces. CxR again gave the lowest cdmax values when MLNS was used (Supplementary Table S2, Table 3). There was a striking similarity in the sLORETA peak cdnet between subspaces for a given forward model (Figs. 6–9). For LORETA, despite the comparable cdmax between subspaces, cdnet results for 3D were between 50

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

1733

Fig. 5. Series of CDR maps taken from BFEC patient 2 (Supplementary Figure S1 IED, averaged) for subspace patch. Note that CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 8.598 ms). Abbreviations: IED, interictal epileptiform discharge; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; CDR, current density reconstruction; patch, cortex with fixed extended sources.

and 100 times lower compared to the CxR, CxN, and patch. 3D fwhm means for LORETA, MNLS, and sLORETA were 100 times wider than CxN and patch means. Forward–L1 fwhm results were especially narrow with means ranging from <0.1 ml (CxR, CxN) to 0.5 ml (3D) (Supplementary Table S2).

3.3.2. Location (x, y, z) and orientation (nx, ny, nz) In BFEC, subspace influenced x location parameter for the nonLORETA approaches (Table 3) – CxN/patch x means differed from 3D/CxR means by 15–25 mm (Supplementary Table S1, Table 3). Relative to CxR and 3D, y mean locations were more anterior with CxN (forward–sLORETA, BEMi/s–MNLS) and with patch (FEMi/ BEMi–LORETA). For the z parameter, CxR led to the most caudal positions for all forward–LORETA set-ups (Supplementary Table S1) where the ny vector was directed posterior-ward. A LORETA–subspace switch from CxR to patch produced a rostral shift in CDR location from the undersurface of the MNI brain to an area proximate to the Sylvian fissure. The converse was seen for forward–sLORETA, where CxR and 3D generated the most rostral z locations. For nz, the most consistent subspace effect was the positive nz polarity for patch (superior-ward projection) versus the negative polarities for 3D and CxR (inferior-ward projection). The most noteworthy qualitative effect of the above subspace-directed location and orientation shifts was for the patch constraint to ‘pull’ the CxR map towards the Sylvian fissure (along a line approximating the Rolandic fissure) or towards the ipsilateral temporal lobe (when the corresponding CxR solution map abutted the upper margin of the Sylvian fissure) (Figs. 3 and 5).

In MTLE, x parameter shifts were driven by patch (Supplementary Table S2 and Table 3). A lateral shift in the peak sLORETA Fdistribution occurred when either CxN or patch replaced CxR as the subspace constraint (Figs. 7–9). Accordingly, the nx vector was the only orientation parameter ever affected by subspace. This effect was only limited to BEMi/s–sLORETA (Supplementary Table S2 and Table 3) but there was polarity concordance. Location parameter z was influenced by subspace across three forward– LORETA set-ups – CxR locations were most caudal and 3D most rostral. Qualitatively, this was mirrored by patch and CxR which produced a caudal shift relative to the 3D subspace localization (Figs. 6–9). Equally notable was the localization afforded by patch relative to the corresponding CxR solution; patch solutions effected a ‘contraction’ of the CxR subspace solution to the basal and/or basolateral temporal lobe surface (Figs. 7 and 9). However, patch solutions were at times accompanied by an island of proximatescale CDR activity encompassing the ipsilateral orbito-frontal cortical surface (Fig. 9). In such cases, the CDR maximum fell within the boundaries of the temporal lobe activation, not the frontal lobe. 3.3.3. snr, rd For the two groups, best-fit snr and the rd parameters were not significantly influenced by subspace (Supplementary Tables S1 and S2, Table 3). 3.4. Reliability For all best-fit CDR output parameters, both intra (A1–A2) and inter (A1–B1, A2–B1)-operator concordance measures (Table 4)

1734

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Fig. 6. Series of CDR maps taken from MTLE patient 4 (Supplementary Figure S2 IED, averaged) for subspace 3D. CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 0.78 ms). Abbreviations: IED, interictal epileptiform discharge; MTLE, mesial temporal lobe epilepsy; CDR, current density reconstruction; 3D, whole grid space.

were generally high with 85% of paired comparisons (28/33) returning kappa > 0.75 (Fleiss rating ‘excellent’). ‘Poor’ inter-operator concordance for the y parameter was best explained by differences for BEMs–MNLS–CxN (64.8 mm), BEMs– MNLS–3D (14 mm), BEMi–L1–CxR (18.4 mm) as distinct from BEMi–sLORETA–CxN (0 mm), FEMi–sLORETA–patch (3 mm), BEMs–LORETA–3D (7 mm), FEMi–LORETA–patch (3.8 mm). ‘Fair’ intra-operator concordance for y was best explained by BEMs– MNLS–3D (14 mm), FEMi–L1–CxN (63.9 mm), BEMi–L1–CxR (16.9 mm) as opposed to BEMi–sLORETA–CxN (0 mm), FEMi–sLORETA–patch (0 mm), BEMs–LORETA–3D (0 mm), FEMi–LORETA– patch (2.2 mm). ‘Fair’ inter-operator concordance for ny (A1–B1) stemmed from MNLS (only inverse model to return opposite polarity) and L1 (most discordant ‘same-polarity’ vector result, 0.28). ‘Fair’ inter-operator concordance for nz (A1–B1, A2–B1) stemmed from L1 (polarity opposite) and MNLS (most discordant ‘samepolarity’, 0.53). By comparison, the most concordant inter-operator vector results arose with sLORETA: ny A1 0.17, B1 0.17, nz A1 0.02, A2 0.02, B1 0.02 (BEMi–sLORETA–CxN), nz A1 0.78, A2 0.78, B1 0.76 (FEMi–sLORETA–patch). 3.5. Software factors Using an Intel Pentium 4Ò processor (2 GB RAM), CDR computation times, CPTs, (Table 5) were sub-second for all forward–subspace set-ups involving sLORETA and MNLS. L1 inverse modelling was next (0.8–9.6 s in BFEC, 1.0–18.7 s in MTLE). Mean L1 CPTs were always longer for MTLE (versus BFEC) and for CxR (versus other subspaces). LORETA CPTs were the longest. Multiple iterations (two to six) were needed, each requiring at least one minute

for 3D and CxR and up to 30 s for CxN and patch. Like L1, LORETA took longer to arrive at a best-fit CDR solution in MTLE and for CxR (mean 7.7 min, range 2.7–17.4 min). CPTs did not differ between forward models. When Profusion EEG files were uploaded to Scan 4.3 for analysis, inversion of the waveform polarity occurred across all channels. This problem was remedied with the Scan montage editor using a linear inversion matrix. No polarity inversion occurred when the corrected Scan file was then uploaded to Curry.

4. Discussion This study has undertaken an analysis of 48 different distributed modelling approaches as applied to the EEG source localization (ESL) of scalp-recorded interictal discharges (IEDs) in two common types of focal epilepsy- the Rolandic variety of benign focal epilepsy of childhood (BFEC), and mesial temporal lobe epilepsy (MTLE) secondary to hippocampal sclerosis. The aim was to identify an optimal distributed modelling set-up by characterizing (quantitatively and qualitatively) the impact of each modelling variable (forward, inverse, and subspace) on the nature of the best-fit solution. Expressed as the proportion of two-way ANOVA comparisons reaching p < 0.05, the forward model was seen to exert the least effect on the best-fit CDR result (BFEC 12%, MTLE 7%), while the inverse and subspace modelling conditions had a greater impact on CDR outcome (BFEC 35%, MTLE 41% for inverse model; BFEC 55%, MTLE 21% for subspace constraint). Such disproportionate effects on CDR outcome were also evidenced qualitatively by way of the changes seen in the CDR display map across the 48 modelling set-ups.

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

1735

Fig. 7. Series of CDR maps taken from MTLE patient 4 (Supplementary Figure S2 IED, averaged) for subspace CxR. CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 0.78 ms). Abbreviations: IED, interictal epileptiform discharge; MTLE, mesial temporal lobe epilepsy; CDR, current density reconstruction; CxR, cortex with rotating sources.

4.1. Forward model The FEMi model accounted for all significant variations in the final fit characteristics for both patient groups, with the BEMs and BEMi models returning very similar results. CDR parameters were, however, infrequently affected and the size of the impact was, in each case, minor in quantitative and qualitative terms. The main effect was seen for cdmax which was lowered by the FEMi. This is probably attributable to the preferential conductivity weighting given by the FEMi to tangentially oriented fields (Fuchs et al., 2007); hence the effect was more often encountered in the BFEC group, a condition classically associated with tangential fields generated by a horizontal dipolar source on the bank of the lower Rolandic sulcus (Wong, 1998). The absence of a forward modelling effect at any level on the sLORETA probability distribution in both patient groups is interesting and suggests that it is the approximation of a realistic head shape rather than any subtler difference between these forward models that underlies the cdmax and cdnet result here. The important upshot is that owing to the present impracticality of generating individualized FEMi forward models by virtue of the heavy computational load, individualizing the forward model via application of the more standardized BEM, or BEMs, to the patient’s MRI scan does not appear to come at a cost in terms of the quality of the final fit, particularly when sLORETA is applied as the inverse model. This result is reminiscent of the findings noted for dipole modelling, in which the fundamental impact of the forward model on the fit quality relies on the shape of the volume conductor (realistic head versus spherical shell) rather than on any subtler difference between realistic models (Plummer

et al., 2008). It remains to be seen whether these observations extend to the ‘more sophisticated’ FEMi models currently being developed which factor into their conductivity set-up white matter anisotropy in addition to skull layer anisotropy (Haueisen et al., 2002). As reflected by the BFEC and MTLE CDR maps, the forward model had minimal impact on best-fit x, y, z location. Indeed, LORETA– CxR was the only set-up to see a significant forward modelling effect quantitatively for both patient groups, but here, maps were heterogeneous and indistinct (Figs. 3 and 7). Simulation studies have also noted this tendency for LORETA to generate spatially blurred CDR maps (Fuchs et al., 1999). A forward modelling effect on best-fit CDR orientation was difficult to characterize. This is not surprising given the millisecond speed with which dipolar and distributed models may undergo a change in their orientation (Ebersole, 2003) in order to ‘keep pace with’ (or to satisfactorily explain) the fluxing surface EEG field topography typically seen in IEDs. Therefore, if orientation-wise comparisons between different CDR models are to carry any validity, the same point or phase of the spike wave complex in BFEC and MTLE should be interrogated. This condition was satisfied in the present study as the peak of the major negative spike in both the BFEC IED (surface negative temporal field, surface positive frontal field) and the MTLE IED (surface negative temporal field) consistently coincided with the best-fit CDR latency by virtue of the peaking signal-to-noise-ratio (SNR) during this phase of the spike. Despite this, the variability in CDR cdmax orientation was high both across and within patients as evidenced by the relatively large standard deviations for all three orientation vectors (nx, ny, and

1736

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Fig. 8. Series of CDR maps taken from MTLE patient 4 (Supplementary Figure S2 IED, averaged) for subspace CxN. CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 0.78 ms). Abbreviations: IED, interictal epileptiform discharge; MTLE, mesial temporal lobe epilepsy; CDR, current density reconstruction; CxN, cortex with fixed sources.

nz). The absence of a significant quantitative difference between the three forward models for the best-fit explained signal variance (residual deviation, rd) and the best-fit signal-to-noise ratio (snr) (Table 1) gives further credence to the overall impression that the forward model has minimal impact on the nature of the best distributed fit in BFEC and MTLE. 4.2. Inverse model The inverse model had a major effect on the distributed ESL result on several fronts in both patient groups. In BFEC, the sLORETA algorithm generated the most consistent source localization results (for orientation and location) that were most compatible with a putative focus at the lower Rolandic sulcus. In MTLE, sLORETA also generated the most consistent ESL results for location and orientation that would be compatible with the well studied electrophysiological behaviour of a propagating source of mesial temporal lobe origin towards the ipsilateral basal and baso-lateral surfaces (Ebersole, 1999, 2003; Ebersole and Hawes-Ebersole, 2007). Of the four inverse algorithms, L1 generated the most ‘seductive’ focal solutions, but at a cost. L1 locations were prone to scatter and the capacity of the algorithm to explain the measured data was found wanting in MTLE, when it returned the highest residual deviation (rd) results. LORETA performed well with the 3D subspace constraint but poorly under the more anatomically informed conditions of the CxR constraint. MNLS solutions were particularly ‘weak’ for current density per se and proved difficult to resolve qualitatively for both the 3D and CxR subspace constraints in both patient groups.

4.3. Subspace constraint For both patient groups, it was reassuring to find that the unexplained variance (rd) did not differ between the 3D, CxR, CxN, and patch constraints for any of the 12 forward–inverse set-ups (Table 3). However, qualitatively, CDR maps were heavily influenced by subspace choice for both groups. CxN maps were fragmented and difficult to interpret (Figs. 4 and 8). Indeed, the 3D maps, while anatomically ‘ignorant’ of the cortical subspace, often generated more robust CDR maps, particularly in the BFEC group (Fig. 2). There was a crucial interplay between the inverse model and the subspace constraint. Under the CxR constraint, both MNLS and LORETA generated implausible solutions. In BFEC, both algorithms consistently localized IED sources to multifocal, bilateral, and subtentorial regions of the cortical subspace (Fig. 3). In MTLE, IEDs were often sourced to the frontal lobes (Fig. 7). Under 3D, CxR, and CxN, L1 returned pin-point solutions. However, the cdnet maxima and near-maxima had a ‘skip’ distribution which was hard to reconcile, especially in view of the time-instant nature of CDR ESL. Of the four subspaces, the patch constraint exercised the most dominant effect, quantitatively (by virtue of the number of CDR parameters influenced directly by it), and qualitatively (by virtue of the way in which the CDR maps for MNLS, LORETA, sLORETA, and to a lesser extent, L1, resembled each other when applied to it) (Fig. 5, Fig. 9). The improved spatial resolution invoked by the patch constraint in both patient groups must be eyed with at least some caution. Under patch, MNLS and LORETA localized BFEC IEDs to the temporal lobe, while L1 localization was bilateral and multilobar (frontal, temporal, occipital). The use of an averaged MRI

1737

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Fig. 9. Series of CDR maps taken from MTLE patient 4 (Supplementary Figure S2 IED, averaged) for subspace patch. CDR maps refer to the same time-interval (within four milliseconds of the spike peak at 0.78 ms). Abbreviations: IED, interictal epileptiform discharge; MTLE, mesial temporal lobe epilepsy; CDR, current density reconstruction; patch, cortex with fixed extended sources.

brain, with its relatively smooth cortical surface, almost certainly influenced ESL accuracy in BFEC, a condition characterized electrographically by a dipolar source in sulcal space (viz. Rolandic fissure). This issue was less problematic in the MTLE group. This is probably because contemporary ESL can be expected, at best, to characterize TLE source behaviour at a sublobar level of accuracy, a scale of resolution that is less likely to be compromised by MRI averaging of the cortex. Of the inverse–subspace combinations, the use of sLORETA constrained by either the CxR or patch subspace appeared to generate the most robust, clinically meaningful results in both patient groups. Significant sLORETA orientation and location parameters by subspace were easiest to reconcile with the qualitative data (or CDR maps) and solutions were clinically plausible. It is quite likely that sLORETA ESL accuracy was adversely influenced by the smoothed cortical surface for the patch constraint. A telling example is the involvement of skip regions in MTLE Patient 4 (Fig. 9) with upper-scale (yellow) CDR mapping taking in two islands of cortex – the basal/basolateral temporal surface and the orbitofrontal surface. In such cases, current density maxima (defin-

ing best-fit, see Fig. 1) were co-incident with the island of activity seated at the ipsilateral temporal lobe. Co-activation of the orbitofrontal region is a likely by-product of the fact that its smoothed cortical surface sits normal to the en-face scalp electrodes. While ‘thresholding’ (as often applied in such settings) further ‘localizes’ the CDR map to the zone immediately surrounding the temporal lobe cdmax, we did not wish to introduce this bias to the qualitative data (we did not ‘threshold’ the quantitative data). While speculative, it is likely that this source of error would be offset (or at least minimized) with the use of a more dedicated temporal electrode array along with the application of individualized forward models that more accurately reflect cortical gyration. Despite the technical limitations associated with the use of a smoothed cortical surface, the ‘contraction’ of the CxR distributed ESL result by the patch set-up (lobar to sublobar in MTLE, proximate Rolandic to low Rolandic/Sylvian fissure in BFEC), did serve as a useful adjunct in teasing out possible IED source behaviour in these two focal epilepsies, particularly when it is recalled that sLORETA’s F-distribution (probability) scale did not differ between subspaces.

Table 4 CDR reliability. Intra (A1–A2) and inter (A1–B1, A2–B1) operator concordance kappa values for CDR outcome parameters. Ten spikes were re-modelled (five from each patient group). Fleiss’ scale was used to interpret kappa scores (>0.75 excellent, >0.4 fair, <0.4 poor). Abbreviations: lat, latency; cd, current density; x, y, and z, location axes; nx, ny, and nz, orientation vectors; fwhm, full width half maximum; snr, signal to noise ratio; rd, residual deviation; CDR, current density reconstruction.

A1–A2 A1–B1 A2–B1

lat

cd

x

y

z

nx

ny

nz

fwhm

snr

rd

0.99 0.97 0.98

0.91 0.96 0.99

0.91 0.99 0.90

0.73 0.77 0.39

0.97 0.92 0.90

0.99 0.78 0.77

0.83 0.63 0.86

0.91 0.69 0.61

0.99 0.95 0.96

0.96 0.89 0.82

0.91 0.89 0.83

1738

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

Table 5 CDR computation times. Means and ranges for computation times (s) for distributed modelling operations in BFEC and MTLE. No differences were noted between forward models. Times were most dependent on the nature of the inverse model. Abbreviations: L1, minimum L1 norm; LORETA, low resolution electromagnetic tomography; MNLS, minimum norm least squares; sLORETA, standardized LORETA; 3D, whole grid space; CxR, cortex with rotating sources; CxN, cortex with fixed sources; patch, cortex with fixed extended sources; BFEC, benign focal epilepsy of childhood with centrotemporal spikes; MTLE, mesial temporal lobe epilepsy. LI

LORETA

BFEC 3D

Mean Range

CxR

2.4 1.7–2.9 7.2 5–9.6

MNLS

sLORETA

MTLE 3.1 1.6–5.9

151.2 80–281

252.2 77–492

0.1 0.04–0.1

0.1 0.08–0.1

0.1 0.09–0.1

0.1 0.1–0.15

10 6–18.7

293.6 136–460

464.1 164–1045

0.2 0.07–0.3

0.3 0.2–0.35

0.3 0.2–0.3

0.3 0.3–0.45

CxN

1.2 0.8–1.6

1.6 1–2.6

35.1 16–54

49.6 19.1–118.

0.1 0.07–0.09

0.1 0.08–0.1

0.1 0.09–0.1

0.1 0.1–0.15

patch

1.4 1.1–1.9

1.8 1–2.9

37.1 16–72

58.2 19.4–124

0.1 0.07–0.3

0.2 0.08–0.4

0.1 0.09–0.3

0.3 0.1–0.5

4.4. Reliability MNLS and L1 inverse algorithms were common denominators when the least concordant comparisons arose. The most reproducible results consistently involved sLORETA. Results for LORETA fell between these extremes. On the other hand, it is somewhat surprising that as many as 85% of all intra- and inter-operator comparisons returned kappa values above 0.75 when it is recalled that the 11 CDR parameters examined were drawn from the final fit properties of only one ‘mini-dipole’ (corresponding to the CDR’s cdmax) of the several thousand mini-dipoles making up each distributed result. Moreover, the largest y location mismatches at 64.8 mm (MNLS inter-operator) and 63.9 mm (L1 intra-operator) should be interpreted in the context of the spatial resolution of scalp EEG that, at best, is in the order of 10 square centimetres. Similarly, the size of the vector differences for ny and nz should be judged knowing that (a) the scalp voltage topography changes on a scale of milliseconds across the course of a spike, and (b) such topography changes are predominantly driven by changes in the orientation of the source rather than by changes in its anatomical point location. 4.5. Software factors These results indicate that the rate determining step for distributed modelling CPT is the nature of the inverse model. The less constraining 3D and CxR subspace conditions would understandably potentiate the number of parameters (for location and orientation) that need to be ‘scanned’ by the inverse algorithm to finalize the optimal fit. MTLE CPTs were consistently longer than corresponding BFEC times. This was probably due to the noisier surface voltage distribution seen for MTLE spikes, a property that would expectedly demand more of the regularization strategy in limiting the best-fit CDR residual deviation. Hence, the use of a more dedicated inferior temporal electrode array in MTLE might serve to shorten CPTs by perhaps reducing the ‘blurring’ of the scalp EEG topography. Despite these caveats, the use of sLORETA as the inverse model would appear to override these factors, rendering such high-demand computations practicable for everyday clinical use. 4.6. Caveats The limitations of using forward models based on an averaged MRI brain scan have been discussed. However, as mentioned, this was the only practical way of directly comparing standard and high resolution boundary element methods with the finite element method. A second issue relates to the decision not to apply a min-

imum threshold setting for the CDR map display. Somewhat like functional MRI studies, arbitrary signal thresholds are commonly set in distributed ESL studies for the purposes of display and statistical analysis. It is likely that the setting of a minimum display threshold (at 50% fwhm for instance), would have ‘sharpened up’ the quality of the CDR images by limiting the contribution of ambiguous or outlying components to the final-fit. However, the fact that such threshold settings are genuinely arbitrary means that (a) selection bias is automatically introduced, (b) complex quantitative data are unlikely to be fully reconciled with corresponding, but preferentially trimmed, qualitative data, and (c) the premise of using single value decomposition and independent component analysis to identify signal that is to be explained over noise is, in effect, ‘down-weighted’ when, at the end of the ESL operation, up to half the data is put aside or ignored. All data, quantitative and qualitative, were carefully examined in tandem in the present study before any conclusions were drawn on the key questions raised, a strength rather than a limitation of the study. Finally it should be noted that the mean rd values across the 48 forward– inverse–subspace modelling set-ups approached 25%, well above the ideal 10% maximum rd targeted by most investigators in the field. This was probably due to a combination of factors – the use of single rather than averaged IEDs, the limited extent of the temporal electrode array in the MTLE patients, the use of a labelmatching system for electrode-anatomical landmark co-registration, and the use of the spike peak for ESL when, at this time-point, much of the scalp EEG signal component can be well removed from the original source due to the effects of cortico-cortical propagation. It was in this setting, and perhaps in spite of such caveats, that the clinical utility of ESL was evidenced. Indeed, it is tempting to speculate on the nature of the results were more covering scalp electrodes used (particularly over the infero-basal temporal regions in the MTLE group) with confirmed three dimensional electrode positions (using a spatial digitizer) and with individualized forward model reconstruction (from the patient’s own MRI brain scan). In our opinion, it is likely that the final results for the two patient groups would have still held with one main difference – the order of spatial resolution arrived at for each inverse algorithm. One would expect a rise in peak current density, a narrowing of the fwhm distributed volume, and a fall in residual deviation. The use of individualized cortical modelling should give rise to improved patch subspace solution accuracy due to the more reliable representation of patient specific, normalized cortex. Indeed, it should be emphasized that the distributed ESL results based on an individualized cortical surface should take precedence over the results obtained from generic (MNI) forward models due to the added risk imposed by the latter in mislocalizing sources when cortical constraints are enforced, particularly those that are normal to surface

C. Plummer et al. / Clinical Neurophysiology 121 (2010) 1726–1739

as with the patch constraint. Future studies are needed to characterize and compare ESL accuracy for distributed source modelling solutions based on individualized versus averaged realistic forward models. Despite the promising ESL results that can be achieved with an appropriate distributed modelling approach, a judgement on its clinical superiority over dipole modelling cannot be extrapolated from the present data. For that, a similar thorough-going comparison needs to be undertaken on the same raw data. Recent findings (Plummer et al., 2010, in press), based on the same BFEC and MTLE IED data studied herein, revealed no difference between the sLORETA–CxR/patch approach and a combined rotating/moving dipolar approach in explaining the surface EEG signal. There were two other key observations. The distributed best-fit consistently fell earlier than the corresponding dipole result during the time-course of the IED upslope – an earlier best-fit places the ESL solution closer to the original IED source. And second, ESL results were complementary, not conflicting – dipolar ESL localization re-enforced distributed ESL characterization at a sublobar level of resolution in both patient groups. To conclude, we have found that an optimal distributed modelling approach for both BFEC and MTLE incorporates sLORETA as the inverse model, CxR and patch as subspace constraints, and BEMi/s or FEMi as the forward model. With this sort of distributed modelling configuration, computation times are sub-second and CDR results are reproducible. While clinical validation of this distributed ESL approach is necessary with the use of simultaneous scalpintracranial recordings and individualized MRI datasets in surgical patient cohorts, our findings indicate that CDR is translatable to the routine clinical work-up of patients with unilateral BFEC of the typical Rolandic variety, and to unilateral MTLE secondary to hippocampal sclerosis. Acknowledgements Compumedics NeuroscanÒ supplied the Scan 4.3Ò and CURRY 5.0 software programs. The first author held an Australian National Health and Medical Research Council Postgraduate Medical Research Scholarship. Drs. Michael Wagner and Manfred Fuchs are employees of Compumedics NeuroscanÒ. Their contributions to the study were restricted to the provision of technical and mathematical advice. Technical assistance also provided by Ms. Linda Seiderer and Ms. Natasha Willems. Ò

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.clinph.2010.04.002. References Baillet S, Garnero L. A Bayesian approach to introducing anatomo-functional priors in the EEG/MEG inverse problem. IEEE Trans Biomed Eng 1997;44:374–85.

1739

Binnie CD, Dekker E, Smit A, Van der Linden G. Practical considerations in the positioning of EEG electrodes. Electroencephalogr Clin Neurophysiol 1982;53:453–8. Ebersole JS. EEG source modelling. The last word. J Clin Neurophysiol 1999;16:297–302. Ebersole JS. Noninvasive localization of epileptogenic foci by EEG source modelling. Epilepsia 2000;41(Suppl. 3):S24–33. Ebersole JS. Cortical generators and EEG voltage fields. In: Ebersole JS, Pedley TA, editors. Current practice of clinical electroencephalography. Philadelphia: Lippincott Williams and Wilkins; 2003. p. 12–31. Ebersole JS, Hawes-Ebersole S. Clinical application of dipole models in the localization of epileptiform activity. J Clin Neurophysiol 2007;24:120–9. Fleiss J. Statistical methods for rates and proportions. 2nd ed. New York: John Wiley and Sons; 1981. Fuchs M, Kastner J, Wagner M, Hawes S, Ebersole JS. A standardized boundary element method volume conductor model. J Clin Neurophysiol 2002;113:702–12. Fuchs M, Wagner M, Kastner J. Development of volume conductor and source models to localize epileptic foci. J Clin Neurophysiol 2007;24:101–19. Fuchs M, Wagner M, Kohler T, Wischmann HA. Linear and nonlinear current density reconstructions. J Clin Neurophysiol 1999;16:267–95. Grova C, Daunizeau J, Lina JM, Benar CG, Benali H, Gotman J. Evaluation of EEG localization methods using realistic simulations of interictal spikes. Neuroimage 2006;29:734–53. Haueisen J, Tuch DS, Ramon C, Schimpf PH, Wedeen VJ, George JS, et al. The influence of brain tissue anisotropy on human EEG and MEG. Neuroimage 2002;15:159–66. Helmholtz H. Über die Methoden, kleinste Zeitteile zu messen, und ihre Anwendung für physiologische Zwecke. Original work translated in Philos Mag 1853;6:313–25. Huiskamp GJ, van der Meij W, van Huffelen A, van Niewenhuizen O. High resolution spatio-temporal EEG–MEG analysis of rolandic spikes. J Clin Neurophysiol 2004;21:84–95. Jasper H. Report of the committee on methods of clinical examination in electroencephalography. Electroencephalogr Clin Neurophysiol 1958;10:370–5. Kim TS, Zhou YX, Kim S, Singh M. EEG distributed source imaging with a realistic finite-element head model. IEEE Trans Nucl Sci 2002;49:745–52. Nummemnaa A, Auranen T, Hamalainen MS, Jaaskelainen IP, Lampinen J, Sams M, et al. Hierarchical Bayesian estimates of distributed MEG sources: theoretical aspects and comparison of variational and MCMC methods. Neuroimage 2007;35:669–85. Pascual-Marqui RD. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find Exp Clin Pharmacol 2002;24(Suppl. D):5–12. Pascual-Marqui RD, Michel CM, Lehmann D. Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain. Int J Psychophysiol 1994;18:49–65. Pataraia E, Lindinger G, Deecke L, Mayer D, Baumgartner C. Combined MEG/EEG analysis of the interictal spike complex in mesial temporal lobe epilepsy. Neuroimage 2005;24:607–14. Plummer C, Harvey AS, Cook M. EEG source localization in focal epilepsy: where are we now? Epilepsia 2008;49:201–18. Plummer C, Litewka L, Farish S, Harvey AS, Cook MJ. Clinical utility of currentgeneration dipole modelling of scalp EEG. Clin Neurophysiol 2007;118:2344–61. Plummer C, Wagner M, Fuchs M, Harvey AS, Cook MJ. Dipole versus distributed EEG source localization for single versus averaged spikes in focal epilepsy. J Clin Neurophysiol 2010;27, in press. Scherg M, Bast T, Berg P. Multiple source analysis of interictal spikes: goals, requirements, and clinical value. J Clin Neurophysiol 1999;16:214–24. Wagner M, Kohler T, Fuchs M, Kastner J. An extended source model for current density reconstructions. In: Nenonen J, Ilmoniemi RJ, Katila T, editors. Biomag 2000: proceedings of the 12th international conference on biomagnetism. Finland: Helsinki University of Technology; 2001. p. 749–52. Wagner M, Wischmann HA, Fuchs M, Kohler T, Drenckhahn R. Current density reconstructions using the L1 norm. Adv Biomagn Res (BIOMAG) 1998;96:393–6. Wong PK. Potential fields, EEG maps, and cortical spike generators. Electroencephalogr Clin Neurophysiol 1998;106:138–41.