Improving the Nulling Beamformer Using Subspace Suppression

Abstract

Magnetoencephalography (MEG) captures the magnetic fields generated by neuronal current sources with sensors outside the head. In MEG analysis these current sources are estimated from the measured data to identify the locations and time courses of neural activity. Since there is no unique solution to this so-called inverse problem, multiple source estimation techniques have been developed. The nulling beamformer (NB), a modified form of the linearly constrained minimum variance (LCMV) beamformer, is specifically used in the process of inferring interregional interactions and is designed to eliminate shared signal contributions, or cross-talk, between regions of interest (ROIs) that would otherwise interfere with the connectivity analyses. The nulling beamformer applies the truncated singular value decomposition (TSVD) to remove small signal contributions from a ROI to the sensor signals. However, ROIs with strong crosstalk will have high separating power in the weaker components, which may be removed by the TSVD operation. To address this issue we propose a new method, the nulling beamformer with subspace suppression (NBSS). This method, controlled by a tuning parameter, reweights the singular values of the gain matrix mapping from source to sensor space such that components with high overlap are reduced. By doing so, we are able to measure signals between nearby source locations with limited cross-talk interference, allowing for reliable cortical connectivity analysis between them. In two simulations, we demonstrated that NBSS reduces cross-talk while retaining ROIs’ signal power, and has higher separating power than both the minimum norm estimate and the nulling beamformer without subspace suppression. We also showed that NBSS successfully localized the auditory M100 event-related field in primary auditory cortex, measured from a subject undergoing an auditory localizer task, and suppressed cross-talk in a nearby region in the superior temporal sulcus.

Introduction

Magnetoencephalography (MEG) is a functional imaging modality that provides accurate temporal (~1 ms) and spatial (~1–2 cm) measures of cortical activity. MEG is extensively used in cognitive neuroscience research (Hari and Salmelin, 2012; Salti et al., 2015; Mamashli et al., 2017; Pantazis et al., 2017) and it is approved for some clinical studies, such as localization of epileptic foci (French et al., 1993; Stufflebeam et al., 2009).

In MEG, a source estimation technique is needed to translate the measured signals to estimates of the underlying neural current sources. MEG source estimation is an ill-posed inverse problem (Wang et al., 1993; Wendel et al., 2009): there is no unique solution and the solutions available are sensitive to noise. Therefore, restricting assumptions about the possible sources are needed and noise sensitivity must be mitigated by regularization. A fast, whole-brain source-localization technique is the minimum norm estimate (MNE), which is a spatially smooth estimate with a relatively wide point spread and strong crosstalk between regions of interest (ROIs) (Dale and Sereno, 1993; Hämäläinen et al., 1993; Wang et al., 1993; Hämäläinen and Ilmoniemi, 1994; Dale et al., 2000). Due to its fast computation time MNE is well-suited to conduct exploratory studies.

Another class of methods, called beamformers, have become increasingly popular for MEG and EEG source estimation, due to their goal to minimize cross-talk. The beamforming method was originally formulated for radar and sonar to pick up signals from specific transmitters while attenuating those from other locations (Van Veen and Buckley, 1988; Hillebrand et al., 2005). Beamformers have since been used in numerous applications including wireless communication, astrophysics, and biomedical signal processing (Baillet et al., 2001; Hillebrand et al., 2005). Broadly, the beamforming involves an array of sensors and combines the signal recorded by each sensor to increase the signal/noise ratio and focuses the entire sensory array on a central spatial location. In brain imaging, MEG and EEG beamformer methods are used to measure the location, magnitude and direction of the magnetic fields resulting from the electric currents flowing inside the brain (MEG) or the electric potential on the scalp surface (EEG).

Although several flavors of beamformer algorithms exist, we are specifically focusing on those algorithms that are closely related to the specific method we are proposing, the nulling beamformer with subspace suppression (NBSS). The Linearly Constrained Minimum Variance (LCMV) beamformer is frequently used in MEG data analysis because it is explicitly aimed at suppressing cross-talk (Van Veen and Buckley, 1988; Spencer et al., 1992; Van Veen et al., 1997; Robinson and Vrba, 1999; Hui et al., 2010). Simulations have shown that the LCMV method results in greater spatial accuracy compared to the MNE method (Hadjipapas et al., 2005). The LCMV method, like other adaptive beamformer methods, makes the assumption that the sources are uncorrelated (Reddy et al., 1987; Hillebrand et al., 2005). However, due to cortico-cortical communication which occurrs even at rest, the signals are invariably correlated across cortical regions (Reddy et al., 1987; Hillebrand et al., 2005).

