International Journal of Bioelectromagnetism Vol. 8, No. 1, pp. 20-30, 2006.

www.ijbem.org

On the Accurate Anatomical Constraints in EEG/MEG Source Imaging

Chang-Hwan Im^a, Bin He^a

^aDepartment of Biomedical Engineering, University of Minnesota, USA.

Correspondence: University of Minnesota, 7-105 BSBE, 312 Church St., Minneapolis, MN 55455, USA

Abstract. The present study introduces source models for EEG/MEG cortical source imaging, which can consider anatomical constraints more accurately. Conventional models have used coarse cortical surface mesh or sampled some vertices from fine surface mesh, and thus they failed to utilize full anatomical information which we can get with sub-millimeter modeling accuracy. Conventional ones placed a single dipolar source on each cortical patch and estimated its intensity by means of various inverse algorithms; whereas the suggested model integrates whole cortical patch area to construct lead field matrix and estimates current density that is assumed to be constant in each patch. We applied the proposed and conventional models to realistic EEG data and compared the results quantitatively. The quantitative comparisons showed that the proposed model can provide more accurate spatial descriptions of the neuronal source distribution.

Keywords: EEG, MEG, Brain mapping, Inverse problem, Lead field construction, Cortical Imaging

1. Introduction

Macrocolumns of tens of thousands of synchronously activated pyramidal cortical neurons are widely believed to be the main EEG and MEG generators because of the coherent distribution of their large dendritic trunks which are locally oriented in parallel, perpendicularly to the cerebral cortical surface [Nunez and Silberstein, 2000; He 2004, 2005]. Nowadays, this physiological phenomenon has been successfully adopted and widely used as a basic anatomical constraint in EEG and MEG source imaging [Dale and Sereno, 1993; Buchner et al., 1997; Kincses et al., 1999; Fuchs et al., 1999; Dale et al., 2000; Babiloni et al., 2001, 2003; Baillet et al., 2001a]. The source imaging with the anatomical constraint, which has been often called cortically distributed source model [Baillet et al., 2001b], resulted in elimination of spurious sources [Baillet et al., 1998] as well as reduction of crosstalk distribution [Liu et al., 1998], compared to conventional voxel (volume pixel) based imaging techniques.

To impose the anatomical constraint, many dipolar sources should be placed on the cortical surface, usually on the interface between white and gray matter of the cerebral cortex extracted from structural MRI, which is relatively easier to be detected than the other borders. We can further constrain each of these dipolar sources to be normal to the surface. Then, the strengths and/or orientations of the dipolar sources are determined using linear (L₂ norm) or nonlinear (L_p norm) estimation methods [Dale and Sereno, 1993; Fuchs et al., 1999]. To determine proper locations and orientations of the scattered sources, the cortical surface is usually tessellated into a huge number of triangular elements, the number of which is often exceeding several hundreds of thousands. Developments of medical image processing techniques and high-resolution structural MRI enabled us to get high-resolution cortical surface with sub-millimeter modeling errors [Dale et al., 1999; Fischl et al., 1999; Fischl and Dale, 2000]. Unfortunately, however, it is computationally inefficient to use whole cortical surface vertices for the source reconstruction because of the increased underdetermined relationship between the limited numbers of measurements and the number of unknown variables to be reconstructed. To reduce the number of possible source locations, some people resampled the fine mesh to a small number of larger triangles. Then, a unitary equivalent current dipole was placed in each node of the triangulated surface, with an orientation parallel to the averaged normal vectors of the surrounding triangles

