The electroencephalographic inverse problem is concerned with the estimation of neural sources underlying measured distributions of scalp potentials. One of widely used source models is the equivalent current dipoles (ECDs) model, where each ECD is completely characterized by three location parameters and three moment parameters specifying the dipole moment magnitude and orientation (Scherg & von Cramon 1985; Fender, 1987; He et al., 1987; Cuffin, 1995; Khosla et al., 1997). The inverse problem is then defined as the estimation of the location and moment parameters of ECDs whose modeled potentials best fit the measured scalp potentials in a least square sense (Fender, 1987; He et al., 1987; He & Mush, 1992; Cuffin, 1995). This least square fitting implies minimizing a cost function, where the cost function is the residue between the measured and modeled scalp potential data. This cost function that depends nonlinearly on the locations and linearly on the moment parameters, can be modeled using the measured potential distribution at a time instant (instantaneous-state ECD model) (Wood, 1982; He et al., 1987; Snyder, 1991; Khosla et al., 1997). Since there are not any constraints placed on the ECD parameters, instantaneous-state modeling (ISM) can allow both moving and rotating ECDs. In ISM, the measured scalp potential distribution at any given time instant is fitted by the source localization method individually. Thus, the positions and moments of the ECDs are allowed to vary as a function of time. This fitting can be done using either a single ECD or multiple ECDs (He & Musha, 1992; Khosla et al., 1997). Obviously, for these ECD-fitting methods this means it must be based on a priori knowledge on the number of sources, e.g., how many ECDs are active at a given time instant.

Recently, Kamijo et al. (2001) used a neural network method to identify the number of ECDs based on ISM. At its heart is a multilayer neural network which takes the measured scalp potential distribution as inputs, uses one hidden layer, and generates one output - the number of ECDs. This method can get a real-time solution of the number of ECDs for one or two ECDs. However, for the two-ECD model, the two dipoles were limited to be located at the different hemispheres, e.g., one in the left hemisphere and another in the right hemisphere. Since brain sources can often be located within the same hemisphere, the above constraint limits the applications of Kamijo et al’s method. Musha and Okamoto (1999) proposed an Akaike Information Criterion (AIC) method to estimate the number of ECDs. Their method showed promise in identifying the number of ECDs with ISM for a one- or two-dipole model. However, the influence of the penalty function (an important coefficient in the Information Criterion method; see Section 2 for details) for the identification accuracy was not considered in their study, and multiple ECDs identification for more than two ECDs was not considered (Musha & Okamoto, 1999).

In the present work, a method integrating the Information Criterion method with the Powell algorithm is proposed to identify the number of ECDs based on instantaneous-state modeling, with the purpose of avoiding the limitation of the active region of ECDs and using the instantaneous scalp EEGs. We first present the procedure to minimize the modified cost functions, which is obtained after dividing the ECD parameters into its linear and nonlinear components in the instantaneous ECD model. Then, we describe the basics of Information Criterion method and show how to implement it. Lastly, we evaluate the proposed method for some electrode systems using computer simulations. The effects of the number of electrodes, the penalty functions and the noise level on the estimation of dipole number are evaluated via computer simulations by the information criterion method.

(1)

where ||·||_F indicates the Frobenius norm, V is a m × 1 vector of measured potentials, G is the transfer matrix, and M is a 3K × 1 vector of ECD moments. Here, m is the number of electrodes and K the number of ECDs at the time instant of interest. The elements of the transfer matrix represent the potential at the pth electrode due to unit strength x , y and z moment components of the kth ECD ( k = 1,2, ··· , K ), respectively. The elements of the transfer matrix are complicated nonlinear functions of the electrode locations, the shape and electrical properties of the head, and the 3K × 1 dipole location vector r . These elements are calculated using either analytical expressions in spherical head models or numerical methods in realistic geometry head models derived from anatomical images obtained via MRIs (He et al., 1987; Cuffin 1995; Khosla et al., 1997; Musha & Okamoto, 1999; He & Lian, 2002).

