Order 1106099: migraine

F

E a

b

U c

#

a

A R R A A

K S M S D M

1

m a [ o o a s l i m t

(

h 1

Biomedical Signal Processing and Control 18 (2015) 263–273

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / b s p c

requency-spatial beamformer for MEG source localization

lizabeth A. Thompson a,∗, Jing Xiang b, Yingying Wang c,1

Electrical Engineering, Purdue University, Fort Wayne (IPFW), 2101 E. Coliseum Blvd., Fort Wayne, IN 46805, United States MEG Research, MEG Center, Division of Neurology, Cincinnati Children’s Hospital Medical Center, 3333 Burnet Avenue, Cincinnati, OH 45229-3039, nited States Division of Developmental Medicine, Laboratories of Cognitive Neuroscience, Gaab Lab, Boston Children’s Hospital, 1 Autumn Street, 6th Floor, Room 640, Boston, MA 02215, United States

r t i c l e i n f o

rticle history: eceived 9 September 2014 eceived in revised form 8 December 2014 ccepted 20 January 2015 vailable online 20 February 2015

eywords: tatistical signal processing EG source localization

patial location ata driven

a b s t r a c t

This article introduces a unique and powerful new method for spatial localization of neuronal sources that exploit the high temporal resolution of magnetoencephalography (MEG) data to locate the originating sources within the brain. A traditional frequency beamforming algorithm was adapted from its conven- tional application to yield information on the spatial location of simulated neuronal signals. The concept is similar to that used in signal source localization in magnetic resonance imaging (MRI) in which spatial location is determined by the frequency of oscillation of the MR signal. Whereas a traditional frequency beamformer uses the time course values of all sensors in the dataset to assign a power value for each possible frequency in the signal, it provides no information on the spatial location of those frequencies. Our approach assigns a power value to each location in the three-dimensional head volume. To compute this power value, the time courses of a subset of sensors closest to that location in space are used rather

odel-free than all the time courses in the dataset. Our novel technique incorporates actual MEG sensor locations of the closest sensors at each location in space. The approach is relatively simple to implement, yields good spatial resolution, and accurately spatially locates a simulated source in low signal-to-noise envi- ronments. In this work, its performance is compared to that of the synthetic aperture magnetometry (SAM) beamformer and shown to exhibit improved spatial resolution.

© 2015 Elsevier Ltd. All rights reserved.

. Introduction

Magnetoencephalography (MEG) measures the extracranial agnetic field generated by the brain’s electrical activity in an

ttempt to obtain information on the underlying neural activity 1]. It is a brain imaging technology of unsurpassed temporal res- lution which produces high resolution (order of milliseconds) ne-dimensional time course data at, typically, 275 locations round the surface of the head. The locations of these 275 sen- ors are shown in Fig. 1. Utilizing this data to deduce the spatial ocation of underlying neuronal sources has proven to be challeng-

ng. In Fig. 1, the source of brain activation is depicted as the red

ass inside the sensors. The group of circles around the perime- er represent the MEG sensor array of the CTF MEG system. The

∗ Corresponding author. Tel.: +1 260 481 6361. E-mail addresses: [email protected] (E.A. Thompson), [email protected]

J. Xiang), [email protected] (Y. Wang). URL: http://www.engr.ipfw.edu/ thompson/ (E.A. Thompson).

1 Tel.: +1 857 218 5479.

ttp://dx.doi.org/10.1016/j.bspc.2015.01.004 746-8094/© 2015 Elsevier Ltd. All rights reserved.

colors of these sensors indicate the signal strength projected from the source. The red and yellow sensors indicate strong detectable signals, while gray indicates a weak or no signal detected at that location.

Recently, a novel MEG source localization technique was intro- duced that adapts a traditional frequency beamforming algorithm to yield information on the spatial location of simulated MEG sources [2]. A traditional frequency beamformer computes the frequencies present in a signal. It does so by computing a power value for each possible frequency in the signal and plotting this power value as a function of frequency. In computing this power value, the traditional frequency beamformere uses the time course values of all sensors in the dataset. Thus, it yields information on the frequencies present in the signal but provides no information on the spatial location of those frequencies. Our novel localization approach adapts this method by assigning a frequency power value to each location in a three-dimensional (3D) head volume.

That is, at each location in the 3D head volume, the time course values of only the closest 21 sensors are extracted and used to compute a power value for a given frequency of interest. In this way, our approach provides a means of quantifying the strength

264 E.A. Thompson et al. / Biomedical Signal Proc

Fig. 1. Locations of 275 MEG sensors in the CTF MEG system. The colors of these sensors indicate the signal strength projected from the source. The red and yellow sensors indicate strong detectable signals, while gray indicates a weak or no signal d l

o l c a s m d

1

p t c B m a

c T m r c f i i E

1

b t [ r E m [

strained minimum variance (LCMV) [49] beamformer and multiple

etected at that location. (For interpretation of the references to color in this figure egend, the reader is referred to the web version of this article.)

f the frequency content of the signal as a function of its spatial ocation. In this research, we more fully develop the approach and ompare its performance to that of the widely recognized synthetic perture magnetometry (SAM) beamformer in accurately locating imulated MEG sources. In addition, in this approach we revise the ethod of inserting simulated activations to the traditional current

ipole.

.1. Background

Whereas functional magnetic resonance imaging (fMRI) [3,4], rovides unparalleled spatial resolution, yielding detailed informa- ion on the anatomical structures of the human brain, it relies on hanges in cerebral blood flow and blood volume, referred to as the OLD (blood-oxygen-level-dependent) effect, which is an indirect easure of brain functionality and limits its temporal resolution to

bout one second [5]. On the other hand, MEG and its electrical counterpart, electroen-

ephalography (EEG), record direct responses from neuron activity. hus, MEG and EEG are capable of tracking brain dynamics with illisecond temporal resolution [5]. This unsurpassed temporal

esolution has led to the widespread use of EEG and MEG by clini- ians and scientists investigating physiologic and pathologic brain unction [6]. While EEG and MEG are superior to all other functional maging methods with regard to temporal resolution, spatial local- zation is impaired by the ill-posedness and non-uniqueness of the EG/MEG inverse problem [5].

.1.1. Sensor time course data techniques To extract information from MEG data regarding interactions

etween spatially disparate brain regions, some current popular echniques rely on the analysis of the sensor time course data 7–10]. Standard methods include coherence analysis and event-

elated synchronization/desynchronization (ERS/ERD). ERD and RS have been demonstrated using a variety of different sensory- otor and cognitive paradigms and at a variety of frequencies

11,12]. Activity has been reported in the alpha [13–16] (8–15 Hz),