[Darvas et al., 2004; Babiloni et al., 2005]. However, this kind of resampling approaches not only requires one more complex image processing procedure, but also has large possibility to lose accurate sulci-gyri structural information of the cortical surface. The other approach is to use a vertex-sampling process which has been frequently referred to as a decimation process [David and Garnero, 2002; Dhond et al., 2003; Lin et al. 2004]. Small number of vertices was downsampled from the cortical surface as regularly as possible and used for source reconstruction purpose; whereas the original mesh information was used only for visualization purpose. This approach is very simple to be applied, but it has some potential problems. First, the decimated vertices may not properly represent orientations of neighboring cortical vertices, especially around highly folded regions. This may result in significant reconstruction errors because neighboring dipolar sources with more appropriate orientations can be overestimated instead of the correct one due to crosstalk effect. Second, if the distribution of the decimated sources is irregular, the current intensity estimate can be distorted because density of the dipolar sources is inversely proportional to the estimated source intensity. Recently, Lin et al. [2005] tackled this problem by incorporating patch areas in the forward model to yield estimates of the surface current density instead of dipole amplitudes at the current locations as well as adopting loose orientation constraint (LOC), which allows some variation of the current direction from the average normal. Their results showed that the use of current density and LOC can improve overall accuracy of the source estimates.

In the present study, we have proposed an alternative source model to easily solve the problems of wrong orientation and irregular patch areas. Our approach does not use a single dipolar source that represents its neighboring vertices but uses all cortical vertices included in many small patches to construct lead field matrix. Since the constructed lead field matrix contains orientation and area information of all cortical vertices, the problems of the conventional approaches can be nicely solved.

2. Methods

2.1. Conventional Source Model

As briefly described in the previous section, conventional source model has utilized reduced number of dipolar sources which were decimated from fine cortical surface structure. Figs. 1(a) and (b) show an example of the tessellated cortical surface and decimated source positions, in which 432,654 original vertices were reduced to 7,866 source positions. As seen in Fig. 1(c), a decimated vertex cannot properly represent orientations of its neighboring vertices, even when the dipole vector was determined by summing all neighboring vectors. We then evaluated the areas of patches and found that the size and shape of each patch are highly irregular¹ [Lin et al., 2005]. Thus, to use single dipolar sources to represent high resolution anatomical images may cause substantial errors in both forward calculation and inverse estimation of bioelectromagnetic sources.

The variable that has been used for the conventional source model is the moment intensity of each dipolar source when orientation constraint is imposed. Then, the relation between the dipole intensity and the measured data can be expressed according to the following system:

x = AQ + n

(1)

where x is a column vector gathering the measurements on N_M sensors at a given time instant; Q is a column vector made of the N corresponding dipole intensities; A is the N_M×N lead field matrix; n is a perturbation or noise vector. The lead field A_ji is defined as electromagnetic quantity of j^th sensor induced by i^th dipolar source with unit intensity. Among various forward calculation methods, in this study, boundary element method (BEM) considering realistic geometry head model was applied [He et al, 1987; Hämäläinen and Sarvas, 1989].

2.2. Proposed Source Model

Contrary to the conventional source model, the proposed approach uses current density instead of the dipole intensity as a variable. First, the cortical surface is divided into N small cortical patches. Then, each cortical patch is assumed to have a constant current density J_i. In this case, the lead field A_ji is defined as electromagnetic quantity of j^th sensor induced by ith cortical patch with unit current density. Suppose that the i^th cortical patch includes k vertices. Each vertex has its own unit normal vector n_m = (n_1m, n_2m, n_3m) and virtual area v_m, where m = 1,…,k. The virtual area was assigned to each vertex as a third of the area of all triangles meeting at a vertex [Chupin et al., 2002].

^[1] The largest patch area is about 2 times larger than the smallest one.

This assumption is valid because the total virtual area remains equal to the actual area of the full tessellation. Then, we can evaluate the lead field a_jm defined as relationship between m^th vertex and j^th sensor by placing a unit dipole vector n_m at the position of the mth vertex. Eventually, the lead field at j^th sensor by a unit current density in i^th cortical patch is evaluated by integrating a_mj over the cortical patch as follows:

Fig. 2 illustrates the conceptual comparison between the conventional and proposed source models. We can see from the figure that the new model can represent size, shape, and orientation of each cortical patch without any loss of anatomical information.

Some previous studies have used similar concept of the constant cortical patch [David and Gareno, 2002; Kincses et al., 1999], but they didn’t aim for compensating anatomical information lost due to downsampling of cortical vertices.

