Anatomical harmonics basis based brain source localization with application to epilepsy.

Brain Source Localization (BSL) using Electroencephalogram (EEG) has been a useful noninvasive modality for the diagnosis of epileptogenic zones, study of evoked related potentials, and brain disorders. The inverse solution of BSL is limited by high computational cost and localization error. The performance is additionally limited by head shape assumption and the corresponding harmonics basis function. In this work, an anatomical harmonics basis (Spherical Harmonics (SH), and more particularly Head Harmonics (H2)) based BSL is presented. The spatio-temporal four shell head model is formulated in SH and H2 domain. The anatomical harmonics domain formulation leads to dimensionality reduction and increased contribution of source eigenvalues, resulting in decreased computation and increased accuracy respectively. The performance of spatial subspace based Multiple Signal Classification (MUSIC) and Recursively Applied and Projected (RAP)-MUSIC method is compared with the proposed SH and H2 counterparts on simulated data. SH and H2 domain processing effectively resolves the problem of high computational cost without sacrificing the inverse source localization accuracy. The proposed H2 MUSIC was additionally validated for epileptogenic zone localization on clinical EEG data. The proposed framework offers an effective solution to clinicians in automated and time efficient seizure localization.

Introduction

The use of an accurate and efficient algorithm for Brain Source Localization (BSL) using Electroencephalography (EEG) measurements has been used in various neuroscience applications such as pre-surgical mapping in patients undergoing resection of epileptogenic zones[1-4], Brain Computer Interface (BCI) based prosthetic limbs[5,6], and attention deficit hyperactive disorder[7,8]. Among these diseases, epilepsy is the most important and common neurological disorder as of the world population is suffering from it. To solve this problem, EEG is regarded as the most commonly used non-invasive diagnosis tool due to its high temporal resolution, portability and cost-effectiveness. The process of brain source localization is carried out in two phases, which are forward modeling and inverse modeling.

Forward modeling involves estimation of the scalp potentials, given head model and current source. The volume conductor head model is chosen based on application, accuracy and computational complexity. The two kinds of volume conductor head models include numerical and analytical model. The numerical models are based on realistic modeling of human head with more complexity and low localization error. Boundary Element Method (BEM)[9], Finite Element Method (FEM)[10], and Finite Difference Method (FDM)[11] belong to numerical model category. In analytical modeling, head is approximated as a set of concentric spheres, each with homogeneous conductivity. This approximation results in low computational complexity. In literature, a single-shell head model[12] with a homogeneous conductivity was reported first. Later, three-shell[13] volume conductor head model was introduced where conductivity of skull was found to be lower than that of scalp/brain tissues. A more accurate four shell head model[14] with additional CerebroSpinal Fluid (CSF) layer between brain and skull, was introduced thereafter that forms the basis of the current work.

The recorded scalp potentials are further utilized to solve the inverse problem for active source localization. BSL methods are widely classified into distributed source (dipole-imaging) model and the dipole-fitting model. The distributed source model assumes that there is a large number of sources confined in an active region of the brain. The BSL, in this case, refers to the estimation of active source amplitudes and orientations using linear optimization techniques such as Bayesian methods, Minimum Norm Estimates (MNE), Weighted Minimum Norm Estimates (WMNE), Low Resolution Electrical Tomography (LORETA), and Local Auto Regressive Average (LAURA). The distributed source model based BSL being highly underdetermined, requires additional constraints like least energy, smooth variation and sparse focal nature of the sources are required to obtain better source estimation[15]. However, the localization of a limited number of equivalent dipoles is the most classical approach to solve the inverse problem. The a priori assumption in this solution is that only one or a few active areas in the brain generated the scalp potential field. Under this constraint, non-linear multidimensional optimization procedures allow to determine the dipole parameters that best explain the observed scalp potential measurements. This includes subspace based Multiple Signal Classification (MUSIC)[16], Recursively Applied and Projected (RAP) MUSIC[17], Truncated RAP-MUSIC (TRAP-MUSIC)[18]), beamforming, and genetic algorithms.