essing and Control 18 (2015) 263–273

beta [17–20] (15–30 Hz), gamma [21–23] (30–50 Hz), and theta [24–26] (6–10 Hz) bands. Several researchers [8] have used time course data in conjunction with wavelet-based techniques.

However, these techniques do not fully exploit the capabil- ity of MEG to localize sources of neuronal activity. Thus, many researchers are attempting to use the MEG time course data to deduce the origin of the underlying sources of brain activ- ity. Task-related changes in spectral power in specific frequency bands are localized using conventional dipole source modeling [27–31], sometimes used in conjunction with contour maps, and beamformers [32–37], to place brain activity in the context of anatomy.

1.1.2. Dipole source models Postsynaptic activity is often modeled as a single current dipole,

which represents an idealized point source [38]. In this case to esti- mate the neural sources, the forward problem must be solved which maps a dipole source of known location, strength, and orientation, to an array of MEG sensors [39].

Given a forward model, that is, an assumed shape and conductiv- ity of the head, and a set of MEG sensor time course measurements, the inverse problem is to find the current density distribution that generated the observed measurements. The solution to the inverse problem is an estimate of the number, location, and time series of each of the dipole sources. Solving the inverse problem involves making inferences regarding the distribution of underlying neu- ronal sources based on the observed MEG values at the surface, and there is no unique result [40].

Since the human evoked response is generated by multiple simultaneously active brain sources producing fields that overlap both in space and time, modeling this activity with a single current dipole can result in erroneous conclusions concerning its neu- ronal generation [41]. Thus, multiple current dipoles have also been applied to spatiotemporal analysis [40]. Dipole fitting [42,43] is a variant used to solve the inverse problem. Imaging methods using the Bayesian perspective have also been used to solve the inverse problem [44–46]. Nonparametric distributed source modeling is an alternative, non-dipole approach to deal with the indeterminacy of the inverse problem [47]. Magnetic Field Tomography (MFT) uses continuous probabilistic estimates of the primary current density to construct tomographic displays of brain activity. It solves the inverse problem using arbitrarily distributed non-dipole sources [48].

1.1.3. Beamformers A beamformer is a processor used in conjunction with an array

of sensors to provide a form of spatial filtering [49]. Spatial filter- ing refers to discrimination of signals on the basis of their spatial location. This concept is analogous to the more familiar temporal filtering, where one discriminates between signals based on their temporal frequency content [37]. The beamformer sensor array is designed to receive signals radiating from a specific location and attenuate those from other locations. The MEG beamformers developed to date rely on a lead field matrix which contains the computed MEG sensor values (forward solution) that would result for all possible locations and orientations of potential sources of neuronal activity. Furthermore, most of the MEG beamforming lit- erature uses a dipole as the source model in calculating this lead field matrix.

The most common MEG beamformers are the linearly con-

signal classification (MUSIC) [35,42,50] and variations thereof. A limitation of the beamforming approach using dipole source mod- els is that they are highly dependent on the underlying assumptions forming the basis of the lead field matrix.

l Proc

2

2

R

w t b m [ v (

X

m d

� R

o i o

t X m fi r i

2

o b b

� P

w

X

i o f t

E.A. Thompson et al. / Biomedical Signa

. Theory

.1. Correlation matrix

The correlation matrix, Rx, of a random vector x is defined by

x = E{xxH } (1)

here H denotes Hermitian transpose and E denotes the expecta- ion operator [51]. When the correlation matrix is unknown, it can e estimated using various methods. The covariance method of esti- ating the correlation matrix is unbiased and positive semidefinite

51] and is composed as follows. Suppose there are N data points in ector x: x(0), x(1),. . .x(N − 1). To estimate a P × P correlation matrix P ≤ N), a data matrix, X, of size (N − P + 1) × P is defined as in Eq. (2):

=

⎡ ⎢⎢⎢⎣

x(P − 1) x(P − 2) · · · x(0) x(P) x(P − 1) · · · x(1)

... ...

...

x(N − 1) x(N − 2) · · · x(N − P)

⎤ ⎥⎥⎥⎦ (2)

The estimate of the correlation matrix, � Rx, using the covariance

ethod is defined by Eq. (3), where X is as defined in Eq. (2) and H enotes Hermitian transpose:

x = 1

N − P + 1 (X) H X (3)

Sometimes the scaling factor 1/(N − P + 1) is omitted [51]. When P = N, the data matrix X reduces to the vector x in reverse

rder, and the scaling factor reduces to 1. When there are M real- zations of vector x, all consisting of N data points (each perhaps riginating from a different sensor), for example,

x0(0), x0(1), · · ·x0(N − 1) x1(0), x1(1), · · ·x1(N − 1) x2(0), x2(1), · · ·x2(N − 1)

· · · xM−1(0), xM−1(1), · · ·xM−1(N − 1)