Figure 1. An example of tessellated cortical surface and decimated sources: (a) cortical surface segmented and tessellated from an MRI T1 images (MNI standard brain); (b) original and decimated cortical vertices. 432,654 original vertices were reduced to 7,866 source positions; (c) area-of-influence around a decimated source (red vector). The orientation was determined by the vector sum of all vertices inside the patch.

2.3. Forward Calculation and Inverse Estimation

In the present study, realistic geometry head model was considered to calculate electric potential induced by a point dipolar source [Hämäläinen and Sarvas, 1989; He et al., 1987]. A three-layer boundary element model which consists of scalp, outer skull and inner skull was adopted [Hämäläinen and Sarvas, 1989; He et al., 2002], which will be described again in the next section.

We used a linear estimation approach [Dale and Sereno 1993, Dale et al. 2000] to reconstruct cortically distributed brain sources. The expression for the inverse operator W is

Figure 2. Conceptual comparison between conventional and new source models: (a) conventional model uses a single dipolar source in a cortical patch to evaluate lead field matrix; (b) proposed model uses all vertices inside the cortical patch to evaluate the lead field matrix. Each vertex in the patch has its own unit normal vector n_m = (n_1m, n_2m, n_3m) and virtual area v_m, where m = 1,…,k.

where A is the lead field matrix, R is a source covariance matrix, and C is a noise covariance matrix. The source distribution can be estimated by multiplying the measured signal at a specific instant x by W. If we assume that both R and C are scalar multiples of identity matrix, this approach becomes identical to minimum norm estimation [Liu et al., 2002]. In this study, the source covariance matrix R was assumed to be a diagonal matrix, which means that we ignored relationships between neighboring sources. The lead field weightings [Gorodnistky et al., 1995] were imposed to each diagonal entry of R. In this study, pre-stimulus time window was used to calculate C. λ² is a regularization parameter and was determined systematically using the following equation [Lin et al., 2004]:

where trace(.) and SNR represent sum of diagonal terms and signal to noise ratio, respectively.

2.4. Computer Simulation

Neuroelectromagnetic inverse problems are hard to be verified by in-vivo experiments because exact source locations inside of the real human brain cannot be estimated a priori. For that reason, artificially constructed forward data are widely used to validate MEG and EEG inverse algorithms [Kincses et al., 1999; Sekihara et al., 2001]. Hence, we applied the proposed approach to artificially constructed EEG data.

We adopted realistic conditions to construct artificial EEG data. We assumed 128 electrodes that were attached on a subject’s scalp according to the extended 10-20 electrode system. To utilize anatomical information, interface between white and gray matter was extracted from MRI T1 images of an MNI standard brain (www.mrc-cbu.cam.ac.uk/Imaging/Common/mnispace.shtml#evans_proc) and tessellated into 865,712 triangular elements and 432,654 vertices. To extract and tessellate the cortical surface, we applied BrainSuite developed in the University of Southern California, CA, USA [Shattuck and Leahy 2002]. For the accurate forward calculation, full head structures were taken into account and BEM was applied [Liu et al. 2002]. In the present study, the three-layer model, consisting of inner and outer skull boundary and scalp surface, was used. 5,372 boundary elements and 2,748 surface nodes were generated from the same MRI data. The relative conductivity values of brain, skull, and scalp were assumed to be 1, 1/16, and 1, respectively [Haueisen et al., 1997; Oostendorp et al., 2000]. Fig. 3 shows the boundary element model used in the present simulation.

Nowadays, for the forward simulations, generating artificial activation patches on a brain cortical surface has been popularized instead of activating some point sources [Im et al. 2003]. To generate activation patches and construct forward data set, the concept of virtual area was adopted. The activation patch was generated using the following process: 1) A point is selected as a seed of an activation patch area. 2) The patch area is extended by including neighboring vertices around the patch. 3) If the total virtual area of the cortical patch exceeds an aimed surface area, the extension of the activation patch is terminated.

In the present study, we generated one activation patch for each simulation. The patch was made of a set of dipoles with constant current density and orientations perpendicular to the cortical surface. Then, the current dipole moment at each vertex was calculated by the product of the current density and the virtual area. The temporal variation of current density J was assumed as follows:

