Application of multi-source minimum variance beamformers for reconstruction of correlated neural activity

NeuroImage 58 (2011) 481–496

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / y n i m g

Application of multi-source minimum variance beamformers for reconstruction of correlated neural activity Alexander Moiseev a,⁎, John M. Gaspar b, Jennifer A. Schneider b, Anthony T. Herdman b, c a b c

Down Syndrome Research Foundation, Burnaby, Canada BC V5B 4J8 Simon Fraser University, Burnaby, Canada BC V5A 1S6 University of British Columbia, Vancouver, Canada BC V6T 1Z4

a r t i c l e

i n f o

Article history: Received 7 December 2010 Revised 18 May 2011 Accepted 21 May 2011 Available online 16 June 2011 Keywords: Magnetoencephalography (MEG) Inverse solutions Source analysis Minimum variance beamformers Correlated sources

a b s t r a c t Linearly constrained minimum variance beamformers are highly effective for analysis of weakly correlated brain activity, but their performance degrades when correlations become significant. Multiple constrained minimum variance (MCMV) beamformers are insensitive to source correlations but require a priori information about the source locations. Besides the question whether unbiased estimates of source positions and orientations can be obtained remained unanswered. In this work, we derive MCMV-based source localizers that can be applied to both induced and evoked brain activity. They may be regarded as a generalization of scalar minimum-variance beamformers for the case of multiple correlated sources. We show that for arbitrary noise covariance these beamformers provide simultaneous unbiased estimates of multiple source positions and orientations and remain bounded at singular points. We also propose an iterative search algorithm that makes it possible to find sources approximately without a priori assumptions about their locations and orientations. Simulations and analyses of real MEG data demonstrate that presented approach is superior to traditional single-source beamformers in situations where correlations between the sources are significant. © 2011 Elsevier Inc. All rights reserved.

Introduction Linearly constrained minimum variance (LCMV) beamformers (Van Veen et al., 1997; Robinson and Vrba, 1999; Sekihara and Nagarajan, 2008; Greenblatt et al., 2005; Huang et al., 2004; Herdman and Cheyne, 2009) are widely used for bioelectromagnetic source reconstruction. In magnetoencephalography (MEG) these adaptive spatial filters reconstruct signals from a given brain location using the output of a magnetic sensor array positioned outside the head. Each sensor registers tiny magnetic fields produced by post-synaptic currents in populations of neurons that are synchronously activated (Williamson and Kaufman, 1981; Hamalainen et al., 1993; Herdman and Cheyne, 2009). Mathematically, the inverse problem of reconstructing the sources of neural activity based on observed magnetic fields outside the brain is ill posed and has no unique solution (Baillet et al., 2001; Greenblatt et al., 2005). To find an inverse solution in the traditional minimum variance beamformer approach, the following two crucial assumptions are made. First, it should be possible to represent the brain activity as a superposition of elementary sources with known fields (the lead fields or forward solutions). The elementary sources used most often are point current dipoles although other choices are also possible (Limpiti et al.,

⁎ Corresponding author at: Down Syndrome Research Foundation, 1409 Sperling Avenue, Burnaby, Canada BC V5B 4J8. Fax: + 1 604 431 9248. E-mail address: [email protected] (A. Moiseev). 1053-8119/$ – see front matter © 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2011.05.081

2006; Jerbi et al., 2002). Second, the time courses of these elementary sources should be uncorrelated or statistically orthogonal. In practice the second assumption may often be violated. Indeed, one should expect the brain regions involved in certain tasks to be connected and to exhibit correlated activity, as is the case, with bilateral sensory areas being simultaneously activated. For example, bilaterally symmetric activations are observed in lateral visual cortices to visual stimulation (Quraan and Cheyne, 2010). The problem becomes even more salient when highly synchronized bilateral activations are sustained for prolonged intervals, as seen with auditory steady state response (ASSR) measurements (Herdman et al., 2003, 2004; Popescu et al., 2008). We should also mention a frequently used approach where covariance matrix of the measured field is calculated using short “active” time intervals around the stimuli rather than over the whole recording. In this case, correlated brain responses time locked to the stimuli may dominate the covariance and therefore affect the beamformer performance. When correlations are significant they result in the partial or complete cancelation of correlated sources as well as the distortion of their time courses (Widrow et al., 1982; Van Veen et al., 1997; Sekihara et al., 2002; Hillebrand and Barnes, 2005). A number of approaches have been proposed to overcome this problem (Brookes et al., 2007; Dalal et al., 2006; Popescu et al., 2008; Quraan and Cheyne, 2010; Hui et al., 2010). For bilateral activations, the most popular method is to use only sensors that are close enough to one of the correlated sources, or to use just one hemisphere of the sensor array. Then the contribution to the covariance from another source is greatly reduced and does not affect the

482

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

beamformer performance (see for ex. Herdman et al., 2003; Popescu et al., 2008; Herdman and Cheyne, 2009). The disadvantage of this approach is that using only a subset of sensors decreases the signal to noise ratio (SNR) and the spatial resolution for non-coherent activations. Another way of dealing with source correlations is described by Kimura et al. (2007). The idea is to construct the sensor covariance matrix that would be measured if the underlying brain sources were uncorrelated. First, the brain activity covariance matrix is estimated by the minimum norm method, then sources are “decorrelated” by discarding the non-diagonal terms of that matrix. Finally, a new version of the physical sensors' covariance is generated which is used in the beamformer analysis. In the approach taken by Brookes et al. (2007), instead of trying to get rid of the source correlations, they are incorporated into the lead field. The new lead field is composed as a linear combination of the dipolar forward solutions of participating sources. This is an example of a multi-source spatial filter, where filter weights for a target location explicitly depend on locations and orientations of all other sources. Most traditional beamformers are single-source, because their weights only depend on the target location. The method suggested by Brookes et al. allowed accurate localization of a pair of highly correlated sources, and in principle (with much computational effort) could be applied to a larger number of sources. As pointed out by the authors, the most significant problem here is that the reconstructed time course is a linear combination of the time courses of the participating sources. Finding individual time courses may be difficult, especially with an increasing number of sources. In publications by both Dalal et al. (2006) and by Hui et al. (2010), additional constraints are imposed on beamformer weights so that the beamformer gain from certain brain regions are (approximately) zero. In fact, this “nulling beamformer” (a term suggested by Hui et al., 2010) is another example of the multi-source approach, because the weights not only depend on the target voxel but also on the parameters of the “null” region. The constraints reduce the number of degrees of freedom available for the spatial filter which is similar to reducing the number of sensors in the array. However, with a priori knowledge of the locations of likely interfering sources, the null-region is small and additional zerogain constraints practically do not affect the beamformer performance. However, with no a priori information this region is large and the decrease in SNR and spatial resolution may become significant. A comparison of this method versus using a subset of the array sensors may be found in Popescu et al. (2008). In recent work by Quraan and Cheyne (2010), the same approach was further studied using both simulations and real MEG data. It was shown that indeed adverse effects of source correlations are strongly reduced, but the method may introduce distortions in source locations and amplitudes. These occur when correlated sources are close to each other, or when the “null”region is located near an existing source. Such distortions should be expected when single source localizers (e.g. pseudo-Z) are constructed using weights obtained in a multi-source setup which is what happens when null-gain constraints are imposed. Rather than looking at source localization, Hui et al. (2010) address the application of the “nulling beamformer” technique for functional connectivity analysis and show that this method effectively eliminates the “cross-talk” (linear mixing) between the sources. In this way, chances to observe spurious connections in the brain are significantly reduced. Recently Diwakar et al. (2011) proposed another spatial filtering technique for dealing with a pair of highly correlated sources, named Dual-Core Beamformer (DCBF). Instead of linearly combining source lead fields (Brookes et al., 2007), the authors start with a composite dual core “vector” beamformer whose lead field matrix includes the forward solutions of both sources. At this point their beamformer does not differ from that of previously mentioned authors (Dalal et al., 2006; Popescu et al., 2008; Quraan and Cheyne, 2010; Hui et al., 2010) nor the beamformers derived in this paper, and belongs to a general family of multiple constrained minimum variance (MCMV) filters to be discussed in a moment. New in this work are both the localizer

function and the method of how source amplitudes and orientations are found. Namely, the latter are obtained as solutions of a generalized eigenvalue problem — a method developed by Sekihara et al., 2004 for a traditional vector beamformer. In DCBF case the same expressions are used, with a single source lead field matrix replaced by the “dual core” lead field matrix. Thus computationally expensive search over possible orientations and relative amplitudes of the two sources that is a part of Brookes et al.'s (2007) procedure are no longer needed. While no mathematical justification is given regarding when such generalization of Sekihara's solution would be valid, computer simulations show that the method is successful for relatively strong source correlations (N0.5) and moderate to high SNRs. Here, we present a new multi-source beamformer method suitable for reconstructing both correlated and uncorrelated sources, assuming that observed brain activity is generated by a limited number of sources, possibly correlated, with known lead fields, plus noise. We additionally postulate that the sources are current dipoles; however the approach itself does not rely on any specific choice of the source. If locations and orientations of all sources are known then a beamformer solution is available which does not suffer from source cancelation effects. This is called a multiple constrained minimum variance (MCMV) beamformer, and it has been known in the sonar and radar fields for decades (Frost, 1972). Mathematically the MCMV beamformer also belongs to the family of linearly constrained minimum variance (LCMV) filters. The acronym MCMV is used here because in MEG and EEG literature, the term LCMV beamformer now has a narrower meaning and is synonymous with the so called “vector” beamformers (Sekihara et al., 2004; Sekihara and Nagarajan, 2008). We follow this tradition here to avoid confusion. As previously mentioned, the MCMV filter serves as a starting point for the approach developed by Dalal et al., 2006; Hui et al., 2010 (see also Popescu et al., 2008 and Quraan and Cheyne, 2010). Source cancelation does not occur with the MCMV beamformer because filter weights for each source are orthogonal to the lead fields of all the other sources. Of course in practice the brain sources are not known a priori and the question is how they can be found. In single-source spatial filters this is done by looking for the maxima of the beamformer-reconstructed activity distribution (see section Source localization: single source beamformers below). A similar procedure with a multi-source filter is not directly applicable in practice because of the extremely large number of combinations of possible source locations and orientations. Furthermore, the choice of the activity measure itself is not obvious. Therefore, this paper serves to address two related issues. First, we derive several measures (or localizers) and prove that they yield unbiased estimates of the source parameters. Thus, all the sources may be accurately determined by finding a global maximum of a localizer function defined on a parameter space of high dimension. Still finding this global maximum in practice turns out to be a challenging task that requires development of special techniques, and this is a second issue we address. Some sort of a sequential procedure looks like a natural approach, especially because it has proved successful in other EEG/MEG beamforming methods like RAP MUSIC (Mosher and Leahy, 1998, 1999) and in the radar field (Stoica and Sharman, 1990). In this work, we develop an iterative algorithm (see section Source search procedure) which allows finding the sources with reasonable computational effort. Finally, we illustrate the results by applying proposed methods to both simulated and real MEG data. Methods MCMV beamformer formulation In this section we briefly overview the derivation of the classic MCMV beamformer (Frost, 1972; Sekihara and Nagarajan, 2008) for the MEG inverse problem. Let b(t) denote the M-dimensional column vector of the MEG sensor readings at time t, where M is the number of

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

