NeuroImage 188 (2019) 161–170

Contents lists available at ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/neuroimage

Evaluation of a dual signal subspace projection algorithm inmagnetoencephalographic recordings from patients with intractableepilepsy and vagus nerve stimulators
Chang Cai a, Jiajing Xu a, Jayabal Velmurugan a, b, c, d, Robert Knowlton g, Kensuke Sekihara e, f,
Srikantan S. Nagarajan a, *, Heidi Kirsch a, g
aDepartment of Radiology and Biomedical Imaging, University of California, San Francisco, CA, 94143-0628, USA
Department of Clinical Neurosciences, National Institute of Mental Health and Neurosciences, Bangalore, Indiac
MEG Research Center, National Institute of Mental Health and Neurosciences, Bangalore, Indiad
Department of Neurology, National Institute of Mental Health and Neurosciences, Bangalore, Indiae
Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45 Yushima, Bunkyo-ku, Tokyo 113-8519, Japanf
Signal Analysis Inc., Hachioji, Tokyo, Japang
Department of Neurology, University of California, San Francisco, CA 94143-0628, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Brain mapping
Magnetoencephalography
DSSP
Intractable epilepsy
Vagus nerve stimulators

Magnetoencephalography (MEG) data is subject to many sources of environmental noise, and interference
rejection is a necessary step in the processing of MEG data. Large amplitude interference caused by sources near
the brain have been common in clinical settings and are difﬁcult to reject. Artifact from vagal nerve stimulators(VNS) is a prototypical example. In this study, we describe a novel MEG interference rejection algorithm called
dual signal subspace projection (DSSP), and evaluate its performance in clinical MEG data from people withepilepsy and implanted VNS. The performance of DSSP was evaluated in a retrospective cohort study of patientswith epilepsy and VNS who had MEG scans for source localization of interictal epileptiform discharges. DSSP wasapplied to the MEG data and compared with benchmark for performance. We evaluated the clinical impact ofinterference rejection based on human expert detection and estimation of the location and time-course ofinterictal spikes, using an empirical Bayesian source reconstruction algorithm (Champagne). Clinical recordings,
after DSSP processing, became more readable and a greater number of interictal epileptic spikes could be clearly
identiﬁed. Source localization results of interictal spikes also signiﬁcantly improved from those achieved before
DSSP processing, including meaningful estimates of activity time courses. Therefore, DSSP is a valuable novel
interference rejection algorithm that can be successfully deployed for the removal of strong artifacts and interferences in MEG.

1. Introduction
From the time of its ﬁrst introduction, magnetoencephalography(MEG) has been used to map functional brain activity noninvasively withgood spatial and excellent temporal resolution, and thus to offer valuable
information for use in clinical neurology and basic neuroscience. However, MEG has suffered from an important shortcoming: it is prone tocontamination from signals other than the signals of interest - including
inevitable non-biological sources like power lines and trains, and biological sources outside of the brain like the heart. Though most of thisinterference is of similar magnitude to brain activity, some of it is highamplitude and needs special attention - including artifact from dental
work, and especially interference from vagal nerve stimulators (VNS),
relatively common in people with intractable epilepsy, that makes it very
difﬁcult for us to see and then to model activity of interest (Nagarajan
et al., 2006a; Taulu and Simola, 2006; Sekihara et al., 2016).
A variety of methods have been used to minimize artifact in magnetoencephalographic recordings with varying degrees of success. Averaging responses over trials is one method commonly used; this takes
advantage of the idea that interference in different trials is statisticallyindependent, whereas evoked signals are not. However, this method
requires a large number of trials, and evoked signals must be relatively* Corresponding author.
E-mail addresses: sri@ucsf.edu (S.S. Nagarajan), Heidi.Kirsch@ucsf.edu (H. Kirsch).
https://doi.org/10.1016/j.neuroimage.2018.11.025
Received 30 July 2018; Received in revised form 22 October 2018; Accepted 16 November 2018
Available online 29 November 20181053-8119/© 2018 Published by Elsevier Inc.

C. Cai et al.

NeuroImage 188 (2019) 161–170

Fig. 1. Schematic showing the processing steps of DSSP.

