Temporal Autocorrelation-Based Beamforming With MEG Neuroimaging Data

Abstract

Characterizing the brain source activity using magnetoencephalography (MEG) requires solving an ill-posed inverse problem. Most source reconstruction procedures are performed in terms of power comparison. However, in the presence of voxel-specific noises, the direct power analysis can be misleading due to the power distortion as suggested by our multiple trial MEG study on a face-perception experiment. To tackle the issue, we propose a temporal autocorrelation-based method for the above analysis. The new method improves the face-perception analysis and identifies several differences between neuronal responses to face and scrambled-face stimuli. By the simulated and real data analyses, we demonstrate that compared to the existing methods, the new proposal can be more robust to voxel-specific noises without compromising on its accuracy in source localization. We further establish the consistency for estimating the proposed index when the number of sensors and the number of time instants are sufficiently large. In particular, we show that the proposed procedure can make a better focus on true sources than its precedents in terms of peak segregation coefficient. Supplementary materials for this article are available online.

• Beamforming
• MEG neuroimaging
• Source localization and reconstruction
• Temporal autocorrelations

1. INTRODUCTION

When stimulus information is being processed and transmitted in the brain, electrical currents, reflecting neuronal response to a stimulus, will be produced and flow through the brain because of neuron firing, which in turn generates orthogonally oriented electromagnetic fields around them in accordance with Maxwell’s equations. The electrical currents are often approximated by neural dipoles. Each dipole is determined by its location, orientation, and time-varying amplitudes (called time-courses). The magnetic fields induced by the neuronal current flow can be measured at the scalp via noninvasive brain imaging devices such as magnetoencephalography (Hämäläinen et al. 1993). The magnitudes of these measured fields are directly related to neuronal current strength, allowing characterization of the amplitudes of electrical brain activity on a millisecond time scale. There are many uses for the magnetoencephalography (MEG), including assisting surgeons in localizing a pathology and assisting researchers in determining the function of various parts of the brain, among others. For example, Henson et al. (2011) conducted a MEG-based experiment, where a healthy young adult undergoes a series of perception judgments on faces or scrambled faces. They aimed at localizing brain regions for human face-perception. While MEG offers a direct measurement of neural activity with very high temporal resolution, its spatial resolution is low, requiring to solve an ill-posed problem of source localization in which we reconstruct source dipoles that underpin the observed magnetic fields. The problem is challenging since the observed magnetic field could result from an infinite number of possible neuronal dipoles. The purpose of this article is to address the challenge within the context of adaptive spatial filtering and to improve the analysis of the above face-perception experiment.

Various approaches can be used to tackle the above challenge. These approaches roughly fall into two categories, namely (spatially and temporarily) global methods and local (i.e., spatially local but temporally global) methods. The former includes Bayesian parametric methods while the latter takes sensor covariance-based spatial filters as special cases. The local methods can be or not be data-adaptive. In nonadaptive settings, a regularized L₂ estimator of the time-course at each time-instance is used. Examples of nonadaptive methods include minimum-norm estimator and its variations (Lin et al. 2006). They are not data-adaptive in the sense that they can be converted to sensor covariance-based spatial filters, where the sensor covariance matrix is estimated by a forward model alone without using the data. The data-adaptive filters are essentially projection-based filters, where the goal is to look for a set of weights that project the sensor data, subjected to the minimum variance constraints (van Veen et al. 1997; Robinson and Vrba 1998; Huang et al. 2004). By the construction, the minimum variance beamformers are scalable to big MEG data via parallel computing, in particular when multiple trials and multiple subjects are involved. Other advantages of the beamformer methods over global methods include that there are no priori assumptions about the number of the underlying sources and that no parametric assumptions are required about background noise distributions (Hillebrand et al. 2005; Friston et al. 2006). Minimum variance beamformers (including linearly constrained minimum variance (LCMV) and synthetic aperture magnetometry (SAM) as special cases) and minimum-norm estimator (MNE) have been widely used in source localization and implemented in main neuroimaging software. See Sekihara and Nagarajan (2010) for a review.

Although a number of studies have shown that the minimum variance beamformers are very useful in practice (e.g., Brookes et al. 2011), the following practically important issues remain to be addressed. First, in the face-perception MEG experiment, the sensor measurements are often broken into distinct blocks (termed epochs), forming multiple trial data for each stimulus. This gives rise to two different approaches to estimating the sensor covariance: (i) using the covariance of the average sensor measurements, where we average sensor measurements over trials first and (ii) using the average trial covariance, where we calculate individual trial covariances first and then average them. The applicability of these two approaches depends on types of neuronal activities to be detected (Friston et al. 2006). Neuronal responses to stimuli can be categorized as evoked or induced. Evoked and induced responses differ in their phase-relationships to the stimulus. Evoked responses are phase locked to the stimulus, whereas induced responses are not. In the presence of induced responses, the first approach suffers from loss of signals due to averaging sensor measurements over trials as shown in Section 3. In this case, one prefers to calculate individual trial covariances first and then average them. Unfortunately, the resulting sensor covariance estimate is often contaminated by background noises, giving a spatially blurred map of neuronal activity as shown by the first two rows of the plots in Figure 1, where we applied the SAM procedure to a face-perception dataset. The high noise level has obscured the distinction between sources and noises on the map. Second, most beamforming analysis is performed in terms of power (i.e., marginal variance of the projected MEG data) comparison. However, the power index can miss temporal-correlation patterns of source time-courses as pointed out by Gross et al. (2001). In the presence of strong voxel-specific noises, the power index can also give a noise-contaminated map at scalp. Under certain regularity conditions, Zhang, Liu, and Green (2014) showed that the SAM can make an error in discriminating the true source from the noisy voxel when the signal-to-noise ratio (SNR) is too small to be estimated accurately.

Figure 1 Neuronal activity maps with the face-perception data in session 5: The plots of the log-contrasts between the faces and scrambled faces on three orthogonal slices through the global peak locations and on 20 transverse slices under the scheme of sensor covariance estimation indexed by ma. The rows 1 and 2, 3 and 4, 5 and 6, 7 and 8, and 9 and 10 are for the SAM, LCMV, LCMV4, depth-weighted MNE, and TAB, respectively. In each subplot, the scale of the log-contrasts between the faces and scrambled faces are normalized to their maximum value and the log-contrasts have been thresholded by zero below. The highlighted yellow colored areas show neuronal activity increases for the faces relative to the scrambled faces. All the subplots are overlaid on the anatomical MRI scan of the brain.

To tackle the above issues, we need to move away our focus from a simple power analysis to a nonpower analysis. Here, we present a family of temporal autocorrelation-based beamformers (short for TABs) as a complement of the SAM approach, and develop the related theory. We propose the use of Ljung–Box statistic as a tool for summarizing the temporal information at each voxel. The Ljung–Box statistic was originally developed for measuring the deviation of a time series from white noises (Ljung and Box 1978). The basic premise behind the new proposal is that any interesting source signal arising from stimulus-related neuronal activities is at least nonwhite. Therefore, a good neuronal activity index should be able to capture temporal deviation from white noises, while bypassing the interference of background noises. The proposed procedure aims at identifying neuronal responses with high autocorrelations. We demonstrate that the proposed procedure can perform favorably in terms of the focus on peaks in the face-perception MEG data analysis. In particular, we identify several new patterns for neuronal responses to face and scrambled-face stimuli. Simulation studies suggest that for multiple trial data, when there exist induced responses, the average-first covariance estimation scheme can substantially deteriorate the accuracy of beamforming. The studies also show that the proposed procedure can work very well even for moderate numbers of sensors and time-instances when the MEG data are noisy. The proposed approach gives better results, in terms of localization bias, compared to standard SAM, LCMV, high-order sensor covariance-based LCMV, and depth-weighted MNE. We carry out a theoretical comparison of the proposed index with the SAM index in terms of estimation consistency and segregation coefficient in the mapping, where the segregation coefficient of a neuronal activity index at voxel r is defined as the index ratio at r and at a peak. The smaller the segregation coefficient, the better resolution (i.e., separation between sources and nonsources) is provided. Our theoretical analysis presents a clear picture of how the spatial and temporal dimensions affect the performance of the proposed procedure. We show that the estimation of the new index is consistent as the numbers of sensors and time-instances, n and J are sufficiently large (say, n ⩾ 90, J ⩾ 500, and $\sqrt{log\left(n\right)/J}\le 10\%$ as demonstrated in the simulations). The large J provides a better estimate of covariance matrix of n sensors, while the large n entails a more accurate source mapping based on the estimated sensor covariance. The proposed index is immunized against noises as its value is mainly determined by the SNR at each voxel. Compared to the SAM map, the new activity map can have a better segregation coefficient in some regions of interest.

The article is organized as follows. The details of the proposed beamformer estimation are provided in Section 2. The simulation studies and real data applications are conducted in Section 3. A theoretical analysis of the proposed inference procedure is provided in the Appendix. The conclusions are made in Section 4. The proofs of the theoretical results and more results on the simulations and the data analysis can be found in the online supplementary materials.

2 METHODOLOGY

Let Y_i(t_j) be the measurement recorded by the MEG sensor i at time t_j, and  Y(t_j) = (Y₁(t_j), …, Y_n(t_j))^T be the measurements from all n sensors at time t_j, where the time points t_j = j/Δ, 1 ⩽ j ⩽ J, the number of the time instants J = bΔ is determined by the time window b and the sampling rate Δ per second, and the number of the sensors n is of order hundreds. Sarvas (1987) showed that the contribution of an individual source to  Y(t_j) can be numerically calculated by the use of a Maxwell’s equation-based forward model and that the contributions of multiple sources can be summed up linearly. Accordingly,  Y(t_j) can be written as (2.1) $$\begin{array}{c}\hfill \mathbf{\text{Y}}\left(t_j\right)={\int }_{\Omega }\mathbf{\text{x}}\left(r\right)\beta (r,t_j)dr+ϵ\left(t_j\right),\end{array}$$(2.1) where Ω is the source space (i.e., the space inside the brain), β(r, t_j) is the source magnitude at location r with unknown orientation η(r), and  x(r) = l(r)η(r) is a linear function of the orientation $\eta \left(r\right)\in {\mathbb{R}}^3$ with l(r) being an n × 3 matrix (called lead field matrix) at location r. The columns l^x(r), l^y(r), and l^z(r) in l(r) are the noiseless output of n sensors when a unit magnitude source at location r is directed in the directions of the x, y, and z axes, respectively. The lead field matrix is known in the sense that it can be calculated by solving a set of Maxwell’s equations (Sarvas 1987).