sensors. Assume that b(t) is generated by n0 b M possibly correlated sources si(θ i, t), i = 1,…, n0 plus noise ν(t): n0

 i i i bðt Þ = ∑ si θ ; t h θ + νðt Þ

ð1Þ

i=1

In this expression, si(θ i, t) is instantaneous amplitude of i-th source, while vector θ i describes all other source parameters. For current dipoles, θ i = {r i, u i} where r i is the source position and u i is the orientation unit vector. The M-dimensional column vectors h i(θ i), ν(t) define the source lead fields and the noise measured by the array, respectively. We assume that source parameters θ i do not change with time, random processes si(θ i, t),ν(t) are stationary and mutually uncorrelated, and 〈ν(t)〉 = 0, where the angle brackets denote statistical averaging. In addition, we assume that 〈si(θ i, t)〉 = 0 (and therefore also 〈b(t)〉 = 0) unless explicitly stated otherwise. In the linear spatial filter approach, the goal is to find an estimate ^si for the unknown amplitude of source i as a weighted sum of the sensor array readings:
T M i i sˆi ðt Þ = ∑ wm bm ðt Þ≡ w bðt Þ

ð2Þ

m=1

Here w i = {w1,…, wM}T is the M-dimensional column vector of the beamformer weights for source i, and superscript “T ” denotes transposition. For the single-source case (n0 = 1) formulation (2) corresponds to the “scalar” beamformer, where the weight vector w depends on both the location and orientation of the source (Sekihara et al., 2004). For the MCMV filter, we additionally require vector wi to satisfy a unit-gain condition: wiT hi = 1; i = 1; …; n , and zero-gain conditions for all other sources: wiT h j = 0; i≠j; i; j = 1; …; n. Here n denotes the MCMV beamformer order and generally it differs from the true number of sources n0 that can be simply unknown. The unit-gain and zero-gain conditions may be written as one matrix equation: T

W H = In

ð3Þ

where M×n matrices H and W are defined as H(Θ)={h1,…, hn}, W(Θ)= {w1,…, wn}, Θ is a vector of parameters of all the sources: Θ={θ1,…, θn} and In is the n-dimensional identity matrix. Explicit expressions for the MCMV beamformer weights are found by minimizing the total reconstructed source power P,
 n 2 T P = ∑ bsˆi N = Tr W RW

ð4Þ

i=1

subject to conditions (3). The well known result (Frost, 1972; Van Veen et al., 1997; Sekihara and Nagarajan, 2008) is −1

W ðΘÞ = R


−1 T −1 H H R H

ð5Þ

where R=b bbT Nis the M×M sensor covariance matrix which includes both signals and noise. Then for the power (Eq. (4)) one gets
 −1 P = Tr S

ð6Þ

with S(Θ)=HTR− 1H. Note that matrix H depends on both locations and orientations of the sources; therefore, these expressions may be regarded as a generalization of the “scalar” beamformer to the multi-source case. For multi-source beamformers, the filter weights for any source i depend not only on its own location and orientation θi but also on the parameters of all other sources: wi =wi(Θ), i=1,…, n. The zero-gain conditions make MCMV beamformer insensitive to correlations between the sources if the zeros are placed exactly where needed. In this case it is clear from Eqs. (1), (2), and (3) that sources sj, j ≠ i do not contribute to reconstructed activity sˆi ðt Þ. Note also that

483

vector LCMV beamformer (Van Veen et al., 1997; Sekihara et al., 2001) is a special case of the MCMV filter. It reconstructs three Cartesian components of the impressed current at given location r. Therefore, this is just an MCMV filter with n = 3, θ i = {r, u i} and predefined orientations u 1 ⊥ u 2 ⊥ u 3. When the source is a current dipole with fixed orientation, its components si, i = 1,…, 3 are 100% correlated, yet the vector beamformer may be successfully applied in this situation and no source cancelation occurs. Source localization: single source beamformers In practice the number of sources n0, their locations and their orientations are not known, so a single-source beamformer (n=1) may be used to find the sources. For each voxel r of a discreet grid covering the region of interest (ROI) some measure of brain activity is calculated. Ideally, its maxima should indicate the source locations. The source orientation u is usually unknown unless anatomical constraints are applied but assume for the moment that u(r) is set. The simplest measure is of course the source power (Eq. (6)): P(r)=1/(h(r, u)TR− 1h(r, u)), but it turns out that maxima of the power are significantly shifted towards the center of the head due to a strong non-linear dependence of the forward solutions h on the depth, and do not correspond to the true source locations. In order to overcome this bias, the localizer should somehow correct for the depth dependence of the lead field (Greenblatt et al., 2005). Among possible power-based measures that satisfy this requirement (for example, see: Robinson and Vrba, 1999; Sekihara and Nagarajan, 2008; Greenblatt et al., 2005; Huang et al., 2004) we'll mention only those related to the multi-source localizers developed in this paper. The first one is Van Veen's neural activity index (AI) ζ (Van Veen et al., 1997):

 T −1 T −1 ζ = PAI ðθÞ = h N h = h R h

ð7Þ

where N stands for the covariance matrix of the noise registered by the senor array. For uncorrelated sensor noise N = b ννT N = σ 2IM and the activity index PAI equals (up to a constant multiplier) to the power ˜ = h = ‖h‖: calculated using normalized forward solutions (NFS) h

T  ˜ R−1 h ˜ = ‖h‖2 = hT R−1 h (Sekihara and Nagarajan, PNFS ðθÞ = 1 = h 2008; Quraan and Cheyne, 2010). The second localizer called “pseudo-Z” Z̶(θ) (Robinson and Vrba, 1999; Sekihara and Nagarajan, 2008) is defined as a ratio of reconstructed source power P to the noise projected by the filter:



 T T T −1 T −1 −1 Z ð̶ θÞ = w Rw = w Nw = h R h = h R NR h

ð8Þ

It should be added that there is a number of other measures that are not power-based but rather utilize spatial distributions of certain statistics of the reconstructed signals (so called Statistical Parametric Maps [SPMs]). As an example, dual state beamformers use distributions of the pseudo-T, pseudo-F statistics of power differences between two conditions (Robinson and Vrba, 1999; Vrba and Robinson, 2001; Barnes and Hillebrand, 2003). Spatial maps of excess kurtosis of the source signals are used to identify zones of epileptic activity in the brain (Kirsch et al., 2006). An SSI (Source Stability Index) estimator (Hymers et al., 2010) identifies source locations by the consistency of the response to repeating stimuli. In the current paper we'll touch the dual state beamformers in the Discussion, however analyses will be limited to the power-based localizers. When source orientation u(r) is not known, a search over orientations is required. Orientation umax that maximizes the selected activity measure should be used in expressions (7) and (8). An effective way to find umax by solving a generalized eigenvalue problem is presented in Sekihara and Nagarajan (2008) and Sekihara et al. (2004). A different type of localizer is the event-related beamformer (Robinson, 2004; Cheyne et al., 2007). In many situations event-

484

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