Several methods have been recently developed for addressing correlated sources. For example, the dual-core beamformer specifically isolates pairs of correlated sources (Brookes et al., 2007; Diwakar et al., 2011). The nulling beamformer (NB) method introduces additional constraints on the LCMV beamformer to suppress contributions from more than two cortical sites (Hui et al., 2010). In the NB method, a predefined set of ROIs is used and it is assumed that all the sources originate from these ROIs. This assumption is suitable for hypothesis-driven studies, where the goal is to compute activity in the ROIs hypothesized to be involved in the neural implementation of a specific task.

The nulling beamformer requires an additional step to avoid overconstraining the solution: the gain matrices mapping the sources in the ROIs to sensor signals are reduced in rank through truncated singular value decomposition (TSVD), which retains a specified number of the largest singular values. However, the TSVD does not address the problem of source separation, and thus the resulting components may be strongly correlated with those of neighboring ROIs.

We propose a novel beamformer method, NBSS, which aims to modify the mapping from sensor space to source space in order to suppress high-crosstalk components. By doing so, we are able to reduce cross-talk between nearby sources for cortical connectivity analysis that has reduced contamination from signal interference. The NBSS method uses a tuning parameter to control the trade-off between the maximal amount of cross-talk and the signal power of an ROI. In two simulations we compared the performance of MNE, NB, and NBSS to measure cross-talk among predefined ROIs. The outcome of the simulations show that NBSS outperformed NB and MNE in reducing cross-talk while retaining signal power.

Methods

In this section we will briefly define the MNE, the nulling beamformer (NB), and the NBSS method that we propose.

Minimum Norm Estimate

The measured MEG signal vector y(t) at time t is assumed to be the sum of contributions of all neural sources:

y(t)=∑gixi(t)+n(t)=Gx(t)+n(t), (1) y(t)=∑gixi(t)+n(t)=Gx(t)+n(t), (1)

where xi(t) is the time course of the ith source, x(t) is the collection of time courses into a vector, gi is the gain vector (column vector of the forward matrix G) for that source, and n(t) is additive noise. The MNE, xMNE(t), is the linear inverse estimate:

xMNE(t)=RGT(GRGT+C)−1y(t), (2) xMNE(t)=RGT(GRGT+C)-1y(t), (2)

where R and C are the source and noise covariance matrices, respectively (Dale et al., 2000).

Nulling Beamformer

The nulling beamformer (NB) method is based on the LCMV beamformer (Van Veen et al., 1997). This model assumes N distinct sources, which can be discrete current dipoles or cortical patches with uniform activity. The aim is to find a linear operator wi that, when applied to y(t), yields unit gain for each source i while the contribution from other sources is minimized. The nulling beamformer introduces the additional constraint that attempts to eliminate the effects of all other confounding sources from the source location of interest (Hui et al., 2010). Thus, for each source, a constraint is applied in order to remove the signal contribution from other sources. The set of constraints for large number of ROIs will overconstrain the system. Therefore, the sources must be limited to a small number P of ROIs. To further reduce the number of constraints, a TSVD is used to obtain a rank L approximation for the gain matrix Gp for each ROI patch, indexed by p:

GLp=ULpSLp(VLp)T, (3) GpL=UpLSpL(VpL)T, (3)

where SLp SpL = diag(σp, 1, …, σp, L) is a diagonal matrix of largest singular values of Gp and the columns of ULp UpL and VLp VpL contain the first L left and right singular vectors, respectively.

With the TSVD applied to each cortical ROI we implemented the nulling beamformer method as defined in (Hui et al., 2010). The following matrices for P cortical patches were constructed:

H=[UL1,…,ULP] (4) H=[U1L,…,UPL] (4)

Bp=[r1dVL1SL1−1,…,rNdVLPSLP−1] (5) Bp=[rd1V1LS1L-1,…,rdNVPLSPL-1] (5)

