1. Introduction

The inverse problem in EEG/MEG is defined as the estimation of the distribution of electromotive force (EMF) in the brain from EEG/MEG by mathematical manipulations. As is well known, the solution to the inverse problem is not unique, since there exist silent EMF distributions that do not generate any electric and/or magnetic field at all outside the closed surface involving them [Rush, 1975]. This difficulty is usually circumvented by making use of some simplified models for the EMFs such as multipoles, moving dipoles, multiple fixed dipoles, continuous source distribution and so on: the parameters of these models can be determined uniquely by fitting the forward solution to the measured EEG/MEG.

Depending on the models for the EMF distribution, various methods have been proposed to solve the inverse problems: equivalent dipole method [Musha and Okamoto, 1999], BESA [Scherg and Picton, 1991], MUSIC [Mosher and Leahy, 1998], LORETA [Pascual-Marqui et al., 1994] and CDR [Wang et al., 1992] to name a few of the major ones. The equivalent dipole method is based on the moving dipole model. In this method EMF sources in the brain are approximated by a small number of current dipoles, and their locations and moments are estimated by fitting the EEG/MEG generated by them to the measured ones. In BESA (Brain Electric Source Analysis) and MUSIC (Multiple Signal Classification), the locations of dipoles are assumed to be fixed during some time interval, and they are determined from the potential distributions measured repeatedly during that time interval. So they are categorized as the spatio-temporal method. LORETA (Low Resolution Electromagnetic Tomography) and CDR (Current Density Reconstructions) are based on the continuous source distribution model. Continuous distribution is discretized as an array of numerous dipoles sited at known locations. Since the model allows a large degree of freedom, there are infinite numbers of solutions that can reconstruct the measured potential distribution on the scalp. In order to make the solution unique, an additional condition is imposed on the solution: namely some kind of norm defined on a huge vector composed of the vector components of all the dipoles must be minimized. So, these methods are categorized as the minimum norm method.

