Journal of Neural Engineering

PAPER

Beamspace dual signal space projection (bDSSP): a method forselective detection of deep sources in MEG measurements
To cite this article: Kensuke Sekihara et al 2018 J. Neural Eng. 15 036026

View the article online for updates and enhancements.

This content was downloaded from IP address 132.174.255.215 on 17/04/2018 at 04:39

Journal of Neural Engineering
J. Neural Eng. 15 (2018) 036026 (17pp)

https://doi.org/10.1088/1741-2552/aab5bd

Beamspace dual signal space projection(bDSSP): a method for selective detectionof deep sources in MEG measurements
Kensuke Sekihara1,2 , Yoshiaki Adachi3, Hiroshi K Kubota4, Chang Cai5and Srikantan S Nagarajan51

Signal Analysis Inc., Hachioji, Tokyo, Japan
Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45
Yushima, Bunkyo-ku, Tokyo 113-8519, Japan
3
Applied Electronics Laboratory, Kanazawa Institute of Technology, Kanazawa, Japan
4
Healthcare Business Center, Business Development Division, Ricoh Ltd., Tokyo, Japan5
Biomagnetic Imaging Laboratory, University of California, San Francisco, 513 Parnassus Avenue, S362,
San Francisco, CA 94143, United States of America
2E-mail: k-sekihara@nifty.com
Received 5 October 2017, revised 9 February 2018
Accepted for publication 12 March 2018
Published 16 April 2018
Abstract

Objective. Magnetoencephalography (MEG) has a well-recognized weakness at detecting
deeper brain activities. This paper proposes a novel algorithm for selective detection of deep
sources by suppressing interference signals from superficial sources in MEG measurements.
Approach. The proposed algorithm combines the beamspace preprocessing method withthe dual signal space projection (DSSP) interference suppression method. A prerequisiteof the proposed algorithm is prior knowledge of the location of the deep sources. The
proposed algorithm first derives the basis vectors that span a local region just covering thelocations of the deep sources. It then estimates the time-domain signal subspace of thesuperficial sources by using the projector composed of these basis vectors. Signals from thedeep sources are extracted by projecting the row space of the data matrix onto the direction
orthogonal to the signal subspace of the superficial sources. Main results. Compared withthe previously proposed beamspace signal space separation (SSS) method, the proposed
algorithm is capable of suppressing much stronger interference from superficial sources. This
capability is demonstrated in our computer simulation as well as experiments using phantom
data. Significance. The proposed bDSSP algorithm can be a powerful tool in studies ofphysiological functions of midbrain and deep brain structures.
Keywords: neuromagnetism, deep sources, electromagnetic brain imaging,magnetoencephalography (MEG), source localization, MEG source space analysis(Some figures may appear in colour only in the online journal)

1. Introduction

sensor arrays with advanced signal processing algorithms, arenow capable of imaging dynamic brain activities with the temporal resolution on the order of submilli-seconds.
In spite of its success, MEG has the well-recognizedweakness that it is not generally good at detecting deepbrain activities, although there have been many attempts onrecording MEG signals from deep brain regions [4] such as

Electrophysiological activity of neurons in the cerebral cortex
generates weak but detectable magnetic fields outside thescalp, and direct non-invasive measurements of these neuronal
activities can be achieved with magnetoencephalography(MEG) [1–3]. Modern MEG systems, combining large-scale
1741-2552/18/036026+17$33.00

1© 2018 IOP Publishing Ltd

Printed in the UKK Sekihara et al

J. Neural Eng. 15 (2018) 036026

2. Dual signal subspace projection algorithm

the thalamus [5–9], amygdala [10, 11], and cerebellum [12,13]. This is primarily because MEG sensor arrays generally
have very low sensitivity near deep brain regions, due to thephysical nature of magnetic signals. That is, since magn­
etic signals attenuate with the inverse-square-distance law,
MEG signals arising from deep sources become significantly
weaker than signals from superficial cortical sources. Also,
it is generally true that when a deeper brain region is active,
superficial cortical regions are simultaneously active. In suchcases, the activity at a deep brain region is masked by overlapped cortical activities, and only these cortical activities aredetected by MEG.
To overcome this weakness, Özkurt et al proposed a methodof decomposing MEG signals into components corresponding
to deep and superficial sources [14]. Their method is based on
the signal space separation (SSS) algorithm, which utilizes thevector spherical harmonics expansions of the measured magn­
etic field [15, 16]. Their method combines the SSS algorithmwith so-called beamspace processing [17, 18]. The method
proposed by Özkurt et al is thus referred to as the beamspace
SSS method herein.
The beamspace processing was originally developed toreduce the data dimensionality. Many popular source localization methods, such as the MUSIC algorithm [19, 20] andadaptive beamforming [21, 22], make use of second-order statistics of the measured data. To obtain stable data statistics, alarge number of time samples are generally needed, and thenumber of required time samples linearly increases as the data
dimension increases. The beamspace processing was developed to reduce this data-sample requirements by reducing
the data dimensionality. The data-dimensionality reduction isachieved by projecting the data vector onto a low-dimensional
subspace [22, 23]. In this paper, the beamspace idea is used
for preprocessing the data by selective interference cancellation. A brief introduction on the beamspace preprocessing ispresented in section 3.1 of this paper.
Inspired by Özkurt et al, this paper proposes a novel algorithm for selective detection of a deep source by suppressing
interference signals from superficial sources. The proposed
method is a combination of the beamspace processing withthe previously-proposed interference suppression method
called dual signal space projection (DSSP) [24]. Thus, theproposed method is here called beamspace DSSP (bDSSP).
Compared with the beamspace SSS, the proposed bDSSP
algorithm is capable of suppressing much stronger interference from superficial sources. This capability is demonstrated
in our computer simulation and experiments using phantom
data.
This paper is organized as follows: following concise
explanations of the DSSP algorithm in section 2, the proposed beamspace DSSP algorithm is described in section 3.
We present computer simulations that test the effectiveness ofthe bDSSP algorithm. In section 5.1, results of experiments
using an MEG phantom are presented to further demonstrate
the effectiveness of the proposed method. In sections 6 and 7,
additional supplementary arguments are presented, including
a comparison between the beamspace SSS and bDSSP
methods.

2.1. Data model

Biomagnetic measurements are usually conducted using asensor array, which simultaneously measures the biomagn­
etic signal with multiple sensors. Let us define the measurement of the mth sensor at time t as ym(t). The measurementfrom the whole sensor array is expressed as a column
vector y(t): y(t) = [y1 (t), y2 (t), . . . , yM (t)]T , which is called
the data vector. Here, M is the number of sensors, and thesuperscript T indicates the matrix transpose. Throughout
this paper, plain italics indicate scalars, lower-case boldface
italics indicate vectors, and upper-case boldface italics indicate matrices.
Let us assume that a unit-magnitude source exists at r
( r = (x, y, z)). When this unit-magnitude source is directed in
the x, y, and z directions, the outputs of the m-th sensor arexyz
respectively denoted as lm(r), lm(r), and lm(r). Let us define
an M × 3 matrix L(r) whose mth row is equal to a 1 × 3 rowyzx
vector [lm(r)]. This matrix L(r), referred to as(r), lm(r), lm
the lead field matrix [1, 25], represents the sensitivity of thesensor array at r.
The DSSP algorithm was proposed in order to remove
interfering magnetic fields overlapped onto signal magnetic
fields. The algorithm assumes the data model:y(t) = yS (t) + yI (t) + ε,(1)where yS (t), (called the signal vector), represents the signal ofinterest, yI (t), (called the interference vector), represents theinterference magnetic field, and ε, (called the random vector),
represents additive sensor noise. We denote the time series
outputs of a sensor array y(t1 ), . . . , y(tK ), where K is the total
number of measured time points. The measured data matrix B
is thus defined as:
B = [y(t1 ), . . . , y(tK )].
(2)The signal matrix is defined as
BS = [yS (t1 ), . . . , yS (tK )],(3)and the interference matrix as
BI = [yI (t1 ), . . . , yI (tK )].
(4)Then, the data model in equation (1) is expressed in a matrix
form as:B = BS + BI + Bε ,(5)where Bε is the noise matrix whose j-th column is equal to thenoise vector ε at time tj.
2.2. Pseudo-signal subspace projector