where rid=1 rdi=1 for i = p and rid=0 rdi=0 for i ≠ p. Using these matrices, the following weight vector was constructed for each patch p:

wp=C−1H(HTC−1H)−1BTp (6) wp=C-1H(HTC-1H)-1BpT (6)

where C is the covariance matrix of the measurement data. The nulling beamformer source timecourses x^NBp(t) at the pth ROI is then:

xNBp(t)=wTpy(t). (7) xpNB(t)=wpTy(t). (7)

TSVD of Nulling Beamformer May Remove Separable ROI Activity

The nulling beamformer works best when the projections between any pair of ROI gain matrices Gp and Gq are small, that is they correspond to low cross talk. However, if this condition does not hold, such as when the distance between the ROIs is small, the nulling beamformer's nulling constraint may mutually nullify large contributions of each ROI's signal. Distant ROIs will have small magnitude components that overlap which are truncated through the TSVD.

Using “Subspace Suppression” in Place of TSVD

We propose a method, subspace suppression (SubS), which will suppress singular values of the gain matrices whose associated basis vectors have high cross talk. To provide greater suppression of crosstalk we couple this method with the nulling beamformer. In the NBSS method, we first compute the projections onto a particular singular value in an ROI's gain matrix from all other basis vectors from the other ROI gain matrices. We compute for the amplitude of the projections from the other ROIs (indexed by q) onto the pth ROI's singular value's basis vectors:

σq→p,l=∑k=1K∣∣∣∣uTp,kGq,lvp,k∣∣∣ (8) σq→p,l=∑k=1K|up,kTGq,lvp,k| (8)

K refers to the number of singular values in G while k indexes the rows in Up and Vp, which are the left and right singular value matrices of the pth ROI, up, k and vp, k are the rows of U and V respectively, and Gq,l=ULqSl,Lq(VLq)T Gq,l=UqLSql,L(VqL)T where Sl,Lq Sql,L is the singular value matrix SLq SqL with all elements set to zero except for the lth singular value along the diagonal. The amount of crosstalk contribution from all ROIs q to ROI p is the summation of the following variables:

σCT,p,l=∑q=1,q≠pPσq→p,l=∑q=1,q≠pP∑k=1K∣∣uTp,kGq,lvp,k∣∣ (9) σCT,p,l=∑q=1,q≠pPσq→p,l=∑q=1,q≠pP∑k=1K|up,kTGq,lvp,k| (9)

This value is normalized through dividing by the summation of cross-talk contributions among all ROIs:

σˆCT,p,l=σCT,p,l∑Ll=1σCT,p,l (10) σ^CT,p,l=σCT,p,l∑l=1LσCT,p,l (10)

The goal of developing the NBSS method is to control the amount of cross-talk through a tuning parameter. To have a tuning parameter that does not have arbitrary limits, we placed the following design constraints on the construction of the NBSS method: the tuning parameter set to zero should not have any effect on the gain matrix, the tuning parameter set to one should completely nullify all crosstalk contributions, and the amount of suppression of a singular value should increase as the amount of cross-talk increases for a fixed value of the tuning parameter. A function that can satisfy the aforementioned properties is:

σNBSSp,l=σp,l1−α1+αβ (11) σp,lNBSS=σp,l1-α1+αβ (11)

where σp, l is the lth singular value of gain matrix Gp, α is the tuning parameter, and β is chosen based on the amount of suppression of the singular desired values. The β value is set such that when α is 0.5, the suppression will equal to one minus proportional amount of cross-talk in the area reflecting the amount of suppression desired at that level. The rule we defined is summarized in Equation (12):

1−α1+αβ=1−σˆCT,p,l∣∣α=0.5 (12) 1-α1+αβ=1-σ^CT,p,l|α=0.5 (12)

We obtain the following expression for the singular values of the new gain matrix, by solving (12) for β and inserting the β into (11):

σNBSSp,l=σp,l1−α1+α⎛⎝⎜1−2σˆCT,p,lσˆCT,p,l−1⎞⎠⎟ (13) σp,lNBSS=σp,l1-α1+α(1-2σ^CT,p,lσ^CT,p,l-1) (13)