In the process of source localization, the inverse solutions is evaluated at each grid point i.e. nearly thousands of times to obtain convergent results. Practical computations of these inverse solutions are considerably expensive. The high computational cost involved thus becomes one of the major obstacles in solving the inverse problem of the EEG. In this work, we proposed a novel framework for active brain source localization by utilizing anatomical harmonics (Spherical Harmonics (SH) and Head Harmonics (H2)) domain processing on measured EEG data. SH have been widely used for acoustic source localization due to ease of array processing in SH domain[19,20]. It requires sensor array to be arranged on a rigid spherical baffle. In literature, the human head is approximated by spherical shape[12-14]. Hence, SH basis functions have been natural choice for accurate representation of EEG data over head[21,22]. Spherical harmonics were utilized in[21] for source reconstruction using dipole imaging method. Surface complexity of brain is estimated in[22] using SH to differentiate autistics, dyslexics, and controls subjects. The inherent assumption in application of SH basis function is that the data is present over entire sphere. However, the cap where the sensors are placed for data acquisition, assumes the shape between hemisphere and sphere. Hence, a more realistic basis function, Head Harmonics (H) were utilized in[23] for Infinite Homogeneous Isotropic Conductor (IHIC) head model and in[24] for Three Layer Concentric Spherical (TLCS) head model. In this work, more realistic four layer head model is first formulated in SH domain followed by H2 domain. The transformed model is then used for BSL using the dipole fitting approach. In particular, MUSIC and RAP-MUSIC BSL methods are presented in SH and H2 domain.

Four shell based EEG forward model in spatio-temporal domain

Spatio-temporal domain data model accounts for change in EEG potential w.r.t. space and time. Forward model is combination of head model and source model. The head model utilized herein is four layer concentric shell head model as shown in Fig. 1a, with the coordinate system in Fig. 1b. Brain, CSF, skull and scalp constitute the four concentric spherical layers with conductivity as , , and respectively. The brain comprises of homogeneous neural tissues and forms the innermost concentric sphere of radius br. This is surrounded by a concentric outer spherical shell of radius cr representing the CSF. This is followed by a third layer with radius dr for the skull. The outermost layer is the scalp of radius R. The BSL problem under consideration, utilizes Equivalent Current Dipole (ECD) source model where each dipole is parameterized by its location and dipole moment .

Figure 1

(a) Four shell head model, and (b) coordinate system.

I Sensors are placed over the scalp to record brain electrical activity due to P active dipole sources. The location of the ith sensor and pth source is given byrespectively, where is elevation angle measured downward from positive Z axis, and is azimuth angle measured anticlockwise from positive X axis. All the dipole sources are assumed to be inside the brain compartment of the four shell model. The EEG potential measured at ith sensor location generated by pth active source dipole at a time instant t is given by[14]and represents the vector dot product. Here denote the relative radii of brain, CSF and skull with respect to the scalp radius. Eccentricity of the dipole is and . Radial unit vector is given by and represents tangential unit vector defined in terms of vector cross product as is the legendre polynomial and is the associated legendre polynomial of order . The terms corresponding to have zero contribution in (2) and is therefore not considered. is the angle between ith sensor and pth dipole with given asUnder the assumption of fixed dipole source location and orientation, dipole moment vector can be written as [16], where is unit orientation vector and is magnitude of the dipole moment. Substituting in (2) and rearranging, the EEG potential at ith sensor due to pth dipole source can now be written asFrom (5), the total EEG potential at the ith electrode due to all the P active dipoles can be written aswhere is gain vector, expressed asfor four shell model. For I electrodes and P dipoles, the EEG potential in (6), can be written in matrix form asRe-writing (8) in compact form, we havewhere elements of is given by (7), and is referred as orientation moment matrix. In the presence of spatially and temporally zero mean white Gaussian noise with variance , the four shell spatio-temporal data model in (9) for time snapshots can be written aswhere manifold matrix is collection of manifold vectors given byIt is to be noted that matrix is order independent and hence would remain same in the order limited scenario.

Manifold matrix