evoked responses from primary sensory cortices, thus diminishing theusefulness of MEG for these patients and, thereby, their hope for recovery(Sekihara et al., 2016). Therefore, developing and testing algorithms forinterference rejection in MEG data is important, especially new algorithms that speciﬁcally address the kind of interference that is not wellhandled by currently available options but that is clinically important(e.g. VNS implant interference). Ideally such an algorithm would berobust and broadly capable of rejection of as many types of interference
as possible. Given that many source localization platforms include leadﬁelds, it would be ideal to offer a tool that is also based on lead ﬁelds.
Right now options are restricted to speciﬁc hardware platforms. For
example, the temporally extended signal space separation method (tSSS)
developed by one MEG manufacturer offers a potential solution (Taulu
and Hari, 2009) but this tool has only been demonstrated for the Elekta
platform and has not been shown for other platforms. In contrast, here weshow a MEG hardware platform independent algorithm for large interference rejection.
Dual signal subspace projection (DSSP) is a newly proposed algorithm
for removal of large interference in biomagnetic measurements, and has
the potential to handle many different kinds of interference (Sekihara
et al., 2016). DSSP is based on the fact that MEG signal has both spatial
and temporal properties. This allows us to deﬁne a signal subspace in thespace domain, and another signal subspace in the time domain. We assume that the interference signal is present all the time across the whole
signal subspace, either inside or outside the spatial-domain signal subspace, or both. In contrast, activity from the brain is presumed to exist
only inside the spatial domain signal subspace. The DSSP algorithm ﬁrst
projects the columns of the measured data matrix onto the inside andoutside of the spatial-domain signal subspace, creating two ‘projected’
data matrices. The intersection of the row spans of these two ‘projected’matrices is then taken to be an estimate of the time-domain interferencesubspace, and artifact removal is carried out on the basis of this estimatedinterference subspace. Details of the DSSP algorithm have been publishedrecently, but the performance of DSSP in assisting the identiﬁcation andlocalization of epileptiform discharges has not been determined. In thispaper we evaluate its ability in these arenas, using subject speciﬁc leadﬁelds and selecting parameters, exploring its capability to handle various

similar and robust (Nagarajan et al., 2007a). Filtering is another widely
applied method, but requires prior knowledge about the interference.
Recently, data-driven approaches such as principal component analysis(PCA), independent component analysis (ICA) have been popular.
However, these methods ask users to make subjective choices during
application (e.g. choice of threshold in PCA and of interference component in ICA), and the methods cannot exploit pre-/poststimulus partitioning of the data (Ossadtchi et al., 2004; Nagarajan et al., 2006b). Joint
decorrelation is another method commonly supposed to be robust to
many types of interference problems, but its use requires the design ofdifferent bias ﬁlters for different interference types, and thus to someextent requires prior knowledge of the interference (de Cheveigne andParra, 2014). Algorithms based on statistical properties of the interference are a class of automated interference algorithm method hailed asboth reliable and robust. The partitioned factor analysis (PFA) algorithm(Nagarajan et al., 2006a, 2007b; Zumer et al., 2008) is implemented byobtaining a probabilistic model from the data distributions in the
pre-stimulus period (when the interference exists) and the post-stimulus
period (when both interference and true signal exist), and then inferringmodel parameters from these distributions. This method handles mosttypes of interference well, but since it relies on the availability of separate
measurements that capture the statistical properties of the interference,
its use is limited to situations where such separate measurements areappropriate, and it is not effective for removing overlapped interference(Nagarajan et al., 2007a). Also, these algorithms may not be effective forinterference of extremely large magnitude relative to the signals being
estimated, which is often seen in MEG data in patients with VNS
implants.
Artifacts of signiﬁcant magnitude are not rare in MEG recordings, andresolving MEG data from distorted recordings is often of great clinicalsigniﬁcance. Particularly in the case of people with intractable epilepsy
who have received VNS implants and have continued refractory focal
onset seizures, MEG studies are an important part of the evaluation forand the planning of resective surgery. Without interference rejection,
MEG data in many people with VNS implants will be completely distorted
by signiﬁcant artifact from the stimulator and the lead-wires, making it
extremely difﬁcult to see interictal epileptiform activity or stimulus

162

Durationof
Epilepsy

MR abnormality

Ictal EEG

Interictal EEG

PET CT

Presumed EZ

Interictal MEG
spikesNum ofspikesbeforeDSSP

Num ofspikesafter
DSSP

Max ﬁeldstrength
BeforeDSSP (pT)Max ﬁeldstrengthafter DSSP
(pT)

1

22

18

Left lateral frontal lobe
cortical dysplasia

Poorly localized;
left frontocentral
regionLeft frontocentral
spikes or polyspikes

Left frontal
onset

Left frontotemporal

0

39

70.2

2.3

2

25

20

Primary read as
normal, secondary
read as bilateralposterior pachygyria

Seizures arising
independently
from eachhemisphere; poorly
localized

Independent
bitemporal spikes;
generalizedparoxysmal fast
actvity

Normal but PET
fusion with MRI
corresponding
hypometabolism
Negative

Unknown to
date

44

33

5.3

2.8

3

44

32

UnremarkableNot available

N/A

Unknown to
date

107

168

148.3

4.8

4

22

20

Encephalomalacia ofthe left temporal lobe,
volume loss of left
hippocampus

Left parietal region

N/A

Left temporoparietaloccipital

Posterior medial left
temporal lobe

3

100

48.3

2.5

5

17

Not
available

Left hipopcampal
atrophy, left
hemispheric corticaldysplasia

