Beamformer reconstruction of correlated sources using a modified source model

www.elsevier.com/locate/ynimg NeuroImage 34 (2007) 1454 – 1465

Beamformer reconstruction of correlated sources using a modified source model Matthew J. Brookes, a,⁎ Claire M. Stevenson, a Gareth R. Barnes, b Arjan Hillebrand, b Michael I.G. Simpson, b Susan T. Francis, a and Peter G. Morris a a

The Sir Peter Mansfield Magnetic Resonance Centre, School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD, UK b The Wellcome Trust Laboratory for MEG Studies, Neuroscience Research Institute, Aston University, Birmingham, B4 7ET, UK Received 17 May 2006; revised 20 October 2006; accepted 2 November 2006 Available online 29 December 2006 This paper introduces a lead field formulation for use in beamformer analysis of MEG data. This ‘dual source beamformer’ is a technique to image two temporally correlated sources using beamformer methodology. We show that while the standard, single source beamformer suppresses the reconstructed power of two spatially separate but temporally correlated sources, the dual source beamformer allows for their accurate reconstruction. The technique is proven to be accurate using simulations. We also show that it can be used to image accurately the auditory steady state response, which is correlated between the left and right auditory cortices. We suggest that this technique represents a useful way of locating correlated sources, particularly if a seed location can be defined a priori for one of the two sources. Such a priori information could be based on previous studies using similar paradigms, or from other functional neuroimaging techniques. © 2006 Elsevier Inc. All rights reserved. Keywords: MEG; Beamformer; Spatial filter; Source modelling; Correlated sources; Dual source model

Introduction Magnetoencephalography (MEG) is a direct, non-invasive technique for the measurement of human brain activity. Postsynaptic current flow within the dendrites of active neurons induces a small magnetic field that can be measured close to the scalp surface using superconducting quantum interference devices (SQUIDS) (Hamalainen et al., 1993). The magnitude of these measured fields is directly related to neuronal current strength, and hence their measurement allows characterisation of the amplitude of electrical brain activity on a millisecond time scale. The major challenge, however, involves mapping spatially the neuronal ⁎ Corresponding author. E-mail address: [email protected] (M.J. Brookes). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2006.11.012

current distribution that underlies the measured magnetic fields. The MEG inverse problem states: given a known magnetic field distribution outside the head, how do we reconstruct an accurate representation of the original neuronal current distribution inside the head? This is an ill-defined problem since any single magnetic field distribution could be caused by an infinite number of current distributions (von Helmholtz, 1853). Mathematically, the problem is ill defined in a Hadamard (1952) sense. For these reasons, while the temporal resolution of MEG is unparalleled, its spatial resolution is dependent on the approach that is used to solve the inverse problem. In recent years, the spatial accuracy of MEG has been significantly improved by the introduction of a number of source imaging techniques: minimum norm (Hamalainen and Ilmoniemi, 1984; Hamalainen et al., 1993), weighted minimum norm (Hamalainen and Ilmoniemi, 1984; Hamalainen et al., 1993) and LORETA (Pascal-Marqui et al., 1994) to name but a few. Each of these techniques relies on its own specific assumption set in order to solve the MEG inverse problem. One particular class of inverse solution techniques is known collectively as MEG beamforming (Robinson and Vrba, 1998; Sekihara et al., 2002a; Van Drongelen et al., 1996; Van Veen et al., 1997). The beamformer methodology represents a spatial filtering approach to source localisation and can be thought of as being somewhat analogous to a traditional frequency filter. In much the same way as a frequency filter can be used to extract only those components of a signal containing some particular temporal characteristic (i.e. oscillating at a predefined frequency), the beamformer spatial filter can be used to extract those signals with a specific spatial characteristic (i.e. it retains only signals originating at some pre-determined spatial location). Such measurements result in an estimate of the electrical activity (the magnitude and orientation of neuronal current) at the predefined location. By sequential application of the beamformer to a number of locations placed on a regular grid spanning the source space (i.e. the brain), a volumetric image of source power can be created (Robinson and Vrba, 1998; Van Veen et al., 1997). Furthermore, by applying the beamformer in both active and

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