The DSSP algorithm assumes the data model in equation (5)with the assumption that the interference sources are locatedoutside the source space6. The DSSP algorithm uses the6

The source space indicates a region in which signal sources can exist. For
example, the entire brain (or entire cortex) is considered as the source spacein most MEG measurements.

2

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

so-called pseudo-signal subspace projector, and to derive it,voxels are defined over the source space, in which the voxel
locations are denoted r1 , . . . , rN . The augmented leadfield
matrix over these voxel locations is defined asrsp(P̆S B) ⊂ KS + KI + K̆ε ,(14)rsp((I − P̆S )B) ⊂ KI + K̆ε ,(15)where KS and KI respectively indicate the time-domain
signal and interference subspaces. Here, we use the notations,rsp(P̆S Bε ) = K̆ε and rsp((I − P̆S )Bε ) = K̆ε .
Using equations (14) and (15), we can finally derive therelationship [27]:F = [L(r1 ), . . . , L(rN )],(6)and the pseudo-signal subspace ĔS is defined such thatĔS = csp(F),(7)where the notation csp(X) indicates the column space of amatrix X . If the voxel interval is sufficiently small and voxel
discretization errors are negligible, we have the relationship ĔS ⊃ ES where ES indicates the true signal subspace.
Therefore, a vector contained in the signal subspace is also
contained in the pseudo-signal subspace.
Let us derive the orthonormal basis vectors of the pseudosignal subspace. To do so, we compute the singular value
decomposition of F :KI ⊃ rsp(P̆S B) ∩ rsp((I − P̆S )B).
(16)The equation above shows that the intersection betweenrsp(P̆S B) and rsp((I − P̆S )B) forms a subset of the interference subspace KI. Once the orthonormal basis vectors of theintersection ψ 1 , . . . , ψ r are obtained, we can compute the projector onto the intersection Πisc such that
Πisc = [ψ 1 , . . . , ψ r ][ψ 1 , . . . , ψ r ]T .
(17)Using this Πisc as the projector onto the interference subspace KI, the interference removal is achieved and the signal
matrix is estimated such thatMF=
λj ej f Tj .
(8)
j=1

 S = B(I − Πisc ) = B(I − [ψ 1 , . . . , ψ r ][ψ 1 , . . . , ψ r ]T ).
B
(18)The method of removing the interference in a manner
described above is called DSSP [24]. Note that since the basis
vectors of the intersection, ψ 1 , . . . , ψ r , span only a subsetof the interference subspace KI, this method cannot perfectly remove interferences. However, when the intersection
rsp(P̆S B) ∩ rsp((I − P̆S )B) is a reasonable approximationof KI, interferences can effectively be removed by the DSSP
algorithm.

Note that F is an M × N matrix, so ej is an M × 1 column
vector and f j is an N × 1 column vector. In equation (8),we assume the relationship M  <  N, and the singular values
are numbered in decreasing order. If the singular values
λ1 , . . . , λτ are distinctively large and other singular values
λτ +1 , . . . , λM are nearly equal to zero, the leading τ singular
vectors e1 , . . . , eτ form orthonormal basis vectors of thepseudo-signal subspace ĔS [26]. Thus, the projector onto ĔS
is obtained usingP̆S = [e1 , . . . , eτ ][e1 , . . . , eτ ]T .
(9)3. Beamspace dual signal subspace projection(bDSSP) algorithm

Note that, since ĔS ⊃ ES, the orthogonal projector (I − P̆S )removes the signal vector, i.e. (I − P̆S )yS (t) = (I − P̆S )BS = 0.

3.1. Beamspace processing
2.3. DSSP algorithm

Beamspace processing refers to a signal processing algorithm
used for data-dimensionality reduction. Such data-dimensionality reduction is achieved by projecting the data vector onto
a low-dimensional subspace. In other words, an M × 1 data
vector y(t) is represented with known basis vectors u1 , . . . , uP,
where the number of basis vectors P is smaller than the dimension of the data vector M. That is, if y(t) is expressed using alinear combination of a set of known P basis vectors such that

A detailed explanation of the DSSP algorithm in the context
of the time-domain signal subspace can be found in [27]. The
DSSP algorithm applies P̆S and I − P̆S to the data matrix B tocreate two kinds of data matrices:(10)P̆S B = BS + P̆S BI + P̆S Bε ,(I − P̆S )B = (I − P̆S )BI + (I − P̆S )Bε .
(11)

P

y(t) ≈
cj (t)uj ,(19)To derive equations (10) and (11), we use P̆S BS = BS and(I − P̆S )BS = 0. Let us use the notation rsp(X) to indicate therow space of a matrix X . Then, according to the arguments in[27], the following relationships hold:j=1

the sensor measurements y1 (t), y2 (t), . . . , yM (t) can be represented by only P coefficients c1 (t), . . . , cP (t). Since weassume P  <  M, the data dimension is reduced from M to P in
equation (19).
The problem here is how to find basis vectors u1 , . . . , uP
which satisfy the relationship in equation (19). A methodof deriving the basis vectors based on the prior knowledgeof signal source locations has been proposed in [23]. In thisrsp(P̆S B) ⊂ rsp(BS ) + rsp(P̆S BI ) + rsp(P̆S Bε ),(12)rsp((I − P̆S )B) ⊂ rsp((I − P̆S )BI ) + rsp((I − P̆S )Bε ).
(13)Since the relationships, rsp(P̆S BI ) = KI and rsp(I − P̆S )BI =KI, hold, equations (12) and (13) lead to the relationships:
3

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 1. Locations of sensors and three sources assumed inour computer simulations. Filled circles indicate the locations of
sensors, the blank circle indicates the location of the deep source,
and the square and triangle indicate the locations of the two
superficial sources. The cross mark indicates the sphere origin. Topleft, top-right, and bottom panels, respectively, show the coronal,
sagittal, and axial projections. The source locations, orientations
and distances to the sphere origin are shown in table 1 in which
Source#2 is the deep source, and the other two sources are called
the superficial sources in our computer simulations.

Figure 2. Time courses assigned to the three sources. Time t isexpressed with the unit of time points. Top panel shows the timecourse assigned to the deep source, and the bottom two panels show
the time courses assigned to the two superficial sources. Each timecourse is normalized between  −1 and 1.

signals from a deep source by suppressing interference fromsuperficial sources. The algorithm is called bDSSP, which is
described as follows: The data model for the bDSSP algorithm
is expressed as

Table 1. Source properties assumed in computer simulation.

SourcenumberOrientation