Left hemisphere
onset

Right temporal sharp
waves, generalized
spike and polyspike
discharges
Left TIRDA, frequent
broad spikes over left
temporo-parietal
region, occasional left
anterior temporal
predominanceIntermittent left
frontotemporaldischarges

Bilateral slow waves
model bilaterally inthe suprasylvian
frontal and infrasylvian temporal
lobes
Right temporal,
right frontal

Lefthemisphere,
probable left
temporal lobe

Left temporal region

47

143

60.5

2.1

6

38

17

Left parietal, left
temporalIndependent
bilateral
frontotemporal

Independent right andleft temporaldischarges

Hypometabolism ofleft temporal lobe, leftparietal lobe, left
posterior occiptal
lobe
Bilateral temporalhypometabolism

None

5

14

21.1

2.0

7

31

25

Unremarkable

Bilateral central/paracentral regions

N/A

54

61.9

1.6

37

19

Left anterior
temporal; also rareright temporal spikes

Bilateral temporalhypometabolism

Left mesial
temporal

24

32

32.0

1.8

9

26

Since
young

T2/FLAIR
hyperintensity andatrophy of bilateral
temporal lobes, L > R
Expected changes frommedial left frontal lobe
corticectomy;
otherwise
unremarkable

Right cingulate
gyrus; L > R
perirolandic regions
Right temporal
spikes, rare left
temporal spikes

35

8

Vertex spike
followed by diffuse
fast activity
Left frontotemporal

Frontal or
temporal;laterality
unknown
Unknown to
date

No interictal

N/A

Frontal;lateralizationunclear but
more likely left

None

43

72

56.6

1.1

10

28

27

Bifrontal sharp waves,left frontal spikes

Increased metabolic
activity in high left
posterior pariental
sulcusUnknown to
date

Right suprasylvian
frontal lobe

3

53

116.7

2.6

Left parietal cavernous
malformation

Poorly localized
and lateralized;
some withpreceding leftparasagittal sharp
waves
Suggestive offrontal onset but
poorly lateralized

Notes

Prior
posteriortemporalresection

NeuroImage 188 (2019) 161–170

Age

163

ID

C. Cai et al.

Table 1
Clinical characteristics of ten subjects. Note that some patients were referred from outside institutions and thus their information was limited to that available at the time of the MEG scan.

C. Cai et al.

NeuroImage 188 (2019) 161–170

Fig. 2. Each panel shows the power spectral density (PSD) for each channel before and after either DSSP (red dots) or ANC (blue dots), for each of the ten subjects inthis study. Points along the diagonal line in each panel indicates lack of interference rejection.

interest are assumed to spatially arise from within the brain. This spatialsubspace is therefore deﬁned by the lead-ﬁeld matrix describing themagnetic ﬁeld distributions arising from activity from all sources within
the brain. The lead-ﬁeld matrix is taken to construct a projector whichallows to split the measured data into two components Bin ﬁelds given by
PS B which stem from brain sources, and Bout ﬁelds given by ðI  PS ÞB
which are orthogonal and do not arise from brain sources. The second
deﬁnition is based on temporal considerations. The temporal representations of Bin and Bout may not be orthogonal to each other and it is likelythat strong interference components that are present in both spatial
components Bin and Bout . These interference components can be identiﬁed by constructing a temporal subspace projector which is orthogonal tothe intersection of the temporal subspaces of Bin and Bout . The ﬁnal
cleaned DSSP data is obtained by projecting the original signal throughthis temporal subspace operator. The difference between the DSSP and

artifacts as part of processing of clinical datasets. In particular, we will
evaluate whether it will be helpful to solve the problem of spike detection
in patients with VNS implants. We will also test its ability to improvesource localization of spikes using Champagne, an empirical Bayesiansource reconstruction algorithm described previously (Julia P Owen
et al., 2012).
2. Method
2.1. The DSSP algorithm
This section introduces the DSSP algorithm brieﬂy; Fig. 1 shows thesteps of the DSSP algorithm and details of the derivation can be found in
Appendix A. In brief, the DSSP algorithm relies on two subspace deﬁnitions. The ﬁrst deﬁnition is based on spatial considerations. Sources of164

C. Cai et al.

NeuroImage 188 (2019) 161–170