hese are stacked together to form the rows of an M × N data matrix , which is subsequently used in calculating the covariance esti- ate of the correlation matrix per Eq. (3). This is equivalent to

nding the covariance estimate of the correlation matrix for each ealization and then summing them all together to yield a compos- te matrix.

.2. Traditional frequency beamformer

Spectrum estimation is used to determine the frequency content f a random signal. The periodogram,

� Px(ω), is a standard, Fourier

ased technique often used for this purpose. It is defined in Eq. (4) elow [51]:

x(ω) = 1 N

∣∣X(ω)∣∣2 (4) here

(ω) = N−1∑ n=0

x(n)e−jωn (5)

s the Discrete Time Fourier transform of the data sequence, x(n), f length N, and ω is the digital frequency in radians. The digital requency is related to the analog frequency, ̋ = 2�f, by the rela- ionship ω = ˝T, where T is the sampling period, in seconds. An

essing and Control 18 (2015) 263–273 265

equivalent, alternate, representation of Eq. (5) is depicted in Eq. (6) below.

X(ω) = wH x (6) where x is a column vector [x(0), x(1),. . .x(N − 1)]T, w is the column vector of Eq. (7) containing the Fourier coefficients:

w = [e−j0ω , e−j1ω , e−j2ω , · · ·e−j(N−1)ω ]T (7) and H denotes Hermitian transpose. Using this alternate represen- tation, the Periodogram of Eq. (4) can also be expressed in alternate form [51]:

� Px(ω) =

1 N

wH xxH w (8)

and its expected value is

E{P̂x(ω)} = 1 N

wH Rxw (9)

where Rx is as defined in Eq. (1). Eq. (9) can be used to implement a traditional frequency beam-

former, with the estimate of the correlation matrix, � Rx, defined per

Eq. (3), used in place of Rx. For each frequency within the range of all possible frequencies, a unique weight vector is defined per Eq. (7). The range of possible frequencies is governed by the Nyquist criteria. Thus, the maximum analog frequency value, fmax, is half the sampling frequency, and the frequency resolution is fs/N where fs is the sampling frequency in samples/sec or Hertz, and N is the number of samples in the data sequence. The analog frequencies are converted to their digital counterparts using the relation ω = 2�fT, which is subsequently inserted into Eq. (7) to yield a unique weight vector for each frequency. Then, using Eq. (9), a power value is com- puted corresponding to each frequency, and the resulting power values are plotted as a function of frequency. The peak value(s) of this spectrum will occur at the frequency or frequencies present in the signal.

For our research, this traditional beamformer was modified to provide a novel technique yielding information on the spatial loca- tion of simulated neuronal sources.

2.3. Adaptation of conventional frequency beamformer

In this research, the traditional frequency beamformer of Eq. (9) was adapted to yield information on the spatial location of sim- ulated neuronal signals. If all the MEG time course sensor data is stacked together to form the rows of a data matrix, X, and subse- quently used to calculate the covariance estimate of the correlation matrix per Eq. (3), subsequent application of Eq. (9) for a range of frequencies would yield information on the frequency content of the MEG data. The result is the traditional frequency beamformer.

To utilize this information in spatial source localization, a three- dimensional (3D) volume in the CTF MEG head frame of reference is defined. For each of the spatial locations in this grid, the locations of the 21 closest MEG sensors are determined. Then, a subset of the full time course data is extracted, consisting only of the time course data from the 21 sensors closest to that spatial location. These 21 time courses are stacked together to form the rows of a data matrix, X, and subsequently used to calculate the covariance estimate of the correlation matrix per Eq. (3). Define this estimate of the correlation matrix as

� Rsubset . Using this correlation matrix in

the calculation of Eq. (9) in conjunction with a weight vector of Eq. (7) defined for the frequency of interest results in a power value that is a function of spatial location, per Eq. (10):

P(x, y, z) = 1 N

wH R̂subset w (10)

The power value resulting from Eq. (10) is subsequently plotted as a function of the spatial coordinates (x, y, z). The process must

2 l Proc

b ( c f r s

2

s

F i l

66 E.A. Thompson et al. / Biomedical Signa

e repeated for each frequency of interest. The power value of Eq. 10) provides a means of quantifying the strength of the frequency ontent of a signal as a function of its spatial location. This concept or signal source localization is similar to that used in magnetic esonance imaging (MRI) in which the frequency content of the ignal enables its localization in space.

.4. Rationale for use of 21 closest sensors

Having selected a search grid in 3D space, the number of closest ensors to each spatial location of this grid was tabulated. In this

ig. 2. Frequency-spatial beamformer output SNR = −2 dB. The site of the single simulate ndicated by the color bar, the red regions identify the highest output values of the algo ocation of the dipole. (For interpretation of the references to color in this figure legend, t

essing and Control 18 (2015) 263–273

research, the 3D search grid consisted of 142,814 spatial locations. None of these 142,814 locations shared the same 20 closest sen- sors, in the same order. That is, using this grid, the sequence of 20 closest sensors, sorted by distance, was unique for each location in the grid. It was subsequently experimentally determined that 21 closest sensors provided optimum results.

3. Methods

A simulated MEG signal was modeled as a dipole with a sinusoidal time course of 10 Hz frequency and 3.54 × 10−8

d dipole at (0, 0, .07) m is identified with an ‘X’ in the frame labeled z = 0.07 m. As rithm. The algorithm attains its maximum output value in close proximity to the he reader is referred to the web version of this article.)

E.A. Thompson et al. / Biomedical Signal Processing and Control 18 (2015) 263–273 267

Fig. 3. Frequency-spatial beamformer output SNR = −2 dB. The sites of the three correlated dipole sources at (−.002, 0, .07), (0, 0, .07), (.002, 0, .07) m are identified with ‘X’ in t mpose o close c

A b h m m s b p

S