In order to reduce the complexity of dipole source localization, we express the moments of dipoles by the so-called “normal equations,” M = G ⁺ V , where G ⁺ = (G^T G)⁻¹ G^T is known as the ”Moore-Penrose pseudo inverse” of the transfer matrix G . The cost function J can then be written as

(2)

The above least squares cost function J depends only on the transfer matrix G , a function of only 3K nonlinear location parameters (He et al., 1987; Khosla et al., 1997). Many methods have been proposed to minimize the cost function J , such as the Levenberg–Marquardt method, the downhill simplex method, the Powell method, the simulated annealing method and the genetic method (He et al., 1987; He & Musha, 1992; Khosla et al, 1997; Yamazaki et al, 2000; Kamijo et al, 2001). In the present study, the Powell algorithm is selected to minimize the cost function J because of its good performance in the small-residual case. In addition, the Powell algorithm uses the successive value to build up a numerical approximation to the derivative matrix (Powell, 1970) instead of using explicit expressions for the matrix. In our implementation of the algorithm we extended the one used by Press (1994) after adding some new features. For the dipole source localization, it is important to have a good guess of the initial value of dipole location parameters. In the present study, the initial guess method derived by Lypchuk T (1991) was used.

(3)

where log(•) is the natural logarithm and 6· n is the parameter number of assumed n ECDs. According to the rule of the IC method, the ECD model with minimum IC _n is selected as the optimal model, and the ECD number n in this model is the optimal ECD number.

The procedure of the proposed method involves two steps: (1) calculating the error functions ( E1 ,···, E_n) for the model set with the assumed ECD number. Each error function, can be obtained from the cost function (2) by the Powell algorithm, where the ECD number n of forward models is assumed to 1,2 , ···. (2) calculating the IC value in equation (3) with these error functions, and choosing the optimal ECD number n which minimizes the IC value.

In 1969, Akaike proposed the AIC, a statistic method incorporating Kullback-Leibler information with the use of maximum likelihood principles and negative entropy. Different penalty functions based on AIC have also been reported (Kwek, 2001). In the present study, we used the penalty functions C(m) as follows:

• C1 = 2 (Akaike, 1969);
• C2 = 2log(log(m)) (Hannan and Quinn, 1979);
• C3 = log(m) (Schwarz, 1978);
• C4 = 2log(m) (Broman, 1997);
• C5 = 3log(m) (Broman, 1997).

It can be seen from equation (3) that if more ECDs are assumed, the first term would become smaller since more ECDs would lead to a smaller fitting error, whereas the second term would become larger which is linearly proportional to the number of ECDs. Therefore, it is important to balance the goodness of fit and the number of parameters (Ljung, 1990). On the other hand, the penalty function will influence this balance. It can lead to different performance in the selection of parameter number. If a penalty function of larger value is used, the second term of Eq (3) would become larger, which would lead to a smaller estimate of the number of ECDs, as compared with that if a penalty function of smaller value would be used. The five penalty functions we evaluated in the present study are listed by their values from small to large. Therefore, the penalty function C₁ tends to give larger estimate of number of ECDs, while C₅ tends to give smaller estimate of number of ECDs, for a given set of instantaneous potentials over the scalp.

• The locations were generated randomly under the constraint that the distance between any two dipoles is larger than 2.0 cm in Case A and larger than 1.0 cm in Case B.
• The orientations were generated randomly.
• The moment magnitudes were less than 0.8.

Within each case, five groups of noisy data were obtained by adding normally distributed noise of different levels (5%, 10%, 20%, 30%, and 40%), where the noise level is defined as a percentage ratio of the root mean square value of the noise to that of the potential data at a given time instant. In the present study, the source localization method is used to calculate the potential error E. For evaluating the accuracy of our source localization method, the average location error and location error are defined as follows:

where is the location error for the i-th sample and and are the real and calculated dipole location vectors for the j-th dipole.

