Revised manuscript

Comparison of DSSP and tSSS algorithms for removing artifacts fromvagus nerve stimulators in magnetoencephalography data

Chang Cai1 , Huicong Kang1,2 , Heidi E. Kirsch3 , Danielle Mizuiri1 , Jessie Chen1 ,Abhishek Bhutada1 , Kensuke Sekihara4,5 , and Srikantan S. Nagarajan1

1

Biomagnetic Imaging Laboratory, University of California, San Francisco, 513 Parnassus Avenue,

S362, San Francisco, CA 94143, USA
2

Department of Neurology, Tongji Hospital affiliated to Tongji Medical College, Huazhong

University of Science and Technology, Wuhan, 430040, P.R.China
3

The UCSF Epilepsy Center and the Departments of Neurology and of Radiology and Biomedical

Imaging, University of California, San Francisco, CA 94143-0138.
4

Signal Analysis Inc., Hachioji, Tokyo, Japan

5

Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45

Yushima, Bunkyo-ku, Tokyo 113-8519, Japan

Corresponding author: Srikantan S. Nagarajan, E-mail: sri@ucsf.edu

Abstract
Objective: Large amplitude artifacts from vagus nerve stimulator (VNS) implants for refractory
epilepsy affect magnetoencephalography (MEG) recordings and are difficult to reject, resulting in
unusable data from this important population of patients who are frequently evaluated for surgicaltreatment of epilepsy. Here we compare the performance of two artifact removal algorithms for
MEG data: the dual signal subspace projection (DSSP) and temporally extended signal space
separation method (tSSS) algorithms. Approach: Each algorithm’s performance was first evaluated in simulations. We then tested the performance of each algorithm on resting-state MEG
data from patients with VNS implants. We also examined how each algorithm improved source
localization of somatosensory evoked fields (SEFs) in patients with VNS implants. Main results:
DSSP and tSSS algorithms have similar ability to reject interference both in simulated and inreal MEG data, if the origin location for tSSS is appropriately set. If the origin set for tSSS is
inappropriate, the signal after tSSS can be distorted, due to a mismatch between the internal
region and the actual source space. Both DSSP and tSSS are able to remove large amplitudeartifacts from outside the brain. DSSP might be a better choice than tSSS when choice of origin
location for tSSS is difficult. Significance: Both DSSP and tSSS algorithms can recover distorted
MEG recordings from people with intractable epilepsy and VNS implants, improving epileptic
spike identification and source localization of both functional activity and epileptiform activity.

1

Introduction

Magnetoencephalography (MEG) has made a tremendous impact on noninvasive brain mapping because it provides good spatial and excellent temporal resolution and can offer valuable information foruse in clinical neurology and in basic neuroscience. Unfortunately, the strength of brain MEG activity
is typically 50-1000 fT, and MEG recordings are prone to contamination arising from inevitable nonbiological sources like power lines and trains, and from biological sources outside of the brain like theheart. Though most of this interference is of similar magnitude to brain activity, some can be quite
large, arising from sources of magnetic field in the environment or even from the patient’s own bodyor head, including, for example, artifact from dental work. A prototypical example the interference
arising from vagal nerve stimulator (VNS) implants in patients with intractable epilepsy; VNS artifact
makes it very difficult for us to consistently and reliably detect and localize epileptiform activity inthese patients [1, 2, 3].
2

The preferred way to deal with interference from the environment is to prevent it. If this is
not possible, a variety of methods have been applied to minimize artifact in MEG recordings. One
commonly used method is to apply high-pass or band-pass filters. However, the brain signal of interest
may be in the same frequency range as the artifact; filtering may not be effective. One can also rejecttrials that contain large artifacts; however, this may cause unacceptable loss of signal from the brain[4]. Another commonly used method is averaging responses over trials: the fundamental idea forthis is that interference in different trials is statistically independent, whereas evoked signals are not.
However, this method requires a large number of trials, and evoked signals must be relatively similar
and robust [5].
Recently, data-driven approaches such as principal component analysis (PCA) and independent
component analysis (ICA) have been used for artifact removal. However, these methods require users
to make subjective choices during application (e.g. choice of threshold in PCA; selection of interferencecomponents for removal in ICA) [6, 7]. Joint decorrelation is another method commonly supposed to
be robust to many types of interference problems, but its use requires the design of different bias filters
for different interference types, and thus to some extent requires prior knowledge of the interference[8].
Ideally, any interference removal algorithm would be robust and broadly capable of rejection of as
many types of interference as possible. Currently, options are restricted to specific MEG hardware
platforms. For example, in the Elekta platform, the temporally extended signal space separation
method (tSSS) was developed based on spherical harmonic decomposition, and offers a potential solution [4]. However, this validity of this method has only been demonstrated for Elekta and has notbeen shown for other hardware platforms. Recently, we proposed an MEG hardware platform independent algorithm named Dual Signal Subspace Projection (DSSP) for removal of large interferencein biomagnetic measurements. In a retrospective clinical study we have shown the potential of DSSP
to deal with high amplitude, periodic interference. DSSP helps to recover distorted MEG recordings
from people with intractable epilepsy and VNS implants, making epileptic spike identification easierand spike mapping better [3, 9].
The performance of DSSP and tSSS have not been directly compared. In this study, we compare
the performance of these two algorithms using simulated data and real MEG recordings from patients
implanted with VNS. We test the algorithms’ ability to remove interference from resting-state MEG

3data in patients with VNS implants and also to improve source localization of somatosensory evoked
fields (SEF) using Champagne, an empirical Bayesian source reconstruction algorithm described previously [10]. Results show that DSSP and tSSS algorithms have similar performance for interference
rejection. However, for the tSSS algorithm to have comparable performance to DSSP, the origin location for the spherical harmonic decompositions must be set properly. Section 2 briefly describes
the DSSP and tSSS algorithms. Section 3 describes the computer simulation used for performancecomparison, followed by a description of their performance on real data analysis in Section 4. Section 5
provides a short discussion followed by a conclusion in Section 6.

2

Methods

2.1

Structure of DSSP and tSSS algorithms