To simplify model (2.1(2.1) $$\begin{array}{c}\hfill \mathbf{\text{Y}}\left(t_j\right)={\int }_{\Omega }\mathbf{\text{x}}\left(r\right)\beta (r,t_j)dr+ϵ\left(t_j\right),\end{array}$$(2.1) ), we discretize the continuous source space by a grid {r₁, …, r_p}, which is distributed throughout the brain. We assume that the true sources are approximately located on the sieve when the sieve is sufficiently dense (i.e., p is sufficiently large). Let β(t_j) = (β₁(t_j), …, β_p(t_j))^T = (β(r₁, t_j), …, β(r_p, t_j))^T be the magnitude vector of the candidate sources at {r₁, …, r_p} and {β(r_k, t_j): 1 ⩽ j ⩽ J} the source time-course at r_k, where the superscript T indicates the matrix transpose. Letting  X = ( x₁, …,  x_p), the model (2.1(2.1) $$\begin{array}{c}\hfill \mathbf{\text{Y}}\left(t_j\right)={\int }_{\Omega }\mathbf{\text{x}}\left(r\right)\beta (r,t_j)dr+ϵ\left(t_j\right),\end{array}$$(2.1) ) can be discretized as follows: (2.2) $$\begin{array}{ccc}\hfill \mathbf{\text{Y}}\left(t_j\right)& =& \sum _{k=1}^p{\mathbf{\text{x}}}_k\beta (r_k,t_j)+ϵ\left(t_j\right)=\mathbf{\text{X}}\beta \left(t_j\right)+ϵ\left(t_j\right),\hfill \\ & & 1\le j\le J,\hfill \end{array}$$(2.2) where 1 ⩽ p < ∞, ϵ(t_j) is the background noise vector of the n sensors at time t_j. The strength of the signal in the time-course β_k( · ) at voxel r_k is often measured by its variability ${lim}_{J\to \infty }{\sum }_{j=1}^J{({\beta }_k\left(t_j\right)-{\bar{\beta }}_k)}^2/J$ (called source power), where ${\bar{\beta }}_k={\sum }_{j=1}^J{\beta }_k\left(t_j\right)/J$ is a temporal average of time-course β_k( · ). In practice, the number of sensors is of order hundreds while the number of candidate sources p, around a few thousands or more, is much larger than n. When p is much larger than n, the model estimation becomes challenging as there are a diverging number of candidate models that can fit to the data.

To localize unknown sources, we propose a family of beamformers by covariance thresholding and shrinkage, which are computed in the following two steps: We first work out a sensor autocovariance estimate by thresholding and shrinkage, and then decompose it and calculate the neuronal activity index. The proposed beamformers can be extended to the setting with two stimuli: we apply the procedure to the MEG data for each stimulus, followed by calculating a log-contrast (the logarithm of the ratio of the indices under two stimuli, an F-type statistic) at each voxel. This will create a log-contrast map. The global peak on the map indicates where the maximum of the log-index increase attains for one stimulus relative to the other. The details are given below.

2.1 Estimation of Sensor Autocovariances

The sensor autocovariance matrix is often estimated by the sample version, $$\begin{array}{ccc}\hfill \widehat{C}\left(l\right)& =& \left({\widehat{c}}_{ij}\left(l\right)\right)=\frac{1}{J}\sum _{j=1}^{J-l}\mathbf{\text{Y}}\left(t_j\right)\mathbf{\text{Y}}{\left(t_{j+l}\right)}^T-\overline{\mathbf{\text{Y}}}{\overline{\mathbf{\text{Y}}}}^T,\hfill \\ & & l=0,1,2,...,\hfill \end{array}$$where $\overline{\mathbf{\text{Y}}}$ is the sample mean { Y(t_j)}. For simplicity, we denote $\widehat{C}=\widehat{C}\left(0\right).$ Let the time-course β_k(t_j) be decomposed to two uncorrelated parts, uncontaminated time-course α_k(t_j) and noise e_k(t_j). Then, under certain conditions, we can show that $\widehat{C}\left(l\right)=C\left(l\right)+o_p\left(1\right)$ with (2.3) $$\begin{array}{ccc}\hfill C\left(l\right)& =& \sum _{k=1}^p{\gamma }_k\left(l\right){\mathbf{\text{x}}}_k{\mathbf{\text{x}}}_k^T+{\sigma }_0^2\left(l\right)I_n,0\le l\le J_0,\hfill \end{array}$$(2.3) where J₀ is the number of lags considered with a default value of 20, I_n is the n × n unit matrix, and $$\begin{array}{ccc}\hfill {\gamma }_k\left(l\right)& =& \underset{J\to \infty }{lim}\sum _{j=1}^{J-l}({\alpha }_k\left(t_j\right)-{\bar{\alpha }}_k)({\alpha }_k\left(t_{j+l}\right)-{\bar{\alpha }}_k)/J+{\sigma }_e^2\left(l\right),\hfill \\ \hfill {\sigma }_e^2\left(l\right)& =& \underset{J\to \infty }{lim}\sum _{j=1}^{J-l}(e_k\left(t_j\right)-{\bar{e}}_k)(e_k\left(t_{j+l}\right)-{\bar{e}}_k)/J,\hfill \\ \hfill {\sigma }_0^2\left(l\right)& =& \left\{\begin{array}{cc}{\sigma }_0^2,\hfill & l=0,\hfill \\ 0,\hfill & l\ne 0\hfill \end{array}\right.\hfill \end{array}$$are the contaminated autocovariance of the kth time-source at lag l, the autocovariance of the source noises at lag l, and the autocovariance of the sensor noises at lag l, respectively.

It has been known that the sample covariance is not a good estimator of the population covariance if its dimension n is large or if the sample covariance is degenerate (Bickel and Levina 2008). In MEG neuroimaging, the sensor sample covariance matrices can be nearly singular due to collinearity between nearby voxels. Here, we apply the thresholding procedure of Bickel and Levina (2008) to the sensor sample covariance $\widehat{C}$, obtaining $\widehat{C}(0,{\tau }_{nJ})=\left({\widehat{c}}_{ij}\left({\tau }_{nJ}\right)\right),$ where ${\widehat{c}}_{ij}\left({\tau }_{nJ}\right)={\widehat{c}}_{ij}I\left(\right|{\widehat{c}}_{ij}|\ge {\tau }_{nJ})$ and the threshold level τ_nJ is changing in n and J. Similarly, we define the thresholded autocovariance $\widehat{C}(l,{\tau }_{nJ})$ for l ⩾ 1. The key assumption behind covariance thresholding is that many of off-diagonal entries are zeros and the number of nonzero off-diagonal entries is growing slowly as the matrix dimension is increasing. It can be seen from the proof of Lemma 2 in the online supplemental materials that the above choice of the thresholding level gives the optimal convergence rate of $\widehat{C}(0,{\tau }_{nJ})$. In the field of sensor array processing, people often adopt a shrinkage approach, that is, by artificially add noise to the data to improve the mapping (Sekihara and Nagarajan 2010). In light of this, like in the software FieldTrip, in our implementation, if the smallest eigenvalue of the sample sensor covariance is too small, then the thresholded covariance may not be positive definite. This can be corrected by adding a small amount of noise (determined by the smallest eigenvalue of the noise covariance matrix or by the Backus-Gibert regularization) to the data (i.e., adding ε × I_p to the sample sensor covariance) after thresholding. See Huang et al. (2004) for the details.

Although we mainly focus on the thresholding approach, for a comparison, we also consider the optimal shrinkage covariance estimator of Ledoit and Wolf (2004), which is defined as $$\begin{array}{c}\hfill {\widehat{C}}_{\text{opt}}=\frac{b_n^2}{d_n^2}{\mu }_n{I}_n+\frac{d_n^2-b_n^2}{d_n^2}\widehat{C},\end{array}$$where $$\begin{array}{ccc}\hfill {\mu }_n& =& ⟨\widehat{C},I_n⟩,d_n^2=⟨\widehat{C}-{\mu }_n{I}_n,\widehat{C}-{\mu }_n{I}_n⟩,\hfill \\ \hfill {\bar{b}}_n^2& =& \frac{1}{J^2}\sum _{j=1}^J⟨{\mathbf{\text{Y}}}_j{\mathbf{\text{Y}}}_j^T-\widehat{C},{\mathbf{\text{Y}}}_j{\mathbf{\text{Y}}}_j^T-\widehat{C}⟩,b_n^2=min({\bar{b}}_n^2,d_n^2),\hfill \end{array}$$and the trace operator < A, B > =tr(AB^T)/n for any n × n matrices A and B. The idea behind the above shrinkage estimator is to find the optimal weighted average of the sample covariance matrix $\widehat{C}$ and the identity matrix via minimizing the expected squared loss. Under certain conditions ${\widehat{C}}_{\text{opt}}$ converges to the true covariance C as n tends to infinity, implying that ${\widehat{C}}_{\text{opt}}$ can be degenerate if C is (Ledoit and Wolf 2004). The autocovariance of the time-course at each grid is estimated by the autocovariance of the projected data along the weighting vector derived from the ${\widehat{C}}_{\text{opt}}$-based SAM.

2.2 Neuronal Activity Indices

We now describe the proposed beamformers in the following two steps: We first approximate the brain space by a regular three-dimensional grid and then work out a neuronal activity index at each grid point. The index distribution is overlaid on a structural image of the subject’s brain, creating a brain source map.

2.2.1 Projection Step

For any grid point r and orientation η, let l(r) be the n × 3 lead field matrix and  x = l(r)η the composite lead field vector. Given η, we derive the weighting vector $w\left(r\right)\widehat{=}w(r,\eta )$ by minimizing the projected power $w^T\widehat{C}(0,{\tau }_{nJ})w$, subject to w^T x = 1. This gives $w(r,\eta )=\widehat{C}{(0,{\tau }_{nJ})}^{-1}\mathbf{\text{x}}/{\mathbf{\text{x}}}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}\mathbf{\text{x}}$ and the estimated source time-course {w(r, η)^T Y(t_j): 1 ⩽ j ⩽ J} at r with the power $1/{\mathbf{\text{x}}}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}\mathbf{\text{x}}$ and SAM index $$\begin{array}{ccc}\hfill \text{SAM}(r,\eta )& =& \frac{w{(r,\eta )}^T\widehat{C}(0,{\tau }_{nJ})w(r,\eta )}{w{(r,\eta )}^{T}w(r,\eta )},\hfill \end{array}$$the estimated autocovariance $w{(r,\eta )}^T\widehat{C}(l,{\tau }_{nJ})w(r,\eta )$, and the estimated autocorrelation at lag l $$\begin{array}{c}\hfill \frac{w{(r,\eta )}^T\widehat{C}(l,{\tau }_{nJ})w(r,\eta )}{w{(r,\eta )}^T\widehat{C}(0,{\tau }_{nJ})w(r,\eta )}=\frac{{\mathbf{\text{x}}}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}\widehat{C}(l,{\tau }_{nJ})\widehat{C}{(0,{\tau }_{nJ})}^{-1}\mathbf{\text{x}}}{{\mathbf{\text{x}}}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}\mathbf{\text{x}}}.\end{array}$$Finally, the optimal orientation η(r) can be estimated via maximizing the SAM index, which is equivalent to calculating the eigenvector associated with the maximum eigenvalue of the matrix $l{\left(r\right)}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}l\left(r\right)$ relative to $l{\left(r\right)}^T\widehat{C}{(0,{\tau }_{nJ})}^{-2}l\left(r\right).$ The above maximum eigenvalue is taken as the value of the conventional SAM index ${\text{SAM}}_{c_0}\left(r\right)$ at r.