In this Section, the spatial domain manifold matrix is expressed in terms of spherical harmonics basis functions () corresponding to sensor and source coordinates. Spherical harmonic addition theorem suggests[25]Substituting (14) in the gain vector in (7), we haveRadial and tangential orientation component of a dipole are expressed asrespectively. For a radially oriented dipole, and . For a tangential oriented dipole and . A mixed oriented dipole has both radial and tangential component as non-zero. As the EEG reflects the electrical activity of the post-synaptic potentials of pyramidal neuron cells oriented perpendicular to the cortical surface, dipoles with radial orientation is considered herein. For a radially oriented dipole, the element of a manifold vector can be written utilizing (15) and (16) asThe expression in (17), can not be computed in practice due to the infinite series summation. EEG mode strength determines the relative weighting of sensor and source spherical harmonics with the EEG signal order n. The parameter is a function of head model and the source radial position and is plotted in Fig. 2(a). The magnitude of diminishes significantly after order 40 and therefore the EEG mode strength becomes insignificant for higher order. This suggests that the summation in (17) can be truncated to a finite order upto without any significant error.

Figure 2

(a) EEG mode strength with order n of different source eccentricity. (b) The 2-D axial termed as ROI was taken at the depth of 3.2 cm from the top of the scalp surface. (c) Four active sources were placed pseudo-randomly on the ROI.

Under the assumption of finite order, n takes value from , where will be referred as EEG signal order. For a given order n, the degree m varies from . Thus, there are total distinct spherical harmonics . Hence, rewriting (17) in the matrix form,where is defined as is defined by replacing the argument to in (19). EEG mode strength is a diagonal matrix given byThe manifold matrix for a total of I electrodes and P dipoles can now be written asand , contain all sensors and dipoles position (,)’s respectively. It may be noted from (21) that the manifold matrix is the product of SH basis functions corresponding to sensor coordinates, EEG mode strength and the SH basis functions corresponding to source coordinates.

Four shell based EEG forward model in spherical harmonics domain

In this Section, four shell based EEG forward model is presented in spherical harmonics domain. The spherical harmonic basis functions is first introduced followed by a detailed derivation of EEG signal decomposition in SH domain.

The spherical harmonics basis functions

The head model considered herein is four shell spherical model. Hence, the SH basis functions becomes a natural choice to a function defined on the head surface. In this work, SH basis functions are explored for forward problem formulation and subsequent localization, followed by a more accurate head harmonics[23,24] basis functions based approach.

As EEG signals are real and discrete in nature, the real Spherical Fourier Transform (SFT) is applied instead of complex SFT. The real SFT pair for the discrete time EEG signal is given by[25]Here, dependency of V on (R, t) is omitted for notation simplicity. The sensor array order is dependent on spatial sampling scheme. Assuming nearly uniform sampling scheme with zero or negligible aliasing error, the number of sensors I should be at least , where assumes value in [1, 1.5][26,27]. Equivalently, the order up to which the EEG potential can be recorded without aliasing is limited by number of sensors asIt is to note that (26) presents an upper limit on array order that is function of number of sensors and the value of . The mathematical formulation of is beyond the scope of the present work.

The real valued spherical harmonics basis functions of order n and degree m is defined aswhere, the normalization constant is given by

Spherical harmonics decomposition of EEG signal

In this Section, spatio-temporal forward model in (12) is reformulated in computationally efficient spherical harmonics domain. In practice, the continuous potential on the scalp surface is spatially sampled by placing the EEG electrodes as per the international electrode placement system. Hence, rewriting the forward SFT in (25) by replacing the integral, we havewhere the quadrature weights are chosen such that the approximation error involved in (29) is minimized. Substituting for from inverse SFT (25) in (29), we haveTo ensure error-free sampling, following orthonormality condition must be satisfied.Here denotes the Kronecker delta function. (31) can be rewritten in matrix form asand . The quadrature weight matrix is estimated in the weighted least squares sense, given byA common choice for the weighting matrix is an identity matrix. An approximated spherical fourier transform in (30) can now be written in matrix from asand the corresponding inverse SFT is given in matrix form asThe quadrature weight matrix transforms the spatial domain forward data model to spherical harmonics domain data model. Multiplying from left to (12) results in SH domain data model aswhere the SH domain manifold matrix is written using (21) as . It is to be noted that the array manifold matrix of dimension I in spatio-temporal domain is reduced to dimension , where . This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain.