As σˆCT,p,l σ^CT,p,l increases, the amount of suppression for the same α also increases. The scaling factor α is a parameter that can range from 0 to 1, with 0 applying no cross-talk suppression, and 1 applying maximal suppression of all overlapping cross-talk components. How scaling varies as a function of α and σˆCT,p,l σ^CT,p,l is illustrated in Figure 1 which shows the plot of the scaling factor on σp, l as a function of α for varying σˆCT,p,l σ^CT,p,l .

FIGURE 1

Figure 1. Scaling factor vs. Tuning parameter and normalized cross-talk. Scaling factor of singular values of the gain matrix of ROI p, σp, l, as a function of the tuning parameter, α, and the normalized cross-talk contribution, σˆCT,p,l σ^CT,p,l . There is no suppression when α is 0 and complete suppression when α is 1.

Nulling Beamformer With Subspace Suppression

We introduced the subspace suppression method into the nulling beamformer procedure by replacing the original singular value matrix in the TSVD in Equation (3) with the reweighted singular value matrix whose singular values are computed in Equation (13) to obtain:

GNBSSp=ULpSNBSSpVL,Tp (14) GpNBSS=UpLSpNBSSVpL,T (14)

The singular values are preweighted to minimize cross-talk between ROIs, thus allowing higher separation ability through the nulling beamformer solution without mutual signal nullification for high degrees of overlapping singular vectors.

Simulations

In each of the simulations we defined a set of ROIs used to test two factors: the effect of various factors on the nulling constraint between two ROIs (Simulation 1), and the effect of multiple ROIs in a realistic MEG scenario (Simulation 2).

General Method for Simulation of Cortical Signals

An MEG measurement with the 306-channel Elekta-Neuromag Vectorview system was assumed with a realistic location and alignment of the sensor array with the subject's head, obtained from an actual measurement (Vaina et al., 2010). The forward solution was computed using a single-compartment boundary-element model (Hamalainen and Sarvas, 1989). A spatial noise-covariance matrix was estimated from the MEG data in the magnetically shielded room without a subject and was used to spatially whiten both the data and the gain matrices. The noise covariance was used to add spatially colored noise to the simulated signals. We generated 1,000 samples of random Gaussian noise in the sensor space as a baseline level with no signal. The amplitude of signals from each ROI was chosen so that it would have a relative noise amplitude (RNA) of 0.1. The RNA is the ratio between the standard deviation of the noise divided by the standard deviation of the signal (Dannhauer et al., 2013). Additional 1,000 samples were generated and added to simulated activation mapped from the ROI to the sensor space.

Using the procedures detailed in the Theory section, we estimated the ROI signals from the MNE, nulling beamformer, and NBSS methods. In the computation of the MNE and beamformer methods the orientations of the sources were not constrained. The ROI signal estimate was obtained by averaging the power of the dipole amplitudes across the three dipole components and locations within each ROI. We chose the power measured from the ROI that generated the signal to be 10 times the noise level from the MNE estimate to provide a baseline of comparison between methods. For the TSVD, we vary L from 1 to 5.

To assess the amount of cross-talk between two ROIs, we chose a criterion often used in wireless communications to determine the quality of cross-talk suppression between channels, the signal-to-interference ratio (SIR) (Andersin et al., 1998; Jeske and Sampath, 2001). We expressed SIR as the log ratio of the power of the source at a location to the power of interfering sources at that location:

SIRp=10log10Pp→p∑Pq≠pPq→p (15) SIRp=10log10Pp→p∑q≠pPPq→p (15)

where Pq→p is the power in ROI p due to signal generated in ROI q, which we estimate through simulations. At a level of 0, the SIR has the same power as other ROIs. We will use the SIR measure to compare the MNE, NB, and NBSS methods.

Comparisons of the Minimum Norm Estimate, Nulling Beamformer, and Nulling Beamformer With Subspace Supression Methods