he frame labeled z = 0.07 m. Due to their close proximity, the X’s are partially superi utput values of the algorithm. The algorithm attains its maximum output value in olor in this figure legend, the reader is referred to the web version of this article.)

mpere-meters amplitude. The dipole source was assumed to e originating at (x, y, z) coordinates (0, 0, .07) m in the ead frame of reference. The forward problem was solved to ap the dipole source to the array of MEG sensors. Zero ean Gaussian noise was then added to create a desired

ignal-to-noise ratio (SNR). SNR was defined per Eq. (11) elow, where PS is the signal power, and PN is the noise ower.

NR = 10log10 (

PS PN

) (11)

d on one another. As indicated by the color bar, the red regions identify the highest proximity to the locations of the dipoles. (For interpretation of the references to

The signal power of a sinusoid, PS, was defined per Eq. (12) below, where A is its peak amplitude. The noise power, PN, was defined as the variance of the random noise, �2.

PS = A2

2 (12)

A 3D volume in the head frame of reference was defined in which x and y values ranged from −0.1 to 0.1 m at increments of

.002 m, and z ranged from .018 to 0.07 m, in increments of .004 m. For each of these locations in 3D space, the locations of the 21 clos- est sensors were determined. The simulated time courses for these 21 closest sensors were stacked together to form a data matrix, X,

268 E.A. Thompson et al. / Biomedical Signal Processing and Control 18 (2015) 263–273

Table 1 Comparison of frequency-spatial beamformer results to those of SAM, 99% threshold, single dipole of amplitude 35.4 nAm and frequency 10 Hz at (0, 0, .07) m.

Trial SNR (dB) Frequency-spatial beamformer SAM

No. true positives No. false positives

Max distance from true location (mm)

No. true positives

No. false positives

Max distance from true location (mm)

1 10 0 2 5.7 0 1 8.2 2 10 0 1 5.7 0 1 11.3 3 10 0 2 5.7 0 1 8 4 10 0 1 4.5 0 1 17.9 5 10 0 1 4.5 0 1 17.2

1 2 0 1 4.5 0 1 10 2 2 0 1 5.7 0 1 8.2 3 2 0 1 5.7 0 1 8.9 4 2 0 1 4.5 0 1 17.9 5 2 0 1 4.5 0 1 8

1 −1 0 1 4.5 0 1 17.2 2 −1 0 2 5.7 0 1 13.4 3 −1 0 2 4.5 0 1 25.1 4 −1 0 2 5.7 0 1 13.4 5 −1 0 2 5.7 0 1 10.2 1 −2 0 1 4.5 0 1 4 2 −2 0 1 5.7 0 1 8.2 3 −2 0 1 5.7 0 1 23.3 4 −2 0 1 4.5 0 1 17.9 5 −2 0 1 4.5 0 1 17.2 1 −3 0 2 4.5 0 1 4 2 −3 0 2 5.7 0 1 17.2 3 −3 0 1 4.5 0 1 17.9 4 −3 0 1 4.5 0 1 18.9 5 −3 0 1 4.5 0 1 14.4 1 −4 0 1 4.5 0 1 2 2 −4 0 2 4.5 0 1 17.9 3 −4 0 1 4.5 0 1 8.5 4 −4 0 1 4.5 0 1 18.9 5 −4 0 1 5.7 0 1 12.6 1 −5 0 1 4.5 0 1 4 2 −5 0 1 2.8 0 1 8 3 −5 0 2 5.7 0 1 6.3 4 −5 0 1 4.5 0 3 20.6 5 −5 0 1 5.7 0 1 23.3 1 −10 0 1 4.5 0 1 11.3 2 −10 0 1 5.7 0 2 12.8

o n p f � R

T C a

3 −10 0 1 5.7 4 −10 0 1 5.7 5 −10 0 1 5.7

f size M × N, where M is the number of sensors (21), and N is the

umber of time points. For computational efficiency, only 200 time oints were utilized, which was adequate to provide the desired requency resolution. At each location in space, a correlation matrix, subset , was then computed per Eq. (3). The analog frequency value,

able 2 omparison of frequency-spatial beamformer results to those of SAM, 95% threshold, th mplitude 35.4 nAm and frequency 10 Hz.

Trial SNR (dB) Frequency-spatial beamformer

No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distance to one of dipoles (mm)

1 +5 0 5 5.7 2 2 +5 0 2 4.5 2.8 3 +5 0 4 5.7 2.8 4 +5 0 4 5.7 2.8 5 +5 0 5 5.7 2

1 −2 0 1 4.5 4 2 −2 0 4 5.7 2 3 −2 0 2 4.5 2 4 −2 0 2 4.5 4 5 −2 0 2 5.7 4

0 2 12.8 0 1 10.8 0 1 6.3

f, of 10 Hz was converted to its digital equivalent, ω, per Eq. (13)

below

ω = 2�f fs

(13)

ree correlated dipole sources at (−.002, 0, .07), (0, 0, .07), (.002, 0, .07) m, each of

SAM

No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distance to one of dipoles (mm)

0 1 10 8 0 2 23.3 20.6 0 1 14.1 12.8 0 2 16.1 11.7 0 1 10 8.9

0 1 12.8 11.7 0 1 10 8.9 0 1 17.2 16.1 0 1 8.2 8 0 1 15.6 14.4

E.A. Thompson et al. / Biomedical Signal Processing and Control 18 (2015) 263–273 269

Table 3 Comparison of Frequency-Spatial Beamformer results to those of SAM, 95% threshold, three uncorrelated dipole sources at (−.002, 0, .07), (0, 0, .07), (.002, 0, .07) m, each of amplitude 35.4 nAm and frequency 10 Hz.

Trial SNR (dB) Frequency-spatial beamformer SAM

No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distance to one of dipoles (mm)

No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distance to one of dipoles (mm)

1 +5 0 4 5.7 2.8 0 1 14.4 13.4 2 +5 0 4 5.7 2.8 0 1 7.2 5.7 3 +5 0 4 5.7 2.8 0 1 8.9 8.2 4 +5 0 4 5.7 2.8 0 1 17.9 17.1 5 +5 0 3 5.7 4 0 1 10 8.9

1 −2 0 4 5.7 2.8 0 1 7.2 6.3 2 −2 0 2 4.5 2.8 0 2 54.6 20.6 3 −2 0 3 5.7 2.8 0 1 8.9 8.2

u v u a

p d a

1 i + f n a

t ( a x r t t 1 w E r o a

T C e

4 −2 0 5 5.7 2 5 −2 0 1 4.5 4

sing a sampling frequency, fs, of 1200 samples/s. The resulting alue of ω was used in composing the weight vector of Eq. (7) for se in Eq. (10). For each location in space, a power value was then ssigned per Eq. (10).

Using Eq. (10), power values over the entire 3D space were com- uted and then the resultant array was scaled to 256 levels for isplay. The maximum value of the array was expected to occur t the location of the simulated activation.

Having simulated a signal consisting of a dipole source with a 0 Hz sinusoidal time course, with added noise, the beamform-

ng algorithm of Eq. (10) was applied, while varying SNR from 10 to −10 dB. Actual MEG sensor array locations were obtained rom test data of a real human subject scanned in a 275 chan- el MEG system and these locations were incorporated into the lgorithm.

Performance of the algorithm in locating the simulated activa- ion was compared to that of the synthetic aperture magnetometry SAM) beamformer [35]. When implementing the SAM algorithm,

3D volume in the head frame of reference was defined in which values ranged from −0.1 to 0.1 m at increments of .002 m, and z anged from .02 to 0.07 m, in increments of .002 m. The length of the ime course data for the SAM simulations consisted of the entire 30 rials in the dataset, each consisting of 4096 samples, for a total of 22,880 samples. For each SNR value from +10 to −10 dB, five trials ere run using different realizations of zero mean Gaussian noise.

ach of these realizations was obtained by setting the seed of the andom number generator to a specific value, enabling execution f the trial using identical noise realizations for both the proposed lgorithm and the SAM beamformer. In each trial, the number of

able 4 omparison of frequency-spatial beamformer results to those of SAM, 95% threshold, th ach of amplitude 35.4 nAm and frequency 10 Hz.

Trial SNR (dB) Frequency-spatial beamformer

No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distanc one of dipoles (mm)

1 +5 0 4 5.7 2.8 2 +5 0 4 5.7 2.8 3 +5 0 4 5.7 2 4 +5 0 4 5.7 2.8 5 +5 0 5 5.7 2

1 −2 0 4 5.7 2.8 2 −2 0 3 5.7 4 3 −2 0 3 5.7 4 4 −2 0 2 5.7 4 5 −2 0 4 5.7 2.8

0 1 10.8 10.2 0 1 6 4

true and false activations was recorded as well as the maximum distance of the predicted from the actual dipole location. Thresh- olds of 90%, 95%, and 99% of the maximum value of the algorithm output were used in this analysis.

Subsequently, a cluster of three dipoles were inserted in the dataset. These were assumed to be originating at (x, y, z) coordi- nates (−.002, 0, .07) m, (0, 0, .07) m, and (.002, 0, .07) m in the head frame of reference. A sinusoidal time course of 10 Hz frequency and 3.54 × 10−8 Ampere-meters amplitude was assumed for all three of these dipoles. For this test, the time courses of all three dipoles were correlated, that is, their phases were identical. Both algorithms were again executed using SNR values of +5 and −2 dB and five different realizations of random Gaussian noise at each SNR value. The number of true and false positives located by each of the two algorithms was recorded as well as the maximum dis- tance from a predicted dipole location to the center dipole located at (0, 0, .07) m. In addition, the closest distance from the predicted dipole locations to any of the three true dipole locations was noted. Next, the same cluster of three dipoles was modified to make them uncorrelated. The phases of the sinusoidal time courses of the three dipoles were modified to be 0◦, 45◦, and 90◦, respectively, while keeping their amplitudes and frequencies unaltered. The experi- ment was repeated using the same testing scenario as that of the correlated sources. A partially correlated scenario was also created. For this, the two outermost dipoles were correlated with each other,

but the center dipole was composed of a 200 sample inactive period followed by a 200 sample active sinusoidal period. Since the samp- ling rate was 1200 samples/s, 200 samples represented 0.1667 s. For this test, 400 time points were utilized for the Frequency-Spatial

ree partially correlated dipole sources at (−.002, 0, .07), (0, 0, .07), (.002, 0, .07) m,

SAM

e to No. true positives

No. false positives

Max distance from center dipole (mm)

Closest distance to one of dipoles (mm)

0 1 13.4 12.6 0 1 10 8.2 0 1 7.2 6.3 0 2 29.5 24.2 0 1 6 6

0 1 10.2 8.2 0 1 14.4 13.4 0 1 14.4 13.4 0 2 29.5 24.2 0 2 23.3 8.5

270 E.A. Thompson et al. / Biomedical Signal Processing and Control 18 (2015) 263–273

Fig. 4. Frequency-spatial beamformer output SNR = −2 dB. The site of the single simulated dipole at (0, 0, .046) m is identified with an ‘X’ in the frame labeled z = 0.046 m. As indicated by the color bar, the red regions identify the highest output values of the algorithm. In this case, the maximum output is not localized as a pinpoint. However, the p For in t

B a

t t

4

d p T

lots of Fig. 5 can be used to reveal the location of the activation for this situation. ( he web version of this article.)

eamformer; 30 trials (122,880 data points) were used in the SAM nalysis.

Finally, the simulated dipole was located at (0, 0, .046) m, and he ability of the two algorithms to correctly locate the dipole in his new location was recorded.

. Results

A typical output for the frequency-spatial beamformer, single ipole scenario, is displayed in Fig. 2. The algorithm provides out- ut for each of the slices in the z-dimension that were evaluated. he location of the single dipole at (0, 0, .07) m is identified with

terpretation of the references to color in this figure legend, the reader is referred to

an ‘X’ in the z = .07 m frame of Fig. 2. As indicated in the asso- ciated colorbar, voxels of highest output are identified in red in Fig. 2; those with lowest output values are blue. Using a thresh- old of 95% of the maximum algorithm output, there are three false positives within 5.7 mm of the true dipole location at (0, 0, .07) m. Using a threshold of 99%, there is a single false positive at (−.004, −.004, .07) m, which is within 5.7 mm of the true dipole location.

A typical output for the frequency-spatial beamformer, three dipole scenario, is displayed in Fig. 3. The locations of the three dipoles at (−.002, 0, .07) m, (0, 0, .07) m, and (.002, 0, .07) m are identified with ‘X’ in the z = .07 m frame of Fig. 3. Due to their close

l Proc

p a t c A r c

f w w t a 9 s d

s A f f f o

F z

E.A. Thompson et al. / Biomedical Signa

roximity, the X’s are partially superimposed on one another. Using threshold of 95% of the maximum algorithm output, there are wo false positives within a maximum distance of 5.7 mm from the enter dipole and within 4 mm of the dipole at (−.002, 0, .07) m. lthough Fig. 3 was obtained using correlated dipole sources, the esults were also typical of those using uncorrelated and partially orrelated sources.

A comparison of results for the novel frequency-spatial beam- ormer algorithm to those of SAM is provided in Tables 1–4. Table 1 as created from simulations using a single dipole at (0, 0, .07) m hile varying the SNR and the individual noise realizations. While

hresholds of 90%, 95%, and 99% of the maximum value of the lgorithm output were used in the analysis, only the results of 9% threshold are displayed in Table 1. On average, the frequency- patial beamformer exhibits greater accuracy in locating the single ipole.

Tables 2–4 provide comparisons of results of the frequency- patial beamformer to those of SAM using three dipole sources. ll were obtained using a 95% threshold. Table 2 displays results

or correlated dipole sources, Table 3 for uncorrelated, and Table 4 or partially correlated dipole sources. The frequency-spatial beam- ormer was consistently more accurate in identifying the locations f the sources than those of SAM for a variety of SNR and

ig. 5. Mesh plots of selected frames of Fig. 4: (a) .07 m; (b) .066 m; (c) .046 m; (d) .042 −0.046 m as the depth of the dipole source.

essing and Control 18 (2015) 263–273 271

correlation scenarios. It was also more accurate in distinguishing between distinct sources in close proximity.

Fig. 4 displays the frequency-spatial beamformer output when the single dipole is located at (0, 0, .046) m. In this case, the output must be displayed in a different fashion to determine the location of the dipole. Fig. 5 displays mesh plots for selected frames of Fig. 4. To conserve space, only the z = .07 m, .066 m, .046 m, and .042 m frames are displayed. The disappearance of the local maxima below z = 0.046 m correctly identifies z = 0.046 m as the depth of the dipole source.

5. Discussion

The Fourier basis of the frequency-spatial beamformer pro- vides an explanation as to why performance of the algorithm is not adversely affected by the correlation condition of the sources. Fourier analysis provides a robust measure of the frequency content of a signal, regardless of its phase. In addition, the frequency- spatial beamformer is not influenced by the orientation of the

dipole sources and should work equally well if the sources are assumed to be in a form other than dipoles, as long as the time course is at least partially periodic. Furthermore, the new algorithm should yield results consistent to those reported here for a wide

m. The disappearance of the local maxima below z = 0.046 m correctly identifies

2 l Proc

r i

i s t s c e t b b 2 t p a

f T t t p f

s p t t S m o n e s i F d a i

a c u o o p s

6

t t n b a s c p t i s p r a

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

72 E.A. Thompson et al. / Biomedical Signa

ange of frequency values. Indeed, the preliminary study reported n Thompson et al. [2] assumed a source frequency of 200 Hz.

The spatial resolution of the frequency-spatial beamformer s based on utilizing information from a large quantity of sen- ors at varying distances from the source. In its present form, he algorithm does not explicitly make a distinction between a ensor close to a source and one farther away. All 21 sensors losest to a given location in space are treated equally How- ver, if fewer sensors were available, as in, for example, EEG, he algorithm could be adapted to include a weighting factor ased on the distance of the sensor from the source. That is, efore computing

� Rsubset of Eq. (10), the data from each of the

1 closest sensors could be weighted in proportion to its dis- ance from the specified location in 3D space. This technique could erhaps improve spatial resolution when fewer sensors are avail- ble.

The frequency-spatial beamformer presented here is suited or detecting activations that are periodic or partially periodic. his is due to the use of the weight vector of Eq. (7) containing he Fourier coefficients. For nonperiodic or spontaneous signals, his weight vector could be modified to another form, for exam- le, wavelet-based, to enable detection of signals of a variety of orms.

In this study, the use of known activations in the form of uperimposed simulations enabled an objective and accurate com- arison to the widely accepted method of SAM. Having established he frequency-spatial beamformer’s merit in accurately locating hese simulated activations and demonstrated its potential over AM, this work represents the next stage in the process of deter- ining its viability for analysis of real MEG data. Implementation

n real data will necessitate a thorough investigation of the sig- ificance of indications where discrepancies between techniques xist. It will also entail some trial and error in the investigation tage to discern the optimum length of the time course to include n the computations and to learn how to interpret the results. uture work will involve testing the algorithm on real human ata containing known activation locations, and confirming its bilty to accurately locate them, in a variety of activation scenar- os.

It is anticipated that the frequency-spatial beamformer will have pplication in a wide variety of research and clinical settings. It ould be a powerful tool in functional connectivity analysis. To se it for this purpose, the power values of Eq. (10) would be btained at various snapshots in time to establish a time-sequence f activations. In clinical settings it could potentially be used for reoperative functional brain mapping and localization of epileptic eizures.

. Conclusions

The concept for signal source localization used in MRI in which he frequency of oscillation of the MR signal enables its localiza- ion in space is applied here to provide a unique and powerful ew method of MEG source localization. A traditional frequency eamforming algorithm has been adapted from its conventional pplication to yield information on the spatial location of neuronal ignals. It is based on the premise that the strength of the frequency omponents of a signal will be strongest at the sensors in closest roximity to the signal, and is not influenced by the orientation of he dipole source. The algorithm performs well on simulated data n low signal-to-noise environments and provides relatively good

patial resolution. This data-driven approach is comparatively sim- le to implement and requires no signal model when applied to eal data. It relies solely on the actual sensor time course values of

subset of sensors closest to a given spatial location.

[

[

essing and Control 18 (2015) 263–273

References

[1] E. Ducla-Soares, Modeling in magnetoencephalography, in: S. Sato (Ed.), Advances in Neurology, Magnetoencephalography, vol. 54, Raven Press, New York, 1990, p. 95.

[2] E. Thompson, S.K. Holland, J. Xiang, Y. Wang, MEG source localization using a frequency beamformer, in: Proceedings of the 2010 IEEE 36th Annual Northeast Bioengineering Conference (NEBEC), 27–28 March, New York, NY, 2010, pp. 1–2.

[3] J. Belliveau, et al., Functional mapping of the human visual cortex by magnetic resonance imaging, Science 254 (1991) 716–719.

[4] S. Ogawa, et al., Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging, Proc. Natl. Acad. Sci. U.S.A. 89 (1992) 5951–5955.

[5] R. Konig, C. Sieluzycki, P. Durks, Tiny signals from the human brain: acquisition and processing of biomagnetic fields in magnetoencephalography, J. Low Temp. Phys. 146 (5–6) (2007) 697–718.

[6] M.S. Hu, G. Stead, Worrell, Automatic identification and removal of scalp ref- erence signal for intracranial EEGs based on independent component analysis, IEEE Trans. Biomed. Eng. 54 (9) (2007) 1560–1572.

[7] A. Ukil, Practical denoising of MEG data using wavelet transform, in: Proceedings of 13th International Conferance on Neural Information Processing, vol. 4233, 2006, pp. 578–585.

[8] K. Jerbi, J. Lachaux, S. Baillet, L. Garnero, Imaging cortical oscillations during sustained visuomotor coordination in MEG, in: IEEE ISBI’04, 2004, pp. 380–383.

[9] G. Marrelec, et al., Conditional correlation as a measure of mediated inter- activity in fMRI and MEG/EEG, IEEE Trans. Signal Proc. 53 (9) (2005) 3503–3516.

10] J. de Munck, et al., A maximumlikelihood estimator for trial-to-trial variations in noisy MEG/EEG data sets, IEEE Trans. Biomed. Eng. 51 (12) (2004) 2123–2128.

11] K. Singh, et al., Task-related changes in cortical synchronization are spatially coincident with the hemodynamic response, Neuroimage 16 (2002) 103–114.

12] E. Basar, et al., Gamma, alpha, delta, and theta oscillations govern cognitive processes, Int. J. Psychophysiol. 39 (2–3) (2001) 241–248.