Location (cm)Distance from thesphere origin (cm)12
3

(1, −1, 0)(0, 1, 1)(1, 0, 1)(0.4, 0.1, 13.3)(2.9, −1.4, 9.3)
(2.4, 2.6, 11.3)

7.0
4.4
6.1

B = Bdeep + Bsup + Bε ,(22)where Bdeep indicates the signal caused from a deep source7,
and Bsup the signal from superficial sources. A prerequisite ofthis algorithm is that the location of the deep source be known.
We then compute the beamspace basis vectors u1 , . . . , uP by
setting the local source space as a small region just covering
the location of the deep source.
The beamspace projector Pdeep is then derived as

proposed method, the augmented lead field matrix F̄ is definedover a local region that just contains the signal sources. The
voxels are defined over this local region and the voxel locations are denoted r̄1 , . . . , r̄N̄ . The augmented leadfield matrix
over these voxel locations is expressed as

Pdeep = [u1 , . . . , uP ][u1 , . . . , uP ]T .
(23)By multiplying Pdeep and I − Pdeep with the data matrix B , weobtain:F̄ = [L(r̄1 ), . . . , L(r̄N̄ )],(20)Pdeep B = Bdeep + Pdeep Bsup + Pdeep Bε ,(24)and its singular value decomposition is given by(I − Pdeep )B = (I − Pdeep )Bsup + (I − Pdeep )Bε .
(25)

R

T
F̄ =
λ̄j ēj f̄ j
(21)Here, we use Pdeep Bdeep = Bdeep . Then, by using exactly thesame derivations as in [27], we can derive

j=1

where R = min{M, N̄}. Let us assume that the leading τ̄ singular values, λ̄1 , . . . , λ̄τ̄ are distinctively large, compared to
the rest of the singular values λ̄τ̄ +1 , . . . , λ̄R. Then, the beamspace basis vectors u1 , . . . , uP are obtained as the leading τ̄singular vectors ē1 , . . . , ēτ̄ where P is equal to τ̄ .

rsp(Pdeep B) ⊂ Ksup + Kdeep + K̄ε ,(26)rsp((I − Pdeep )B) ⊂ Ksup + K̄ε ,(27)where Kdeep and Ksup are the time-domain signal subspaces ofthe deep and superficial sources, respectively. We also use thenotations: rsp(Pdeep Bε ) = K̄ε and rsp((I − Pdeep )Bε ) = K̄ε .
Using equations (26) and (27), We can finally derive,

3.2. bDSSP algorithm

Use of the beamspace basis vectors with the DSSP algorithm leads to a novel algorithm that can selectively detect

74

We assume that the target deep source is a single source for simplicity.

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Ksup ⊃ rsp(Pdeep B) ∩ rsp((I − Pdeep )B).
(28)The equation above indicates that the intersection between therow spaces of Pdeep B and (I − Pdeep )B forms a subset of Ksup .
The orthonormal basis set of the intersection,rsp(Pdeep B) ∩ rsp((I − Pdeep )B), can be obtained using theprocedure presented in [24, 27, 28]. Denoting these orthonormal basis vectors by φ1 , . . . , φr , the projector onto theintersection is obtained asΠisc = [φ1 , . . . , φr ][φ1 , . . . , φr ]T .

Using this Πisc as the projector onto the signal subspace Ksup ,the signal from the deep source is estimated by projecting therows of the data matrix onto the direction orthogonal to Ksup ,such that deep = B(I − Πisc ) = B(I − [φ1 , . . . , φr ][φ1 , . . . , φr ]T ).
B
(29)Note that since the basis vectors φ1 , . . . , φr span only a part ofKsup , the orthogonal projection on the right-hand side of equation (29) cannot perfectly remove Bsup . Nonetheless, whenrsp(Pdeep B) ∩ rsp((I − Pdeep )B) ≈ Ksup holds, a significant
reduction of Bsup can be attained. One such situation is thecase in which the superficial signal Bsup is dominated in B .
The procedure of the bDSSP algorithm can be summarized
as follows:(1) Set the local source space so that it covers the locationof the target deep source, and compute the beamspace
projector Pdeep using equation (23).
(2) Apply the projector Pdeep to the data matrix B to create
two data sets Pdeep B and (I − Pdeep )B.
(3) Apply the singular value decomposition to Pdeep B toderive the orthonormal basis set of rsp(Pdeep B); the basis
set is denoted by {x1 , . . . , xµ } where μ is the dimension
of rsp(Pdeep B). Also, apply the singular value decomposition to (I − Pdeep )B to derive the orthonormal basis set ofrsp((I − Pdeep )B); the basis set is denoted by {y1 , . . . , yν }where ν is the dimension of rsp((I − Pdeep )B).
(4) Define matrices X = [xT1 , . . . , xTµ ] and Y = [yT1 , . . . , yTν ],and compute the singular value decomposition of XT Y ,such thatγ1 · · · 0..  T T ,..
XT Y = Q  ...
(30)
.
.
0

...

Figure 3. (a) Results of experiments when the signal tointerference ratio (SIR), γ, was set to one. Top panel: simulated
sensor time courses generated when one deep and two superficial
sources are active. Middle panel: sensor time courses generated
when only the deep source is active with the same amount of thesensor noise added. Bottom panel: results of applying the bDSSP
algorithm to the sensor data in the top panel. (b) Results of thesame experiments except that SIR (γ) was set to 0.05. Top panel:
simulated sensor time courses. Bottom panel: results of applying
the bDSSP algorithm to the sensor data in the top panel. Sensor
time courses in each panel are normalized to each maximum value,
and the ordinate of the figure indicates the values normalized
between  −1 and 1.

(8) Apply any of the source localization methods to obtainsource images.

γν		where Q and T are matrices whose columns consist of thesingular vectors.
(5) Observe the relation:

4. Computer simulations

γ1 = γ2 = · · · = γr ≈ 1 > γr+1  · · ·  γν ,A series of computer simulations were carried out to test theeffectiveness of our proposed algorithm. A sensor alignmentof the 275-channel whole-head sensor array from the Omega(VMS Medtech, Coquitlam, Canada) neuromagnetometer wasused. The sensor lead field was computed using the Sarvas
formula [25] with the sphere origin set at (0.1, 0, 6.3). In our
computer simulation, three sources, (two superficial sourcesand one deep source), were assumed. The locations of the		and determine the dimension of the intersection
rsp(Pdeep B) ∩ rsp((I − Pdeep )B) to be r.
(6) Obtain the orthonormal basis set of the intersection
rsp(Pdeep B) ∩ rsp((I − Pdeep )B) as the first r columns of
the matrix XQ, which are denoted φ1 , . . . , φr .
(7) Estimate the sensor time courses of the deep source byapplying equation (29).
5

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 5. (a) Results of applying the bDSSP algorithm tosensor data with SNR equal to 1.5. (The SNR is defined asBdeep F /Bε F) These sensor time courses are normalized
between  −1 and 1, and the ordinate indicates the normalized
values. (b) Source reconstruction results obtained from the bDSSPprocessed sensor data in (a). In these results, the top left, top rightand bottom panels, respectively, show the axial, coronal and sagittal
projections of the 3D source reconstruction results.

