Multiple Dipole Source Position and Orientation Estimation Using Non-Invasive EEG-like Signals

The problem of precisely estimating the position and orientation of multiple dipoles using synthetic EEG signals is considered in this paper. After determining a proper forward model, a nonlinear constrained optimization problem with regularization is solved, and the results are compared with a widely used research code, namely EEGLAB. A thorough sensitivity analysis of the estimation algorithm to the parameters (such as the number of samples and sensors) in the assumed signal measurement model is conducted. To confirm the efficacy of the proposed source identification algorithm on any category of data sets, three different kinds of data-synthetic model data, visually evoked clinical EEG data, and seizure clinical EEG data are used. Furthermore, the algorithm is tested on both the spherical head model and the realistic head model based on the MNI coordinates. The numerical results and comparisons with the EEGLAB show very good agreement, with little pre-processing required for the acquired data.


Introduction
The human brain composed of neurons that connect with each other via electrical signals. One can record and measure these activities using an electroencephalogram (EEG). An essential use of the EEG is in locating the generating source of these signals, usually approximated by dipoles. This is important because, in some particular circumstances, neurons may not function optimally and could make the equivalent dipole generate abnormal signals. This could be the result of seizures or other brain disorders. In order to isolate such disorders, the challenge is to find a non-invasive way to locate the anomalous source. In [1], the authors addressed the location of abnormality for mildly depressed patients. In this case, only few regions were associated with depression. Thus, to treat these disorders, source localization is crucial and vital in clinical subjects exhibiting such neural activity [2]. EEG signals' source localization has been extensively studied. Cohen et al. in [3] additionally compared the accuracy of using EEG signals versus magnetoencephalogram (MEG) signals and showed that EEG signals are as useful as MEG signals for source localization problems. Furthermore, in [4], the authors described trend source localization methods using the finite element method for modeling the human head, as well as defined a time-slices approach.
Among all the studies, there are two important steps for source localization: (1) The forward model for EEG signal approximation; and (2) The inverse problem for locating the generating source. In the forward problem, the electrode potentials are calculated based on the given source properties. Many review articles address different forward problem approaches pertaining to source localization as in [5]. Other studies such as [6,7] focused on a specific forward problem approach such as implementing the boundary element method (BEM) and its effect on source localization error. Moreover, in [8], the effect of forward

Head Model
In this study, two different head models are evaluated: (1) a hemi-spherical (half sphere) head model, as shown in Figure 1; and (2) a realistic head model based on the Montreal Neurological Institute (MNI) coordinates. A comparison of these two head models is presented in Figure 2. Note that in both head models, +y axis passes through the nasion, and the +x axis passes through the right ear. First, let us consider a half sphere as a brain model which is shown in Figure 1.  According to Figure 1a, the location of any ith dipole with respect to the origin by utilizing spherical coordinates can be presented in a vector l i = [l x i l y i l z i ] T , given as: In this study, the radius of the brain spherical model is considered as r = 10 cm, and thus, r i ∈ [0, 10). Moreover, the two angles are defined as θ i ∈ [0, π 2 ) and φ i ∈ [0, 2π) to cover the whole half-sphere head model.
Similarly, the jth sensor location is described as w j = [w x j w y j w z j ] T , obtained using Figure 1b illustrates the relation of the ith dipole and two selected sensors. The vector s i = s i µ i ∈ R 3×1 is the dipole signal strength (with units as Coulombmeter), µ i ∈ R 3×1 is the unit vector denoting the orientation of the dipole, and s i ∈ R is the magnitude. An alternate representation could be s i = s x i s y i s z i T . Note, the latter requires three parameters to specify the i-th dipole as opposed to four parameters for the former. However, for this study, the 4 parameter representation is used since it provided better estimates of the dipole orientation. Figure 1b also shows the distance between a dipole and two different sensors, which is denoted by d ji where j is the sensor number, and i the dipole number.