international 10–20 system that includes subtemporal electrodes. Thirty
to forty minutes of spontaneous data were obtained in intervals of
10–15 min with the patient asleep and awake. The position of the patient's head in the dewar relative to the MEG sensors was determinedusing indicator coils before and after each recording interval to verify
adequate sampling of the entire ﬁeld. The data were then bandpass
ﬁltered ofﬂine, initially at 1–70 Hz. More details of the recording
methods have been previously described (Nagarajan et al., 2007a). As
artifact commonly distorted MEG recordings from the patients with VNS
implants, in order to enable for visual analysis and dipole ﬁtting of raw
data, additional bandpass ﬁlters (typically 10–70 Hz or 20–70 Hz) wereapplied as needed during analysis of MEG data. After the application of
DSSP for artifact removal, all data were bandpass ﬁltered at 1–70 Hz.
2.4. Epileptic spike analysis
Spikes were visually identiﬁed by experts - a certiﬁed EEG technologist (MM) and a clinical neurophysiologist (JV) - and their results wereconﬁrmed by a board-certiﬁed clinical neurophysiologist and epileptologist (HEK). To ensure bias-free analysis of data, all reviewers wereblinded to the identity of the subject and their clinical data during spike
identiﬁcation and localization procedures. EEG spikes were identiﬁed
based on the criteria deﬁned by the International Federation of Clinical
Neurophysiology (IFCN) (Deuschl Eisenet al., 1999) and the ACMEGS (I
Bagic et al., 2011) for EEG epileptiform discharges. MEG spikes werechosen for analysis based on duration (< 80ms), morphology, ﬁeld map,and lack of associated artifact. The onset of each spike was marked as therising deﬂection of the ﬁrst sharp negativity from the baseline andequivalent current dipoles were ﬁt using commercial software provided
by CTF Systems (VSM MedTech, Port Coquitlam, British Columbia). Only
localized spikes with a goodness of ﬁt higher than 90% were accepted.
Co-registration of dipoles to MRI scans was performed using ﬁducials(nasion and preauricular points) to produce magnetic source images(MSI) of dipoles superimposed on anatomic images. The ﬁtted spike dipoles were then inspected and validated according to their location.
Simultaneous EEG during MEG was used to deﬁne and conﬁrm spikes on
MEG, ascertaining that a signal was not an artifact or another physiologic
feature, and also to identify spikes when MEG recordings were heavily
contaminated by VNS artifact (ie when MEG data were signiﬁcantly
distorted, spike identiﬁcation relied heavily on EEG).

Fig. 3. Averaged power spectral density (PSD) across all channels and subjects,
before and after DSSP and ANC.

Fig. 4. Maximum absolute value of the magnetic ﬁelds across channels and time
before and after DSSP for all ten subjects. Note that the y-axis is plotted in log
units showing one order of magnitude reduction by DSSP.

2.5. DSSP performance evaluation
In the absence of ground truth data, it is somewhat difﬁcult to evaluate performance in absolute terms. We ﬁrst examine the average power
spectrum of the sensor data before and after DSSP processing. To get asense of DSSP algorithm performance we compare power spectram obtained following DSSP with those obtained following another interference removal method (Adaptive Noise Canceling, or ANC). Traditionally,
ANC makes use of data from reference sensors that collect data containing only the interference but not the signal of interest (Widrow et al.,
1975; Adachi et al., 2001). However, in the case of data with VNS
implant artifacts reference sensor data is not available. Therefore, as an
approximation, for this study we deﬁne the reference sensor as the MEG
channel with the highest power where artifact dominates, and can bedeemed to be an “interference only” reference channel. ANC also uses theidea of subspace projection, where the reference sensor is used to create
the temporal span of interference, then data from each sensor is thenprojected onto the subspace orthogonal to the span of interference. This
subspace projection ensures that the contributions from interference toeach sensor is eliminated, leaving intact sources of interest. We also
compare the interference rejection performance of DSSP and of ANC, by
comparing the input and output signal power in each sensor of the array.
After DSSP implementation, cleaned MEG recordings were analyzed
as described above (without additional band pass ﬁltering) by three individuals with expertise in interictal spike detection (MM, JV, HEK) who

the tSSS algorithm is the way of splitting the data into spatial componentsof Bin and Bout . DSSP uses the lead-ﬁeld to perform this whereas the tSSS
uses spherical harmonic expansions of the magnetic ﬁeld data.
2.2. Subjects
We selected 10 epilepsy patients with VNS who underwent a clinical
MEG study as part of epilepsy surgery evaluation at the University ofCalifornia, San Francisco (UCSF) Biomagnetic Imaging Laboratory (BIL)
between November 24th, 2004 and May 6th, 2016. Prior to MEG, allpatients had high-resolution epilepsy protocol 3T T1-MRI scans for coregistration of localized spikes. Table 1 summarizes clinical characteristicsof these subjects.
2.3. MEG recordings
Simultaneous EEG and MEG recordings were performed inside amagnetically shielded room with a 275 channel whole-head axial
gradiometer system (VSM MedTech, Port Coquitlam, British Columbia).
MEG data were recorded from each patient in a passband of 0–75 Hz
using a CTF 275 channel whole cortex MEG helmet while simultaneous
twenty-one channel scalp EEG data were recorded using a modiﬁed165

C. Cai et al.

NeuroImage 188 (2019) 161–170