We applied the three methods, MNE, NB, and NBSS methods, to two simulations. In the first simulation, we tested the effects of different parameters on the SIR in the three methods. We sampled a vertex (the reference) in the intraparietal sulcus (IPS) and a random vertex (the test vertex) at a random distance away. The test vertex was allowed to be from 2 to 10 cm away at a 2 cm interval with a random jitter of 0.1 cm. For each vertex, we defined a spherical region of 1 cm radius. Using the MNE, nulling beamformer (NB), and the NBSS method we proposed, we calculated SIR between the reference vertex and the test vertices. In the NBSS method, the tuning parameter alpha varied from 0.05 to 0.95 in 0.05 steps. Subsequently, we chose the alpha parameter with maximal SIR. We repeated the simulation 100 times, randomly choosing new reference vertices within IPS. This produced 100 simulated “trials” for each sampled vertex point.

In the second simulation we used a set of cortical ROIs activated in one of our recent MEG study on a cognitive task on visual and auditory motion (Vaina et al., 2010). We investigated four right-hemisphere ROIs: MT+, superior temporal polysensory area (STP), primary auditory area (aud), and the dorso-lateral prefrontal cortex (DLPFC) (Figure 2). We compared the cross-talk with SIR as resulted from using the MNE, NB, and NBSS methods. We computed 1000 simulated “trials” of independent single samples generated from each of the regions.

FIGURE 2

Figure 2. Regions of Interest on right hemisphere cortical surface. Four regions-of-interest (ROIs, green) on the right hemisphere cortical surface of a single subject (adapted from [17]): MT+, superior temporal sulcus (STP), primary auditory area (aud), and dorsolateral prefrontal cortex (DLPFC).

In the (Vaina et al., 2010) study, we measured MEG data from a subject undergoing an auditory localization task. We presented a pure 300 Hz auditory tone to both ears at 69 dB for 1 s. The auditory tone induces a peak in the primary auditory cortex at 100 ms from auditory onset, labeled the M100 evoked field (Reite et al., 1978). Our interest was to localize the signal in the primary auditory area (aud) and minimize cross-talk in the neighboring area STP. Therefore, we reconstructed the signal in the primary auditory (aud) and STP regions using the MNE, NB, and NBSS methods and compared the amplitude at the M100 peak.

Results and Discussion

Simulation 1

Our first simulation investigated how the SIR is affected by distance between two ROIs as distance between them is increased. Figure 3 shows that the nulling beamformer (NB) is significantly greater in SIR than the MNE, even for the 2 cm distance [t(198) = 54.48, p < 0.001]. The NBSS provides higher SIR values at all distances [t-test at 2 cm between NB and NBSS, t(198) = 3.36, p = 0.001].

FIGURE 3

Figure 3. Signal to Interference Ratio (SIR) between 1 cm radius regions of varying distance. Error bars represent 1 standard error over 100 repetitions.

For the shortest distance that we tested (2 cm), we varied the effect of the truncation amount of the TSVD for 1 through 5 dimensions (Figure 4). With L = 1, there is no difference between NB and NBSS methods, since NBSS relies on modulating the gain between components. However, in this situation, both NB and NBSS (M = 10.26, SD = 2.59) have higher SIR than MNE (M = 5.47, SD = 1.77) [t(198) = 15.15, p < 0.001]. However, as the number of components increased, the effectiveness of the NB and NBSS methods decreased after L = 2. At L = 5, SIR with NBSS (M = 5.45, SD = 2.7) was lower than in MNE (M = 4.43, SD = 4.46) [t(198) = 1.98, p = 0.049]. The decrease of SIR for L > 2 may be due to increased signal overlap between the two ROIs.

FIGURE 4

Figure 4. Signal to Interference Ratio (SIR) between 1 cm radius regions of varying truncation levels in NB and NBSS methods. Error bars represent 1 standard error over 100 repetitions.

We fixed the truncation amount to L = 3 and varied the residual noise amplitude (Figure 5). The SIR difference [t(198) = 6.12, p < 0.001]. However, with an RNA level of 10 where the signal is below the noise level on each trial, the NB (M = 1.71, SD = 1.39) and NBSS (M = 2.44, SD = 3.22) methods have significantly higher SIR than MNE (M = 0.81, SD = 1.40) [between NB and MNE t(198) = 2.54, p = 0.012]. NBSS has slightly higher SIR than NB, but is not significant [t(198) = 1.61, p = 0.110].

FIGURE 5

Figure 5. Signal to Interference Ratio (SIR) between 1 cm radius regions of varying Residual Noise Amplitudes. Error bars represent 1 standard error over 100 repetitions.