by projecting the source time courses in figure 2 through thesensor lead field.
We first performed simulations with a moderately highsignal-to-interference ratio. Simulated sensor measurements
were computed by adding the sensor time courses generated
from the deep source and those from the superficial source.
We set the ratio γ = Bdeep F /Bsup F to be one at which
XF indicates the Frobenius norm of a matrix X 8. This ratio γ
is called the signal to interference ratio (SIR), because, in our
experiments, the magnetic field from the superficial sources
is considered the interference, and that from the deep source
considered the signal. The simulated sensor measurements
given by B = Bdeep + Bsup + Bε are shown in the top panel of
figure 3(a). Here, a moderare amount of white Gaussian noise
was added to simulate the sensor noise and the SNR defined asBdeep F /Bε F was set at 4.5. The middle panel of figure 3(a)shows Bdeep + Bε, which are the sensor time courses generated from only the deep source with the same amount of thesensor noise added. These time courses work as the groundtruth for our experiments described below.
We applied the proposed bDSSP algorithm to the simulated measurements shown in the top panel, and the results

Figure 4. (a) Source reconstruction results obtained using thesimulated sensor time courses in the top panel of figure 3(b). The
cross marks indicate the location of the deep source. (b) Sourcereconstruction results from the bDSSP-processed sensor data in thebottom panel of figure 3(b). In these results, the top left, top rightand bottom panels, respectively, show the axial, coronal and sagittal
projections of the 3D source reconstruction results.

three sources are shown with sensor locations in figure 1. In
this figure, the filled circles indicate the locations of the sensors; the location of the deep source is indicated by the blankcircle; and the locations of the two superficial sources areindicated by the square and the triangle. The sphere origin isindicated by a cross mark.
The source properties (the locations, orientations, and distances from the sphere origin) are summarized in table 1. Note
that, in our computer simulations, Source#2 is considered
the deep source and the other two sources are considered thesuperficial sources. The distance of Source#2 to the sphereorigin was set to 4.4 cm, and distances of the other sources
to the sphere origin were set between 6 and 7 cm. The timecourses assigned to these three sources are shown in figure 2.
The time t is expressed with the unit of time points ranging
from t  =  −1200 to 1200. Sensor time courses were generated

8
As defined in equation (22), Bsup and Bdeep are the signal matrices corre­
sponding to superficial and deep sources.

6

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 6. Results of Monte Carlo experiments on the mismatch
between the position of the local source space and the location of
the target deep source. Top panel: plots of the sensor-time course
correlation Φ versus the mismatch distance D. Bottom panel: plotsof the localization error of the deep source versus D. The position of
the local source space was set at fifty different locations having thesame D. A filled dot represents the results of each Monte Carlo trial,and blank circles indicate the average of the Monte Carlo results ateach D.

are shown in the bottom panel of figure 3(a). In the application, the local source space was chosen as a 2 cm cubic region
whose center was equal to the location of the deep source. The
results show that by using the bDSSP algorithm, time courses
nearly identical to those of the ground truth (in the middle
panel) were obtained. These results show that the bDSSP
algorithm effectively suppresses the interference from thesuperficial source and the signal sensor time courses from thedeep-source can be retrieved.
In our computer simulations, the correlation between the
groun truth (the sensor time courses Bdeep + Bε in the middle
panel of figure 3(a)) and the bDSSP-processed sensor time deep (in the bottom panel of figure 3(a)) is computed
courses B
such that

Figure 7. Results of Monte Carlo experiments on the robustness ofthe bDSSP algorithm to the source configuration. (a) One hundred
randomly-selected locations of the deep source generated in thisMonte Carlo study. Each location is shown by a cross mark. Twoblank circles show the locations of the superficial sources, which
are the same as those used in the previous experiments. (b) Top
panel: histogram of the one-hundred trials with respect to the sensor
time course correlation Φ . Bottom panel: histogram of the onehundred trials with respect to the source localization errors. In thisstudy, SIR was set at 0.05, and SNR was set at 4.5.

M
1 
Φ
=
ρm ,
(31)
M

We next conducted computer simulations for a case ofsignificantly low SIR of 0.05. That is, the interference fromsuperficial sources was 20 times stronger than the signal fromthe deep source. The generated sensor time courses in thiscase are shown in the top panel of figure 3(b). The results ofapplying the bDSSP algorithm are shown in the bottom panel.
The results here also attain Φ of 0.85, and this fact indicatesthat the proposed bDSSP algorithm still effectively removed
the interference from the superficial sources even when SIR
was so low.

m=1

where ρm is the correlation coefficient between the mth row
 deep and that of Bdeep + Bε. Here, Φ is computed by averof B
aging ρm across all sensor channels. This Φ can be a measureindicating how well the bDSSP algorithm can retrieve thesensor time courses from the target deep source, and is called
the sensor time-course correlation. The results in the bottom
of figure 3(a) attain Φ of 0.85. We use this Φ to assess theresults in our Monte Carlo experiments below.
7

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 8. Configuration of a dry phantom used in our experiments. (a) The squares show the locations of dipole sources. Dipole pairs
annotated by ‘dipole pair#1’ and ‘dipole pair#2’ were those used in our experiments. Dipole pair#1 consists of dipole 4-1 and dipole
4-2, and the orientations of these dipoles are perpendicular to each other. Dipole pair#2 consists of dipole 2-1 and dipole 1-2, and theorientations of these dipoles are parallel. The center of the disc matches the sphere origin of the spherical conductor. (b) Isoscelestriangular-shaped coils that simulate dipole sources in the spherical homogeneous conductor can be seen. (c) Depiction showing how thedisc-shaped phantom was installed inside the sensor helmet.

of the sensor noise on the algorithm performance. The same
computer simulation as for the low SIR experiments mentioned above was used. However, in this case, the sensor noise
was added such that the SNR (defined as Bdeep F /Bε F)
was approximately equal to 1.5. The bDSSP-processed sensor
time courses are shown in figure 5(a). These results showthat the sensor time courses of the deep-source still can be
retrieved. The results of source reconstruction obtained fromthese sensor data are shown in figure 5(b). Again, the deep
source is clearly reconstructed, although these results have asmall localization error of 0.7 cm, which was obtained as the
distance between the assumed and reconstructed locations of
the deep source9.
In computer simulations so far, the local source space
was chosen as a 2 cm cubic region whose center was exactly
equal to the location of the deep source. However, in real-life
applications, the estimated location of the target deep source
may have some uncertainty. Therefore, we conducted Monte
Carlo computer simulation to see the robustness of the bDSSP
algorithm to the mismatch between the position of the local
source space and the location of the target deep source. In
this simulation, the distance between the deep source and theWe performed source reconstruction experiments using
these simulated sensor data with SIR of 0.05. The adaptive
beamformer reconstruction algorithm [22] was applied to thedata in figures 3(b). The results of the reconstruction experiments are shown in figure 4. In these experiments, voxels
were defined with 0.5 cm interval. The locations of the threesources were set at one of voxel locations. In these resultsin figure 4, the sources were reconstructed exactly at theirassumed locations.
The source image from the simulated sensor measurementsin the top panel of figure 3(b) is shown in figure 4(a). Here, thetop left, top right and bottom panels, respectively, show theaxial, coronal and sagittal projections of the 3D source distribution. In these results, two superficial sources are reconstructed, but the deep source was not clearly reconstructedbecause of a large intensity difference between the superficial
and deep sources. Note that the cross marks in figure 4(a) indicate the location of the deep source. The source reconstruction
results obtained using the bDSSP-processed sensor data in thebottom panel of figure 3(b) are shown in figure 4(b). Here,the target deep source is successfully reconstructed, demonstrating the effectiveness of the proposed algorithm.
We also performed source reconstruction experiments
using low SIR and low SNR data, in order to see the effect9