EEG Signal Measurement
The mathematical model for the measured EEG signal is adapted from [9,21], together with a measurement noise model (n j ), as shown in Equations (3) and (4) reflects the inverse square field from the Biot-Savart law as in [9]. This model represents the values of measured signals (in µV) that would have been obtained using electrodes located on the patient's scalp and is shown below: n ∈ R m×1 represents the Gaussian white noise process with some specified covariance. Thus, f j,i ∈ R is the signal strength of the ith dipole received at the jth sensor and g j,i is given as where d ji shows the distance between the ith dipole and the jth sensor.
Furthermore, denote G i ∈ R m×3 as the gain matrix for the ith dipole with m sensors, and constant brain conductivity ζ (µS/cm) [11]. Thus: The EEG signals, as received by the m sensors, are then compactly represented as This model illustrates the relationship between the dipole properties and the collected EEG signals, considering the noise of the sample collecting process. Note that F ∈ R m×n .
To account for the fact that the conductivity of the brain material is non-uniform, we propose a piece-wise constant conductivity to account for the soft matter, as well as the skeletal tissue before the dipole signal is received at the sensor. This is modeled as follows to present a more realistic conductivity: In what follows, each of the terms introduced in Equation (7) such as ζ j,i and ρ j,i will be elaborated upon.