2.2.2 Mapping Step

The value of the TAB index at grid point r can be calculated by using the formula $$\begin{array}{c}\hfill {\text{TAB}}_{c_0}\left(r\right)=(J+2)J\sum _{l=1}^{J_0}\frac{\widehat{\rho }{(r,\eta \left(r\right),l)}^2}{J-l},\end{array}$$where $\widehat{\rho }(r,\eta \left(r\right),l)$ is the autocorrelation at lag l, that is, $$\begin{array}{c}\hfill \widehat{\rho }(r,\eta \left(r\right),l)=\frac{\eta {\left(r\right)}^{T}l{\left(r\right)}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}\widehat{C}(l,{\tau }_{nJ})\widehat{C}{(0,{\tau }_{nJ})}^{-1}l\left(r\right)\eta \left(r\right)}{\eta {\left(r\right)}^{T}l{\left(r\right)}^T\widehat{C}{(0,{\tau }_{nJ})}^{-1}l\left(r\right)\eta \left(r\right)}.\end{array}$$When r is running over the grid, ${\text{TAB}}_{c_0}\left(r\right)$ creates a map in the brain. The maximum peak of the map gives a location estimator of one of latent sources and the corresponding latent time-course at the peak can be estimated by projecting the sensor data along the optimal beamforming weight w(r, η). We also calculate the local peaks on the transverse slices of the brain, identifying multiple sources.

2.3 Choosing Tuning Constants

Note that the tuning constant J₀ in the TAB index depends on the temporal structure of the underlying time-course at each voxel. If we impose certain parametric models on time-courses, then the Akaike information criterion can be used to determine the optimal J₀. However, the computation would be very expensive when we run the order selection for thousands of voxels. Here, we simply choose a large J₀ so that TAB can capture temporal patterns in time-courses. For evoked and induced response, we show by simulations that the performance of TAB is not sensitive to the choice of J₀ when J₀ ⩾ 20. See Figure S14 in the online supplemental materials for the details.

To choose the thresholding level, we first estimate the background noise level by using the data. In practice, the MEG imaging is often run on a subject first without stimulus and then with stimulus. This allows us to calculate the sample covariance $\widehat{C}$ for the MEG data under stimulus as well as the sample covariance ${\widehat{C}}_0$ for the prestimulus data. The latter can provide an estimator of the background noise level σ²₀. In light of the convergence rate of the sensor covariance matrix estimator shown in the Appendix, we set the thresholding level ${\tau }_{nJ}=c_0{\widehat{\sigma }}_0^2\sqrt{log\left(n\right)/J}$ with c₀ being optimized over a range of values. We threshold $\widehat{C}$ by τ_nJ, where ${\widehat{\sigma }}_0^2$ is the minimum diagonal element in ${\widehat{C}}_0$. There are other estimators for σ²₀ such as the maximum and median diagonal elements. However, any scaling effects between these estimators will be compensated for by the rescaling of c₀ as t_nJ depends on c₀ and ${\widehat{\sigma }}_0^2$ only through their product. When c₀ = 0, the proposed covariance estimator reduces to the standard sample covariance implemented in the software FieldTrip. To optimize c₀, for each value of c₀, we apply the proposed procedure to the data and obtain the maximum SAM and TAB indices: (2.4) $$\begin{array}{ccc}\hfill {\text{SAM}}_{c_0}& =& max\{{\text{SAM}}_{c_0}\left(r\right):r\text{in}\text{the}\text{grid}\},\hfill \\ \hfill {\text{TAB}}_{c_0}& =& max\{{\text{TAB}}_{c_0}\left(r\right):r\text{in}\text{the}\text{grid}\}\hfill \end{array}$$(2.4) In both simulations and a real data analysis, we will show that c₀ ∈ D₀ = {0, 0.5, 1, 1.5, 2} has covered its useful range. We choose c₀ in which the corresponding index attains maximum or minimum, which are called ma and mi, respectively. By choosing c₀, ma intends to increase the maximum index value, while mi tries to reduce source interference. In the two stimuli case, we choose c₀ in which the corresponding log-contrast attains maximum or minimum. The similar schemes of ma and mi can be defined for other beamformers.

3 NUMERICAL RESULTS

We begin with the following notations. For any estimator $\widehat{r}$ of an source location r, let the localization bias be defined as $E|\widehat{r}-r|$, where $|\widehat{r}-r|$ is the L₁ distance between $\widehat{r}$ and r. Let the maximum coherence between two locations r₁ and r₂, ρ_max be defined by $$\begin{array}{c}\hfill {\rho }_{max}(r_1,r_2)=\underset{\left|\right|{\eta }_1\left|\right|=1,\left|\right|{\eta }_2\left|\right|=1}{max}\frac{{\left(l\left(r_1\right){\eta }_1\right)}^{T}l\left(r_2\right){\eta }_2}{\left|\right|l\left(r_1\right){\eta }_1\left|\right||·|l\left(r_2\right){\eta }_2\left|\right|},\end{array}$$which is a measure of collinearity between r₁ and r₂. In this section, we will examine the performance of the proposed beamformers under various scenarios by simulations and a real data analysis. We will consider both the thresholding-based and the optimal shrinkage-based sensor covariance estimators in calculating the neuronal activity indices for the SAM, TAB, LCMV, LCMV4 (the LCMV based on the sensor covariance matrix of order 3), and the depth-weighted MNE with weight ||l(r)||^{− 1} at voxel r, where ||l(r)|| is the Forbenius norm of the leads field matrix at r. See Huang et al. (2004) and Lin et al. (2006) for the definitions of the LCMV, LCMV4, and depth-weighted MNE procedures. In simulation studies, we will compare the TAB to the procedures SAM, LCMV, LCMV4, and depth-weighted MNE in terms of localization biases of the peak locations to the underlying sources and demonstrate that the TAB is more robust to noises than the other four procedures. We will show that the accuracy of the localization is affected by source location coherence as well as temporal correlations between sources.

3.1 Simulated Data

In the following, we will focus on finite sample scenarios with moderate values of $n^2\sqrt{log\left(n\right)/J}$, while a theoretical analysis will be presented for large $n^2\sqrt{log\left(n\right)/J}$ in the Appendix. By simulations, we will answer the following questions arising in the analysis of the face-perception experiment:

1. Does the TAB procedure outperform the existing procedures when heavy noises are presented in the data?

2. Has the TAB procedure been improved by using thresholded or shrinkage covariance estimator?

3. To what extent will the performance of the proposed procedure deteriorate by source location coherence and source correlations?

4. In the presence of induced responses, to what extent will the performance of the proposed procedure be affected by averaging multiple trial data?

For this purpose, in the simulations we considered two head models (one was produced by a real human subject while the other was based on a mathematical model) with different numbers of MEG sensors and a broad range of sources. We compared the TAB to the SAM, LCMV, LCMV4, and depth-weighted MNE. We presented only the results for J = 500 and J = 1000 due to space limitation and that the results are similar when J = 2000, 3000, and 4000. We considered four settings as follows.

Setting 1 (Single Trial Data With Evoked Responses): We created a 184-sensor MEG system (CTF/VSM) using a real adult human subject head shape, which was approximated by multiple local spheres and downloaded at http://fieldtrip.fcdonders.nl/download (see Oostenveld et al. 2011). We constructed 2905 regular three-dimensional grid points of resolution 1 cm within the head. These candidate source positions were aligned with the axes of the head coordinate system. A lead field matrix L (with dimension 184 × 8715) between the sensors and the grid points was then calculated by using the open software FieldTrip.