related (ER) brain responses to repeatedly presented stimuli are of interest. All other ongoing brain activity is regarded as noise. ERsignals that have stable latencies and consistently repeating time courses are called phase-locked or evoked responses. As opposed to oscillatory activity whose mean time course is zero, phase-locked responses may be effectively separated from the background by averaging MEG signals over trials. This idea applied to a beamformer led to ER-localizers, which typically have higher signal to noise ratios than beamformers based solely on power (i.e. PAI, PNFS or Z̶ ). Sekihara et al. (2005, and Greenblatt et al. (2005), showed that for an ideal case of a single dipolar source and diagonal white noise both PNFS(θ) and Z̶ðθÞ maximize for θ = θ0 where θ0 = {r0, u0} denotes the true source location and orientation. In other words, both these localizers are unbiased estimators of the source parameters. To our knowledge, the location bias of ER-localizers has not yet been studied. Multiple source beamformers To be able to find sources in a similar way with a MCMV beamformer (nN 1), at least two considerations need to be addressed. First, the activity measure is to be determined. Obviously total power P is not an option as it yields a location bias even in a single source case. In addition, for nN 1 matrix S=HTR− 1H in Eq. (6) becomes degenerate when θi →θ j for any pair of sources i, j resulting in singularity for power. This problem is dealt ˜ = H = ‖H‖ (Dalal et al., 2006) or by using weights with by using H   ˜ i = 1 = σ 2 wi = ‖wi ‖ (Quraan and Cheyne, 2010). However, authors w present no explanation why beamformer maxima found in this manner should coincide with true source locations if nN 1, even when the noise is diagonal. The second consideration addresses how to find the maxima of the activity distribution in the multi-source case. The straight “brute force” ROI scan requires approximately (V ×K)n weights calculations, where V is the total number of voxels and K is the number of tested orientations. For n≥2 this usually amounts to an ROI scan that makes the brute force approach infeasible with current research computing resources. In this section we address the first question; the second question is dealt with later (see “Source search procedure”). Table 1 lists four localizer functions that are applicable in a multi-source case for an arbitrary noise covariance matrix N. Matrix definitions are listed in Table 2. Each localizer function is a distribution in a 5n-dimensional Θ-space where n is the number of current dipoles. The first two localizers in Table 1 are power based and may be used for both evoked (phase-locked) and oscillatory activity. The other two apply when the evoked component of the source activity sðt Þ = bfsi gT N;i = 1; :::;n0 (and therefore the averaged field bðt Þ) does not average to zero. While all four localizers have important common properties, we will precede by examining each one in greater detail. Power-based multiple source localizers The first power-based estimator for n sources, is the “multi-source activity index” (MAI) defined by:

 1=2 −1 1=2 −1 PMAI ðΘÞ = Tr G S G −I n ≡Tr GS −n

Table 2 Common matrix definitions. Expression

Matrix size; comments

S G G0

HTR− 1H HTN− 1H G|Θ = Θ0

(n × n) (n × n) (n0 × n0); value of G for true source coordinates

T

HTR− 1NR− 1H D E

R

where n×n matrix G=HTN− 1H provides forward-solution normalization.

bb

t

H T R−1 RR−1 H bssT N −s sT; s≡bs N bssT Nt

E C C

(n × n) (n0 × n0); the true source covariance matrix (n0 × n0); 2nd moments of the averaged source amplitudes

For n = 1(Eq. (9)) reduces to PMAI = ζ − 1 where ζ is the “neural activity index” (Eq. (7)). ζ is a “(signal + noise)/noise”-type ratio, because in expression (7) (h TR − 1h) − 1 is the total reconstructed power that includes both source and noise contributions, while (h TN − 1h) − 1 is the power due solely to the noise. Therefore ζ − 1 is a “pure” signal-to-noise-type ratio, which equals zero when there are no sources. Form (9) for PMAI is a generalization of expression ζ − 1 to a multi-source case. Both matrices S(Θ)and G(Θ)become degenerate if any of the sources coincide (θ i → θ j), however Tr(GS − 1) remains finite (see Appendix A). Consequently function PMAI(Θ) has no singularities. It is also non-negative as the trace of a symmetric non-negative matrix. The second estimator in Table 1 is “multi-source pseudo-Z” (MPZ) — a generalization of expression (8) to the multiple source case. We set

 −1 = 2 −1 = 2 −1 PMPZ ðΘÞ = Tr T −n ST −I n ≡Tr ST

ð10Þ

where matrix T is specified in Table 2. This definition reduces to Z̶−1 for n = 1. Again, 1 is subtracted to obtain a pure signal-to-noise-type ratio that equals zero when there are no sources. Converting Eq. (8) to a matrix form leads to Eq. (10). Expression (10) is not equivalent to the ðI Þ “type I” SNR definition for the vector beamformer: Z̶V = TrðS Þ = TrðT Þ given in Sekihara and Nagarajan, 2008, sec 6.4 or variations of the latter found in Huang et al., 2004, nor to the localizer used by Quraan and Cheyne, 2010. Similarly to PMAI(Eq. (9)), PMPZ(Θ) has no singularities and is non-negative because it is the trace of a symmetric non-negative matrix. Event-related multiple source localizers To construct a localizer based on the phase-locked activity, allow the source time courses si(θ i, t) in Eq. (1) to have non-zero means. Then si(θ i, t) may be represented as a sum of the event-average time course si and a fluctuating part s˜i , where si is the phase-locked   response: si = si + s˜i ; si = hsi i; s˜i ≡0. It follows that b(t) is no longer a zero-mean process, Dand E can be represented in the same way: bi = bi + b˜ i ; bi = hbi i; b˜ i ≡0. The expression for the measurement covariance matrix R should now be written as def

ð9Þ

(n × n) (M × M); 2nd moments matrix of the averaged field

T

R=

D
D TE 
T E = b˜ b˜ b−b b−b

In the single-source case, one can form a ratio of power of the phaselocked response sˆðθ; t Þ = wT ðθÞbðt Þ to projected noise: Z E̶ R ðθ; t Þ =

Table 1 MCMV-based source localizers. Localizer name

Base matrix expression (“ME”) 1/2 − 1 1/2

Localizer function (Tr(“ME”)) −1

Multi-source activity index (MAI) Multi-source pseudo-Z (MPZ)

G S G − In T − 1/2ST − 1/2 − In

PMAI = Tr(GS ) − n 1 PMPZ = Tr(ST− ) − n

Multi-source event related beamformer (MER)

T − 1/2ET − 1/2

PMER = Tr(ET− )

Reduced multi-source event related beamformer (rMER)

S

− 1/2

ES

− 1/2

1

PrMER = Tr(ES

−1

)

Global maximum value max PMAI = Tr(CG0) max PMPZ = Tr(CG0)
 max PMER = Tr CG0  

−1 max PrMER = Tr C G−1 0 + C

A. Moiseev et al. / NeuroImage 58 (2011) 481–496



 T wT b b w = wT Nw and use it as a localizer (Cheyne et al., 2007). For more statistically robust results, this phase-locked power may be additionally averaged across a time interval (Robinson, 2004) (we denote this operation as 〈〉t in Table 2). Minding that weights w do not depend on time and using expression (5) we get
D 

 E 
T T T T −1 −1 T −1 −1 Z ̶ER = w bb w = w Nw = h R RR h = h R NR h t

D

T

where R = bb

E t

. Thus in a single source case we have a ratio of two

numbers Z E̶ R = E = T with E = h R−1 RR−1 h and T = h TR − 1NR − 1h. To obtain a multi-source event-related beamformer (MER), we replace the ratio E/T with its matrix version, the same way it was done for power-based beamformers PMAI(Eq. (9)) and PMPZ(Eq. (10)), and get: T



 −1 = 2 −1 = 2 −1 PMER = Tr T ≡Tr ET ET

ð11Þ

where E and T are now matrices specified in Table 2. PMER represents a “pure” SNR because by definition noise has no contribution to the phaselocked response so there is no need to subtract Tr(In) in Eq. (11). A variation of this localizer may be derived based on the following considerations. In many cases evoked activity only occurs in short time intervals and makes small contribution to the overall field covariance R. In other words matrices R and N do not differ much. For this reason, consider a “reduced” version of the PMER localizer where N is replaced with R. Then in Eq. (11) matrix T is substituted with S = H TR − 1H yielding

 −1 = 2 −1 = 2 −1 ≡Tr ES PrMER = Tr S ES

ð12Þ

While we derived Eq. (12) assuming R ~ N, it turns out that important properties that are discussed in the next section hold for PrMER without this assumption. An attractive feature of the “reduced” localizer PrMER is that noise covariance N does not enter expression Eq. (12). This may be a significant advantage because getting a statistically reliable noise estimate may be not easy. The price paid for this advantage is lower spatial resolution compared to the “full” PMER beamformer, as discussed below. Common properties and differences All multi-source localizers P* * (where P* * refers to any of PMAI, PMPZ, PMER, PrMER) have common properties listed below. Proofs may be found in Appendix A. 1. P* * is a non-negative function of source parameters Θ = {θ 1,…, θ n}. P* *(Θ) remains finite at singular points where θ i → θ j, i ≠ j. 2. P* *(Θ) reaches a global maximum when a) the beamformer order becomes equal to the actual number of sources n0, and b) coordinates Θ match those of true sources Θ0: max

max P ðΘÞ = P ðΘ0 Þ≡P n;Θ

In other words, P* *(Θ) provides an unbiased estimate of the true source locations and orientations. 3. Saturation: Once the global maximum is reached, a further increase in the beamformer of the localizer. n order does not change o the
value  ˜ = θ1 ; …; θn0 ; θn0 + 1 ; … P Θ ˜ = P ðΘ0 Þ = P max . I.e. for any Θ 0  0 Property 2 is crucial because it guarantees that when the forward model and covariance matrices R, N are known exactly, the beamformer maximum corresponds to exact locations and orientations of all the sources. This statement holds irrespective to any specific form of the noise covariance N, for arbitrary large correlations between the sources and for any signal to noise ratio. As a side result,

485

property 2 provides mathematical justification for use of eventrelated localizers in a single-source case. The saturation property places a natural limit on the beamformer order and has clear physical meaning. The expressions in Table 1 are various forms of signal to noise ratios. If some of true sources have not been found, the “signal” power may be increased by adding another source. Consequently the global maximum will not be reached if the beamformer order is smaller than the actual number of sources. When all sources are found, the “signal” power cannot grow any more no matter how many additional fictitious sources are added, therefore the ratio saturates (see Figs. 2, 6B). This property gives a principal means to determine the actual number of sources. However this reasoning applies only if all of the true sources are found correctly. If sources are missed or found inaccurately, then the beamformer value will continue to grow with the beamformer order even if the order becomes larger than the actual number of sources n and the global maximum will not be reached (see Fig. 6A). While possessing common properties listed above, the localizers in Table 1 are not equivalent. Consider PMAI and PMPZ first. Although both achieve the same global maximum value (see Table 1), it holds that always PMPZ(Θ) ≤ PMAI(Θ) and equality is only achieved when all true sources are found (Appendix A). In other words, PMPZ has a sharper peak compared to PMAI. This difference is most significant for high SNR situations. On the other hand, for this very reason PMAI is less sensitive to the source search procedure errors than PMPZ. Similarly, it can be shown that “reduced” evoked response beamformer PrMER has lower spatial resolution compared to PMER and for high SNR situations the latter is preferable. For very low SNR cases the difference becomes negligible and the advantage of not having to specify (potentially inaccurate) noise model makes PrMER preferable. Source search procedure With estimators in Table 1, all sources may be localized provided a global maximum of P* *(Θ) can be found. For a single source case (n = 1) there is no problem because P* *(Θ) simply needs to be calculated once for each voxel in the region of interest. This requires approximately V×K beamformer computations (K is the number of source orientations at each voxel). In the multi-source case (n N 1), this number becomes (V × K) n and grows exponentially with n because all combinations of locations and orientations are calculated. Typically V ~ 104 and the number of calculations becomes too large even for n = 2 and K = 1 (two sources with constrained orientation). This problem can be addressed in several ways. For example, Diwakar et al., 2011 used a modified Powell search procedure. Other algorithms for finding a global maximum of a function in highly dimensional spaces (i.e. simulated annealing) may also be considered. This way, all sources are found “at once” when the global maximum is identified. However, as was already mentioned, sequential algorithms, where sources are found one by one in a lower-dimensional parameter space, proved to be productive in similar situations (Mosher and Leahy, 1998, 1999; Stoica and Sharman, 1990). The implementation of this idea is usually task specific and concrete mathematical details should be worked out on a case by case basis. In this section we develop a sequential search procedure for the multi-source beamformers that requires only V × n beamformer calculations for localization of n sources. In practice, even strongly correlated sources are unlikely to cancel each other completely (see Sekihara et al., 2002; Hillebrand and Barnes, 2005; Quraan and Cheyne, 2010 and Numerical examples). For full cancelation sources should have approximately equal and relatively high SNR and very high correlation (≥0.8). Usually one of those conditions is not met, therefore the algorithm starts with a single source beamformer. The largest activation detected is designated as a first source. Its amplitude may have been reduced by coherent interference, but at the moment one is only interested in the location. On step k of the algorithm we fix coordinates of previously found (k-1) sources and label

486

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

these as “references”. The H-matrix may be written as H = {HR, hk(θ k)}, where columns of a “reference” matrix HR are forward solutions of already found sources: HR = {h1,…, hk − 1}. We have hk(θk) = H k(rk)uk where columns of matrix H k = {hxk, hyk, hzk} are lead fields for sources at rk, oriented along axes x, y and z. To find source k, vector rk is run over the entire ROI (V voxels). In Appendix B we report that the dipole orientation u k(rk) maximizing the localizer value can be found as an eigenvector corresponding to the largest eigenvalue λ of the following generalized eigenvalue problem:

 k k k k D r u = λF  r u

ð13Þ

Here D* *, F* * are 3 × 3 matrices. A search over the source orientations is not required. Expressions for D* *, F* * are
−1
−1
−1
−1 D = AkR ARR BRR ARR ARk −AkR ARR BRk −BkR ARR ARk + Bkk ; ð14Þ
−1 kk kR RR Rk F  = A −A A A