This section briefly describes the DSSP and tSSS algorithms. Detailed explanations of these algorithms
have been published [3, 11, 12]. These algorithms share a common structure outlined in Fig. 1. As
shown in this figure, these algorithms consist of three steps: signal nulling, estimation of time-domaininterference subspace, and time-domain signal space projection (SSP) [11]. The signal nulling step
involves projecting out the signal components from the measured sensor data by applying the spatialdomain SSP algorithm [13, 14], i.e., by projecting the data matrix onto the subspace orthogonal to
the signal subspace, which is called the noise subspace. This procedure is referred to as the signal
nulling in this paper, and the operator used for it is called the signal nulling operator.
Theoretically, an ideal signal nulling operator is the noise subspace projector. However, since thesignal and noise subspace projectors are unknown, we must use something that can substitute forthese projectors. For this, the DSSP algorithm uses the pseudo-signal subspace projector P̆ S definedin Eq. (11) in Appendix. The tSSS algorithm uses the SSS internal subspace projector P int definedin Eq. (30), of which derivation is described in detail in Section B.
When using P int , the internal region needs to be matched with the source space1 , and this matching
between the internal region and the source space can be attained by setting the origin for the sphericalharmonics decomposition at a proper location. Note that the location of this origin is a parameter
1 The source space indicates the region in which sources can exist. The whole brain, in general, is considered as thesource space in MEG.

4

a user should specify when implementing the tSSS algorithm. This is because the tSSS algorithm
computes the SSS basis vectors using the spherical polar coordinate with this origin, and then derives
P int using the SSS basis vectors. Thus, this origin of the polar coordinate is called the SSS origin
in this paper. If the location of the SSS origin is not properly set, a mismatch between the source
space and the internal region arises that can cause a significant distortion in tSSS interference removal
results, as shown in our computer simulation in Section 3.
Note that the only practical difference between the DSSP and tSSS algorithms is in the mannerof approximating the signal subspace projector in the first step; the following steps are exactly thesame for these algorithms. The second step estimates the interference subspace by computing theintersection between the row spaces of the two modified data matrices obtained with and without
signal nulling. The third step implements the time-domain SSP by using the interference subspace
projector estimated in the second step. These two steps are further explained below.

2.2

Estimation of interference subspace and interference removal