We assumed that there were two nonzero sources θ₁(t) and θ₂(t) located at r₁ = (5, −5, 6)^T and r₂ = (5, 5, 6) in the front lobe of the brain (which is known to be important for syntactic processing) with the dumping cosine patterns θ₁(t) = η₁β₁(t) and θ₂(t) = η₂β₂(t), respectively. The L₁ distance and the maximum coherence between the two sources are equal to 10 and 0.2181, respectively. Here, $$\begin{array}{ccc}\hfill {\eta }_1& =& {\left(\frac{2}{\sqrt{6}},\frac{1}{\sqrt{6}},\frac{1}{\sqrt{6}}\right)}^T,{\eta }_2=(1,0,0),\hfill \\ \hfill a_1& =& \sqrt{6},a_2=8/5,\hfill \\ \hfill {\beta }_k\left(t\right)/a_k& =& z_k\left(t\right)\hfill \\ & & +\left\{\begin{array}{cc}exp\left(-{\left({\displaystyle t-\frac{m_k}{m_0}-\pi }\right)}^2\right)\hfill & \\ \times \left({\displaystyle 1-\frac{{(t-m_k/m_0-m_0/60,000)}^2}{{(1-1/m_0+m_0/60,000)}^2}}\right)\hfill & \\ \times cos\left(2\pi f_k\left({\displaystyle t-\frac{m_k}{m_0}}\right)-\pi \right),\hfill & \\ {\displaystyle \frac{m_k}{m_0}\le t\le 1+\frac{(m_k-1)}{m_0}}\hfill & \\ 0,\hfill \end{array}\right.\hfill \end{array}$$otherwise, where m₀, m_k and f_k are two factors related to time-shifts and frequencies of the cosine waves, and z_k(t) ∼ AR(1), that is, z_k(t) = 0.2z_k(t − 1) + e(t), {e(t)} is a source-specific (white) noise process with mean 0 and variance 0.1². We generated a pair of time-courses for each of the following combinations of m₀, (m₁, m₂) and (f₁, f₂): (1) m₀ = 250, (m₁, m₂) = m₀ × (1/8, 1/4) and (f₁, f₂) = (1, 3); (2) m₀ = 250, (m₁, m₂) = m₀ × (1/8, 1/8) and (f₁, f₂) = (1, 1).

Setting 2 (Single Trial Data With Phase-Locked Oscillations): It is similar to Setting 1. The details can be found in the online supplemental materials.

Setting 3 (Multiple Trial Data With Randomly Phased Oscillations and 91 MEG Sensors): Motivated by Davies-Thompson and Andrews (2012), we assumed that there were two nonzero sources θ₁(t) and θ₂(t) located at the precuneus region r₁ = ( − 5, 5, 5)^T and r₂ = ( − 4, −4, 8) with the random-phased oscillatory patterns θ₁(t) = η₁β₁(t) and θ₂(t) = η₂β₂(t), respectively. We used the same head model as in Setting 2. The L₁ distance and the maximum coherence between r₁ and r₂ are 13 and 0.0489, respectively. The two locations have been shown to be in a face-perception region of the brain (Zhang, Liu, and Green 2014). Here, $$\begin{array}{ccc}\hfill {\eta }_1& =& {\left(\frac{10}{\sqrt{102}},\frac{1}{\sqrt{102}},\frac{1}{\sqrt{102}}\right)}^T,\hfill \\ \hfill {\beta }_1\left(t\right)& =& \sqrt{102}cos\left(20t\pi \right),0\le t\le 1.\hfill \\ \hfill {\eta }_2& =& {\left(1,0,0\right)}^T,{\beta }_2\left(t\right)=8cos\left(20t\pi \right),0\le t\le 1.\hfill \end{array}$$

Setting 4 (Three Sources With 184 MEG Sensors): It is similar to Setting 1. The details can be found in the online supplemental materials.

The pairs of signal curves β₁(t) and β₂(t) in these settings were plotted in Figure S1 in the online supplemental materials. From the left to right in each row, the plots were made for Setting 1(1), Setting 1(2), Setting 2, Setting 3 with a single pair of randomly phased oscillatory functions, and Setting 3 with the average of 40 randomly phased oscillatory functions, respectively. The top and bottom rows were made for the sample rates J = 500 and J = 1000, respectively. In particular, the last two plots in each row show that after averaging, the amplitudes of randomly phased oscillatory functions were reduced by 85%. Note that the temporal correlation coefficients between the paired curves in the top row are, respectively, 0.066, 0.912, 1, 1, and 1 in the top row and 0.028, 0.912, 1, 1, and 1 in the bottom row. The 10 scenarios can be categorized in two groups in terms of their temporal correlation coefficients, namely, weakly temporally correlated group and strongly temporally correlated group.

In each of Settings 1, 2, and 3, we simulated the sensor measurements at t from the model (3.1) $$\begin{array}{c}\hfill Y\left(t\right)=l\left(r_1\right){\eta }_1{\beta }_1\left(t\right)+l\left(r_2\right){\eta }_2{\beta }_2\left(t\right)+ϵ\left(t\right),\end{array}$$(3.1) where ϵ(t) is the sensor noise vector. In Setting 4, the model is (3.2) $$\begin{array}{c}\hfill Y\left(t\right)=l\left(r_1\right){\eta }_1{\beta }_1\left(t\right)+l\left(r_2\right){\eta }_2{\beta }_2\left(t\right)+l\left(r_3\right){\eta }_3{\beta }_3\left(t\right)+ϵ\left(t\right).\end{array}$$(3.2) The sensors were measured at the time instants t_k = bk/J, k = 0, 1, 2, …, J − 1. The signal strength (SS) in the sensor space was defined by $$\begin{array}{c}\hfill \text{SS}=\sqrt{{\sum }_v{\sum }_{k=0}^{J-1}\left|\right|l\left(r_v\right){\eta }_v{\beta }_v\left(t_k\right){\left|\right|}^2/J}.\end{array}$$For each k, we sampled N_nk from an n-dimensional standard Normal N_nk and set $ϵ\left(t_k\right)=\text{SS}\times N_{nk}/\sqrt{\text{SNR}}.$ We considered two values of SNR, 1/20² and 1/30² in Settings 1 and 4, 1/25 and 1/0.64 in Setting 2, and 1/10² and 1/15² in Setting 3.

In each scenario of Settings 1, 2, and 4, we generated 30 independent datasets of {Y(t_k), 0 ⩽ k ⩽ J − 1} paired with {ϵ(t_k), 0 ⩽ k ⩽ J − 1} by models (3.1(3.1) $$\begin{array}{c}\hfill Y\left(t\right)=l\left(r_1\right){\eta }_1{\beta }_1\left(t\right)+l\left(r_2\right){\eta }_2{\beta }_2\left(t\right)+ϵ\left(t\right),\end{array}$$(3.1) ) and (3.2(3.2) $$\begin{array}{c}\hfill Y\left(t\right)=l\left(r_1\right){\eta }_1{\beta }_1\left(t\right)+l\left(r_2\right){\eta }_2{\beta }_2\left(t\right)+l\left(r_3\right){\eta }_3{\beta }_3\left(t\right)+ϵ\left(t\right).\end{array}$$(3.2) ), respectively. In Setting 3, we generated 60 × 40 independent datasets of {Y(t_k), 0 ⩽ k ⩽ J − 1} paired with {ϵ(t_k), 0 ⩽ k ⩽ J − 1} by model (3.1(3.1) $$\begin{array}{c}\hfill Y\left(t\right)=l\left(r_1\right){\eta }_1{\beta }_1\left(t\right)+l\left(r_2\right){\eta }_2{\beta }_2\left(t\right)+ϵ\left(t\right),\end{array}$$(3.1) ), forming 60 multiple trial datasets, with 40 trials each. For the first 30 multiple trial datasets, we estimated the sensor covariances by the strategy of calculating individual trial covariances first, followed by an average (short for “first-covariance-then-average”). For the remaining 30 multiple trial datasets, we estimated the sensor covariances by averaging the sensor measurements over trials first and then calculating the covariance of the averaged data (short for “first-average-then-covariance”). Here, we imitated a practical scenario, where the MEG imaging was run on a subject first without stimulus and then with stimulus. The former provides an estimator of the background noise level. The width b of time-window for sampling is 2 in Setting 1 and 1 in Settings 2, 3, and 4.

For each dataset, we calculated the sample covariance $\widehat{C}$ of {Y(t_k), 0 ⩽ k ⩽ J − 1} derived from the stimulus data and the corresponding sample covariance ${\widehat{C}}_0$ of the background noises derived from the prestimulus data. Then, for each of TAB, SAM, LCMV, LCMV4, and depth-weighted MNE, we calculated their sensor covariance estimates for each stimulus, which were indexed by c₀ = 0, 0.5, 1, 1.5, 2,  sh, ma, and mi. Note that, when c₀ = 0, the resulting SAM, LCMV, and depth-weighted MNE reduce to their standard versions implemented in the FieldTrip. We applied the corresponding TAB, SAM, LCMV, LCMV4, and MNE procedures to each of the datasets, obtaining the TAB, SAM, LCMV, LCMV4, and MNE-based maximum location estimates, respectively. We then calculated the localization biases of these maximum location estimates to the nearest true sources over 30 datasets for these five procedures, respectively. The details of the results are displayed by the box-and-whisker plots in Figures 2 and 3, and in Figures S2– S13 in the online supplemental materials.

Figure 2 Setting 1 (Comparison of the TAB to the other procedures when J = 1000): Two sources are located at CTF coordinates (5, −5, 6)^T cm and (5, 5, 6)^T cm, respectively, with temporal correlation 0.066 in Setting 1(1) and 0.912 in Setting 1(2). The first four rows show the box-and-whisker plots of the global peak-derived minimum localization bias differences, say TAB-SAM, TAB-LCMV, TAB-LCMV4, and TAB-MNE against the sensor covariance estimates indexed by c₀ = 0, 0.5, 1, 1.5, 2,  sh, ma, and mi for the combinations of n = 184, SNR = 1/20², 1/30², and J = 1000, respectively, in Setting 1(1). The remaining four rows are for Setting 1(2). The red colored lines in the boxes are the medians. Note that the optimal shrinkage sensor covariance estimate is indexed by “sh.” The negative values of the differences indicate that the TAB can achieve less localization biases than its competitors.

Figure 3 Setting 3 (TAB performances under two schemes of handling multiple trial data): Two sources are located at CTF coordinates (3, −1, 4)^T cm and ( − 5, 2, 6)^T cm, respectively, with temporal correlation 1. Comparison of the performances of the TAB under two schemes of handling multiple trial data, one is averaging sensor measurements first and the other is calculating individual sensor covariances first and then averaging them. The left and right two columns are for the first scheme and the second scheme, respectively. The rows show, from the left to right, the box-and-whisker plots of the TAB index values and the TAB global peak-derived minimum localization biases to the underlying two sources against the sensor covariance estimates used indexed by c₀ = 0, 0.5, 1, 1.5, 2, sh, ma, and mi for the combinations of n = 91, the number of trials is 40, SNR = 1/10², 1/15², and J = 500, 1000, respectively. The red colored lines in the boxes are the medians. Comparing the left two columns to the right two columns, we can see that the localization biases will increase if we perform the TAB on average sensor measurements in the presence of randomly phased oscillatory signals. This is due to that averaging randomly phased sensor measurements may reduce the signal-to-noise ratios of the sources in the data.

The results can be summarized as follows:

• The results, which are displayed in the columns 2 and 4 of Figure 3 and also Figures S2, S5, and S11 in the online supplemental materials, demonstrate that using the sensor covariance estimators indexed by ma and  sh in the TAB can offer superior source localization than that using the conventional sample covariance estimator. Similarly, ma (or mi) and  sh can improve the performances of the conventional sample covariance estimator-based SAM, LCMV, LCMV4, and depth-weighted MNE. The details are omitted due to space limitations.

• The results, which are shown in the right two columns of the above figures, confirm that the TAB procedure can outperform the SAM, LCMV, LCMV4, and depth-weighted MNE in terms of localization bias when the data are noisy. For instance, the first rows in Figure 2 show that compared to the SAM, the TAB can reduce the localization bias by around 85% when the scheme ma is used and by around 75% when the scheme  sh is used.

• The plots in Figures S3, S6, S8, and S12 in the online supplemental materials confirm that the source location coherence does have severe effects on localizing the true sources. For example, r₂ was masked in Settings 1(1) and 1(2), whereas r₁ was masked in Settings 2 and 3. The effects of source correlations are also visible. For instance, comparing the rows 2 and 4 in Figure S2 in the online supplemental materials, we can find that the localization biases are almost doubled when SNR= 1/900 and the source correlation coefficient changed from 0.066 to 0.912.

• By comparing the localization bias plots in the columns 2 and 4 in Figure 3, we reveal that for randomly phased multiple trial data, using the first-covariance-then-average strategy for sensor covariance estimation in the TAB can offer a less biased localization than using the first-average-then-covariance strategy. For instance, the plots at row 1 and column 2 and at row 1 and column 4 in Figure 3 show that the localization biases reduced to zeros when we changed the strategy from using the first-average-then-covariance to using the first-covariance-then-average.

• Figure S10 in the online supplemental materials demonstrates that the TAB can perform better than the SAM, LCMV, LCMV4, and depth-weighted MNE even when adopting the first-average-then-covariance strategy in estimating sensor covariance matrix. This is not surprising because the TAB should be better than the SAM, LCMV, LCMV4, and depth-weighted MNE in detecting temporal deviation from background noises.

3.2 Face-Perception Data

We applied the proposed beamformers to six multiple trial MEG datasets acquired by Wakeman and Henson (Henson et al. 2011). The datasets were generated by six sessions with 146, 148, 150, 147, 148, and 148 face or scrambled face trials on a single subject, respectively. During each of the trials, a healthy young adult was presented with one exemplar of two types of visual stimuli: a face or a scrambled face. Each trial consisted of a central fixation cross (presented for a random duration of 400 to 600 ms), followed by a face or scrambled face (presented for a random duration of 800 to 1000 ms), followed by a central circle for 1700 ms. The subject used either his/her left or right index finger to report whether he/she thought the stimulus was symmetrical or asymmetrical vertically through its center. The MEG data were collected with a Neuromag VectorView system, containing a magnetometer and two orthogonal, planar gradiometers located at each of 102 positions within a hemispherical array situated in a light, magnetically shielded room. The sampling rate was 1100 Hz. Here, we analyzed the data from the 102 magnetometers only. We focused on identifying locations that showed activity increases for the faces relative to the scrambled faces.

To this end, we first created a grid system of 1 cm resolution with 1487 grid points, using the subject’s anatomical magnetic resonance imaging (MRI) data. Then, for each session we applied the neuroimaging software SPM8 to read and preprocess the recorded data, and to epoch the data for the face stimulus and the scrambled face stimulus, respectively. This gives rise to 146, 148, 150, 147, 148, and 148 epochs (i.e., trials) of 700 ms (770 time instants) with 200 ms prestimulus and 500 ms post-stimulus, which were corresponding to the six sessions, respectively. For each of the two stimuli, we calculated the averaged sample covariance $\widehat{C}$ over the trials by using the post-stimulus data and the corresponding averaged noise covariance ${\widehat{C}}_0$ by using prestimulus data. We thresholded $\widehat{C}$ by $c_0{\widehat{\sigma }}_0^2\sqrt{log\left(n\right)/J}$, where n = 102, J = 551, and ${\widehat{\sigma }}_0^2$ is the minimum diagonal element in ${\widehat{C}}_0$. We also calculated the optimal shrinkage estimator. We reported the corresponding results to the sensor covariance estimates indexed by the three schemes c₀ = 0, ma, and  sh, respectively.

After the preprocessing, for each session, we applied both the TAB, SAM, LCMV, LCMV4, and depth-weighted MNE to the face dataset and the scrambled face dataset, respectively, followed by calculating the corresponding log-contrast for each of 1487 grid points. We interpolated and aligned the log-contrasts with the structural MRI of the subject. We also averaged the log-index ratios at each voxel over six sessions, interpolated, and relayed them on the structural MRI of the subject.

Taking Session 5 as an example, to compare the results derived from these procedures, we plotted the log-contrasts for each procedure on 20 traverse slices of the brain and on the three orthogonal slices through the peak location (i.e., maximum location estimate) in Figure 1. The plots for other sessions can be found in the online supplementary materials. The plots in Figure 1 show a clear distinction between the SAM map (the first two rows) and the TAB map (the last two rows): The SAM gives rise to a blurred neuronal activity map due to contamination from noises, whereas the TAB produces a clean and more sensible neuronal activity map. The SAM attains the peak at the CTF coordinates ( − 4.05, −3.95, −0.05)cm, while the TAB attained the peak at the CTF coordinates ( − 1.05, 5.95, 4.05)cm. The TAB highlights much less activated areas (with light yellow colored) than does the SAM (with heavily yellow colored), suggesting a large part of the activated areas derived from the SAM may be artifacts generated by voxel-specified noises as their temporal autocorrelations are low. This agrees with the theoretical analysis in the next section. The SAM-derived global peak has failed in identifying a meaning location as the location ( − 4.05, −3.95, −0.05)cm is located in the cerebellum and not close to any known face-perception regions in the brain. In contrast, the TAB peak location ( − 1.05, 5.95, 4.05)cm is in the fusiform face area (FFA), which is well known to be a face-perception region. The LCMV4 is able to provide a focused map whereas the standard LCMV and depth-weighted MNE give slightly smeared maps. A similar phenomenon was observed in the median-nerve MEG response by Huang et al. (2004). However, unlike the TAB, the LCMV4 may achieve its focus at the cost of its accuracy in terms of localization bias as implied by our simulations. See Figure 2. The plots in Figures S17, S19, S21, S23, S25, S26 in the online supplemental materials show that the thresholding and shrinkage schemes perform similarly. So we focus only on the thresholding-based procedures below. To reveal the types of the estimated sources, in Figures S30– S35 in the online supplemental materials, we further estimated and plotted the time-courses at the first four local peaks by projecting the sensor measurements along the optimal beamforming vectors for the ma-based SAM and TAB, respectively. The TAB revealed that two evoked response patterns are different from those shown by the SAM. For example, for Session 5, at the location ( − 1.05, 6.95, 4.05)cm, the estimated time-course for the scrambled face is almost flat, whereas the one for the face has a dumping S-pattern in the time interval (100, 200). In particular, the TAB identified a monotone response patterns at location (7, 3, 3)cm. The time-courses derived from LCMV, LCMV4, and MNE are not presented here due to space limitations.

Similar conclusions can also be made for the other sessions by use of Table 1 and its extension in the online supplemental materials. Interestingly, Table 1 shows that the peak locations are varying across the sessions, due to a session-specific random effect. There are two possible explanations for the aforementioned phenomenon: i) there might be a dynamic network in the brain for face-perception. Only part of it would appear in each session and (ii) there might be a so-called habituation effect for the brain processing the stimuli, that is, in the first time when you saw a face, the brain responded in one way, but in the second, the third, or the fourth times, as you have seen the same face before, your brain lost interests and therefore responded differently. So there was a nonstationary response pattern over time, which increased the variance over sessions. However, the voxels unrelated to the network showed the same (low) level of variability over time.