The reconstructed source location was defined as the location of the peakmaximum.

8

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 9. (a) Sensor time courses measured when only the superficial dipole (marked ‘4-1’ in figure 8(a)) was active. The superficial
dipole was driven by an 11 Hz sinusoid. (b) Sensor time courses measured when only the deep dipole (marked ‘4-2’ in figure 8(a)) was
active. The deep dipole was driven with an amplitude-modulated sinusoid whose carrier frequency was 15 Hz and modulation frequency
was 1 Hz. These sensor time courses in (a) and (b) are normalized to each maximum field intensity. The abscissa of these figures indicatestime in ms. (c) Source reconstruction results of the superficial dipole. (d) Source reconstruction results of the deep dipole. In each figure,the left, middle, and right panels, respectively, show the axial, coronal and sagittal projections of the 3D source reconstruction results.

center of the local source space was set at several non-zerovalues, keeping the local source space to be the same 2 cm
cubic volume. This distance is called the mismatch distanceand denoted by D. Here, to assess the algorithm performance,
we used two measures: the sensor time-course correlation Φin equation (31) and the localization error of the deep source.
Results of the experiments are shown in figure 6. Here, thetop panel shows the plots of the sensor-time course correlation
Φ versus the mismatch distance D. The bottom panel indicatesthe plots of the localization error of the deep source versus
D. In these Monte Carlo experiments, at each value of D, theposition of the local source space was set to fifty random locations having the same D. In the plots in figure 6, a filled dot
indicates the results of each Monte Carlo trial, and a blankcircle indicates the average across the results at each valueof D. These plots show that when D is less than 1 cm, there isnearly no influence and even when D  =  1.5 cm, the influenceis still very small, demonstrating the algorithm’s robustness
to the mismatch between the local source space and the target
source location.
Finally, we conducted Monte Carlo simulation to see therobustness of the bDSSP algorithm to the source configuration.

9

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 10. (a) The sensor time courses measured when the superficial and deep dipoles were simultaneously active. The magnetic

field from the superficial dipole was sixteen times stronger than the magnetic field from the deep dipole. (b) The sensor time coursesof the bDSSP results. These sensor time courses in (a) and (b) are normalized to each maximum field intensity. The abscissa of thesefigures indicates time in ms. (c) Source reconstruction results from the sensor data in (a). (d) Source reconstruction results from the bDSSPprocessed sensor data in (b). The left, middle, and right panels, respectively, show the axial, coronal and sagittal projections of the 3D
source reconstruction results.

equal to 0.05. In these results, the 96% of the Monte Carlotrials attain the values of Φ between 0.8 and 0.9. Also, 97% of
the Monte Carlo results attain the localization errors less than
0.8 cm. Note that, in this Monte Carlo study, since the deep
source location was randomly chosen, the source was not necessarily located at one of voxel locations, and the localization
error also contains contributions from the voxel discretization,
which is 0.25 cm in average for this study that uses the voxel
interval of 0.5 cm. The results of these Monte Carlo experiments demonstrate that the performance of the bDSSP algorithm is not much affected by the configuration of the deep
source with respect to superficial sources.

In this simulation, one hundred Monte Carlo trials were generated, and in each trial, the location and orientation of the
deep source were randomly chosen. One hundred locationsof the deep source generated in this Monte Carlo study are
shown by cross marks in figure 7(a). In this study, a random
value between 0 and 1 cm was assigned to the mismatch distance D. The results of these one hundred Monte Carlo trials
were assessed by using the sensor time course correlation Φ
and the source localization error.
In figure 7(b), the histogram of Φ is shown at the top panel,and the histogram of source localization errors is shown at
the bottom panel. Here, SNR was set equal to 4.5 and SIR (γ)

10K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 11. Results of experiments using a pair of dipoles denoted ‘Dipole pair#2’; the deep dipole is denoted ‘2-1’, and the superficial

dipole ‘1-2’. (a) The reconstruction results of the superficial dipole. (b) The reconstruction results of the deep dipole. The left, middle, andright panels, respectively, show the axial, coronal and sagittal projections of the 3D source reconstruction results.

strength of 1.42 mA. The deep dipole, annotated as ‘4-2’,
was driven by an amplitude-modulated sinusoid in which thecarrier frequency was 15 Hz and the modulation frequency
was 1 Hz. The current strength to drive the deep dipole was
0.225 mA. The current values of the two dipoles were chosen in order for the magnetic field of the superficial dipole tohave a 16-times stronger intensity than that of the deep dipole.
Namely, the signal-to-interference ratio (SIR), γ, was set to
1/16 in the first experiments. The data were acquired for 2 s
using a sampling frequency of 1 kHz.
The 160-sensor time courses measured when only thesuperficial dipole was active are shown in figure 9(a), andthe time courses measured when only the deep dipole wasactive are shown in figure 9(b). Adaptive beamformer sourcereconstruction was applied to these data sets. The image of thesuperficial dipole is shown in figure 9(c), and the image of the
deep dipole is shown in figure 9(d). The sensor time courses
and the dipole location of the deep source in figures 9(b) and(d) work as the ground truth in the following experiments.
The 160-sensor time courses measured when the superficial and deep dipoles were simultaneously active are shownin figure 10(a). In these sensor data, since the signal from thesuperficial dipole was sixteen times stronger than the signal
from the deep dipole, the sensor time courses were dominatedby the signal from the superficial source. Source reconstruction results from these sensor data are shown in figure 10(c).
Although the sensor data only show the dominated superficial
dipole activity, the reconstruction results show both the superficial and deep dipoles.
We then applied the bDSSP algorithm to detect the signal
from the deep source. We set the local source space at a
1 cm-cubic region whose center was equal to the location of the
deep dipole. The bDSSP algorithm was applied to the sensor

5. Experiments
5.1. Experiments using phantom data
5.1.1. Phantom data measurements. Experiments using datafrom an MEG phantom were performed to test the usefulnessof the proposed algorithm. The phantom used in our experiments is shown in figure 8. In this phantom, dipole sources
consist of isosceles-triangular-shaped coils; these triangularcoils generate magnetic fields expressed by the Sarvas formula[25]. Thus, the coils behave like dipole sources in the spherical homogeneous conductor [29]. A phantom using such triangular coils is called a dry phantom [30]. The triangular coils
installed on the surface of a disc-shaped phantom can be seenin figures 8(a) and (b), in which the squares with a pale color
show the locations of the coils. The center of the disc matches
the sphere origin of the spherical conductor.
A whole-head MEG system with a 160-channel sensor
array [31], installed at Applied Electronics Laboratory,
Kanazawa Institute of Technology, Amaike, Kanazawa,
Japan, was used to measure the phantom data. Figure 8(c)shows how the phantom was installed within the bore of the
MEG sensor helmet. The phantom is equipped with marker
coils, which were used to match between the phantom andsensor coordinates.
5.1.2. Experiments using dipole pair#1. We performed twocases of experiments using different dipole pairs. The firstcase used a pair of dipole sources shown with the annotation‘Dipole pair#1’ in figure 8(a). These dipoles are perpend­
icular to each other and 2 cm apart. The pair was placed near
the parietal-lobe region. The superficial dipole, which is annotated as ‘4-1’, was driven by an 11 Hz sinusoid with a current
11

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 12. Results of experiments using a pair of dipoles denoted ‘Dipole pair#2’; the deep dipole is denoted ‘2-1’, and the superficial