In the proposed method, the error functions En used in the IC method were obtained by performing source localization with the assumed ECD number n. Before discussing the identification of the ECD number, we firstly investigated the features of source localization in this section. Fig. 1 shows the average location errors ē_pos for all dipole cases, where the potential data (Case A) with 10% noise were generated in the 128-electrode configuration. Fig. 1 indicates that the average location error increases when the noise level or the dipole number increases. For multi-dipole cases, the noise has significant influence on the average location error.

Fig. 1 Average location errors for five cases with one to five dipoles. 128-electrode configuration was used and the minimum distance among dipoles was 2 cm (Case A). The average location error is the mean of location errors of 500 samples, e.g. , where the (more ...)

In order to illustrate the performance of our method, results from the potential data (Case A) with 10% noise in a 128 electrode configuration were examined. Fig. 2 shows the En and ICn of one sample of the three dipoles case. From Fig. 2(a), it is noted that the error functions decrease gradually while the assumed ECD number increases. And E1 – E2 is larger than E2 – E3 , E3 – E4 and E4 – E5 . Fig. 2(b) shows the identification results and effect of the penalty function in the IC method on the identification results. When the penalty functions C2 , C3 and C4 are employed, the minimum value of each of the three IC_n curves is 3, which shows that the correct ECD number can be obtained with the proposed method. When the largest C5 and the smallest C1 in all five penalty functions are used, the minimum IC value in these IC_n curves are 4 and 2, respectively. According to the rule of the IC method [See Method Section], the detected ECD number with five penalty functions ( C1,··· ,C5 ) are 4, 3, 3, 3 and 2 respectively (corresponding to the minimum in each of the IC_n curve). It indicates that the penalty function can affect the result of identification. Fig. 2(b) suggests that it is important to select the suitable penalty function. Fig. 3 show the distribution of the identification results in the IC method with different penalty functions. The results show that the accuracies of identification using the five penalty functions for one, two, three, four, and five dipoles. The IC method with the third penalty function ( C3 = log(m) ) provides the best performance in estimating the ECD number. The accuracies (in %) of from one to five dipoles cases were 97%, 94%, 88%, 49%, and 23%, respectively, where the accuracy defined as the percentage of the correct samples in all the samples.

Fig. 2 Error functions E and the IC values IC of one sample for the three dipoles case. The five penalty functions are used in this example. 10% noise and 128-electrode configuration, and Case A were used. In the upper panel, the error function E is obtained (more ...)

Fig. 3 Distribution of estimated ECDs numbers with five penalty functions for 10% noise, 128-electrode configuration, and Case A. (a) one dipole case, (b) two dipoles case, (c) three dipoles case, (d) four dipoles case and (e) five dipoles case.

Fig. 4 shows the effect of noise level on the performance of estimation of ECD number. Up to five dipoles were considered in the simulation. Fig. 4 suggests that the accuracies of C3 are still higher than other IC functions, and provides the best performance regardless of the noise level. This result suggests that, for a given electrode configuration and the head model, we can determine an unique optimal penalty function. For example, 128 electrodes and Case A, the optimal penalty function is C3 . Results also indicate that the noise can influence accuracies of identification of multiple ECDs cases significantly. For example, accuracies of three dipoles cases with C3 are 93%, 88%, 65%, 40% and 18% for five noise levels (from 5% to 40%). According to the IC method, if the sharp decrease in the error function is significant, the IC method can estimate the correct number of ECDs in most samples. It can therefore be inferred that the higher the noise level, the smoother the IC curve. Such smoothed IC curve reduces its performance in estimating the number of ECDs.

Fig. 4 Accuracies of five cases with the proposed method for 128-electrode configuration and Case A. The accuracy is a percent ratio between the number of correct samples and the number of all of samples, where the correct sample refers to the situation where (more ...)