Here matrices A, B should be substituted with S, G for the MAI localizer, with T, S for the MPZ, with T, E for the MER or with S, E for the rMER localizers, respectively. Matrices S, G, T, E with upper indices coincide with those in Table 2, with H replaced by either H Ror H k. For example, S kR = (H k) TR − 1H R, etc. The largest maximum of P* *(θ1,…, θk − 1 ;rk, uk(rk)) as a function of rk determines the coordinates and orientation of the source number k. Finally the weight matrix W is calculated using Eq. (5). Note that when a new source is added parameters θi, ib k of the sources found earlier do not change, but their weights wi are updated. As a result, source cancelation effects are removed and previously hidden correlated activity is revealed. Ideally, the iterations should be stopped when the localizer value reaches a plateau. In practice this only happens in high SNR cases. The saturation effect turns out to be very sensitive to small errors, especially in orientation, and may not be observed with the approximate algorithm described above. We found that a more reliable way is to monitor singlesource “power” pseudo-Z values z k̶ = wkT Rwk = wkT Nwk of each new pffiffiffiffiffiffi source, or equivalently “amplitude” pseudo-Z's z ̶k . For event-related beamformers, the corresponding quantity is z k̶ = wkT Rwk = wkT Nwk where N is an estimate of the covariance of the averaged noise field: N = N = ðnumber of trialsÞ. Iterations should be terminated when z k̶ or pffiffiffiffiffiffi z ̶k becomes small and starts fluctuating around a baseline level. Note that in practice this baseline level is usually not equal to one and must be determined experimentally. This is because the true noise covariance is not known and its sample estimates or simple analytical models are used to calculate z k̶ . Numerical examples Simulations All modeling were done for a 151-channel radial gradiometer MEG system (CTF/VSM) using a real adult human subject head shape which

was approximated by multiple local spheres. First, we describe in detail three concrete examples (simulations I–III), then describe statistical results obtained when source locations and orientations were selected randomly. For simulations I–III, 300 trials of 2 s duration were created with a sample rate of 200/s and a 50-Hz low pass filter, simulating several point current dipoles plus noise. Either white diagonal noise (uncorrelated) or human brain noise collected from an adult subject was used. In all three simulations sources were positioned on the same axial plane and had identical orientations. Their time courses consisted of one 0.3-s long 20Hz sinusoidal pulse with a Hanning envelope per trial, starting at t = 1 s after the start of each trial. When the brain noise was used, covariances R and N were estimated over [1.0, 1.3] s “active” and [0.7, 1.0] s “baseline” time intervals, respectively. For white noise, analytical (infinite integration time) expressions for R and N were used. Parameters of the simulations are detailed in Table 3. A whole head search for maxima of the beamformer distribution on a 2 mm step grid with the total number of voxels V ~ 105 was performed using the algorithm described above. In Simulation I we considered a case where two relatively strong (20 nA m) parallel correlated dipoles (c = 0.8) were positioned far apart (see Table 3). Fig. 1A presents cross-sections of these distributions through the source plane for four iterations of the algorithm. True source locations were marked with crosses and circles show the positions picked up by the beamformer. The distributions were normalized on their max values (with reference locations excluded). In this case, all four localizers successfully found both sources in the first two iterations. It can be seen that PMPZ has slightly sharper maxima than PMAI. Likewise, PMER is marginally sharper than PrMER. The sources were so far apart that this difference in spatial resolution did not affect the results. The 1st iteration corresponded to conventional single-source beamformers. Even in this high SNR case, the second largest maxima of the singlesource power localizers were shifted from the true position of the second source by 2-3 mm due to coherent interference from the source 1. On iteration 2, correct locations were found. Fig. 1B demonstrates that dipole amplitudes obtained by a single-source beamformer (the first iteration) were reduced by more than 40% compared to the true values while multi-source beamformers reconstruct them almost perfectly. pffiffiffiffiffiffi Fig. 1C presents “amplitude” single source pseudo-Z z k̶ (i.e. the square root of the “power” pseudo-Z z ̶k ) for locations found in iteration k. For “sources” 3 and 4 the pseudo-Z is close to 1 meaning that we are simply sampling the noise field, and that these sources are not real. In other words it indicates that true sources are recovered and iterations may be stopped. The saturation observed in the P*max plot (see Fig. 2) suggests * the same. Imagine now that the true noise covariance matrix N is not known. In this situation it is quite common to resort to a white noise approximation for N. Fig. 3A shows PMPZ and PMER distributions for the same simulated data, but where diagonal sensor noise of 3 fT/√Hz was used in the beamformer calculations. Even in such a high-SNR case, location errors for the power-based beamformer (MPZ) for sources 1 and 2 are 1 cm and 2.5 cm, respectively. Their reconstructed amplitudes

Table 3 Setup for simulations I–III. For all simulations sources belong to one horizontal plane and have the same orientation u = {1, 0, 0}. For simulation III, source 1 is uncorrelated with other three; correlation between sources 2, 3, 4 equals to 0.8. Simulation #

Noise

Source #

{x, y, z}, cm

Amplitude, nA*m

RMS, nA m

Correlation

I

Brain Brain

III

White

{.5, 4.5, 8.5} {.5, − 5.1, 8.5} {.5, 4.5, 8.5} {.5, − 6.0, 8.5} {.0, − 4.5, 8.0} {2.0, 2.5, 8.0} {− 1.0, 3.5, 8.0} {.0, 5.5, 8.0}

20 20 6 6 51.7 163.6 136.6 51.7

8.6 8.6 2.6 2.6 22.2 70.2 58.7 22.2

.8

II

1 2 1 2 1 2 3 4

1 .8 or 0

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

487

Fig. 1. Simulation I. (A) — Cross sections of beamformer distributions through the source plane z = 8.5 cm; X and Y axes units are cm. The distributions are normalized on their maximum values over the head volume, with reference locations excluded. True source locations are marked with crosses, circles show projections of the locations picked up by a beamformer on the source plane. (B) — Reconstructed amplitudes of sources found in first two iterations after subtracting projected noise: light gray — iteration 1, dark gray — iteration 2, and black — the true pffiffiffiffiffiffi amplitude. (C) — Single source pseudo-Z z k̶ for beamformer maxima found in four iterations of the algorithm.

Fig. 2. Simulation I: Maximum values of the localizers over the whole head volume P*max * found in iterations 1–4, normalized on the P*max * in iteration #4. All curves saturate after both true sources are found.

are ~1.7 times larger than the true amplitudes (Fig. 3B). The singlesource pseudo-Z plot (Fig. 3C) indicates that other sources do exist. Both these effects occurred because the diagonal noise model explains only a small part of the true noise field. In particular, noise sources that were not included in the model contribute to the reconstructed source amplitudes because of leakage. For event-related beamformer (MER) the location errors are smaller (6 mm and 3 mm) but still significant and source amplitudes are also overestimated. At the same time, as we saw the noise-independent rMER beamformer localizes both sources perfectly (Fig. 1). This example demonstrates that even in high-SNR cases supplying a correct noise model is essential. Without such a model application of rMER localizer is beneficial. In simulation II the sources are still far apart, but are much weaker and the correlation coefficient has been increased to 1. Again, spatial resolution does not play a major role here so both power-based and event-related localizers yield similar results. Power-based beamformers completely missed source 2 on the 1st step of the algorithm (the single-source beamformer), and could only localize source 2 with about 3 mm error on subsequent steps (Fig. 4A). Because of this error, the amplitude of the 2nd source was not fully recovered (Fig. 4B) and the pseudo-Z value (Fig. 4C) for source 2 is as low as for the false sources 3, 4. Generally source 2 should be discarded as unreliable, unless other considerations (for example, the source time course or its anatomical location) indicate otherwise. Event-related beamformers

488

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

Fig. 3. Simulation I. (A) — Cross sections of PMPZ and PMER distributions through the source plane z = 8.5 cm normalized on their maximum values over the head volume, with reference locations excluded. X and Y axes units are cm. Diagonal sensor noise with density 3 fT/rtHz was used in the calculations instead of real brain noise. True source locations are marked with crosses, circles show projections of the locations found by a beamformer. (B) — Reconstructed source amplitudes in first two iterations after subtracting projected pffiffiffiffiffiffi noise: light gray — iteration 1, dark gray — iteration 2, and black — the true amplitude. (C) — Single source pseudo-Z Z k̶ for beamformer maxima found in four iterations.

Fig. 4. Simulation II. (A) — Cross sections of beamformer distributions through the source plane z = 8.5 cm; normalized on their maximum values over the head volume, with reference locations excluded. True source locations are marked with crosses, circles show projections of the locations found by a beamformer. (B) — Reconstructed amplitudes of sources pffiffiffiffiffiffi found in first two iterations after subtracting projected noise: light gray — iteration 1, dark gray — iteration 2, and black — the true amplitude. (C) — Single source pseudo-Z Z ̶k for beamformer maxima found in four iterations.

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

489

Fig. 5. Simulation III. (A) — Cross sections of beamformer distributions through the source plane z = 8 cm; normalized on their global maximum values, with reference locations excluded. True source locations are marked with crosses, circles show projections of the locations found by a beamformer. (B) — Reconstructed amplitudes of sources found in first six iterations after pffiffiffiffiffiffi subtracting projected noise: light gray — iteration 1, dark gray — iteration 6, black — the true amplitude. (C) — Single source pseudo-Z z k̶ for beamformer maxima found in six iterations.

performed much better than power based MAI and MPZ. Although the single-source versions (i.e. iteration 1) misplaced source 2 by ~ 8 mm, multi-source beamformers successfully recovered true dipole positions and amplitudes (Figs. 4A and B). Simulation III is designed as a difficult case that demonstrates differences in performance of various localizers. Three parallel correlated sources (c = 0.8) are close to each other. The fourth source is uncorrelated with the other three. Diagonal white noise of 15 fT/√Hz density in the 0.5-50 Hz frequency band was assumed. We start with the two power-based localizers. On the first step MAI identified the three-source cluster as one big blob whose maximum does not correspond to any of the true source positions (Fig. 5A). Further steps eventually recovered all correlated sources with about 1 cm location error. Finally the uncorrelated source was also found with no location error. The pseudo-Z plot (Fig. 5C) indicated that only three sources out of the four should be taken as reliable and large errors in reconstructed amplitudes of correlated sources were observed (Fig. 5B). Because of its better spatial resolution, the MPZ localizer correctly picked up the uncorrelated source and detected another three sources almost exactly right on the first step, although with amplitudes ~ 2 times smaller than real ones due to the source cancelation (Figs. 5A and B). Subsequent iterations yielded correct positions of all sources and almost perfect reconstruction of their amplitudes. It is also clear from the pseudo-Z plot (Fig. 5C) that all sources were recovered and no further iterations were needed. Similarly, the MER localizer performed much better than rMER (Fig. 5A). Unlike MPZ, MER located correlated sources with small (several millimeters) errors. This was because signal to noise ratio for event-related beamformers is higher and therefore source cancelation