Four shell based EEG forward model in head harmonics domain

In this Section, four shell based EEG forward model is presented in head harmonics domain. The head harmonic basis functions is first introduced followed by a detailed derivation of EEG signal decomposition in H domain.

The head harmonics basis function

The concept of spherical harmonics has been extended to develop a set of HemiSpherical Harmonics (HSH) basis functions that are defined over the unit hemisphere[28]. HSH basis functions have been found useful for representation of surface reflectance functions[29], brain source localization[30] and surface description[31]. The HSH basis functions considers the EEG potential field to be defined on the hemisphere and therefore the human head is modeled as hemisphere[30]. However, the EEG signal acquired over head assumes the shape between hemisphere and sphere. Hence, based on realistic head dimension and EEG sensor placement system, another set of basis functions (called Head Harmonics (H)) was proposed in[23,24]. It is to note that for EEG acquisition, the electrodes are placed over the head surface with an elevation angle in the range range[27]. Therefore, it is more appropriate to utilize the H basis functions (where ) instead of SH basis functions (where ). The real valued orthonormal H basis functions () are expressed as[23]where is a normalization constant and is given byand is the shifted ALPs.

Head harmonics decomposition of EEG signal

The spatio-temporal forward model in (12), is reformulated in head harmonics domain. The SH decomposition of EEG signal in (38), assumes the potential field to be sampled on the sphere. However, the anatomical brain structure and the EEG sensor array placement suggest for the potential belonging to two third (i.e ) of the sphere only. Hence, the data model in (38) is modified by a conversion matrix () that transforms the SH domain data model to H domain. This is accomplished by changing the basis function aswhere the understanding of forward SFT may be used to write the element of the basis conversion matrix asRewriting (41) in matrix form we have,Multiplying to (43) and then utilizing the spherical array orthogonality condition in (32), the basis conversion matrix can be written asThe SH domain data model in (38) can now be transferred into the H domain data model with multiplication of the conversion matrix asIt may be noted that the transformation matrix transforms the spatial domain data model to the head harmonics counterpart. Additionally, the SH transformation matrix for converting the spatial domain data model to the spherical harmonics counterpart is . In case of no transformation from the original spatial domain, the spatial transformation matrix . The transformation matrices and , as detailed in Section ‘Head Harmonics Decomposition of EEG Signal’, converts the spatio-temporal data model to SH and H domain data model, respectively. In this process, the dimensionality of the data model is changed from I to , where [26,27]. This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain.

Spatial, SH and H2 based inverse brain source localization algorithms

Following the development of four shell forward data model in domain, dipole fitting algorithms are presented herein for EEG inverse problem. In particular, subspace based MUSIC[16] and RAP-MUSIC[17] methods are formulated in the domain. The subspace based methods utilize the estimated noise or signal subspace to compute the corresponding cost functions. A peak in the cost function is attributed to the neural activity. It is to note that spatial MUSIC and RAP-MUSIC algorithms utilize the spatial domain data model in (12) while the proposed SH and H based methods utilize the data model in (38) and (48) respectively. The spatial domain EEG signal is first transformed to the desired domain counterpart by multiplying the corresponding transformation matrix to obtain .

In the MUSIC algorithm, the estimates of active brain source location are obtained by projecting the domain array manifold vector over the estimated noise subspace. The domain MUSIC spectrum is expressed aswhere is the noise subspace projection matrix obtained from eigen value decomposition of covariance matrix [16]. The MUSIC spectrum in (49) can also be written in a simplified manner as[23]where is the minimum eigenvalue of (.), and is the Left Singular Vector (LSV) of . One of the drawback of MUSIC algorithm is that it searches for multiple local peaks in the entire head volume, which makes the task time consuming and susceptible to false source localization. To overcome the limitation, RAP-MUSIC was proposed that extracts the global maxima at each recursion steps. For the kth recursion step, the domain RAP-MUSIC spectrum is defined aswhere is the LSV of and . The orthogonal projector projects out the topography of the already found sources aswith . Matrix concatenates the array manifold vectors of already found sources, and is given by

