760

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

Reconstructing Spatio-Temporal Activities of
Neural Sources Using an MEG Vector Beamformer
TechniqueKensuke Sekihara*, Member, IEEE, Srikantan S. Nagarajan, David Poeppel, Alec Marantz, and Yasushi Miyashita

Abstract—We have developed a method suitable for reconstructing spatio-temporal activities of neural sources by using
magnetoencephalogram (MEG) data. The method extends theadaptive beamformer technique originally proposed by Borgiotti
and Kaplan to incorporate the vector beamformer formulation inwhich a set of three weight vectors are used to detect the source
activity in three orthogonal directions. The weight vectors of thevector-extended version of the Borgiotti–Kaplan beamformer
are then projected onto the signal subspace of the measurement
covariance matrix to obtain the final form of the proposed
beamformer’s weight vectors. Our numerical experiments showthat both spatial resolution and output signal-to-noise ratio of the
proposed beamformer are significantly higher than those of theminimum-variance-based vector beamformer used in previous
investigations. We also applied the proposed beamformer to two
sets of auditory-evoked MEG data, and the results clearly demonstrated the method’s capability of reconstructing spatio-temporal
activities of neural sources.
Index Terms—Beamformer, biomagnetism, functional neuroimaging, magnetoencephalography, MEG inverse problems,
neuromagnetic signal processing.

I. INTRODUCTIONAMONG the various kinds of functional neuroimaging
methodologies, the major advantage of magnetoencephalography (MEG) is its ability to provide fine time
resolution of the millisecond order [1]. Neuromagnetic imaging
can thus be used to visualize neural activities with such a finetime resolution, and to provide functional information aboutbrain dynamics [2]. Toward this goal, a number of algorithms
for reconstructing spatio-temporal source activities have beendeveloped. Well-known approaches for this reconstruction
employ the model of the equivalent current dipole (ECD) [3],which assumes a highly localized source. Although this ECD
Manuscript received March 17, 2000; revised April 6, 2001. The work of
S. S. Nagarajan was supported by a grant from the Whitaker Foundation. This
work was carried out as part of the MIT-JST International Cooperative Research
Project “Mind Articulation.” Asterisk indicates corresponding author.
*K. Sekihara is with the Department of Electronic Systems and Engineering,Tokyo Metropolitan Institute of Technology, Asahigaoka 6-6, Hino, Tokyo 1910065, Japan (e-mail: ksekiha@cc.tmit.ac.jp).
S. S. Nagarajan is with the Department of Bioengineering, University of Utah,
Salt Lake City, UT 84112-9202 USA.
D. Poeppel is with the Department of Linguistics and Biology, University of
Maryland, College Park, MD 20742 USA.
A. Marantz is with the Department of Linguistics and Philosophy, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
Y. Miyashita is with the Department of Physiology, The University of Tokyo,
School of Medicine, Hongo, Tokyo 113-0033, Japan.
Publisher Item Identifier S 0018-9294(01)05134-5.

model has successfully been applied to neuromagnetic data, wecannot rely on ECD modeling when the source is distributedor when no information on the spatial extent of the source isavailable.
Another approach for spatio-temporal reconstruction is basedon the linear estimation method [4]; it assumes voxels in the reconstruction region and attempts to estimate the moment of asource assigned to each voxel by using least-squares fitting. Although this approach does not impose any models on the neuromagnetic source, a naive form of such an approach has a serious
problem in that the estimation becomes severely ill posed. Thisis because the number of voxels generally reaches, at least, a
few thousand, so several thousand parameters need to be estimated from the measured data obtained at only one to two hundred points on the scalp surface. To reduce the influence fromthis ill-posed condition, an efficient method of constraining theleast-squares solution should be developed. This has been an
area of active research, and many kinds of investigations in thisdirection have been reported [5]–[7].
In this paper, we explore the possibility of applying a class oftechniques called the adaptive beamformer to this reconstruction problem. The adaptive beamformer provides a versatileform of spatial filtering suitable for processing data from anarray of sensors. Adaptive-beamformer-type techniques wereoriginally developed in the fields of array signal processing, including radar, sonar, and seismic exploration [8], and they havebeen already applied to the MEG/EEG source-reconstruction
problem [9]–[12]. In these investigations, the minimum-variance beamformer, which is one of most popular adaptive
beamformer techniques, was modified to incorporate the detection of three-dimensional (3-D) vector sources. Particularly,
in [9], [10], [12], a vector beamformer technique has beendeveloped on the basis of the minimum-variance beamformer;
the vector beamformer uses a set of three weight vectors for
detecting the source activity in three orthogonal directions
such as , , and , thereby reconstructing not only the source
magnitude but also the source orientation.
This paper develops a vector beamformer technique on thebasis of the Borgiotti–Kaplan beamformer [13]. The developed
beamformer performs significantly better than the minimumvariance beamformer used in the previous investigations, withrespect to the spatial resolution and the output signal-to-noise
ratio (SNR). In Section II, after a brief introduction to the minimum-variance-based beamformer technique, we formulate our
proposed vector beamformer. In Section III, a series of numerical experiments verify the effectiveness of the proposed beam-