effects were more pronounced. For example, on the 1st step amplitudes of correlated sources found by MER were 5–10 times smaller than the real ones, as opposed to only 50% reduction in the case of MPZ (Fig. 5B). The rMER yielded much larger location errors, and only three of the four sources were recovered reliably according to the pseudo-Z plot. Thus with the good knowledge of the noise field, using the MER beamformer is preferable. Beamformer global maximum values found as a function of the iteration number (or equivalently, the beamformer order) are shown in Fig. 6A. No beamformer demonstrated clear saturation except for MPZ, and even that one continued ramping after all true sources were found. Interestingly enough the MER beamformer that performed quite well for all practical purposes was farthest from the theoretical saturation limit. Evidently, this was a consequence of higher SNR and sharper maxima in the Θ-space: the MER beamformer became more sensitive to localization errors that occurred when an approximate search procedure was applied. Fig. 6B illustrates that if we were able to find the exact locations and orientations, perfect saturation would indeed be achieved. To evaluate beamformer performance for a larger number of possible situations, we ran a set of simulations, where source locations and orientations were selected randomly. Among numerous possible choices of parameters, we decided to use configurations that should be most difficult for traditional beamformers — namely, where source correlations are significant and with SNRs equal for all sources. Again, multi-sphere head model and real human brain noise were used. In each run there were three sources, with correlation coefficients c = 0.8 for sources 1 and 2, c = 0.5 for sources 2 and 3 and c = 0.2 for sources 1 and 3. The SNR was defined as a square root of the ratio of the signal power, registered by the array, to total noise power: SNR =

490

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

Fig. 6. Simulation III: Maximum values of the localizers over the whole head volume P*max * as a function of iteration number, normalized on the theoretical (true) global maximum of corresponding beamformer. (A) — Results obtained using the search procedure; (B) — Results obtained if exact source locations and orientations were found on each step.

s(‖h‖ 2/tr(N)) 1/2. For a randomly selected position/orientation θ, h(θ) was calculated and source amplitude s was adjusted to obtain required SNR value. Realizations where the amplitude of any source exceeded a physiologically plausible range of 3 to 70 nA m were discarded. For every beamformer type, three sets of 200 simulations were run — one set for each of three possible SNR values: 0.5 (strong sources), 0.25 (moderate sources) and 0.1 (weak sources). To make it more intuitive, these SNR values roughly corresponded to the tangential source located 2 cm below the head surface, with amplitudes 28, 14 and 6 nA m, respectively. To locate the sources, a whole head search on the 5 mm-step grid was performed using analytically calculated covariance matrix R. Note that generated source positions did not coincide with the grid voxels, which is more likely the case in a real life scenario. Otherwise, over-optimistic estimates of the beamformer performance were obtained. (For example, average localization errors of 1 to 2 mm were obtained even for the weak sources.) The simulation results for multi- and single-source beamformers are presented in Tables 4 and 5. For each beamformer type, column ‘Success rate’ shows percent of successful localizations. Localization was considered successful only if each of the three sources were found with location error not exceeding 1 cm. Column ‘d, cm’ contains location error averaged over the sources and over all runs. Column ‘Rel. amp’ presents averaged ratio (in percents) of reconstructed source amplitude to the true source

amplitude. Depending on the SNR, different performance results for different beamformer types were observed. When sources were strong (SNR = 0.5) single and multi-source beamformers had comparable success rates in the range 78–88%, except for MPZ, MER beamformers. These failed more frequently simply because their maxima were missed when true sources fell in between the voxels. For comparison, when uncorrelated sources coinciding with the grid voxels were simulated for the same SNR value, the success rate was slightly above 90%, while in 10% of the cases at least one of the sources landed close to the center of the head or in peripheral areas with poor sensor coverage, and could not be located with 1 cm accuracy. Thus, 90% rate should be regarded as the maximum achievable success rate for our setup. Average location errors for single- and multi-source beamformers did not differ. The biggest difference was in reconstructed source amplitudes. Single source beamformers lost more than 40% of the source amplitudes (or ~60% in power) due to the source cancelation. Multi-source beamformers lost 8–12%, partly due to the discrete grid that prevented nulls from being placed exactly where needed. When sources were weak (SNR=0.1), the situation was different. While multi-source beamformers were mostly successful (72–77%), single source beamformers failed in more than half of the runs (success rates 26–46%). Among multi-source localizers, MAI beamformer was the worst. At the same time, the source cancelation effect was greatly reduced,

Table 4 Random simulations results for multi-source beamformers. SNR

0.5 0.25 0.1

MAI

MPZ

MER

rMER

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

83 80 72

0.3 0.34 0.37

92 96 95

72 86 76

0.28 0.31 0.36

88 94 96

72 86 76

0.28 0.31 0.36

88 94 96

88 86 77

0.27 0.3 0.35

90 95 95

Table 5 Random simulations results for single-source beamformers. AI, PZ, ER and rER denote the single-source versions of MAI, MPZ. MER, rMER, respectively. SNR

0.5 0.25 0.1

AI

PZ

ER

rER

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

Success rate

d, cm

Rel. amp

78 56 26

0.3 0.36 0.4

58 72 93

86 73 46

0.28 0.32 0.4

57 69 93

86 73 46

0.28 0.32 0.4

57 69 93

86 73 46

0.28 0.32 0.4

57 69 93

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

and both single source beamformers (when they did not fail) and multisource ones reconstructed 93–95% of the true source amplitude. For moderate sources (SNR=0.25) the situation was between those two extremes, with 80–86% success rates for multi-source beamformers and 56–73% rates for their single source counterparts. For both moderate and weak sources multi-source localizers exhibited smaller location errors than the single source ones the difference being strongly statistically significant (p-valueb 0.01 by a pair wise T-test). We should finally note that all beamformers (whether single or multi-source) could not identify the source orientation better than with 30 to 40° of error, independent of SNR. Moreover, orientation was determined with 2 to 3° of accuracy even for the low SNR case, or with better than 0.5 degree of accuracy for the high SNR case if all sources coincide with the grid voxels. As a note, sources likely do not coincide with grid voxels in practice, and thus we suggest that source orientations found by the beamformers in question would not be reliable. MEG AEF data analysis To test the methodology described above with real measurements, data obtained in a group study of an auditory steady state response (ASSR) were analyzed. It is well known that in such experiments sustained highly correlated bilateral neural activations are observed making it difficult for traditional single source beamformers (Herdman et al., 2003; Popescu et al., 2008). The data from 11 typically developed human subjects were collected using a 151-channel MEG scanner (CTF/VSM) with a frequency range of 0–150 Hz and a sample rate of 600/s. Tones of 500 ms duration were randomly presented at 1000 and 1200 Hz with a sound level of 80 dB SPL (sound pressure level) modulated by 100% depth at 40 Hz. After preprocessing and artifact removal, about 200 epochs of the MEG data

491

for each subject were left for analysis. Sample covariance matrices R were calculated using 300 ms segments of data starting 40 ms after the stimulus onset. To estimate the noise covariance matrix N for each participant, 500 ms “baseline” windows preceding the stimulus were used. For the purposes of group analysis, the following source search procedure was applied. On each step of the algorithm, individual subject's beamformer distributions were normalized on their mean over the corresponding ROI, transformed to a standard MNI brain space using SPM99 software package, and then averaged (Singh et al., 2002). The normalization prevented individual distributions with large dynamic range from over-influencing the group results. The strongest maximum of the group-averaged data was selected as the next source location. It was then transformed back to each subject's brain space. If an individual functional image had maxima within 2 cm distance from the projected group location, the closest one was used as the next source position for this participant. If no such maxima were found, the projected group location itself was used. Such an approach ensured that beamformer distributions with anatomically aligned null-gain constraints were averaged on each step of the algorithm. All four beamformers discussed above were applied to the data. However, power-based ones, especially MAI, showed poor performance and were able to reliably localize only one source in right auditory area (Talairach coordinates ({55, −24, 6} mm)). This was not surprising considering that overall power of the steady-state signal was relatively small. On the other hand, event-related beamformers were successful in locating both right and left auditory sources (positioned at {58, −25, 1} and {−56, −24, 10} mm, respectively), as shown in Fig. 7A. On the first iteration (for single-source event-related beamformer) the left source was not detected, but it became visible on the second step when coherent interference from the right auditory activation was removed. On individual subjects, the left

Fig. 7. ASSR experiment. (A) — Coronal and axial cross sections of the MER beamformer distribution through the right auditory source location; normalized on their global maximum values, with reference locations excluded. Black circles mark the right source and the projection of the left source on the section plane. (B) — Reconstructed amplitudes of sources pffiffiffiffiffiffi found in the first two iterations: light gray — iteration 1, and dark gray — iteration 2. (C) — Single source pseudo-Z z k̶ for beamformer maxima found in four iterations.

492

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

auditory source was either missed or misplaced by ~ 1 cm when a single source beamformer was used, but was found when multisource filter was applied. In other words, the behavior demonstrated in simulations was also observed in the results using the real data. Fig. 7B shows average source amplitudes for two steps of the algorithm and Fig. 7C shows single-source pseudo-Z's for 1st four iterations. Iteration 2 showed an increase of ~ 50% in amplitude of the ASSR signals with the removal of cancelation effects. Fig. 7C shows that after four steps all significant phase-locked activity was extracted and iterations could be stopped. Discussion The beamformers introduced above provide means for reconstructing multi-source activity because they take into account source correlations and (at least, ideally) give accurate estimates of source locations and orientations. They differ from the (also MCMV-based) beamformers previously used (Dalal et al., 2006; Popescu et al., 2008; Quraan and Cheyne, 2010; Hui et al., 2010; Diwakar et al., 2011) in both theoretical and practical implications, as detailed below. Theoretical implications First, beamformers introduced in this paper are unbiased multisource estimators of source locations and orientations. Localizers applied by previous researches (Dalal et al., 2006; Popescu et al., 2008; Quraan and Cheyne, 2010) are in fact single-source ones, but are based on beamformer weights derived with additional null constraints. Generally these estimators are biased. Diwakar et al., 2011 used a multi-source localizer but its unbiased property was not established. We should note here that for traditional single source spatial filters the unbiased property albeit theoretically attractive might not be of practical advantage, because these beamformers still yield biased estimates in the presence of other (even uncorrelated) sources. This fact cannot be changed by say using a finer grid or more accurate estimates of covariance matrices. For multiple source beamformers the situation is different. We showed that in principle one must get true source locations and orientations irrespective of SNR, assuming that covariances R and N are known exactly and provided that an exhaustive search is performed. When conditions are not ideal, we can still count on getting closer to the correct solution with more accurately estimated R, N matrices and/or using more elaborate multidimensional search procedures. Second, in our approach we determine optimal orientation for each source, which is analogous to the scalar beamformer in the single source scenario. Other MCMV-based approaches utilized the vector beamformer (Dalal et al., 2006; Popescu et al., 2008; Quraan and Cheyne, 2010; Hui et al., 2010). The latter is known to have lower signal to noise ratio in an important practical case where the source orientation is fixed (Sekihara et al., 2004; Sekihara and Nagarajan, 2008; Robinson and Vrba, 1999). Dual core beamformer (Diwakar et al., 2011) identifed source orientations but only for the special case of two highly correlated sources and large enough SNRs. Third, we proved the saturation property of the new localizers which has no analogs in single source beamformers. In the multisource approach, the number of sources n becomes an additional independent parameter that the weights and the localizer function depend on. When the saturation is reached we know that there are no more sources to search for (see Figs. 2 and 6B). The results presented in here have three implications for traditional single source beamformers (a special case with n = 1). First, if only one source exists, both
neural  activity index ζ =
T T (h TN − 1h)/(h TR − 1h) and pseudo-Z Z ̶ = h R−1 h = h R−1 NR−1 h provide unbiased estimates of location and orientation for an arbitrary noise covariance matrix N, which was proven for the pseudo-Z localizer in the case of diagonal white noise (Sekihara et al., 2005;

Greenblatt et al., 2005). Our result extends this finding to any noise matrices. To our knowledge, the same property for ζ is established here for the first time. In particular for diagonal white noise ζ becomes a scaled version of the PNFS estimator, which under white noise assumption has no bias (Sekihara et al., 2005). However for arbitrary noise covariance ζ and PNFS are no longer equivalent and the unbiased property for PNFS does not hold. At the same time, approximating true noise covariance with a diagonal one can lead to large localization errors, as discussed earlier (see simulation I, Fig. 3). Second, we present a new pseudo-Z definition for the LCMV beamformer which is the PMPZ expression (10) calculated by using lead field matrices HLCMV = {hx, hy, hz}. The LCMV (“vector”) beamformer (Van Veen et al., 1997) has lower spatial resolution than the “scalar” one (Robinson and Vrba, 1999; Sekihara et al., 2004) but is often used because knowledge of the source orientation is not required. Several definitions of pseudo-Z type expressions for the vector beamformer were suggested in previous publications (Sekihara and Nagarajan, 2008; Huang et al., 2004; Quraan and Cheyne, 2010), however to our knowledge there is no analytical proof that these definitions provide unbiased estimates of the source location. Simulations confirm that each definition typically allows locating sources with accuracy compatible with the beamformer peak spread but does not necessarily yield precise locations, as one would expect when all model parameters are known exactly. On the other hand, expression (10) does provide a localizer definition that is unbiased because vector beamformer is just a special case of the MCMV filter. Third, it follows from our results that the source location found by an event-related beamformer is unbiased provided that the source orientation is found based on the evoked field — that is by maximizing a single
 T source version of expression (11): PER = E = T = h R−1 RR−1 h =