Simulated data analysis

Various experiments were conducted to illustrate the advantages of SH and H domain based BSL algorithms. Both simulated and real EEG data were utilized for this purpose. In simulation, four layer concentric shell head model with radius of brain, CSF, skull, and scalp was considered to be 8.0 cm, 8.2 cm, 8.7 cm, and 9.2 cm respectively. The brain and scalp conductivities were set to 0.33 . The conductivities of skull and CSFs were taken to be 1/80 and 5 times of that of the brain. A total of Sensors were placed over head scalp utilizing electrode placement system. The order was set to 60 considering the diminishing nature of EEG mode strength for higher order. The number of discrete time samples was 200. The inter-grid gap was chosen to be 1 mm. The true source grid (i.e., the Region Of Interest (ROI)) was a 2-D axial slice at the depth of 3.2 cm from the top of the scalp surface as shown in Fig. 2b. The simulated sources were placed pseudo-randomly on ROI so that they were at least 3 cm from each other and at least 2 cm away from the center of ROI. Example of such a source distribution for is shown in Fig. 2c. Active source location and orientation is assumed to be fixed. Potential data were generated considering the radial orientation of active brain sources. All experiments were conducted with Monte-Carlo repetitions. It may be noted that the peak of the spectrum was explored only at ROI to reduce the computation load of the search over the entire volume.

Spatial resolution

The capability of BSL algorithms in estimating the spatially closed active source location is analysed herein by plotting the respective cost function spectrum. Two dipole sources were placed at and on ROI. The inter source correlation and distance between the two sources was taken to be 0.2 and 2 cm respectively. The Signal to Noise Ratio (SNR) was set to 5dB. The order was taken to be 3 for both SH and H domain. The 3-D view of spatial MUSIC, SH-MUSIC and H-MUSIC cost function is plotted in Fig. 3 (a)-(c). It may be noted that the spatial MUSIC fails to localize closely spaced active dipoles at low SNR, providing a single false location. However, the proposed SH and H counterparts localizes the two sources well with high resolution. This may be attributed to the increased contribution of the initial two eigenvalues of the covariance matrix corresponding to the two sources. The contribution of the initial two eigen values is in spatial domain. This contribution increases to and in the SH and H domain respectively.

Figure 3

3D view of (a) spatial MUSIC, (b) SH MUSIC and (c) H MUSIC for two active dipole source.

Additionally, the minimum distance between the two sources that have been localized via different BSL algorithms are quantified. Two sources were placed pseudo-randomly on ROI so that they were at least 2 cm away from origin. The inter source distance was varied from 0.5 cm to 2.0 cm. The minimum inter source distance is the shortest possible distance between two sources that can be localized. The minimum inter source distance and the corresponding localization error is illustrated in Fig. 4. It is to note that H based MUSIC and RAP-MUSIC attains least localization error and minimum inter source distance when compared to spatial and SH counterparts. In particular, the H RAP-MUSIC has the best performance.

Figure 4

Spatial resolution.

Localization error

In this Section, the Localization Error (LE) of spatial, SH and H based BSL methods is presented. The LE for L iteration and P dipole source is defined aswhere (x, y) is the true source location and is the estimated source location. Simulation experiments were conducted to compute the LE with varying SNRs, array orders, number of active dipoles, and inter source correlation.

LE with varying SNR and array order for single source