13] G. Pfurtscheller, Graphical display and statistical evaluation of event-related desynchronization (ERD), Electroencephalogr. Clin. Neurophys. 43 (5) (1977) 757–760.

14] E. Basar, et al., Alpha oscillations in brain functioning: an integrative theory, Int. J. Psychophysiol. 26 (1–3) (1997) 5–29.

15] W. Klimesch, et al., Simultaneous desynchronization and synchronization of different alpha responses in the human electroencephalograph: a neglected paradox, Neurosci. Lett. 284 (1–2) (2000) 97–100.

16] C.M. Krauss, et al., Test-retest consistency of the event-related desynchronization/event-related synchronization of the 4–6, 8–10, and 10–12 Hz frequency bands during a memory task, Clin. Neurophysiol. 112 (5) (2001) 750–757.

17] G. Pfurtscheller, A. Stancak, G. Edlinger, On the existence of different types of central beta rhythms below 30 Hz, Electroencephalogr. Clin. Neurophysiol. 102 (4) (1997) 316–325.

18] C. Andrew, G. Pfurtscheller, Lack of bilateral coherence of post-movement cen- tral beta oscillations in the human electroencephalogram, Neurosci. Lett. 273 (2) (1999) 89–92.

19] M. van Burik, G. Pfurtscheller, Functional imaging of postmovement beta event- related synchronization, J. Clin. Neurophysiol. 16 (4) (1999) 383–390.