T h R−1 NR−1 h . To our knowledge, the unbiased nature of event-related beamformers has not been analytically proven before. Mind also that common approaches are to find source orientations using power-based localizers (Robinson, 2004; Cheyne et al., 2007) and then applying those to evoked fields to identify the maxima. Assuming no other sources exist, the power-based beamformer should find the same (true) source orientation, but in practice the SNR of power-based localizer can be quite poor (this is the reason why the event-related beamformers were developed). If other sources exist orientation found by the power beamformer is affected by the oscillatory activations that do not have a phase-locked component; therefore it might be determined inaccurately and result in a localization error. On the contrary, direct maximization of PMER via solution of the eigenvalue problem fully utilizes higher SNR of the event-related beamformer. The same considerations hold for PrMER beamformer.

Practical implications Several things must be mentioned regarding application of the suggested multi-source beamformers in practice. Although theoretically all sources can be found using exhaustive search in the 5n-dimensional parameter space, it is computationally infeasible with current computing resources found in most laboratories. We introduced a source search algorithm that makes the MCMV approach practically possible without making a priori assumptions about the number of sources and their locations. However this procedure is approximate. As a result, even though the unbiased property of all localizers holds regardless of the signal to noise ratio, SNR does matter. It determines the spatial resolution and strongly affects the performance of the search procedure. For poor SNR, the location of the strongest maximum of a single source beamformer distribution might be significantly off the true source position which will adversely affect the performance of subsequent steps for all localizers, as we saw in simulation III (Fig. 5). This is also reflected in lower success rates and larger average localization errors obtained in random simulations for weaker sources. Moreover, for low SNRs the

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

procedure is more vulnerable to errors in estimates of matrices R and N (Brookes et al., 2010). Another practical consideration relates to the choice of the noise covariance matrix N. When used in the beamformer equations, an incorrect noise covariance matrix results in localization errors. This happens because anything that is not included in the noise field is automatically treated as a source signal by the beamformer. For example, if matrix N is only based on the sensor noise, then spontaneous neural activations are effectively added to the sources in the search space. In low SNR cases these activations might be randomly located by the search algorithm and seen as real sources instead of the true sources of interest (see Fig. 3). We suggest using rMER localizer for evoked activity for low SNR situations when the true noise covariance cannot be estimated accurately. For any iterative algorithm, an important question is at what point the iterations should be stopped. At first glance, the saturation effect looks like a good criterion (Fig. 2) but it turns out to be sensitive even to small localization errors. In other words, all sources might already be found accurately for all practical purposes, with the saturation still not being established (see i.e. MER plot in Fig. 6A). At the same time, the beamformer order n cannot be increased indefinitely. In practice, it should be small compared to the number of channels M. Larger n reduces the number of degrees of freedom (M–n) available for suppression of interfering sources, eventually leading to performance degradation. In our experience, the most reliable way of controlling the pffiffiffiffiffiffi iterations is to estimate the single-source pseudo-Z ( z ̶k ) for each new localizer maximum. We currently suggest that experimenter reviews this quantity for several iterations to identify which source to include, based on fluctuations around sources with the lowest SNR. For example, in Fig. 5B sources 1, 3 and 6 in MAI would be considered real sources while 2 and 4 are not, and 5 is possible source (especially if there are any additional considerations to include it). The search should be stopped when new sources are not significantly above the established baseline level. The practical aspects discussed above are illustrated by the real ASSR group data analysis. In this case, the low SNR of the steady state response due in part to the source correlations made it impossible for the power-based localizers to reliably identify even the strongest (i.e. the right auditory) activation therefore the search algorithm failed. Note that a more sophisticated approach using narrow band filtering could make power-based beamformers succeed (Herdman et al., 2003) but this possibility was not explored in the present paper. Application of event-related beamformers succeeded even using the broad band data and the results are generally consistent with those obtained using event-related SAM (Cheyne et al., 2007; Herdman and Cheyne, 2009). We also observed an increase in signal amplitude by about 50% (Fig. 7B) when coherent interference was removed by the multi-source beamformer — as should be expected both based on our simulations and on the results from the literature (Popescu et al., 2008; Quraan and Cheyne, 2010; Herdman and Cheyne, 2009). In a single-source case, so called dual state beamformer (Robinson and Vrba, 1999; Vrba and Robinson, 2001; Barnes and Hillebrand, 2003) proved to be successful, for locating oscillatory sources, when eventrelated beamformers cannot be applied. In a dual state method, spatial distributions (SPMs) of a certain statistical measure of the differences in the brain activity between the conditions are constructed. Typical examples are pseudo-T and pseudo-F statistics (Robinson and Vrba, 1999; Vrba and Robinson, 2001). Maxima and minima of these SPMs, ideally, should point to the sources appearing or missing in each condition compared to the other. It would be beneficial to extend the dual state approach to a multi-source case. However, there are at least two problems complicating such an extension. First, differences in the spatial distributions of multi-source localizers are only meaningful if distributions for each condition correspond to equal beamformer orders n and identical sets of constraints (3). Only in this case the dual state SPM will reflect changes in the brain activity rather than the differences

493

between localizer functions themselves. Therefore, either identical constraints should be enforced for each condition, or subtraction of the distributions should be avoided. Second, the unbiased property of existing dual state localizers was never established. As we show below, pseudo-F estimator is unbiased if additional assumptions about the sources are made. However simulations show that generally the single source dual state localizers are biased. At the same time, the unbiased property is a crucial feature of the beamformers developed in this work and one would like to preserve it in a dual state approach. Both problems outlined above cannot be properly addressed within the limits of this paper. Nevertheless, here we demonstrate one important case where direct extension of MAI and MPZ localizers to a dual-state scenario is straightforward. Assume that brain sources observed in condition 2 consist of those in condition 1 and also new sources uncorrelated with the sources in condition 1. A typical practical example would be when a resting state condition is compared to an “active” condition exhibiting additional brain activity in response to the stimuli. One can write the following expressions for sensor covariance matrices observed in conditions 1 and 2: R1 = N + H01C1H01T; R2 = N + H01C1H01T + H02C2H02T, with C1, 2, H01, 2 denoting source covariances and lead fields for respective conditions, and N being the common noise covariance matrix. Therefore R2 = R1 + H02C2H02T. In this expression matrix R1 plays a role of the noise which is uncorrelated with the sources C2. Consequently one can construct localizers for condition 2 using covariance matrix R1 in place of the noise covariance N (see expression (A2) in Appendix A). Specifically, for the MAI beamformer one gets P MAI dual = Tr(GS − 1 ) − n = Tr[(H T R 1− 1 H)(H T R 2− 1 H) − 1 ] − n (see Table 1). The derivation in Appendix A holds, and it follows that PMAI dual constitutes an unbiased localizer for the “new” sources that are observed only in condition 2. Note that for a single source case the last expression is simplified to PMAI dual = P2/P1 − 1 where P1, 2 are the reconstructed source powers for respective conditions, and coincides with the pseudo-F statistic introduced in Robinson and Vrba, 1999. Thus PMAI dual is a generalization of the dual-state pseudo-F beamformer for a multiple source case. At the same time, this proves the unbiased property of a single source pseudo-F beamformer provided the assumptions made above are satisfied. The “dual state” variant of the MPZ beamformer has better spatial resolution than MAI and is also unbiased: PMPZ dual = Tr(ST − 1) − n = Tr[(H TR2− 1H)(H TR2− 1R1R2− 1H) − 1] − n. Note that in a single source case this expression is not reduced to pseudo-F or pseudo-T statistic. Both PMAI dual, PMPZ dual only see sources specific to condition 2 and ignore those that are common for conditions 1 and 2. This is illustrated in Fig. 8 for the MPZ dual beamformer as an example. The same foursource setup as in Simulation 3 is used (see Fig. 5 and Table 3), but this time source #2 is uncorrelated with all others, while sources #1, 3, 4 are correlated with c = 0.8. In condition 1 only source #2 is present. In condition 2 sources #1, 3, 4 also become activated. Fig. 8 demonstrates that MPZ-dual beamformer successfully finds sources #1, 3, 4 specific to the condition 2 while being completely blind to strong source #2, even though the latter is quite close to sources #3 and #4. Although dual state evoked beamformers were not considered in the literature, a logical question is — would it be possible to construct similar dual state versions of MER or rMER using the same assumptions? The answer is negative, because while oscillatory part of the signals in condition 1 is “hidden” (see Fig. 8), both MER and rMER beamformers can still see their phase-locked part. Therefore extrema of MER or rMER “dual” distributions will not only be due to the differences between the conditions. For event-related beamformers (including a single source case) other approaches should be looked for, possibly based on direct subtractions of the distributions in conditions 1 and 2. Our final note relates to application of the methods developed in this work to brain connectivity analysis (see Schoffelen and Gross, 2009 and references there in). Although single source minimum variance filters were applied for this purpose (Gross et al., 2001), the major concern here is linear mixing of the source time courses (Hui et al., 2010) leading