0018–9294/01$10.00 © 2001 IEEE

SEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES

former. In Section IV, we apply the proposed beamformer totwo sets of auditory MEG data. The results of these applications
demonstrate its capability of reconstructing spatio-temporal activities of neural sources. Throughout this paper, plain italics indicate scalars, lower-case boldface italics indicate vectors, andupper-case boldface italics indicate matrices.
II. METHOD
A. Definitions and Problem Formulation
Let us define the magnetic field measured by the thas, and a column vector
detector coil at time
as a set of measured
is the total number of detector coils and thedata where
superscript indicates the matrix transpose. A spatial locationis represented by a 3-D vector :.
current sources are assumed to generate the
A total ofneuromagnetic field, and the locations of these sources are. The moment magnitude of the thdenoted as, and the source magnitude
source at time is defined as.
vector is defined as
To express the orientation of the th source, we define the an,gles between its moment vector and the , , and axes as, and, respectively. The orientation of the th source, whereis defined as a vector
and is equal to, or throughout thismatrix that expresses the orientations
paper. We define aof all sources as such that
..
.
..
.

..

.

We assume in this paper that the orientation of each source istime independent.
The lead field vector for the component of a source atis defined as. Here,expresses the th sensor output induced by the unit-magnitude
source that is located at and directed in the direction. We,
define the lead field matrix as
which represents the sensitivity of the sensor array at . The
lead field vector representing the sensitivity of a sensor array, which is calculated
in the direction at is denoted as. The composite lead field matrix for thefromentire set of sources is defined as

761

magnitude of a source located at and directed in thetion using the following linear spatial filter operation:(3)is the estimated moment magnitude. In (3), thewhererepresents a set of weights that charactercolumn vector
izes the property of the beamformer. It should be pointed out that
according to (3), the power of the output noise due to the addiis proportional totive white noise
where
is assumed to be the white Gaussian noise andis
its variance. That is, the power of the output noise is proportional, which is called the white noise gain for this reason.
toB. Existing Beamformer Techniques for Reconstructing Neural
Source Activities
1) Minimum-Variance Distortionless Beamformer: Onewell-known beamformer technique is the minimum-variancedistortionless beamformer [14], in which the weight vector
is obtained by minimizing
with the constraint of
1. (Although the weight vector depends on the
pointing location , we omit the explicit notation of unless
this omission causes ambiguity.) This weight vector is(4)This minimum-variance beamformer is widely used in various
signal-processing fields. One problem arises when we apply it
to reconstructing neural sources. That is, to calculate the beamformer weight by using (4), we must first determine the source
orientation at each . This determination is not straightforward, although a method have been proposed for this purpose[15].
2) Vector-ExtensionofMinimum-VarianceBeamformer: Instead of estimating the source orientation andmagnitude separately, the minimum-variance beamformer canbe modified to estimate not only the source magnitude but also
the source orientation. Such a beamformer simultaneously estimates the source moment in three orthogonal directions such
as , , and . Let us denote the unit vectors in the , , anddirections, as , , and , respectively, i.e.,,, and. Let us also denote the
,
, and
as
,
weight vectors that estimate, and. Then, these weight vectors can be derived by thefollowing minimizations with multiple constraints:(1)
The relationship between

andis then expressed as(2)

is the additive noise. We define, for later use, the
where
such
covariance matrix of the measured magnetic field as
where
indicatesthat
the ensemble average.
To estimate the source moment from the measured magnetic
field, we focus on the class of techniques referred to as a beamformer [8]. The beamformer technique estimates the moment

direc-

subject toandsubject to
and

762

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

subject to
and

(5)