The number of active dipole is taken to be to study the effect of varying SNR and array order on LE. The source was placed in the pseudo-random manner. The LE is presented in Table 1 for MUSIC method in spatial, SH and H domain. As expected, LE decreases consistently with increasing SNR. LE variation is additionally presented with varying array order in accordance with (26). High localization error of SH and H MUSIC methods was observed for array order . This can be attributed to poor signal representation at lower order. The error further decreases with increase in order to give an optimum array order for a given number of sources. It is interesting to note that for single source, array order has the least localization error as highlighted in the Table 1 for SH/H based localization. An increase in the localization error was observed thereafter thus giving an upper bound on the array order. This is consistent with the array signal processing[20]. In general, it may be observed that the H MUSIC method outperforms the SH counterpart. For multiple source localization, the optimal array order with the least localization error is presented next.

Table 1

LE of (a) MUSIC, (b) SH MUSIC, and (c) H MUSIC with different SNR values and array order .

N_a	SNR
	-5 dB			0 dB			5 dB			10 dB			15 dB
	(a)	(b)	(c)	(a)	(b)	(c)	(a)	(b)	(c)	(a)	(b)	(c)	(a)	(b)	(c)
1	11.702	51.891	44.653	5.123	36.770	37.933	1.606	44.382	44.405	0.557	37.388	36.441	0.097	34.258	37.006
2		7.141	6.379		3.311	2.936		1.028	0.738		0.177	0.147		0.044	0.014
3		8.175	7.995		3.052	3.231		0.930	0.885		0.128	0.123		0.024	0.024
4		9.841	9.722		3.705	3.614		1.073	1.133		0.244	0.244		0.048	0.038
5		10.224	10.350		4.168	4.093		1.274	1.280		0.373	0.367		0.062	0.80
6		11.109	10.823		4.426	4.394		1.424	1.300		0.415	0.379		0.082	0.062
7		11.528	11.435		4.644	4.631		1.448	1.447		0.514	0.433		0.072	0.082
8		11.690	11.569		4.833	4.825		1.535	1.469		0.503	0.513		0.095	0.082
9		12.033	11.708		5.018	4.950		1.592	1.550		0.513	0.523		0.111	0.115
10		12.246	12.158		5.161	5.022		1.613	1.574		0.537	0.619		0.127	0.121

LE with varying array order for two sources

Performance of spatial, SH and H based MUSIC only method was explored in previous Section. In this Section, performance of MUSIC and RAP-MUSIC is presented for two active dipole sources with inter-source correlation as 0.2, placed pseudo randomly at 5 dB SNR. Localization error of MUSIC and RAP-MUSIC methods in the three domain with varying order is presented in Fig. 5a. As the spatial methods are independent of the array order, the localization error remains constant with order . It may be noted that H, SH domain methods outperformed spatial counterparts. A valley in localization error plot may be observed at order 3 for the two active dipole sources. In all the three domains, the performance of RAP-MUSIC is better than MUSIC.

Figure 5

Performance of Spatial, SH and H based MUSIC, RAP-MUSIC for two concurrent active sources at 5 dB SNR, (a) Localization error with varying array order (b) Localization error with varying inter source correlation (c) Computation time with varying array order .

LE analysis for RAP-MUSIC with varying array order and number of active dipoles

Optimum array order is experimentally estimated in this Section with varying number of active dipole sources. With an increase in the number of active sources, finding the local maxima in the MUSIC spectrum is susceptible to false source localization. Hence, LE results are presented in Table 2 for RAP MUSIC method only. The number of sources were varied from 1 to 5 at 5 dB SNR. The inter-source correlation between each adjacent pair was set to [0.20, 0.17, 0.17, 0.20]. It may be noted that the spatial domain RAP-MUSIC method localizes sources upto 3. On the other hand, the proposed SH and H methods successively localize all the five sources. An upward shift in optimal array order with an increase in the number of active sources was observed. This may be because of that appears to be dependent on the number of active dipole sources. Variation in the optimum array order with varying number of sources is highlighted in Table 2.

Table 2

Localization error analysis of RAP-MUSIC in different domain with varying number of sources.

N_a	Domain	Number of Sources
		1	2	3	4	5
N.A.	Spatial	1.606	2.650	6.247	-	-
2	SH	1.028	2.931	6.096	12.892	24.390
	H^2	0.738	2.562	4.949	12.126	23.441