494

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

Fig. 8. Dual state MPZ beamformer, using the same sources as in Simulation III. Only source #2 exists in condition 1; all four sources are present in condition 2. Sources #1, 3, 4 (see Table 3) are correlated with c = 0.8 and are uncorrelated with source #2. (A) — Cross sections of beamformer distributions through the source plane z = 8 cm; normalized on their global maximum values, with reference locations excluded. True source locations are marked with crosses, circles show projections of the locations found by a beamformer. (B) — Reconstructed amplitudes pffiffiffiffiffiffi of sources found in first four iterations after subtracting projected noise: light gray — iteration 1, dark gray — iteration 4, and black — the true amplitude. (C) — Single source pseudo-Z z k̶ for the beamformer maxima.

to spurious connectivity. Namely, as a consequence of finite spatial resolution of a filter, the reconstructed source signal always contains a certain amount of other source signals mixed in. In other words, there is always a “leakage” of one source signal into another. As a result, false connections between the sources might be observed. Although such leakage is not specific to minimum variance beamformers, the latter are especially vulnerable to it if the sources are statistically dependent, because by design the minimum variance filter exploits source correlations to minimize the output power and therefore to cancel the sources. As a consequence larger source correlations result in worse beamformer SNR, diminished spatial resolution, stronger source leakage and greater signal mixing. More details regarding this problem can be found in Sekihara et al., 2002. We want to note that the application of MCMV beamformers to this area looks quite natural for two reasons. First they are not based on the source independence assumption. Second, the linear mixing is now greatly reduced (ideally-excluded) due to the null-gain constraints imposed on the interferer's locations (see Hui et al., 2010). For connectivity analysis purposes, source positions are usually preset but orientations are not known. The latter can be obtained by repeatedly solving eigenvalue problems (13), starting with the strongest source. For a small number of sources, for example when pair-wise connectivity is looked at, even an exhaustive brute force search for orientations can be applied.

Acknowledgment The Natural Sciences and Engineering Research Council of Canada (grant #341602-2007) and the Down Syndrome Research Foundation funded this research. The authors also wish to thank Dr. Harold Wilson for his helpful advice and a very fruitful discussion. Appendix A. Proof of properties of the localizers in Table 1 Table 1 summarizes the MCMV-based source localizers introduced in this work. Table 2 and the derivations of expressions (3)–(6) above define the matrices that appear in Table 1 and in this Appendix. Matrices without subscript 0 are calculated with the trial source parameters Θ; matrices with subscript 0 are calculated using the true (but unknown) source forward solutions H0. The object of the analysis will be to show that localizers calculated with H0 achieve maximum values. We start with expression (9) for the multi-source activity index:

 1 = 2 −1 1 = 2 −1 −n PMAI ðΘÞ = Tr G S G −I n ≡Tr GS

It is straightforward to show that covariance R and its inverse can be written as:

Conclusions In this work, we derived new MCMV-based source localizers of two kinds: a) those that can be applied to both induced and evoked activity; and b) those targeting only the evoked component. They can be regarded as a generalization of scalar minimum-variance beamformers to the case of multiple correlated sources. We showed that for arbitrary noise covariance both type of localizers yield simultaneous unbiased estimates of multiple source positions and orientations and remain bounded in singular points. The saturation property of the localizers places a natural limit on the beamformer order and at least ideally provides means to determine the true number of sources. In addition we proved the unbiased property of several traditional (single-source) power and event-related localizers for arbitrary noise covariance matrix N, and derived an unbiased pseudo-Z expression for a vector LCMV beamformer. To apply the multi-source localizers in practice, we proposed an iterative source search algorithm that makes it possible to estimate multiple source locations without a priori assumptions about their locations or degree of correlation. The approach presented was shown to yield good results in simulations and with real MEG data, where correlations between the sources are significant. Moreover the multi-source beamformers provided more accurate locations and higher SNR of reconstructed activity compared to traditional single-source beamformers.

ðA1Þ

R = N + H 0 CH T0 ;


−1 R−1 = N −1 −N −1 H 0 C −1 + G0 H T0 N −1

ðA2Þ

Using the expressions above, the matrix inversion lemma (A + BCD) − 1 = A − 1 − A − 1B(C − 1 + DA − 1B) − 1DA − 1 (Hager, 1989) and Table 2 we get:
−1 −1 T H = G−K C + G0 K ;
−1 ; = G−1 + G−1 K C −1 + G0 −K T G−1 K K T G−1 T

−1

S=H R S −1

ðA3Þ

with matrix K = HTN− 1H0 describing the overlap between the n-source trial forward solution and that of the n0 true sources. Then from (A1) we have 0 0 1−1 1
 1 = 2 −1 1 = 2 T −1 C T −1 = 2 C B −1 = 2 B −1 PMAI ðΘÞ = Tr G S G −I n = Tr@G K @C + G0 −K G K A K G A |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} X

ðA4Þ In Eq. (A4), consider now n0 × n0 matrix X = G0 − K TG − 1K. Using explicit formulae for matrices G, G0, K and the SVD decompositions for

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

noise-weighted forward solutions: N − 1/2H = UDV T , N − 1/2H0 = U0D0V0T leads to the following expression for X:


T T T X = V 0 D0 I n 0 − U 0 U U U 0 D0 V 0

ðA5Þ

Here U is (M×n) matrix with orthonormal columns: U TU = In; V is an n×n orthogonal matrix, D is an n×n diagonal matrix and similarly for U0, D0, V0. Clearly the positive-definite character of X is fully determined by a symmetric square matrix In0 − (U0TU)(U TU0) which is: a) Positive definite if the column spaces of U, U0 (and consequently of H, H0) do not have a common subspace. That is no forward solution hi, i = 1,…,n matches one of the true forward solutions h0j , j = 1,…, n0 b) Positive semi-definite if some of the h i match h0j c) Zero if all sources are found and H0 ∈ H. Eq. (A5) shows that the same statements can be made regarding matrix X. Let A, B be positive (semi) definite matrices. We use notation A ≥ B if matrix (A–B) is positive semi definite (and similarly for other inequalities). Then we have 1) 2) 3) 4)

For any matrix M, M TAM ≥ 0; If A − 1exists, then (A + B) − 1 ≤ A − 1; If A ≥ B then Tr(A) ≥ Tr(B); Tr(AB) = Tr(BA). This gives the following expression for PMAI:



−1 −1 = 2 −1 T −1 = 2 −1 = 2 T −1 = 2 ≤Tr G PMAI ðΘÞ = Tr G K C +X K G KCK G
 1 = 2 T −1 1=2 = Tr C K G KC

Recalling that K TG − 1K = G0 − X we get


Tr C

1=2

T

K G

−1

KC

1=2





 1=2 1=2 1=2 1=2 = Tr C ðG0 −X ÞC ≤Tr C G0 C

= Tr ðCG0 Þ max Thus we proved that P MAI (Θ) ≤ Tr(CG 0 ), P MAI (Θ) = Tr(CG 0 ) (Table 1) and the equality is only possible when X is zero. That is, when all true sources are found, the localizer achieves its maximum. Note that the limit on PMAI(Θ) holds even when sources are close (θ i → θ j) and that PMAI(Θ) is bounded at singular points. Consider now multi-source pseudo-Z PMPZ(Θ)=Tr(T− 1/2ST− 1/2 −In). Using Eq. (A2), one gets for matrix T


−1
−1
−1 −1 T −1 −1 T T = G−2K C + G0 K + K C + G0 G0 C + G0 K ðA6Þ Using Eq. (A3) to derive an expression for SG − 1S shows that T can be written as: T = SG

−1


−1
−1 
−1 T −1 −1 T S + K C + G0 G0 −K G K C + G0 K ðA7Þ

As we know matrix X=G0 −KTG− 1K is non-negative, so is the 2nd term in Eq. (A7) and we have T ≥ SG − 1S ⇔ T − 1 ≤ (SG − 1S) − 1 = S− 1GS− 1. Consequently the following estimate holds: PMPZ



 −1 −1 = Tr ST −n≤Tr GS −n = PMAI :

In other words, always PMPZ(Θ) ≤ PMAI(Θ). Equality is only achieved when X = 0 — that is when all true sources are found. This completes the proof that the global maximum of PMPZ is the true source

495

description, that PMPZ has higher spatial resolution than PMAI and that max max PMPZ = PMAI (see Table 1). For the MER localizer, PMER =Tr(ET− 1) we start by examining matrix E = H T R−1 RR−1 H. It is straightforward to show using Eq. (A2) that E = YCY T with n×n0 matrix Y given by Y=K−K(C− 1 +G0)− 1G0. Substituting Y into Eq. (A6) for T gives: T=G−KG0− 1KT +YG0− 1YT. Introducing another matrix X1 =G−KG0− 1KT gives T=X1 +YG0− 1YT. X1 has the same structure as matrix X above, and we see that X1 ≥0. Assume for the moment that n≤n0. Then matrix Y is wide or square, matrix YG0− 1YT is non-singular and consequently (YG0− 1YT)− 1exists. Therefore −1 T=X1 +YG0− 1YT ≥YG0− 1YT ⇒T ≤(YG0− 1YT
)− 1,
and we get
 the folT −1 lowing inequality: PMER = Tr YCY T T −1 ≤Tr CY T YG−1 Y . 0 Y Next use a singular value decomposition: G0− 1/2YT = UyDyVyT. Since T

Y is “tall”, Dy, Vy are square (n×n) invertible matrices, while Uy is n0 × n −1 T −1 1/2 and we can write YT(YG U U TG1/2 ≤ G . That proves
0 Y ) Y = G  0 y y 0 max 0 inequality PMER = Tr YCY T T −1 ≤Tr CG0 = PMER when n ≤ n 0 . Equality is achieved if n = n0 and matrix X1 = 0 which happens only when Θ is the true source distribution. To prove the inequality for n N n0 note that one can always add fictitious sources with zero amplitude so that n becomes smaller than n0 again and previous logic applies. At the same
time,  sources with zero amplitudes do not contribute to the trace