dipole ‘1-2’. (a) The sensor time courses measured when the superficial and deep dipoles were simultaneously active. (b) The bDSSPprocessed sensor time courses. These sensor time courses in (a) and (b) are normalized to each maximum field intensity. The abscissa ofthese figures indicates time in ms. (a) Source image reconstructed using the sensor data in (a). (b) Source image reconstructed using thebDSSP-processed sensor data in (b). The left, middle, and right panels, respectively, show the axial, coronal and sagittal projections of the
3D source reconstruction results.

data in figure 10(a) to extract the signal from the deep source.
The resultant sensor time courses are shown in figure 10(b),and source reconstruction results are shown in figure 10(d).
Comparisons between these results and those in figures 9(b)and (d) demonstrate that the proposed bDSSP algorithm cansuccessfully extract the activity of the deep dipole from thesensor data dominated by a large interference from the superficial source. Note that the difference between the source locations in figures 9(d) and 10(d) is 0.5 cm.

parallel to each other and 4 cm apart, and they were placed near
the temporal-lobe region. The superficial dipole, which is annotated as ‘2-1’, was driven by an 11 Hz sinusoid with a current
strength of 0.71 mA. The deep dipole, annotated as ‘1-2’, wasdriven with the same amplitude-modulated sinusoid as in theprevious experiments with a current strength of 0.5 mA. These
current values were chosen in order for the intensity of the magn­
etic field from the superficial dipole to be eight times stronger
than that from the deep dipole. The reconstruction results of thesuperficial dipole are shown in figure 11(a). The reconstruction
results of the deep dipole are shown in figure 11(b).
The sensor time courses obtained when the superficial
and deep dipoles were simultaneously active are shown in

5.1.3. Experiments using dipole pair#2. We conducted thesame experiments using a pair of dipole sources shown withthe annotation ‘Dipole pair#2’ in figure 8(a). These dipoles are
12

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 13. Source reconstruction results obtained from synthetic sensor data ysy (t) obtained such that ysy (t) = ydeep (t) + αysup (t). Here,ydeep (t) is the sensor data measured when only the deep dipole (dipole‘1-2’) was active, and ysup (t) is the sensor data measured when only
the superficial dipole (dipole‘2-1’) was active. (a) Source reconstruction results when α = 0.04. (b) Source reconstruction results when
α = 0.02. (c) Source reconstruction results when α = 0.01.

figure 12(a). Again, the sensor data were dominated by thesignal from the superficial dipole. Results of source reconstruction applied to these sensor data are shown in figure 12(c).
In these results, since the superficial source is dominated, thedeep dipole can hardly be seen. We applied the bDSSP algorithm for the removal of the interference from the superficial
source activity. The resultant sensor time courses are shownin figure 12(b). These results show that the interference signal
from the superficial dipole was successfully removed. The
source reconstruction results obtained using these bDSSPprocessed sensor data are shown in figure 12(d). Althoughconsiderable blur was introduced and the difference between
the source locations in figures 11(b) and 12(d) is 2.5 cm in thiscase, the deep dipole can still be detected in figure 12(d).
In figure 12(d), the resultant source image contains a large
blur and a fairly large localization error was caused. This isdue to the high spatial correlation between the lead fields of thetwo sources [22]. (Since the two sources are parallel dipoles,
the similarity of their lead fields are considerably high, eventhough they are 4 cm apart.) By making use of the blur in thebDSSP results, we can make a rough estimate on how muchsignal from the superficial source remains in the bDSSP-processed results. Let us denote the sensor data measured when

only the superficial dipole (dipole‘2-1’) was active by ysup (t),and the data measured when only the deep dipole (dipole‘1-2’)
was active by ydeep (t). Then, the synthetic sensor data ysy (t)is computed such that ysy (t) = ydeep (t) + αysup (t) where thepositive value α controlled the mixture ratio. The sourcereconstruction results using ysy (t) with α = 0.04, α = 0.02,
and α = 0.01 are shown in figures 13(a)–(c), respectively.
Let us compare these reconstruction results with thebDSSP interference removal results in figure 12(d), and we
can see that the results in figure 13(b) are closest to those infigure 12(d), This observation leads to the estimation thatthe value of α could be around 0.02 in the bDSSP results infigure 12(d). Since the value of α was equal to 8 in the original sensor data shown in figure 12(a), the attenuation of theinterference signal from the superficial dipole is estimated as
0.02/8  =  1/400. That is, we can estimate that the bDSSP algorithm has reduced the intensity of the superficial interference
by a factor of 400 in this experiment.
5.2. Experiments using real MEG data

The bDSSP algorithm was also applied to somatosensory
MEG data to further demonstrate its usefulness. The data set
13

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 14. (a) Top panel: somatosensory MEG data averaged across 240 trials collected using the 275-channel CTF MEG scanner. Bottom

panel: the sensor time courses of the bDSSP results. (b) Top panel: results of source reconstruction obtained from the sensor data in the top
panel of (a). The blank squares show the boundary of the local source space used when applying the bDSSP algorithm. A single source,labelled by (A), was detected at the primary somatosensory area. Bottom panel: results of source reconstruction obtained from the bDSSPprocessed sensor data shown in the bottom panel of (a). Three sourses, labelled as (B)–(D), are detected. The cross sectional MR images
were chosen as those at the voxel having the maximum intensity. (c) Reconstructed source time courses. The top, second top, third top, andbottom panels respectively show the time courses of voxels indicated by (A)–(D). The all four time courses are normalized to a common
value equal to the maximum value of the time course in the top panel.

image (in the upper panel of figure 14(b)), can be observed.
(They are labelled as (B)–(D).)
The reconstructed source time courses are shown infigure 14(c). These time courses indicate that the three sources
observed in the results in the bottom panel of figure 14(b)were stimulus elicited. The results in the bottom panel of
figure 14(b) show that the two sources labeled as (B) and (D)are located on the posterior bank of the central sulcus and considered as somatosensory activities. The third source, labeledas (C), is located near the supplementary motor area (SMA)
and can also be considered as a somatosensory-related activity.
These results demonstrate the usefulness of the bDSSP algorithm for investigating weak brain activities usually hidden bysimultaneously-existing, much stronger activities.

was collected from a healthy subject using the 275-channel
CTF MEG scanner. Vibrotactile stimuli were applied to a subject’s right index finger for a total of 240 trials, and sensor
time courses were averaged across these trials time-aligned
to the onset of the tactile stimulus. The averaged sensor timecourses are shown in the upper panel in figure 14(a) in whicha peak due to the activity from the hand area in the primarysomatosensory cortex (S1) can be seen 50–60 ms after stimulation. The results of source reconstruction from these sensor
data are shown in the upper panel of figure 14(b). Here, weused the sparse Bayesian (Champagne) source reconstruction
algorithm developed by our group [32, 33]. In these results,
a clear and focussed activity (labelled as (A)) arising fromthe primary somatosensory cortex can be observed. However,no other somatosensory activities that are weaker than theactivity in the primary somatosensory cortex were detected.
We applied the bDSSP algorithm in order to detect weakersomatosensory cortical activity by suppressing the strong S1
activity. A local source space was set so that it covers more
ventral regions of the somatosensory cortex, as indicated by
the rectangles in figure 14(b). The sensor time courses of thebDSSP results are shown in the bottom panel of figure 14(a).
Here, it can be seen that the large peak near 50 ms is removed.
The results of source reconstruction obtained from these timecourses are shown in the bottom panel of figure 14(b). Here,three sources, not observed in the original reconstruction