20] M. Taniguchi, et al., Movement-related desynchronization of the cerebral cor- tex studied with spatially filtered magnetoencephalography, Neuroimage 12 (3) (2000) 298–306.

21] U. Ribary, Magnetic field tomography of coherent thalamocortical 40-Hz oscil- lations in humans, Proc. Natl. Acad. Sci. U. S. A. 88 (240) (1991) 11037–11041.

22] M. Joliot, U. Ribary, R. Llinas, Human oscillatory brain activity near 40 Hz coex- ists with cognitive temporal binding, Proc. Natl. Acad. Sci. U. S. A. 91 (24) (1994) 11748–11751.

23] C. Tallon-Baudry, et al., Stimulus specificity of phase-locked and non-phase- locked 40 Hz visual responses in human, J. Neurosci. 16 (13) (1996) 4240–4349.

24] C. Tesche, J. Karhu, Theta oscillations index human hippocampal activation during a working memory task, Proc. Natl. Acad. Sci. U. S. A. 97 (2) (2000) 919–924.

25] A. Burgess, J. Gruzalier, Short duration power changes in the EEG during recog- nition memory for words and faces, Psychophysiology 37 (5) (2000) 596–606.

26] W. Klimesch, et al., Episodic retrieval is reflected by a process specific increase in human electroencephalographic theta activity, Neurosci. Lett. 302 (1) (2001) 49–52.