After calculating electric potential at the 128-channel electrodes assuming 200 Hz-sampling rate, we added real brain noise, which was obtained from a pre-stimulus period of a practical EEG experiment. The original signal without noise was scaled in order for the signal-to-noise ratio to be approximately 10 dB and 7 dB. Fig. 4 shows an example of the artificial EEG signals with respect to time.

Figure 3. Boundary element model for EEG forward calculations. 5372 elements and 2748 nodes were generated. Note that the cortical surface meshes were not included in the EEG forward calculation. They were used only for positioning dipolar sources.

Figure 4. An example of simulated EEG signals with real brain noise (SNR = 7 dB).

Although we used constant current patches to generate artificial EEG signals, they were not corresponding to the cortical patches used for the inverse estimation because they were generated independently using full tessellation of the cortical surface. In other words, the ideal patches of constant activity do not recover exactly the patches considered in the source modeling.

3. Simulations Results

We applied three different source models to the artificial EEG data. The three cases are as follows:

(Case 1) One dipolar source was placed on each patch, which is the conventional source model. The dipole intensity was used as a variable and the dipole orientation was determined by summing up normal vectors of neighboring vertices inside the patch.

(Case 2) One dipole source was placed on each patch as in the conventional source model, but area information of each patch was considered. Current density was used as a variable and area of each patch was multiplied by the conventional lead field. The dipole orientation was determined by summing up normal vectors of neighboring vertices inside the patch. This case was simulated to investigate the influence of different patch sizes on the solution accuracy.

(Case 3) Lead field matrix was evaluated by integrating all cortical vertices inside each patch, which is the present source model.

The three cases were tested for 50 activation patches of which the positions and sizes were randomly chosen. Two different SNRs, 10 dB and 7 dB, were simulated for each patch location. The same inverse method given in Eq. (3) was applied to the three cases. Fig. 5 and Fig. 6 show examples of results for two different patch locations when SNR was 7 dB. Each figure shows exact source location and source distributions reconstructed at 100 ms. The magnitude of the variables was normalized with respect to maximum value and sources that exceeded 0.1 were visualized. From the results, we can see the followings:

- When comparing the results of (Case 1) and (Case 2), we can see that the irregular distribution of patch sizes has just small influence on the solution accuracy. We can also see intuitively that the emergence of small oscillatory sources was reduced slightly by using current density as a variable, which coincides well with a previous study [Lin et al., 2005].

- When the proposed source model was applied, the resultant distributions were more focalized compared to those of (Case 1) and (Case 2). Moreover, the shapes of the patches were more clearly reconstructed. Since the maximum magnitude was increased, many noisy sources had smaller normalized values than the cutoff magnitude (0.1) and removed from the visualization.

For more quantitative comparison, an assessment criterion named DF (degree of focalization) was introduced to measure the accuracy of the reconstructed source [Im et al., 2003]. The DF is defined as energy reconstructed in a patch area divided by overall energy^[2]. If the energy in the assumed patch is exactly the same as that of whole source space, DF becomes 1. On the contrary, if the source is perfectly mislocalized, DF will be close to 0. Thus, higher DF implies that the method can reconstruct more accurate and focalized source distribution. This assessment criterion is similar to the ROC curve analysis which investigates the relationship between true positive fraction and false positive fraction [Darvas et al., 2004]. Fig. 7 shows the comparison of the DF values averaged for 50 activation patch simulations with different patch sizes and locations. The DF values decreased as the SNR increased, but the overall tendencies remained unchanged. We applied t-test to verify statistical difference of DF values between new and conventional source models. The t-test applied from (Case 1) and (Case 2) against (Case 3) stated statistical differences (p < 0.0005), while the (Case 1) and (Case 2) does not differ in their mean DF values (p > 0.1). It can be seen from the results that the proposed source model could result in more focalized and accurate source distribution very consistently. These results demonstrate that the orientation errors affect solution accuracy much more than irregular patch size does.

^[2] In this study, the integration was performed only for limited sources of which the magnitude exceeded 30% of maximum value because the general minimum norm solution results in small DF values.