Fig. 5. A representative case showing the effect of the application of DSSP (Subject 5 from Table 1). (a) EEG epoch corresponding to MEG epoch (selected channels)(b) Raw MEG recordings (selected channels) (c) DSSP-processed MEG data. The red line marks a spike not identiﬁed in the raw data but seen in the DSSP-processed
data, the green line in (a) and (c) mark a different spike. (d) Sensors topographical maps for EEG and MEG at the red line time-point corresponding to the spike (before
and after DSSP). (e) Spike localization using Champagne on DSSP-processed data. (f) Normalized intensity for the spike of interest reconstructed through Champagne.

Fig. 6. Another representative case (Subject 4 from Table 1). Legends are identical to the previous ﬁgure.

each case studied, ten representative spikes seen well on EEG but poorly
on MEG in the unprocessed recordings were selected for analysis. For
each of the spike selections, Champagne was run for 300 ms (i.e. using
180 data points) on truncated MEG epochs each centered on the selected
spike; the required noise estimate input for Champagne was obtained byrunning VBFA on the 1s of MEG preceding the truncated epoch. The spike
source reconstruction map obtained after implementing Champagne wasused to judge the performance of DSSP: we observed the activation valueof the localized spike activity and the recovered activity time-series, andcompared these with standard clinical spike mapping as described above,
as well as with correlative clinical data. Across all subjects, we evaluate
the consistency of spike localization results with Champagne following

were blinded to the results of the initial (pre-DSSP) analysis. The results
were then compared with the original analysis and included quantiﬁcation of the number of spikes identiﬁed and localized, and concordancewith other clinical information (EEG, semiology, MRI lesion if present).
We also use t-tests to evaluate whether the number of spikes identiﬁed
and localized were different before and after DSSP.
Finally, DSSP was integrated into a newer source localization pipeline: a united Bayesian framework for MEG/EEG source imaging that
includes Variational Bayes Factor Analysis (VBFA) for noise approximation and a Sparse Bayesian Algorithm (Champagne) for source localization (Wipf and Nagarajan, 2009; Wipf et al., 2010), to see whether
localization improved upon the incorporation of DSSP algorithm. For

166

C. Cai et al.

NeuroImage 188 (2019) 161–170

reduction in PSD after ANC processing. Upon examination of the PSD as afunction of frequency, it is found that the PSD decreases after DSSP and
ANC, occur mainly in the low frequency bands from 0 to 30 Hz, withnegligible reductions above 30 Hz. However, the reductions in low frequencies are much greater for DSSP when compared to ANC (Fig. 3).
Fig. 4 shows maximum absolute value of the magnetic ﬁeld strength
across channels and time before and after DSSP for each of the ten subjects. Artifact levels were consistently reduced by an order of magnitude,
before DSSP the artifact levels were 62.1 pT (62.1  13.55) and werereduced to 2.4 pT (2.4  0.32) after DSSP processing (Fig. 4). This
reduction of 96.13% was highly statistically signiﬁcant (two-tailed
paired t-test, t-value ¼ 4.40777, p-value ¼ .00034).

3.2. Expert detection and localization of interictal spikes
After the application of the DSSP algorithm to MEG data from people
with intractable epilepsy and VNS, both spikes that were well seen on
EEG, and those seen primarily on MEG could be visually identiﬁed byexperts from MEG background at a high rate. Fig. 5 shows EEG data (a)and MEG data (b) from a patient with intractable epilepsy and a VNS
implant. It features VNS artifact with partial periodicity, low frequency
and high amplitude. After DSSP is applied, as shown in (c), this periodic
feature is greatly diminished and the background looks similar to that ofmost people with epilepsy who do not have VNS. Fig. 5(d) shows sensors
topographical maps for MEG after DSSP and EEG maps. Fig. 6 shows
similar data for another patient. Fig. 7 shows the average number of
spikes that could be identiﬁed by visual inspection of the MEG and theaverage number of spikes that could be localized by topographical inspection before and after DSSP. Over twice as many spikes could be
identiﬁed after DSSP when compared to before DSSP (p ¼ 0:037), andover four times as many spikes could be localized after DSSP whencompared to before DSSP (p ¼ 0:007).
Source localization of identiﬁed interictal spikes in MEG data before
DSSP with the Champagne algorithm resulted in localization failure innine out of the ten cases; in all these cases either no strong activation
could be found, or the activity was localized to unusual positions (e.g.
near or outside of the skull), and the activity time-series recovered did
not resemble that of a spike. However, localization results with Cham-

Fig. 7. The averaged number of spikes that could be identiﬁed by visual inspection and localized by topographical inspection before and after the application of DSSP for ten subjects with standard error bars.