In Figure 6, we show that the SIR is affected by the size of the ROI (Figure 6). We fixed the distance to 2 cm from the center of the ROI and compare a 0.5 cm radius to a 1 cm radius. MNE signal power dropped significantly with a larger ROI radius (radius = 0.5 cm, M = 10.67, SD = 0.22; radius = 1 cm, M = 6.83, SD = 2.17), since the ROIs were physically closer and therefore spreading more signal into the neighboring area [t(198) = 18.20, p < 0.001]. However, since NB and NBSS methods selectively choose signal components that best suppress cross-talk, SIR was improved with larger radius, due to having more signal components that can be manipulated by reweighting singular values using NBSS [e.g., for NBSS at radius = 0.5 cm, M = 11.38, SD = 6.51; radius = 1 cm, M = 25.43, SD = 3.01, t(198) = 19.64, p < 0.001].

FIGURE 6

Figure 6. Signal to Interference Ratio (SIR) between regions of fixed 2 cm distance of 0.5 cm and 1 cm radii. Error bars represent 1 standard error over 100 repetitions.

Simulation 2

The second simulation modeled a more realistic scenario using multiple ROIs. Figure 7, shows the SIR between the MNE, NB, and NBSS methods. The NBSS method outperformed NB and MNE in all cases. The NB method, for most ROIs resulted in higher SIR than the MNE method, except when applied to the aud ROI. This may be due to strong cross-talk from STP, which is close to the aud area. As alpha varied, the SIR increased or decreased in non-continuously. The reweighting procedure in the NBSS method reorders singular values, which are ranked by their corresponding power. Therefore, the reordering drastically changed the projections between sensor and source space when varying alpha.

FIGURE 7

Figure 7. Signal to Interference Ratio (SIR) between four regions. Four regions of interest are defined in simulation 2 (MT+, aud, STP, DLPFC). SIR is measured in each region and with the minimum norm estimate (MNE), nulling beamformer (NB), and nulling beamformer with subspace suppression (NBSS) with alpha tuning parameter ranging from 0.05 to 0.95 at 0.05 intervals.

Using the optimal parameters showed in Figure 7, we estimated activation in the primary auditory cortex (aud) and the STP area over all trials from −0.1 to 0.3 s using each of the 3 methods using the auditory localizer MEG data from one subject in (Vaina et al., 2010). We were interested in measuring the auditory M100 peak, occurring at approximately 100 ms. The results are shown in Figure 8. Although the amplitude in aud using MNE is higher than the others (z = 8.55, p < 0.001), there is significant signal spread into STP (z = 4.51, p < 0.001). NB effectively suppresses the signal in STP (z = 2.78, p = 0.003), but also causes significant loss of signal in aud (z = 2.91, p = 0.002). NBSS recovered the high amplitude signal in aud, similar to the MNE solution (z = 5.96, p < 0.001). The amplitude of the signal in STP, however, was suppressed, but was above significance (z = 2.94, p = 0.002).

FIGURE 8

Figure 8. Signal Amplitude in auditory cortex and STP. Amplitude is normalized such that noise amplitude has consistent height in all figures. Time is relative to the auditory tone onset (at 0 s) therefore including a silent period (−0.1 to 0 s) and auditory tone period (0 to 0.3 s).

Conclusion

The nulling beamformer effectively suppresses cross-talk between distant cortical ROIs. However, for ROIs that are close together, the preliminary step of applying the TSVD may remove components that represent the underlying signal. We showed that that the nulling beamformer significantly suppressed signal power in closely neighboring ROIs, thus allowing the measurement of cortical connectivity between closer sources with reduced interference. Importantly, the subspace suppression method reduces the contributions of highly overlapping gain components and thus it prevents the removal of components that significantly separate the signals between ROIs.

In the simulations described in this paper, we optimized SIR through a discrete sampling over the range of values between 0 and 1. We have also shown that while SIR provides a means of comparing signal strength to the combined cross-talk between regions, instead it may be preferred to minimize the maximal cross-talk. Through this optimization, the largest component of cross-talk would be minimized. In future work, we will formalize the optimization of the alpha parameter through SIR and explore alternative optimization measures.