27] R. Grasman, H. Huizenga, L. Waldorp, K. Bocker, P. Molenaar, Frequency domain simultaneous source and source coherence estimation with an application to MEG, IEEE Trans. Biomed. Eng. 51 (1) (2004) 45–55.

28] K. Iramina, et al., Measurement of somatosensory evoked response using func- tional MR images and MEG, IEEE Trans. Magn. 33 (5) (1997) 4260–4262.

29] A. Theiben, et al., A first estimate for EEG/MEG dipole reconstructions, in: IEEE 18th Conference on EMBS, 1996, pp. 810–811.

30] J. Mosher, P. Lewis, R. Leahy, Multiple dipole modeling and localization from spatio-temporal MEG data, IEEE Trans. Biomed. Eng. 39 (6) (1992) 541–557.

31] T. Knosche, Transformation of whole-head MEG recordings between different

sensor positions, Biomed. Eng. 47 (3) (2002) 59–62.

32] A. Herdman, et al., Determination of activation areas in the human auditory cortex by means of synthetic aperture magnetometry, Neuroimage 20 (2003) 995–1005.

l Proc

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

[

E.A. Thompson et al. / Biomedical Signa

33] A. Hadjipapas, et al., Assessing interactions of linear and nonlinear neuronal sources using MEG beamformers: a proof of concept, Clin. Neurophys. 116 (6) (2005) 1300–1313.