The minimum-variance beamformer with multiple linear constraints, referred to as the linearly constrained minimum-variance beamformer, is known to have the following solution [9],
[10]:
(6)
.
where is defined as
Equation (5) indicates that, when estimating one of thethree orthogonal components of the source moment, weneed to suppress the other two components. This is because
when estimating one of the three components, the other two
components behave like perfectly correlated virtual sources.
Therefore, without this suppression, considerable amountof signal cancellation should arise. By applying these null
constraints, however, we can avoid this signal cancellation, andthe beamformer can detect the source moment projected inthree orthogonal directions. Such a beamformer is referred toas a vector beamformer.
There are two problems when applying (6) to actual
MEG/EEG source reconstruction problems. First, the beamformer output has erroneously large values near the center ofthe sphere used for the forward calculation. This is because
becomes very small when approaches the center of-dependent artifacts, the use
the sphere. To avoid thesehas beenof the normalized lead field matrix
suggested [10], [12].
Second, the performance of the beamformer in (6) is very
sensitive to errors in calculating the lead field matrix, when applied to the spatio-temporal reconstruction of the source activities [15]. This problem is known to be partly solved by replacingin (6) with its regularized inverse[11], [12],[16], [17] where the parameter is the regularization parameter. Such replacement, however, is known to degrade the spatial resolution, providing a tradeoff between the SNR and thespatial resolution [11]. In the following section, we propose abeamformer that is free from this tradeoff and can attain an inherently higher spatial resolution.

Fig. 1. The coordinate system and the source-detector configuration used in
the numerical experiments. The coordinate origin was set at the center of thedetector coil located at the center of the coil array. The cross section at x 1.0
is shown. The circle shows the cross section of the sphere used for the forward
calculation. The square shows the reconstruction region for the experiments
whose results are shown in Figs. 3–6.

=

TABLE I
SOURCE PARAMETER VALUES USED FOR THE NUMERICAL EXPERIMENTS
IN SECTION III

1) Vector-Type Borgiotti–Kaplan Beamformer: The vectorextended Borgiotti–Kaplan beamformer is obtained by using thefollowing constrained minimizations:subject toandsubject toandC. Formulation of the Proposed Method
The proposed beamformer is formulated on the basis of the
beamformer developed by Borgiotti and Kaplan [13]; its weightwith the constraint ofvector is derived by minimizing
1. Because the Borgiotti–Kaplan beamformer has theunit white noise gain, the output power of the Borgiotti–Kaplan
beamformer is equal to the power of the signal normalized
by the power of the noise [13]. The weight vectors of the
proposed beamformer are derived by a two-step procedure.
The first step extends the Borgiotti–Kaplan beamformer to thevector-type beamformer. The second step further extends thisvector-extended version of the Borgiotti–Kaplan beamformer
to an eigenspace-projected beamformer.

subject to
and

(7)

. Let us introduce a scalar
We first derive the expression forwhere can be determined
constant such that
1. Then, the constrained optifrom the relationshipmization problem in (7) is changed tosubject to(8)

SEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES

763

(a)

(b)
Fig. 2. (a) Time courses of the three sources assumed in the numerical experiments. Time courses from the first to the third sources are shown from the top to
the bottom, respectively. The time courses are shown for the time window between 0 and 400 ms. Each time course is normalized by its maximum value. Threevertical broken lines indicate the time instants 220, 268, and 300 ms at which the source-moment magnitude is displayed in Figs. 3(a)–6(a). (b) The generated
magnetic field used for the numerical experiments.

The solution of this optimization problem is known to havethe formThus, we getfrom the relationship
1. Using exactly the same derivation, the weights
andcan be derived, and a set of the weights is expressed as

(9)
(11)
Then, we have
(10)
where

It can be shown that the above beamformer retains the property of the Borgiotti–Kaplan beamformer, and its output poweris equal to the power of the source activity normalized by thepower of the output noise due to the additive sensor noise.
2) Extension to an Eigenspace-Projection Beamformer: The extension to an eigenspace-projection beamformer is attained by projecting the weight vectors in (11) onto

764IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

the signal subspace of the measurement covariance matrix.
The eigenspace projection improves the output SNR without
sacrificing the spatial resolution. The general analysis regarding
how this eigenspace projection improves the output SNR has
been reported [18].
Unless the source activities are perfectly correlated with eachhas eigenvalues greater thanandeigenother,values equal to where is the variance of the additive noise.
Assuming that the eigenvalues are numbered in decreasingas,order, let us define the matrix
withare the eigenvectors of .
whereis the maximum-likelihood estimate
The column span of[19]. The weight vectors of theof the signal subspace of
eigenspace-projected beamformer is obtained by using [20]
(a)

(12)
The projection onto the signal subspace using the above equation, however, invalidates the null constraints imposed on the orthogonal components. This can be understood by considering,. The null constraints in this casefor example, the case of
0 and
0. However, let us conshould be
sider