Figure 5. Results of a realistic EEG simulation – a patch was assumed around right inferior temporal lobe: (a) Exact patch location; (b) Source distribution of (Case 1); (c) Source distribution of (Case 2); (d) Source distribution of (Case 3). All quantities were normalized with respect to their own maximum value. Sources that exceed 0.1 are visualized. SNR = 7 dB. DF values of (b), (c), and (d) are 0.081, 0.083, and 0.184, respectively.

Figure 6. Results of a realistic EEG simulation – a patch was assumed around left frontal lobe: (a) Exact patch location; (b) Source distribution of (Case 1); (c) Source distribution of (Case 2); (d) Source distribution of (Case 3). All quantities were normalized with respect to their own maximum value. Sources that exceed 0.1 are visualized. SNR = 7 dB. DF values of (b), (c), and (d) are 0.065, 0.069, and 0.153, respectively.

Figure 7. Comparison of DF values averaged for 50 cortical patch simulations. Values in parentheses are standard deviations.

4. Discussion

Many researchers have been trying to make the spatial resolution of EEG and MEG comparable to that of fMRI or PET, which is believed to have higher spatial resolution than EEG or MEG. As for hardware aspect, more and more researchers are moving to a higher number of channels, and EEG acquisition systems with 128 channels are not uncommon anymore. Even 256 channel EEG systems are commercially available now [Suarez et al., 2000]. In the case of MEG, over 250 channel systems have been already popularized (e.g. CTF 275 MEG system) and being widely used. Combination of different modalities is now becoming a new alternative to enhance the limited spatial resolution of the electromagnetic-based imaging techniques [Dale et al., 2000]. As for software aspect, a lot of source imaging algorithms have been developed to get more accurate and focalized source estimate. Many methods are now adopting brain anatomy as a basic constraint because high-resolution structural MRI technology and image processing techniques enable us to get more accurate information on the human brain anatomy [Dale et al., 1999; Fischl and Dale, 2000].

The present study deals with a problem arising when imposing the anatomical constraint. Conventional studies have generated coarse cortical surface mesh or sampled some vertices from fine cortical surface mesh, to place the dipolar sources, but either model could not take full advantage of the accurate anatomical information which currently has sub-millimeter accuracy. The conventional model placed one dipolar source in each cortical patch and estimated the dipole intensity; whereas the proposed model utilizes every cortical vertex inside a cortical patch in evaluating a lead field matrix and estimates the current density of the patch. Since the proposed model can be regarded as a kind of numerical integrations over each source patch, it can consider more accurate shapes and orientations of the cortical patches.

If one tries to use all cortical vertices over several hundreds of thousands for the source reconstruction, two problems can arise. First, the use of more variables strengthens the underdetermined relationship between the number of measurements and that of unknowns, and thus it becomes more difficult to get a unique and accurate solution without additional a priori information. The second problem is the increment of computational cost. Supposing the numbers of decimated and original vertices are N_d and lN_d, respectively, the use of the whole source space takes more time in the order of O(l²) to calculate the matrix multiplications in Eq. (3). Besides, additional time required to construct lead field matrix is relatively negligible compared to that for the matrix multiplications, because the process requires relatively less increment of time in the order of O(l) when inversion of boundary element transfer matrix is stored in a computer memory. In the case of the proposed model, however, additional computation time is not significant because the number of variables does not change.

In the present study, the proposed model has been applied to realistic EEG simulations and the results have been compared to those of conventional ones. Boundary element method considering the realistic geometry head model was used for the forward calculations. For the verification, we generated 50 cortical patches with different sizes and locations, and simulated time-varying EEG signals. To be more realistic, real brain noise extracted from a pre-stimulus period of a practical experiment was added to the artificial EEG signals. For the quantitative comparisons, we adopted a criterion named DF, which can measure how well the reconstructed sources recover the original patches. We could see from the simulations that the consideration of area information did not improve the results as much as we expected. However, the application of the proposed model resulted in more accurate source estimate.