3	SH	0.930	1.860	2.426	3.421	4.363
	H^2	0.885	1.744	2.348	3.466	4.585
4	SH	1.073	2.089	2.443	3.119	3.805
	H^2	1.133	2.008	2.393	3.120	3.798
5	SH	1.274	2.188	2.604	3.154	3.846
	H^2	1.280	2.107	2.575	3.224	3.851
6	SH	1.424	2.253	2.721	3.450	3.947
	H^2	1.300	2.264	2.688	3.301	3.902
7	SH	1.448	2.367	2.873	3.624	4.195
	H^2	1.447	2.419	2.837	3.570	4.174

LE with varying inter source correlation

The subspace based source localization methods work on the prior assumption that the active sources are uncorrelated. When inter source correlation is significant (synchronous source case), BSL methods provide erroneous estimation of dipole sources. Therefore, to study the applicability of SH and H based methods in the realistic scenario of synchronous source activation is of great significance. In the present simulation, the pearson correlation coefficient was varied from 0 to 0.8 for the case of two simultaneously active dipole sources at 5 dB SNR. The harmonics order for the SH and H based methods was taken to be 3 as discussed in Fig. 5b. Effect of inter source correlation on the localization error is presented in Fig. 5b. It can be observed that the SH and H domain based methods perform relatively better than their spatial domain counterparts in synchronous source localization. In particular the H based method outperforms the rest.

Computation time

The computational efficiency of SH and H domain based BSL methods such as MUSIC and RAP-MUSIC is presented herein. The transformation matrices and , as detailed in Section ‘Head Harmonics Decomposition of EEG Signal’, converts the spatio-temporal data model to SH and H domain, respectively. In this process, the dimensionality of the data model is changed from I to , where [26,27]. This reduction in dimensionality of the data is responsible for reduced computational cost in the transformed domain. The advantage of data model formulation in SH and H domain in terms of dimensionality reduction is illustrated as a reduction in computation time. A 64-bit type system having Intel Core i7-6700 processor with 4 cores @ 3.40 GHz, and 16 GB RAM was utilized. Windows 10 Pro Version 1803 with MATLAB Version: R2019b was used for the evaluation. Computation time was measured using tic, toc MATLAB function at 5dB SNR. In the simulation, number of dipole sources were taken to be . The average computation time was evaluated with varying sensor array order , as presented in Fig. 5c. It was observed that the computation time for SH and H based methods is same. Present simulation demonstrate the effectiveness of SH and H domain methods when compared with spatial domain methods. A significant reduction in the computation time for harmonics order is observed. It may be noted from Sect. 7.2 that the optimum harmonic order for two sources is 3. The high computational cost involved in solving the inverse problem using conventional spatial domain methods is a major obstacle in real-time EEG source localization.

Success rate

At a particular SNR value, the success rate of an algorithm is defined as the ratio of the number of localizations achieved with a 0 localization error to the total number of attempts. We further extend our performance evaluation by examining the success rate in for spatial, SH and H domain MUSIC methods. The number of Monte-Carlo trials were taken to be L = 1000. The order of SH and H harmonics was set to 2. The success rate in for a single source is presented in Fig. 6. At low SNR, significant difference in success rate can be observed between the SH, H and spatial domain methods. In particular, the H MUSIC is seen as the robust and most accurate source localization algorithm.

Figure 6

Success rate of MUSIC in the spatial, SH and H domain.

Real data analysis

For epilepsy surgery, localization of the seizure focus is time consuming and requires evaluation of neurological tests by experienced clinicians. In this Section, additional validation of the proposed head harmonic-based BSL is presented for epileptic seizure location using real clinical EEG data.

Subjects