The approach described in this paper is the improving of source localization during real-time applications. By estimating noise and data covariances, this method can be used to separate the signal in multiple ROIs in real-time. Wu et al. (2012) developed a near maximum-likelihood source estimation that relies on nonlinear optimization through search. However, due to the use of a parameter search, this method cannot be used in real-time application of source estimation.

Author Contributions

KR developed the proposed algorithm, carried out the simulations, and co-wrote the manuscript. MH assisted in development of the algorithm, advised on MEG principles, and assisted with writing the manuscript. LV was involved and oversaw all stages of the paper (development of algorithm, simulations, co-writing manuscript).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Seppo Ahlfors for his insightful comments and suggestions in an earlier version of this manuscript. This research was supported in part by the NSF grant 1545668 to LV.

References

Andersin, M., Mandayam, N. B., and Yates, R. D. (1998). Subspace based estimation of the signal to interference ratio for TDMA cellular systems. Wireless Netw. 4, 241–247. doi: 10.1023/A:1019160207631

Baillet, S., Mosher, J. C., and Leahy, R. M. (2001). Electromagnetic brain mapping. Signal Process. Mag. IEEE 18, 14–30. doi: 10.1109/79.962275

Brookes, M. J., Stevenson, C. M., Barnes, G. R., Hillebrand, A., Simpson, M. I., Francis, S. T., et al. (2007). Beamformer reconstruction of correlated sources using a modified source model. NeuroImage 34, 1454–1465. doi: 10.1016/j.neuroimage.2006.11.012

Dale, A. M., Liu, A. K., Fischl, B. R., Buckner, R. L., Belliveau, J. W., Lewine, J. D., et al. (2000). Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron 26, 55–67. doi: 10.1016/S0896-6273(00)81138-1

Dale, A. M., and Sereno, M. I. (1993). Improved localizadon of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach. J. Cogn. Neurosci. 5, 162–176. doi: 10.1162/jocn.1993.5.2.162

Dannhauer, M., Lämmel, E., Wolters, C. H., and Knösche, T. R. (2013). Spatio-temporal regularization in linear distributed source reconstruction from EEG/MEG: a critical evaluation. Brain Topogr. 26, 229–246. doi: 10.1007/s10548-012-0263-9

Diwakar, M., Huang, M.-X., Srinivasan, R., Harrington, D. L., Robb, A., Angeles, A., et al. (2011). Dual-core Beamformer for obtaining highly correlated neuronal networks in MEG. NeuroImage 54, 253–263. doi: 10.1016/j.neuroimage.2010.07.023

French, J. A., Williamson, P. D., Thadani, V. M., Darcey, T. M., Mattson, R. H., Spencer, S. S., et al. (1993). Characteristics of medial temporal lobe epilepsy: I. Results of history and physical examination. Ann Neurol. 34, 774–780. doi: 10.1002/ana.410340604

Hadjipapas, A., Hillebrand, A., Holliday, I. E., Singh, K. D., and Barnes, G. R. (2005). Assessing interactions of linear and nonlinear neuronal sources using MEG beamformers: a proof of concept. Clin. Neurophysiol. 116, 1300–1313. doi: 10.1016/j.clinph.2005.01.014

Hämäläinen, M. S., Hari, R., Ilmoniemi, R. J., Knuutila, J., and Lounasmaa, O. V. (1993). Magnetoencephalography: theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65, 413–497. doi: 10.1103/RevModPhys.65.413

Hämäläinen, M. S., and Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. doi: 10.1007/BF02512476

Hamalainen, M. S., and Sarvas, J. (1989). Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. Biomed. Eng. IEEE Trans. 36, 165–171. doi: 10.1109/10.16463

Hari, R., and Salmelin, R. (2012). Magnetoencephalography: from SQUIDs to neuroscience: neuroimage 20th anniversary special edition. Neuroimage 61, 386–396. doi: 10.1016/j.neuroimage.2011.11.074

Hillebrand, A., Singh, K. D., Holliday, I. E., Furlong, P. L., and Barnes, G. R. (2005). A new approach to neuroimaging with magnetoencephalography. Hum. Brain Mapp. 25, 199–211. doi: 10.1002/hbm.20102

