A new multimodal cortical source imaging algorithm for integrating simultaneously recorded EEG and MEG

Abstract

In this study, we propose a new cortical source imaging algorithm for integrating simultaneously recorded electroencephalography (EEG) and magnetoencephalography (MEG), which takes into account the different sensitivity characteristics of the two modalities with respect to cortical source orientations. It is generally accepted that MEG can record neuronal electrical activities with higher spatial resolution than EEG because magnetic field is less affected by the inhomogeneous conductivity profile of the human head than electric field. Nevertheless, it is also well known that MEG cannot reliably detect neuronal sources with radial orientation, whereas EEG is relatively less dependent on the source orientations than MEG. However, this intrinsic difference has not previously been taken into account in the integrative cortical source imaging using simultaneously recorded EEG and MEG data. To consider the different directional sensitivities of the two modalities, we conceptually decomposed the sensor and source spaces into radial and tangential components and developed an explicit formulation described as a general type of weighted minimum norm estimation. Using extensive computer simulations, we demonstrated that our proposed algorithm substantially enhances localization accuracy over almost the entire cortical surface compared to that of the conventional integration algorithm and results in a significant improvement especially when either a radial or tangential source component is dominant.

Keywords:

• cortical source imaging
• inverse problem
• multimodal neuroimaging
• electroencephalography
• magnetoencephalography

1. Introduction

Electroencephalography (EEG) and magnetoencephalography (MEG) have been widely used in clinical and cognitive neuroscience as powerful neuroimaging modalities that can estimate neuronal electrical activities with millisecond temporal resolutions. Despite their excellent temporal resolution, however, the spatial resolutions provided by EEG and MEG are not comparable to that provided by functional magnetic resonance imaging (fMRI) or positron emission tomography (PET) due to limited spatial samplings, uncertainties in the forward modelling and additive noise/artifacts 1–3. The spatial resolutions of EEG or MEG can be substantially improved by performing source imaging or by solving an inverse problem to estimate the EEG or MEG sources 4,5. In particular, EEG and MEG source imaging plays major roles in presurgical evaluation and surgical planning for patients with intractable drug-resistant epilepsy, because epileptogenic zones are only infrequently identified as lesions on structural MR images 6. In these applications, accurate estimation of neuronal electrical sources is of particular importance to reduce the size of the intracranial EEG grids and to avoid misplacement of the grid electrode locations.

Recently developed MEG instruments allow for simultaneous recording of magnetic and electrical fields originating from brain electrical activities, and simultaneous EEG and MEG data are routinely recorded in several clinical applications 7,8. Although both EEG and MEG signals are generated from identical neuronal current sources, their properties are intrinsically different. That is, MEG cannot reliably detect the radial component of neuronal currents, whereas EEG measures both radial and tangential source activities with slightly higher sensitivity to radial activities 9–15. Because of this intrinsic difference between EEG and MEG, some interictal spikes appear only in MEG but not in EEG, and sometimes vice versa 16–20. For example, in a study 20 with 43 patients suffering from temporal lobe epilepsy, eight patients had interictal spikes seen only in MEG recordings, while one patient had spikes that appeared only in EEG recordings. Moreover, some case studies have reported considerable spatial disagreements between the neuronal sources estimated using EEG and MEG 11,17,18,21,22. Therefore, some researchers have suggested that integrating simultaneously recorded EEG and MEG data may be more beneficial than using a single modality 23–30.

Several new approaches have been proposed to integrate simultaneously recorded EEG and MEG signals since Wood 31 first used single-channel MEG together with simultaneously recorded EEG to identify underlying neuronal sources in the somatosensory cortex. Using the equivalent current dipole (ECD) model, which assumes small numbers of current dipoles to approximate the flow of electrical current in a small brain area, Stok et al. 32 attempted to localize ECDs based on a dataset composed of 21-channel EEG signals and 21-channel MEG signals. To consider the different directional sensitivities of EEG and MEG, Cohen and Cuffin 10 proposed a three-step ECD localization algorithm, involving the following steps: (1) tangential components of each ECD are estimated using MEG data only; (2) electric potentials at each EEG electrode are calculated by solving an EEG forward problem using the estimated tangential ECD moments. Residual EEG data are then calculated by subtracting the calculated forward EEG data from the measured EEG data and (3) the radial component of each ECD is estimated using the residual EEG. Based on this idea, Goncalves et al. 33 and, more recently, Huang et al. 28 proposed alternative ECD localization algorithms that first estimate the locations and tangential components of ECDs with MEG and then determine the radial components of the ECDs and electrical conductivities of the brain, scalp, and skull using both EEG and MEG data sets. In these approaches, a concentric sphere head model was adopted based on the assumption that the MEG signal is generated only by the tangential components of ECDs. In a realistic geometry head model, however, this assumption could result in localization errors because the radial components of ECDs can also contribute to the generation of MEG 13. Moreover, when some brain activities are detected only using EEG, the stepwise approach described above cannot reliably localize the source locations. More basically, the ECD localization approach requires a priori assumptions about the number of ECDs and their initial locations, which are often difficult to determine due to a lack of preliminary information.