Regarding interference removal, the model for the measured data is defined in a matrix form such that
B = BS + BI + Bε,(1)where B is the data matrix, B S the signal matrix, B I the interference matrix, and B ε a sensor noiseb S , which is equal to P̆ S in the DSSP
matrix. The approximate signal subspace projector is denoted P
b S satisfies the relationship: P
b S B S ≈ B S andalgorithm and P int in the tSSS algorithm. This P
b S . The both DSSP and tSSS algorithms modify thethe signal nulling operator is obtained as I − P
b S and I − P
b S to B, resulting inmeassured data matrix in two ways by applying P
b S B = BS + P
b S BI + P
b S Bε,B in = P(2)b S )B = (I − P
b S )B I + (I − P
b S )B ε ,B out = (I − P(3)where B in and B out are two kinds of modified data matrices obtained from B.
In Eqs. (2) and (3), we can observe that B I -related terms exist as common components in B inb S B I and (I − P
b S )B I are both equal to the row space of B I , which
and B out . The row spaces of P
is equal to the time-domain interference subspace KI . Also, the row space of B I can be extracted
as the intersection between the row spaces of B in and B out . Therefore, the time-domain interference5subspace is obtained as
KI = row space of B I ≈ row space of B in ∩ row space of B out .

(4)

The basis vectors of the intersection can be derived using the algorithm described in [11, 15].
Once the orthonormal basis vectors of the intersection, ψ 1 , . . . , ψ r , are obtained, we can compute theprojector onto the intersection ΠX such that
ΠX = [ψ 1 , . . . , ψ r ][ψ 1 , . . . , ψ r ]T .

(5)

Using this ΠX as the projector onto the interference subspace KI , the interference removal (andresultantly the estimation of the signal matrix) is achieved by applying the time-domain SSP:b S = B(I − ΠX ) = B(I − [ψ 1 , . . . , ψ r ][ψ 1 , . . . , ψ r ]T ),B(6)b S is the estimated signal matrix.
where B

3

Computer simulations

As mentioned in Section 2.1, the DSSP and tSSS algorithms differ in their signal-nulling operators.
The tSSS algorithm uses the projector onto the SSS internal subspace as the signal-nulling operator.
This SSS internal subspace projector, however, requires that the location of the SSS origin, which isa user-specified parameter, must be properly chosen. If it is not, a mismatch between the internal
region and the source space could arise and this mismatch may cause signal distortion in the tSSS
interference-removed results. Computer simulation in this section explores this point, and presents
two kinds of examples: one in which the location of the SSS origin is properly chosen and the otherin which it is not properly chosen. We then compare results of the DSSP algorithm with those of thetSSS algorithm obtained with the proper location of the SSS origin.
In our computer simulation, a sensor alignment of the 275-channel whole-head sensor array fromthe OmegaTM (CTF MEG Technology, Coquitlam, BC,Canada) neuromagnetometer was used. The
coordinate system and source-sensor configuration used in the computer simulation are depicted inFig. 2(a). We assumed a single source at a location shown by a square in this figure; its location waspresumably in the right parietal area near the primary somatosensory cortex, representing a source
located near the border of the source space.
6

We assumed two locations of the SSS origin. These locations are marked as ”A” and ”B” in
Fig. 2(a), and referred to as location A and location B in this section. The location A is 10 cm below
the vertex of the sensor array, and approximately 4 cm above a typical location of the origin of asingle sphere fitted to the inner border of a head MRI. The origin of this sphere is called the local
sphere origin, and abbreviated as LSO in this paper. The location B is 16 cm below the vertex of thesensor array, and 2 cm below LSO. These two locations are shown by red filled circles in Fig. 2(a).
In order to see if these locations of the SSS origin are proper choices or not, we need to visualize
where the signal nulling region exists and how widely it is extended, depending on the location of the
SSS origin. For this visualization, the signal gain by the signal nulling operator was plotted. The
signal gain GS (r) is defined such that
GS (r) =b S )y S (r)kk(I − P,ky S (r)k(7)

where y S (r) is the signal vector due to a single source at r. The distribution of GS (r) was computedby scanning a single source over the whole source space. The plots of GS (r) are shown in Fig. 2(b).
In this figure, the left-hand, middle, and right-hand columns, respectively, show maximum-intensity
projections of the three-dimensional distribution of GS (r) onto the y − z, x − z, and x − y planes.
The colors of these figures are in accordance with the values of the signal gain GS (r), and the colorof the inner most region (dark blue) indicates the value of GS (r) less than 1/100.
In Fig. 2(b) top row in this figure shows the the signal gain mapping obtained using the pseudo
signal subspace projector P̆ S (Eq. (11)) used in the DSSP algorithms. The middle and bottom rows
show the signal gain mappings from the SSS internal subspace projector P int . The middle row shows
the case in which the SSS origin was set at location A, and the bottom shows the case in which the
SSS origin was set at location B. These plots indicate that, for the tSSS method, the signal nulling
region depends directly on the SSS origin location. When the origin was set at location A, the signal
nulling region extends to cover the location of the source near the border of the source space, whichsuggests that location A can be an appropriate choice of the SSS origin location. However, when theorigin was set at location B, the source was located slightly outside the signal nulling region, whichsuggests that location B may be an inappropriate choice that causes a mismatch between the internal
region and the source space.
We then computed the internal and the external components of the signal magnetic field generated
from that source located near the source-space border (marked by a square in Fig. 2(a)). To do so,7

we assign the time course expressed by an exponentially-dumped sinusoid to this source. The signalsource time course was then projected to the sensor time courses through the sensor lead field, which
is obtained using the homogeneous spherical head model [16]. The generated signal data B S with
2400 time points was shown in Fig. 3(a).
b S B S and B out = (I −The top and the bottom panels in Fig. 3(b), respectively, show B in = P
b S )B S for the DSSP case where P̆ S is used as P
b S . Here, the bottom panel show the results ofP
the signal nulling, and these results show that signal nulling was almost attained, although a small
amount of signal remained.
Figure 3(c) shows the results of B in (top panel) and B out (bottom panel) for the tSSS cases inwhich the operator P int was computed with the SSS origin set at location A. We can see that signal
nulling was nearly effectively attained in this case. Figure 3(d) shows again the results of B in (toppanel) and B out (bottom panel) for the tSSS cases in which the operator P int was computed with the
SSS origin set at location B. In this case, a considerable amount of the signal remained in the bottom
panel. These results can be predicted from the signal-gain plots in the bottom panel of Fig. 2(b)where the source was outside the signal nulling region.
Next, we performed interference removal experiments. We put a single interference source justoutside the sensor array (50 cm distant from the center of the array). This interference simulated
some types of brain stimulator artifacts. The interference source was assumed to have a low-passfiltered random time course. The measured data was generated by adding the interference data to thesignal data (shown in Fig. 3(a)) with the interference to signal ratio (ISR) equal to 10. The generated
interference-overlapped data are shown in Fig. 4(a). Here, since the interference magnetic field wasten times stronger than the signal magnetic field, the data were dominated by the interference.
The DSSP algorithm was applied to these generated data, and the results are shown in the top
panel in Fig. 4(b). The field contour map at 1200 ms of these results is shown in the top-right panelof Fig. 4(c). Note that the ground truth of the sensor time courses is shown in Fig. 3(a). The groundtruth of the field contour map (obtained from the signal magnetic field) is shown in the top-left panelof Fig. 4(c). Comparison between the DSSP results and these ground truth results indicates that the
DSSP algorithm almost effectively removed the interference.
The second panel from the top in Fig. 4(b) shows the results of the tSSS algorithm with the SSS
origin set at location A. The field contour map of these tSSS results is shown in the bottom-left panel

8

of Fig. 4(c). These tSSS results are almost the same as the ground truth results, indicating that tSSS
algorithm effectively removed the interference.
The bottom panel in Fig. 4(b) shows the results of the tSSS algorithm with the origin set at
location B. The field map of these results is shown in the bottom-right panel of Fig. 4(c). Althoughthese time courses suggest that the interference was removed pretty well, the field map indicates thatsome amount of distortion was caused in the interference-removed results. These results in Figs. 3
and 4 confirm that location A represents a proper location of the SSS origin, and location B is not a
proper one.
The cause of the signal distortion can be explained in the following manner. When the signal
magnetic field has components belonging to the external subspace, these external components are
considered as a part of the interference and they are removed through the tSSS interference removal
process, resulting in the distortion of the signal magnetic field.
In summary, from the results of our computer simulation, we can draw conclusions that DSSP
and tSSS algorithms have almost the same performance for interference removal. However, in orderfor the tSSS algorithm to have a performance comparable to that of the DSSP algorithm, the SSS
origin should be properly set so that the internal region matches the source space. If the origin is notproperly chosen, a mismatch between the internal region and the source space can distort the signal
magnetic field.

4

Experiments using real MEG data from patients with a implanted VNS device

4.1

Collecting MEG data and general aspects of analysis

We chose 10 epilepsy patients with VNS implants who underwent a clinical MEG study as partof epilepsy surgery evaluation at the University of California, San Francisco (UCSF) Biomagnetic
Imaging Laboratory (BIL) between November 24th, 2004 and May 6th, 2016. Prior to the MEG scan,
all patients had high-resolution epilepsy protocol 3T T1-MRI scans for coregistration of localized
spikes. MEG recordings were performed inside a magnetically shielded room with a CTF 275 channel
whole-head axial gradiometer system.

9

We collected two kinds of MEG data from each patient: somatosensory evoked MEG and restingstate spontaneous MEG. Somatosensory evoked data sets (sensory evoked field paradigm) were collected with the sample rate set at 1.2 kHz. The patients’ right index finger was tapped for a total of
256 trials. A peak should typically exist around 50 ms after stimulation in the contralateral (in thiscase, the left) somatosensory cortical area for the hand, i.e., dorsal region of the postcentral gyrus.
The leadfield for each subject was calculated in NUTMEG 4.1 [17] with an 8 mm voxel grid
resulting in approximately 5300 voxels generated. We chose a pre-stimulus time window from -100 ms
to 0 ms (120 sample points) as a baseline period and the post-stimulus time window was set from 0
ms to 200 ms (240 sample points) where the stimulus onset time was at 0 ms time point.
Source localization of the somatosensory data was conducted using a united Bayesian frameworkfor MEG/EEG source imaging that includes Variational Bayes Factor Analysis (VBFA) for noise
learning [18] and a sparse Bayesian algorithm (Champagne) for source localization [10, 19, 20]. The
localization results were overlaid on the standardized brain atlas provided by the NUTMEG-SPM8
pipeline.
The resting-state MEG data were collected at 1.2 kHz then downsampled to 600 Hz in a passbandof 0-75 Hz. Thirty to forty minutes of spontaneous data were obtained in intervals of 10-15 min withthe patient asleep and awake. The position of the patient’s head in the dewar relative to the MEG
sensors was determined using indicator coils before and after each recording interval to verify adequate
sampling of the entire field. The data were then bandpass filtered offline, at 1-70 Hz. More details of
the recording methods have been previously described [5, 9].
These resting-state data were then used for the epilepsy spike analysis. Spikes were visually
identified by an expert clinical neurophysiologist (HK) and results were confirmed by a board-certified
clinical neurophysiologist and epileptologist (HEK) using the same criteria as used in our prior work
[9].
For SEF data, we ran DSSP and tSSS on single trials and then averaged the results. For each trial,
the time-window for analysis was 0.3 second obtained at a sampling rate of 1.2 KHz. For epileptic
spike data, we used a 10 second time-window for DSSP and tSSS analysis at sampling rate for ourdata is 600 Hz.

10

4.2

VNS artifact removal in somatosensory evoked MEG

Somatosensory MEG recordings from selected channels are shown for four representative patients in
Fig. 5. The original MEG recordings (before artifact removal) are shown in Fig. 5 (a). These originalrecordings feature time courses with partial periodicity, low frequency and high amplitude, typical of
VNS artifact. The results of DSSP artifact removal are shown in Fig. 5 (b). In these results, theperiodic feature is greatly diminished and a peak around the latency of 50 ms is clearly observed.
To use the tSSS application, we need to provide the location of the SSS origin. In order to evaluate
the impact of origin selection, we applied the tSSS algorithm while setting the SSS origin to each offour different locations: 4 cm below LSO, LSO, 4 cm above LSO, and 8 cm above LSO2 . The resultsof tSSS artifact removal using each of these four SSS origin locations are shown in Fig. 5(c)–(f). It
can be seen that tSSS and DSSP give similar artifact-removal results when the SSS origin was set at
4 cm or 8 cm above LSO. For the other two origin settings, however, tSSS apparently failed to remove
most of the overlapped artifacts, and some artifacts remained in the processed results.
The source localization was conducted to see whether the localization improved upon the incorporation of the DSSP and tSSS algorithms and to compare their performance. The results for the same
four subject are shown in Fig. 6. Results from the original data are shown in Fig. 6 (a), in which thesource localization obviously failed for all four patient cases. Results from the DSSP processed data
are shown in Fig. 6 (b). We can see that the sources were consistently localized near the primarysomatosensory area in these results. Results from tSSS processed data are shown from Fig. 6 (c) to(f). The results similar to those from DSSP application was obtained using the tSSS algorithm with
the SSS origin set at 4 cm above LSO (Fig. 6(e)). The tSSS application with other settings of the
SSS origin, however, gave no clear results where source locations were inconsistent across these four
patients.
The results of source localization experiments for ten patients are summarized in Fig. 7 where thebar graphs indicate the number of successful localizations. The success or failure of source localization
results was judged by the locations of local peaks in the source reconstruction results. if a peakmaximum was within a region whose label was ”Postcentral gyrus” on the standardized brain atlas,
we considered the source localization successful.
2 As mentioned in Section 3, LSO is the location of the origin of a single sphere fitted to the inner border of thesubject’s head MRI.

11

We can see that using the DSSP-processed data, the source localization was successful for all ten
subjects. In the tSSS results, eight of ten attained the successful localization when the location of the
SSS origin was set at 4 cm above LSO. However, for other settings of the SSS origin, the number ofsuccessful localizations was no more than eight. For only one subject were results localized using theoriginal VNS-overlapped MEG sensor data.
In summary, the results from these somatosensory experiments indicate that DSSP and tSSS algorithms have nearly the same performances for VNS artifact removal. However, when implementing thetSSS algorithm, the location of the SSS origin should be properly set to yield performance comparable
to that of the DSSP algorithm. If the SSS origin location is not properly set, the performance of thetSSS algorithm can be significantly worse than that of the DSSP algorithm. Also, the results in Fig. 7
clearly indicate that the proper choice of the SSS origin location is in the vicinity of 4 cm above LSO.

4.3

Artifact removal in resting-state MEG

We next performed VNS artifact removal experiments using resting-state MEG. In Fig. 8(a), original
sensor data for a representative patient are shown. A large VNS artifact with a partial periodicity
can be seen. The DSSP processed sensor data are shown in Fig. 8(b) in which the periodic featureis greatly diminished and the background looks similar to that of most patients with epilepsy andwithout a VNS device.
According to the results of somatosensory experiments, the proper location of the SSS origin isaround 4 cm above LSO. In order to confirm if this is also true for resting state MEG in epilepsy, weagain tested the same four locations of the SSS origin when implementing the tSSS algorithm. The
sensor data after tSSS artifact removal with the SSS origin set at these four locations are shown in
Fig 8(c)–(f).
As can be seen in Fig 8(c) and (d), when the SSS origin was set at 4 cm below LSO or at LSO, thetSSS algorithm failed to remove all components of the artifact. When the SSS origin was set at 4 cm
above LSO or at 8 cm above LSO, results similar to the DSSP results were obtained. The third rowof Fig. 8 shows the sensor topographical maps. The topographical maps also confirm that tSSS with
the SSS origin set at 4 cm above LSO or at 8 cm above LSO can give results very similar to thosefrom the DSSP algorithm.
We performed visual identification of epileptiform activity (spikes) in these data sets. The vertical12

red and green lines indicate spikes that were identified in the DSSP-processed data but not identified
in the original sensor data. These spikes were also identified in the results from tSSS with the SSS
origin set at 4 cm above LSO or at 8 cm above LSO.
To see how close the DSSP and tSSS results are, we computed the correlation between DSSP- andtSSS-processed resting-state sensor time courses. In these resting-state measurements, each subject
data contained approximately 100 trials with 275 sensor channels. Correlation coefficients were averaged across trials and sensor channels, and then averaged across ten subjects. Averaged correlation
coefficients are plotted in Fig. 9 with respect to the SSS origin location. As is shown here, the correlation coefficient attains its maximum value when the SSS origin location is set at 4 cm above theLSO, and the maximum correlation reaches 0.85.
From experiments using resting-state MEG data sets, we can draw the conclusion that the DSSP
and tSSS algorithms yield most similar results when the SSS origin is set near 4 cm above LSO.
However, for the other four locations of the SSS origin, the performance of the tSSS algorithm is not
as good as the performance of the DSSP algorithm.

5

Discussion

In this paper we compare the DSSP and tSSS algorithms for MEG artifact removal. We conductedcomputer simulation and experiments using real MEG data, which were collected from patients with
an implanted VNS device. The two algorithms have similar performances for VNS artifact removal
both in simulations and in real MEG data.
However, unlike the DSSP algorithm, the tSSS algorithm has a user-defined parameter, the SSS
origin, which is the origin of the spherical polar coordinate for computing the SSS basis vectors.
We found that, in order for the tSSS algorithm to have performance comparable to that of the DSSP
algorithm, the location of the SSS origin must be properly chosen. If not, there is a mismatch between
the internal region and the source space, which may cause a distortion of the signal magnetic field.
According to our experiments using MEG data from ten patients, the proper location of the SSS
origin is in the vicinity of 4 cm above the origin of a sphere fitted to the inner border of subject’s
head MRI. Naturally, the question may arise if this proper location depends on the shape of thesensor helmet and sensor arrangement, implying that it might change across different MEG hardware

13systems. Answering this question was not easy for us in isolation, because our access to MEG hardware
is currently limited to the CTF MEG system installed in our laboratory. Nonetheless, we intend in
future to explore whether our results apply to different MEG hardware, by collaborating with users
of different systems.
Recent algorithms based on learning statistical properties of the interference have been hailed
as reliable and robust. One such example is the partitioned factor analysis (PFA) algorithm [1,21, 22]. The PFA involves learning of a probabilistic model of the interference from the data in a
separate interval (where only the interference exists) and applying this to an active period (when bothinterference and true signals of interest exist). This method handles most types of interference well,
but since it relies on the availability of separate measurements that capture the statistical propertiesof the interference, its use is limited to situations where such separate measurements are appropriate[5]. However, these methods would not be effective for interference that is constant in the data, as in
the case of VNS implant artifacts, or if the interference is present only in the active dataset of interest,
as in artifacts that occur in active period but not during a baseline measurement period. Also, thesealgorithms may not be effective for interference of extremely large magnitude relative to the signals
being estimated, as is often seen in MEG data from patients with VNS implants. Nevertheless, future
performance comparison with these unsupervised machine learning benchmark algorithms for artifact
removal are warranted.
In this paper, we have tested the performances of the DSSP and tSSS algorithms for VNS artifact
removal. It may be interesting to explore the possibilities of applying these algorithms to removalof some types of physiological artifacts common in MEG. Indeed, DSSP and tSSS may be used forremoval of cardiac, ocular or dental artifacts. In our research group, we have successfully deployed
DSSP and tSSS for all of these sources of artifacts in MEG data. Algorithm performance of both
DSSP and tSSS requires some degree of spatial separation of the sources of artifacts from the sourcesof interest [3, 11]. Therefore, performance is expected to be better for artifactual sources that may
be further from the brain. Consistent with this idea, in our experience, we get better performance forremoval of cardiac and dental sources in comparison to ocular artifacts (data not shown). Some typesof independent component analysis (ICA) algorithms are known to be effective for the revoal of such
artifacts. Nevertheless, our future plan includes the performance comparison of these two algorithms
with the ICA algorithms, because the ICA algorithms require subjective selection of interference

14components but DSSP/tSSS algorithms do not require such procedure.

6

Conclusion

Both DSSP and tSSS are able to remove large amplitude artifacts from outside the brain. DSSP can bea better choice than tSSS when the origin location for tSSS is not easy to choose. Both algorithms canrecover distorted MEG recordings from people with intractable epilepsy and VNS implants, improving
epileptic spike identification and source localization of functional activity and epileptiform activity.

Acknowledgment
The authors would like to thank Susanne Honma and Anne Findlay for collecting much of the MEG
data, and all members and collaborators of the Biomagnetic Imaging Laboratory for their support.
This work was supported by NIH grants R01EB022717, R01DC013979, R01NS100440, R01DC017696,
and UCOP-MRP-17-454755.

A

Pseudo-signal subspace projector

The DSSP algorithm uses the pseudo-signal subspace projector to approximate the signal subspace
projector. This section describes how to derive this projector. We define voxels over the source space,in which the voxel locations are denoted r 1 , . . . , r N . The augmented leadfield matrix over these voxel
locations is defined as
F = [L(r 1 ), . . . , L(r N )],(8)and the pseudo-signal subspace ĔS is defined such thatĔS = column space of F = span{L(r 1 ), . . . , L(r N )}.

(9)

The true signal subspace is defined as
ES = span{L(r q1 ), . . . , L(r qq )}where r q1 , . . . , r qq are the source locations assuming that there are q sources in the source space. If thevoxel interval is sufficiently small and voxel discretization errors are negligible, each of source locations15

r q1 , . . . , r qq is equal to one of voxel locations, and the space spanned by {{L(r q1 ), . . . , L(r qq )} is a subsetof the space spanned by {L(r 1 ), . . . , L(r N )}, resulting in the relationship ĔS ⊃ ES .
Let us derive the orthonormal basis vectors of the pseudo-signal subspace. To do so, we computethe singular value decomposition of F :
F =

M
X

λj ej f Tj .

(10)

j=1

If the singular values λ1 , . . . , λτ are distinctively large and other singular values λτ +1 , . . . , λM arenearly equal to zero, the leading τ singular vectors e1 , . . . , eτ form orthonormal basis vectors of thepseudo-signal subspace ĔS [23]. Thus, the projector onto ĔS is obtained using
P̆ S = [e1 , . . . , eτ ][e1 , . . . , eτ ]T .

(11)

Note that, since ĔS ⊃ ES holds, the orthogonal projector (I − P̆ S ) removes the signal components,
resulting in (I − P̆ S )B S = 0.

B

SSS internal subspace projector

B.1

Spherical harmonics expansion of the magnetic field

This section describes the derivation of the SSS internal subspace projector used for approximating
the signal subspace projector in the tSSS algorithm. One fundamental assumption of the SSS method
is that the sensors are installed in a source-free region, which is referred to as the sensor region. Then,the magnetic field at r, B(r), is expressed using the spherical polar coordinate r = (r, θ, φ) by
B(r) = −µ0∞ X`
X`=1 m=−`

∞

α`,m

`

X X
ν `,m (θ, φ)− µ0β`,m r`−1 ω `,m (θ, φ),
`+2
r

(12)

`=1 m=−`

where µ0 indicates the magnetic permeability of free space. In Eq. (12), ν `,m (θ, φ) and ω `,m (θ, φ) are
the modified vector spherical harmonics [24, 25].
In the right-hand side of Eq. (12), the first term represents the magnetic field generated fromsources located closer to the origin than the sensors are. The second term represents the magnetic
field from sources located farther from the origin than the sensors are. The region closer to the origin
than the sensors is referred to as the internal region, and the region farther from the origin than thesensors is referred to as the external region. Note here that the origin location of the spherical polar
16

coordinate r = (r, θ, φ) is a user-defined parameter, and a user should carefully choose the origin.
This is because a proper choice of the origin can result in the signal sources to be located within theinternal region but an improper choice may cause a situation where some signal sources are locatedoutside the internal region, as shown in our computer simulation and experiments using the real MEG
data.

B.2

Derivation of SSS basis vectors

Let us derive the SSS basis vectors. To do so, the magnetic signal detected by the j-th sensor isdenoted by bj , and the location and the normal vector of the j-th sensor by r j and ζ j . Then, we havebj = B(r j ) · ζ j = bjint + bjext ,(13)

where the notation “·” indicates taking the inner product between two vectors. Here, bjint and bjext
respectively represent magnetic components originating from the internal and external regions. These
components are expressed as
bjint = −

∞ X
`
X

α`,m



ν `,m (θj , φj ) · ζ jrj`+2`=1 m=−`

bjext = −

∞ X
`
X

,

(14)



β`,m rj`−1 ω `,m (θj , φj ) · ζ j ,(15)`=1 m=−`

where we set µ0 = 1 for simplicity. Let us define the internal and external components of the vector
M TT1b as bint = [b1int , . . . , bM
int ] and bext = [bext , . . . , bext ] , which are expressed such thatbint =∞ X`
X

α`,m c`,m ,(16)

β`,m d`,m ,(17)`=1 m=−`

bext =∞ X`
X`=1 m=−`

where column vectors c`,m and d`,m are given by
1c`,m = d`,m = 



r1`+2

[ν `,m (θ1 , φ1 ) · ζ 1 ]..
.

1`+2 [ν `,m (θM , φM ) · ζ M ]rM