DSSP processing of data with the presumed epileptogenic zone (based on
lesion, ictal EEG, semiology and other clinical data) using a chi-squaretest.
3. Results
3.1. DSSP vs ANC
We compare the power spectral densities (PSD) before and after
artifact rejection with DSSP and ANC algorithms. Fig. 2 shows the PSD for
all channels for all ten subjects included in this study. A diagonal line ineach panel indicates lack of interference rejection. Although both algorithms show reductions in PSD after processing, in all subjects DSSP
shows greater reduction in all channels regardless of the power level ineach channel (most red dots are consistently below the diagonal). In
contrast, the performance of the ANC algorithm is neither consistentacross subjects nor across channels, some channels show negligible

Fig. 8. Source localization results for all ten subjects using Champagne after DSSP.
167

C. Cai et al.

NeuroImage 188 (2019) 161–170

included all patients who had VNS implants in our study. Second, this is aretrospective and non-randomized study. Currently, for new cases with
VNS artifact, we are undertaking a prospective study applying DSSP prior
to initial analysis, and results of this prospective study will be publishedin the future. Finally, although better localization was achieved with theaddition of DSSP to Champagne, whether this truly improved outcomes
as a results of epileptogenic zone mapping is unknown. Additional information about the true epileptogenic zone (e.g. from follow-up after
resective surgery) will be needed to make such judgments. Therefore, we
can only conclude that DSSP helps achieve spike mapping, but cannot
evaluate the impact of localization accuracy on outcome. We are collecting follow-up information from patients who went through surgery
after MEG recording, and in future this information could be used as agold standard to judge the performance of DSSP.
Indeed, DSSP will not perform good artifact rejection if the sources
are located close to the brain and if the artifacts are correlated with brain
activity. The DSSP algorithm is suited to remove interferences fromsources located fairly close to the boundary of the source space. The eyeblink and cardiac artifacts encountered in MEG measurements are also
such interferences.

pagne on MEG data after DSSP resulted in marked improvement.
Figs. 5(e) and 6(e) show source localization, and Figs. 5(f) and 6(f) showrespective Champagne time series for the spike of interest in the twocases previously discussed. In Fig. 8, all cases could be localized correctly
after running DSSP. A Chi-Square test comparing the number of localizations before and after DSSP that are consistent or inconsistent with thepresumed epileptogenic zone demonstrates that Champagne localizations were clearer and more consistent after DSSP (χ 2 ¼ 16:364 and p ¼
0:0001).
In order to quantify the accuracy of spatial localization of spikes with
the Champagne algorithm after DSSP, we evaluated the cluster radius forten representative spikes from each subject (to prevent biases for subjects
with a greater number of spikes). The averaged radius of the spatial
clusters of interictal spikes after DSSP is 13.3 mm (SEM ¼ 1.3) (note thatthe resolution of the forward model is 8 mm), thus indicating the highaccuracy and precision of spike localization with Champagne after DSSP.
4. Discussion
In this study, DSSP was evaluated using typical clinical data frompeople with epilepsy and VNS, showing its potential to diminish the inﬂuence of interference. DSSP processing of MEG data enabled better visual identiﬁcation of spikes, making meaningful the MEG recordings thatwere contaminated and previously of limited value. Finally, when integrated with the Champagne source reconstruction algorithm, DSSP did
help to achieve more reasonable spike localizations and meaningful
recovered spike activity time series. The successful rejection of VNS
artifact using DSSP should therefore improve treatment (including surgical planning for resection or other localized therapies) of people withintractable epilepsy and VNS. Given these results, as we gain further
experience with DSSP, its potential use in the setting of other interference
types can also be explored.
Analysis of the power spectrum of sensor data revealed that DSSP
resulted in greater reductions in lower frequency bands when compared
to an adaptive noise cancelation benchmark. However, in future a moredetailed comparison of performance with other benchmark algorithms
including the tSSS algorithm is warranted (Taulu and Hari, 2009).
There are several limitations to this study. First, small sample size
could potentially limit the generalizability of our results. However, we

5. Conclusion
In short, DSSP is a novel interference rejection algorithm worth
exploration. The retrospective clinical study has shown its potential todeal with high amplitude, periodic interference currently not handled
well by other algorithms. DSSP helped to recover distorted MEG recordings from people with intractable epilepsy and VNS implants, making epileptic spike identiﬁcation easier and spike mapping better. The
speciﬁcity of this improved spike mapping is still unknown.
Acknowledgment
The authors would like to thank Susanne Honma, Danielle Mizuiri
and Anne Findlay for collecting much of the MEG data in the Biomagnetic
Imaging Laboratory. We would also like to thank Mary Mantle for blinded analysis of some of the MEG data. This work was supported by NIH
grants R01EB022717, R01DC013979, UCOP MRP-17-454755, and a giftfrom Ricoh Company Ltd..