•
Variable conductivity (ζ j,i ) within the head model: It is assumed that a dipole can be located anywhere in the brain (half-sphere head model), and the sensors are located on the patient's scalp. This means the generated signals from dipoles pass through different parts of the head, such as the brain's soft tissue, the spongy bone, the compact bone part, and the skin, to reach the sensors. Since each of these materials have different conductivities, the signal conduction is affected accordingly. In order to model this changing conductivity, the ζ j,i is modeled as a uniformly distributed random variable between 0.1 and 0.9 all the way through the head (including the soft parts to the hard parts such as bone). • Sensor distribution: Several methods are available to locate the electrodes on the patient's scalp, usually named by numbers indicating the standard locations. For instance, the two most popular sensor distributions are 10-10 and 10-20 [22]. This study aimed to find a precise estimation of the dipole properties, regardless of how the sensors are distributed. In order to fulfill this goal, a random distribution of sensors is considered in this paper. Note that the sensors are uniformly randomly distributed over the scalp. Figure 3 illustrates a sample of uniform random distribution for 128 sensors. • The adjacency of the dipoles to the sensors and the varying received signal strength (RSS) as a function of the location of the dipoles with respect to the sensors (ρ j,i ): the accuracy of the collected signals is directly related to the distance between the source and the sensor. In this paper, the location of sensors is assumed to be uniformly randomly distributed. Thus, the distance between the dipole to each sensor can be varied from 0 (co-located source and sensor) to the brain's diameter (i.e., 2r = 20 cm in this paper). Thus, greater proximity to the sensor translates into a better RSS value of the EEG signal. In order to consider this, a distance-based signal strength attenuation term ρ j,i is modeled as follows where α and β are constants and chosen such that the coefficient ρ j,i lies within 0.01 and 1.0. Note that the units of ρ j,i are the same as those of the charge, i.e., Coulomb. Note that these coefficients eliminate the effect of weak and noisy data and keep the high strength signals. In order to illustrate this concept, consider the following example. For the closest distance (the dipole is right under the sensor), the collected value will be the same as the real value generated by the dipole, which means that the coefficient is equal to exactly 1.0. With the same approach, for the longest distance (the brain diameter), the collected value will be 0.01 times the real value, which means that one can neglect it. The synthetic signals generated using the aforementioned EEG measurement model are shown in Figure 4.
In this paper, to illustrate the effectiveness of the measurement model as well as the estimation algorithm, it is assumed that the dipoles are fixed in orientation and magnitude. The head is modeled as a half-sphere with a diameter of 20 cm as it is suggested in [23]. Figure 5 shows a comparison of BESA (https://www.besa.de/products/besa-research/ features/head-model-selection/, accessed on 15 January 2023) data, with specific head model parameters (conductivities) shown in Figure 5a against that synthesized using our method assuming that ζ is uniformly distributed between 0.2 and 0.4 S/m. The results are in very good agreement. Henceforth, all synthetic data were generated using the proposed mathematical model following this verification.

Problem Statement
Succinctly stated, the problem that is solved in this paper is as follows. "Assuming multiple dipole sources (n) with fixed orientations and strengths, and given m noisy synthetic EEG signal measurements, we seek to find an accurate estimate of the unknown parameters that characterize dipoles, such as their distinct locations, orientations, and the strengths of the dipoles using a simple phenomenological measurement model provided in Equation (6)".

Solution Methodology
The parameters to be estimated for the i-th dipole are denoted compactly as the vector p i ∈ R 7×1 : The above-mentioned unknowns all appear in the right side of Equation (6). Thus, the inverse problem is solved to obtain the unknown components in the p vector.
Among the different approaches to solving an inverse problem, this paper chooses an optimization method where the L 2 norm of the estimation error is minimized. The estimation error is defined as the difference in the real measured signals and those predicted by the estimated values of the unknown parameters using the measurement model in Equation (6).p i = l x il y il z iμ x iμ y iμ z iŝ i T The predicted value of the EEG signal from the parameter estimates is given bŷ The objective of the inverse problem is to minimize the weighted L 2 norm of the measurement residuals, shown as: where W is a symmetric weighting matrix. In this paper, W is chosen to be an identity matrix. The cost function is augmented with the unit norm constraint of the dipole orientation, i.e., μ i 2 = 1.
where λ i is the Lagrange Multiplier corresponding to the ith dipole orientation constraint in the equation above. Note that the procedure for solving this problem would be to set up the necessary conditions (from the gradient of the cost), and determine an update for the parameters from one iteration to the next using the gradient and the Hessian (second derivative of the cost function with respect to the decision variables), and the application of the Karush-Kuhn-Tucker (KKT) conditions. These procedures are built into a nonlinear solver such as 'fmincon' in the Optimization Toolbox of MATLAB, which allows for an explicit specification of the cost to be minimized, nonlinear constraints that the decision variables satisfy, as well as any bounds on the decision values that need to be respected.
The results obtained are discussed in the next section.

Simulation Results
Before presenting the detailed simulation results, a thorough sensitivity and characteristics analysis of the solution to the inverse problem with regard to the number of sensors to be employed and the number of samples required to reliably provide a converged solution was performed. The signals generated for the analysis have an average signal-to-noise ratio (SNR) of approximately 20. Table 1 provides the results of the simulation in this study. As shown in the table, the source localization process uses an initial guess and iteratively converges to the actual values. As shown in Table 1, the estimated value is very accurate for all three selected dipoles, indicating that the source localization algorithm is consistent. To study the sensitivity of the estimation algorithm to the number of sensors and samples, the simulation is performed for different numbers of sensors, as well as samples for three different dipoles, and the results are shown in Figures 6-8. These 3D plots show the changes in the total error percentage of the estimated values in terms of the increasing number of sensors and samples simultaneously. As expected, by growing the number of samples and sensors, the amount of collected data increases. As a result, the estimation of unknown values is more accurate. Moreover, one can determine the least value of error and the corresponding number of samples and sensors. For instance, in this simulation, 48 sensors and 350 samples gave the least error percentage.

Model Sensitivity and Characteristics Analysis Using Synthetic Data and the Spherical Head Model
In order to study the effect of the distance between dipoles and sensors, three different locations are chosen for the simulation, one deeper in the brain, one closer to the scalp and sensors, and one dipole is chosen somewhere in between the locations of the other two. Based on the estimation algorithm described in the previous section, sensors receive weaker signals from the deeper dipoles. This leads to a higher amount of estimation error for the deeper dipoles, i.e., dipole #1 in this study. Considering this information regarding the location of the dipoles, one can interpret Figures 6-8 more accurately. This case can be observed by comparing the magnitude of the estimation error for each dipole from these three figures.    9-13 present the error percentage for estimating each variable in the parameter vector p i separately, where the error bar is the confidence interval value. It was mentioned previously that the electrodes are randomly located on the patient's scalp to make sure this simulation is working correctly regardless of the sensors' location. However, in some cases, the sensors could be located somewhere far from the dipole, and as a result, the collected data are noisier and weaker. Clearly, this happens when the number of sensors is excessively low and they cannot cover the scalp adequately. As a result, the confidence intervals for the low number of sensors are more significant. In other words, if the small number of sensors are located in some place near the dipoles, the result is still acceptable. Figures 9-11 show the error percentage for each variable. There is a lower error for each variable in terms of the number of sensors. Thus, it is not possible to choose the optimum number of sensors only using the error in each variable. To address this problem, an average error of all the parameters is studied, as shown in Figure 12 which provides a more robust estimate for the optimum number of sensors required. We note that (a) the minimum error percentage occurs when around 19 sensors are used; and (b) the error percentage reduces as the number of sensors increase (in general). This occurs because a higher number of sensors cover more of the scalp area, and as a result, the collected EEG data have more information. To summarize, Figure 12 represents a converged and consistent value of error after a specific number of sensors (20 sensors) are used. This means a specific amount of collected data are enough to make the model equations well defined, and the estimation problem well posed to find a consistent estimate. Thus, one can use the minimum number of sensors, as indicated by these numerical experiments, and conduct a simpler actual experiment (in terms of cost and computational burden). In order to further test this hypothesis, this simulation was performed using up to 100 sensors, and the result is provided in Figure 13.     It was previously discussed that 19 sensors would be good enough for this study to estimate the dipole properties. Similarly, a sensitivity study with respect to the number of samples was carried out. Figure 14 shows the changes in error percentage by increasing the number of samples. It was determined that 400 samples were sufficient to estimate all the dipoles' parameters. The effect of distance between the sensor and a dipole is one of the important considerations in this study. In order to see this, two different sensors are considered, one closer to a specific dipole, and the other one far from the same dipole, as was previously shown in Figure 1. The collected EEG signal from these two dipoles is shown in Figure 15. As can be seen from this figure, the sensor collected a stronger EEG signal from the closer dipole (shown in blue) rather than the dipole far from the sensor (shown in red). In conclusion, the proposed mathematical model for the dipole is very general and captures the key features of signal propagation from the source to the sensor. This simple model allows significant flexibility in terms of medium dependent conduction, signal strength variation due to distance, as well as other unmodeled noises.

Model Performance for Source Localization Using EEGLAB
In this study, the half-sphere head model is used to simplify the simulation. Furthermore, different parts of the brain's conductivity are assumed to be random numbers between 0.1 and 0.9 as mentioned before in Section 2.
In order to check how precise these assumptions are, one can compare the derived results provided by the present algorithm with the generated results from available open source software such as BESA, Cartool, EEGLAB, etc. [24]. In this study, the EEGLAB [25] software is chosen, with several options to choose the head model. The MNI head model is one of the most accurate ones selected in this study to compare with the spherical head model. The first step in using EEGLAB is to import the same generated signal used in the previous section. The location of sensors in the MATLAB-generated EEG signal are the standard 64 locations of BioSemi (known as 10-20 standard system), and the same locations will be considered in EEGLAB. The MNI head model in EEGLAB has a head radius of 85 mm. Note that, for having a more accurate comparison, the head model in the previous sections that considered a half sphere is also modified to an 85 mm radius.
Note that prior to the source localization step using clinical data in EEGLAB, the user has to reject certain data by visual inspection. The proposed algorithm (discussed previously) on the other hand does not require this manual step. Furthermore, since only dipole activities generate the available EEG data, there is no need to denoise or cancel any other 'artifact' in this case. The DIPFIT toolbox is used in EEGLAB to generate the source localization results. This MNE-based toolbox is designed by Robert Oostenveld and is available at Fieldtrip [26]. Applying the DIPFIT function in EEGLAB gives the dipole's location in Talairach coordinates, as shown in Figure 16a. Furthermore, Figure 16a,b are located next to each other to show the similarity of the source localization from the proposed algorithm and the result obtained from EEGLAB. Given the differences between these two approaches and those associated with the representations (Talairach coordinates vs. the Cartesian coordinates), the results are in close agreement. Besides the visual presentation of the results, Table 2 provides the numerical results for the three dipoles located in the head model. The error mentioned in this table is the distance between the actual source and the source localization result in mm, which is defined in Equation (11). Note that for all three dipoles, the error is significantly lower for the present algorithm compared to the EEGLAB result. We thus verify that the constrained nonlinear leastsquares-based source localization method's performance agrees well with that produced by EEGLAB.

Characteristics Analysis Using the Clinical Data and the MNI Coordinates
Now that the efficacy of the proposed algorithm has been proven in the previous sections, one can assess this algorithm using clinical sample data. Therefore, a suggested dataset from EEGLAB is used in this section [20,27]. This dataset belongs to a visual attention task. In this task, each event is a three-second time interval where the subject should press a button right after seeing a square on a screen in front of them. Given the quality of the dataset and the presence of bad quality signals, it is suggested to select good quality time intervals rather than using the whole dataset. Thus, in this study, two different time intervals are selected to quantify the source localization algorithm. As illustrated in Figure 17, one of the selected events is from 113 to 116 s, and the other is from 146 to 149 s. The selected events are utilized to generate the source localization results on both the EEGLAB software and the introduced algorithm. The EEGLAB source localization algorithm suggests the number of dipole locations to be as many as the channel numbers. In other words, for this specific dataset that was recorded by 30 EEG signals, EEGLAB found 30 possible dipoles that generate this signal. However, only the results with residual variances (RVs) less than 15% are acceptable according to EEGLAB [28]. The residual variance in EEGLAB software is defined based on matching the dipole projection to the EEG electrodes. Based on this definition, it is also important to mention that smaller numbers of RV indicate the most accuracy, and the results are more reliable. Figures 18 and 19 represent the source localization results for the experiment that started at 113 s and finished at 116 s. Figure 18 illustrates all the dipole estimations simultaneously. Given this general presentation of results, one can see the similarity between the EEGLAB result and the introduced algorithm result. Note that, as mentioned earlier, Figure 18b is the general result by EEGLAB where all the dipole estimations are not necessarily correct. Figure 18c is the final result where any estimation with the RV over 15% is filtered out.  EEGLAB usually utilizes an MRI scan to show the source localization result, making the comparison less convenient. This suggests providing another method to compare the results, and Figure 19 is provided to compare the dipole estimations individually in the same brain map using the Talairach coordinates.
Among   Finally, the algorithm is also tested on clinical data collected from a patient with an active seizure case. The Temple University Hospital provided the EEG data used in this section as their open source database [29]. These data were collected by 19 EEG sensors with the international 10-20 EEG electrodes, as shown in Figure 22b.
After applying the EEG source localization to the seizure EEG data, the comparison in Figure 23 shows that the EEGLAB provides 19 estimated dipoles, and only four have an RV > 15%. However, the proposed algorithm results in 35 dipole estimations, including the four dipole results from EEGLAB. Figure 24 presents the four similar results compared individually.
Given a more significant number of dipole estimations from the MATLAB simulation compared to the EEGLAB, the accuracy of this result is illustrated by providing the channel data. As shown in Figure 22, the highest brain activity occurs around sensors Fp1, Fp2, and F3, which matches both the EEGLAB and MATLAB simulation result. On the other hand, there are other channels with noticeable brain activity, such as Cz, C3, Pz, and P4. This means there could be other sources of brain activity around the central area, even though they are not as strong as the frontal head area. These sources are also discerned by the proposed algorithm's result, as shown in Figure 23a. In other words, the introduced algorithm is more sensitive to all the brain signals and covers a broader range of source localization.

Conclusions
This study introduces a constrained nonlinear least-squares algorithm to estimate the dipole properties based on the collected EEG signals. The mathematical model for a signal from the dipole is shown to be very flexible and includes features that model piece-wise conductivity, varying received signal strength (RSS), and random sensor distribution on the scalp. This model flexibility enables one to solve the inverse problem for many different head types and/or numbers or sensors and samples. To give an illustration of this fact, this algorithm is applied to two different types of EEG data: (1) synthetic data generated by the forward model; and (2) clinical data. The results emphasize the accuracy and performance of the present algorithm for both EEG signal categories. The source localization error for the synthetic data is less than 0.1 percent. In the case of testing the clinical data, this algorithm works as accurately as EEGLAB, which is illustrated in the figures. The results also show that the presented algorithm performs extremely well for the localization of multiple dipoles and is verified with the results obtained from EEGLAB. Notably, the introduced algorithm can be more sensitive to different sources in the head model. As a result, for some cases, such as the seizure data utilized in this study, the proposed source localization algorithm can find more active dipoles compared to EEGLAB. Future work will focus on further optimizing the number of sensors, decreasing the computation load of source localization algorithms where the FEM and BEM are used, and prototyping the algorithm on a chip that can be embedded in an EEG helmet.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.