Table 1 The peak locations and face-perception regions they belong to for the six runs and for the averaged log-index over the six runs. TAB_ma and SAM_ma denote the TAB and the SAM when the thresholding scheme ma is used. By N/A we mean that the location is not close to any face-selective regions

	SAMma		TABma
Sess.	Loc. (cm)	Reg.	Loc. (cm)	Reg.
1	(−0.05, 5.95,6.05)	PCu	(−5.45, 4.45,1.95)	Close to OFA
2	(−5.05,5.45,6.05)	PCu	(−5.05,0.95,4.95)	OFA
3	(7.05,2.45,1.55)	N/A	(7.95,−0.05,3.95)	N/A
4	(−2.95,−6.35,6.05)	N/A	(−6.95,2.05,6.95)	PCu
5	(−4.05,−3.95,−0.05)	N/A	(−1.05,5.95,4.05)	FFA
6	(6.05,5.45,2.55)	N/A	(4.45,3.05,0.55)	AMG
Average	(−6.95,2.05,7.05)	PCu	(−6.95,2.05,7.05)	PCu

Neuropsychological evidences supported that the face-perception process can involve multiple face-selective regions in human cortex such as the FFA, the occipital face area (OFA), the superior temporal sulcus (STS), the precuneus (PCu), and the amygdala (AMG) (e.g., Davies-Thompson and Andrews 2012). The process proceeds via a series of stages in which the OFA computes an early structural description of a face while higher-level face-selective regions such as the FFA and the anterior temporal lobe compute the invariant aspects of a face. From Table 1, we note that most of the peak locations derived from the TAB procedure are in one of these regions, whereas most of the peak locations derived from the SAM are not close to these regions. This indicates that the TAB can be a very useful tool to improve the SAM-based beamforming in practice. Compared to Henson et al. (2011), we revealed several more locations with different responses to the face and the scrambled face stimuli, reflecting the neuropsychological fact that there is a dynamic network in the brain for the face-perception (Luckhoo et al. 2012).

