IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997 1075 Spatio-Temporal EEG Source Localization Using Simulated Annealing Deepak Khosla,* Manbir Singh, and Manuel Don Abstract— The estimation of multiple dipole parameters in spatio-temporal source modeling (STSM) of electroencephalo- graphic (EEG) data is a difficult nonlinear optimization problem due to multiple local minima in the cost function. A straightfor- ward iterative optimization approach to such a problem is very susceptible to being trapped in a local minimum, thereby resulting in incorrect estimates of the dipole parameters. In this paper, we present and evaluate a more robust optimization approach based on the simulated annealing algorithm. The complexity of this approach for the STSM problem was reduced by separating the dipole parameters into linear (moment) and nonlinear (location) components. The effectiveness of the proposed method and its superiority over the traditional nonlinear simplex technique in escaping local minima were tested and demonstrated through computer simulations. The annealing algorithm and its imple- mentation for multidipole estimation are also discussed. We found the simulated annealing approach to be 7–31% more effective than the simplex method at converging to the true global min- imum for a number of different kinds of three-dipole problems simulated in this work. In addition, the computational cost of the proposed approach was only marginally higher than its simplex counterpart. The annealing method also yielded similar solutions irrespective of the initial guesses used. The proposed simulated annealing method is an attractive alternative to the simplex method that is currently more common in dipole estimation applications. Index Terms— Brain mapping, dipoles, electroencephalogra- phy, estimation, functional imaging, human, inverse problem, nonlinear optimization. I. INTRODUCTION T HE ELECTROENCEPHALOGRAPHIC (EEG) inverse problem is the estimation of the neural current sources underlying a measured distribution of scalp potentials. It is common practice to model neural activity by one or more equivalent current dipoles, where each dipole is completely characterized by three location parameters and three values specifying dipole moment magnitude and orientation [1]–[3]. Physiological justifications for the current dipole model are given elsewhere [4]. The inverse problem is then defined as the estimation of the location and moment parameters of one Manuscript received September 6, 1996; revised April 17, 1997. This work was supported in part by Grants NIA-NIH P50 AG05142, NIH NS-53213 by the University of Southern California, Los Angeles, CA, and an internal grant from the House Ear Institute, Los Angeles, CA. This work was initiated and carried out in part when D. Khosla was with the Department of Biomedical Engineering, University of Southern California, Los Angeles, CA. Asterisk indicates corresponding author. *D. Khosla is with the House Ear Institute, Los Angeles, CA 90057 USA (e-mail: dkholsa@hei.org). M. Singh is with the Departments of Radiology and Biomedical Engineer- ing, University of Southern California, Los Angeles, CA 90007 USA. M. Don is with the House Ear Institute, Los Angeles, CA 90057 USA. Publisher Item Identifier S 0018-9294(97)07415-6. or more dipoles whose modeled potentials best fit the actual measurements in a least-squares sense [2], [3], [5]. This least- squares fitting implies minimizing a cost function, where the cost function is the residual variance between the measured and modeled scalp signals. This function, which depends nonlinearly on the dipole locations and linearly on the moment parameters, can be modeled using either single time points of the measured potential distribution (instantaneous-state dipole model [3], [6]) or the entire EEG data set (spatio-temporal source model [1], [7]–[10]). In instantaneous-state modeling (ISM), single time points are fitted separately. Thus, the position and orientation of the dipole(s) are allowed to vary as a function of latency [2]. This fitting can be done using either a single dipole (when the activ- ity at that time point is restricted to a small portion of the cor- tex [9]) or multiple dipoles (when the scalp topography is due to two or more simultaneously active brain regions [11], [12]). It is expected that single-dipole sources can be localized fairly for a proper choice of the head model and good signal-to-noise ratio (SNR) of the measurements [13], [14]. However, when multiple-dipole are fitted to a single scalp topography, the lo- calizations can be unreliable because the number of parameters to be estimated starts approaching the number of independent recording channels [14]. An effective approach then is to fit multiple sources using multiple scalp topographies from a number of consecutive time points, i.e., STSM. This modeling approach was first applied to EEG analysis by Scherg et al. [1]. As the name suggests, STSM takes into account the spatio-temporal course of the signals as a whole, instead of considering each time point separately. This increases the range of measurements that can be used in model fitting. The dipole sources are assumed to be fixed in location throughout the measurement interval, while their orientations are either fixed or varying [1], [5], [7]. Other details and methodological considerations of STSM are discussed in [5] and [8]–[10]. In general, the performance of both of the above modeling approaches strongly depends on the optimization technique used to find the unknown dipole parameters. This occurs because the error or cost function is generally not globally con- cave and may contain multiple equivalent solutions [14]. The most widely used optimization methods for solving the EEG inverse problem can be classified into two groups: gradient methods, which use function and derivative information (e.g., Levenberg–Marquardt [15]), and search methods (nongradient techniques) which use only function values (e.g., Nelder–Mead downhill simplex [16]). Both of these methods minimize the cost function by iteratively adjusting the parameters of the dipoles. The final solution with these methods often depends 0018–9294/97$10.00  1997 IEEE