In contrast to the ECD model, the distributed source model assumes numerous current dipoles scattered in a limited source space, usually a tessellated cerebral cortex, and the dipole moments are determined using linear or nonlinear estimation methods 34. The distributed source approach does not require any preliminary information on the number or initial locations of brain activations, allowing inexperienced users to more easily localize MEG and EEG sources. Some studies have reported that the spatial localization accuracy of cortical source estimation can be enhanced in the distributed source model by using EEG and MEG data simultaneously 23,25,30,35,36. In these studies, both EEG and MEG recordings were incorporated into a single matrix equation expressing the linear relationship between the strength of each dipole moment and the simultaneously recorded EEG and MEG data sets. The use of more physical recordings is expected to enhance the overall localization accuracy compared to single-modality-based localization 30. However, our simulation results show that the conventional fusion imaging method does not always guarantee enhanced localization accuracy, particularly when a cortical source has a dominant radial component. Therefore, to successfully integrate EEG and MEG data, a new source imaging algorithm that can accurately estimate neuronal current distributions when neuronal sources are dominated by a radial component is required. To achieve the maximum synergy effect from the multimodal integration of EEG and MEG, we considered the different directional sensitivity characteristics of EEG and MEG. As mentioned above, these characteristics are fairly well understood and have already been considered in ECD localization algorithms 28,33; however, such information has not been effectively applied to EEG-MEG fusion imaging algorithms based on the distributed source model.

Our goal in the present study was to develop a new multimodal source imaging algorithm that can integrate simultaneously recorded EEG and MEG data to enhance overall localization accuracy, even under circumstances in which a specific directional component of the neuronal source is dominant. To account for the different directional sensitivities of EEG and MEG, we decomposed the sensor and source spaces into radial and tangential components conceptually and developed a new explicit formulation to solve the inverse problem. We then verified the performance of the proposed method through extensive computer simulations using a realistic head model.

2. Methods

2.1. Decomposition of source and sensor spaces

While trivial in a spherical head model, the radial directions of cortical sources in a realistic geometric head model need to be defined differently. In the present study, the radial direction r of a cortical source was defined as the orientation along which the total magnetic flux density generated by a unit dipole placed at a source location is minimized. To identify the radial direction, singular value decomposition was applied to the MEG leadfield matrix and the singular vector corresponding to the weakest singular value was assigned to the radial directional vector 28. Because the cortical current is generally assumed to be oriented perpendicularly to the cortical surface 35, the unit normal directional vector n can be explicitly defined at every location on the cortical surface. Using the Gram-Schmidt orthogonal process, the tangential direction t at the i-th cortical vertex (1 ≤ i ≤ m) can be uniquely determined as(1) where ⟨x, y⟩ denotes the inner product of the two vectors x and y, |x| denotes the Euclidian norm of the vector x, indices in parentheses represent the vertex number, and m is the number of cortical vertices. After evaluating the tangential, radial and normal directions at every cortical vertex, we define m by m diagonal matrices P_t and P_r, whose (i, i)-th elements represent the ratios of tangential and radial components to the normal component of the i-th cortical source, respectively, that is(2) and(3) for each 1 ≤ i, j ≤ m. Then, P_t and P_r satisfy(4) where I denotes an identity matrix of order m.

In the source space, the cortical source j =  [j₁, j₂, … , j_i, … , j_m]^T, oriented perpendicularly to the cortical surface, can be decomposed into tangential source components j_t =  [j_t,1, j_t,2, … , j_t,i, … , j_t,m]^T and radial source components j_r =  [j_r,1, j_r,2, … , j_r,i, … , j_r,m]^T using P_t and P_r, in the form of(5)

Using (4) and (5), the cortical sources can be rewritten as(6) and a schematic drawing of the source space is shown in Figure 1.

Figure 1. A schematic drawing presenting the decomposed radial and tangential components of a cortical source defined on i-th cortical vertex.