(13)
and
are not necessarily in the signal
Becauseandsubspace, we generally have
and, therefore,and, leading to the relationships
and. Consequently, we conclude that the signal
does not preserve the null consubspace projector
straints.
It can, however, be shown that the eigenspace-projection beamformer in (12) can detect the three orthogonalcomponents of the source moment even though the null
constraints are not preserved. Omitting the time notation,
into two
let us decompose the measured magnetic field
where
is the magnetic fieldparts,
is the congenerated from the target source at , andcantribution from other sources. The magnetic field,
be expressed as
expresses the orientation of the target
where
is its moment magnitude. Then, the estimatedsource, andcomponent of the source moment,, is expressed as. Since thedoes not pass the signal other than that from ,weight0 and, consequently
we have

Because the vector
subspace, we get(b)Fig. 3. (a) Results of the spatio-temporal reconstruction obtained using
the minimum-variance-based vector beamformer in (6). The upper-left,
upper-right, and lower-left maps respectively show the snapshots of thesource-moment magnitude at 220, 268, and 300 ms. The lower-right map
shows the time-averaged reconstruction. (b) Estimated time courses from thefirst to the third sources are shown from the top to the bottom, respectively. The
three vertical broken lines indicate the time instants of 220, 268, and 300 ms.

Therefore, we finally get

(14)

(16)

is in the signal

and
Similarly, we can also obtain. Thus, the eigenspace-projection
beamformer in (12) can detect the three orthogonal componentsof the source moment, even though the null constraints are notpreserved.

(15)

SEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES

765

(a)

(a)

(b)
Fig. 4. (a) Results of the spatio-temporal reconstruction obtained using theminimum-variance-based vector beamformer in (6) together with the use of theregularized inverse (R + I ) . The parameter was set at 0:003 , where is the largest eigenvalue of R . The upper-left, upper-right, and lower-leftmaps respectively show the snapshots of the source-moment magnitude at 220,
268, and 300 ms. The lower-right map shows the time-averaged reconstruction.
(b) Estimated time courses from the first to the third sources are shown from the
top to the bottom, respectively. The three vertical broken lines indicate the time
instants of 220, 268, and 300 ms.

Fig. 5. (a) Results of the spatio-temporal reconstruction with the weightvectors obtained using (11) alone. The upper-left, upper-right, and lower-leftmaps respectively show the snapshots of the source-moment magnitude at
220 ms, 268 ms, and 300 ms. The lower-right map shows the time-averaged
reconstruction. (b) Estimated time courses from the first to the third sources
are shown from the top to the bottom, respectively. The three vertical brokenlines indicate the time instants of 220, 268, and 300 ms.

III. NUMERICAL EXPERIMENTS
A. Data Generation
We conducted a series of numerical experiments to test theeffectiveness of the proposed method. A coil alignment of the
37-channel Magnes biomagnetic measurement system (Biomagnetic Technologies Inc., San Diego, CA) was used in theseexperiments. The coordinate system used in our numericalexperiments is illustrated in Fig. 1. The values of the spatial
were expressed in centimeters. Three
coordinates
signal sources were assumed to exist on a plane defined as(b)1.0. The locations as well as the orientations of the sources arelisted in Table I. Because a spherical homogeneous conductor
was used, we express[3] with the origin set atthe source-moment vector using the two tangential components.
The magnetic field was generated at a 1-ms interval from
300 to 400 ms. The moment time courses of the three sources
are shown in Fig. 2(a) for the time window between 0 and
400 ms. The white Gaussian noise was added to the generated
magnetic field, and the SNR, defined as the ratio of the Frobenius norm of the signal-magnetic-field data matrix to that of

766

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

Fig. 7. The x, y , and z coordinates used to express the reconstruction resultsin Section IV. The midpoint between the left and right preauricular points was
defined as the coordinate origin. The axis directed away from the origin toward
the left preauricular point was defined as the y axis, and that from the origin to
the nasion was the x axis. The z axis was defined as the axis perpendicular
to both these axes and was directed from the origin to the vertex.

(a)

+

+

+

Note that since these weight vectors are calculated for any spatial location , the source-moment distribution at any location
can be reconstructed in a perfectly post-processing manner. The
reconstruction region was set as an area defined byandon the plane
1 (as indicated by thesquare in Fig. 1), and the reconstruction interval was 1 mm in
the and directions.
andwere obtained, an angle repreOnce
plane, , wassenting the mean source direction in thecalculated using
if

if

(b)
Fig. 6. (a) Results of the spatio-temporal reconstruction obtained using
the proposed vector beamformer technique [(11) and (12)]. The upper-left,
upper-right, and lower-left maps respectively show the snapshots of thesource-moment magnitude at 220, 268, and 300 ms. The lower-right map
shows the time-averaged reconstruction. (b) Estimated time courses from thefirst to the third sources are shown from the top to the bottom, respectively. The
three vertical broken lines indicate the time instants of 220, 268, and 300 ms.