An electrophysiological and imaging data of five patients with drug resistant epilepsy of long duration were included for study. The mean age of the cohort was 24.4 ± 7.33 with standard deviation range as 15 to 34 year. The patients were admitted to the Deenanath Mangeshkar Hospital & Research Center at Pune, Maharashtra, India. All of these patients were on multiple Anti-Seizure Medications (ASM). All patients underwent thorough pre-surgical evaluation with long term Video EEG, 3 Tesla Magnetic Resonance Imaging (MRI) scan, Positron Emission Tomography (PET) scan and neuropsychological assessments. All the patients had lesion on the MRI along with corresponding hypo metabolism on the PET scan with good electrophysiological concordance. All of these patients were discussed in weekly multidisciplinary meetings and were advised surgical resection. Imaging data was blinded and researchers did not know the site or side of the lesion before EEG source analysis. Informed Consent is obtained from all subjects and the study is conducted in accordance with relevant guidelines and regulations approved by the institutional ethical committee of Deenanath Mangeshkar Hospital & Research Centre (Ref: DMHRC Code- IHR/2021/Apr/NK/403).

EEG recordings

EEG recordings were done using standard 10-20 international system (19 electrodes along with ground and reference electrodes). Additional two electrodes T1 and T2 were added to detect anterior temporal discharges along with Electromyography (EMG) and Electrocardiography (ECG) electrodes. All the electrodes were pasted using standard gel and impedance was kept below 5k. Sampling frequency of EEG was set 2 KHz. Signals were passed through band pass filter of 1 Hz -70 Hz, followed by Notch filter to remove ’50Hz’ power line interference during recording. An epoch of 20 seconds was selected for source analysis by expert lab technician from the long duration recorded EEG, where first five seconds were pre-ictal recordings and another 15 seconds were seizure onset and propagation across the channels. Data pre-processing was performed offline on the selected epoch of data using the EEGLAB toolbox. The EMG and ECG channels present in the dataset were not considered.

Source localization

Head harmonic based MUSIC algorithm is evaluated on the real epileptic seizure data as it provides all the peaks corresponding to active dipoles simultaneously. The data recorded from 21 channels is interpolated to 64 channel data using EEGLab toolbox ‘spherical’ interpolation method. The number of source was assumed to be two as eigen value decomposition of covariance matrix had 95% of its energy confined to a rank-2 subspace. Five epileptic patients data were analyzed. Epileptogenic zone was observed in the left temporal lobe for the patients P1, P2, P3 and P5. The patient P4 had epileptogenic zone in the right temporal lobe. Fig. 7 illustrates the H MUSIC cost function plot of patient P1 having epiletogenic zone on the left temporal lobe. A pronounced cortical activation is focused in the corresponding epileptic region. The results are validated by the neurosurgeon who performed standard temporal lobe resection surgeries on left side in 4 patients (P1, P2, P3, P5) and on right side in 1 patient (P4). Our result demonstrates the effectiveness of H domain algorithm in successfully localizing the true epileptic activation region on real EEG data.

Figure 7

H MUSIC cost function plot of a subject P1 having seizure on left temporal lobe.

Ethical approval

All methods were carried out in accordance with relevant guidelines and regulations approved by the institutional ethical committee of Deenanath Mangeshkar Hospital and Research Centre (Ref: DMHRC Code- IHR_2021_Apr_NK_403 , Dated: 2nd Apr 2021). A generic patient consent is taken but specific consent is waived off by ethics committee because this study did not alter the patient management and it’s only involves studying of data acquired for other reason.

Conclusion

In this work, we proposed Spherical Harmonics (SH) and Head Harmonics (H) domain processing framework for active dipole localization. The four shell head model in spatio-temporal domain is formulated in computationally efficient SH and H domain. A qualitative and quantitative performance comparison of Multiple Signal Classification (MUSIC) and Recursively Applied and Projected (RAP)-MUSIC methods in spatial, SH and H domain is presented on simulated data. SH and H domain processing effectively solves the problem of high computational cost without sacrificing the inverse source localization accuracy. The performance of H based RAP-MUSIC was best when compared with SH and spatial domain counterparts. The performance dependency of SH and H domain methods on array order was experimentally established. The proposed H MUSIC was additionally validated for epileptogenic zone localization using clinical EEG data. The proposed framework offers an effective solution for clinicians in automated and time efficient seizure localization.