A great variety of methods including the major ones mentioned above had been installed already in some systems on the market. CURRY (http://www.neuro.com/) supports all major methods and provides many kinds of signal processing facilities. CURRY can be utilized for modalities such as MRI, fMRI, CT and PET as well as EEG/MEG. BESA (http://www.besa.de/) is a system designed to perform BESA method. BESA supports both EEG and MEG and can be used with MRI, fMRI and so on. BrainStorm (http://neuroimage.usc.edu/brainstorm/) is free software that realizes the equivalent dipole method, MUSIC and some kind of minimum norm method for both EEG and MEG. There is also free software called LORETA-Key that performs LORETA (http://www.unizh.ch/keyinst/). These softwares are free but they are not complete by themselves. For example, they usually lack ability to construct head models. On the other hand, CURRY and other systems on the market are, so to speak, all-round systems: they provide all the functions. However, a considerable amount of professional knowledge is required to master these systems. And they are too expensive for most researchers and/or clinicians who are interested in the inverse analyses.

Therefore, we are developing a complete but simple and user-friendly system named "Brain Space Navigator" or "BS-navi" for short. This Window-based system is specialized for the equivalent dipole method. BS-navi provides a utility to construct real shaped multi-shell head models from MRIs or X-ray CTs, and to estimate equivalent dipoles from measured EEGs efficiently based on the head model thus constructed. It also enables us to display estimated dipoles on an appropriate slice of MRI, or plot them in a 3-dimensional image of the head obtained by means of the volume rendering.

We focus on the equivalent dipole method because it can be easily and intuitively understood by clinicians. Moreover, since no restriction is imposed on the dipoles, this method seems adequate so long as electromotive forces are expected to be localized. On the other hand, in the spatio-temporal method such as BESA and MUSIC, dipole locations are supposed to be fixed during some specified time interval. Also in the minimum norm method such as LORETA or CDR, the additional condition implied to the solution in order to make it unique has no connection with the electrophysiological phenomena, and so it is quite difficult to evaluate the adequacy of this condition.

In section 2, the construction processes of a real shaped multi-shell head model are explained step by step. Then, an efficient forward algorithm specialized for multi-shell head model is introduced in section 3, and the actual EEG analyzing procedures using BS-navi are presented in section 4.

2.Construction of the Head Model

2.1. Boundary Extraction and Prototypical Head Model

The real shaped multi-shell head model is constructed based on the X-ray CTs or MRIs. As an example, an MRI is shown in Fig.1(a). The bone appears as dark region with small pixel values, while the brain and the skin appear as gray regions where the pixel values are a little bit larger. By selecting adequate threshold for the pixel value, gray scale MRI is converted into a binary image as shown in Fig.1(b). Then the boundaries between white and black regions are extracted as shown in Fig.1(c). These processes are repeated for all the slices. From these cross sections, the 3-dimensional shapes of the tissue boundaries are constructed (see Fig.2). They are the scalp surface, outer surface of the skull and the brain surface. Each boundary is represented by thousands of nodal points and contains thousands of triangles.

Figure 1. Original MRI (a), binary image (b) and extracted boundaries(c): the scalp surface, outer surface of the skull and the brain surface.

Figure 2. Prototypical head model (scalp, skull and brain) constructed from the boundary curves.

2.2. Reduction of Nodal Points

Each shell of the prototypical head model thus constructed is triangulated with many nodal points, and hence it takes long time to calculate potential distributions using this head model. We have to reduce the number of nodes in order to calculate potential distribution within practical time. To this end, we make use of a quasi-regular polyhedron that is, roughly speaking, a sphere triangulated homogeneously. An example with 912 nodes and 1820 triangles is shown in Fig.3(a). This quasi-regular polyhedron is placed in the scalp surface (b), and nodal points on the sphere are projected onto the scalp surface. Thus we get scalp surface (c) triangulated with the same number of nodes on the sphere, 912 in this case.

Figure 3. A quasi-regular polyhedron with 912 nodes and 1820 triangles (a). It is placed in the scalp surface (b), and nodal points on the sphere are projected onto the scalp surface (c)

Figure 4. Nodal points of the quasi-regular polyhedron are projected onto the skull surface, but the front part of the skull is deformed considerably.

However, this method does not always work well. For example, if the method is applied to the skull surface, the front part of the skull is deformed considerably as shown in Fig.4 because it is concave around there.

In order to overcome this difficulty, the skull surface, a concave polyhedron, is transformed into a convex polyhedron by inflation as shown in Fig.5. Inflation of a polyhedron is performed by supposing that some amount of electric charge is placed at each nodal point of the polyhedron, and all the sides of the polyhedron are made of elastic strings. The polyhedron is allowed to expand freely in a viscous liquid that damps down the inflation to a stationary state. Indeed, dynamic motion of the polyhedron is solved until it reaches stationary state. Then nodal points of a quasi-regular polyhedron are projected onto this inflated polyhedron. Finally, these projected points are brought back to the original polyhedron.

A quasi-regular polyhedron is specified by two integers a and b. The number of nodes on a quasi-regular polyhedron of type (a,b) is given by N=10(a²+ab+b²)+2 . And so, the number of nodes can be adjusted easily as occasion demands: a large number of nodes for accurate computation, or a small number of nodes for fast computation. Fig.6 shows the original skull surface triangulated by 2592 nodes, and the ones with reduced nodal points: 1962 and 912 nodes, respectively. As can be seen from these figures, number of nodes can be reduced without spoiling the shape of models.

Figure 5. Transformation of a concave polyhedron into a convex one by inflation

Figure 6. The original skull surface triangulated by 2592 nodes (a), and the ones with reduced nodal points: 1962 (b) and 912 nodes (c), respectively

2.3. Electrode Arrangement

Electrodes are arranged on the scalp surface following the standard 10-20 method, and their locations are determined based on 5 anatomical markers specified on X-ray CTs or MRIs. The markers are nasion, inion, both ears and the top of the head as indicated by yellow points in Fig.7(a). For example, the locations of the electrodes on the median plane are determined in the following way. The median plane is defined as the plane that contains nasion, inion and the top of the head. At first, the cross line of the scalp surface with the median plane is determined, and its length from nasion to inion is computed. Then, an electrode is placed at 10% of the total length from nasion. In the same way, the other electrodes are placed at 30%, 50%, 70% and 90% from nasion. All the 21 electrodes are arranged in the same way as shown in Fig.7(c).

Figure 7. Electrode locations are determined based on 5 anatomical markers specified on X-ray CTs or MRIs. They are nasion, inion, both ears and the top of the head (a). Based on these marker location, 21 electrodes are arranged following the standard 10-20 method (b),(c).

3.     An Efficient Forward Algorithm Specialized for Multi-Shell Head Models

Since our head model, real shaped multi-shell head model, is piecewise homogeneous, potential distribution can be calculated effectively by the boundary element method (BEM). In the direct BEM, each compartment is treated separately, and finally all the compartments are combined together through the boundary conditions in order to get the equations for the whole system.

Figure 8. Multi-shell head model composed of N compartments immersed in an infinite homogeneous region (a). After discretization, each compartment is represented by a block of resistance network with multiple input and output (b).

Suppose that a human head is approximated by a multi-shell model shown in Fig.8(a). The EMF source is contained only in the innermost compartment, namely the brain. After discretization, each compartment is represented by a block of resistance network with multiple input and output as shown in Fig.8(b) [Musha and Okamoto, 1999]. If the number of nodes on the boundary Γ₁ is N₁ , the innermost compartment Ω₁ , the brain, is represented by a block with N₁ output. And if the number of nodes on the boundary Γ₂ is N₂ , the second compartment Ω₂ is represented by a block with N₁ input and N₂ output, and so on. Vectors u₁ , j₁ refers to the potential distribution and the normal component of current density on the boundary Γ₁ , and u₂ , j₂ to those on the boundary Γ₂ . The input-output relation of the n-th compartment Ω_n is given in a matrix form as Eq.(1) for n=2~N .

(1)

where Tⁿ_n-1.n-1 , Tⁿ_n-1.n , Tⁿ_n.n-1 and Tⁿ_n.n are matrices obtained by BEM: Tⁿ_k.l is a N_k×N_l matrix. Replacement of n in Eq.(1) by n+1 yields an expression for Ω_n+1 , and substituting it in Eq.(1) eliminates u_n and j_n as

(2)

where

(3)

Note that Eq.(2) is nothing but the input-output relation for the compartment that is obtained by merging the two adjacent compartments Ω_n and Ω_n+1 . In this way, the adjacent compartments can be merged into one compartment. Repeating this process can combine Ω₂ through Ω_N together and all the parameters related to boundaries Γ₂
through Γ_N-1 are eliminated to get

(4)

For Ω₁ , which includes all the EMF sources, we have

(5)

where H₁¹ and G₁¹ are N₁×N₁ matrices obtained by BEM, σ₁ is the conductivity of the brain and u_∞ is the brain surface potential distribution that would be generated by the EMF source in the brain if the whole space were filled with the medium with the same conductivity as the brain. Elimination of j₁ from Eqs.(4) and (5) leads to

(6)

where H⁺ is the Moore-Penrose generalized inverse of the square matrix H . When no current flows into the head, we have j_N=0 , and when it does the values of j_N are known, while u_∞ is easily calculated from the EMF in the brain. From Eqs.(4) and (6), the potential distribution u_N on the scalp is determined as

(7)

Let be the M-dimensional vector composed of the potentials at M electrode sites on the scalp, which are determined from the N_N -dimensional vector by means of the linear interpolation. Then, is given as a product or a M×N_N matrix P and the vector u_N : = Pu_N. Since j_N is zero except for the current injection experiment, we get the final form as

(8)

Potential distribution on the scalp can be calculated just by operating the M×N_N transfer matrix T on the vector u_∞ . As mentioned before, u_∞ is the brain surface potential distribution that would be generated by the EMF source in the brain if the whole space were filled with the medium with the same conductivity as the brain. And so, it is quite easy to calculate u_∞ . Note that the transfer matrix T does not depend on the EMF source. It depends only on the head model and the electrode locations. So, we need to calculate the transfer matrix only once for each subject.

    4. BS-navi

BS-navi is specialized for the equivalent dipole method that can be easily and intuitively understood by clinicians. A dipole analyzing system can constructed just by installing BS-navi software into an existing Window-based EEG system. Therefore, it may be introduced not only by small research institutes but also by practicing doctors for clinical use.

As is explained in order in this section, BS-navi provides utilities to construct real shaped multi-shell head models from MRIs or X-ray CTs, to estimate equivalent dipoles from measured EEGs quite efficiently and to display estimated dipoles on an appropriate slice of MRI, or plot them in a 3-dimensional image of the head obtained by means of the volume rendering. Usually, a three-shell head model called the "SSB head model", which is composed of the scalp, skull and the brain compartments, is used. However, if necessary, the compartment of the cerebrospinal fluid can be taken into account by making use of the "SSLB head model". "L" in "SSLB" stands for the liquor or the cerebrospinal fluid. BS-navi handles single moving dipole (SMD) and two moving dipole (TMD) models. In principle, the equivalent dipole method is applicable regardless of the number of dipoles as long as the number of freedom of the measurement is larger than that of the dipoles. However, because of the ill-posedness of the inverse problem, it is practically impossible to estimate more than two moving dipoles accurately enough to the clinical applications [Okamoto et al., 1983].

4.1. Electrode Arrangement

Fig.9 shows the workspace of BS-navi. All the procedures such as the construction of head model, computation of transfer matrix, dipole estimation, and so on, are performed on this workspace. In Fig.9 a MRI is selected, and its histogram of the pixel value is shown. Based on this histogram, one can specify the threshold for the pixel value to make the binary image, or have the computer determine an adequate threshold automatically. Once the binary image is obtained, one can initiate edge extraction procedure by clicking the "Edge" button.

Figure 9. The workspace of BS-navi with the "MRI panel" and the "Edge extraction panel".

4.2. Head Model Construction and Computation of the Transfer Matrix

It is a little bit tedious task to extract boundaries for all the tens of slices, but once the boundary data are ready, it takes only a few seconds to construct a prototypical 3-shell head model. Head models with reduced number of nodal points are created automatically. Then, the head models thus constructed can be inspected visually by clicking "View" button after selecting one of "Scalp", "Skull" or "Brain" radio buttons. Fig.10(a) shows the interactive 3D display of the scalp model. One can rotate this model freely by dragging some point on it.

Figure 10. Interactive 3D display of the head model (a), and the panel for computing the transfer matrix (b).

In order to calculate the transfer matrix, one has to specify the number of nodal points in each boundary, the conductivity of each compartment, electrode arrangement and a file in which the resultant transfer matrix is stored as shown in Fig.10(b). Then one can initiate the computation by clicking "Matrix" button.

4.3. Estimation of the Equivalent Dipoles

EEG data are shown as in Fig.11(a). Two frames, frame no.250 and 500 in this example indicated by two magenta lines in this figure, are specified in order to eliminate the offset potential. The mean value of EEG during this time interval is subtracted from original EEG. Check boxes on the right side of the panel are used to indicate bad channels. EEG data specified as bad channels are not used in the dipole estimation. Topographies or the potential distributions on the scalp are displayed on lower left part of the panel. Positive potential regions are painted in red, and negative potential regions in blue.

In order to estimate the equivalent dipoles, you need a transfer matrix file, an initial location file, an EEG file and a file in which the results will be stored. One can specify all these files by clicking the buttons on the panel. Equivalent dipole method is applied to each frame of the potential distribution starting from the initial frame to the final frame indicated by two blue lines that are specified by clicking on the EEG panel. Computation starts when the "start" button is clicked. The progress of the computation is indicated as shown in Fig.11(b). Dipole locations are plotted as they are obtained in frontal, left sagittal and horizontal view together with the cross sections of the head model.

Figure 11. The panel for EEG/Topo (a), and the one for showing the progress of dipole estimation (b).

4.4. Estimation of the Equivalent Dipoles

By default, resultant equivalent dipoles are plotted together with the cross sections of the head model in frontal, left sagittal and horizontal view as shown in Fig.12(a). In this case, only three dipoles are plotted, which are estimated during the time interval indicated by two blue lines with the dipolarity or the goodness of fit index larger than the specified value of 98%. EEG curves, time course of the dipolarity and RMS (root mean square value of the potentials measured at 21 or 19 electrodes on the scalp) are plotted on lower right part of the panel. Topography at a specified time is also shown in the panel. Clicking this "CT/MRI" button, one can plot dipoles on MRI as shown in Fig.12(b).

Figure 12. Equivalent dipoles plotted with cross sections of the head model (a) and MRI (b).

If one clicks a point plotted on the sagittal plane for example, the frontal and horizontal planes are moved so that they include the specified point. And the new cross sections of the head with these new planes are shown as in Fig.13(a). Now, the specified point is on these frontal and horizontal planes, but it is not necessarily on the sagittal plane. Then, by clicking the corresponding point on the frontal or horizontal plane, one can move the sagittal plane so that the specified point is on all the three mutually orthogonal planes. In this way, it is possible to confirm the anatomical location of the dipole.

Figure 13. Confirmation of the anatomical location of a dipole. Clicking a point on the sagittal plane results in moving the frontal and horizontal planes so that they include the specified point (a). Then, clicking the corresponding point on the frontal or horizontal plane confirms that the specified point is on all the three mutually orthogonal planes.

It is also possible to plot the dipole on a truncated head as shown in Fig.14(a) by clicking the "VR" button, which stands for the volume rendering. Moreover, dipoles are plotted as dots in a transparent brain model as shown in Fig.14(b). This is an interactive 3D display, and one can rotate this model freely by dragging some point on it.

Figure 14. Dipoles plotted on a truncated head (a), or in a transparent brain model (b).

Figure 15. Dipole plotted with all the EEG curves (a). An example of the two-dipole approximation (b).

Many other functions are available. For example, it is possible to check all the EEG curves like Fig.15(a). Results from the two-dipole approximation are plotted as shown in Fig.15(b). One dipole of each pair is plotted in red and the other in green.

5. Conclusions

BS-navi is a software designed to perform the equivalent dipole method. In BS-navi, the electromotive forces are approximated by one or two dipoles. The locations, as well as their vector components at a point of time are determined by minimizing the square deviation between the potential distribution measured at that point of time, and the one generated by the dipoles. Since no restriction is imposed on the dipoles, this method can be applied so long as electromotive forces are expected to be localized, but it is practically impossible to estimate three or more dipoles.

BS-navi provides utilities to construct real shaped multi-shell head models (SSB and/or SSLB) from MRIs or X-ray CTs, to estimate equivalent dipoles (single and/or two-dipole approximation) from measured EEGs quite efficiently and to display estimated dipoles on an appropriate slice of MRI, or plot them on a truncated head obtained by means of the volume rendering.

Because of the efficient algorithm to calculate potential distribution in a real shaped multi-shell head model, an equivalent dipole can be estimated within 200ms by an ordinary personal computer in the case of single dipole approximation. This means that it is possible to trace a dipole in real time using BS-navi.

References

Rush S: On the independence of magnetic and electric body surface recordings, IEEE Trans., BME-22:157-167, 1975

Musha T and Okamoto Y: Forward and inverse problems of EEG dipole localization, Critical Reviews in Biomedical Eng., 27[3-5]189-239, 1999

Scherg M and Picton TW: Separation and identification of event-related potential components by brain electric source analysis, Event-Related Brain Research (EEG Suppl.42), 24-37, 1991

Mosher JC and Leahy RM: Recursive MUSIC: A framework for EEG and MEG source Localization, IEEE Trans., 45: 1342-1354, 1998

Pascual-Marqui RM, Michel C and Lehman D: Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain, Int. Jour. of Psychophysiol., 18:49-65, 1994

Wang J, Williamson SJ and Kaufman L: Magnetic source images determined by a lead-field analysis: the unique minimum-norm least-square estimation, IEEE Trans., BME-39, 665-675, 1992

Okamoto Y, Teramachi Y and Musha T: Limitation of the Inverse Problem in Body Surface Potential Mapping, IEEE Trans. BME-30:749-754, 1983