the noise matrix, was set to 18. The generated magnetic field
is shown in Fig. 2(b). This SNR is higher than that in typicalcases of actual MEG measurements. We used such a high
SNR value because the differences in the source estimation
results obtained using tested beamformer techniques can be
more easily observed under such high SNR conditions; suchdifferences might otherwise be obscured by noise effects.
B. Spatio-Temporal Reconstruction Experiments
The spatio-temporal reconstruction was performed by using
and

(17)

(18)

where indicates the average over the time window with which
was calculated. Then, the time course expressed in the mean
, and that in its orthogonal direction,
source direction,
were given by

(19)
and
In the following experiments, we used
when displaying the time course of a source activity.
To display the results of the spatio-temporal reconstruction,
three time points at 220, 268, and 300 ms [marked in Fig. 2(a)]were selected. The amplitude of the second source happened to
be zero at 220 ms, all the sources had nonzero amplitudes at
268 ms, and only the second source had a nonzero amplitudeat 300 ms. The snapshots of the source magnitude distribution
at these three time points, andthe time averaged reconstruction
the following experiments.

were displayed inSEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES

767

Fig. 8. The auditory-somatosensory combined response measured by simultaneously applying an auditory stimulus and a somatosensory stimulus. The auditory
stimulus was a 1-kHz pure tone delivered to the subject’s right ear and the somatosensory stimulus was a 30-ms-duration tactile pulse delivered to the distal segmentof the right index finger. A total of 256 epochs were averaged.

Fig. 9. Results of the spatio-temporal reconstruction from the auditory-somatosensory combined response in Fig. 8. The reconstruction was obtained using the
proposed beamformer. The maximum intensity projection of the source-moment magnitude onto the axial (left), coronal (middle), and sagittal (right) slices are
displayed. The source-magnitude distributions are shown at latencies of (a) 65 ms, (b) 138 ms, and (c) 194 ms. These time instants are shown by the vertical brokenlines in Fig. 10.

C. Results from Minimum-Variance Vector Beamformer
The results of the spatio-temporal reconstruction obtainedusing the minimum-variance vector beamformer in (6) with
the normalized lead field matrix are shown in Fig. 3(a).
The estimated time courses at the pixels nearest to the threesource locations are shown in Fig. 3(b). These results showthat the reconstruction at each instant in time was fairly
noisy: the snapshot at 220 ms showed some influence fromthe second source, and the snapshot at 300 ms contained
the activities of the first and third sources. The time-averaged reconstruction, however, clearly resolved three active
sources. Note that this time-averaged reconstruction is equal
to the map of the neural activity index proposed in [10],768IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

Fig. 10. Time courses of the points nearest to (a) the primary somatosensory cortex and (b) the primary auditory cortex. The solid and broken plotted lines
correspond, respectively, to s^ (r ; t) and s^ (r ; t). Three vertical broken lines indicate the time instants of 65, 138, and 194 ms.

Fig. 11. The auditory-button-press response. The data was measured with a subject who pressed a response button with his left index finger when he heard the
1-kHz pure tone delivered to his right ear. A total of 256 epochs are averaged.

because the noise was white Gaussian and the noise covariance matrix was expressed as the unit matrix in thesenumerical experiments.
It is known that this poor output SNR is due to the use of[11], [12], [16]. Thus, we next
the direct matrix inversion
tested the minimum-variance beamformer together with the use
instead of. Theof the regularized inverse, whereis theregularization parameter was set at. The results in Fig. 4(a) show that alargest eigenvalue ofconsiderable amount of blur was introduced. The estimated timecourses are shown in Fig. 4(b). Fig. 4 shows that the SNR of
the beamformer output was considerably increased in this case,although each time course shows some influence from neighboring sources. The results here demonstrated that the regularization leads to a tradeoff between the spatial resolution and the
SNR of the beamformer output.
D. Results from Proposed Vector Beamformer
We first show the reconstruction results from weight vectors
obtained using (11) alone. This is equivalent to the vector-extended Borgiotti–Kaplan beamformer without the eigenspace
projection. The results are shown in Fig. 5. Comparison between the time-averaged reconstruction in Fig. 3(a) and that inFig. 5(a) confirms that the Borgiotti–Kaplan-type beamformer
has a spatial resolution much higher than the minimum-variance beamformer. The spatio-temporal reconstruction, however,is very noisy for both cases.