The proposed model is more plausible than the conventional ones because it considers the whole cortical surface area without any loss of anatomical information. The present study only applied a linear inverse operator to the inverse estimation, but the proposed model is applicable to other inverse techniques that use same anatomical constraint because this model does not deal with forward or inverse calculation methods but just construction of lead field matrix containing more accurate anatomical information. Although the simulations were performed only for EEG, the model can be applied to MEG source estimation in the same manner.

In summary, we have used a source patch model that can consider accurate anatomical information in EEG/MEG cortical source imaging process. The proposed model has been applied to realistic EEG simulations, and compared quantitatively with conventional ones. The present simulation results show that the proposed model provides enhanced performance in reconstructing cortical activations, as compared with conventional cortical source models. It is expected that the present model will serve as a useful means to get high resolution cortical source images in various EEG/MEG applications.

Acknowledgements

This work was supported in part by NIH EB00178, NSF BES-0411898 and the Biomedical Engineering Institute of the University of Minnesota . CH Im was supported in part by the Korea Research Foundation Grant funded by Korea Government (MOEHRD, Basic Research Promotion Fund) (M01-2005-000-10132-0).

References

Babiloni F, Carducci F, Cincotti F, Del Gratta C, Pizzella V, Romani GL, Rossini PM, Tecchio F, Babiloni C. Linear inverse source estimate of combined EEG and MEG data related to voluntary movements. Human Brain Mapping, 14: 197-209, 2001.

Babiloni F, Babiloni C, Carducci F, Romani GL, Rossini PM, Angelone LM, Cincotti F. Multimodal integration of high-resolution EEG and functional magnetic resonance imaging data: a simulation study. Neuroimage, 19(1): 1-15, 2003.

Babiloni F, Cinocotti F, Babiloni C, Carducci F, Mattia D, Astolfi L, Basilisco A, Rossini PM, Ding L, Ni Y, Cheng J, Christine K, Sweeney J, He B. Estimation of the cortical functional connectivity with the multimodal integration of high-resolution EEG and fMRI data by directed transfer function. Neuroimage, 24: 118-131, 2005.

Baillet S. Toward Functional Brain Imaging of Cortical Electrophysiology Markovian Models for Magneto and Electroencephalogram Source Estimation and Experimental Assessments. Ph. D. dissertation: University of Paris XI, 1998.

Baillet S, Riera JJ, Marin G, Mangin JF, Aubert J, Garnero L. Evaluation of inverse methods and head models for EEG source localization using a human skull phantom. Physics in Medicine and Biology, 46: 77-96, 2001a.

Baillet S, Mosher JC, Leahy RM. Electromagnetic Brain Mapping. IEEE Signal Processing Magazine, Nov: 14-30, 2001b.

Buchner H, Knoll G, Fuchs M, Rienäcker A, Beckmann R, Wagner M, Silny J, Jorg P. Inverse localization of electric dipole current sources in finite element models of the human head. Electroencephalography and Clinical Neurophysiology, 102: 267-278, 1997.

Chupin M, Baillet S, Okada C, Hasboun D, Garnero L. On the detection of hippocampus activity with MEG. In proceedings of the International Conference on Biomagnetism (BIOMAG2002), 2002.

Dale AM, Sereno MI. Improved localization of cortical activity by combining EEG and MEG with MRI surface reconstruction: a linear approach. Journal of Cognitive Neuroscience, 5: 162-176, 1993.

Dale AM, Fischl B, Sereno MI. Cortical Surface-Based Analysis I. Segmentation and Surface Reconstruction. Neuroimage, 9: 179-194, 1999.

Dale AM, Liu AK, Fischl BR, Buckner RL, Belliveau JW, Lewine JD, Halgren E. Dynamic Statistical Parametric Mapping: Combining fMRI and MEG for High-Resolution Imaging of Cortical Activity. Neuron, 26: 55-67, 2000.

Darvas F, Pantazis D, Kucukaltun-Yildirim E, Leahy RM. Mapping human brain function with MEG and EEG: methods and validation. Neuroimage, 23: S289-S299, 2004.

David O, Garnero L. Time-Coherent Expansion of MEG/EEG Cortical Sources. Neuroimage, 17: 1277-1289, 2002.

