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

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 surgical treatment of epilepsy. Here we compare the performance of two artifact removal algorithms for MEG data: dual signal subspace projection (DSSP) and temporally extended signal space separation (tSSS). 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 in patients with VNS implants. Main results. DSSP and tSSS algorithms have a similar ability to reject interference in both simulated and real 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-amplitude artifacts from outside the brain. DSSP might be a better choice than tSSS when the 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.


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 for use 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 non-biological sources such as power lines and trains, and from biological sources outside the brain, for example the heart. Although most of this interference is of a similar magnitude to brain activity, some can be quite large, arising from magnetic fields in the environment or even from the patient's own body or head, including, for example, artifacts from dental work. A prototypical example is the interference arising from vagal nerve stimulator (VNS) implants in patients with intractable epilepsy-VNS artifacts make it very difficult for us to consistently and reliably detect and localize epileptiform activity in such patients [1][2][3].
The preferred way to deal with interference from the environ ment is to prevent it. If this is not possible a variety of methods have been applied to minimize artifacts 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, so filtering may not be effective. One can also reject trials 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 behind this 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 the removal of artifacts. However, these methods require users to make subjective choices during application (e.g. choice of threshold in PCA; selection of interference components 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 (tSSS) method was developed based on spherical harmonic decomposition, and offers a potential solution [4]. However, the validity of this method has only been demonstrated for Elekta and has not been shown for other hardware platforms. Recently, we proposed a MEG hardware platform-independent algorithm named dual signal subspace projection (DSSP) for removal of large interference in 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 identification of epileptic spikes easier and improving spike mapping [3,9].
The performance of DSSP and tSSS has not been directly compared. In this study we compare the performance of these two algorithms using simulated data and real MEG recordings from patients with VNS implants. We test the algorithms' ability to remove interference from resting-state MEG data in patients with VNS implants and also to improve source localization of somatosensory evoked fields (SEFs) using Champagne, an empirical Bayesian source reconstruction algorithm described previously [10]. Results show that DSSP and tSSS algorithms have a 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 performance comparison, 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.

Structure of the DSSP and tSSS algorithms
This section briefly describes the DSSP and tSSS algorithms. Detailed explanations of these algorithms have been published [2,3,11]. These algorithms share a common structure outlined in figure 1. As shown in this figure, these algorithms consist of three steps: signal nulling, estimation of time-domain interference 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 spatial-domain SSP algorithm [12,13], 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 signal nulling in this paper, and the operator used for it is called the signal nulling operator.
Theoretically, the noise subspace projector is an ideal signal nulling operator. However, since the signal and noise subspace projectors are unknown, we must use something that can substitute for these projectors. For this, the DSSP algorithm uses the pseudo-signal subspace projector P S defined in equation (A.4) in the appendix. The tSSS algorithm uses the SSS internal subspace projector P int defined in equation (B.19), the derivation of which is described in detail in appendix B.
When using P int , the internal region needs to be matched with the source space 7 , and this matching between the internal region and the source space can be attained by setting the origin for the spherical harmonics decomposition at a proper location. Note that the location of this origin is a parameter that should be specified by the user 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 the interference removal results of the tSSS algorithm, as shown in our computer simulation in section 3.
Note that the only practical difference between the DSSP and tSSS algorithms is in the manner of approximating the signal subspace projector in the first step; the subsequent steps are exactly the same for these algorithms. The second step estimates the interference subspace by computing the intersection 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.

Estimation of interference subspace and interference removal
Regarding interference removal, the model for the measured data is defined in a matrix form such that where B is the data matrix, B S is the signal matrix, B I is the interference matrix and B ε is a sensor noise matrix. The approximate signal subspace projector is denoted by P S , which is equal to P S in the DSSP algorithm and P int in the tSSS algorithm. This P S satisfies the relationship P S B S ≈ B S and the signal nulling operator is obtained as I − P S . The both DSSP and tSSS algorithms modify the measured data matrix in two ways by applying P S and I − P S to B, resulting in where B in and B out are two kinds of modified data matrices obtained from B. In equations (2) and (3), we can observe that B I -related terms exist as common components in B in and B out . The row spaces of P S B I and (I − P S )B I are both equal to the row space of B I , which is equal to the time-domain interference subspace K I . 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 interference subspace is obtained as The basis vectors of the intersection can be derived using the algorithm described in [11,14]. Once the orthonormal basis vectors of the intersection, ψ 1 , . . . , ψ r , are obtained, we can compute the projector onto the intersection Π X such that Using Π X as the projector onto the interference subspace K I , the interference removal (and as a resultant the estimation of the signal matrix) is achieved by applying the time-domain SSP: where B S is the estimated signal matrix.

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 is a 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 other in which it is not properly chosen. We then compare the results of the DSSP algorithm with those of the tSSS algorithm obtained with proper location of the SSS origin. In our computer simulation, a sensor alignment of the 275-channel whole-head sensor array from the Omega TM (CTF MEG Technology, Coquitlam, BC, Canada) neuromagnetometer was used. The coordinate system and source-sensor configuration used in the computer simulation are depicted in figure 2(a). We assumed a single source at a location shown by a (a black square in (a) and an orange square in (b)) in this figure; its location was presumably in the right parietal area near the primary somatosensory cortex, representing a source located near the border of the source space.
We assumed two locations of the SSS origin. These locations are marked as 'A' and 'B' in figure 2(a), and referred to as location A and location B in this section. Location A is 10 cm below the vertex of the sensor array, and approximately 4 cm above a typical location of the origin of a single 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. Location B is 16 cm below the vertex of the sensor array, and 2 cm below LSO. These two locations are shown by red filled circles in figure 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 G S (r) is defined such that where y S (r) is the signal vector due to a single source at r. The distribution of G S (r) was computed by scanning a single source over the whole source space. The plots of G S (r) are shown in figure 2(b). In this figure, the left-hand, middle and right-hand columns, respectively, show maximum-intensity projections of the three-dimensional distribution of G S (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 G S (r), and the color of the innermost region (dark blue) indicates the value of G S (r) less than 1/100. In figure 2(b) the top row shows the the signal gain mapping obtained using the pseudo signal subspace projector P S (equation (A.4)) used in the DSSP algorithm. 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, which suggests that location A can be an appropriate choice for the SSS origin location. However, when the origin was set at location B, the source was located slightly outside the signal nulling region, which suggests 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 figure 2(a)). To do this, 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 [15]. The generated signal data B S with 2400 time points are shown in figure 3(a).
The top and the bottom panels in figure 3(b), respectively, show B in = P S B S and B out = (I − P S )B S for the DSSP case where P S is used as P S . Here, the bottom panel shows the results of 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 in which 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 (top panel) 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 remains in the bottom panel. These results can be predicted from the signalgain plots in the bottom panel of figure 2(b) where the source was outside the signal nulling region.
Next, we performed interference removal experiments. We put a single interference source just outside 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-pass-filtered random time course. The measured data were generated by adding the interference data to the signal data (shown in figure 3(a)) with an interference to signal ratio (ISR) equal to 10. The generated interference-overlapped data are shown in figure 4(a). Here, since the interference magnetic field was ten   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 figure 4(b). The field contour map of these results at 1200 ms is shown in the top right panel of figure 4(c). Note that the ground truth of the sensor time courses is shown in figure 3(a). The ground truth of the field contour map (obtained from the signal magnetic field) is shown in the top left panel of figure 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 figure 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 of figure 4(c). These tSSS results are almost the same as the ground truth results, indicating that the tSSS algorithm effectively removed the interference.
The bottom panel in figure 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 figure 4(c). Although these time courses suggest that the interference was removed pretty well, the field map indicates that some amount of distortion was caused in the interferenceremoved results. These results in figures 3 and 4 confirm that location A represents a proper location of the SSS origin while 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 comp onents 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 conclude that DSSP and tSSS algorithms have almost the same performance for interference removal. However, in order for 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 not properly chosen, a mismatch between the internal region and the source space can distort the signal magnetic field.

Collecting MEG data and general aspects of analysis
We chose ten epilepsy patients with VNS implants who underwent a clinical MEG study as part of evaluation for epilepsy surgery at the University of California, San Francisco (UCSF) Biomagnetic Imaging Laboratory (BIL) between 24 November 2004 and 6 May 2016. Prior to the MEG scan, all patients had high-resolution epilepsy protocol 3T T1-MRI scans for co-registration of localized spikes. MEG recordings were performed inside a magnetically shielded room with a CTF 275-channel whole-head axial gradiometer system.
We collected two kinds of MEG data from each patient: somatosensory evoked MEG and resting-state spontaneous MEG. Somatosensory evoked data sets (sensory evoked field paradigm) were collected with the sample rate set at 1.2 kHz. A patient's 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 this case, the left) somatosensory cortical area for the hand, i.e. the dorsal region of the postcentral gyrus.
The leadfield for each subject was calculated in NUTMEG 4.1 [16] with an 8 mm voxel grid resulting in the generation of approximately 5300 voxels. 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 the 0 ms time point.
Source localization of the somatosensory data was conducted using a united Bayesian framework for MEG/EEG source imaging that includes variational Bayes factor analysis (VBFA) for noise learning [17] and a sparse Bayesian algorithm (Champagne) for source localization [10,18,19]. 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 passband of 0-75 Hz. 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 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 in our previous 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 s obtained at a sampling rate of 1.2 kHz. For epileptic spike data, we used a 10 s time-window for DSSP and tSSS analysis at a sampling rate for our data of 600 Hz.

VNS artifact removal in somatosensory evoked MEG
Somatosensory MEG recordings from selected channels are shown for four representative patients in figure 5. The original MEG recordings (before artifact removal) are shown in figure 5(a). These original recordings feature time courses with partial periodicity, low frequency and high amplitude, typical of VNS artifacts. The results of DSSP artifact removal are shown in figure 5(b). In these results, the periodic 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 of four different locations: 4 cm below LSO, LSO, 4 cm above LSO and 8 cm above LSO 8 . The results of tSSS artifact removal using each of these four SSS origin locations are shown in figures 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.
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 figure 6. Results from the original data are shown in figure 6(a), in which the source localization obviously failed for all four patients. Results from the DSSP processed data are shown in figure 6(b). We can see that the sources were consistently localized near the primary somatosensory area in these results. Results from tSSS processed data are shown in figures 6(c)-(f). Results similar to those from DSSP application were obtained using the tSSS algorithm with the SSS origin set at 4 cm above LSO (figure 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.   figure 5. Each row 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.
The results of source localization experiments for ten patients are summarized in figure 7 where the bar 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 peak maximum was within a region with a label of 'Postcentral gyrus' on the standardized brain atlas, we considered the source localization successful.
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 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 of successful localizations was no more than eight. The results for only one subject were localized using the original VNSoverlapped MEG sensor data.
In summary, the results from these somatosensory experiments indicate that DSSP and tSSS algorithms have nearly the same performance for VNS artifact removal. However, when implementing the tSSS 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 the tSSS algorithm can be significantly worse than that of the DSSP algorithm. Also, the results in figure 7 clearly indicate that the proper choice of the SSS origin location is in the vicinity of 4 cm above LSO.

Artifact removal in resting-state MEG
We next performed VNS artifact removal experiments using resting-state MEG. In figure 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 figure 8(b) in which the periodic feature is greatly diminished and the background looks similar to that of most patients with epilepsy and without a VNS device.
According to the results of somatosensory experiments, the proper location of the SSS origin is around 4 cm above LSO. In order to confirm if this is also true for resting-state MEG in epilepsy we again 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 figures 8(c)-(f).
As can be seen in figures 8(c) and (d), when the SSS origin was set at 4 cm below LSO or at LSO, the tSSS 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 row of figure 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 those from the DSSP algorithm.
We performed visual identification of epileptiform activity (spikes) in these data sets. The vertical red lines indicate spikes that were identified in the DSSP-processed data but not 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-and tSSS-processed resting-state sensor time courses. In these resting-state measurements, each subject's 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 figure 9 with respect to the SSS origin location. As shown here, the correlation coefficient attains its maximum value when the SSS origin location is set at 4 cm above the LSO, 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 their 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.

Discussion
In this paper we compare the DSSP and tSSS algorithms for MEG artifact removal. We conducted computer simulation and experiments using real MEG data 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 it is 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 the subject's head MRI. Naturally, the question may arise whether this proper location depends on the shape of the sensor helmet and sensor arrangement, implying that it might change across different MEG hardware 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,5,20]. The PFA algorithm involves learning 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 both interference 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 properties of 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 an active period but not during a baseline measurement period. Also, these algorithms 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 removal of some types of physiological artifacts common in MEG. Indeed, DSSP and tSSS may be used for removal of cardiac, ocular or dental artifacts. In our research group, we have successfully deployed DSSP and tSSS for all these sources of artifacts in MEG data. The algorithm performance of both DSSP and tSSS requires some degree of spatial separation of the sources of artifacts from the sources of 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 for removal of cardiac and dental sources than for ocular artifacts (data not shown). Some types of ICA algorithms are known to be effective for the removal of such artifacts. Nevertheless, our future plan includes the performance comparison of these two algorithms with ICA algorithms, because ICA algorithms require subjective selection of interference components but DSSP/tSSS algorithms do not require such procedure.

Conclusion
Both DSSP and tSSS are able to remove large-amplitude artifacts from outside the brain. DSSP can be a better option than tSSS when the origin location for tSSS is not easy to choose. Both algorithms can recover distorted MEG recordings from people with intractable epilepsy and VNS implants, improving epileptic spike identification and source localization of functional activity and epileptiform activity.
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 and the pseudo-signal subspace Ȇ S is defined such that 2) The true signal subspace is defined as where r q 1 , . . . , r q q are the source locations assuming that there are q sources in the source space. If the voxel interval is sufficiently small and voxel discretization errors are negligible, each of the source locations r q 1 , . . . , r q q is equal to one of the voxel locations, and the space spanned by {{L(r q 1 ), . . . , L(r q q )} is a subset of the space spanned by {L(r 1 ), . . . , L(r N )}, resulting in the relationship Ȇ S ⊃ E S .
Let us derive the orthonormal basis vectors of the pseudosignal subspace. To do this, we compute the singular value decomposition of F: If the singular values λ 1 , . . . , λ τ are distinctively large and other singular values λ τ +1 , . . . , λ M are nearly equal to zero, the leading τ singular vectors e 1 , . . . , e τ form orthonormal basis vectors of the pseudo-signal subspace Ȇ S [21]. Thus, the projector onto Ȇ S is obtained using 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 sourcefree 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 where µ 0 indicates the magnetic permeability of free space. In equation (B.1), ν ,m (θ, φ) and ω ,m (θ, φ) are the modified vector spherical harmonics [22,23]. In the right-hand side of equation (B.1), the first term represents the magnetic field generated from sources 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 the sensors is referred to as the external region. Note here that the origin location of the spherical polar coordinate r = (r, θ, φ) is a user-defined parameter, and the 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 the internal region but an improper choice may cause a situation where some signal sources are located outside 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 this, the magnetic signal detected by the j th sensor is denoted by b j , and the location and the normal vector of the j th sensor by r j and ζ j . Then, we have b j = B(r j ) · ζ j = b j int + b j ext , (B.2) where the notation '·' indicates taking the inner product between two vectors. Here, b j int and b j ext respectively represent magnetic components originating from the internal and external regions. These components are expressed as where column vectors c ,m and d ,m are given by . . . with the data vector b. That is, the matrix P int acts as a projector that passes the internal components and blocks the external ones 9 . Therefore, we call P int the SSS internal subspace projector in this paper. In exactly the same manner, the SSS external subspace projector P ext is derived as This P ext passes the external components but blocks the internal ones. It can be seen from equations (B.19) and (B.20) 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 the matrix C, and P int does not change the vector if it is already in the internal subspace. Therefore, the application 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 the internal region and the source space and some signal sources are located outside the internal region, the signal nulling operator I − P int only removes that part of B S that belongs to the column space of C.