We then applied the proposed vector beamformer obtainedusing (11) and (12) to the same computer-generated data set.
The reconstructed source distributions are shown in Fig. 6(a),and the estimated time courses are shown in Fig. 6(b). Comparison between Figs. 5 and 6 confirms that the eigenspace projection can improve the SNR with almost no sacrifice of thespatial resolution. Comparing the results in Fig. 6 with the minimum-variance results in Fig. 3, we can clearly see that the proposed beamformer technique significantly improved both spatial resolution and output SNR.
IV. APPLICATION TO AUDITORY-EVOKED MEG DATA
We applied the proposed beamformer technique to two setsof auditory-evoked MEG data to demonstrate its spatio-temporal reconstruction capability. The auditory-evoked fields were
measured using the 37-channel Magnes biomagnetometer installed at the Biomagnetic Imaging Laboratory, University ofCalifornia, San Francisco. The auditory stimulus was presented
to the subject’s right ear. The sensor array was placed above thesubject’s left hemisphere with the position adjusted to optimallyrecord the N1m auditory-evoked field. The average inter-stimulus interval was 2 s, with the interval randomly varied between
1.75 s and 2.25 s. The sampling frequency was set at 1 kHz. An
on-line filter with a bandwidth from 1 to 400 Hz was used, andno post-processing digital filter was applied. To express the results of reconstructing source activities in this section, we used
the head coordinate system illustrated in Fig. 7.

SEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES

769

Fig. 12. Results of the spatio-temporal reconstruction from the auditory-button-press response in Fig. 11. The maximum intensity projection of the source-moment
magnitude onto the axial (left), coronal (middle), and sagittal (right) slices are displayed. The source-magnitude distributions are shown at latencies of (a) 100 ms,(b) 170 ms, and (c) 200 ms. These time instants are shown by the vertical broken lines in Fig. 13.

Fig. 13. Time courses of the points nearest to (a) the primary auditory cortex and (b) the center of the motor activities. The solid and broken plotted lines
correspond, respectively, to s^ (r ; t) and s^ (r ; t). Three vertical broken lines indicate the time instants of 100, 170, and 200 ms.

The first data we tested was the auditory-somatosensory
combined response measured by simultaneously applying anauditory stimulus and a somatosensory stimulus to a malesubject. The auditory stimulus was a 1-kHz pure tone witha 200-ms duration, and the somatosensory stimulus was a
30-ms-duration tactile pulse (17 psi) delivered to the distalsegment of the right index finger. These two stimuli started
at the same time. A total of 256 epochs were measured,and the response averaged over all the epochs is shown in
Fig. 8. We applied the proposed beamformer to this aver-aged data. The data in the time window ranging from 0 to.
300 ms was used for calculating the covariance matrix
was set at two because
The signal subspace dimension
showed two distinctly large
the eigenvalue spectrum of
eigenvalues.
The reconstructed source-magnitude maps at three latencies,
65, 138, and 194 ms, are shown in Fig. 9. The source magnitude
map at 138 ms [Fig. 9(b)] contains a source activity presumably at the primary somatosensory cortex. The source magnitude map at 194 ms [Fig. 9(c)] shows a source activity at the770

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 48, NO. 7, JULY 2001

primary auditory cortex. The map at 65 ms [Fig. 9(a)] contains
both of these activities.
The time courses of the points at the primary somatosensory
and auditory cortices are shown in Fig. 10(a) and (b), respectively. The coordinates of these cortices were determined fromandthe maximum points in Fig. 9(b) and (c). Bothwere plotted in Fig. 10 so that we could check whether the orientation of the source was fixed or changed during the time
was calculated. In Fig. 10(a) the P50
window with whichpeak, which is known to represent activity of the primary somatosensory cortex, is observed near the latency of 50 ms. In
Fig. 10(b), the auditory N1m peak is observed near the latencyof 100 ms. These time course plots also show that both sourceswere active at 65 ms, that only the primary somatosensory area
was active at 138 ms, and that only the primary auditory area
was active at 194 ms. These observations are consistent withthe source-activity behavior shown in Fig. 9(a)–(c).
We next applied the proposed beamformer to auditory-button-press data. When we measured this data set, thesubject pressed a response button with his left index finger
when he heard the 1-kHz pure tone. Thus, the data should
contain motor activities fairly closely time-locked to theauditory stimulus. The results obtained by averaging a total of
256 epochs are shown in Fig. 11. The proposed beamformer
was applied to this averaged auditory-button-press data. The
data from a time window ranging from 0 to 300 ms was used
and the signal subspace dimension was set
for calculating
at two. The results of the reconstruction are shown in Fig. 12.
Here, the reconstructed source-magnitude distribution at 100
ms is shown in Fig. 12(a), that at 170 ms is shown in Fig. 12(b),and that at 200 ms is shown in Fig. 12(c). The reconstructedresults at 100 ms contained clear activity near the primary
auditory cortex in the left temporal area. The results at 170 ms
showed that the activity at the primary auditory area was still
dominant. The results at 200 ms showed wide-spread activities,
presumably around motor and premotor areas.
We determined the location of the primary auditory area bychoosing the maximum point in Fig. 12(a), and the location of
the center of the motor activities by choosing the maximum
point in Fig. 12(c). The time course of the activity in the primary
auditory area is shown in Fig. 13(a). The time course forms apeak near the latency around 100 ms, clearly showing the auditory N1m component. The auditory activity was relatively weak,
but still dominant at the latency of 170 ms. The time course of
the motor activities is shown in Fig. 13(b), which indicates thatthe motor activity formed a peak near 200 ms. These plots arealso consistent with the results in Fig. 12.
V. CONCLUSION
We have developed a novel MEG vector beamformer technique suitable for reconstructing spatio-temporal activities ofneural sources. The developed beamformer is formulated in a
two-step procedure: the first step extends the Borgiotti–Kaplan
beamformer to the vector-type beamformer and the second step
projects its weight vectors onto the signal subspace of the measurement covariance matrix. The proposed beamformer has thespatial resolution and the output SNR, both significantly higher