Hui, H. B., Pantazis, D., Bressler, S. L., and Leahy, R. M. (2010). Identifying true cortical interactions in MEG using the nulling beamformer. NeuroImage 49, 3161–3174. doi: 10.1016/j.neuroimage.2009.10.078

Jeske, D. R., and Sampath, A. (2001). “Signal-to-interference ratio estimation based on decision feedback,” in Vehicular Technology Conference, VTC 2001 Spring., IEEE VTS 53rd: IEEE (Rhodes), 2484–2488.

Mamashli, F., Khan, S., Bharadwaj, H., Michmizos, K., Ganesan, S., Garel, K. A., et al. (2017). Auditory processing in noise is associated with complex patterns of disrupted functional connectivity in autism spectrum disorder. Autism Res. 10, 631–647. doi: 10.1002/aur.1714

Pantazis, D., Fang, M., Qin, S., Mohsenzadeh, Y., Li, Q., and Cichy, R. M. (2017). Decoding the orientation of contrast edges from MEG evoked and induced responses. Neuroimage doi: 10.1016/j.neuroimage.2017.07.022. [Epub ahead of print].

Reddy, V., Paulraj, A., and Kailath, T. (1987). Performance analysis of the optimum beamformer in the presence of correlated sources and its behavior under spatial smoothing. Acoust. Speech Signal Process. IEEE Trans. 35, 927–936. doi: 10.1109/TASSP.1987.1165239

Reite, M., Zimmerman, J. T., and Zimmerman, J. E. (1978). Human magnetic auditory evoked fields. Electroenceph. Clin. Neurophysiol. 45, 114–117. doi: 10.1016/0013-4694(78)90349-8

Robinson, S., and Vrba, J. (1999). “Functional neuroimaging by synthetic aperture magnetometry,” in Recent Advances in Biomagnetism, ed T. Yoshimoto (Sendai: Tohoku University Press), 302–305.

Salti, M., Monto, S., Charles, L., King, J. R., Parkkonen, L., and Dehaene, S. (2015). Distinct cortical codes and temporal dynamics for conscious and unconscious percepts. Elife 4:e05652. doi: 10.7554/eLife.05652

Spencer, M. E., Leahy, R. M., Mosher, J. C., and Lewis, P. S. (1992). “Adaptive filters for monitoring localized brain activity from surface potential time series,” in Signals, Systems and Computers, 1992 Conference Record of The Twenty-Sixth Asilomar Conference on (Pacific Grove, CA), 156–161.

Stufflebeam, S. M., Tanaka, N., and Ahlfors, S. P. (2009). Clinical applications of magnetoencephalography. Hum. Brain Mapp. 30, 1813–1823. doi: 10.1002/hbm.20792

Vaina, L. M., Calabro, F., Lin, F.-H., and Hämäläinen, M. S. (2010). “Long-range coupling of prefrontal cortex and visual (mt) or polysensory (stp) cortical areas in motion perception,” in 17th International Conference on Biomagnetism Advances in Biomagnetism – Biomag2010, eds S. Supek and A. Sušac (Berlin; Heidelberg: Springer), 298–301.

Van Veen, B. D., and Buckley, K. M. (1988). Beamforming: a versatile approach to spatial filtering. ASSP Mag. IEEE 5, 4–24. doi: 10.1109/53.665

Van Veen, B. D., Van Drongelen, W., Yuchtman, M., and Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. Biomed. Eng. IEEE Trans. 44, 867–880. doi: 10.1109/10.623056

Wang, J.-Z., Williamson, S., and Kaufman, L. (1993). Magnetic source imaging based on the minimum-norm least-squares inverse. Brain Topogr. 5, 365–371. doi: 10.1007/BF01128692

Wendel, K., Väisänen, O., Malmivuo, J., Gencer, N. G., Vanrumste, B., Durka, P., et al. (2009). EEG/MEG source imaging: methods, challenges, and open issues. Comput. Intell. Neurosci. 2009, 1–12. doi: 10.1155/2009/656092

Wu, S. C., Swindlehurst, A. L., Wang, P. T., and Nenadic, Z. (2012). Efficient dipole parameter estimation in EEG systems with near-ML performance. IEEE Trans. Biomed. Eng. 59, 1339–1348. doi: 10.1109/TBME.2012.2187336