4 DISCUSSION AND CONCLUSION

We have proposed a family of beamformers for inverting electromagnetic models and detecting temporally autocorrelated sources in MEG. These beamformers are local methods in the spatial domain but not in the temporal domain. The proposed neuronal activity index is different from a source stability index proposed by Hymers et al. (2010), in which one needs to randomly split the source-courses into two datasets each containing half the epochs. One performs beamforming on each dataset and then checks the consistency of the results derived from the two datasets. Compared to the existing one, the new index is much easier for a theoretical analysis as it has skipped the random splitting step. The proposed beamformers are based on the thresholded sensor covariance estimator optimized over grid points $c_0{\widehat{\sigma }}_0^2\sqrt{log\left(n\right)/J}$, c₀ = 0, 0.5, 1, 1, 1.5, 2. These points are derived from the asymptotic theory. In Figures S15 and S16 in the online supplemental materials, we have demonstrated that arbitrary chosen but dense grids perform very similar to and sometimes worse than the scheme $\mathit{max}$. This is because the grid points in $\mathit{max}$ are chosen under the guidance of the theory.

In a standard beamforming analysis, active sources are often constrained to remain the same over the short time interval of interest (Gramfort et al. 2013). Here, to detect the evoked and induced responses, we have restricted our analysis in a short period. The implicit assumption is that the source is stationary in the analysis period. While this assumption is reasonable for the face-perception experiment, it is not a good model for nonstationary settings where multiple transient sources activate sequentially or simultaneously during a long time interval. In these settings, a varying index in time and frequency is required (Gramfort et al. 2013).

We have examined the performance of the proposed procedure by simulations and a real data analysis. We have also investigated the theoretical behavior of the proposed procedure in the Appendix. We have showed that the index estimation is consistent. Compared to its precedents, the proposed procedure has several advantages: It is more robust to voxel-specific noises, has a better segregation coefficient, can help reduce source localization bias, and is scalable to big MEG data. Our simulation studies have indicated that the proposed procedure can outperform the standard SAM, LCMV, high-order covariance-based procedure LCMV4, and depth-weighted MNE when the data are noisy. In particular, we have shown by using real MEG datasets that the proposed procedure can be a powerful tool for analyzing multiple trial MEG data. The main strength of TAB over traditional approaches lies in its robustness to background noises and its accuracy in detecting true sources. Unlike TAB, LCMV4 can be significantly biased to artifacts although it is robust to background noises.

There is a limitation in our theoretical analysis, where we assume the underlying sources are uncorrelated. Extending the current theory to a general setting along the lines discussed in Sekihara and Nagarajan (2010) is possible but very tedious. It is beyond the scope of the current article.

The proposed methodology can be used to generate an F-type test for locations in the brain that segregates two stimulus conditions. The conventional way to estimate the null distributions is to shuffle two stimulus conditions randomly to generate a few hundred permuted datasets. We then apply the proposed procedure to each of them, calculating the corresponding permuted log-contrast at each location in the brain. At each location, we count the percentage of times that the permuted log-contrasts are larger than or equal to the observed log-contrast, which can be taken as the estimated p-values of these F-type tests. Taking Session 5 in the above face-perception experiment as an example, we obtained 500 permuted datasets by performing the Matlab command “randperm” on the stimuli of 148 trials. We applied the TAB procedure to these permuted datasets, obtaining an estimated p-value at each grid point. For instance, the estimated p-values at the first five peaks on the TAB map are 0.062, 0.072, 0.024, 0.114, and 0.044, respectively. To highlight the significant (1 − p)-values, we thresholded them by 0.9 and overlaid them on the anatomical MRI scan of the brain in Figure S36 in the online supplemental materials. The maximum (1 − p)-value 0.973 (thus the minimum p-value 0.017) attained at the location with CTF ( − 0.95, 4.95, 1.55)cm in the FFA region. An open issue with the permutation-derived null distribution is that it may be substantially biased when the face and scrambled face data were not drawn from the same distribution. This hampers multiple testing adjusting.

The proposed procedure can be thought as performing nonorthogonal principal component decomposition on the sensor covariance matrix by exploiting the second-order structure of the sample covariance space. If the sample covariance space is highly nonlinear itself, then neither orthogonal principal component analysis (PCA) nor nonorthogonal PCA is useful. However, the concept of nonlinear component decompositions such as independent component analysis (Hyvärinen, Karhunen, and Oja 2001) can be modified to meet the challenge.

Finally, we note that the source localization problem has arisen from other scientific research areas, including biomedicine, radar, sonar, and geophysics. See Morris et al. (2011), Minsley, Sogade, and Morgan (2007), and the references therein. The proposed procedure can be applicable to these areas after slight modifications.

SUPPLEMENTARY MATERIALS

The detailed proofs of the theoretical results can be found in the online supplementary materials.

Additional information

Notes on contributors

Jian Zhang

This is an Open Access article. Non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly attributed, cited, and is not altered, transformed, or built upon in any way, is permitted. The moral rights of the named author(s) have been asserted.Jian Zhang is Professor, School of Mathematics, Statistics and Actuarial Science, University of Kent, Canterbury, Kent CT2 7NF, United Kingdom (E-mail: [email protected]). Li Su is Senior Research Associate, Department of Psychiatry, School of Clinical Medicine, and Department of Psychology, University of Cambridge, Cambridge, UK (E-mail: [email protected]). The authors thank Professor Richard Henson, MRC Cognition and Brain Sciences Unit, Cambridge, for sharing the MEG neuroimaging data with them and for his comments on the earlier draft of the article. The authors are also grateful to the editor, an associate editor, and two reviewers for their constructive comments that have helped to improve a number of important issues. The participation of Dr. Li Su in this project was partly supported by the European Research Council under project code 230570-NEUROLEX, by Medical Research Council Cognition and Brain Sciences Unit (U.1055.04.002.00001.01), and by the NIHR Biomedical Research Centre and Biomedical Research Unit in Dementia based at Cambridge University Hospitals NHS Foundation Trust.

Li Su

This is an Open Access article. Non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly attributed, cited, and is not altered, transformed, or built upon in any way, is permitted. The moral rights of the named author(s) have been asserted.Jian Zhang is Professor, School of Mathematics, Statistics and Actuarial Science, University of Kent, Canterbury, Kent CT2 7NF, United Kingdom (E-mail: [email protected]). Li Su is Senior Research Associate, Department of Psychiatry, School of Clinical Medicine, and Department of Psychology, University of Cambridge, Cambridge, UK (E-mail: [email protected]). The authors thank Professor Richard Henson, MRC Cognition and Brain Sciences Unit, Cambridge, for sharing the MEG neuroimaging data with them and for his comments on the earlier draft of the article. The authors are also grateful to the editor, an associate editor, and two reviewers for their constructive comments that have helped to improve a number of important issues. The participation of Dr. Li Su in this project was partly supported by the European Research Council under project code 230570-NEUROLEX, by Medical Research Council Cognition and Brain Sciences Unit (U.1055.04.002.00001.01), and by the NIHR Biomedical Research Centre and Biomedical Research Unit in Dementia based at Cambridge University Hospitals NHS Foundation Trust.

APPENDIX: THEORETICAL ANALYSIS

In this Appendix, we present a theoretical analysis for the proposed beamformers when both n and J are large enough. Here, allowing n to tend to infinity is just an analytic device for finding the spatial factors that affect the performances of the proposed beamformers. We assume that there are q( ⩽ p) unknown true sources among p candidate sources in model (2.2(2.2) $$\begin{array}{ccc}\hfill \mathbf{\text{Y}}\left(t_j\right)& =& \sum _{k=1}^p{\mathbf{\text{x}}}_k\beta (r_k,t_j)+ϵ\left(t_j\right)=\mathbf{\text{X}}\beta \left(t_j\right)+ϵ\left(t_j\right),\hfill \\ & & 1\le j\le J,\hfill \end{array}$$(2.2) ). We aim at screening for the true sources from the candidate sources by using the TAB. To make the model (2.2(2.2) $$\begin{array}{ccc}\hfill \mathbf{\text{Y}}\left(t_j\right)& =& \sum _{k=1}^p{\mathbf{\text{x}}}_k\beta (r_k,t_j)+ϵ\left(t_j\right)=\mathbf{\text{X}}\beta \left(t_j\right)+ϵ\left(t_j\right),\hfill \\ & & 1\le j\le J,\hfill \end{array}$$(2.2) ) identifiable, we assume the following condition:

(A1): The source processes β( · ) are temporally uncorrelated with the noise process ϵ(t). The trends α( · ) are temporally uncorrelated with the noise process  e( · ) and the limits of ${\sum }_{j=1}^J{({\alpha }_k\left(t_j\right)-{\bar{\alpha }}_k)}^2/J$ and ${\sum }_{j=1}^{J-l}(e_k\left(t_j\right)-{\bar{e}}_k)(e_k\left(t_{j+l}\right)-{\bar{e}}_k)/J$ exist in probability as J → ∞. The source time-courses at different voxels are asymptotically uncorrelated with each others. That is, for any fixed J₀, as J → ∞, $$\begin{array}{ccc}\hfill \sum _{j=1}^{J-l}(ϵ\left(t_j\right)-\overline{ϵ}){(\beta \left(t_{j+l}\right)-\bar{\beta })}^T/J& =& o_p\left(1\right),0\le l\le J_0,\hfill \\ \hfill \sum _{j=1}^{J-l}({\beta }_{k_1}\left(t_j\right)-{\bar{\beta }}_{k_1}){({\beta }_{k_2}\left(t_{j+l}\right)-{\bar{\beta }}_{k_2})}^T/J& =& o_p\left(1\right),k_1\ne k_2,\hfill \\ & & 0\le l\le J_0,\hfill \\ \hfill \sum _{j=1}^{J-l}(\mathbf{\text{e}}\left(t_j\right)-\overline{\mathbf{\text{e}}}){(\alpha \left(t_{j+l}\right)-\bar{\alpha })}^T/J& =& o_p\left(1\right),0\le l\le J_0,\hfill \end{array}$$and ${lim}_{J\to \infty }{\sum }_{j=1}^J{({\beta }_k\left(t_j\right)-{\bar{\beta }}_k)}^2/J$ exists in probability. Note that Sekihara and Nagarajan (2010) gave a heuristic discussion on the performance of the SAM when the second equation does not hold.