than those of the minimum-variance vector beamformer usedin the previous investigations. Our numerical experiments verified the superiority of the proposed method, and its application
to two sets of auditory MEG data demonstrated the method’s
spatio-temporal reconstruction capability.
ACKNOWLEDGMENT
The authors would like to thank Dr. T. Roberts and S. Honma
for their help in performing the auditory MEG measurements.
REFERENCES[1] M. Hämäläinen, R. Hari, R. J. IImoniemi, J. Knuutila, and O. V.
Lounasmaa, “Magnetoencephalography-theory, instrumentation, andapplications to noninvasive studies of the working human brain,” Rev.
Mod. Phys., vol. 65, pp. 413–497, 1993.
[2] T. P. L. Roberts, D. Poeppel, and H. A. Rowley, “Magnetoencephalography and magnetic source imaging,” Neuropsychiatry, Neuropsych..,Behavioral Neurol., vol. 11, pp. 49–64, 1998.
[3] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[4] M. S. Hämäläinen and R. J. Ilmoniemi, “Interpreting measured magnetic
fields of the brain: Estimates of current distributions,” Helsinki Univ.
Technol., Helsinki, Finland, Tech. Rep. TKK-F-A559, 1984.
[5] R. D. Pascual-Marqui and C. M. Michel, “Low resolution electromagnetic tomography: A new method for localizing electrical activity in thebrain,” Int. J. Psychophysiol., vol. 18, pp. 49–65, 1994.
[6] K. Uutela, M. Hämäläinen, and E. Somersalo, “Visualization ofmagnetoencephalographic data using minimum current estimate,”
NeuroImage, vol. 10, pp. 173–180, 1999.
[7] I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm,”
IEEE Trans. Signal Processing, vol. 45, pp. 600–616, Mar. 1997.
[8] B. D. Van Veen and K. M. Buckley, “Beamforming: A versatile approachto spatial filtering,” IEEE Acoust. Speech, Signal Processing Mag., vol.
5, pp. 4–24, Apr. 1988.
[9] M. E. Spencer, R. M. Leahy, J. C. Mosher, and P. S. Lewis, “Adaptive
filters for monitoring localized brain activity from surface potential timeseries,” in Proc. 26th Annu. Asilomer Conf. Signals, Systems, and Computers, Nov. 1992, pp. 156–161.
[10] B. D. Van Veen, W. van Drongelen, M. Yuchtman, and A. Suzuki,“Localization of brain electrical activity via linearly constrained
minimum variance spatial filtering,” IEEE Trans. Biomed. Eng., vol.
44, pp. 867–880, Sept. 1997.
[11] S. E. Robinson and J. Vrba, “Functional neuroimaging by synthetic aperture magnetometry (SAM),” in Recent Advances in Biomagnetism, T.
Yoshimoto, et al., Eds. Sendai, Japan: Tohoku Univ. Press, 1999, pp.
302–305.
[12] J. Gross and A. A. Ioannides, “Linear transformations of data space in
MEG,” Phys. Med. Biol., vol. 44, pp. 2081–2097, 1999.
[13] G. Borgiotti and L. J. Kaplan, “Superresolution of uncorrelated interference sources by using adaptive array technique,” IEEE Trans. Antennas
Propagat., vol. AP-27, pp. 842–845, 1979.
[14] J. Capon, “High-resolution frequency wavenumber spectrum analysis,”
Proc. IEEE, vol. 57, pp. 1408–1419, 1969.
[15] K. Sekihara, D. Poeppel, and Y. Miyashita, “Application of eigenspace
beamformer to virtual depth-electrode measurement using MEG,” inProc. 2nd Int. Symp. Noninvasive Functional Source Imaging Within the
Human Brain and Heart (Biomedizinische Technik), S. Supek, Ed., Zagreb, Croatia, Sept. 1999, pp. 127–130.
[16] H. Cox, R. M. Zeskind, and M. M. Owen, “Robust adaptive beamforming,” IEEE Trans. Signal Process., vol. SP–35, pp. 1365–1376,
1987.
[17] B. D. Carlson, “Covariance matrix estimation errors and diagonalloading in adaptive arrays,” IEEE Trans. Aerosp. Electron. Syst., vol.
24, pp. 397–401, July 1988.
[18] L. Chang and C. C. Yeh, “Performance of DMI and eigenspace-based
beamformers,” IEEE Trans. Antennas Propagat., vol. 40, pp.
1336–1347, Nov. 1992.
[19] L. L. Scharf, Statistical Signal Processing: Detection, Estimation, andTime Series Analysis. New York: Addison-Wesley, 1991.