The leadfield matrices corresponding to EEG and MEG project the source space to the sensor space so that(7) where v represents the electric potential recorded at scalp electrodes, b represents the magnetic flux density recorded at SQUID sensors, and K_eeg and K_meg represent the leadfield matrices of EEG and MEG, respectively, each of which can be obtained by solving EEG and MEG forward problems.

In the sensor space, the electric potential v at the scalp EEG sensors is conceptually decomposed into v_t and v_r, which are generated by j_t and j_r, respectively, and the magnetic flux density b passing through MEG sensors is also decomposed into b_t and b_r, generated by j_t and j_r, respectively; that is(8)

Then, the relations between sources and the directional components of v and b can be described as(9) where the directional leadfield matrices, K_t,eeg, K_r,eeg, K_t,meg and K_r,meg have the following relationships:(10)

When linear inverse estimation is applied, the tangential component of cortical sources can be estimated using either(11) or(12) where G_t,eeg and G_t,meg represent the tangential inverse operators with respect to EEG and MEG, respectively. Similarly, the radial components of cortical sources can be estimated using either(13) or(14) where G_r,eeg and G_r,meg represent the radial inverse operators of EEG and MEG, respectively. The derivations of the directional inverse operators are provided in the Appendix.

2.2. A proposed algorithm

First, we consider the decomposition of the estimated source . Similarly to (5) and (6), the estimated source can also be decomposed into tangential and radial components; that is,(15)

The tangential component was estimated from (12) and the radial component was estimated from (13). Substituting (12) and (13) into (15) gives(16) and applying (5) and (9) to (16) gives(17)

When we define the error between the exact and reconstructed sources e as(18) we can rewrite (17) as(19)

Rearranging (19),(20) where

Now, we consider a general constrained minimization problem for weighted minimum norm estimation (WMNE) 37,(21) where W represents the source weighting matrix and(22)

Then, the solution of (21) is known to be(23)

In our problem, we need to minimize the error between the exact and reconstructed sources e, defined in (18); thus the definition of our problem becomes(24)

In our problem, y represents the combined EEG and MEG data K represents the combined leadfield matrix and j is the source strength that we need to estimate. In the linear equation y = Kj, each of the leadfield matrices of EEG and MEG as well as each of the EEG and MEG signal datasets are normalized so that the elements in the combined leadfield matrix equation have equivalent orders. To accomplish this, the leadfield matrices of EEG and MEG (K_eeg and K_meg) are first normalized by the matrix norms of EEG and MEG, respectively. Then, the same scale factors are applied to the EEG and MEG signal data (v and b), respectively. The normalized leadfield matrices and data sets are then stacked into a single leadfield matrix K and a data vector y.

If we set the source weighting matrix W in (21) to(25) we can transform the general WMNE problem in (21) into our problem in (24), because the equation (22) becomes(26)

Then, we can directly utilize the known solution of WMNE written in (23). Thus, the solution of our problem given in (24) can be found by substituting (25) into (23):(27) where(28)

When additive noise is present in the signals, a regularization term needs to be introduced. Then the expression for the solution 1 becomes(29) where λ is a regularization parameter and was determined using the generalized cross validation method 38. Note that the notation is a dummy variable used only to construct the source weighting matrix W during the formulation, and denoted in (29) is the final source estimate.

3. Simulation studies

3.1. Simulation set-ups

Neuroelectromagnetic inverse problems are hard to verify using in vivo experiments because exact source locations inside a human brain are not known a priori. Therefore, artificially constructed forward data have been widely used to validate MEG and EEG inverse algorithms 38,39. Hence, we applied the new inverse method introduced in the previous section to artificially constructed EEG and MEG data sets. The MEG sensor layout used for the simulation was adopted from a commercial 148-channel whole-head magnetometer system (Magnets 2500 WH; Biomagnetic Technologies, San Diego, CA; see Figure 1(A) in 40). EEG sensors were generated by projecting the 148 MEG sensors to the nearest points on the scalp surface to allow for a direct comparison of the performances of the two modalities 27.

Because synchronously activated pyramidal cortical neurons, located perpendicularly to the cortical surface, are believed to be the main EEG and MEG generators, many studies have adopted this physiological phenomenon as a basic anatomical constraint in EEG or MEG source imaging 36,37,41,42. To use this anatomical information, we extracted the interface between the white and gray matter from structural MRI images of a standard brain atlas (180*217*180 pixels, 1*1*1 mm) provided by the Montreal Neurological Institute. To extract and tessellate the cortical surface, we used CURRY6 for Windows® (Compumedics, Inc., El Paso, TX). Although advances in medical image processing and high resolution structural MRI allow high resolution cortical surfaces with sub-millimeter modelling errors to be obtained 43,44, it is computationally inefficient to use whole cortical surface vertices for source reconstruction purposes because of the underdetermined relationship between a limited number of sensors and a larger number of source locations. To reduce the number of possible source locations, a smaller number of vertices was downsampled from the cortical surface as regularly as possible and used only for source reconstruction purposes, whereas the original mesh information was used only for visualization purposes 45,46. In our simulation study, we downsampled 11,373 vertices and 22,774 triangular elements from the original dense cortical vertices.