34] K. Sekihara, et al., Application of an MEG eigenspace beamformer to reconstruc- ting spatio-temporal activities of neural sources, Hum. Brain Mapp. 15 (2002) 199–215.

35] J. Vrba, S. Robinson, Linearly constrained minimum variance beamformers, syn- thetic aperture magnetometry, and MUSIC in MEG applications, in: Proceedings of the 34th Asilomar Conference, vol. 1, 2000, pp. 313–317.

36] W. Gaetz, D. Cheyne, Localization of human somatosensory cortex using spa- tially filtered magnetoencephalography, Neurosci. Lett. 340 (2003) 161–164.

37] B. Van Veen, W. van Drongelen, M. Yuchtman, A. Suzuki, Localization of brain electrical activity via linearly constrained minimum variance spatial filtering, IEEE Trans. Biomed. Eng. 44 (9) (1997) 867–880.

38] V. Pizzella, G. Romani, Principles of magnetoencephalography, in: S. Sato (Ed.), Advances in Neurology, Magnetoencephalography, vol. 54, Raven Press, New York, 1990, pp. 1–6.

39] F. Darvas, D. Pantazis, E. Yildirim, R. Leahy, Mapping human brain func- tion with MEG and EEG: methods and validation, Neuroimage 23 (2004) 5289–5299.

40] M. Balish, R. Muratore, The inverse problem in electroencephalography and magnetoencephalography, in: S. Sato (Ed.), Advances in Neurology, Magne- toencephalography, vol. 54, Raven Press, New York, 1990, pp. 79–88.

41] C. Baumgartner, Clinical Electrophysiology of the Somatosensory Cortex, Springer-Verlag, New York, 1993, pp. 18–19.

[

essing and Control 18 (2015) 263–273 273

42] J. Mosher, R. Leahy, Recursive MUSIC: a framework for EEG and MEG source localization, IEEE Trans. Biomed. Eng. 45 (11) (1998) 1342–1355.

43] M. Scherg, Fundamentals of dipole source potential analysis, Adv. Audiol. 6 (1990) 40–69.

44] S. Baillet, I. Garnero, A Bayesian approach to introducing anatomo-functional priors in the EEG/MEG inverse problem, IEEE Trans. Biomed. Eng. 44 (5) (1997) 374–385.

45] G. Adde, M. Clerc, R. Keriven, Imaging methods for MEG/EEG inverse problem, Int. J. Bioelectromagn. 7 (2) (2005) 111–114.

46] R. Pascual-Marqui, C. Michel, D. Lehman, Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain, Int. J. Psychophysiol. 18 (1994) 49–65.

47] P. Lewis, J. Mosher, R. Leahy, Neuromagnetic source reconstruction, in: Proceedings of the IEEE 20th International Conference on Acoustics, Speech, and Signal Processing (ICASSP), Detroit, MI, 1995, pp. 2911–2914.

48] A.A. Ioannides, J.P.R. Bolton, C.J.S. Clarke, Continuous probabilistic solutions to the biomagnetic inverse problem, Inverse Probl. 6 (1990) 523–542.

49] B. Van Veen, K. Buckley, Beamforming: a versatile approach to spatial filtering, IEEE ASSP Mag. 5 (2) (1988) 4–24.

50] T. Knosche, B. Maess, A. Nakaamura, A. Friederici, Human communication

investigated with magnetoencephalography: speech, music, and gestures, in: H. Preissl (Ed.), Magnetoencephalography, Elsevier Academic Press, Amster- dam, 2005, pp. 79–120.

51] C.W. Therrien, Discrete Random Signals and Statistical Signal Processing, Pren- tice Hall, Englewood Cliffs, NJ, 1992, pp. 307–310, 586–596, 610.