We proceed in the following two steps. First, we consider the ideal situation where the sensor autocovariance matrices are known. Then, we investigate the asymptotic performance of the beamformer when the sensor autocovariance matrices are estimated based on the sensor measurements at J time instants. To simplify the derivation, we reparameterize the model (2.2(2.2) $$\begin{array}{ccc}\hfill \mathbf{\text{Y}}\left(t_j\right)& =& \sum _{k=1}^p{\mathbf{\text{x}}}_k\beta (r_k,t_j)+ϵ\left(t_j\right)=\mathbf{\text{X}}\beta \left(t_j\right)+ϵ\left(t_j\right),\hfill \\ & & 1\le j\le J,\hfill \end{array}$$(2.2) ) as in Zhang, Liu, and Green (2014). For the notation simplicity, we still let  x_k and β_k( · ) stand for the lead field vector and the time-course after the reparameterization. Both the SAM index and the TAB index are invariant under this reparameterization. For any A⊆1: p = {1, 2, …, p}, define C_A(l) = ∑_{k ∈ A} x_k x^T_kγ_k(l) + σ²₀(l)I_n, $C_{[-j]}={\sum }_{k\ne j}{\mathbf{\text{x}}}_k{\mathbf{\text{x}}}_k^T{\gamma }_k\left(l\right)+{\sigma }_0^2\left(l\right)I_n$. Let C_k(l) denote C_{1: k}(l).

A.1 Beamforming With Known Sensor Autocovariances

When the sensor processes are observed over an infinite number of time instants, under the assumption (A1), the sensor autocovariance matrix C_p(l), 1 ⩽ l ⩽ h can be fully recovered. In this ideal situation, we perform beamforming directly on C_p(l), 1 ⩽ l ⩽ J₀ and then reconstruct the unknown sources by using Equation (2.3(2.3) $$\begin{array}{ccc}\hfill C\left(l\right)& =& \sum _{k=1}^p{\gamma }_k\left(l\right){\mathbf{\text{x}}}_k{\mathbf{\text{x}}}_k^T+{\sigma }_0^2\left(l\right)I_n,0\le l\le J_0,\hfill \end{array}$$(2.3) ). For this purpose, for any grid point r in the brain, let  x =  x(r, η) denote $\sqrt{n}l\left(r\right)\eta /\left|\right|l\left(r\right)\eta \left|\right|$, the composite lead field vector at r with orientation η. For any two locations r and r₁, define ρ( x(r, η),  x₁) =  x(r, η)^T x₁/n = 1 − || x(r, η) −  x₁||²/(2n), which measures how close (r, η) is to (r₁, η₁). In fact, we can show that $$\begin{array}{c}\hfill \left|\right|\mathbf{\text{x}}(r,\eta )-{\mathbf{\text{x}}}_1{\left|\right|}^2/\left(2n\right)\le 1-\rho {(\mathbf{\text{x}}(r,\eta ),{\mathbf{\text{x}}}_1)}^2\le \left|\right|\mathbf{\text{x}}(r,\eta )-{\mathbf{\text{x}}}_1{\left|\right|}^2/n,\end{array}$$where || x(r, η) −  x₁||²/n is the average of the squared coordinate Euclidean-distances between  x(r, η) and  x₁. Let ρ( x,  x₁,  y) = ρ( x,  y) − ρ( x,  x₁)ρ( x₁,  x₂) and LT(r) denote the theoretical TAB index at the voxel r $$\begin{array}{c}\hfill \text{LT}\left(r\right)=J(J+2)\sum _{l=1}^{J_0}\frac{1}{J-l}{\left(\frac{{\gamma }_r\left(l\right)}{{\gamma }_r}\right)}^2,\end{array}$$where γ_r(l) denotes the autocovariance function at r, γ_r = γ_r(0) and set γ_r(l)/γ_r = 0 when γ_r = 0. Note that γ_r(l), l ⩾ 1 are not affected by voxel-specific noises and that under condition (A1), power γ_r can be written as the uncontaminated power γ_r0 plus the voxel-specific noise level σ²_e. We have $$\begin{array}{c}\hfill \frac{{\gamma }_r\left(l\right)}{{\gamma }_r}=\frac{{\gamma }_r\left(l\right)/{\sigma }_e^2}{1+{\gamma }_{r0}/{\sigma }_e^2},\end{array}$$which implies that LT(r) depends on the voxel-specific SNR γ_r0/σ²_e and that LT(r) tends to zero when σ²_e is very large. This is in contrast to the SAM index, where the SAM index will tend to infinity when σ²_e tends to infinity. In the following, to make the sources identifiable, we assume that

(A2): For any $r\ne r_1,$ l(r) and l(r₁) are linearly independent in the sense that for any orientations η and η₁ and nonzero constant c₁ with $max\left\{\right|\left|l\left(r\right)\eta \right||,|\left|l\left(r_1\right){\eta }_1\right|\left|\right\}\ne 0$, we have $l\left(r\right)\eta \ne c_{1}l\left(r_1\right){\eta }_1$.

A.1.1 Single Source

We first investigate the single source case (q = 1), in which a single source of nonzero power is located at r₁ with orientation η₁. Note that in this setting LT(r) = 0 when $r\ne r_1.$ The following result shows that when the sensor covariance is known, the TAB-based source map can accurately recover the true source location and orientation, and the index estimation is consistent. In addition, the TAB-based map has a smaller segregation coefficient than does the SAM when n is sufficiently large.

Proposition A.1.

Under Conditions (A1) and (A2), TAB(r) attains the maximum uniquely at r₁ and tends to zero in a rate of (n(1 − ρ( x(r, η(r)),  x₁)²))^{− 2} for $r\ne r_1$. And for r = r₁, $$\begin{array}{c}\hfill \frac{\text{TAB}\left(r_1\right)}{\text{LT}\left(r_1\right)}=\frac{1}{{\left(1+{\sigma }_0^2/\left(n{\gamma }_1\right)\right)}^2}=1-\frac{2{\sigma }_0^2}{n{\gamma }_1}+O(1/n^2),\end{array}$$when n is large enough. Furthermore, the segregation coefficient at location r admits $$\begin{array}{c}\hfill \frac{\text{TAB}\left(r\right)}{\text{TAB}\left(r_1\right)}=\frac{\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^4}{{\left(1+{\sigma }_0^{-2}{\gamma }_{1}n(1-\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^2)\right)}^2}.\end{array}$$

Remark A.1.

When n(1 − ρ( x(r, η(r)),  x₁)²) is large, the segregation coefficient at r has the rate of O((n(1 − ρ( x(r, η(r)),  x₁)²))^{− 2}). Note that Zhang, Liu, and Green (2014) showed that the SAM index also attains the maximum at the true source location r₁ with the segregation coefficient at location r, $$\begin{array}{c}\hfill \frac{\text{SAM}\left(r\right)}{\text{SAM}\left(r_1\right)}=\frac{1+{\gamma }_1{\sigma }_0^{-2}n(1-\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^2)}{1+(2+{\gamma }_1{\sigma }_0^{-2}n){\gamma }_1{\sigma }_0^{-2}n(1-\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^2)}.\end{array}$$In particular, when n(1 − ρ( x(r, η(r)),  x₁)²) is large, the corresponding segregation coefficient has the rate of O(1/n). Therefore, the TAB-based map can outperform the SAM in segregating the true source r₁ from r when n(1 − ρ( x(r, η(r)),  x₁)²) is large enough. Furthermore, the SAM map is less smoother than the TAB map in terms of their full width half maximums (Barnes and Hillebrand 2003). In fact, the full width half maximums (FWHMs) around r₁ for the SAM and the TAB can be, respectively, calculated by solving the equations SAM(r)/SAM(r₁) = 1/2 and TAB(r)/TAB(r₁) = 1/2, which are equivalent to solving $$\begin{array}{ccc}\hfill n(1-\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^2)& =& \frac{{\sigma }_0^4}{n{\gamma }_1^2},\hfill \\ \hfill n(1-\rho {(\mathbf{\text{x}}(r,\eta \left(r\right)),{\mathbf{\text{x}}}_1)}^2)& =& \frac{\sqrt{2}-1}{{\gamma }_1/{\sigma }_0^2+\sqrt{2}/n},\hfill \end{array}$$respectively. This proves that the SAM has a much smaller FWHM than does the TAB if we note that n(1 − ρ( x(r, η(r)),  x₁)²) is the squared distance between  x(r, η(r)) and  x₁.

A.1.2 Multiple Sources

We assume that there exist q unknown sources located at r_k, 1 ⩽ k ⩽ q with orientations η_k, 1 ⩽ k ⩽ q, respectively. We investigate the case where the underlying source locations are asymptotically separated in the so-called lead field distance. Let  x and  y denote the composite lead field vectors at locations r_x and r_y with orientation η_x and η_y, respectively. Let ${\rho }_{j_1{j}_2}=\rho ({\mathbf{\text{x}}}_{j_1},{\mathbf{\text{x}}}_{j_2}).$ We define the notations of partial correlations of sources a_yx|k, b_yx|k, c_yx|k, 1 ⩽ k ⩽ q in the online supplemental materials and introduce a regularity condition for identifying the underlying sources below.

For any 1 ⩽ j ⩽ q, we rearrange the sequences ( x₁, γ₁), …, ( x_q, γ_q) into a new order, namely, ( x₁, γ₁), …, ( x_{j − 1}, γ_{j − 1}), ( x_{j + 1}, γ_{j + 1}), …, ( x_q, γ_q), ( x_j, γ_j) and recalculate the corresponding values of a_{qq|(q − 1)} and b_{qq|(q − 1)} for this new sequence. These two new values are denoted by a_{jj|[ − j]} and b_{jj|[ − j]}, respectively. When a_{(k + 1)(k + 1)|k}’s are positive, a_{jj|[ − j]} is also positive. We call a_{jj|[ − j]} the correlation factor of the jth source with the remaining sources. Let a_nq = nmin _{1 ⩽ k ⩽ q − 1}a_{(k + 1)(k + 1)|k}.