For the accurate forward calculation, we applied a first-order node-based boundary element method (BEM) to calculate the forward magnetic field and electric potential distributions 47. We obtained EEG and MEG leadfield matrices by applying BEM to three-layer tessellated boundary surfaces, consisting of the inner and outer skull boundaries and scalp surface, which were generated from the same MRI data using CURRY6. A total of 3393 nodes were used for the node-based BEM computation. The relative conductivity values of the brain, skull, and scalp were assumed to be 1, 1/16 and 1 (S/m), respectively 48,49.

We performed extensive computer simulations to quantitatively compare the performance of our proposed method with that of the conventional method in which the source weighting matrix W in (21) is the identity matrix. We assumed that current sources were constant cortical patches composed of a set of dipoles with constant dipole moments and orientations perpendicular to the cortical surface. To generate activation patches and construct a forward data set, we adopted the concept of a virtual area. The activation patch was generated using the following process: (1) a point was selected as a seed of an activation patch; (2) the patch was then extended to include neighbouring vertices around the patch and (3) if the total virtual area of the cortical patch exceeded the targeted surface area, the extension of the activation patch was terminated 50. Because some cortical surface regions were too distant from sensors to generate detectable EEG and MEG signals in a noisy environment, a limited numbers of source patches was chosen, from which the distance to the scalp surface did not exceed 50 mm. Source patches on the cerebellum were also excluded. Consequently, our computer simulations used 4568 cortical patches with a mean area of 68.4 ± 17.6 mm².

Whereas, there are many evaluation metrics available to assess the source localization accuracy of ECD models (Euclidean distance, for example), quantifying the reconstructed cortical sources in the distributed source model is far from trivial, and there is no consensus as to which metric should be used 27,50,51. Because we used an extended patch source in our simulations, we assessed the accuracy of source estimation using the criterion called degrees of focalization (DF), which quantifies how much of the reconstructed source is contained in the reference source patch 50. This validation metric was defined as follows and ranged in value from 0 to 100:(30) where Ω denotes the whole source space and Π denotes the reference source patch.

To further investigate the influence of the orientations of the cortical sources on localization accuracy, the 4568 source patches were classified on the basis of the proportion of the radial component γ, defined as(31) where Π denotes the source patch area, N_v represents the number of vertices in Π, i represents the i-th cortical vertex, j_r,i is the amplitude of the radial component of a cortical source at the i-th vertex, and γ ranges from 0 to 1. A value of γ close to 0 indicates that a cortical source patch is oriented in the tangential direction, whereas a value of γ close to 1 indicates that a source patch is oriented in the radial direction. The histogram depicted in Figure 2 shows the distribution of the number of source patches with respect to γ with a bin size of 0.05. From this figure, it is clear that more cortical sources are oriented in a tangential direction than are oriented in a radial direction, reflecting the well-known fact that the area of the sulci is larger than that of the gyri. Because the number of source patches with a γ value exceeding 0.8 was not sufficient to estimate average localization accuracies, which was smaller than a third of the smallest number of source patches in a single bin in which the γ value was less than 0.8, we excluded those patches when analyzing the simulation results shown in Figure 3.

Figure 2. Distribution of cortical source patches with respect to the proportion of the radial component of sources. Because the number of patches in which the proportion was greater than 0.8 was not sufficient to estimate average localization accuracies, we excluded these patches when analysing the simulation results shown in Figure 3.

Figure 3. Average localization accuracies in four cases with respect to the proportion of the radial component: (case A) cortical sources estimated with EEG data alone; (case B) cortical sources estimated with MEG data alone; (case C) the conventional combined EEG-MEG source estimation method was applied and (case D) the proposed EEG-MEG integration method was applied.

3.2. Simulation results