control time windows (e.g. during and prior to stimulation), then the task-related change in electrical power can also be assessed (Vrba and Robinson, 2001). Current implementations of the MEG beamformer are based on two main assumptions: (i) current sources in the brain are dipolar, and (ii) time courses of two or more spatially distinct neuronal sources are orthogonal. These assumptions lead to two main limitations for beamformer reconstructions. Firstly, the quality of the source reconstruction is dependent on the extent of the source (as the area of active cortex is increased, the beamformer reconstructed source power is suppressed (Hillebrand et al., 2005a; Vrba, 2002). Secondly, if two spatially distinct sources become correlated in time, then the beamformer reconstructed source power is attenuated and for highly correlated sources, little or none of the original source power will be recovered (Van Veen et al., 1997; Sekihara et al., 2002b). In practice, it is unlikely that a high degree of source correlation would persist for more than a few seconds (Hadjipapas et al., 2005; Hillebrand and Barnes, 2005), and therefore the effect of correlated sources can be minimised by ensuring that analyses are done using large time/frequency windows. In addition, the use of short time/frequency windows leads to inaccurate estimation of the data covariance matrix. However, complete suppression can still remain for some cortical areas with high levels of connectivity, particularly in the case of driven or steady state responses (Herdman et al., 2003a). The suppression of source power in this way is problematic, particularly in the case where accurate reconstruction of source amplitude is required. For example, many recent studies (see Brookes et al., 2005; Singh et al., 2002; or Preissl, 2005 for a review) have looked into correlations between MEG and functional magnetic resonance imaging (fMRI). The aim of such studies is to compare spatially and temporally the recorded blood oxygenation level dependent (BOLD) response in fMRI and the detectable neuromagnetic effects in the MEG. However, if quantitative information is to be gained from the MEG (for example, if signal amplitude is to be extracted and used as an input model to predict the BOLD effect), we have to take into account the effects of amplitude reduction due to correlated sources. Theoretical methods have been proposed to reconstruct the true source amplitudes (Sekihara et al., 2001); however, these methods do not also allow for the localisation of such sources. In the present paper, we introduce a simple but effective technique to reconstruct two (or more) correlated sources using a variation to the beamformer lead-field model. This variation is similar to that described previously in the context of MUSIC source localisation (Mosher and Leahy, 1998). The theory of the beamformer spatial filter is first summarised, with the emphasis placed on the so-called ‘non-linear beamformer’ otherwise known as synthetic aperture magnetometry (SAM) (Robinson and Vrba, 1998). Following this, an explanation of the correlated source problem is put forward, and a solution introduced in the form of a dual source lead field model. Both the single and dual source beamformers are tested using a simulation, and the dual source model is shown to be successful in locating two correlated sources with equal amplitude. Subsequently, source space images of correlated sources are derived based on both single and dual source models for a variety of different source correlations, in order to characterise the dual source beamformer performance in terms of both localisation and time course reconstruction. Simulations are also undertaken to show that the dual source

1455

beamformer works in the case where correlated sources have different amplitudes. Finally, the new dual source beamformer is applied to experimental data, and we show that the auditory steady state response (Kuwada et al., 1986; Rees et al., 1986; Ross et al., 2000), correlated across hemispheres, can be imaged accurately using this approach. Furthermore, we show that these results are consistent with those obtained using a modified single source beamformer approach. Theory The non-linear beamformer The beamformer spatial filter involves the calculation of a weighted sum of magnetic field measurements made at each of the n MEG sensors. If m(t) represents an n dimensional column vector of magnetic field measurements made at time t, and Wθ is an n dimensional column vector of weighting parameters for some predefined source space location and orientation, then the reconstructed source strength estimate Qˆ θ (t) is given by: Qˆ θ ðtÞ ¼ W Tθ mðtÞ

ð1Þ

Where the 6-element vector θ represents both the source location and orientation and the superscript T indicates a matrix transpose. The weighting parameters (Wθ) are derived based on power minimisation. The overall amount of power emanating from within the source space is minimised with the linear constraint that any power originating from the location of interest (θ) must remain in the output signal: j 2k minW θ Qˆ θ subject to W Tθ Lθ ¼ 1 ð2Þ Where Qˆ 2θ represents source power for location and orientation θ and is given by Qˆ 2θ = WTθ CWθ. Here, C represents the n × n data covariance matrix calculated over a time frequency window of interest (Robinson and Vrba, 1998; Van Veen et al., 1997; Barnes and Hillebrand, 2003; Hillebrand et al., 2005b), Lθ is the lead field vector, which is an n dimensional column vector containing the magnetic fields that would be measured at each of the n MEG sensors in response to a source of unit amplitude with location and orientation specified by θ. In most cases, Lθ is based on a dipolar model of neuronal current and computed using an appropriate model for the volume conductor (Sarvas, 1987). The solution to the beamformer problem is given by (Robinson and Vrba, 1998; Van Veen et al., 1997; Hillebrand et al., 2005b): h i1 W Tθ ¼ LTθ fC þ lΣg1 Lθ LTθ fC þ lΣg1 ð3Þ Σ is a diagonal matrix representing the white noise at each of the MEG channels and μ is a Backus–Gilbert regularisation parameter and is used to adjust the trade off between the full width at half maximum of the point spread function of the beamformer image (i.e. spatial resolution), and the magnitude of uncorrelated noise in the time course reconstruction. Assuming that the magnitude of uncorrelated sensor noise is equivalent across MEG sensors, we can write that μΣ = μσ2I, where σ2 is the white noise variance and I represents the identity matrix. Eq. (3) shows that the weighting parameters are dependent on the lead field model, the data covariance, the signal-to-noise ratio of the recorded data and the matrix regularisation parameter μ. Using Eq. (3), the

1456

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

mean squared source power at location and orientation θ can be shown to be given by: 2 Qˆ θ

¼

W Tθ CW θ

¼

LTθ fC þ lΣg1 CfC þ lΣg1 Lθ 2

fLTθ fC þ lΣg1 Lθ g

ð4Þ

In the unregularised case, i.e. μ = 0, this reduces to:  1 2 Qˆ θ ¼ LTθ C1 Lθ

ð5Þ

In principle, Eq. (4) or (5) could be calculated at a number of source locations spanning the source space in order to produce a whole head volumetric image of brain electrical activity. However, the signal-to-noise ratio of source estimates declines with depth into the head, and so close to the centre of the brain, the noise power may dominate over any genuine source power. For this reason, source power estimates must be normalised by an estimate of noise power ˆv 2θ = WTθ ΣWθ. The normalised power estimate is therefore given by: 2 Qˆ W T CW θ : Z¯ θ ¼ 2θ ¼ Tθ vˆ θ W θ ΣW θ

ð6Þ

2 Qˆ θ  T¯ θ ¼ 2 ðaÞ vˆ þ θ

2

Qˆ θ ðcÞ vˆ 2 : θ ðcÞ

A limitation of SAM is that spatially separate, yet temporally covariant, sources are suppressed. This is because for two temporally correlated sources, the forward solution at each MEG sensor will add linearly, and this vector sum may then not correspond to the lead field pattern from either of the two individual source locations. Power minimisation subject to the unity pass band will lead to rejection of both sources, and therefore cause suppression of the reconstructed power at the two source locations (Hillebrand et al., 2005b). Mathematically, this problem can be understood in terms of singular value decomposition of the data covariance matrix. We now restrict our analysis to the unregularised beamformer and consider singular value decomposition of the data covariance matrix thus: C ¼ USVT :

ð7Þ

Where (a)Qˆ 2θ represents the beamformer power estimate during an ‘active’ time window, (c)Qˆ 2θ represents the beamformer power estimate during a ‘control’ window and (a)ˆv 2θ and (c)ˆv 2θ represent noise power estimates during the active and control periods respectively. Again normalisation by a noise metric is used to account for the inherent non-uniform noise projection. To compute volumetric T¯ or Z¯ statistical images, a separate set of weights is computed for a number of possible source locations placed on a regular grid covering the entirety of source space. The orientation of potential sources is derived based on an optimisation of source power. Since, using a spherical volume conductor model, MEG is insensitive to radial dipoles, the source orientation is rotated in a plane tangential to the radial direction. The angle for which the pseudo-Z-statistic is maximised is then used as the optimum orientation for the source strength estimate at each location. Mathematically, this can be expressed as: 2 ! Qˆ r;d Z¯ opt ¼ max 2 ; 0o ≤δ ≤180o ð8Þ d vˆ r;d where r represents source location and δ is the angle of the source orientation in the tangential plane with respect to the azimuthal direction (i.e. θ = (r,δ)). (Notice that Z¯(δ) is entirely equivalent for 0° ≤ δ ≤ 180° and 180° ≤ δ ≤ 360°.) Assessments of cortical poweror stimulus-related change in cortical power are made at each location in source space, for the optimum source orientation defined using Eq. (8). The sampling of source space in this way is known as non-linear beamforming or synthetic aperture magnetometry (SAM) (Robinson and Vrba, 1998).

ð9Þ

Where, since C is both square, and symmetric, the left and right singular vectors are equivalent (i.e. U = V). S is the diagonal matrix containing the singular values. The inverse of C can be written as: C1 ¼ US1 UT ¼

This is termed either the pseudo-Z-statistic, or neural activity index (Robinson and Vrba, 1998; Vrba and Robinson, 2001). In order to measure the change in cortical power related to the performance during some specific task, the pseudo-T-statistic (Vrba and Robinson, 2001) can be employed: ðaÞ

Failure of the SAM beamformer and the dual source extension

n X

1

T

Sði;iÞ UðiÞ UðiÞ

ð10Þ

i¼1

Where S(i,i) corresponds to the ith diagonal element of S and U is the ith column of U. Put another way, this means that the covariance matrix, and its inverse, effectively represents a weighted sum of its individual singular vectors. The singular vectors themselves correspond to the field pattern produced by each of the temporally orthogonal sources active during the recording. The singular values give relative weightings to those sources. Considering only the first few singular values in the single value decomposition, it is simple to show that the beamformer power is given by (Gross et al., 2003): (i)

" #1 T T T LTθ Uð1Þ Uð1Þ Lθ LTθ Uð2Þ Uð2Þ Lθ LTθ Uð3Þ Uð3Þ Lθ 2 ˆ Qθ ¼ þ þ þ ::: Sð1;1Þ Sð2;2Þ Sð3;3Þ ð11Þ Eq. (11) informs us as to how the beamformer can act to suppress two correlated sources. Let us consider for the moment the situation in which two spatially distinct sources, Qθ1 and Qθ2, exist. If these sources are orthogonal (i.e. uncorrelated in time and giving rise to orthogonal field distributions) then (assuming a high signal-to-noise ratio) the first singular vector, U(1), will be indicative of the field pattern produced by Qθ1 and the second singular vector, U(2), will be indicative of the field produced by Qθ2. Subsequent singular vectors will be representative of the measured noise. It is clear that in the case where Lθ ⇒ Lθ1 (i.e. the field produced by a test source at θ matches Lθ1 and therefore U(1)), then the first term in Eq. (11) will be maximised, and a high output power would be observed in the final image. Similarly, as Lθ ⇒ Lθ2, d output power is again high. For cases where the singular vectors are not representative of the dipolar field patterns described by the lead fields Lθ, the output power remains small. For this reason, sources that do not originate from within the brain, or that do not produce dipolar field patterns, are suppressed. In the correlated source case, the time courses of Qθ1 and Qθ2 are no longer orthogonal, and the first singular vector U(1) will incorporate field patterns from both Qθ1 and Qθ2. Since the sum of

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

the fields produced by the two correlated sources (now represented by U(1)) does not match either Lθ1 or Lθ2, output power in the final image will always be small, and therefore the two sources will be suppressed in the final source power estimate. Eqs. (6) and (7) show that this remains the case for pseudo-Z- and pseudo-Tstatistical images. In this situation, there exist several techniques by which to retrieve the original source power. The simplest technique involves the regularisation parameter. A heavily regularised beamformer will be less susceptible to cancellation than the unregularised beamformer. However, such techniques cause a loss of spatial resolution, and may not result in an accurate estimation of source power. An attractive method would involve the modification of the covariance matrix such that Qθ1 and Qθ2 are once again orthogonal in their singular vector representation. This has been attempted by the use of higher order covariance matrices (Huang et al., 2004) and has shown promise in localising correlated activity both in simulation and in a median nerve experiment. A simpler approach, however, is to modify the lead field vector such that it is based on 2 dipoles instead of one. In this ‘dual source model,’ we assume that two correlated sources, Qθ1 and Qθ2, exist where θ1 and θ2 represent the position and orientation of sources 1 and 2, respectively, with lead field patterns Lθ1 and Lθ2, respectively. In our dual source beamformer implementation, the lead field pattern required should be computed based on both θ1 and θ2 such that: 1 1 LDUAL ¼ Lθ1 þ Lθ2 2 2

ð12Þ

By analogy with Eq. (3), the beamformer weights can be calculated as: h i1 1 1 W Tθ1 ;θ2 ¼ LTDUAL fC þ lr2 Ig LDUAL LTDUAL fC þ lr2 Ig ð13Þ Once calculated, these weights can be used to derive a source power estimate, pseudo-Z-statistic, or pseudo-T-statistic in exactly the same way as that described above by replacing Wθ with Wθ1,θ2. In this case, the dual source beamformer minimises energy originating from locations other than the two defined by θ1 and θ2 simultaneously. The optimum angle for the dipole orientation can again be found by rotation of the dipoles in the tangential plane, and nonlinear optimisation of the pseudo-Z-statistic. However, in this case, since two sources are to be modelled, unless source orientations are known a priori, two separate rotations are required and the search for an optimal angles becomes two dimensional. This can be expressed mathematically as:   o 0 ≤δ ≤ 180o Z¯ DUAL opt ¼ max Z¯ θ ;θ ; o 1 ð14Þ 1 2 0 ≤δ2 ≤ 360o d1 ;d2 Where δ1 represents the orientation angle of source one with respect to the azimuthal direction in the tangential plane, and likewise δ2 represents the orientation angle of source two with respect to the azimuthal direction the tangential plane. Determination of the optimal source orientation in the dual source beamformer is different from that in the single source case. For a single source, the dipole orientation and the polarity of the reconstructed source amplitude are interchangeable. This is the case using any (single source) solution to the MEG inverse problem. For example, a source described by Q = 1nAm and δ = 0° is also

1457

equivalently well described by Q = − 1nAm and δ = 180°. This means that for a single source non-linear beamformer, Z¯(δ) is entirely equivalent for 0° ≤ δ ≤ 180° and 180° ≤ δ ≤ 360°. However, using the dual source beamformer, two sources are modelled simultaneously and at any single instant in time, the reconstructed source amplitude must represent the sum of the individual amplitudes from each source. For this reason, the relative orientation difference between the two sources must be correctly estimated. For example, suppose that source one can be described by Q1 = 1nAm and δ1 = 0° and source two can be described by Q2 = 1nAm and δ2 = 90°. In the beamformer reconstruction, an accurate image will be obtained if source estimates are given by Qˆ 1 = 1nAm, ˆδ 1 = 0°; and Qˆ 2 = 1nAm, δˆ 2 = 90° or equivalently by Qˆ 1 = − 1nAm, δˆ 1 = 180°; and Qˆ 2 = −1nAm, δˆ 2 = 270°. However, in the case where Qˆ 1 = 1nAm, δˆ 1 = 0°; and Qˆ = −1nAm, δˆ 2 = 270°, or equivalently Qˆ 1 = − 1nAm, δˆ 1 = 180°; and Qˆ 2 = 1nAm, δˆ 2 = 90°, the relative orientation between the two sources is incorrectly estimated and therefore the source amplitudes (and source power) will subtract, meaning that a correct reconstruction will not be obtained. For the above reasons, the quantity δ1 should be varied between 0° and 180° and for each value of δ1, δ2 must be varied between 0° and 360°. This technique ensures that the relative orientation difference between the two sources is correct. Used in this way, Eq. (14) results in optimal orientation estimates for both sources at any given source locations. These optimised angles should be used in order to make estimates of source power- or stimulus-related change in source power using the dual source beamformer. An estimate of the time course of neuronal activity can also be obtained using the dual source beamformer in the same way as that described by Eq. (1) (i.e. Qˆθ1,θ2(t) = WθT1 ,θ2m(t)). The quantity Qˆθ1,θ2(t) is based on two source locations and, as described above, represents the sum of activity originating at both locations (i.e. Qˆθ1,θ2(t) ≈ Qθ1(t) + Qθ2(t)). In the case where Qθ1 and Qθ2 are perfectly correlated, the correlated time course may then be accurately expressed as Qˆθ1,θ2(t)/2. A possible problem with the dual source beamformer is that if two spatially distinct but temporally correlated sources are not equal in source amplitude then their spatial localisation will be biased. In order to explain this confound, consider again Eq. (11) but this time rewritten specifically for the dual source beamformer such that: " Qˆ θ1 ;θ2 ¼

T

LTθ1 ;θ2 Uð1Þ Uð1Þ Lθ1 ;θ2 Sð1;1Þ

#

T

þ

LTθ1 ;θ2 Uð2Þ Uð2Þ Lθ1 ;θ2 Sð2;2Þ

þ: : : ð15Þ

Here, the singular vector U(1) is indicative of the two correlated sources Qθ1 and Qθ2, and if Qθ1 ≈ Qθ2 then U(1) will be well described by the lead fields (Lθ1,θ2) such that the beamformer output power is an accurate reflection of the true source power. However, if Qθ1 ≠ Qθ2 then Lθ1,θ2 will not accurately reflect U(1). This is because Lθ1,θ2 (= Lθ1 + Lθ2) is based on dipoles of equal amplitude placed at θ1 and θ2. If for example Qθ1 << Qθ2, Lθ1,θ2 will not be accurately scaled and therefore the beamformer estimates of the source locations are likely to be incorrect, and the pseudoZ-statistics (or pseudo-T-statistics) greatly reduced. To solve the problem, we introduce a simple ‘amplitude optimisation’ routine, such that the lead field strength of the two sources can be modulated until the optimum relative source

1458

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

strengths are found. The relative lead field strengths for the two modelled sources are varied such that: LDUAL ¼ a1 Lθ1 þ a2 Lθ2

ð16Þ

In order to ensure that the correct summed source amplitude is reconstructed, these parameters must be derived such that a1 + a2 = 1, and therefore Eq. (16) is better expressed LDUAL = aLθ1 + (1 − a)Lθ2. The optimal relative source amplitude parameter, a, is found by optimisation of the pseudo-Z-statistic: Z¯ DUALopt ¼ maxðZ¯ a Þ; a

0VaV1

ð17Þ

Where Z¯a is given by substitution of Eq. (16) into Eqs. (13) and (5). The use of this amplitude optimisation allows for accurate spatial reconstruction of the two sources and derivation of the relative source strengths. This non-linear search adds an extra degree of freedom to those already represented by θ1 and θ2. The dual source beamformer, as described above, represents a ‘brute force’ approach to solving the correlated source problem and while attractive due to its simplicity, practically it is computationally intensive. Given enough computing power, the dual source beamformer would ideally allow for all combinations and permutations of both θ1 and θ2. This means that for each possible θ1, one computes the required parameter (i.e. pseudo-T-statistic or pseudo-Z-statistic), using optimised orientations, for all θ2 (and if required all a). This technique therefore represents a way of identifying and accurately locating two correlated sources anywhere in the source space with no a priori information. Furthermore, one might envisage the addition of further correlated sources if so required to explain the measured data. In practice, however, such an approach has a limited practical benefit due to the huge processing time required. However, as described below, in order to cut down on computational time, it is often possible to obtain the approximate position of one (or in some cases all) of the correlated sources either empirically, through the use of a single source beamformer, or by a second neuroimaging modality, for example fMRI. Methodology Simulations

used as a volume conductor model. The forward solution itself was based on the derivation by Sarvas (1987). White noise was added to the simulated p data ffiffiffiffiffiffi at each MEG sensor with an r.m.s. amplitude of 31.6 ft (10ft= Hz with a 10-Hz bandwidth (35–45 Hz)). Reconstruction of correlated sources Functional images of the two simulated correlated sources were made using a MATLAB implementation of both the single and dual source beamformer (μ = 0 in both cases). Covariance matrices were generated based on the simulated data spanning all time, and the 35- to 45-Hz frequency window. Within this time frequency window, the two sources were correlated with a correlation coefficient of 0.94. In order to reduce computational time, the known location and orientation of the first source was used as a seed, meaning that (for the purposes of this simulation) the lead fields Lθ1 remained fixed at the known position. Pseudo-Zstatistical images were then constructed as per the theory section in order to find the location, orientation and amplitude of the second source Qθ2. Having found Qθ2 using Qθ1 as a seed position, the two sources were swapped and the location and orientation of Qθ1 found using Qθ2 as a seed position. Effect of variation in source correlation on beamformer reconstruction A second set of simulations was undertaken in order to characterise the behaviour of both the single source, and dual source beamformer approaches at varying degrees of source correlation. The source configuration from the initial simulations was maintained, with two 5-nAm sources positioned approximately in left and right auditory cortex. The additive source noise was, however, reduced to zero in order to ensure maximum variability of the correlation coefficient between the two sources. Again, both sources were defined to include an initial 5 s of baseline activity, followed by 5 s of 40 Hz oscillation, and then a further 5 s of rest. The degree of correlation between the two sources was varied by changing the phase of the second source such that: Qθ1 ¼ Asinð2pftÞ and Qθ2 ¼ Asinð2pft þ vÞ

Simulations were used in order to assess the performance both of the single and dual source beamformers in reconstructing two spatially distinct, but correlated sources. For all simulations, the third order gradiometer configuration of a 275-channel whole head MEG system (CTF Systems Inc., Port Coquitlam, Canada) was used. The source space was based on a subject’s head shape measured from an anatomical MRI scan and its location based on a real experimental recording session. The two sources were positioned approximately in the left and right auditory cortices; however, for ease of visualisation, both sources were constrained to the same coronal plane. The orientation of each source was defined perpendicular to the radial direction. The sources were dipolar, of strength 5 nAm and both sources were defined to include an initial 5 s of baseline activity, followed by 5 s of 40 Hz oscillation, and then a further 5 s of rest. Additive uncorrelated source noise of amplitude 0.2 nAm was included. Ten trials were simulated at a sampling frequency of 600 Hz. For the solution to the forward problem, the best fitting sphere to the subject's head shape was

ð18Þ

ð19Þ

Where χ represents the phase difference between the sources and was varied between 0° and 90° in 10° steps. This resulted in 10 different source correlation coefficients in the range zero to one. As outlined in the theory section, the reconstructed source activity using a dual source beamformer should represent the sum of the two correlated sources. However, this relationship should also hold true even if the sources become uncorrelated. For this reason, the dual source beamformer reconstructed time course should be given by: Qθ1 þ Qθ2 ¼ A½sinð2pftÞ þ sinð2pft þ vÞ

ð20Þ

which reduces to:

 v v Qθ1 þ Qθ2 ¼ 2Asin 2pft þ cos 2 2

ð21Þ

For all ten values of χ, simulated data were analysed using both a single and dual source beamformer. Functional images and source time courses were reconstructed.

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

Effect of variation in relative source strengths In order to look into the effect of the dual source beamformer when the two sources are not equal in magnitude, a third simulation was undertaken. Here, the simulation parameters were as described above, however, while the amplitude of the source in the right hand auditory cortex remained at 5 nAm, the amplitude of the source in the left auditory cortex was increased to 9 nAm. This gave a simulated relative source amplitude parameter a = 0.643. Functional images and source time courses were reconstructed using the dual source beamformer both with and without amplitude optimisation (again μ = 0). In both cases, the source in the left hand auditory cortex was used as a seed position and the dual-source beamformer used to image the right hand source. A cubic ROI (5 cm side length) was used with the true source location at the centre. In the first case (i.e. without amplitude modulation), the parameter a remained fixed at 0.5. In the second case, for each image voxel, a was varied between 0 and 1. The value of the final pseudo-Z-statistical image was taken to be the optimised value of Z¯ as per Eq. (17). Auditory experiment The experiment was carried out on a 275-channel CTF system using the third order gradiometer configuration and a sampling rate of 600 Hz. The paradigm involved simple auditory stimulation and a single subject took part in the study in order to provide a proof of

1459

principle for the dual source beamformer technique. The subject listened passively to a 500-Hz pure tone that was amplitude modulated using a 40-Hz sinusoid. Each epoch was 3.8 s long and comprised 2 s of auditory stimulation followed by 1.8 s rest. This was repeated 75 times. The paradigm was known to induce bilateral steady state activity that was correlated across the left and right auditory cortices oscillating at 40 Hz (Herdman et al., 2003b; Ross et al., 2005; Simpson et al., 2005). Following data acquisition, a 3D digitiser (Polhemus Isotrack) was used to digitise the shape of the subject's head. A 150-Hz antialiasing filter and a 50-Hz powerline filter were applied to the MEG data that had also been DC corrected. Co-registration of the MEG data to MP-RAGE anatomical MR images (1 mm3 resolution taken on a Siemens Trio 3 T system) was then performed by matching (Fitzgibbon, 2001) the digitised head shape to the head surface extracted from the subjects MRI using custom made software. All beamformer analyses (both single source and dual source) were based on the formulation of pseudo-T-statistical images as per Eq. (6). In all cases, the active time window was defined as 1 s ≤ t ≤ 2 s with trial onset at t = 0, the passive window was defined as 2.5 s ≤ t ≤ 3.5 s. A frequency window of 38 Hz ≤ f ≤ 42 Hz was used to gain maximum sensitivity to the auditory steady state response. The data were initially analysed using a single source beamformer. Source estimates were made at each vertex of a regularly spaced grid spanning the entire source space with a spatial resolution of 5 mm.

Fig. 1. Source localisation and time course reconstruction of two correlated sources using a single (left) and dual (right) source beamformer implementation. The time course for the single source beamformer was estimated for the voxel where the original source was placed in the left hemisphere (though the reconstruction for the source in the right hand hemisphere is of a similar amplitude). The original source locations are (60 , 50) mm and (50, −50 mm). Note that the dual source beamformer reconstructs accurately the locations and summed amplitude of these two correlated sources.

1460

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

Fig. 2. Images showing the single source beamformer reconstruction (top) and dual source beamformer reconstruction (bottom) of one of the two simulated sources at varying degrees of correlation. Notice that whereas the single source beamformer fails to accurately reconstruct sources with a high degree of correlation, the dual source beamformer only fails when the two sources become orthogonal. The original source location is central to the field of view.

An unregularised dual source beamformer was applied to the recorded data. Analysis comprised two stages. Stage one

For this proof of principle study, the auditory cortex has been specifically chosen since the correlated sources lie in opposite hemispheres. For this reason, accurate reconstructions of correlated source power can be made using a single source beamformer by

Initially a constrained search for the optimal source positions and orientations of the two correlated sources was undertaken. For every θ1, a pseudo-T-statistic was computed for all possible locations and orientations θ2. However, in order to save computational time, potential locations were restricted to two 5 cm3 volumes spanning approximately the left and right auditory areas with a 0.5-cm resolution. This means that for each of the (113=) 1331 locations of the first source, we consider all of the potential 1331 locations of the second source giving (13312=) 1771561 possible combinations of the two source locations. For each combination, all possible orientations of both sources were considered as per Eq. (14). δ1 was varied from 0 to 180° in steps of 4°. δ2 was varied between 0° and 360° also in steps of 4°. This two-dimensional non-linear search resulted in optimised orientations for each potential pair of source locations, which was used in subsequent calculations. The pseudo-T-statistic was computed on each iteration of the algorithm. Note that for this example, the relative source amplitude parameter, a, was kept constant at 0.5 since sources were expected to be of approximately equal amplitude. The pair of locations for which the pseudo-T-statistic was a maximum was used as the optimised dipolar locations. Stage two Having found the optimum locations and orientations using step one, the seeded approach to dual source beamforming was adopted in order to produce an image of the two correlated sources. In the same way as for the simulations, the location of source one, as defined by stage one, was used as the seed location. A dual source beamformer image was created by allowing all possible locations of source two (this time no constraint was placed on the possible locations of source two). Having obtained this first image, a second image was obtained by switching the sources such that source two became the seed and source one was allowed to roam. The final correlated source image was then constructed by combining (summing) the two individual images. Relative source amplitudes were also checked using the technique outlined by Eq. (17).

Fig. 3. The angular distribution of the pseudo-Z-statistic at all of the simulated phase differences. For sources with a high degree of correlation, a single peak dominates and is an accurate representation of the simulated source orientation. As correlation is reduced, a second peak is observed 180° out of phase with the first. This is representative of the source orthogonalisation at low correlation coefficients.

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

basing the beamformer solution on only half of the available channels (i.e. MEG sensors in either the right or left hemisphere are switched off). Obviously, in general this is not a viable approach (for correlated activity in primary and lateral areas of the visual cortex for example), but here it allows for a simple confirmation of results. For this reason, two further single source beamformer images were created, the first using only channels in the left hemisphere, the second using only channels in the right hemisphere. These images were used as a standard against which to check the dual source localisation accuracy. In addition, further verification of the accuracy of the dual source beamformer technique was achieved by use of a 2-dipole fit to the averaged data. A single spatiotemporal dipole fit was performed (latency 1– 2 s) using a 2-dipole model and standard commercial software (CTF Systems Inc.). Errors in the fitted parameters were assessed using a Monte Carlo technique (Medvick et al., 1989). Results Fig. 1 shows the results for the simulation with 2 strongly correlated sources. Source localisation and time course reconstruction using a single source beamformer are shown on the left, source localisation and time course reconstruction using the dual source beamformer are on the right. Note that for the single source beamformer almost all of the power is suppressed in the

1461

final image, whereas in the dual source approach, the power is maintained. Also note that in the time course reconstruction, using a single source beamformer the reconstructed amplitude is much lower than the simulated level (5 nAm), whereas using the dual source beamformer, the full amplitude of each source has been summed giving a combined amplitude of 10 nAm. In this case, an accurate temporal reconstruction of source activity for each individual source can be derived by simply dividing the reconstructed time course for the dual source by two. Fig. 2 shows both the single source and dual source reconstruction of one of the two simulated sources at varying degrees of correlation (relative phase). For low phase difference (high correlation), the single source beamformer fails to reconstruct accurately the source, whereas the method works well as the sources become temporally orthogonal. In contrast, the dual source beamformer works well in the case where a high correlation exists between the two sources; however, spatial resolution is lost when the two sources become orthogonal. The reason for the loss in spatial resolution at low correlation using a dual source beamformer is (at least in part) explained in Fig. 3. Here, the angular distribution of the pseudo-Z-statistic is shown for all of the simulated phase differences. For sources with a high degree of correlation, a single peak dominates the distribution and is an accurate representation of the simulated source orientation. As correlation is reduced, a second peak is observed 180° out of phase

Fig. 4. Simulated source time courses reconstructed using the single and dual source beamformers. (A) Crossed data points show the amplitude of the dual source reconstruction as a function of phase change. Circled data points show the amplitude of the single source reconstruction as a function of phase change. The dual source reconstructed signal should represent the sum of the two individual source time courses, and the theoretically derived amplitude (Eq. (21)) is shown by the dashed line. The single source reconstruction shows how signal amplitude is suppressed with increasing correlation between sources. (B) Raw signals reconstructed using the dual source beamformer for all 10 phase differences. (C) Raw signals reconstructed using the single source beamformer for all 10 phase differences.

1462

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

with the first. This second peak shows that at low correlation coefficients, the dual source beamformer cannot accurately derive the orientations of the two simulated sources. Fig. 4A shows the single and dual source beamformer reconstructed source amplitudes as a function of the phase angle between the simulated sources. These amplitudes were measured in the time window 6 s ≤ t ≤ 9 s, where t = 0 is defined as the start of each trial. The dual source reconstructed time course accurately represents the sum of the two individual source time courses. The single source reconstruction shows how signal amplitude is suppressed with increasing correlation between sources. In support of this, Fig. 4B shows the reconstructed signals using the dual source beamformer in the time window 6.00 s ≤ t ≤ 6.05 s, again for all 10 phase differences. These signals are in agreement with those theoretically derived in Eq. (21). In contrast, Fig. 4C shows the raw signals recons-

tructed using the single source beamformer for all 10 phase differences. Here we observe a loss of amplitude with increasing correlation. Fig. 5 shows the results of the amplitude optimisation. Fig. 5A shows a dual source beamformer image derived both with and without amplitude optimisation for two non-equal sources. As shown, the amplitude optimisation corrects for the spatial distortion caused by the sources being of different amplitude. Notice also that the pseudo-Z-statistic is increased from 4.5 to 13 in the case where amplitude optimisation was used. Fig. 5B shows the pseudo-Z-statistic as a function of a, Notice the sharpness of the optimised peak indicating that amplitude optimisation can identify accurately the relative source amplitude parameter. The figure also shows that the derived value for a was close to the simulated value of 0.643. Fig. 5C shows the reconstructed time course. Once again, the dual source beamformer has accurately reconstructed the sum of the two individual time courses. Fig. 6A shows the experimental results from the auditory data. We see that whereas the dual source beamformer allows accurate reconstruction of the two correlated sources in the left and right auditory areas, source power remains suppressed in the single source beamformer reconstruction. The peak locations as obtained with the dual source beamformer were compared to those derived using a single source beamformer (where half of the MEG sensors were eliminated in order to prevent the correlated source problem). This revealed a discrepancy of 1 image voxel in the case of the left hand peak, and zero in the case of the right hand peak. Fig. 6B shows the spatial relationship between the peaks in the dual source beamformer image, and a 2-dipole fit. In the case of the right hand peak, the discrepancy between the peak maxima and dipole is 4 mm. In the case of the left hand peak, the discrepancy between the peak maxima and dipole is 4.5 mm. In both cases the discrepancy is less than the voxel dimension in the beamformer image. The Monte Carlo derived error estimate in the dipole fits were 2 mm in the case of the left hand dipole and 4 mm in the case of the right hand dipole. In order to show the angular sensitivity of the technique, the result from the 2-dimensional non-linear search for optimum source orientations (Eq. (14)) is shown in Fig. 7 for the peak locations in the dual source image of auditory cortex. As shown, the technique estimates source orientation with reasonable angular resolution. Finally, amplitude optimisation showed that the two sources were approximately equal in magnitude (result not shown). Discussion

Fig. 5. Results of the amplitude optimisation. (A) The dual source beamformer image derived both with and without amplitude optimisation. Note that the pseudo-Z-statistic is increased from 4.5 to 13 in the case where amplitude optimisation has been used and that the peak value is in the correct location (green dot) in this case. (B) The measured pseudo-Z-statistic as a function of a. Note that the derived value of a is close to the simulated value of 0.643. (C) The reconstructed time course showing the sum of the two individual source time courses.

This paper has shown that while the standard, single source, beamformer suppresses the reconstructed power of two spatially separate but temporally correlated sources, a simple modification to the source model can allow for accurate spatial reconstruction of those two sources. We have shown that this technique can be used to locate temporally correlated sources, and further that in all cases, the reconstructed time course is an accurate representation of the summed activity in those areas. Also, for strongly correlated sources, an accurate temporal representation of the source activity in each individual source can be obtained using amplitude optimisation. The technique has been validated using simulations and experimental data and has shown that accurate reconstruction of correlated steady state auditory activity is achievable. The main limitation of the dual source beamformer technique is that the reconstructed time course represents the summed source

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

1463

Fig. 6. (A) Experimental results showing the localisation of correlated sources in left and right auditory cortex. The dual source reconstruction is shown on the left and is thresholded at T ¯ = 1. The single source reconstruction is shown on the right and is not thresholded. Again note that for two temporally correlated sources, the dual source approach is successful whereas the standard beamformer fails. (B) The spatial relationship between the peaks in the dual source beamformer image (peak voxel shown in red) and a 2-dipole fit (shown as the yellow and green circles).

activity of both dipoles. We have shown that, in the case where two dipoles are perfectly correlated, the reconstructed time course is an accurate representation of the two original time courses and the individual source time courses can be obtained by amplitude optimisation. However, in the case where only a partial correlation

Fig. 7. The result of the non-linear search of optimum orientations of two correlated sources. In this case, the search was performed for the optimum locations of the left and right auditory sources.

exists, in some cases the correct amplitudes of the individual sources cannot be regained, and the temporal signature will be incorrect. This is apparent from Fig. 4, in which we notice that for non-zero phase differences, although the reconstructed time course using the dual source beamformer is an accurate reflection of the summed activity, it is not reflective of the time course from either one of the two original sources. A previously described technique (Sekihara et al., 2002b) has, however, been shown to accurately reconstruct the source amplitudes of partially correlated sources. We therefore propose that, in such cases, the two methods be used in conjunction. The dual source technique can be used to locate temporally correlated sources, and the technique of Sekihara et al. (2002b) can be used to accurately reconstruct time courses in the case of a partial correlation. (Also since the dual source beamformer does accurately reconstruct summed activity, this might be used as a validation of the technique described in Sekihara et al., 2002b.) The dual source beamformer approach is computationally demanding. Given a reasonable set of starting parameters (i.e. an anatomically derived constraint), it is possible to localise accurately two correlated sources by considering all possible combinations and permutations of source orientation and location in a reasonable time frame (approximately 12 h using our MATLAB implementation on a high-spec PC). However, while extension of the search areas, or an increase in the number of correlated sources is possible, these would extend computational time considerably. Techniques exist by which to speed up the process, for example moving from a non-linear to a linear beamformer would remove the requirement for the search for optimum orientations. In this case the spatial resolution of the technique is likely to be reduced (Robinson and Vrba, 1998), but it would be a useful technique by which to obtain an initial guess of the source locations. In a recently described technique (Sekihara et al., 2004), the optimum orientation is estimated from a singular value decomposition of the reconstructed SNR matrix

1464

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465

for 3 orthogonal source orientations. Application of this methodology may speed up the process since it also removes the requirement for the non-linear angular optimisation. Symmetrically constrained dual source models are also routinely used in the BESA software, and this too cuts down the computational burden of dual source modelling. However, this is reliant on the anatomical structure of the brain being perfectly symmetric and so should not be considered an optimal solution to the problem. A beamformer solution is obtained independently for each voxel, or combination of voxels in case of the dual-source beamformer; hence, distributed computing may also be used in order to speed up the computation. Finally, the beamformer calculation itself is independent across voxels (i.e. the calculation at one voxel is not related to that at any other). This makes the algorithm well suited to parallelisation. For this reason, the use of large computer clusters would also be extremely advantageous in reducing computational time. The dual source beamformer represents a valuable addition to beamformer analysis if there are suspicions of correlated sources in the recorded MEG data. For two spatially distinct sources, if the calculated pseudo-T-statistic (or equivalently pseudo-Z-statistic) is higher using a dual source model than it is using a single source model at either of the two source locations, then we conclude that correlations exist between those two sources. Here, we have used the example of auditory steady state activity to prove this, and we show that two correlated sources can be imaged accurately using the new technique. Other potential applications can be conceived, for example, in the case of a visual experiment involving moving stimuli, one might expect correlated activity in both primary and lateral (V5/MT) visual areas. The dual source beamformer could be used to locate this activity whereas a single source beamformer would be of limited use. The technique itself is speeded up considerably in the case where the location of one of the correlated sources can be estimated a priori. This might be achieved based on anatomical information, MEG/EEG results using similar paradigms, or results from other functional imaging techniques. Again using the example of a visual stimulus, a seed position could be constructed for primary visual cortex activity by using a static visual stimulus and reconstruction of the response to the static stimulus, or by performing an equivalent fMRI experiment. Given a situation in which the initial guess at a seed location is reasonably accurate, the dual source beamformer will then efficiently locate the second correlated source and reconstruct the summed time course activity. In some cases, that initial guess could then be modified using several iterations of the dual source algorithm. Conclusion In conclusion, the dual source beamformer represents a way to image two correlated sources using beamformer methodology. We have shown that while the standard, single source beamformer suppresses the reconstructed power of two spatially separate but temporally correlated sources, the dual source beamformer allows for accurate reconstruction. The technique has been proven in simulations and using real experimental data. The technique is computationally demanding, although the use of large computer clusters and use of a priori information may shorten computational time. Finally, we suggest that this technique represents a useful way of locating correlated sources, particularly if a seed location can be defined a priori.

References Barnes, G.R., Hillebrand, A., 2003. Statistical flattening of beamformer images. Hum. Brain Mapp. 18, 1–12. Brookes, M.J., Gibson, A.M., Hall, S.D., Furlong, P.L., Barnes, G.R., Hillebrand, A., Singh, K.D., Holliday, I.E., Francis, S.T., Morris, P.G., 2005. GLM-beamformer method demonstrates stationary field, alpha ERD and gamma ERS co-localisation with fMRI BOLD response in visual cortex. NeuroImage 26, 302–308. Fitzgibbon, A.W., 2001. Robust Registration of 2D and 3D Point 662–670. Gross, J., Timmermann, L., Kujala, J., Salmelin, R., Schnitzler, A., 2003. Properties of MEG tomographic maps obtained with spatial filtering. NeuroImage 19, 1329–1336. Hadamard, J., 1952. Lectures On Cauchy's Problem In Linear Partial Differential Equations. Hadjipapas, A., Hillebrand, A., Holliday, I.E., Singh, K.D., Barnes, G.R., 2005. Assessing intersections of linear and non-linear neuronal sources using MEG beamformers: a proof of concept. Clin. Neurophysiol. 116, 1300–1313. Hamalainen, M.S., Ilmoniemi, R.J., 1984. Interpreting measured magnetic fields of the brain: estimates of current distributions. (Helsinki University of Technology, Finland) Technical Report TKKK-F-A599. Hamalainen, M., Hari, R., Ilmoniemi, R.J., Knuutila, J., Lounasmaa, O.V., 1993. Magnetoencephalography—Theory, instrumentation and applications to non-invasive studies of the working human brain. Rev. Modern Phys. 21 (2), 413–460. Herdman, A.T., Wollbrink, A., Chau, W., Ishii, R., Ross, B., Pantev, C., 2003a. Determination of activation areas in the human auditory cortex by means of synthetic aperture magnetometry. NeuroImage 20 (2), 995–1005. Herdman, A.T., Wollbrink, A., Chau, W., Ishii, R., Ross, B., Pantev, C., 2003b. Determination of activation areas in the human auditory cortex by means of synthetic aperture magnetometry. NeuroImage 20, 995–1005. Hillebrand, A., Barnes, G.R., 2005. Beamformer analysis of MEG data. Int. Rev. Neurobiol. 68, 149–171 (Special volume on Magnetoencephalography). Hillebrand, A., Brookes, M.J., Singh, K.D., Furlong, P.L., Barnes, G.R., 2005a. Spatial extent of neuronal activity measured with MEG beamformers. Proceedings of the 11th Annual Meeting of the Organisation for Human Brain Mapping, Toronto. Hillebrand, A., Singh, K.D., Holliday, I.E., Furlong, P.L., Barnes, G.R., 2005b. A new approach to neuroimaging with magnetoencephalography. Hum. Brain Mapp. 25, 199–211. Huang, M.-X., Shih, J.J., Lee, R.R., Harrington, D.L., Thoma, R.J., Weisend, M.P., Hanlon, F., Paulson, K.M., Martin, K., Miller, G.A., Canive, J.M., 2004. Commonalities and differences among vectorized beamformers in electromagnetic source mapping. Brain Topogr. 16 (3), 139–158. Kuwada, S., Batra, R., Maher, V.L., 1986. Scalp potentials of normal and hearing-impaired subjects in response to sinusoidally amplitudemodulated tones. Hear. Res. 21 (2), 179–192. Medvick, P.A., Lewis, P.S., Aine, C., Flynn, E.R., 1989. Monte Carlo analysis of localization errors in magnetoencephalography. In: Williamson, S.J., Hoke, M., Stroink, G., Kotani, M. (Eds.), Advances in Biomagnetism. Plenum Press, New York, pp. 543–546. Mosher, J.C., Leahy, R.M., 1998. Recursive MUSIC: a framework for EEG and MEG source localisation. IEEE Trans. Biomed. Eng. 45 (11), 1342–1354. Pascal-Marqui, R.D., Michel, C.M., Lehmann, D., 1994. Low resolution electromagnetic tomography, a new method for localizing electrical activity in the brain. Int. J. Psychophysiol. 18, 49–65. Preissl, H., 2005. Magnetoencephalography. Acedemic Press. Rees, A., Green, G., G.R., Kay, R.H., 1986. Steady-state evoked responses to sinusoidally amplitude-modulated sounds recorded in man. Hear. Res. 23 (2), 123–133. Robinson, S., Vrba, J., 1998. Functional neuroimaging by synthetic aperture magnetometry. In: Yoshimoto, T., Kotani, M., Kuriki, S., Karibe, H.,

M.J. Brookes et al. / NeuroImage 34 (2007) 1454–1465 Nakasato, N. (Eds.), Recent Advances in Biomagnetism. Tohoku Univ. Press, Sendai, pp. 302–305. Ross, B., Borgmann, C., Draganova, R., 2000. A high-precision magnetoencephalography study of human auditory steady-state response to amplitude-modulated tones. J. Acoust. Soc. Am. 108 (2), 679–691. Ross, B., Herdman, A.T., Pantev, C., 2005. Right hemispheric laterality of human 40 Hz auditory steady-state responses. Cereb. Cortex 15 (12), 2029–2039. Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol. 32 (1), 11–22. Sekihara, K., Nagarajan, S.S., Peoppel, D., Marantz, A., Miyashita, Y., 2001. Reconstructing spatio-temporal activities of neural sources using an MEG vector beamformer technique. IEEE Trans. Biomed. Eng. 48 (7), 760–771. Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., 2002a. 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., Miyashita, Y.M., 2002b. Application of an MEG eigenspace beamformer to reconstructing spatiotemporal activities of neural sources. Hum. Brain Mapp. 15, 199–215. Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., 2004. Asymptotic

1465

SNR of scalar and vector minimum-variance beamformers for neuromagnetic source reconstruction. IEEE Trans. Biomed. Eng. 51, 1726–1734. Simpson, M.I.G., Hadjipapas, A., Barnes, G.R., Furlong, P.L., Witton, C., 2005. Imaging the dynamics of the auditory steady-state evoked response. Neurosci. Lett. 16 (3), 195–197. Singh, K.D., Barnes, G.R., Hillebrand, A., Forde, E.M., Williams, A.L., 2002. Task related changes in cortical synchrony are spatially coincident with the haemodynamic response. NeuroImage 16, 103–114. Van Drongelen, W., Yuchtman, M., Van Veen, B.D., Van Huffelen, A.C., 1996. A spatial filtering technique to detect and localize multiple sources in the brain. Brain Topogr. 9 (1), 39–49. Van Veen, B.D., Van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localisation of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9). von Helmholtz, H., 1853. U¨ber einige Gesetze der Vertheilung elektrischer Stro¨me in ko¨rperlichen Leitern, mit Anwendung auf die thierischelektrischen Versuche. Ann. Phys. Chem. 89, 353–377. Vrba, J., 2002. Magnetoencephalography: the art of finding a needle in a haystack. Physica C 368, 1–9. Vrba, J., Robinson, S.E., 2001. Signal processing in Magnetoencephalography. Methods 25, 249–271.