Dhond RP, Marinkovic K, Dale AM, Witzel T, Halgren E. Spatiotemporal maps of past-tense verb inflection. Neuroimage, 19: 91-100, 2003.

Fischl B, Sereno MI, Dale AM. Cortical Surface-Based Analysis II. Inflation, Flattening, and a Surface-Based Coordinate System. Neuroimage, 9: 195-207, 1999.

Fischl B, Dale AM. Measuring the thickness of the human cerebral cortex from magnetic resonance images. Proceedings of the National Academy of Sciences of the United States of America, 97: 11050-11055, 2000.

Fuchs M, Wagner M, Kohler T, Wischmann H-A. Linear and nonlinear current density reconstructions. Journal of Clinical Neurophysiology. 16: 267-295, 1997.

Gorodnistky IF, George JS, Rao BD. Neuromagnetic imaging with FOCUSS: a recursive weighted minimum norm algorithm. Electroencephalography and Clinical Neurophysiology, 95: 231-251, 1995.

Hämäläinen MS, Sarvas J. Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. IEEE Transactions on Biomedical Engineering, 36:165-171, 1989.

Haueisen J, Ramon C, Eiselt M, Brauer H, Nowak H. Influence of tissue resistivities on neuromagnetic fields and electric potentials studied with a finite element model of the head. IEEE Transactions on Biomedical Engineering, 44: 727-735, 1997.

He B, Musha T, Okamoto Y, Homma S, Nakajima Y, Sato T. Electric dipole tracing in the brain by means of the boundary element method and its accuracy. IEEE Transactions on Biomedical Engineering, 34: 406-414, 1987.

He B, Zhang X, Lian J, Sasaki H, Wu D, Towle VL. Boundary element method-based cortical potential imaging of somatosensory evoked potentials using subjects' magnetic resonance images. Neuroimage, 16: 564-576, 2002.

He B. Modeling and Imaging of Bioelectric Activity: Principles and Applications. Kluwer Academic/Plenum Publishers, 2004.

He B. Neural Engineering. Kluwer Academic/Plenum Publishers, 2005.

Im C-H, An K-O, Jung H-K, Kwon H, Lee Y-H. Assessment criteria for MEG/EEG cortical patch tests. Physics in Medicine and Biology, 48: 2561-2573, 2003.

Kincses WE, Braun C, Kaiser S, Elbert T. Modeling extended sources of event-related potentials using anatomical and physiological constraints. Human Brain Mapping, 8: 182-193, 1999.

Lin F-H, Witzel T, Hämäläinen MS, Dale AM, Belliveau JW, Stufflebeam SM. Spectral spatiotemporal imaging of cortical oscillations and interations in the human brain. Neuroimage, 23: 582-595, 2004.

Lin F-H, Belliveau JW, Dale AM, Hämäläinen MS. Distributed current estimates using cortical orientation constraints. Human Brain Mapping, in press, 2006.

Liu AK, Belliveau JW, Dale AM. Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proceedings of the National Academy of Sciences of the United States of America, 95: 8945-8950, 1998.

Liu AK, Dale AM, Belliveau JW. Monte Carlo simulation studies of EEG and MEG localization accuracy. Human Brain Mapping, 16: 47-62, 2002.

Nunez PL, Silberstein RB. On the relationship of synaptic activity to macroscopic measurements: Does co-registration of EEG with fMRI make sense? Brain Topography, 13:79-96, 2000.

Sekihara K, Nagarajan SS, Poeppel D, Marantz A, Miyashita Y. Reconstructing spatio-temporal activities of neural sources using an MEG vector beamformer technique. IEEE Transactions on Biomedical Engineering, 48: 760-771, 2001.

Shattuck DW, Leahy RM. BrainSuite: An automated cortical surface identification tool. Medical Image Analysis, 6: 129-142, 2002.

Suarez E, Viegas MD, Adjouadi M, Barreto A. Relating induced changes in EEG signals to orientation of visual stimuli using the ESI-256 machine. Biomedical Sciences Instrumentation, 36: 33-38, 2000.