We then evaluated source localization accuracies for four different cases: (1) source imaging with only EEG data (EEG alone case); (2) source imaging with only MEG data (MEG alone case); (3) source imaging using conventional integration method based on a simple order normalization to combine the EEG and MEG data sets (conventional method) and (4) source imaging using the proposed multimodal integration method (proposed method). We reconstructed cortical source distributions using the EEG/MEG forward data sets simulated for each of the 4568 cortical patches. White Gaussian noise was added to the simulated EEG and MEG signals. We set the signal-to-noise ratio (SNR) values of EEG and MEG to 10 and 30, respectively, considering that EEG data generally contains more noise than MEG data in practice when the SNR values were defined as ten times of the log-scaled square root of the ratio of the simulated signal power to the noise power. The variations in the DF values averaged in each bin with respect to the proportion of the radial component (γ) values are shown in Figure 3. The localization accuracy of the EEG alone case increased and that of the MEG alone case decreased as the γ value increased, demonstrating that EEG and MEG source localization results are dependent upon the source orientations when constant background noise is added to the neural electrical signals. These results indicate that MEG is better than EEG at estimating tangential sources whereas EEG is better than MEG at estimating radial sources.

The conventional method generally yielded more accurate source estimation results than either EEG alone or MEG alone cases, but it did not always enhance localization accuracy, particularly when the radial source component was dominant (γ > 0.6). In contrast, our proposed method enhanced the localization accuracy for every γ value and moreover, resulted in a significant improvement in localization accuracy, particularly when either the radial or tangential components of the cortical sources were dominant.

The spatial distributions of the DF values measured at each cortical source patch for the four cases listed above are shown in Figure 4; the DF value of the source patch was assigned to the centre of the patch. Figure 4(A) and (B) show that the EEG alone and MEG alone cases provided relatively more accurate source estimates only around some specific cortical areas. Interestingly, the cortical areas with high DF values in the EEG alone case (Figure 4A) and the MEG alone case (Figure 4B) were separated and complementary to each other. In contrast, the conventional method enhanced the localization accuracy in cortical areas in which either the EEG alone or MEG alone cases could not accurately reconstruct source distributions (Figure 4C). However, the enhancement in the absolute DF values was comparable to those obtained from either of the two modalities, showing that an additional synergy effect was not obtainable using the conventional approach. In contrast, our proposed method showed significantly enhanced localization accuracy compared to that of the conventional integration method (Figure 4D).

Figure 4. Localization accuracy of four cases mapped on the cortical surface: (A) EEG method; (B) MEG method; (C) conventional integration method and (D) proposed integration method. The colour code indicates the DF value assigned at the centre of each reference source patch. The colour map of the DF value was thresholded at 3%.

4. Conclusion

Although both EEG and MEG measure neuronal electrical activities, their sensitivities to neuronal source orientations are intrinsically different. The different directional sensitivity characteristics were considered in dipolar source localization for enhanced localization accuracy, but have not been adopted in the cortical source imaging. In this study, we proposed a new multimodal cortical source imaging method to integrate simultaneously recorded EEG and MEG recordings, considering the different directional sensitivities of EEG and MEG. Our simulation studies showed that the proposed method can enhance the source localization accuracy significantly, regardless of the cortical source locations.

The different sensitivities of EEG and MEG to source orientations have been investigated by several research groups 13,52–54. However, consideration of the directional sensitivity characteristics in the integrated EEG/MEG source imaging has not been reported prior to our study. Previous studies concerning integrated EEG and MEG source imaging have focused only on how to preprocess EEG and MEG leadfield matrices to reduce the condition number of the combined linear system 36,55. We developed an explicit formulation based on the weighted minimum norm estimation with a source weighting matrix reflecting the decomposed directional components of EEG and MEG, yielding results robust to the cortical source orientations.

As mentioned in the introduction, we expect that the most useful application of our method will be the localization of epileptogenic sources, because precise estimation of epileptogenic foci is of critical importance in surgical planning for epilepsy patients. We plan to apply our method to practical EEG and MEG data simultaneously recorded from successfully treated partial epilepsy patients in future studies.

Acknowledgements

This work was supported in part by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. 2010-0015604) and in part by the research fund of Hanyang University (HY-2011 -N).

Appendix

The directional inverse operator, G_t,eeg, can be computed as follows. Applying (5) and (11), the estimated source can be written as(32)

Using (9), (A.1) can be rewritten as(33) and from (5)(34)

When we define the error between the exact and reconstructed sources ϵ as(35) we can rewrite (A.3) as(36)

Rearranging (A.5), we have the following relationship:(37)

To find G_t,meg that minimizes the error ϵ for any j, we consider the following minimization problem:(38) where tr{·} is the trace of a matrix and is defined as the sum of the main diagonal entries. This expression can be explicitly minimized by taking the derivative with respect to G_t,meg using the following properties:(39) and the derivative set to zero. Then, we have(40)

The explicit formulations for G_t,eeg, G_r,eeg, and G_r,meg can be derived in a similar way. Then, we have(41)