Appendix A. Derivation of DSSPAppendix A.1. Data modelThis section brieﬂy describes the DSSP algorithm. A full explanation of the algorithm is available in (Sekihara et al., 2016). Also, a detailed
explanation of the DSSP algorithm in the context of the time-domain signal subspace can be found in (Sekihara and Nagarajan, 2017). Let us deﬁne themeasurement of the m-th sensor at time t as ym ðtÞ. The measurement from 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 the superscript T indicates the matrix transpose. 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 ofthe m-th sensor are respectively denoted as lxm ðrÞ, lym ðrÞ, and lzm ðrÞ. Let us deﬁne an M  3 matrix LðrÞ whose m-th row is equal to a 1  3 row vector ½lxm ðrÞ;lym ðrÞ; lzm ðrÞ. This matrix LðrÞ, referred to as the lead ﬁeld matrix, represents the sensitivity of the sensor array at r.
The DSSP algorithm was proposed in order to remove interfering magnetic ﬁelds overlapped onto signal magnetic ﬁelds. The algorithm assumes thedata model:yðtÞ ¼ y S ðtÞ þ yI ðtÞ þ ε;

(A.1)

where yS ðtÞ, (called the signal vector), represents the signal of interest, y I ðtÞ, (called the interference vector), represents the interference magnetic ﬁeld,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 deﬁned as: B ¼ ½yðt1 Þ; …; yðtK Þ The signal matrix is deﬁned as BS ¼ ½yS ðt1 Þ; …;
y S ðtK Þ; and the interference matrix as BI ¼ ½y I ðt1 Þ; …; y I ðtK Þ Then, the data model in Eq. (A.1) is expressed in a matrix form as:
B ¼ BS þ BI þ Bε ;

(A.2)

168

C. Cai et al.

NeuroImage 188 (2019) 161–170

where Bε is the noise matrix whose j-th column is equal to the noise vector ε at time tj .
Appendix A.2. Pseudo-signal subspace projector
The dual signal space projection (DSSP) algorithm assumes that the interference sources are located outside the source space which indicates aregion in which signal sources can exist. The DSSP algorithm uses the so-called pseudo-signal subspace projector, and to derive it, voxels are deﬁned
over the source space, in which the voxel locations are denoted r 1 ; …; r N . The augmented leadﬁeld matrix over these voxel locations is deﬁned as
F ¼ ½Lðr 1 Þ; …; Lðr N Þ;(A.3)

and the pseudo-signal subspace E S is deﬁned such thatE S ¼ cspðFÞ;(A.4)

where the notation cspðXÞ indicates the column space of a matrix X. If the voxel interval is sufﬁciently small and voxel discretization errors are
negligible, we have the relationship E S ⊃E S where E S 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 pseudo-signal subspace. To do so, we compute the singular value decomposition of F:
F¼

M
X
j¼1

λj ej f Tj ;(A.5)

where ej and f j are left and right singular vectors. In Eq. (A.5), we assume the relationship M < N, and the singular values are numbered in a decreasingorder. 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 the pseudo-signal subspace E S . Thus, the projector onto E S is obtained using
PS ¼ ½e1 ; …; eτ ½e1 ; …; eτ T :(A.6)

Note that ðI  PS ÞyS ðtÞ ¼ ðI  PS ÞBS ¼ 0 holds.
Appendix A.3. DSSP algorithm
The DSSP algorithm applies PS and I  PS to the data matrix B to create two kinds of data matrices:
PS B ¼ BS þ PS BI þ PS Bε ;(A.7)

ðI  PS ÞB ¼ ðI  PS ÞBI þ ðI  PS ÞBε :(A.8)

Let us use the notation rspðXÞ to indicate the row space of a matrix X. Then, the relationships, rspðPS BI Þ ¼ K I , rspððI  PS ÞBI Þ ¼ K I , and rspðBS Þ ¼
K S hold, where K S and K I respectively indicate the time-domain signal and interference subspaces. According to arguments in (Sekihara andNagarajan, 2017), we can ﬁnally derive the relationship:
K I ⊃rspðPS BÞ \ rspððI  PS ÞBÞ:(A.9)

The equation above shows that the intersection between rspðPS BÞ and rspððI  PS ÞBÞ forms a subset of the interference subspace K I . The basis
vectors of the intersection can be derived using the algorithm described in (Gene et al., 2012). Once the orthonormal basis vectors of the intersection
ψ 1 ; …; ψ r are obtained, we can compute the projector onto the intersection ΠI such that
Π I ¼ ½ψ 1 ; …; ψ r ½ψ 1 ; …; ψ r T :(A.10)

Using this ΠI as the projector onto the (time-domain) interference subspace K I , the interference removal is achieved and the signal matrix isestimated by the time-domain signal space projection (Sekihara and Nagarajan, 2017), which isb S ¼ BðI  ΠI Þ ¼ B I  ½ψ 1 ; …; ψ r ½ψ 1 ; …; ψ r T :B

(A.11)The method of removing the interference in a manner described above is called dual signal space projection (DSSP). Note that since the basis vectorsof the intersection, ψ 1 ;…; ψ r , span only a subset of the interference subspace K I , this method cannot perfectly remove interferences. However, when theintersection rspðPS BÞ \ rspððI  PS ÞBÞ is a reasonable approximation of K I , interferences can effectively be removed by the DSSP algorithm.