Finally, the effect of different electrode configurations on the identification results for Case A and B were investigated. Table 1 and Table 2 summarize the highest accuracy by the IC method with the optimal penalty functions. For different electrode configurations, we suggest the following choices of penalty functions, e.g. 1) C4 for 32 and 64 electrode configurations, and C3 for 128 and 256 electrode configurations (Case A); 2) C4 for 32 electrode configuration, and C3 for 64, 128 and 256 electrode configurations (Case B). For two ECDs case, it is observed that the accuracies of Case A (Table 1) are close to Case B (Table 2). This result suggests that when the distance between any two dipoles is larger than 1 cm the identification accuracy does not change much with the dipole distances.

Table 1 Accuracies (%) of estimation by the IC method for different electrode configurations. Source configuration: Case A.

Table 2 Accuracies (%) of estimation by the IC method for different electrode configurations. Source configuration: Case B

In the present study, we have proposed a method to estimate the ECD number from instantaneous scalp EEG measurements by combining the Powell algorithm with the Information Criterion method. It overcomes the limitations in previously reported methods, such as the division of the dipole active region in the neural network method and the assumption of fixed dipole locations in the spatio-temporal model. Our main results of numerical simulations are as follows. Firstly, the choice of penalty functions in the IC method can affect the identification accuracy substantially. The optimal penalty functions in the IC method are determined, given the forward model and the electrode configuration. Secondly, noise is an important factor reducing the identification accuracy. In most cases, the sharp decrease existing in the case of low noise level disappears in the case of high noise level. The selection of the penalty function depends on the noise level. The optimal penalty functions for different noise levels are suggested according to the computer simulation results. Thirdly, the increased electrode number can increase the accuracy of identification markedly. The optimal penalty functions for different electrode configurations are suggested according to the computer simulation. Lastly, it is found that the accuracies of identification are similar when the minimum distance between dipoles is larger than 1 or 2 cm for 2 ECDs. When there are 3 ECDs, identification accuracies for Case A and B using the 256 electrode configuration were 82% and 65%, respectively. It suggests that the more the ECDs, the lower the spatial resolution of the present method. So for the purpose of estimating the number of ECDs, the minimum distance between dipoles would need to be increased, in order to obtain reliable estimation.

Dipole source localization has been pursued for decades (Wood, 1982; Scherg M, and von Cramon, 1985; He et al., 1987; Fender, 1987; He & Musha, 1992; Cuffin, 1995; Musha & Okamoto, 1999). One of the challenges of dipole source localization or other parametric source estimation methods lies in the need of knowing the number of equivalent current dipoles a priori. So when the number of source dipoles is appropriately estimated then the method works well in some cases. On the other hand, if the number of ECDs used in the dipole source localization is not correct, then misleading results maybe obtained. For this reason, distributed source imaging approaches that do not require a priori knowledge on the number of source dipoles have been investigated (Hamalainen & Ilmoniemi, 1984; Dale & Sereno, 1993; Gevins et al, 1994; Pascual-Marqui et al., 1994; Babiloni et al, 1997, 2001; He et al, 1999, 2001, 2002a, 2002b; Zhang et al, 2003). Therefore, development of methods being able to estimate reliably the number of source dipoles is of importance for further enhancing dipole source localization and ultimate improving our ability to image and localize brain electrical sources from noninvasive electromagnetic measurements.

While we have explicitly considered the equivalent moving dipole solutions, the proposed method should be of interest for sub-space source localization approaches (Mosher et al, 1992; Xu et al, 2004). The number of source dipoles is also a parameter which needs to be determined in sub-space source localization, although not as severe as the moving dipole localization problem.

Note that the present results are based on the use of simplified spherical head model. Further investigations using realistic geometry head models would provide further tests of the proposed method. Nonetheless, the present promising results suggest that the proposed method can play an important role in the identification of the ECD number for EEG dipole source localization. The principles of the proposed method should also be applicable to MEG dipole source localization.

We are grateful to anonymous reviewers for constructive comments, and Yuan Xu for proofreading the manuscript. This work was supported in part by NIH R01EB00178, NSF BES-0411898, and NSF BES-0411480.