6. Complementary bDSSP algorithm

14

There can be a complementary way to implement the bDSSP
algorithm. In the complementary version, the prerequisite of
the algorithm is the information on the location of the superficial (interference) sources. The local source space is set ata region including the superficial sources, and beamspacebasis vectors are computed by following exactly the sameprocedure described in section 3. Denoting the basis vectors
u1 , . . . , uP, we can compute the beamspace projector Psup:
as 

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

Figure 16. (a) Results of applying the beamspace SSS method tosimulated sensor measurements in the top panel of figure 3(a). (a)Results obtained with LT set at 1. (b) Results obtained with LT set at
2. (c) Results obtained with LT set at 3. These sensor time courses in(a)–(c) are normalized to each maximum field intensity.

(I − Psup )Bsup = 0.
(34)The equation above indicates that the intersection,rsp(Psup B) ∩ rsp((I − Psup )B), forms a subset of Kdeep, andthe projector onto the intersection, Πisc, can be used as the
projector onto Kdeep. Thus, the signal from the deep source is deep = BΠisc. The algorithm is called theextracted by using B‘complementary bDSSP’, (cbDSSP) algorithm.
We applied the cbDSSP algorithm to the computer-generated data in the top panel of figure 3(a), and the results areshown in figure 15(a). The algorithm effectively removed theinterference in this computer simulation. We then applied
the cbDSSP algorithm to the phantom data in figure 10(a)obtained using dipole pair#1, and the results are shown infigure 15(b). Here, the algorithm fails to remove the interference from the superficial dipole.
Figure 15(c) shows Psup B and (I − Psup )B in these phantom
experiments. The data set Psup B is shown in the upper paneland (I − Psup )B, in the lower panel. Here, the two data sets
are normalized by the same value, (the maximum value of
Psup B). It can be seen in figure 15(c) that although most of Bsup
comp­onents were removed in (I − Psup )B, a small amount of
11 Hz sinusoid, the component from the superficial source,
still remains. This is probably because Bsup is much stronger
than Bdeep in the phantom data in figure 10(a)10.
The failure of the cbDSSP algorithm in the phantom experiments indicates one of its limitations. While the bDSSP algorithm uses the projector Pdeep to suppress the signal from thedeep source, the cbDSSP algorithm uses Psup to suppress thesignal from the superficial source. Since it is generally truethat the signal from the superficial source is much stronger
than the signal from the deep source, the cbDSSP algorithm

Kdeep ⊃ rsp(Psup B) ∩ rsp((I − Psup )B).
(35)

10
As mentioned, Bsup is sixteen times stronger than Bdeep in these phantom
data.

Figure 15. (a) Results of applying the cbDSSP algorithm tothe computer generated data in the top panel of figure 3(a). (b)Results of applying the cbDSSP algorithm to the phantom datain figure 10(a) obtained using dipole pair#1. The sensor timecourses in (a) and (b) are normalized to each maximum value, andthe ordinate of these figures indicates the normalized values of themagnetic field intensity. (c) Two data sets, Psup B and (I − Psup )B,from the phantom data, obtained using dipole pair#1. The data
set Psup B is shown in the upper panel and the data set (I − Psup )Bis in the lower panel. These data sets are normalized to the same
maximum value of Psup B.

Psup = [uP ][
uP ]T . Note that Psup is a projector
u1 , . . . , 
u1 , . . . , onto the pseudo-signal subspace of the superficial source. (Inthe original bDSSP algorithm, the beamspace projector isdenoted Pdeep because it is a projector onto the pseudo signal
subspace of the deep source.)
By multiplying Psup and I − Psup with the data matrix B ,we obtain
Psup B = Bsup + Psup Bdeep + Psup Bε ,(32)(I − Psup )B = (I − Psup )Bdeep + (I − Psup )Bε ,(33)where we use

Therefore, we can finally derive [27],15

K Sekihara et al

J. Neural Eng. 15 (2018) 036026

to be used with the bDSSP algorithm. We will investigate on
this point, and publish the results in a future occasion.
This paper proposes a novel algorithm for selective detection of a deep source by suppressing interference signals fromsuperficial sources in MEG measurements. The proposed algorithm combines the beamspace preprocessing method with thedual signal space projection (DSSP) interference suppression
algorithm. A prerequisite of the proposed algorithm is priorknowledge on the location of the deep source. After presenting
a concise review on the DSSP interference removal algorithm,
this paper introduces a method we call the beamspace DSSP
algorithm. Compared with the previously proposed beamspace
signal space separation (SSS) method, the proposed algorithm iscapable of suppressing much stronger interference from superficial sources. This capability is demonstrated in our comp­uter
simulations as well as experiments using phantom data.
The proposed bDSSP algorithm can be a powerful tool inneuroscience studies of physiological functions of midbrain
structures such as the thalamus, amygdala, hippocampus andthe basal ganglia, and of the cerebellum, because there aremany studies that require accurate localization of physiological and pathophysiological activities in deep brain regions.
One example is non-invasive imaging studies of neural oscillations that examine the computational role of these brain
regions. Other examples include studies on pathophysiological activities from deep brain structures. Such pathologic
activities are implicated in a variety of conditions such as
epilepsy and dementias (hippocampus and medial temporallobe), movement disorders (basal ganglia, cerebellum), andneurodevelopmental disorders (thalamus etc). Accurate
reconstruction of such activities in deep brain regions can lead
to better understanding of these pathological conditions andpotentially contributes to improved clinical managements andtreatment strategies.

must be less robust to various types of errors. This is becausesuch errors may force (I − Psup )Bsup to have non-zero values,
while (I − Pdeep )Bdeep approximates zero.
7. Discussion and summary
The beamspace SSS method [14] has been proposed to attainthe selective detection of components from deep sources by
suppressing components from superficial sources. In the original SSS method, signal components of the data vector yS are
estimated using the expansion:∞ yS (t) =α,m c,m ,(36)
=1 m=−

where c,m are the SSS basis vectors for the internal part andα,m , are called the multipole components, and represent thecoefficients of the expansion [15, 16]. Here, the time notation is omitted from α,m for simplicity. The beamspace SSS
method utilizes the fact that components from superficial
sources are represented more by larger  terms, and comp­onents from deep sources are represented more by smaller 
terms. Thus, estimating the data vector using only lower-order
SSS basis vectors provides an estimated data vector containing more deep-source components and less superficialsource components. This can be achieved in simplest form bytruncating the basis vector expansion, as suggested in [14].
That is, to estimate the deep source activity, the beamspace
SSS algorithm uses
LT 



y
(t)
=
α,m c,m ,
(37)
deep
=1 m=−