Deuschl, Günther, Eisen, Andrew, et al., 1999. Recommendations for the Practice ofClinical Neurophysiology: Guidelines of the International Federation of Clinical
Neurophysiology.
Gene, H., Golub, Van Loan, Charles F., 2012. In: Matrix Computations, Volume 3. The
Johns Hopkins University Press.
I Bagic, Anto, Knowlton, Robert C., Rose, Douglas F., Ebersole, John S., ACMEGS Clinical
Practice Guideline (CPG) Committee, et al., 2011. American clinical

References
Adachi, Y., Shimogawara, M., Higuchi, M., Haruta, Y., Ochiai, M., 2001. Reduction of
non-periodic environmental magnetic noise in MEG measurement by continuously
adjusted least squares method. IEEE Trans. Appl. Supercond. 11 (1), 669–672.
de Cheveigne, Alain, Parra, Lucas C., 2014. Joint decorrelation, a versatile tool for
multichannel data analysis. Neuroimage 98, 487–505.

169

C. Cai et al.

NeuroImage 188 (2019) 161–170
Sekihara, K., Nagarajan, S.S., 2017. Subspace-based interference removal methods for a
multichannel biomagnetic sensor array. J. Neural. Eng. 14 (5), 051001.
Sekihara, Kensuke, Kawabata, Yuya, Ushio, Shuta, Sumiya, Satoshi, Kawabata, Shigenori,
Adachi, Yoshiaki, Nagarajan, Srikantan S., 2016. Dual signal subspace projection(DSSP): a novel algorithm for removing large interference in biomagnetic
measurements. J. Neural. Eng. 13 (3), 036007.
Taulu, Samu, Hari, Riitta, 2009. Removal of magnetoencephalographic artifacts withtemporal signal-space separation: demonstration with single-trial auditory-evoked
responses. Hum. Brain Mapp. 30 (5), 1524–1534.
Taulu, Samu, Simola, Juha, 2006. Spatiotemporal signal space separation method forrejecting nearby interference in MEG measurements. Phys. Med. Biol. 51 (7), 1759.
Widrow, Bernard, Glover, John R., McCool, John M., Kaunitz, John, Williams, Charles S.,
Hearn, Robert H., Zeidler, James R., Dong, JR Eugene, Goodlin, Robert C., 1975.
Adaptive noise cancelling: principles and applications. Proc. IEEE 63 (12),1692–1716.
Wipf, David, Nagarajan, Srikantan, 2009. A uniﬁed bayesian framework for MEG/EEG
source imaging. Neuroimage 44 (3), 947–966.
Wipf, David P., Owen, Julia P., Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan
S., 2010. Robust bayesian estimation of the location, orientation, and time course ofmultiple correlated neural sources using MEG. Neuroimage 49 (1), 641–655.
Zumer, Johanna M., Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan S., 2008.
Probabilistic algorithms for MEG/EEG source reconstruction using temporal basis
functions learned from data. Neuroimage 41 (3), 924–940.

magnetoencephalography society clinical practice guideline 1: recording and analysisof spontaneous cerebral activity. J. Clin. Neurophysiol. 28 (4), 348–354.
Owen, Julia P., Wipf, David P., Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan
S., 2012. Performance evaluation of the champagne source reconstruction algorithm
on simulated and real M/EEG data. Neuroimage 60 (1), 305–323.
Nagarajan, Srikantan S., Attias, Hagai T., Hild II, Kenneth E., Sekihara, Kensuke, 2006.
A graphical model for estimating stimulus-evoked brain responses frommagnetoencephalography data with large background brain activity. Neuroimage 30(2), 400–416.
Nagarajan, Srikantan S., Attias, Hagai T., Sekihara, Kensuke, Hild, Kenneth E., 2006.
Partitioned factor analysis for interference suppression and source extraction. In:
International Conference on Independent Component Analysis and Signal Separation.
Springer, pp. 189–197.
Nagarajan, Srikantan S., Attias, Hagai T., Hild, Kenneth E., Sekihara, Kensuke, 2007.
A probabilistic algorithm for robust interference suppression in bioelectromagnetic
sensor data. Stat. Med. 26 (21), 3886–3910.
Nagarajan, S.S., Attias, H.T., Hild, K. Ee, Sekihara, K., 2007. A probabilistic algorithm forrobust interference suppression in bioelectromagnetic sensor data. Stat. Med. 26,
3886–3910.
Ossadtchi, A., Baillet, S., Mosher, J.C., Thyerlei, D., Sutherling, W., Leahy, R.M., 2004.
Automated interictal spike detection and source localization inmagnetoencephalography using independent components analysis and spatiotemporal clustering. Clin. Neurophysiol. 115 (3), 508–522.

170