We say that the source locations (and orientations) are asymptotically separable if they satisfy the condition $$\begin{array}{c}\hfill \text{(A3):}{a}_{nq}\to \infty \text{and}{b}_{yx|q}=O\left(1\right)\text{for}x\text{and}y\text{in}\text{the}\text{grid,}\text{as}n\to \infty .\end{array}$$

It follows from the definition that b_yx|q = O(1) if max _{2 ⩽ k ⩽ q}a_{kx|(k − 1)}/a_{kk|(k − 1)} = O(1) and max _{2 ⩽ k ⩽ q}a_{yk|(k − 1)}/a_{kk|(k − 1)} = O(1) (i.e., given { x₁, …,  x_{k − 1}}, the partial regression coefficients of  x and  y with respect to  x_k are bounded). In particular, if max _{2 ⩽ k ⩽ q}a_{kk|(k − 1)} is bounded below from zero as n tends to infinity, then the above partial regression coefficients are bounded and thus b_yx|q = O(1). Let k_m = 0 if a_nq → ∞ and k_m = min {1 ⩽ k ⩽ p − 1: na_{(k + 1)(k + 1)|k} = O(1)} if a_nq = O(1). Under assumption (A3), we can show that c_jj|q, 1 ⩽ j ⩽ q are bounded by the definition.

Imposing the condition that a_nq → ∞ as n → ∞ is equivalent to say that for any λ > 0, the maximum partial coherence of  x_k, 1 ⩽ k ⩽ p is bounded above by 1 − λ/n for a large n. Note that if a_nq → ∞, then all a_{(k + 1)(k + 1)|k} are positive when n is large enough. Therefore,  x_k, 1 ⩽ k ⩽ q are linearly independent for a large n because the inverse of the matrix $\left({\rho }_{j_1{j}_2}\right)$ can be obtained by iteratively performing sweep operations (q − 1) times on the matrix. See Goodnight (1979).

The following theorem shows that if max _{1 ⩽ j ⩽ q}max _{2 ⩽ k ⩽ q}a_{kj|(k − 1)}/a_{kk|(k − 1)} = O(1) (i.e., the partial regression coefficients of { x_j} are bounded), then a_nq → ∞ is sufficient for the theoretical TAB index to be consistently estimated. In particular, the estimation biases are of order (a_nq)^{− 1}. We now describe our general mapping theorem as follows.

Theorem A.1.

Under assumption (A3), for 1 ⩽ j ⩽ q, the TAB index admits $$\begin{array}{ccc}\hfill \frac{\text{TAB}\left(r_j\right)}{\text{LT}\left(r_j\right)}& =& 1+\frac{2{\gamma }_j{c}_{jj|q}}{n}+O\left(a_{nq}^{-2}\right).\hfill \end{array}$$In particular, when q = 2, we have $$\begin{array}{c}\hfill \frac{\text{TAB}\left(r_j\right)}{\text{LT}\left(r_j\right)}=1-\frac{2{\sigma }_0^2}{{\gamma }_{j}n(1-\rho {({\mathbf{\text{x}}}_1,{\mathbf{\text{x}}}_2)}^2)}+O\left(\frac{{\sigma }_0^4}{{\left(n(1-\rho {({\mathbf{\text{x}}}_1,{\mathbf{\text{x}}}_2)}^2)\right)}^2}\right).\end{array}$$For any voxel r_x, if LT(r_x) = 0, then TAB(r_x) = 0. Otherwise, if na_xx|q → ∞ (i.e.,  x is asymptotically separable from the true sources), then $$\begin{array}{ccc}\hfill \text{TAB}\left(r_x\right)& =& J(J+2)\sum _{l=1}^{J_0}\frac{1}{J-l}{\left(\sum _{k=1}^q{\gamma }_k\left(l\right)\frac{b_{kx|q}^2}{n{a}_{xx|q}}\right)}^2\hfill \\ & & \times \left\{1+\frac{2}{n}\left(\frac{2{\sum }_{k=1}^q{\gamma }_k\left(l\right)b_{kx|q}{c}_{kx|q}}{{\sum }_{k=1}^q{\gamma }_k\left(l\right)b_{kx|q}^2}-\frac{b_{xx|q}}{a_{xx|q}}\right)+O\left(a_{nq}^{-2}\right)\right\}.\hfill \end{array}$$And the segregation coefficient at r_x with respect to {r_j: 1 ⩽ j ⩽ q}, $$\begin{array}{c}\hfill \underset{1\le j\le q}{max}\frac{\text{TAB}\left(r_x\right)}{\text{TAB}\left(r_j\right)}=O\left({\left(n{a}_{xx|q}\right)}^{-2}\right).\end{array}$$

Remark A.2.

The above theorem indicates that the proposed index estimation is consistent with the underlying index when n is large enough and can have a better segregation coefficient when min _jna²_xx|q/a_jj|q → ∞ as n → ∞. The last point follows from the fact that under assumption (A3), $$\begin{array}{ccc}\hfill \frac{\text{SAM}\left(r\right)}{\text{SAM}\left(r_j\right)}& =& \frac{{\sigma }_0^2}{n{\gamma }_j{a}_{jj\left|\right[-j]}}\left(1+\frac{1}{n}\left(b_{xx|q}-\frac{{\sigma }_0^2(1+2{\gamma }_j{b}_{jj\left|\right[-j]})}{{\gamma }_j{a}_{jj\left|\right[-j]}}\right)\right.\hfill \\ & & \left.+O\left(a_{nq}^{-2}\right)\right)\hfill \\ & =& O\left({\left(n{a}_{jj|q}\right)}^{-1}\right).\hfill \end{array}$$See Zhang, Liu, and Green (2014).

A.2 Beamforming With Estimated Sensor Autocovariances

To derive a convergence result for ${\widehat{C}}_q\left(l\right)$ as both n and J tend to infinity, we reformulate the model (2.2(2.2) $$\begin{array}{ccc}\hfill \mathbf{\text{Y}}\left(t_j\right)& =& \sum _{k=1}^p{\mathbf{\text{x}}}_k\beta (r_k,t_j)+ϵ\left(t_j\right)=\mathbf{\text{X}}\beta \left(t_j\right)+ϵ\left(t_j\right),\hfill \\ & & 1\le j\le J,\hfill \end{array}$$(2.2) ) as $$\begin{array}{ccc}\hfill {\mathbf{\text{Y}}}_i\left(t\right)& =& {\phi }_i\left(t\right)+z_i\left(t\right),{\phi }_i\left(t\right)=\sum _{k=1}^q{x}_{ik}{\alpha }_k\left(t\right),\hfill \\ \hfill {\mathbf{\text{z}}}_i\left(t\right)& =& \sum _{k=1}^q{x}_{ik}{e}_k\left(t\right)+{ϵ}_i\left(t\right).\hfill \end{array}$$By the definition, we have E[ z_i(t)] = 0. Let ${\mathcal{F}}_{-\infty }^0=\sigma (\mathbf{\text{z}}\left(t\right):0\le t\le 0)$ and ${\mathcal{F}}_k^{\infty }=\sigma (\mathbf{\text{z}}\left(t\right):t\ge k)$. Define the mixing coefficient $a\left(k\right)={sup}_{A\in {\mathcal{F}}_{-\infty }^0,B\in {\mathcal{F}}_k^{\infty }}|P\left(A\right)P\left(B\right)-P(A\cap B)|.$ We need to impose more regularity conditions on the above model as follows.

(A4): For 1 ⩽ i ⩽ n,  z_i(t) is stationary and ergodic. There exist positive constants κ₁ and τ₁ such that for any u > 0, $$\begin{array}{c}\hfill \underset{1\le i\le n}{max}P\left(\left|\right|{\mathbf{\text{z}}}_i\left(t\right)\left|\right|>u\right)\le exp(1-{\tau }_1{u}^{{\kappa }_1}).\end{array}$$

(A5): For k ⩾ 0, $a\left(k\right)\le exp(-{\tau }_2{k}^{{\tau }_2})$.

(A6): ${max}_{1\le i\le q}{sup}_{t\ge 0}\left|{\alpha }_i\left(t\right)\right|<\infty .$ ${lim}_{J\to \infty }\frac{1}{J}{\sum }_{j=1}^J\alpha \left(t_j\right)$ exists and $$\begin{array}{ccc}\hfill \frac{1}{J}\sum _{j=1}^{J-l}(\phi \left(t_j\right)-\bar{\phi }){(\phi \left(t_{j+l}\right)-\bar{\phi })}^T& =& \mathbf{\text{X}}\text{diag}({\gamma }_{{\alpha }_1}\left(l\right),...,{\gamma }_{{\alpha }_q}\left(l\right)){\mathbf{\text{X}}}^T\hfill \\ & & +O\left(n\sqrt{log\left(n\right)/J}\right).\hfill \end{array}$$

Conditions (A4) and (A5) are required to prove the uniform ergodic results for  z_i(t) and  z_i(t) z_j(t). Condition (A6) can be verified when we take φ_i(t) as a trajectory sampled from a stationary process similar to the processes  z_i(t), 1 ⩽ i ⩽ n.

In the following theorem, we let γ_r(l) and γ_r denote the autocovariance function and variance at voxel r, respectively. We provide a convergence result for the estimation of the TAB index as follows. The result can be straightforwardly extended to the index-ratio.

Theorem A.2.

Under conditions (A1)∼(A6), if $n^2\sqrt{log\left(n\right)/J}\to 0$, then for 1 ⩽ j ⩽ q, $$\begin{array}{c}\hfill \text{TAB}\left(r_j\right)/\text{LT}\left(r_j\right)=1+O_p\left(1/n+a_{nq}^{-2}+n^2\sqrt{log\left(n\right)/J}\right)\end{array}$$as n and J tend to infinity. For voxel r_x with na_xx|q → ∞, $$\begin{array}{c}\hfill \text{TAB}\left(r\right)/J=O_p\left(1/n+{\left(n{a}_{xx|q}\right)}^{-2}+n^2\sqrt{log\left(n\right)/J}\right).\end{array}$$