where ydeep (t) is the estimated deep source activity. In theequation above, the summation with respect to the index is truncated at LT. We applied the beamspace SSS algorithm
to the simulated sensor data in the top panel of figure 3(a).
The results obtained by truncating  to 1, 2, and 3 are shownin figure 16. The results show that the beamspace SSS canenhance the signal from the deep source. However, considerable amount of signal from the superficial source still remains
even when  is set at 1, the minimum value of . These results
suggest that the beamspace SSS method has only a limited
capability of suppressing superficial interference.
In our computer simulation and our experiments using
phantom data, the spherical homogeneous conductor model
was used. One limitation with this conductor model is thatthe discrepancy between the actual and computed lead fields
becomes large for midbrain and deep brain regions. This
may affect the performance of the beamspace projector Pdeep
defined in equation (23). That is, one key factor for the successof the algorithm is that the relationship (I − Pdeep )Bdeep ≈ 0
holds. However, due to the inaccuracy of the conductor model,(I − Pdeep )Bdeep could have non-zero components, which
may cause removal of the signal magnetic field from a deep
source, as well as removal of the interference from superficial
sources. Therefore, to avoid such signal cancellation, a better
conductor model, such as the realistic head model, may have

ORCID iDs
Kensuke Sekihara

https://orcid.org/0000-0002-2048-3797

References[1] Hämäläinen M, Hari R, Ilmoniemi R J, Knuutila J andLounasmaa O V 1993 Magnetoencephalography-theory,
instrumentation, and applications to noninvasive studies of
the working human brain Rev. Mod. Phys. 65 413–97[2] Hari R and Salmelin R 2012 Magnetoencephalography: fromSQUIDs to neuroscience: neuroimage 20th anniversary
special edition Neuroimage 61 386–96[3] Supek S and Aine C J 2016 Magnetoencephalography (Berlin:
Springer)[4] Attal Y, Bhattacharjee M, Yelnik J, Cottereau B, Lefèvre J,
Okada Y, Bardinet E, Chupin M and Baillet S 2007
Modeling and detecting deep brain activity with MEG &EEG 29th Annual Int. Conf. of the IEEE Engineering in
Medicine and Biology Society (IEEE) pp 4937–40[5] Jerbi K et al 2007 Coherent neural representation of hand
speed in humans revealed by MEG imaging Proc. Natl
Acad. Sci. 104 7676–81
16

K Sekihara et al

J. Neural Eng. 15 (2018) 036026[20] Mosher J C, Lewis P S and Leahy R M 1992 Multiple dipole
modeling and localization from spatio-temporal MEG data
IEEE Trans. Biomed. Eng. 39 541–57[21] Van Veen B D, Van Drongelen W, Yuchtman M and Suzuki A
1997 Localization of brain electrical activity via linearly
constrained minimum variance spatial filtering IEEE Trans.
Biomed. Eng. 44 867–80[22] Sekihara K and Nagarajan S S 2008 Adaptive Spatial
Filters for Electromagnetic Brain Imaging (Berlin:
Springer)[23] Rodríguez-Rivera A, Baryshnikov B V, Van Veen B D and
Wakai R T 2006 MEG and EEG source localization inbeamspace IEEE Trans. Biomed. Eng. 53 430–41[24] Sekihara K, Kawabata Y, Ushio S, Sumiya S, Kawabata S,
Adachi Y and Nagarajan S S 2016 Dual signal subspace
projection (DSSP): a novel algorithm for removing large
interference in biomagnetic measurements J. Neural Eng.
13 036007[25] Sarvas J 1987 Basic mathematical and electromagnetic
concepts of the biomagnetic inverse problem Phys. Med.
Biol. 32 11–22[26] Ipsen I C 2009 Numerical Matrix Analysis: Linear Systems
and Least Squares (Philadelphia: SIAM)[27] Sekihara K and Nagarajan S S 2017 Subspace-based
interference removal methods for a multichannel
biomagnetic sensor array J. Neural Eng. 14 051001[28] Golub G H and Van Loan C F 2012 Matrix Computations vol
3 (Baltimore, MD: Johns Hopkins University Press)[29] Ilmoniemi R, Hämäläinen M and Knuutila J 1985 The forward
and inverse problems in the spherical model Biomagnetism:
Applications and Theory ed H Weinberg et al (Oxford:
Pergamon)[30] Oyama D, Adachi Y, Yumoto M, Hashimoto I and Uehara G
2015 Dry phantom for magnetoencephalography:
configuration, calibration, and contribution J. Neurosci.
Methods 251 24–36[31] Uehara G, Adachi Y, Kawai J, Shimogawara M, Higuchi M,
Haruta Y, Ogata H and Hisashi K 2003 Multi-channel
SQUID systems for biomagnetic measurement IEICE
Trans. Electron. 86 43–54[32] Wipf D P, Owen J P, Attias H T, Sekihara K and
Nagarajan S S 2010 Robust Bayesian estimation of
the location, orientation, and time course of multiple
correlated neural sources using MEG NeuroImage
49 641–55[33] Sekihara K and Nagarajan S S 2015 Electromagnetic Brain
Imaging: a Bayesian Perspective (Berlin: Springer)[6] Tesche C 1996 Non-invasive imaging of neuronal population
dynamics in human thalamus Brain Res. 729 253–8[7] Roux F, Wibral M, Singer W, Aru J and Uhlhaas P J 2013 The
phase of thalamic alpha activity modulates cortical gammaband activity: evidence from resting-state MEG recordings
J. Neurosci. 33 17827–35[8] Lithari C, Moratti S and Weisz N 2015 Thalamocorticalinteractions underlying visual fear conditioning in humans
Human Brain Mapp. 36 4592–603[9] Jha A et al 2016 Functional connectivity of thepedunculopontine nucleus and surrounding region inparkinson’s disease Cerebral Cortex 27 54–67[10] Luo Q, Holroyd T, Mitchell D, Yu H, Cheng X, Hodgkinson C,
Chen G, McCaffrey D, Goldman D and Blair R J 2017
Heightened amygdala responsiveness in s-carriers of 5-httlpr
genetic polymorphism reflects enhanced cortical rather thansubcortical inputs: an MEG study Human Brain Mapp.
38 4313–21[11] Garvert M M, Friston K J, Dolan R J and Garrido M I
2014 Subcortical amygdala pathways enable rapid faceprocessing Neuroimage 102 309–16[12] Hashimoto I, Kimura T, Tanosaki M, Iguchi Y and Sekihara K
2003 Muscle afferent inputs from the hand activate humancerebellum sequentially through parallel and climbing fiber
systems Clin. Neurophysiol. 114 2107–17[13] Herrojo Ruiz M, Maess B, Altenmüller E, Curio G and
Nikulin V V 2017 Cingulate and cerebellar beta oscillations
are engaged in the acquisition of auditory-motor sequences
Human Brain Mapp. 38 5161–79[14] Özkurt T E, Sun M and Sclabassi R J 2008 Decompositionof magnetoencephalographic data into components
corresponding to deep and superficial sources IEEE Trans.
Biomed. Eng. 55 1716–27[15] Taulu S and Kajola M 2005 Presentation of electromagnetic
multichannel data: the signal space separation method
J. Appl. Phys. 97 124905[16] Taulu S, Simola J and Kajola M 2005 Applications of thesignal space separation method IEEE Trans. Signal Process.
53 3359–72[17] Buckley K M and Xu X-L 1990 Spatial-spectrum estimation
in a location sector IEEE Trans. Acoust. Speech Signal
Process. 38 1842–52[18] Lee H B and Wengrovitz M S 1990 Resolution threshold of
beamspace MUSIC for two closely spaced emitters IEEE
Trans. Acoust. Speech Signal Process. 38 1545–59[19] Schmidt R O 1981 A signal subspace approach to multiple
emitter location and spectral estimation PhD Thesis
Stanford University Stanford, CA

17