SEKIHARA et al.: RECONSTRUCTING SPATIO-TEMPORAL ACTIVITIES OF NEURAL SOURCES[20] D. D. Feldman and L. J. Griffiths, “A constrained projection approachfor robust adaptive beamforming,” in Proc. Int. Conf. Acoust., Speech,
and Signal Processing, Toronto, Canada, May 1991, pp. 1357–1360.

Kensuke Sekihara (M’88) received the M.S. degreein 1976 and the Ph.D. degree in 1987 both from
Tokyo Institute of Technology, Tokyo, Japan.
From 1976 to 2000, he worked with Central
Research Laboratory, Hitachi, Ltd., Tokyo, Japan.
He was a visiting Research Scientist at Stanford
University, Stanford, CA, from 1985 to 1986, and at
Basic Development, Siemens Medical Engineering,
Erlangen, Germany, from 1991 to 1992. From1996 to 2000, He worked with “Mind Articulation"
research project sponsored by Japan Science andTechnology Corporation. He is currently Professor at Tokyo Metropolitan
Institute of Technology, Tokyo, Japan. His research interests include thebiomagnetic inverse problems, and statistical estimation theory, especially its
application to noninvasive functional neuroimaging.
Dr. Sekihara is a member of the IEEE Medicine and Biology Society, and theIEEE Signal Processing Society.

Srikantan S. Nagarajan received the B.S. degree inelectrical engineering from the University of Madras,
Madras, India, and the M.S. and Ph.D degrees fromthe Department of Biomedical Engineering at Case
Western Reserve University, Cleveland, OH, in 1993
and 1995, respectively.
From 1995-1998, he was a Post-doctoral Fellow in
the Keck Center for Integrative Neuroscience at the
University of California, San Francisco (UCSF). In
1999, he was a full-time Research Scientist at Scientific Learning Corporation, Berkeley, CA, and an
Adjunct Assistant Professor in the Department of Otolaryngology at UCSF. He
is currently an Assistant Professor in the Department of Bioengineering at the
University of Utah, Salt Lake City. His research interests in neural engineering
include bioelectromagnetism, systems and computational neuroscience, and statistical signal processing.

771

David Poeppel received the B.S. and Ph.D degrees incognitive neuroscience from Massachusetts Instituteof Technology, Cambridge, in 1990 and 1995, respectively.
He is currently an Assistant Professor in the
Departments of Linguistics and Biology at the
University of Maryland at College Park. His researchuses the functional neuroimaging methods magnetoencephalography (MEG), electroencephalography(EEG, positron emission tomography (PET), andfunctional magnetic resonance imaging (fMRI) to
investigate the neural basis of speech and language processing.

Alec Marantz received the B.A. degree in psycholinguistics from Oberlin College, Oberlin, OH, in 1978
and the Ph.D. degree in linguistics from the Massachusetts Institute of Technology (M.I.T.), Cambridge,
MA, in 1981.
He joined the faculty at M.I.T. in 1990, where heis currently Professor of Linguistics and Head of the
Department of Linguistics and Philosophy. His research interests include the syntax and morphology ofnatural languages, linguistic universals, and the neurobiology of language. He is currently involved inrevising morphological theory within linguistics and in exploring MEG techniques to uncover how the brain process language.

Yasushi Miyashita received the Ph. D. degree inphysiology from the University of Tokyo, Tokyo,
Japan, in 1979.
He is currently Professor and Chairman of the
Physiology Department, the University of Tokyo
School of Medicine, Tokyo, Japan. His reserach interests include neural basis of cognition in primates,
functional imaging, and image processing with MRI,
MEG, and optical devices.