brain sciences Article
Feasibility of Equivalent Dipole Models for Electroencephalogram-Based Brain Computer Interfaces Paul H. Schimpf
ID
Department of Computer Science, Eastern Washington University, Cheney, WA 99004, USA; [email protected]; Tel.: +1-509-359-6908 Received: 5 August 2017; Accepted: 13 September 2017; Published: 15 September 2017
Abstract: This article examines the localization errors of equivalent dipolar sources inverted from the surface electroencephalogram in order to determine the feasibility of using their location as classification parameters for non-invasive brain computer interfaces. Inverse localization errors are examined for two head models: a model represented by four concentric spheres and a realistic model based on medical imagery. It is shown that the spherical model results in localization ambiguity such that a number of dipolar sources, with different azimuths and varying orientations, provide a near match to the electroencephalogram of the best equivalent source. No such ambiguity exists for the elevation of inverted sources, indicating that for spherical head models, only the elevation of inverted sources (and not the azimuth) can be expected to provide meaningful classification parameters for brain–computer interfaces. In a realistic head model, all three parameters of the inverted source location are found to be reliable, providing a more robust set of parameters. In both cases, the residual error hypersurfaces demonstrate local minima, indicating that a search for the best-matching sources should be global. Source localization error vs. signal-to-noise ratio is also demonstrated for both head models. Keywords: brain–computer interface; electroencephalogram; equivalent source model
1. Introduction The vast majority of non-invasive brain computer interfaces (BCI) based on the electroencephalogram (EEG) focus on temporal events such as sensorimotor rhythms and evoked potentials [1]. Few attempts have been made to make use of source localization other than by strategic location of the EEG recording electrodes. This is perhaps due to the fact that temporal information in the EEG is well understood, reproducible, and not overly challenging to extract in terms of the required real-time processing power. In contrast, inverting the EEG to localize brain activity is known to be an ill-posed problem [2], and is highly underdetermined when linearized with a lead-field approach [3]. While there are a number of algorithms to circumvent this problem, they tend to be compute-intensive and might represent some challenges for real-time computation in a portable system. In any event, few attempts to make use of equivalent source localization for EEG-based BCI have been published. Qin, L. et al. reported on a pilot study for BCI based on equivalent source location using off-line EEG data, wherein subjects were asked to imagine right or left hand movement based on a visual cue [4]. The time window and frequency band of the event-related desynchronization (ERD) [5] were hand-selected from the data, pre-processed via independent component analysis [6], and inverted to a source map using a regularized weighted minimum norm [7]. A single best-fit equivalent dipole was also found using the Nelder–Mead Simplex method [8]. A three concentric sphere head model was used to represent the domain [9]. All currents within 80% of the largest were summed over each hemisphere, and if the sum on the ipsilateral side was larger than the contralateral side, the source mapping method Brain Sci. 2017, 7, 118; doi:10.3390/brainsci7090118
www.mdpi.com/journal/brainsci
Brain Sci. 2017, 7, 118
2 of 11
was judged to have correctly classified the imagined hand movement. If the single best-fit equivalent dipole was located in the ipsilateral hemisphere of the brain, the equivalent dipole approach was judged to have correctly classified the imagined hand movement. Single-trial classification accuracy of 78.9% (source map) and 80.6% (equivalent dipole) were reported on a single test subject. Because the time and frequency data were hand-selected, no detection algorithm was discussed and thus no information was available regarding false positives. Nevertheless, the results indicated promise for the idea of using the localization of equivalent dipoles from EEG for BCI. No information was included on processing time, although it was stated that the approach could be implemented in an on-line BCI application with little alteration. While the noise subspace was estimated and discarded, and the source signal-to-noise ratio (SNR) was reported to be better than the scalp SNR, no numerical estimates of the SNR were reported. A follow-up study [10] expanded the analysis to two equivalent dipoles. Manual selection of the time windows was again applied and Nelder–Mead Simplex was again used for the search. The method was judged to have correctly classified the imagined hand movement if the stronger dipole was located in the ipsilateral hemisphere. An average classification rate across four subjects was 80%, a promising result but not higher than the preliminary study using a single dipole. It is important to note that the equivalent dipole locations, orientations, and amplitudes with this approach are not intended to represent an accurate model of brain activity. They instead represent an amalgam of brain activity as represented by a few dipolar sources. When a lead-field approach is used the forward problem is represented by a linear system as [3]: Ls = v
(1)
where L is the (m × 3n) lead-field, or forward transfer, matrix, s is the (3n × 1) vector of dipolar source amplitudes, with each dipolar source represented by three orthogonal components of the dipole moment, and v is the (m × 1) vector of EEG samples. In the most general case, the number of unknown source amplitudes (3n) is far greater than the number of EEG measurements (m), representing a severely underdetermined problem that is also ill-posed [2]. There are a number of algorithms that select from the infinity of solutions via various assumptions. When the number of non-zero sources is known to be, or can be restricted to be, some number less than 3n, the problem becomes underdetermined, with a single solution in the least-squares sense: −1 sq = LT LT q Lq qv
(2)
where q indicates some subset of the n candidate sources along with the corresponding columns of L. One must then search for the subset q that minimizes the residual minkLq sq − vk2 q
(3)
The question in this case is not the accuracy to which the fitted dipoles represent actual brain activity, but rather whether the location of the best fitting dipoles are useful as classification parameters in a BCI algorithm. This will depend on a subject’s ability, perhaps with training, to control the equivalent dipole locations. Because of the ill-posed nature of the problem, a reasonable precursor question is the extent to which the actual best fitting sources can be found given the head model being used. That is the question investigated here. 2. Materials and Methods 2.1. Head Models The reproducibility of a single known equivalent dipole in the face of noisy EEG measurements was investigated in both a four-shell spherical head model [11] and a realistic head model.
Brain Sci. 2017, 7, 118
3 of 11
The four-shell model has been used in several other studies (e.g., [12–15]) and is composed of concentric Brain Sci. 2017,spheres 7, 118 of varying conductivity representing the brain, cerebral-spinal fluid (CSF), 3skull, of 11 and scalp, as shown in Figure 1. The origin of the model is the center of the spheres, and the radii of Brain Sci. 2017, 7, 118 were 3 ofof 11the boundaries 7575 mm respectively. TheThe modeled conductivities the tissue tissue boundaries wereset setatat63, 63,65, 65,71, 71,and and mm respectively. modeled conductivities of layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19]. the layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19]. the tissue boundaries were set at 63, 65, 71, and 75 mm respectively. The modeled conductivities of the layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19].
Figure 1. of the spherical model. Figure 1. Cross-section Cross-section through the center of the thespherical sphericalmodel. model. Figure 1. Cross-sectionthrough throughthe the center center of
10001000 candidate source distributed pseudo-uniformly[20] [20] overthethe upper candidate sourcelocations locationswere were distributed distributed pseudo-uniformly upper 1000 candidate source locations were pseudo-uniformly [20]over over the upper hemisphere of the cortical surface, with approximate spacing of 4.8mm. mm.256 256EEG EEG electrodes were hemisphere cortical surface, withan spacing of were hemisphere of of thethe cortical surface, withan anapproximate approximate spacing4.8 of 4.8 mm. 256electrodes EEG electrodes likewise distributed over thetheupper surface.The Thelocations locationsof of candidate likewise distributed over upperhemisphere hemisphere of of the the scalp scalp surface. candidate were likewise distributed over the upper hemisphere of the scalp surface. The locations of candidate sources electrodes shownininFigure Figure2.2.The The forward forward problem source was sources andand electrodes areare shown problemfor foreach eachcandidate candidate source was sources and electrodes are shown in Figure 2. The forward problem for each candidate source was obtained from the analytic solution in [11] and used to form a 256 × 3000 lead-field, with three obtained from the analytic solution in [11] and used form a 256 × 3000 lead-field, with three obtained from dipole the analytic solution inper [11]candidate and used to form a 256 × 3000 lead-field, with three orthogonal moment directionsper source location. orthogonal dipole moment directions candidate source location. orthogonal dipole moment directions per candidate source location.
Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources in Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources blue, EEG electrodes in red). in blue, EEG electrodes in red).
Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources TheEEG realistic head model is based on classified magnetic resonance images of a human head. The in blue, electrodes in red). tissues were classified according to the method discussed in [21]. The model contained 11 tissue types
The realistic head model is based on classified magnetic resonance images of a human head. The tissues were classified according to the method discussed in [21]. The model contained 11 tissue types
Brain Sci. 2017, 7, 118
Brain Sci. 2017, 7, 118
4 of 11
4 of 11
The realistic head model is based on classified magnetic resonance images of a human head. The tissues were classified according to the method discussed in [21]. The model contained 11 tissue with resistivities obtained from [16]. EEG observation locations were obtained from the 82 points of types with resistivities obtained from [16]. EEG observation locations were obtained from the 82 points an extended EEG 10–20 layout, visually interpolated with an additional 63 points, resulting in 145 of an extended EEG 10–20 layout, visually interpolated with an additional 63 points, resulting in EEG channels. The candidate source space consisted of 3035 locations sampled from the cortical 145 EEG channels. The candidate source space consisted of 3035 locations sampled from the cortical surfacesurface at a resolution of 4.0 × 4.0 × 3.2 resulting inina alead matrixofofdimension dimension × 9105, at a resolution of 4.0 × 4.0 × mm, 3.2 mm, resulting leadfield field matrix 145145 × 9105, with three orthogonal source directions perper location. matrixwas wascalculated calculated using with three orthogonal source directions location.The Thelead-field lead-field matrix using thethe finite element method [15,22]. This model is illustrated in Figure 3. This model, along with a finite finite element method [15,22]. This model is illustrated in Figure This model, along with a finite element for non-conforming meshes, is available fromthe theauthor author at at [23]. element solversolver for non-conforming meshes, is available from
FigureFigure 3. Realistic head model with cutaway showing internal tissues. EEG electrodes in red. 3. Realistic head model with cutaway showing internal tissues. EEG electrodes in red.
2.2. Source Localization ErrorError 2.2. Source Localization As discussed in in the the primary question is the theoretical extent to which As discussed theintroduction, introduction, the primary question is the theoretical extent to which equivalent equivalent dipoles be reliably located the EEG given the mathematically ill-posed nature of dipoles can becan reliably located from the from EEG given the mathematically ill-posed nature of the inverse problem. This question was investigated using random on the twoon head Formodels. each the inverse problem. This question was investigated usingtrials random trials themodels. two head model, ground-truth dipolesdipoles were generated at randomat source locations, with random orientations. For each model, ground-truth were generated random source locations, with random The forward of Equation was solved EEGsolved measurements normally distributed noise orientations. The model forward model of(1)Equation (1)forwas for EEGand measurements and normally was added. The best fitting to the noisy EEG wasnoisy then found by performing distributed noise was added. Thedipole best fitting dipole to the EEG was then foundanbyexhaustive performing search for the location that minimized the residual error in Equation (3). The mean and standard an exhaustive search for the location that minimized the residual error in Equation (3). The mean and deviation of the Euclidean distance between the actual and inverted dipole locations was then plotted standard deviation of the Euclidean distance between the actual and inverted dipole locations was against signal-to-noise (SNR) ratio. Surface plots of the residual error in Equation (3) over all candidate then plotted against were signal-to-noise (SNR) ratio. Surface plots of theinresidual error ininsight Equation source locations also examined for the presence of local minima, order to provide as to (3) over all candidate source locations were also examined for the presence of local minima, in order to whether a global search algorithm is required for this approach. provide insight as to whether a global search algorithm is required for this approach. 3. Results and Discussion 3.1. Spherical Head Model Figure 4 illustrates localization error vs. SNR for the spherical head model for 250 randomized
Brain Sci. 2017, 7, 118
5 of 11
3. Results and Discussion 3.1. Spherical Head Model Figure 4 illustrates localization error vs. SNR for the spherical head model for 250 randomized trials at each SNR. The possible range of localization error depends on the location of the ground truth source in the hemisphere of possibilities (see Figure 2) which was selected randomly for each5trial. Brain Sci. 2017, 7, 118 of 11 Figure 4 shows that the localization degrades dramatically for SNRs lower than 70 dB, which is an unrealistically largelarge SNR SNR for single trial EEG The compression of the standard an unrealistically for single trialrecordings. EEG recordings. The compression of the deviation standard at lower SNRs is explained by the factby that finite and inverted aresources restricted deviation at lower SNRs is explained thethe factdomain that theis domain is finite andsources inverted are to the cortical surface, which creates an upper limit to the possible localization error. As the noise restricted to the cortical surface, which creates an upper limit to the possible localization error. As the level (lower(lower SNR),SNR), the mean localization approaches an asymptote givengiven by the case noiseincreases level increases the mean localization approaches an asymptote byworst the worst scenario of a completely random guess.guess. As theAs center (mean)(mean) of the of histogram increases, the width case scenario of a completely random the center the histogram increases, the of the histogram is compressed because of the upper limit to the error. Thus, when the noise is low, width of the histogram is compressed because of the upper limit to the error. Thus, when the noise is both the mean and standard deviation are small (the localization is reliable), both grow with the noise low, both the mean and standard deviation are small (the localization is reliable), both grow with the level width the histogram begins tobegins encroach on the upper limit on error, causes the noiseuntil levelthe until the of width of the histogram to encroach on the upper limitwhich on error, which histogram to then shrink somewhat the mean as continues to approach worst casethe scenario. causes the width histogram width to then shrinkassomewhat the mean continuesthe to approach worst As aid to interpretation the trend in Figure 4, a completely random guess as to source location caseanscenario. As an aid to of interpretation of the trend in Figure 4, a completely random guess as to when the true location is at the apex (north pole) of the cortical surface has a mean of 53 mm and a source location when the true location is at the apex (north pole) of the cortical surface has a mean of standard deviation of 25 mm. When the true location is anywhere on the equator, the mean error of a 53 mm and a standard deviation of 25 mm. When the true location is anywhere on the equator, the completely guess israndom 85 mm with standard deviation of 25 mm. mean error random of a completely guessa is 85 mm with a standard deviation of 25 mm.
Figure 4. Figure 4. Localization Localization error error vs. vs. signal-to-noise signal-to-noise (SNR) (SNR) ratio ratio for for the the spherical spherical model. model. Mean Mean and and standard standard deviation are shown for for 250 250 randomized randomized trials trials at at each each SNR. SNR. Localization Localization is poor below below 70 70 dB dB (the (the deviation are shown is poor radius of the cortical sphere is 63 mm). radius of the cortical sphere is 63 mm).
A fundamental problem is that constraining the inverse to an underdetermined system does not A fundamental problem is that constraining the inverse to an underdetermined system does not avoid the ill-posed nature, and small perturbations in the measurement can produce widely different avoid the ill-posed nature, and small perturbations in the measurement can produce widely different localizations. In this case, however, the localization ambiguity occurs only for the azimuth of the localizations. In this case, however, the localization ambiguity occurs only for the azimuth of the inverted source location. This is illustrated for noise-free data in Figure 5, which plots the residual inverted source location. This is illustrated for noise-free data in Figure 5, which plots the residual from Equation (3) for all candidate sources against the azimuth (theta) and elevation (phi) of the from Equation (3) for all candidate sources against the azimuth (theta) and elevation (phi) of the source. source. All candidate sources reside at the same radius (63 mm) in this model. This shows that there All candidate sources reside at the same radius (63 mm) in this model. This shows that there are a are a number of good fits to the EEG data (i.e., leaving a small residual) that are at different azimuths, number of good fits to the EEG data (i.e., leaving a small residual) that are at different azimuths, but but only one elevation. only one elevation.
Brain Sci. 2017, 7, 118 BrainSci. Sci.2017, 2017,7,7,118 118 Brain
6 of 11 of11 11 66of
Figure 5. candidate sources for for an noise-free case. All good sources Figure 5. 5. Residual Residualat atall all candidate sources forillustrative an illustrative illustrative noise-free case. Allcandidate good candidate candidate Figure Residual at all candidate sources an noise-free case. All good (low residual) reside at the same elevation (phi), with a range of azimuths (theta) that covers the entire sources(low (lowresidual) residual)reside resideatatthe thesame sameelevation elevation(phi), (phi),with withaarange rangeof ofazimuths azimuths(theta) (theta)that thatcovers covers sources 2π radians. theentire entire2π 2πradians. radians. the
furtherillustrates illustratesthis thispoint point by plotting the elevation error vs. SNR, SNR, which performs Figure 666 further byby plotting thethe elevation error vs. vs. SNR, which performs well Figure illustrates this point plotting elevation error which performs wellover over much larger SNR range than source location. over a much larger SNRSNR range thanthan source location. well aamuch larger range source location.
(SNR)ratio Figure Figure6.6.Elevation Elevationerror errorvs. vs.signal-to-noise signal-to-noise(SNR) ratiofor forthe thespherical sphericalmodel. model.Mean Meanand andstandard standard deviation are shown for 250 randomized trials at each SNR. The elevation range for candidate sources deviation are shown for 250 randomized trials at each SNR. The elevation range for candidate sources to90 90degrees. degrees. isis00to
Figure 77 provides provides additional additional insight insight into into the the azimuth azimuth ambiguity. ambiguity. ItIt plots plots the the location location and and Figure moment of a ground truth dipole source (red) along with the locations and moments of all candidate moment of a ground truth dipole source (red) along with the locations and moments of all candidate sourcesthat thatprovide providean anEEG EEGresidual residualwithin within1% 1%(blue). (blue).The TheEEG EEGisisplotted plottedunderneath underneaththe thedipoles dipolesas as sources
Brain Sci. 2017, 7, 118
7 of 11
Figure 7 provides additional insight into the azimuth ambiguity. It plots the location and moment of a ground truth dipole source (red) along with the locations and moments of all candidate sources Brain Sci. 2017, 7, 118 7 of 11 that provide an EEG residual within 1% (blue). The EEG is plotted underneath the dipoles as an aid to understanding how thesehow different could produce same EEG on theEEG surface of the model. an aid to understanding thesedipoles different dipoles couldthe produce the same on the surface of This shows that there are a number of good candidate sources that can reproduce the surface EEG, the model. This shows that there are a number of good candidate sources that can reproduce the distributed around the entire circumference of the cortex, all differing orientations. only the surface EEG, distributed around the entire circumference ofwith the cortex, all with differingThus orientations. elevation can provide reliable information in this case. It is the spherical symmetry of the model that Thus only the elevation can provide reliable information in this case. It is the spherical symmetry of produces theambiguity azimuth of the modelthis thatambiguity producesin this inthe theinverted azimuthsources. of the inverted sources.
Figure 7. Dipolar sources that produce an Electroencephalogram (EEG) within 1% of the EEG produced Figure 7. Dipolar sources that produce Electroencephalogram (EEG) within at1% the EEG by the true source (red). The surface EEGan produced by the true source is plotted theofbottom of produced the figure. by the true source (red). The surface EEG produced by the true source is plotted at the bottom of the figure.
Taken together, Figures 5–7 indicate that only the elevations of inverted sources are reliably Taken together, Figures 5–7 indicate that only the elevations of inverted sources are reliably invertible from the EEG when a spherical head model is used. Because the azimuth of the location is invertible from the EEG when a spherical head model is used. Because the azimuth of the location is unreliable, including it in the parameter set could confound a BCI algorithm when a spherical head unreliable, including it in the parameter set could confound a BCI algorithm when a spherical head model is used. model is used. 3.2. Realistic Head Model 3.2. Realistic Head Model Figure 8 illustrates localization error vs. SNR for the realistic head model. It performs much better Figure 8 illustrates localization error vs. SNR for the realistic head model. It performs much over a range of SNRs than the spherical head model, with reasonable results down to 10 dB. The use of better over a range of SNRs than the spherical head model, with reasonable results down to 10 dB. a realistic head model thus provides at least three reasonable BCI parameters per inverted source: x, y, The use of a realistic head model thus provides at least three reasonable BCI parameters per inverted and z, or equivalently, r, theta, and phi. source: x, y, and z, or equivalently, r, theta, and phi.
Brain Sci. 2017, 7, 118
8 of 11
Brain Sci. 2017, 7, 118
8 of 11
Figure 8.Figure Localization errorerror vs. vs. signal-to-noise ratioforfor realistic Mean and 8. Localization signal-to-noise (SNR) (SNR) ratio thethe realistic head head model.model. Mean and standard deviation are shown for 250 randomized trials at each SNR. Localization is much better than standard deviation are shown for 250 randomized trials at each SNR. Localization is much better than the spherical head model at lower SNRs. the spherical head model at lower SNRs. 3.3. Local vs. Global Minima
3.3. Local vs. Global Minima
For the spherical head model, Figure 5 also illustrates local minima in the residual error surface
For that thewould spherical model, also illustrates local minima in model, the residual error trap ahead localized searchFigure for the 5 minimum residual. In the realistic head the origin for surface that would trap asources localized searchin for the minimum residual. In the realistic head the origin candidate was placed the Cerebellum in order to compare approximately to the model, coordinate system of the spherical model. A difference is that candidate sources in the realistic model do not for candidate sources was placed in the Cerebellum in order to compare approximately to the all fall at theof same The residual hypersurface thuscandidate has three dimensions and cannot be model coordinate system theradius. spherical model.error A difference is that sources in the realistic shown in its entirety on a single plot. Figure 9 shows the projection of a residual hypersurface onto the do not all fall at the same radius. The residual error hypersurface thus has three dimensions and azimuth (theta) and radial directions, indicating the presence of local minima for the realistic head cannot be shown its aentirety on a for single plot. Figure shows the projection of a residual model as well.inThus localized search the best-fitting dipoles9for either model cannot be expected hypersurface onto the (theta) and radial indicating the presence of local minima to be reliable, andazimuth a global search strategy should directions, be used. A simple a global search is to exhaustively compute the residual for all possibilities, for the realistic head approach model astowell. Thus a localized search for the best-fitting dipoles for either model p ), where n is the number of candidate sources (one-third which has a computational complexity of O(n cannot be expected to be reliable, and a global search strategy should be used. the number columns in the lead-field) and p is the number of sources allowed in the solution. This can A simple approach to a global search is to exhaustively compute the residual for all possibilities, become quickly intractable as p increases. On a Core i3-4360 CPU (3.7 GHz, Intel Inc., Santa Clara, CA, p which has a computational complexity of O(n whereLA, n USA) is therunning number of candidate sources (oneUSA) with 16 GB of DDR3 RAM (Crucial Inc.,),Gretna, Octave 4.2.1 ([24]) under third theMicrosoft numberWindows columns in the lead-field) and p is the number of sources 7 (Microsoft Corp., Redmond, WA USA), this search requiredallowed less than in 0.1the s forsolution. = 1 on thequickly realistic intractable head model. as A value of p = 2 required 202 s. The pseudo-inverses for different This canpbecome p increases. On a Core i3-4360 CPU (3.7 GHz, Intel Inc., Santa column subsets of the lead-field are trivially parallelizable. Thus p = 2 is within the range of real-time Clara, CA, USA) with 16 GB of DDR3 RAM (Crucial Inc., Gretna, LA, USA) running Octave 4.2.1 performance if a speedup on the order of 200 can be achieved. This is currently under investigation ([24]) under Microsoft Windows 7 (Microsoft Corp., Redmond, WA USA), this search required less using GPU parallelization. than 0.1 s for p = 1 on the realistic head model. A value of p = 2 required 202 s. The pseudo-inverses for different column subsets of the lead-field are trivially parallelizable. Thus p = 2 is within the range of real-time performance if a speedup on the order of 200 can be achieved. This is currently under investigation using GPU parallelization.
Brain Sci. 2017, 7, 118 Brain Sci. 2017, 7, 118
9 of 11 9 of 11
Figure 9. Residual at all candidate sources for an illustrative noise-free case in the realistic head model. Figure 9. Residual at all candidate sources for an illustrative noise-free case in the realistic head model.
3.4. Summary and Applicability 3.4. Summary and Applicability The results above show that the following should be taken into account when inverting equivalent The results above show that the following should be taken into account when inverting dipole locations from the EEG: equivalent dipole locations from the EEG: (a) When a spherical head model is used, the inverted azimuth of the dipole is unreliable in the face (a) When a spherical head model is used, the inverted azimuth of the dipole is unreliable in the face of even modest amounts of noise, although the elevation of the dipole is reliable. of even modest amounts of noise, although the elevation of the dipole is reliable. (b) All components components of of the the location location are are reasonably (b) All reasonably reliable reliable when when aa realistic realistic head head model model is is used. used. (c) The residual cost function for locating the dipole exhibits local minima, and thus a global (c) The residual cost function for locating the dipole exhibits local minima, and thus a global search search strategy should should be be used. used. strategy
For an an EEG-based EEG-based BCI, BCI, the the performance performance cannot cannot therefore therefore be be expected expected to to be be optimal optimal ifif aa spherical spherical For head model model is is used used with with the the full full inverted inverted dipole dipole location, location,or orwhen whenusing usingaalocalized localizedsearch searchstrategy. strategy. head It should be noted that it is not necessarily the case that a realistic head model for BCI be It should be noted that it is not necessarily the case that a realistic head model for subjectBCI be specific. Given the somewhat limited success of using a spherical model [4,10], along with subject-specific. Given the somewhat limited success of using a spherical model [4,10], along these with findings on location ambiguity, there is good to expect that any realistic model, model, such assuch [25], these findings on location ambiguity, there is reason good reason to expect that any realistic would improve BCI performance without needing to be subject-specific. as [25], would improve BCI performance without needing to be subject-specific. The findings findings here here are are specific specific to to source source localization localization from from the the EEG. EEG. Similar Similar concerns concerns may may or or may may The not exist for source localization from the magnetoencephalogram (MEG), which is related to the EEG not exist for source localization from the magnetoencephalogram (MEG), which is related to the EEG problem via via the the Biot–Savart Biot–Savart law law [15]: [15]: problem
Z J rˆ µ00 J × rˆ dv B B = 4 r22 dv 4π r
(4) (4)
Brain Sci. 2017, 7, 118
10 of 11
The use of a vector measurement (B) and the fact that J represents both volume and source currents may mitigate the source localization ambiguities found with the EEG, which is something not investigated here. It is worth noting that other modalities provide more direct opportunities to characterize the location of brain activity. As an imaging technology, functional magnetic resonance imaging (fMRI) [26] provides spatial information directly. Functional near-infrared spectroscopy (fNIRS) [27] provides some opportunities for obtaining spatial information through the strategic placement of emitter–detector pairs. The fNIRS technology appears promising and may well prove to supplant EEG as a common BCI approach when inexpensive fNIRS headsets become available as is the case for EEG (e.g., [28]). 4. Conclusions This work shows that equivalent dipole approaches to EEG-based BCI produce ambiguous results in source azimuth when a spherical head model is used, and can thus be expected to provide reliable information only in terms of source elevation. It also shows that a realistic head model does not suffer such localization ambiguity, providing a more robust set of parameters for potential BCI classification. It further shows that the EEG residual hypersurface demonstrates local minima, indicating that a search for best-fitting dipolar sources should be global. Acknowledgments: The cost for open access publication of this work was provided by Eastern Washington University. Conflicts of Interest: The author declares no conflict of interest.
References 1.
2. 3. 4. 5. 6. 7. 8. 9. 10.
11. 12. 13.
Sellers, E.W.; Arbel, Y.; Donchin, E.; Pfurtscheller, G.; McFarland, D.J.; Allison, B.Z.; Faller, J.; Neuper, C. Existing BCIs. In Brain-Computer Interfaces, Principles and Practices; Wolpaw, J.R., Wolpaw, E.W., Eds.; Oxford University Press: New York, NY, USA, 2012; pp. 215–249, ISBN 978-0-19-538885-5. Kabanikhin, S.I. Definitions and examples of inverse and ill-posed problems. J. Inverse Ill-Posed Probl. 2008, 16, 317–357. [CrossRef] Weinstein, D.; Zhukov, L.; Johnson, C. Lead-field bases for electroencephalography source imaging. Ann. Biomed. Eng. 2000, 28, 1059–1065. [CrossRef] [PubMed] Qin, L.; Ding, L.; He, B. Motor Imagery Classification by Means of Source Analysis for Brain Computer Interface Applications. J. Neural Eng. 2005, 2, 65–72. [CrossRef] [PubMed] Pfurtscheller, G. Event-Related Desynchronization; Elsevier: New York, NY, USA, 1999; pp. 303–325. Hyvärinen, A.; Oja, E. A Fast Fixed-Point Algorithm for Independent Component Analysis. Neural Comput. 1997, 9, 1483–1492. [CrossRef] Hansen, P.C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 1992, 34, 561–580. [CrossRef] Kowalik, J.; Osborne, M.R. Methods for Unconstrained Optimization Problems; Elsevier: New York, NY, USA, 1968. Rush, S.; Driscoll, D. Current distribution in the brain from surface electrodes. Anesth. Analg. 1968, 47, 717–723. [CrossRef] [PubMed] Kamousi, B.; Liu, Z.; He, B. Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis. IEEE Trans. Neural Sys. Rehabil. Eng. 2005, 13, 166–171. [CrossRef] [PubMed] Zhou, H. Computation of the potential distribution in a four-layer anisotropic concentric spherical volume conductor. IEEE Trans. Biomed. Eng. 1992, 39, 154–158. [CrossRef] [PubMed] Hoesk, R.S., Jr.; Sances, A.; Larson, S.J. The contributions of intracerebral currents of the EEG and evoked potentials. IEEE Trans. Biomed. Eng. 1978, 25, 405–413. [CrossRef] [PubMed] Stok, C.J. The influence of model parameters on EEG/MEG single dipole source estimation. IEEE Trans. Biomed. Eng. 1987, 34, 289–296. [CrossRef] [PubMed]
Brain Sci. 2017, 7, 118
14. 15. 16. 17. 18. 19.
20. 21.
22. 23. 24. 25.
26.
27. 28.
11 of 11
Schimpf, P.H.; Ramon, C.; Haueisen, J. Dipole models for the EEG and MEG. IEEE Trans. Biomed. Eng. 2002, 49, 409–418. [CrossRef] [PubMed] Schimpf, P.H. Application of Quasi-static Magnetic Reciprocity to Finite Element Models of the MEG Lead-Field. IEEE Trans. Biomed. Eng. 2007, 54, 2082–2088. [CrossRef] [PubMed] Geddes, A.; Baker, L.E. The specific resistance of biological material—A compendium of data for the biomedical engineer and physiologist. Med. Biol. Eng. 1967, 5, 271–293. [CrossRef] [PubMed] Rush, S.; Driscoll, D. EEG-electrode sensitivity—An application of reciprocity. IEEE Trans. Biomed. Eng. 1969, 16, 15–22. [CrossRef] [PubMed] Oostendorp, T.F.; Delbecke, J.; Stegman, D.F. The conductivity of the human skull: Results in vivo and in vitro measurements. IEEE Trans. Biomed. Eng. 2000, 47, 1487–1492. [CrossRef] [PubMed] Hoekema, R.; Huiskamp, G.J.M.; Wieneke, G.H.; Leijten, F.S.S.; van Veelen, C.W.M.; van Rijen, P.C.; van Huffelen, A.C. Measurement of the conductivity of skull, temporarily removed during epilepsy surgery. Brain Topogr. 2004, 16, 29–38. [CrossRef] Deserno, M. How to Generate Equidistributed Points on the Surface of a Sphere. Available online: http://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf (accessed on 22 May 2017). Shrinidhi, N.; Haynor, D.R.; Wang, Y.; Joregenson, D.B.; Bardy, G.H.; Kim, Y. An efficient tissue classifier for building patient-specific finite element models from X-ray CT images. IEEE Trans. Biomed. Eng. 1996, 43, 333–337. [CrossRef] [PubMed] Schimpf, P.; Haueisen, J.; Ramon, C.; Nowak, H. Realistic computer modeling of electric and magnetic fields of human head and torso. Parallel Comput. 1998, 24, 1433–1460. [CrossRef] GalerWin. Available online: https://sites.google.com/site/pschimpf99/home/software/galerwin (accessed on 14 September 2017). GNU Octave. Available online: https://www.gnu.org/software/octave/ (accessed on 14 September 2017). Lancaster, J.L.; Woldorff, M.G.; Parsons, L.M.; Liotti, M.; Freitas, C.S.; Rainey, L.; Kochunov, P.V.; Nickerson, D.; Mikiten, S.A.; Fox, P.T. Automated Talairach Atlas labels for functional brain mapping. Hum. Brain Mapp. 2000, 10, 120–131. [CrossRef] Weiskopf, N.; Mathiak, K.; Bock, S.W.; Scharnowski, F.; Veit, R.; Grodd, W.; Goebel, R.; Birbaumer, N. Principles of a brain-computer interface (BCI) based on real-time functional magnetic resonance imaging (fMRI). IEEE Trans. Biomed. Eng. 2004, 51, 966–970. [CrossRef] [PubMed] Naseer, N.; Hong, K.S. fNIRS-based BCIs: A review. Front. Hum. Neurosci. 2015, 9, 3. [CrossRef] [PubMed] Stytsenko, K.; Jablonskis, E.; Prahm, C. Evaluation of consumer EEG device Emotiv EPOC. In Proceedings of the MEi: CogSci Conference, Ljubljana, Slovenia, 17 June–18 June 2011. © 2017 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Feasibility of Equivalent Dipole Models for Electroencephalogram-Based Brain Computer Interfaces Paul H. Schimpf
ID
Department of Computer Science, Eastern Washington University, Cheney, WA 99004, USA; [email protected]; Tel.: +1-509-359-6908 Received: 5 August 2017; Accepted: 13 September 2017; Published: 15 September 2017
Abstract: This article examines the localization errors of equivalent dipolar sources inverted from the surface electroencephalogram in order to determine the feasibility of using their location as classification parameters for non-invasive brain computer interfaces. Inverse localization errors are examined for two head models: a model represented by four concentric spheres and a realistic model based on medical imagery. It is shown that the spherical model results in localization ambiguity such that a number of dipolar sources, with different azimuths and varying orientations, provide a near match to the electroencephalogram of the best equivalent source. No such ambiguity exists for the elevation of inverted sources, indicating that for spherical head models, only the elevation of inverted sources (and not the azimuth) can be expected to provide meaningful classification parameters for brain–computer interfaces. In a realistic head model, all three parameters of the inverted source location are found to be reliable, providing a more robust set of parameters. In both cases, the residual error hypersurfaces demonstrate local minima, indicating that a search for the best-matching sources should be global. Source localization error vs. signal-to-noise ratio is also demonstrated for both head models. Keywords: brain–computer interface; electroencephalogram; equivalent source model
1. Introduction The vast majority of non-invasive brain computer interfaces (BCI) based on the electroencephalogram (EEG) focus on temporal events such as sensorimotor rhythms and evoked potentials [1]. Few attempts have been made to make use of source localization other than by strategic location of the EEG recording electrodes. This is perhaps due to the fact that temporal information in the EEG is well understood, reproducible, and not overly challenging to extract in terms of the required real-time processing power. In contrast, inverting the EEG to localize brain activity is known to be an ill-posed problem [2], and is highly underdetermined when linearized with a lead-field approach [3]. While there are a number of algorithms to circumvent this problem, they tend to be compute-intensive and might represent some challenges for real-time computation in a portable system. In any event, few attempts to make use of equivalent source localization for EEG-based BCI have been published. Qin, L. et al. reported on a pilot study for BCI based on equivalent source location using off-line EEG data, wherein subjects were asked to imagine right or left hand movement based on a visual cue [4]. The time window and frequency band of the event-related desynchronization (ERD) [5] were hand-selected from the data, pre-processed via independent component analysis [6], and inverted to a source map using a regularized weighted minimum norm [7]. A single best-fit equivalent dipole was also found using the Nelder–Mead Simplex method [8]. A three concentric sphere head model was used to represent the domain [9]. All currents within 80% of the largest were summed over each hemisphere, and if the sum on the ipsilateral side was larger than the contralateral side, the source mapping method Brain Sci. 2017, 7, 118; doi:10.3390/brainsci7090118
www.mdpi.com/journal/brainsci
Brain Sci. 2017, 7, 118
2 of 11
was judged to have correctly classified the imagined hand movement. If the single best-fit equivalent dipole was located in the ipsilateral hemisphere of the brain, the equivalent dipole approach was judged to have correctly classified the imagined hand movement. Single-trial classification accuracy of 78.9% (source map) and 80.6% (equivalent dipole) were reported on a single test subject. Because the time and frequency data were hand-selected, no detection algorithm was discussed and thus no information was available regarding false positives. Nevertheless, the results indicated promise for the idea of using the localization of equivalent dipoles from EEG for BCI. No information was included on processing time, although it was stated that the approach could be implemented in an on-line BCI application with little alteration. While the noise subspace was estimated and discarded, and the source signal-to-noise ratio (SNR) was reported to be better than the scalp SNR, no numerical estimates of the SNR were reported. A follow-up study [10] expanded the analysis to two equivalent dipoles. Manual selection of the time windows was again applied and Nelder–Mead Simplex was again used for the search. The method was judged to have correctly classified the imagined hand movement if the stronger dipole was located in the ipsilateral hemisphere. An average classification rate across four subjects was 80%, a promising result but not higher than the preliminary study using a single dipole. It is important to note that the equivalent dipole locations, orientations, and amplitudes with this approach are not intended to represent an accurate model of brain activity. They instead represent an amalgam of brain activity as represented by a few dipolar sources. When a lead-field approach is used the forward problem is represented by a linear system as [3]: Ls = v
(1)
where L is the (m × 3n) lead-field, or forward transfer, matrix, s is the (3n × 1) vector of dipolar source amplitudes, with each dipolar source represented by three orthogonal components of the dipole moment, and v is the (m × 1) vector of EEG samples. In the most general case, the number of unknown source amplitudes (3n) is far greater than the number of EEG measurements (m), representing a severely underdetermined problem that is also ill-posed [2]. There are a number of algorithms that select from the infinity of solutions via various assumptions. When the number of non-zero sources is known to be, or can be restricted to be, some number less than 3n, the problem becomes underdetermined, with a single solution in the least-squares sense: −1 sq = LT LT q Lq qv
(2)
where q indicates some subset of the n candidate sources along with the corresponding columns of L. One must then search for the subset q that minimizes the residual minkLq sq − vk2 q
(3)
The question in this case is not the accuracy to which the fitted dipoles represent actual brain activity, but rather whether the location of the best fitting dipoles are useful as classification parameters in a BCI algorithm. This will depend on a subject’s ability, perhaps with training, to control the equivalent dipole locations. Because of the ill-posed nature of the problem, a reasonable precursor question is the extent to which the actual best fitting sources can be found given the head model being used. That is the question investigated here. 2. Materials and Methods 2.1. Head Models The reproducibility of a single known equivalent dipole in the face of noisy EEG measurements was investigated in both a four-shell spherical head model [11] and a realistic head model.
Brain Sci. 2017, 7, 118
3 of 11
The four-shell model has been used in several other studies (e.g., [12–15]) and is composed of concentric Brain Sci. 2017,spheres 7, 118 of varying conductivity representing the brain, cerebral-spinal fluid (CSF), 3skull, of 11 and scalp, as shown in Figure 1. The origin of the model is the center of the spheres, and the radii of Brain Sci. 2017, 7, 118 were 3 ofof 11the boundaries 7575 mm respectively. TheThe modeled conductivities the tissue tissue boundaries wereset setatat63, 63,65, 65,71, 71,and and mm respectively. modeled conductivities of layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19]. the layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19]. the tissue boundaries were set at 63, 65, 71, and 75 mm respectively. The modeled conductivities of the layers were 0.33, 1.0, 0.042, and 0.33 S/m, respectively [16–19].
Figure 1. of the spherical model. Figure 1. Cross-section Cross-section through the center of the thespherical sphericalmodel. model. Figure 1. Cross-sectionthrough throughthe the center center of
10001000 candidate source distributed pseudo-uniformly[20] [20] overthethe upper candidate sourcelocations locationswere were distributed distributed pseudo-uniformly upper 1000 candidate source locations were pseudo-uniformly [20]over over the upper hemisphere of the cortical surface, with approximate spacing of 4.8mm. mm.256 256EEG EEG electrodes were hemisphere cortical surface, withan spacing of were hemisphere of of thethe cortical surface, withan anapproximate approximate spacing4.8 of 4.8 mm. 256electrodes EEG electrodes likewise distributed over thetheupper surface.The Thelocations locationsof of candidate likewise distributed over upperhemisphere hemisphere of of the the scalp scalp surface. candidate were likewise distributed over the upper hemisphere of the scalp surface. The locations of candidate sources electrodes shownininFigure Figure2.2.The The forward forward problem source was sources andand electrodes areare shown problemfor foreach eachcandidate candidate source was sources and electrodes are shown in Figure 2. The forward problem for each candidate source was obtained from the analytic solution in [11] and used to form a 256 × 3000 lead-field, with three obtained from the analytic solution in [11] and used form a 256 × 3000 lead-field, with three obtained from dipole the analytic solution inper [11]candidate and used to form a 256 × 3000 lead-field, with three orthogonal moment directionsper source location. orthogonal dipole moment directions candidate source location. orthogonal dipole moment directions per candidate source location.
Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources in Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources blue, EEG electrodes in red). in blue, EEG electrodes in red).
Figure 2. Candidate source and EEG electrodes in spherical head model (upper hemisphere, sources TheEEG realistic head model is based on classified magnetic resonance images of a human head. The in blue, electrodes in red). tissues were classified according to the method discussed in [21]. The model contained 11 tissue types
The realistic head model is based on classified magnetic resonance images of a human head. The tissues were classified according to the method discussed in [21]. The model contained 11 tissue types
Brain Sci. 2017, 7, 118
Brain Sci. 2017, 7, 118
4 of 11
4 of 11
The realistic head model is based on classified magnetic resonance images of a human head. The tissues were classified according to the method discussed in [21]. The model contained 11 tissue with resistivities obtained from [16]. EEG observation locations were obtained from the 82 points of types with resistivities obtained from [16]. EEG observation locations were obtained from the 82 points an extended EEG 10–20 layout, visually interpolated with an additional 63 points, resulting in 145 of an extended EEG 10–20 layout, visually interpolated with an additional 63 points, resulting in EEG channels. The candidate source space consisted of 3035 locations sampled from the cortical 145 EEG channels. The candidate source space consisted of 3035 locations sampled from the cortical surfacesurface at a resolution of 4.0 × 4.0 × 3.2 resulting inina alead matrixofofdimension dimension × 9105, at a resolution of 4.0 × 4.0 × mm, 3.2 mm, resulting leadfield field matrix 145145 × 9105, with three orthogonal source directions perper location. matrixwas wascalculated calculated using with three orthogonal source directions location.The Thelead-field lead-field matrix using thethe finite element method [15,22]. This model is illustrated in Figure 3. This model, along with a finite finite element method [15,22]. This model is illustrated in Figure This model, along with a finite element for non-conforming meshes, is available fromthe theauthor author at at [23]. element solversolver for non-conforming meshes, is available from
FigureFigure 3. Realistic head model with cutaway showing internal tissues. EEG electrodes in red. 3. Realistic head model with cutaway showing internal tissues. EEG electrodes in red.
2.2. Source Localization ErrorError 2.2. Source Localization As discussed in in the the primary question is the theoretical extent to which As discussed theintroduction, introduction, the primary question is the theoretical extent to which equivalent equivalent dipoles be reliably located the EEG given the mathematically ill-posed nature of dipoles can becan reliably located from the from EEG given the mathematically ill-posed nature of the inverse problem. This question was investigated using random on the twoon head Formodels. each the inverse problem. This question was investigated usingtrials random trials themodels. two head model, ground-truth dipolesdipoles were generated at randomat source locations, with random orientations. For each model, ground-truth were generated random source locations, with random The forward of Equation was solved EEGsolved measurements normally distributed noise orientations. The model forward model of(1)Equation (1)forwas for EEGand measurements and normally was added. The best fitting to the noisy EEG wasnoisy then found by performing distributed noise was added. Thedipole best fitting dipole to the EEG was then foundanbyexhaustive performing search for the location that minimized the residual error in Equation (3). The mean and standard an exhaustive search for the location that minimized the residual error in Equation (3). The mean and deviation of the Euclidean distance between the actual and inverted dipole locations was then plotted standard deviation of the Euclidean distance between the actual and inverted dipole locations was against signal-to-noise (SNR) ratio. Surface plots of the residual error in Equation (3) over all candidate then plotted against were signal-to-noise (SNR) ratio. Surface plots of theinresidual error ininsight Equation source locations also examined for the presence of local minima, order to provide as to (3) over all candidate source locations were also examined for the presence of local minima, in order to whether a global search algorithm is required for this approach. provide insight as to whether a global search algorithm is required for this approach. 3. Results and Discussion 3.1. Spherical Head Model Figure 4 illustrates localization error vs. SNR for the spherical head model for 250 randomized
Brain Sci. 2017, 7, 118
5 of 11
3. Results and Discussion 3.1. Spherical Head Model Figure 4 illustrates localization error vs. SNR for the spherical head model for 250 randomized trials at each SNR. The possible range of localization error depends on the location of the ground truth source in the hemisphere of possibilities (see Figure 2) which was selected randomly for each5trial. Brain Sci. 2017, 7, 118 of 11 Figure 4 shows that the localization degrades dramatically for SNRs lower than 70 dB, which is an unrealistically largelarge SNR SNR for single trial EEG The compression of the standard an unrealistically for single trialrecordings. EEG recordings. The compression of the deviation standard at lower SNRs is explained by the factby that finite and inverted aresources restricted deviation at lower SNRs is explained thethe factdomain that theis domain is finite andsources inverted are to the cortical surface, which creates an upper limit to the possible localization error. As the noise restricted to the cortical surface, which creates an upper limit to the possible localization error. As the level (lower(lower SNR),SNR), the mean localization approaches an asymptote givengiven by the case noiseincreases level increases the mean localization approaches an asymptote byworst the worst scenario of a completely random guess.guess. As theAs center (mean)(mean) of the of histogram increases, the width case scenario of a completely random the center the histogram increases, the of the histogram is compressed because of the upper limit to the error. Thus, when the noise is low, width of the histogram is compressed because of the upper limit to the error. Thus, when the noise is both the mean and standard deviation are small (the localization is reliable), both grow with the noise low, both the mean and standard deviation are small (the localization is reliable), both grow with the level width the histogram begins tobegins encroach on the upper limit on error, causes the noiseuntil levelthe until the of width of the histogram to encroach on the upper limitwhich on error, which histogram to then shrink somewhat the mean as continues to approach worst casethe scenario. causes the width histogram width to then shrinkassomewhat the mean continuesthe to approach worst As aid to interpretation the trend in Figure 4, a completely random guess as to source location caseanscenario. As an aid to of interpretation of the trend in Figure 4, a completely random guess as to when the true location is at the apex (north pole) of the cortical surface has a mean of 53 mm and a source location when the true location is at the apex (north pole) of the cortical surface has a mean of standard deviation of 25 mm. When the true location is anywhere on the equator, the mean error of a 53 mm and a standard deviation of 25 mm. When the true location is anywhere on the equator, the completely guess israndom 85 mm with standard deviation of 25 mm. mean error random of a completely guessa is 85 mm with a standard deviation of 25 mm.
Figure 4. Figure 4. Localization Localization error error vs. vs. signal-to-noise signal-to-noise (SNR) (SNR) ratio ratio for for the the spherical spherical model. model. Mean Mean and and standard standard deviation are shown for for 250 250 randomized randomized trials trials at at each each SNR. SNR. Localization Localization is poor below below 70 70 dB dB (the (the deviation are shown is poor radius of the cortical sphere is 63 mm). radius of the cortical sphere is 63 mm).
A fundamental problem is that constraining the inverse to an underdetermined system does not A fundamental problem is that constraining the inverse to an underdetermined system does not avoid the ill-posed nature, and small perturbations in the measurement can produce widely different avoid the ill-posed nature, and small perturbations in the measurement can produce widely different localizations. In this case, however, the localization ambiguity occurs only for the azimuth of the localizations. In this case, however, the localization ambiguity occurs only for the azimuth of the inverted source location. This is illustrated for noise-free data in Figure 5, which plots the residual inverted source location. This is illustrated for noise-free data in Figure 5, which plots the residual from Equation (3) for all candidate sources against the azimuth (theta) and elevation (phi) of the from Equation (3) for all candidate sources against the azimuth (theta) and elevation (phi) of the source. source. All candidate sources reside at the same radius (63 mm) in this model. This shows that there All candidate sources reside at the same radius (63 mm) in this model. This shows that there are a are a number of good fits to the EEG data (i.e., leaving a small residual) that are at different azimuths, number of good fits to the EEG data (i.e., leaving a small residual) that are at different azimuths, but but only one elevation. only one elevation.
Brain Sci. 2017, 7, 118 BrainSci. Sci.2017, 2017,7,7,118 118 Brain
6 of 11 of11 11 66of
Figure 5. candidate sources for for an noise-free case. All good sources Figure 5. 5. Residual Residualat atall all candidate sources forillustrative an illustrative illustrative noise-free case. Allcandidate good candidate candidate Figure Residual at all candidate sources an noise-free case. All good (low residual) reside at the same elevation (phi), with a range of azimuths (theta) that covers the entire sources(low (lowresidual) residual)reside resideatatthe thesame sameelevation elevation(phi), (phi),with withaarange rangeof ofazimuths azimuths(theta) (theta)that thatcovers covers sources 2π radians. theentire entire2π 2πradians. radians. the
furtherillustrates illustratesthis thispoint point by plotting the elevation error vs. SNR, SNR, which performs Figure 666 further byby plotting thethe elevation error vs. vs. SNR, which performs well Figure illustrates this point plotting elevation error which performs wellover over much larger SNR range than source location. over a much larger SNRSNR range thanthan source location. well aamuch larger range source location.
(SNR)ratio Figure Figure6.6.Elevation Elevationerror errorvs. vs.signal-to-noise signal-to-noise(SNR) ratiofor forthe thespherical sphericalmodel. model.Mean Meanand andstandard standard deviation are shown for 250 randomized trials at each SNR. The elevation range for candidate sources deviation are shown for 250 randomized trials at each SNR. The elevation range for candidate sources to90 90degrees. degrees. isis00to
Figure 77 provides provides additional additional insight insight into into the the azimuth azimuth ambiguity. ambiguity. ItIt plots plots the the location location and and Figure moment of a ground truth dipole source (red) along with the locations and moments of all candidate moment of a ground truth dipole source (red) along with the locations and moments of all candidate sourcesthat thatprovide providean anEEG EEGresidual residualwithin within1% 1%(blue). (blue).The TheEEG EEGisisplotted plottedunderneath underneaththe thedipoles dipolesas as sources
Brain Sci. 2017, 7, 118
7 of 11
Figure 7 provides additional insight into the azimuth ambiguity. It plots the location and moment of a ground truth dipole source (red) along with the locations and moments of all candidate sources Brain Sci. 2017, 7, 118 7 of 11 that provide an EEG residual within 1% (blue). The EEG is plotted underneath the dipoles as an aid to understanding how thesehow different could produce same EEG on theEEG surface of the model. an aid to understanding thesedipoles different dipoles couldthe produce the same on the surface of This shows that there are a number of good candidate sources that can reproduce the surface EEG, the model. This shows that there are a number of good candidate sources that can reproduce the distributed around the entire circumference of the cortex, all differing orientations. only the surface EEG, distributed around the entire circumference ofwith the cortex, all with differingThus orientations. elevation can provide reliable information in this case. It is the spherical symmetry of the model that Thus only the elevation can provide reliable information in this case. It is the spherical symmetry of produces theambiguity azimuth of the modelthis thatambiguity producesin this inthe theinverted azimuthsources. of the inverted sources.
Figure 7. Dipolar sources that produce an Electroencephalogram (EEG) within 1% of the EEG produced Figure 7. Dipolar sources that produce Electroencephalogram (EEG) within at1% the EEG by the true source (red). The surface EEGan produced by the true source is plotted theofbottom of produced the figure. by the true source (red). The surface EEG produced by the true source is plotted at the bottom of the figure.
Taken together, Figures 5–7 indicate that only the elevations of inverted sources are reliably Taken together, Figures 5–7 indicate that only the elevations of inverted sources are reliably invertible from the EEG when a spherical head model is used. Because the azimuth of the location is invertible from the EEG when a spherical head model is used. Because the azimuth of the location is unreliable, including it in the parameter set could confound a BCI algorithm when a spherical head unreliable, including it in the parameter set could confound a BCI algorithm when a spherical head model is used. model is used. 3.2. Realistic Head Model 3.2. Realistic Head Model Figure 8 illustrates localization error vs. SNR for the realistic head model. It performs much better Figure 8 illustrates localization error vs. SNR for the realistic head model. It performs much over a range of SNRs than the spherical head model, with reasonable results down to 10 dB. The use of better over a range of SNRs than the spherical head model, with reasonable results down to 10 dB. a realistic head model thus provides at least three reasonable BCI parameters per inverted source: x, y, The use of a realistic head model thus provides at least three reasonable BCI parameters per inverted and z, or equivalently, r, theta, and phi. source: x, y, and z, or equivalently, r, theta, and phi.
Brain Sci. 2017, 7, 118
8 of 11
Brain Sci. 2017, 7, 118
8 of 11
Figure 8.Figure Localization errorerror vs. vs. signal-to-noise ratioforfor realistic Mean and 8. Localization signal-to-noise (SNR) (SNR) ratio thethe realistic head head model.model. Mean and standard deviation are shown for 250 randomized trials at each SNR. Localization is much better than standard deviation are shown for 250 randomized trials at each SNR. Localization is much better than the spherical head model at lower SNRs. the spherical head model at lower SNRs. 3.3. Local vs. Global Minima
3.3. Local vs. Global Minima
For the spherical head model, Figure 5 also illustrates local minima in the residual error surface
For that thewould spherical model, also illustrates local minima in model, the residual error trap ahead localized searchFigure for the 5 minimum residual. In the realistic head the origin for surface that would trap asources localized searchin for the minimum residual. In the realistic head the origin candidate was placed the Cerebellum in order to compare approximately to the model, coordinate system of the spherical model. A difference is that candidate sources in the realistic model do not for candidate sources was placed in the Cerebellum in order to compare approximately to the all fall at theof same The residual hypersurface thuscandidate has three dimensions and cannot be model coordinate system theradius. spherical model.error A difference is that sources in the realistic shown in its entirety on a single plot. Figure 9 shows the projection of a residual hypersurface onto the do not all fall at the same radius. The residual error hypersurface thus has three dimensions and azimuth (theta) and radial directions, indicating the presence of local minima for the realistic head cannot be shown its aentirety on a for single plot. Figure shows the projection of a residual model as well.inThus localized search the best-fitting dipoles9for either model cannot be expected hypersurface onto the (theta) and radial indicating the presence of local minima to be reliable, andazimuth a global search strategy should directions, be used. A simple a global search is to exhaustively compute the residual for all possibilities, for the realistic head approach model astowell. Thus a localized search for the best-fitting dipoles for either model p ), where n is the number of candidate sources (one-third which has a computational complexity of O(n cannot be expected to be reliable, and a global search strategy should be used. the number columns in the lead-field) and p is the number of sources allowed in the solution. This can A simple approach to a global search is to exhaustively compute the residual for all possibilities, become quickly intractable as p increases. On a Core i3-4360 CPU (3.7 GHz, Intel Inc., Santa Clara, CA, p which has a computational complexity of O(n whereLA, n USA) is therunning number of candidate sources (oneUSA) with 16 GB of DDR3 RAM (Crucial Inc.,),Gretna, Octave 4.2.1 ([24]) under third theMicrosoft numberWindows columns in the lead-field) and p is the number of sources 7 (Microsoft Corp., Redmond, WA USA), this search requiredallowed less than in 0.1the s forsolution. = 1 on thequickly realistic intractable head model. as A value of p = 2 required 202 s. The pseudo-inverses for different This canpbecome p increases. On a Core i3-4360 CPU (3.7 GHz, Intel Inc., Santa column subsets of the lead-field are trivially parallelizable. Thus p = 2 is within the range of real-time Clara, CA, USA) with 16 GB of DDR3 RAM (Crucial Inc., Gretna, LA, USA) running Octave 4.2.1 performance if a speedup on the order of 200 can be achieved. This is currently under investigation ([24]) under Microsoft Windows 7 (Microsoft Corp., Redmond, WA USA), this search required less using GPU parallelization. than 0.1 s for p = 1 on the realistic head model. A value of p = 2 required 202 s. The pseudo-inverses for different column subsets of the lead-field are trivially parallelizable. Thus p = 2 is within the range of real-time performance if a speedup on the order of 200 can be achieved. This is currently under investigation using GPU parallelization.
Brain Sci. 2017, 7, 118 Brain Sci. 2017, 7, 118
9 of 11 9 of 11
Figure 9. Residual at all candidate sources for an illustrative noise-free case in the realistic head model. Figure 9. Residual at all candidate sources for an illustrative noise-free case in the realistic head model.
3.4. Summary and Applicability 3.4. Summary and Applicability The results above show that the following should be taken into account when inverting equivalent The results above show that the following should be taken into account when inverting dipole locations from the EEG: equivalent dipole locations from the EEG: (a) When a spherical head model is used, the inverted azimuth of the dipole is unreliable in the face (a) When a spherical head model is used, the inverted azimuth of the dipole is unreliable in the face of even modest amounts of noise, although the elevation of the dipole is reliable. of even modest amounts of noise, although the elevation of the dipole is reliable. (b) All components components of of the the location location are are reasonably (b) All reasonably reliable reliable when when aa realistic realistic head head model model is is used. used. (c) The residual cost function for locating the dipole exhibits local minima, and thus a global (c) The residual cost function for locating the dipole exhibits local minima, and thus a global search search strategy should should be be used. used. strategy
For an an EEG-based EEG-based BCI, BCI, the the performance performance cannot cannot therefore therefore be be expected expected to to be be optimal optimal ifif aa spherical spherical For head model model is is used used with with the the full full inverted inverted dipole dipole location, location,or orwhen whenusing usingaalocalized localizedsearch searchstrategy. strategy. head It should be noted that it is not necessarily the case that a realistic head model for BCI be It should be noted that it is not necessarily the case that a realistic head model for subjectBCI be specific. Given the somewhat limited success of using a spherical model [4,10], along with subject-specific. Given the somewhat limited success of using a spherical model [4,10], along these with findings on location ambiguity, there is good to expect that any realistic model, model, such assuch [25], these findings on location ambiguity, there is reason good reason to expect that any realistic would improve BCI performance without needing to be subject-specific. as [25], would improve BCI performance without needing to be subject-specific. The findings findings here here are are specific specific to to source source localization localization from from the the EEG. EEG. Similar Similar concerns concerns may may or or may may The not exist for source localization from the magnetoencephalogram (MEG), which is related to the EEG not exist for source localization from the magnetoencephalogram (MEG), which is related to the EEG problem via via the the Biot–Savart Biot–Savart law law [15]: [15]: problem
Z J rˆ µ00 J × rˆ dv B B = 4 r22 dv 4π r
(4) (4)
Brain Sci. 2017, 7, 118
10 of 11
The use of a vector measurement (B) and the fact that J represents both volume and source currents may mitigate the source localization ambiguities found with the EEG, which is something not investigated here. It is worth noting that other modalities provide more direct opportunities to characterize the location of brain activity. As an imaging technology, functional magnetic resonance imaging (fMRI) [26] provides spatial information directly. Functional near-infrared spectroscopy (fNIRS) [27] provides some opportunities for obtaining spatial information through the strategic placement of emitter–detector pairs. The fNIRS technology appears promising and may well prove to supplant EEG as a common BCI approach when inexpensive fNIRS headsets become available as is the case for EEG (e.g., [28]). 4. Conclusions This work shows that equivalent dipole approaches to EEG-based BCI produce ambiguous results in source azimuth when a spherical head model is used, and can thus be expected to provide reliable information only in terms of source elevation. It also shows that a realistic head model does not suffer such localization ambiguity, providing a more robust set of parameters for potential BCI classification. It further shows that the EEG residual hypersurface demonstrates local minima, indicating that a search for best-fitting dipolar sources should be global. Acknowledgments: The cost for open access publication of this work was provided by Eastern Washington University. Conflicts of Interest: The author declares no conflict of interest.
References 1.
2. 3. 4. 5. 6. 7. 8. 9. 10.
11. 12. 13.
Sellers, E.W.; Arbel, Y.; Donchin, E.; Pfurtscheller, G.; McFarland, D.J.; Allison, B.Z.; Faller, J.; Neuper, C. Existing BCIs. In Brain-Computer Interfaces, Principles and Practices; Wolpaw, J.R., Wolpaw, E.W., Eds.; Oxford University Press: New York, NY, USA, 2012; pp. 215–249, ISBN 978-0-19-538885-5. Kabanikhin, S.I. Definitions and examples of inverse and ill-posed problems. J. Inverse Ill-Posed Probl. 2008, 16, 317–357. [CrossRef] Weinstein, D.; Zhukov, L.; Johnson, C. Lead-field bases for electroencephalography source imaging. Ann. Biomed. Eng. 2000, 28, 1059–1065. [CrossRef] [PubMed] Qin, L.; Ding, L.; He, B. Motor Imagery Classification by Means of Source Analysis for Brain Computer Interface Applications. J. Neural Eng. 2005, 2, 65–72. [CrossRef] [PubMed] Pfurtscheller, G. Event-Related Desynchronization; Elsevier: New York, NY, USA, 1999; pp. 303–325. Hyvärinen, A.; Oja, E. A Fast Fixed-Point Algorithm for Independent Component Analysis. Neural Comput. 1997, 9, 1483–1492. [CrossRef] Hansen, P.C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 1992, 34, 561–580. [CrossRef] Kowalik, J.; Osborne, M.R. Methods for Unconstrained Optimization Problems; Elsevier: New York, NY, USA, 1968. Rush, S.; Driscoll, D. Current distribution in the brain from surface electrodes. Anesth. Analg. 1968, 47, 717–723. [CrossRef] [PubMed] Kamousi, B.; Liu, Z.; He, B. Classification of motor imagery tasks for brain-computer interface applications by means of two equivalent dipoles analysis. IEEE Trans. Neural Sys. Rehabil. Eng. 2005, 13, 166–171. [CrossRef] [PubMed] Zhou, H. Computation of the potential distribution in a four-layer anisotropic concentric spherical volume conductor. IEEE Trans. Biomed. Eng. 1992, 39, 154–158. [CrossRef] [PubMed] Hoesk, R.S., Jr.; Sances, A.; Larson, S.J. The contributions of intracerebral currents of the EEG and evoked potentials. IEEE Trans. Biomed. Eng. 1978, 25, 405–413. [CrossRef] [PubMed] Stok, C.J. The influence of model parameters on EEG/MEG single dipole source estimation. IEEE Trans. Biomed. Eng. 1987, 34, 289–296. [CrossRef] [PubMed]
Brain Sci. 2017, 7, 118
14. 15. 16. 17. 18. 19.
20. 21.
22. 23. 24. 25.
26.
27. 28.
11 of 11
Schimpf, P.H.; Ramon, C.; Haueisen, J. Dipole models for the EEG and MEG. IEEE Trans. Biomed. Eng. 2002, 49, 409–418. [CrossRef] [PubMed] Schimpf, P.H. Application of Quasi-static Magnetic Reciprocity to Finite Element Models of the MEG Lead-Field. IEEE Trans. Biomed. Eng. 2007, 54, 2082–2088. [CrossRef] [PubMed] Geddes, A.; Baker, L.E. The specific resistance of biological material—A compendium of data for the biomedical engineer and physiologist. Med. Biol. Eng. 1967, 5, 271–293. [CrossRef] [PubMed] Rush, S.; Driscoll, D. EEG-electrode sensitivity—An application of reciprocity. IEEE Trans. Biomed. Eng. 1969, 16, 15–22. [CrossRef] [PubMed] Oostendorp, T.F.; Delbecke, J.; Stegman, D.F. The conductivity of the human skull: Results in vivo and in vitro measurements. IEEE Trans. Biomed. Eng. 2000, 47, 1487–1492. [CrossRef] [PubMed] Hoekema, R.; Huiskamp, G.J.M.; Wieneke, G.H.; Leijten, F.S.S.; van Veelen, C.W.M.; van Rijen, P.C.; van Huffelen, A.C. Measurement of the conductivity of skull, temporarily removed during epilepsy surgery. Brain Topogr. 2004, 16, 29–38. [CrossRef] Deserno, M. How to Generate Equidistributed Points on the Surface of a Sphere. Available online: http://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf (accessed on 22 May 2017). Shrinidhi, N.; Haynor, D.R.; Wang, Y.; Joregenson, D.B.; Bardy, G.H.; Kim, Y. An efficient tissue classifier for building patient-specific finite element models from X-ray CT images. IEEE Trans. Biomed. Eng. 1996, 43, 333–337. [CrossRef] [PubMed] Schimpf, P.; Haueisen, J.; Ramon, C.; Nowak, H. Realistic computer modeling of electric and magnetic fields of human head and torso. Parallel Comput. 1998, 24, 1433–1460. [CrossRef] GalerWin. Available online: https://sites.google.com/site/pschimpf99/home/software/galerwin (accessed on 14 September 2017). GNU Octave. Available online: https://www.gnu.org/software/octave/ (accessed on 14 September 2017). Lancaster, J.L.; Woldorff, M.G.; Parsons, L.M.; Liotti, M.; Freitas, C.S.; Rainey, L.; Kochunov, P.V.; Nickerson, D.; Mikiten, S.A.; Fox, P.T. Automated Talairach Atlas labels for functional brain mapping. Hum. Brain Mapp. 2000, 10, 120–131. [CrossRef] Weiskopf, N.; Mathiak, K.; Bock, S.W.; Scharnowski, F.; Veit, R.; Grodd, W.; Goebel, R.; Birbaumer, N. Principles of a brain-computer interface (BCI) based on real-time functional magnetic resonance imaging (fMRI). IEEE Trans. Biomed. Eng. 2004, 51, 966–970. [CrossRef] [PubMed] Naseer, N.; Hong, K.S. fNIRS-based BCIs: A review. Front. Hum. Neurosci. 2015, 9, 3. [CrossRef] [PubMed] Stytsenko, K.; Jablonskis, E.; Prahm, C. Evaluation of consumer EEG device Emotiv EPOC. In Proceedings of the MEi: CogSci Conference, Ljubljana, Slovenia, 17 June–18 June 2011. © 2017 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