r1`−1 [ω `,m (θ1 , φ1 ) · ζ 1 ]
..
.
`−1
rM
[ω `,m (θM , φM ) · ζ M ]

17




,





.



(18)

Truncating the summation with respect to the multiple order ` to LC for bint and LD for bext , wefinally obtainb = bint + bext =

LC X
`
X

α`,m c`,m +`=1 m=−`LD X`
X

β`,m d`,m .

(19)`=1 m=−`

Thus, defining
C = [c1,−1 , c1,0 , c1,1 , . . . , cLC ,LC ] ,(20)

D = [d1,−1 , d1,0 , d1,1 , . . . , dLD ,LD ] ,

(21)

T

(22)

T

(23)

α = [α1,−1 , α1,0 , α1,1 , . . . , αLC ,LC ] ,β = [β1,−1 , β1,0 , β1,1 , . . . , βLD ,LD ] ,we obtain





 α 
b = bint + bext = Cα + Dβ = [C, D] 
 = Sx,β(24)

where S = [C, D] and x = [αT , β T ]T .

B.3

Derivation of internal and external subspace projector

Equation (24) is the basis for estimating the internal and external components bint and bext fromb T ]T is obtained as:
b = [αbT , βgiven magnetic signal data b. That is, the least squares estimate xb = (S T S)−1 S T b.
x

(25)

b
b
bint = C α,

(26)

b
b
bext = D β.

(27)Then, bint and bext are estimated as

After some derivations, we can finally obtainb ≈ C T (CC T + DD T )−1 b.
α(28)

The derivation of Eq. (28) is presented in [11][26].
Using Eqs. (26) and (28), we obtainbbint ≈ CC T (CC T + DD T )−1 b,
18

(29)

and thus, the internal component bint can be extracted by multiplying
P int = CC T (CC T + DD T )−1(30)

with the data vector b. That is, the matrix P int acts as a projector that passes the internal components
and blocks the external ones3 . Therefore, we call P int the SSS internal subspace projector in thispaper. In exactly the same manner, the SSS external subspace projector P ext is derived as
P ext = DD T (CC T + DD T )−1 .

(31)This P ext passes the external components but blocks the internal ones. It can be seen from Eqs. (30)and (30) that the relationship I − P int = P ext holds.
This operator P int projects a vector into the SSS internal subspace spanned by the columns of thematrix C, and P int does not change the vector if it is already in the internal subspace. Therefore, theapplication of I−P int eliminates the internal components, and we have the relationship (I−P int )B S =0 if the signal B S is generated from the internal region. However, if there is a mismatch between theinternal region and the source space and some signal sources are located outside the internal region,
the signal nulling operator I − P int only removes a part of B S that belongs to the column space ofC.

References[1] S. S. Nagarajan, H. T. Attias, K. E. Hild II, and K. Sekihara, “A graphical model for estimating
stimulus-evoked brain responses from magnetoencephalography data with large background brain
activity,” NeuroImage, vol. 30, no. 2, pp. 400–416, 2006.
[2] S. Taulu and J. Simola, “Spatiotemporal signal space separation method for rejecting nearby
interference in MEG measurements,” Physics in Medicine & Biology, vol. 51, no. 7, p. 1759,
2006.
[3] K. Sekihara, Y. Kawabata, S. Ushio, S. Sumiya, S. Kawabata, Y. Adachi, and S. S. Nagarajan,“Dual signal subspace projection (DSSP): a novel algorithm for removing large interference inbiomagnetic measurements,” Journal of Neural Engineering, vol. 13, no. 3, p. 036007, 2016.
3 Note that since (P

int )

2 = P

T
int and (P int )= P int do not hold, P int is not exactly a projector in strict linear

algebra terminology.

19[4] S. Taulu and R. Hari, “Removal of magnetoencephalographic artifacts with temporal signal-space
separation: demonstration with single-trial auditory-evoked responses,” Human brain mapping,vol. 30, no. 5, pp. 1524–1534, 2009.
[5] S. S. Nagarajan, H. T. Attias, K. E. Hild, and K. Sekihara, “A probabilistic algorithm for robustinterference suppression in bioelectromagnetic sensor data,” Statistics in medicine, vol. 26, no. 21,
pp. 3886–3910, 2007.
[6] A. Ossadtchi, S. Baillet, J. Mosher, D. Thyerlei, W. Sutherling, and R. Leahy, “Automatedinterictal spike detection and source localization in magnetoencephalography using independent
components analysis and spatio-temporal clustering,” Clinical Neurophysiology, vol. 115, no. 3,
pp. 508–522, 2004.
[7] S. S. Nagarajan, H. T. Attias, K. Sekihara, and K. E. Hild, “Partitioned factor analysis for interference suppression and source extraction,” in International Conference on Independent Component Analysis and Signal Separation, pp. 189–197, Springer, 2006.
[8] A. de Cheveigné and L. C. Parra, “Joint decorrelation, a versatile tool for multichannel data
analysis,” Neuroimage, vol. 98, pp. 487–505, 2014.
[9] C. Cai, J. Xu, J. Velmurugan, R. Knowlton, K. Sekihara, S. S. Nagarajan, and H. Kirsch, “Evaluation of a dual signal subspace projection algorithm in magnetoencephalographic recordings frompatients with intractable epilepsy and vagus nerve stimulators,” NeuroImage, vol. 188, pp. 161–
170, 2019.
[10] J. P. Owen, D. P. Wipf, H. T. Attias, K. Sekihara, and S. S. Nagarajan, “Performance evaluation of the champagne source reconstruction algorithm on simulated and real M/EEG data,”
Neuroimage, vol. 60, no. 1, pp. 305–323, 2012.
[11] K. Sekihara and S. S. Nagarajan, “Subspace-based interference removal methods for a multichannel biomagnetic sensor array,” Journal of neural engineering, vol. 14, no. 5, p. 051001, 2017.
[12] S. Taulu and J. Simola, “Spatiotemporal signal space separation method for rejecting nearby
interference in MEG measurements,” Phys. Med. Biol., vol. 51, pp. 1759–1768, 2006.

20

[13] M. Uusitalo and R. Ilmoniemi, “Signal-space projection method for separating MEG or EEG into
components,” Medical and Biological Engineering and Computing, vol. 35, no. 2, pp. 135–140,
1997.
[14] G. Nolte and G. Curio, “The effect of signal-space projection on dipole localization accuracy,”IEEE Trans. Biomed. Eng., vol. 46, pp. 400–408, 1999.
[15] G. H. Golub and C. F. Van Loan, Matrix computations, vol. 3. The Johns Hopkins University
Press, 2012.
[16] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[17] S. S. Dalal, J. M. Zumer, V. Agrawal, K. E. Hild, K. Sekihara, and S. S. Nagarajan, “NUTMEG:
A neuromagnetic source reconstruction toolbox,” Neurol Clin Neurophysiol., vol. 30, p. 52, 2004.
[18] K. Sekihara and S. S. Nagarajan, Electromagnetic brain imaging: A Bayesian perspective. Berlin,
Heidelber: Springer-Verlag, 2015.
[19] D. P. Wipf, J. P. Owen, H. T. Attias, K. Sekihara, and S. S. Nagarajan, “Robust Bayesianestimation of the location, orientation, and time course of multiple correlated neural sources
using MEG,” NeuroImage, vol. 49, pp. 641–655, 2010.
[20] C. Cai, K. Sekihara, and S. S. Nagarajan, “Hierarchical multiscale bayesian algorithm for robust
meg/eeg source reconstruction,” NeuroImage, vol. 183, pp. 698–715, 2018.
[21] S. S. Nagarajan, H. T. Attias, K. E. Hild, and K. Sekihara, “A probabilistic algorithm forrobust interference suppression in bioelectromagnetic sensor data,” Statistics in Medicine, vol. 26,
pp. 3886–910, 2007.
[22] J. M. Zumer, H. T. Attias, K. Sekihara, and S. S. Nagarajan, “Probabilistic algorithms for
MEG/EEG source reconstruction using temporal basis functions learned from data,” NeuroImage,
vol. 41, no. 3, pp. 924–940, 2008.
[23] I. C. Ipsen, Numerical Matrix Analysis: Linear Systems and Least Squares. SIAM, 2009.
[24] S. Taulu and M. Kajola, “Presentation of electromagnetic multichannel data: the signal space
separation method,” Journal of Applied Physics, vol. 97, no. 12, p. 124905, 2005.
21

[25] E. Hill, “The theory of vector spherical harmonics,” American Journal of Physics, vol. 22, no. 4,
pp. 211–214, 1954.
[26] K. Sekihara, “Signal space separation method for a biomagnetic sensor array arranged on aflat plane for magnetocardiographic applications: A computer simulation study,” Journal ofhealthcare engineering, vol. 2018, 2018.

22

DSSP

Interference overlapped data

Spatial subspace projector(from Lead field matrix)
tSSS
Spatial subspace projector(from spherical harmonic
decomposition)

Step#1Signal nulling by spatial-domain SSP
Step#2Estimation of time-domain interferencesubspace
Step#3
Interference removal by time-domain SSP

Interference removed data
Figure 1: Conceptual sketch of the DSSP and tSSS algorithms.

23

B

B
z (cm)Az (cm)Ax (cm)y (cm)y (cm)(a)x (cm)(b)
Figure 2: Simulation configuration and spatial distribution of the signal gain. a) Source-sensor configuration for the computer simulation with the two locations of the SSS origin (red circles) and thesource location (orange square). (b) Plots of the spatial distribution of the signal gain. Top row:
results from the projector P̆ S used in the DSSP algorithm (Eq. (11)). Middle row: results from the
projector P int used in the tSSS algorithm (Eq. 30) with the SSS origin set at location A. Bottom
row: results from the projector P int with the SSS origin set at location B. Left-hand, middle, andright-hand columns, respectively, show maximum-intensity projections of GS (r) onto the y − z, x − z,
and x − y planes.
24normalized

1

-1

(a)

time (ms)normalized

1-1

normalized

1

-1

(b)

time (ms)normalized

1-1

normalized

1

-1

(c)

time (ms)normalized

1-1

normalized

1

-1

(d)

time (ms)Figure 3: The internal and the external components of the signal magnetic field generated by DSSP
and tSSS with different configurations. (a) Generated signal data B S . (b) B in (top panel) and B out(bottom panel) from the DSSP algorithm. (c) B in (top panel) and B out (bottom panel) from thetSSS algorithm with the SSS origin set at location A. (d) B in (top panel) and B out (bottom panel)for the tSSS algorithm with the SSS origin set at location B. The ordinate of the plot in (a) indicatesnormalized values of the field intensity. The ordinates of the plots in (b)-(d) indicate normalized
values in which the normalization is performed by the maximum value in the top panel for (b) and(c), and that in the bottom panel for (d).

25

normalized

1

-1

(a)
normalized

0.1

normalized

-0.1
0.1

-0.1

normalized

0.1

-0.1

(b)

(c)
Figure 4: The results of the interference removal experiments. (a) Sensor time courses of generated
interference-overlapped data. The ordinate indicates normalized values of the field intensity. (b)Sensor time courses of the interference removed results. Top panel: results from the DSSP algorithm.
Middle panel: results from the tSSS algorithm with the SSS origin set at location A. Bottom panel:
results from the tSSS algorithm with the origin set at location B. The ordinates in these plots indicate
the normalized values in which the normalization is performed by the maximum value of the field
intensity for (a). (c) The maps of the magnetic field at t=1200. Top left: map of the original signal
magnetic field (ground truth). Top right: DSSP results. Bottom left: tSSS results with the SSS origin
set at location A. Bottom right: tSSS results with the SSS origin set at location B.
26

amplitude (fT)time (ms)

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Application of DSSP and tSSS on SEF recordings from patients with a VNS implant. MEG
recordings from selected channels are shown for four representative patients. Each row correspond tosensor data from each of four patients. Columns indicate (a) Original sensor data, (b) DSSP-processed
sensor data, (c) tSSS-processed sensor data with the SSS origin set at 4 cm below LSO, (d) tSSSprocessed sensor data with the SSS origin set at LSO, (e) tSSS-processed sensor data with the SSS
origin set at 4 cm above LSO, and (f) tSSS-processed sensor data with the SSS origin set at 8 cm
above LSO.

27

Max

Min

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Source localization results for the same four representative patients as in Fig. 5. Eachrow corresponds to results from each of four patients. Columns indicate: results from (a) original
sensor data, (b) DSSP-processed sensor data, (c) tSSS-processed sensor data with the SSS origin set
at 4 cm below LSO, (d) tSSS-processed sensor data with the SSS origin set at LSO, (e) tSSS-processed
sensor data with the SSS origin set at 4 cm above LSO, and (f) tSSS-processed sensor data with the
SSS origin set at 8 cm above LSO. For the source localization, we applied a pipeline using a united
Bayesian framework for MEG/EEG source imaging that includes Variational Bayes Factor Analysis
and a sparse Bayesian algorithm (Champagne). The color indicates the strength of the localized brain
activity.

28

Number of successful reconstruction

10
9
8
7
6
5
4
3
2
1
0

RAW

DSSP

LSO-4cmLSO

LSO+4cm

LSO+8cm

Figure 7: Number of successful localizations of SEF sources across ten patients. The source localization
is judged successful if a local peak in the reconstruction results is within a region whose label is”Postcentral gyrus” in the standardized brain atlas.

29

L

R

(a)

pT
18
15
12
9
6
3
0
-3
-6
-9
-12
-15
-18

pT
1.8
1.5
1.2
0.9
0.6
0.3
0
-0.3
-0.6
-0.9
-1.2
-1.5
-1.8

(b)

(c)

pT
3.6
3
2.4
1.8
1.2
0.6
0
-0.6
-1.2
-1.8
-2.4
-3
-3.6

L

R

(d)

pT
3.6
3
2.4
1.8
1.2
0.6
0
-0.6
-1.2
-1.8
-2.4
-3
-3.6

pT
1.8
1.5
1.2
0.9
0.6
0.3
0
-0.3
-0.6
-0.9
-1.2
-1.5
-1.8

(e)

(f)

pT
1.8
1.5
1.2
0.9
0.6
0.3
0
-0.3
-0.6
-0.9
-1.2
-1.5
-1.8

Figure 8: Sensor time courses of the resting-state MEG for a representative patient with a implanted
VNS device. (a) Original sensor data before artifact removal. (b) DSSP-processed results. (c) tSSS
results obtained with the SSS origin set at 4 cm below LSO. (d) tSSS results obtained with the SSS
origin set at LSO. (e) tSSS results obtained with the SSS origin set at 4 cm above LSO. (f) tSSS
results obtained with the SSS origin set at 8 cm above LSO. The first and the second rows show theselected MEG sensors time courses from the left and the right hemispheres, respectively. The sensor
topographical maps are shown in the third row. The vertical red and green lines indicate spikes thatwere identified in the DSSP-processed data but not identified in the original data. The third row
shows maps at the same time point indicated by the red line.
30

0.9
0.85
0.8

Correlation

0.75
0.7
0.65
0.6
0.55
0.5

LSO-4cmLSO

LSO+4cm

LSO+8cm

Figure 9: Correlation between DSSP and tSSS processed sensor time courses from resting-state MEG
data. Average correlation with standard error bars are shown across ten patients with respect to thelocation of the SSS origin relative to LSO.

31