Tr CG0 , which proves the case for n N n0. The derivation for the rMER localizer PrMER = Tr(ES − 1) directly follows the one for MER because matrix S in Eq. (A3) can be represented as S = X1 + Y(G0− 1 + C)Y T. The inequality for PrMER is PrMER (Θ) ≤ Tr(C(G0− 1 + C) − 1) and the maximum is realized only when Θ is the true source distribution. Appendix B. Source search procedure: derivation of expressions (13), (14) On the step k of the algorithm, the H-matrix can be written as H = {H R, h k(θ k)}, where columns of “reference” matrix H R are forward solutions of the k–1 sources already found: H R = {h 1,…, h k − 1}. The next source's forward solution h k is expressed as h k (θ k) = H k(r k)u k where columns of matrix H k = {hxk, hyk, hzk} are lead fields for dipoles at location r k, oriented along axes x, y and z, and u k is the unknown source orientation. We demonstrate the derivation using the MAI beamformer as an example. To find the orientation of source k, we need to maximize the localizer functionPMAI = Tr(GS− 1) − n with respect to u k. First write matrices G, S

 g RR g Rk in a block form, namely G = where gRR = H RTN− 1HR, kR kk g g g Rk = (g kR)T = HRTN− 1hk, g kk = hkTN− 1hk, and same expressions for blocks of S, with N− 1replaced by R− 1. The block representation of the inverse matrix Q ≡ S− 1 can be obtained using Shur's decomposition formula: qRR =(sRR)− 1 +(sRR)− 1sRkqkkskR(sRR)− 1; qRk =−(sRR)− 1sRkqkk; qkk =(skk −skR(sRR)− 1sRk)− 1. Noting that only diagonal blocks of the product GS − 1 ≡GQ will contribute to the trace, we get PMAI = Tr (g RRq RR + g Rkq kR +g kRq Rk + g kkq kk) − n. Substituting expressions for q * * and noticing that a) Tr(g RR(s RR) − 1) − n does not depend on u k; b) q kk is simply a number and consequently can be pulled out of the trace; c) Tr(AB) =Tr(BA) — we arrive to the following maximization problem. Find

max




−1
−1 kk kR RR RR RR Rk kR Rk kR Rk kk kk Tr q s s g s s +q g +g q +g q =

j uk j = 1

2
−1 3
−1
−1
−1 kR RR RR RR Rk kR RR Rk kR RR Rk kk s s g s s −g s s −s s g +g 6 7 max 4 5  −1 Rk skk −skR sRR s j uk j = 1

ðB1Þ

In Eq. (B1), we can make dependence on uk explicit by replacing hk with H k(r k)u k and introducing g kR = u kTH kTN − 1H R = u kTG kR with

496

A. Moiseev et al. / NeuroImage 58 (2011) 481–496

GkR =HkTN− 1HR (and similarly for other non-diagonal elements of blockmatrices G* *, S* *). To have a uniform notation, also denote SRR ≡sRR, GRR ≡gRR. Then problem (B1) is reduced to finding a maximum of a ratio of two quadratic forms:  3 2
−1
−1
−1
−1 kT S kR S RR GRR S RR S Rk −GkR S RR S Rk −S kR S RR GRk + Gkk uk 7 6u 7
 max 6   4 5 j uk j = 1 ukT S kk −S kR S RR −1 S Rk uk = max

j uk j = 1

ukT DMAI uk ukT F MAI uk

ðB2Þ

According to a well known mathematical result, the optimum u k in Eq. (B2) is obtained by solving the eigenvalue problem Eq. (13). Expressions for other beamformers are obtained almost identically, leading to the generic formula Eq. (14). References Baillet, S., Mosher, J.C., Leahy, R.M., 2001. Electromagnetic brain mapping. IEEE Signal Process. Mag. 18, 14–30. Barnes, R.B., Hillebrand, A., 2003. Statistical flattening of MEG beamformer images. Hum. Brain Mapp. 18, 1–12. Brookes, M., Stevenson, C., Barnes, G., Hillebrand, A., Simpson, M., Francis, S., Morrisa, P., 2007. Beamformer reconstruction of correlated sources using a modified source model. NeuroImage 34 (4), 1454–1465. Brookes, M.J., Zumer, J.M., Stevenson, C.M., Hale, J.R., Barnes, G.R., Vrba, J., Morris, P.G., 2010. Investigating spatial specificity and data averaging in MEG. NeuroImage 49 (1), 525–538. Cheyne, D., Bostan, A.C., Gaetz, W., Pang, E.W., 2007. Event-related beamforming: a robust method for presurgical functional mapping using MEG. Neurol. Clin. Neurophysiol. 118, 1691–1704. Dalal, S., Sekihara, K., Nagarajan, S., 2006. Modified beamformers for coherent source region suppression. IEEE Trans. Biomed. Eng. 53 (7), 1357–1363. Diwakar, M., Huang, M.-X., Srinivasan, R., Harrington, D.L., Robb, A., Angeles, A., Muzzatti, L., Pakdaman, R., Song, T., Theilmann, R.J., Lee, R.R., 2011. Dual-Core Beamformer for obtaining highly correlated neuronal networks in MEG. NeuroImage 54, 253–263. Frost, O.L., 1972. An algorithm for linearly constrained adaptive array processing. Proc. IEEE, 60, no. 8, pp. 926–935. Greenblatt, R.E., Ossadtchi, A., Pflieger, M.E., 2005. Local linear estimators for the bioelectromagnetic inverse problem. IEEE Trans. Signal Process. 53, 3403–3412. Gross, J., Kujara, J., Hamalainen, M., Timmermann, L., Schnitzler, A., Salmelin, R., 2001. Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc. Nat. Acad. Sci. 98, 694–699. Hager, William W., 1989. Updating the inverse of a matrix. SIAM Rev. 31 (2), 221–239. Hamalainen, M., Hari, R., Ilmoniemi, R.J., Knuutila, J.L., Olli, V., 1993. Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Modern Phys. 65, 413–497. Herdman, A.T., Cheyne, D., 2009. A Practical guide for MEG and beamforming. In: Handy, T.C. (Ed.), Brain Signal Analysis: Advances in Neuroelectric and Neuromagnetic Methods. MIT Press, Cambridge, MA. Herdman, A.T., Wollbrink, A., Chau, W., Ishii, R., Ross, B., Pantev, C., 2003. Determination of activation areas in the human auditory cortex by means of synthetic aperture magnetometry. NeuroImage 20, 995–1005. Herdman, A.T., Wollbrink, A., Chau, W., Ishii, R., Ross, B., Pantev, C., 2004. Localization of transient and steady-state auditory evoked responses using synthetic aperture magnetometry. Brain Cogn. 54 (2), 149–151. Hillebrand, A., Barnes, G.R., 2005. Beamformer analysis of MEG data. Int. Rev. Neurobiol. 68, 149–171 (Special volume on Magnetoencephalography).

Huang, M.-X., Shih, J.J., Lee, R.R., Harrington, D.L., Thoma, R.J., Weisend, M.P., Hanlon, F., Paulson, K.M., Li, T., Martin, K., Miller, G.A., Canive, J.M., 2004. Commonalities and differences among vectorized beamformers in electromagnetic source imaging. Brain Topogr. 16 (3), 139–158. Hui, H.B., Pantazis, D., Bressler, S.L., Leahy, R.M., 2010. Identifying true cortical interactions in MEG using the nulling beamformer. NeuroImage 49, 3161–3174. Hymers, M., Prendergast, G., Johnson, S.R., Green, G., 2010. Source stability index: a novel beamforming based localisation metric. NeuroImage 49, 1385–1397. Jerbi, K., Mosher, J.C., Baillet, S., Leahy, R.M., 2002. On MEG forward modelling using multipolar expansions. Phys. Med. Biol. 47, 523–555. Kimura, T., Kako, M., Kamiyama, H., Ishiyama, A., Kasai, N., Watanabe, Y., 2007. Inverse solution for time-correlated multiple sources using Beamformer method. Int. Congr. Ser. 1300, 417–420. Kirsch, H.E., Robinson, S.E., Mantle, M., Nagarajan, S., 2006. Automated localization of magnetoencephalographic interictal spikes by adaptive spatial filtering. Clin. Neurophysiol. 117, 2264–2271. Limpiti, T., Van Veen, B.D., Wakai, R.T., 2006. Cortical patch basis model for spatially extended neural activity. IEEE Trans. Biomed. Eng. 53 (9), 1740–1754. Mosher, J.C., Leahy, R.M., 1998. Recursive MUSIC: a framework for EEG and MEG source localization. IEEE Trans. Biomed. Eng. 45 (11), 1342–1354. Mosher, J.C., Leahy, R.M., 1999. Source localization using recursively applied and projected (RAP) MUSIC. IEEE Trans. Sign. Proc. 47 (2), 332–340. Popescu, M., Popescu, E., Chan, T., Blunt, S., Lewine, J., 2008. Spatial-temporal reconstruction of bilateral auditory steady state responses using MEG beamformers. IEEE Trans. Biomed. Eng. 55 (3), 1092–1102. Quraan, M., Cheyne, D., 2010. Reconstruction of correlated brain activity with adaptive spatial filters in MEG. NeuroImage 49, 2387–2400. Robinson, S.E., 2004. Localization of event-related activity by SAM(erf). Neurol. Clin. Neurophysiol. 109. Robinson, S., Vrba, J., 1999. Functional neuroimaging by synthetic aperture magnetometry (SAM). Recent Advances in Biomagnetism. Tohoku Univ. Press, Sendai, Japan, pp. 302–305. Schoffelen, J.-M., Gross, J., 2009. Source connectivity analysis with MEG and EEG. Hum. Brain Mapp. 30 (6), 1857–1865. Sekihara, K., Nagarajan, S.S., 2008. Adaptive spatial filters for electromagnetic brain imaging. Springer, Berlin. Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., Miyashita, Y., 2001. Reconstructing spatio-temporal activities of neural sources using an MEG Vector beamformer technique. IEEE Trans. Biomed. Eng., 48, no. 7, pp. 760–771. Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., 2002. Performance of an MEG adaptive-beamformer technique in the presence of correlated neural activities: effects on signal intensity and time-course estimates. IEEE Trans. Biomed. Eng. 49 (12), 1534–1546. Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., 2004. Asymptotic SNR of scalar and vector minimum-variance beamformers for neuromagnetic source reconstruction. IEEE Trans. Biomed. Eng., 51, No 10, pp. 1726–1734. Sekihara, K., Sahani, M., Nagarajan, S., 2005. Localization bias and spatial resolution of adaptive and non-adaptive spatial filters for MEG source reconstruction. NeuroImage 25, 1056–1067. Singh, K.D., Barnes, G.R., Hillebrand, A., Forde, E.M.E., Williams, A.L., 2002. Task-related changes in cortical synchronization are spatially coincident with the hemodynamic response. NeuroImage 16, 103–114. Stoica, P., Sharman, K.C., 1990. Maximum likelihood methods for direction-of-arrival estimation. IEEE Trans. Signal Process. 38, 1132–1143. Van Veen, B., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867–880. Vrba, J., Robinson, S.E., 2001. Signal processing in magnetoencephalography. Methods 25 (2), 249–271. Widrow, B., Duvall, K.M., Gooch, P.P., Newman, W.C., 1982. Signal cancellation phenomena in adaptive arrays: causes and cures. IEEE Trans. Antennas Propag., vol. AP-30, no. 5, pp. 469–478. Williamson, S.J., Kaufman, L., 1981. Biomagnetism. J. Magn. Magn. Mater. 22, 129–201.