www.elsevier.com/locate/ynimg
NeuroImage 27 (2005) 368 – 376

A simple nonparametric statistical thresholding for MEG spatial-filter
source reconstruction images
Kensuke Sekihara,a,* Maneesh Sahani,b and Srikantan S. Nagarajanc
aDepartment of Systems Design and Engineering, Tokyo Metropolitan University, Asahigaoka 6-6, Hino, Tokyo 191-0065, Japan
Gatsby Computational Neuroscience Unit, University College London, 17 Queen Square, London, WClN 3AR, UKc
Department of Radiology, University of California, San Francisco, 513 Parnassus Avenue, S362 San Francisco, CA 94143, USAb

Received 12 October 2004; revised 3 February 2005; accepted 5 April 2005
Available online 29 June 2005

This paper proposes a simple statistical method for extracting target
source activities from spatio-temporal source activities reconstructedfrom MEG measurements. The method requires measurements in acontrol condition, which contains only non-target source activities. The
method derives, at each pixel location, an empirical probability
distribution of the non-target source activity using the time coursereconstruction obtained from the control period. The statisticalthreshold that can extract the target source activities is derived fromthe empirical distributions obtained from all pixel locations. Here, themultiple comparison problem is addressed with a two-step procedure
involving standardizing these empirical distributions and deriving anempirical distribution of the maximum pseudo T value at each pixel
location. The results of applying the proposed method to auditoryevoked measurements are presented to demonstrate the method’s
effectiveness.
D 2005 Elsevier Inc. All rights reserved.
Keywords: Magnetoencephalography; Statistical threshold; Biomagnetism;
Nonparametric statistics; Spatial filter; Spatio-temporal reconstruction

Introduction
Among the various technologies for noninvasive neural
measurement, the major advantage of magnetoencephalography(MEG) is its ability to provide fine temporal resolution on the orderof milliseconds (Hämäläinen et al., 1993). Neuromagnetic imaging
can therefore visualize neural activities with such a fine time
resolution and provide functional information about brain dynamics. One major problem is that the measured MEG signal
generally contains not only a magnetic field associated with thesignal sources of interest but also contains interference magnetic
fields generated from non-target activities. Such non-target
activities include spontaneous brain activities or some evoked* Corresponding author. Fax: +81 42 585 8642.
E-mail address: ksekiha@cc.tmit.ac.jp (K. Sekihara).
Available online on ScienceDirect (www.sciencedirect.com).
1053-8119/$ - see front matter D 2005 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2005.04.009

activities that are not the interest of the current investigation. These
non-target activities generally overlap with the target signal
activities in the source reconstruction, and they often makeinterpreting the reconstruction results difficult.
In most studies using positron emission tomography (PET) or
functional magnetic resonance imaging (fMRI), the experiments
are carefully designed to extract only target activities. A common
example of such experimental design contains two kinds of stimuli:a task stimulus and a control stimulus. The task stimulus generally
elicits the target cortical activities as well as other activities
associated with the target activities. The control stimulus isdesigned to elicit only the latter activities. Then, by calculating
the statistical difference between the images measured with the twokinds of stimuli, the target activities can be revealed. Bothparametric (Friston et al., 1995; Worsley et al., 1996) andnonparametric statistics (Nichols and Holmes, 2001) have beenused to calculate such statistical differences.
This paper proposes a simple method of statistical subtraction
between task and control measurements. The method is applicableto spatio-temporal source reconstruction from MEG/EEG measurements. It assumes neural activities to be quasi-stochastic, and it
uses nonparametric statistics to derive an empirical probability
distribution of these activities using the time course reconstructionin the control period. This empirical distribution is then used forderiving an appropriate value for the statistical thresholding. The
thresholding can extract the target activities that exist only in thetask measurements by eliminating other non-target activities thatexist both in the task and control measurements.
In this paper, we present the proposed statistical thresholding
method using spatial filter source reconstruction (Sekihara andNagarajan, 2004). This is because the formulation of the spatial
filter is relatively simple and the spatial filter techniques have beensuccessfully applied to MEG source analysis (Hashimoto et al.,
2003; Ishii et al., 2003). However, the applicability of the proposed
method is not limited to the spatial filter formulation and it can beused with any type of source estimation method that can provide
the spatio-temporal source reconstruction, i.e., that can reconstructsource time courses at all pixel locations.

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

Spatial filter formulation for MEG source reconstruction
Definitions
We define the magnetic field measured by the mth detector coilat time t as b(t) and a column vector b(t) = [b 1(t), b 2(t),. . ., b M (t)]T
as a set of measured data where M is the total number of sensor coils
and superscript T indicates the matrix transpose. The spatia1
location is represented by a three-dimensional vector r: r = (x, y, z).
The second-order moment matrix of the measurement is denoted
R, i.e., R = bb(t)b T (t) where b  indicates the ensemble average,
which is replaced by the time average over a certain time window
in practice. (When bb(t) å 0 holds, R is also equal to thecovariance matrix of the measurement.) The magnitude of thesource moment is denoted s(r, t). The source orientation isdefined as a three-dimensional column vector h(r) = [g x (r), g y (r)g z (r)]T whose f component (where f equals x, y, or z) is equal tothe cosine of the angle between the direction of the source and thef axis.
f
We define l m(r) as the output of the mth sensor; the output is
induced by the unit-magnitude source located at r and pointing in
the f direction. The column vector l f (r) is defined as l f (r) = [l 1f(r), l 2f (r), . . ., l Mf (r)]T. We define the lead field matrix, which
represents the sensitivity of the whole sensor array at r, as L(r) =[l x (r), l y (r), l z (r)].

I

Typical non-adaptive and adaptive spatial filters
Spatial filter techniques estimate the source current density by
applying a simple linear operation to the measured data, i.e.,Mŝs ðr; t Þ ¼ wT ðrÞbðtÞ ¼~ wm ðrÞbm ðtÞ;

ð1Þ

m ¼ 1

where ŝ(r, t) is the estimated source magnitude. The column vector369

source distributions (Baillet et al., 2001; Dale et al., 2000; PascualMarqui, 2002).
The other type of spatial filter is an adaptive spatial filter inwhich the filter weight depends on the measurements. The bestknown adaptive spatial filter is the minimum variance spatial filter,
which is customarily referred to as the minimum-variance beamformer (Robinson and Vrba, 1999; Sekihara and Seholz, 1996;
Sekihara et al., 2001; van Veen et al., 1997). In this method, theweight vector is given byl T ðrÞR1
w T ðr Þ ¼ T;
l ðrÞR1 l ðrÞ

ð5Þ

where l(r) is defined as l(r) = L(r)h opt(r). (The method for
determining h opt(r) for the adaptive spatial filter is described in(Sekihara and Seholz, 1996). Furthermore, in practice, thenormalized lead field l(r)/||l(r)|| is used in Eq. (5) to avoid asource location bias caused by the variation of the lead field norm||l(r)|| (Sekihara et al., 2005). The minimum-variance beamformer
can be extended to the eigenspace-projection beamformer, which is
known to be tolerant of errors in the forward modeling or in theestimation of the data covariance matrix (Sekihara et al., 2002).
The extension is attained by projecting the weight vector in Eq. (5)onto the signal subspace of the measurement covariance matrix.
That is, redefining the weight vector in Eq. (5) as w (MV) (r), theweight vector for the eigenspace-projection beamformer isobtained usingwðrÞ ¼ ES ETS wðMVÞ ðrÞ:ð6ÞIn this equation, ES is a matrix whose columns consist of thesignal-level eigenvectors of R, and ES EST is the projection matrix
that projects a vector onto the signal subspace of R. Thiseigenspace projection beamformer was used in the experiments
described in Numerical experiments and Experiments.

wðrÞ ¼ ½w1 ðrÞ; . . . . . . ;wM ðrÞT
represents a set of the filter weights. The filter weights w(r) should
only pass the signal from a source with a location r and reject thesignals generated at other locations. Since the source is a threedimensional vector quantity, the weight vector generally depends
on the source directions. In this paper, we define w(r) as the weightvector in the optimum direction h opt(r), which is determined as the
direction that gives the maximum spatial-filter outputs at each r.
There are two types of spatial filter techniques. One is a nonadaptive method in which the filter weight is independent of themeasurements. The best-known non-adaptive spatial filter is theminimum-norm estimate (Hämäläinen and Ilmoniemi, 1984). The
filter weight is expressed aswðrÞ ¼ G 1 LðrÞhopt ðrÞ:ð2ÞThe matrix G is often referred to as the gram matrix, which is given
by calculating the overlap between the lead fields,
G ¼ XLðrÞLT ðrÞdr:ð3Þ

The estimated current density is then expressed asŝs ðr; t Þ ¼ hTopt ðrÞLT ðrÞG1 bðt Þ:

ð4Þ

The minimum-norm spatial filter in its original form is known togive erroneous source reconstruction results. Therefore, to improve
the performance, it is often used with some kind of constraint on

Evaluation of statistical significance using parametric statistics
The evaluation of the statistical significance of the spatial filter
outputs has typically been performed using parametric statistics(Barnes and Hillebrand, 2003; Dale et al., 2000; Gross et al., 2001;
Robinson and Vrba, 1999). The basic assumption of the parametricmethod is that the measurement consists of deterministic signal andGaussian noise, i.e.,bðt Þ ¼ bI ðtÞ þ nðtÞ;ð7Þ

where b I (t) is the signal of interest, i.e., the signal generated frombrain sources that are the target of current investigation.
In Eq. (7), n(t) is the noise vector, and each element of n(t) is
assumed to follow N 0;r20 , which indicates the Gaussian
distribution with zero mean and a variance of r 02. The spatial
filter outputs ŝ(r, t) are expressed asŝs ðr; tÞ ¼ wT ðrÞbðt Þ ¼ wT ðrÞbI ðt Þ þ wT ðrÞnðt Þ:

ð8Þ

Therefore, since these Gaussian processes are assumed to be
uncorrelated between
different sensor recordings,
the outputs ŝ(r, t)follows N wT ðrÞ bI ðt Þ;r20 jjwðrÞjj2 ; which is a Gaussian
distribution with a mean of w T (r)bI (t) and a variance ofr 02||w(r)||2. Actually, since r 02 must be estimated from the measured
data, the distribution of ŝ(r, t) is not exactly represented by theGaussian distribution but by the t distribution.

370

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

The statistical evaluation can be performed by testing the null
hypothesis at each pixel location. (The null hypothesis is that there
is no signal source activity.) That is, the z score under the null
hypothesis, ŝ(r, t)/(r 0||w(r)||), is calculated and compared to z a/2,
which is the two-tailed z score corresponding to the a level ofsignificance, which is equal to the probability of Type I error. This
procedure is performed at each pixel location, and if the calculated
z score is higher than z a/2, the estimated source activity ŝ(r, t) is
considered to be statistically significant. This procedure can beextended to incorporate the multiple comparison problem (Barnes
and Hillebrand, 2003).

Proposed nonparametric statistical significance evaluation
The signal and noise model expressed in Eq. (7) is, in general,
insufficient to express real-world measurements, and the measured
data should be expressed asbðt Þ ¼ bI ðt Þ þ bn ðt Þ þ nðtÞ;ð9Þ

where b n (t) is the magnetic field generated from sources other than
the signal sources, such as spontaneous brain activities or someevoked activities that are not the target of the current investigation.
This b n (t) is often referred to as the brain noise. Here, we propose a
simple nonparametric method that can take such brain noise into
consideration.
The spatial filter outputs obtained from b(t) is expressed as

and then by choosing |ŝc (r, t ( q))| as R(r) where q = [(1 a)K c ] and |[ ] indicates the maximum integer that does notexceed the value in parenthesis. However, the statisticalthreshold obtained in this manner does not take multiple
comparisons into consideration, and instead of implementing
the above mentioned procedure, the following procedure is
performed.
(iii) Statistical thresholding with multiple comparisons. The
proposed method uses maximum statistics (Blair andKarniski, 1994; Pantazis et al., 2003) to address themultiple comparison problem. To utilize maximum statistics, we first standardize the empirical distribution of|ŝc (r, t j )| by calculating T(r, t j ) such that

I

 
 
  jŝs c r; t j j  bjŝs c r; tj jc
T r; t j ¼
r ðr Þ
r̂

ð13Þ

Here,

ŝs ðr;t Þ ¼ wT ðrÞbI ðt Þ þ wT ðrÞbn ðtÞ þ wT ðrÞnðtÞ
¼ ŝs I ðr; t Þ þ ŝs c ðr;t Þ;

is obtained at each pixel location r, F̂(x) is rewritten as
F̂(x|r) below.
(ii) Statistical thresholding without multiple comparisons. UsingF̂(x|r), we could obtain the statistical threshold at r, R(r),such that R(r) = F̂ 1(1  a|r) where a is a level of significance. In practice, the inverse of the empirical distribution
can be calculated by first sorting |ŝc (r, t j )| in increasing order:






jŝs c r; tð2Þ j > jŝs c r; tðKc Þ j;ð12Þjŝs c r; t ð1Þ j

ð10Þ

 2
  2
r 2 ðrÞ ¼ bŝs c r; tj c  bjŝs c r; t j jc ;
r̂

I

where
T

T

ŝs I ðr;tÞ ¼ w ðrÞbI ðtÞ and ŝs c ðr;t Þ ¼ w ðrÞðbn ðt Þ þ nðt ÞÞ:ð11ÞHere, ŝ I (r, t) is the estimated source activity of interest, and ŝ c (r,t) is the estimated background interference plus noise. The problem
with the parametric modeling described in the preceding section isthat it cannot efficiently take the background interference into
account because the Gaussianity assumption may not hold for ŝ c (r,t). The key assumption in the proposed method is that the control
measurement that can provide b c (t) = b n (t) + n(t) is available. In
general, this assumption is approximately fulfilled because the prestimulus measurement can be considered as a control in many cases.
Using this control measurement, the proposed method first derives
an empirical distribution of ŝ c (r, t), and with this empirical
distribution, the method determines the statistical threshold. The
procedures are as follows:(i) Empirical distribution formation. We calculate ŝ c (r, t j ) byapplying the spatial filter to the control measurement b c (t j )where t j is the discrete time point in the control measurement. We then calculate F̂(x), which is the empirical
distribution of the modulus of the time course |ŝ c (r, t j )|,such that F̂(x) = #{|ŝ c (r, t j )| x}/K c where #{|ŝ c (r, t j )| < x}indicates the number of |ŝ c (r, t j )| which is less than orequal to x, and K c is the number of total time points inthe control measurement. The implicit assumption here isthat the time course reconstruction |ŝ c (r, t j )| is identicallyand independently distributed (IID) over the pre-stimulus
time period. This procedure is repeated, and the empirical
distribution is calculated at all pixel locations. Since F̂(x)and b c indicates the time average over the control
period, i.e.,
Kc
 2
 2
1
bŝs c r; tj c ¼
ŝs c r; t j and~Kc j ¼ 1

 
 
1 Kc
~ jŝs c r; tj j:
bjŝs c r; t j jc ¼
Kc j ¼ 1
We then calculate the maximum T value T max (r) at each pixel
location. The maximum T value at the ith pixel location is denotedi
T max, where i = 1,. . ., K N , and K N indicates the total number ofipixels. We next obtain the empirical distribution of T max, Ĥ (x), suchiithat Ĥ (x) = #{T max x}/K N , where #{T max x} is the number ofi
T maxvalues which is less than or equal to x. We can then obtain theiththreshold of the T max
value for the a-significance level, T max, suchththat T max= Ĥ (1  a). The inverse of this empirical distribution canibe calculated by first sorting T max
in increasing order:
ð1Þ
ð2Þ
ðKN Þ
Tmax
< Tmax< > < Tmax
;

ð14Þ

p)
th
as T max
where p = [(1  a)K N ]. We finally obtainand choose T (maxthe statistical threshold for the spatial-filter reconstruction, R(r), bythconverting T max
into the source activity value, that is, thr ðrÞ þ bjŝs c r;tj jc :r̂~ ðrÞ ¼ Tmax
ð15ÞWe evaluate the statistical significance of the spatial filter
outputs by comparing the outputs |ŝ(r, t)| with ~ (r), and when|ŝ(r, t)| > ~ (r), the outputs ŝ(r, t) are considered to be statisticallysignificant. It should be mentioned that, once we take T max (r) at

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

371

each pixel location, the distributional properties of the process do
not matter anymore, and the assumption that the time coursereconstruction is IID is not needed. Instead, this multiple
comparison procedure imposes a new assumption that the shapeof the control source distribution is the same at each pixel, although
its scale may be different.

Numerical experiments
We conducted numerical experiments to show the effectivenessof the proposed statistical thresholding. We use a sensor alignmentof the 37-sensor array from the Magnesi (4D Neuroimaging Inc.,San Diego) neuromagnetometer. The source-sensor configurationand the coordinate system are illustrated in Fig. 1(a). The
coordinate origin is set at the center of the sensor coil located at
the center of the array. The three point-like sources, shown by thesmall filled circles in this figure, are assumed to be located at (0,1, 6) cm, (0, l, 6) cm, and (0, 1.6, 7.2) cm on the same plane(x = 0). These locations of the three sources are denoted r l, r 2, andr 3. The simulated magnetic field is calculated for 400 ms prestimulus and 400 ms post-stimulus time windows with a sampling
rate of 1 kHz. Here, the nearly orthogonal three time courses
shown in Fig. 1(b) are used as the time courses of the three sources.
In these experiments, the pre-stimulus period is considered thecontrol period. The first and second sources are considered the signal
sources of interest because they are only active in the post-stimulus(task) period. The third source is considered the control source
because it is active both in the pre- and post-stimulus periods. To
simulate brain background activity (which is so-called brain noise),spontaneous MEG measured from an awake human subject wasadded to this computer-generated magnetic field to create simulated
magnetic recordings. Here, the spontaneous MEG was measured
with a sampling rate of 1 kHz using the same sensor array andaveraged over 400 trials. These spontaneous MEG data are shown in
the fourth panel (from the top) of Fig. 1(b). The resulting simulated
magnetic recordings are shown in the bottom panel of Fig. 1(b).
The eigenspace-projected adaptive spatial filter (Sekihara et al.,2002), mentioned in Typical non-adaptive and adaptive spatial
filters, was applied to these simulated recordings. The data between
0 and 400 ms were used for calculating the covariance matrix, andthe weight vector of the spatial filter was obtained with thiscovariance matrix. The dimension of the signal subspace for the
eigenspace projection was set at three. The spatial filter wasapplied to both the pre- and post-stimulus data, and thereconstructed results are shown in Fig. 2. Fig. 2(a) shows thereconstructed time courses at the three source locations, ŝ(r 1, t),ŝ(r 2, t), and ŝ(r 3, t). Three snapshot reconstructions at 220, 265,
and 300 ms, |ŝ(r, 220)|, |ŝ(r, 265)|, and |ŝ(r, 300)| are shown in Fig.
2(b). The reconstruction – averaged over q
thewhole post-stimulus
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
time window between 0 and 400 ms, bŝs ðr;tÞ2 post (where b post
indicates the time average over the whole post-stimulus period)—isshown in Fig. 2(c). The reconstruction region was defined such that2.5 < y < 2.5 cm and 8 < z < 4 cm with the pixel interval of
0.5 cm, resulting in the total number of pixels equal to 7373.
To apply the proposed statistical thresholding, we first calculate
the modulus of the reconstructed time courses |ŝ(r, t)| and derive an
empirical null distribution from the pre-stimulus portion of |ŝ(r, t)|
at each pixel location. The magnitude time courses |ŝ(r 1, t)|, |ŝ(r 2,t)|, and |ŝ(r 3, t)| are shown in Fig. 3(a). The resultant empirical
distributions at r 1, r 2 , and r 3, expressed as histograms of the pre-I

Fig. 1. (a) The coordinate system and source-sensor configuration used in
the numerical experiments. The coordinate origin was set at the center ofthe sensor coil located at the center of the array. The three point-likesources, shown by the small filled circles, were assumed to be located at (0,1, 6) cm, (0, l, 6) cm, and (0, 1.6, 7.2) cm on the plane of x = 0. The
large circle indicates the projection of the sphere used for the forward
calculation. (b) The first three panels from top to bottom show the timecourses assumed as the time courses of the first, second, and third sources,
respectively. The fourth panel shows the spontaneous MEG added to thegenerated signal magnetic field. The bottom panel shows the simulated
magnetic recordings used for the reconstruction experiments.

372

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

The results of this thresholding with a 5% significance level applied
to the reconstruction results in Fig. 2(b) are shown in Fig. 4(b). The
thresholded time-averaged reconstruction is shown in Fig. 4(c). The
results indicate that the third source, active during both the pre- andpost-stimulus periods, is removed from the post-stimulus reconstruction, verifying the effectiveness of the proposed method.

Experiments
We applied the proposed method to auditory-evoked MEG data
to test its effectiveness. The auditory-evoked fields were measured
using the 275-channel Omega-275i (VSM MedTech Ltd., Port
Coquitlam) whole-cortex biomagnetometer installed at the Biomagnetic Imaging Laboratory, University of California, San
Francisco. The auditory stimulus (1-kHz pure tone) was presented
to the subject’s right ear. The average inter-stimulus interval was 2s, with the interval randomly varied between 1.75 s and 2.25 s. The
sampling frequency was set at 4 kHz, and an on-line filter with abandwidth from 1 to 2 kHz was used. A total of 400 epochs were
measured, and these 400 epochs were averaged to produce theauditory-evoked recordings shown in Fig. 5. Here, although clear
P50m and N100m peaks can be observed, we can see that theseFig. 2. (a) The reconstructed time courses for the first source ŝ(r 1, t)(upper), the second source ŝ(r 2, t) (middle), and the third source ŝ(r 3, t)(bottom). (b) The snapshot reconstruction at 220 ms |ŝ(r, 220)| (upper left),
265 ms |ŝ(r, 265)| (upper right), and 300 ms |ŝ(r, 300)| (lower left). The time
instants at 220, 265, and 300 ms are shown by the three broken verticallines in panel q(a).
(c) The reconstruction averaged over the post-stimulus
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
time window, bŝs ðr;tÞ2 post :stimulus values of |ŝ(r 1, t)|, |ŝ(r 2, t)|, and |ŝ(r 3, t)|, are shown in Fig.
3(b). The empirical distributions are then standardized using Eq.
i(13), and the T max, values are derived from the standardizedi
distributions at all pixel locations. The distribution of T max
expressed as a histogram is shown in Fig. 4(a). Using thisthdistribution, the value of T max
is determined to be 3.40 for asignificance level a of 0.05. Finally, using Eq. (15), we derive thethreshold ~ (r), and when |ŝ(r, t)| < ~ (r), |ŝ(r, t)| is set equal to zero.

Fig. 3. (a) Reconstructed magnitude time courses, |ŝ(r 1, t)|, |ŝ(r 2, t)|, and|ŝ(r 3, t)|, are shown from top to bottom, respectively. The horizontal brokenlines in the upper two panels show the threshold values at r 1 and r 2. (b)Histograms of the pre-stimulus values of |ŝ(r 1, t)| (upper left), |ŝ(r 2, t)|(upper right), and |ŝ(r 3, t)| (lower left).

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

373

averaged recordings contain a considerable amount of periodic
background activity.
The eigenspace-projected adaptive beamformer (Sekihara et al.,2002) was applied to these averaged recordings. The dimension of
the signal subspace was set at 2. The data between 0 to 200 ms
were used to calculate the covariance matrix R, and the weightvector was obtained with this covariance matrix. Two latencies, 44
ms and 86 ms, are selected. One is near the peak of P50m, and theother is near the peak of N100m. These time points are shown in
the two vertical broken lines in Fig. 5. The snapshots of the sourcereconstruction results at these latencies are shown in Figs. 6(a) and(b). Both sets of the results contain a clear localized source in theleft temporal lobe probably near the primary auditory area.
However, the reconstruction results at 44 ms also contain otherdiffuse activity. The time-averaged reconstruction
obtained fromqﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
the whole pre-stimulus period (400 0 ms), bŝs ðr;tÞ2 pre ; is shownin Fig. 6(c) (where ( )pre indicates the time average over this prestimulus period). These results also contain a diffuse source similar
to that found in the snapshot at 44 ms, suggesting that these diffuse
sources were caused by the periodic background activities
observed in the waveforms in Fig. 5.
To apply the proposed statistical thresholding, we use themeasurements taken during the pre-stimulus period as a control. We
reconstruct the time course |ŝ(r, t j )| of the pre-stimulus period, andthe empirical distribution of |ŝ(r, t j )| is calculated at each pixel
location. The empirical distributions are then standardized using Eq.
i(13) to determine T max
for the ith pixel, and the same procedure isirepeated at all pixel locations to derive the distribution of T max. Theiempirical distribution of T max is shown as a histogram in Fig. 7(a).
Here, the number of pixels in the three-dimensional reconstructionthgrid was 6800. The value of T max
was determined to be 3.34 using
this empirical distribution. Then, using Eq. (15), we derive thethreshold R(r). When |ŝ(r, t)| < R (r), |ŝ(r, t)| is set to zero. The
results of this thresholding applied to the snapshot images in Figs.
6(a) and (b) are shown in Figs. 7(b) and (c), respectively. Here, thesignificance level a was set at 0.05. The diffuse source activities
contained in the results in Fig. 6(a) have been removed in Fig. 7(b).
These results demonstrate the effectiveness of the proposed
statistical thresholding in removing the influence of backgroundsource activities in the reconstruction results.

IithFig. 4. (a) Histogram of T max. The value of T max
is determined to be 3.40
for a significance level a of 0.05 using this distribution. (b) Thresholded
results with a = 0.05 for |ŝ(r, 220)| (upper left), |ŝ(r, 265)| (upper
right), andqﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ|ŝ(r, 300)| (lower left). (c) The thresholded reconstruction for bŝs ðr;tÞ2 post :Discussion and conclusion
We have developed a simple and novel nonparametric statisticalthresholding procedure for tomographic reconstruction results from

Fig. 5. (a) The 400-epoch-averaged auditory-evoked fields measured using the 275-channel sensor array. Among the 275 sensor recordings, the recordings from132 sensors covering the subject’s left hemisphere are displayed.

374

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

Fig. 6. The maximum-intensity projections of the source reconstruction obtained using the eigenspace-projected adaptive spatial filter. (a) Snapshot at 44 ms oflatency, (b) snapshot at 86 ms of latency, and (c) time averaged reconstruction obtained from the whole pre-stimulus period (400 and 0 ms). The left column
shows the maximum intensity projections of the three-dimensional reconstruction onto the axial plane. The middle column shows those onto the coronal plane.
The left column shows those onto the sagittal plane. The upper ease letters L and R show the left and the right hemispheres, respectively.

MEG data. Several methods with their aims similar to the one forthe proposed method have recently been proposed (Greenblatt andPflieger, 2004; Pantazis et al., 2003; Singh et al., 2003).
Particularly, the methods in (Greenblatt and Pflieger, 2004;
Pantazis et al., 2003) can be alternative methods to the one
proposed in this paper. These methods use permutation tests toassess the statistical significance of the source reconstruction
results. Thus, it not only requires intensive computer time but also
requires raw epoch data to be stored. On the contrary, the method
proposed in this paper does not use such computer intensive resampling methods as the permutation tests or the bootstrap and canthus be implemented with much less computer time. In addition,
our method does not require raw epoch data to be stored but uses
rather only the averaged data.
The proposed method uses the maximum statistics to address
the multiple comparison problems. Without such a multiple
comparison procedure, 100a% of pixels may exhibit falsepositive activations in a single snapshot image (where a is thelevel of significance). The multiple comparison procedure,
however, reduces this false-positive probability to the level atwhich only 100a% of snapshot images may contain the falsepositive activations. However, since the proposed procedure does
not perform the multiple comparisons in the temporal dimension,a reconstructed time course may still contain 100a% of time
points that exhibit false activations. This fact can be seen in Fig.
3(a) in which the derived threshold values are indicated by thehorizontal broken lines in the upper two panels. Particularly in the
middle panel, several time points in the pre-stimulus period

exceed the threshold value, resulting in the false-positive
activations at these time points. To avoid such false activations,the multiple comparisons should also be performed in thetemporal dimension, and for such space – time multiple comparisons, an empirical distribution should be derived not only at eachspatial point but also at each time point. Such a method forspace – time multiple comparisons is currently under investigation,with the results to be published in the near future.
We pointed out in Proposed nonparametric statistical significance evaluation that we impose the assumption that the timecourse reconstruction |ŝc (r, t j )| is identically and independently
distributed (IID) to derive an empirical distribution of |ŝc (r, t j )|.
This assumption is essential for the procedure without multiple
comparisons. However, once we take T max (r) at each pixel
location, the distributional properties of the process do not matter
anymore, and the assumption that the time course reconstruction
is IID is not needed for our procedure with multiple comparison.
In other words, even if the pre-stimulus process, ŝ c (r, t j ), hasstrong temporal autocorrelations, we can still compare poststimulus values to the maximal values of the pre-stimulus
process. However, instead of the IID assumption, our multiple
comparison procedure imposes a new assumption that the shapeof the control source distribution is the same at each pixel,
although its scale may be different. This assumption is moreattainable than the IID assumption because brain source timecourses are known to have significant autocorrelations, but suchautocorrelations should be more or less similar at each pixel
location.

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

375

iFig. 7. (a) Histogram of T max
from all pixel locations. (b) The results of the proposed statistical thresholding applied to the snapshot shown in Fig. 6(a). (c) The
results of the proposed statistical thresholding applied to the snapshot shown in Fig. 6(b).

The other implicit assumption needed for our proposed
procedure is that the source activities in the control (pre-stimulus)
period are not correlated with those in the task (post-stimulus)
period, that is, the control source activities are stationary
throughout the task and control periods. This assumption is notexactly hold in such a situation that background (control source)
activities are changed before and after the stimulus. Actually, we
can often observe that the alpha activity is suppressed by a visual
stimulus. In such situations, however, we can assume that thestimulus will also evoke non-phase locked activity that should be
accounted for separately from the ‘‘pure’’ control source activities
that are unchanged by the stimulus. The proposed method can
only remove the influence of the brain sources that are
uninfluenced by the stimulus, i.e., ‘‘pure’’ control source activities,
and it cannot remove the influence of the stimulus-evoked nonphase locked activity. The proposed method should work foraverage evoked responses because, in such averaged responses,
the contribution of the stimulus-evoked non-phase locked activity
is much reduced by averaging, and its influence generally can beneglected.
In Proposed nonparametric statistical significance evaluation,
the proposed method is described using the one-tailed test.
However, the two-tailed test may be more appropriate when target
source activities can vary in the both directions, i.e., when target
sources can be deactivated by stimuli. The extension of the

proposed method to include the two-tailed test can be performedby first calculating the minimum T value T min(r) at each pixel
location where the value of T is calculated using Eq. (13). Denotingithe minimum T value at the ith pixel location as T min, we next
obtain the empirical distribution of T min, Ĥ min(x) such thatiĤ min(x) = #{T minx} / K N and then obtain the threshold of thethth
T min value for the a-significance level, T min, such that T min=
1
Ĥ min (1  a / 2). We can obtain the statistical threshold corresponding to T thmin, R min(r), by using thr ðrÞ þ bjŝs c r;tj jc :r̂
~ min ðrÞ ¼ Tmin

ð16Þ

When |ŝ(r, t)| R min(r), we can conclude that there is a statisticallysignificant decrease in the source activity at r and time t. The otherside of the threshold can be obtained in exactly the same manner as
described in Proposed nonparametric statistical significance evaluathtion except using T max= Ĥ 1(1  a / 2).
Although the proposed method has been described as a method
for assessing the statistical significance of the source reconstruction
results, the method can also be applied to determine the statisticalsignificance of b I (t) in the sensor recordings b(t). The null
hypothesis here is b I (t) = 0 in the task measurements. Let us
denote the kth component of the vectors b(t), b c (t), and b I (t) asb k (t), b ck (t), and bkI (t), respectively. To determine the statistical

376

K. Sekihara et al. / NeuroImage 27 (2005) 368 – 376

significance of bkI (t), the empirical distribution of |b ck (t j )| is firststandardized by calculating    jbkc tj j  bjbkc tj jck
T c tj ¼
;
ð17Þ
r kb
r̂
where
 k 2 2
  2
r b ¼ bbkc tj c  bjbkc tj jc ;r̂It j is a discrete time point in the control period, and b c again
indicates the time average over the control period. We thencalculate the maximum T value for the kth channel recording,
k
k
T max
, and obtain the empirical distribution of T max, Ĥ b (x), suchk1that Ĥ b (x) = {T max x} / M . We then obtain the threshold of thethth
T max value for the a-significance level T max
such that T max=
1Ĥ b (1  a) and finally obtain the statistical threshold for the kth
sensor recording, R bk , by using thr kb þ bjbkc tj jc :r̂
~ kb ¼ Tmax
ð18Þ
Therefore, when b k (t) > R bk , the null hypothesis is rejected and we
can conclude that a statistically significant b I (t) exists.
In conclusion, this paper proposes a simple and efficient methodof statistical thresholding for MEG spatio-temporal source reconstruction. The method assumes the neural source activities to bequasi-stochastic, and it derives an empirical distribution of nontarget sources from the time course reconstruction in the control
period. It then derives the value of statistical threshold based on theempirical distribution with the multiple comparisons taken into
account. In summary, the method provides an effective means ofremoving source activities not of interest to the current measurements and of extracting the source activity of interest.

Acknowledgments
This work has been supported by Grants-in-Aid from the
Ministry of Education, Science, Culture and Sports in Japan(Cl3680948 and Cl6500296) for KS. This work has been supported
in part by the Whitaker Foundation and by the National Institute of
Health (P4lRR12553-03 and R0l-DC004855-0lAl) for SN.

References
Baillet, S., Mosher, J.C., Leahy, R.M., 2001. Electromagnetic brain
mapping. IEEE Signal Process. Mag. 18, 14 – 30.
Barnes, G.R., Hillebrand, A., 2003. Statistical flattening of MEG beamformer images. Hum. Brain Mapp. 18, 1 – 12.
Blair, B.C., Karniski, W., 1994. Distribution-free statistical analysis ofsurface and volumetric maps. In: Thatcher, R.W. (Ed.), Functional
Neuroimaging. Academic Press.
Dale, A.M., Liu, A.K., Fisehl, B.R., Buekner, R.L., Belliveau, J.W.,
Lewine, J.D., Halgren, E., 2000. Dynamic statistical parametricmapping: combining fMR1 and MEG for high-resolution imaging ofcortical activity. Neuron 26, 55 – 67.
Friston, K.J., Holmes, A.P., Worsley, K.J., Polin, J.-P., Frith, C.D.,

1k#{T maxto x.

kx} is the number of T max
values which are less than or equal

Frackowiak, R.S.J., 1995. Statistical parametric maps in functional
imaging: a general linear approach. Hum. Brain Mapp. 2, 189 – 210.
Greenblatt, R.E., Pflieger, M.E., 2004. Randomization-based hypothesis
testing from event-related data. Brain Topogr. 16, 225 – 232.
Gross, J., Kujara, J., Hämäläinen, M., Timmermann, L., Sehnitzler, A.,
Salmelin, R., 2001. Dynamic imaging of coherent sources: studying
neural interactions in the human brain. Proc. Natl. Acad. Sci. 98,
694 – 699.
Hämäläinen, M.S. Ilmoniemi, R.J., 1984. ‘‘Interpreting Measured Magnetic
Fields of the Brain: Estimates of Current Distributions,’’ Tech. Rep.
TKK-F-A559, Helsinki University of Technology.
Hämäläinen, M., Hari, R., IImoniemi, R.J., Knuutila, J., Lounasmaa, O.V.,
1993. Magnetoencephalography—Theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod.
Phys. 65, 413 – 497.
Hashimoto, I., Kimura, T., Tanosaki, M., Iguehi, Y., Sekihara, K., 2003.
Muscle afferent inputs from the hand activate human cerebellum
sequentially through parallel and climbing fiber systems. Clin. Neurophysiol. 114, 2107 – 2117.
Ishii, R., Dziewas, R., Chau, W., Soros, P., Okamoto, H., Gunji, A.,
Pantev, C., 2003. Current source density distribution of sleep spindles
in humans as found by synthetic aperture magnetometry. Neurosci.
Lett. 340, 25 – 28.
Nichols, T.E., Holmes, A.P., 2001. Nonparametric permutation tests forfunctional neuroimaging: a primer with examples. Hum. Brain Mapp.
15, 1 – 25.
Pantazis, D., Nichols, T.E., Baillet, S., Leahy, R.M., 2003 (July).
Spatiotemporal localization of significant activation in MEG using
permutation tests. In: Taylor, C.J., Noble, J.A. (Eds.), Proceedings of the
18th Conference on Information Processing in Medical Imaging.
Springer Series in Medical Imaging, Ambleside, UK, pp. 512 – 523.
Pascual-Marqui, R.D., 2002. Standardized low resolution brain electromagnetic tomography (sloreta): technical details. Methods Find. Exp.
Clin. Pharmacol. 24, 5 – 12.
Robinson, S.E., Vrba, J., 1999. Functional neuroimaging by synthetic
aperture magnetometry (SAM). In: Yoshimoto, T., et al., (Eds.), Recent
Advances in Biomagnetism. Tohoku Univ. Press, Sendai, pp. 302 – 305.
Sekihara, K., Nagarajan, S.S., 2004. Neuromagnetic source reconstruction and inverse modeling. In: He, B. (Ed.), Modeling and Imaging
of Bioelectric Activity—Principles and Applications. Kluwer Academic/Plenum Publishers, pp. 213 – 250.
Sekihara, K., Seholz, B., 1996. Generalized Wiener estimation of threedimensional current distribution from biomagnetie measurements. In:
Aine, C.I., et al., (Eds.), Biomag 96: Proceedings of the TenthInternational Conference on Biomagnetism. Springer-Verlag, New
York, pp. 338341.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., Miyashita, Y.,
2001. Reconstructing spatio-temporal activities of neural sources using
an MEG vector beam-former technique. IEEE Trans. Biomed. Eng. 48,
760 – 771.
Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., Miyashita, Y., 2002.
Application of an MEG eigenspace beamformer to reconstructing
spatio-temporal activities of neural sources. Hum. Brain Mapp. 15,
199 – 215.
Sekihara, K., Sahani, M., Nagarajan, S.S., 2005. Location bias and spatial
resolution of adaptive and non-adaptive spatial filters for MEG sourcereconstruction. Neurolmage 25, 1056 – 1067.
Singh, K.D., Barnes, G.R., Hillebrand, A., 2003. Group imaging of taskrelated changes in cortical synchronisation using nonparametricpermutation testing. NeuroImage 19, 1589 – 1601.
van Veen, B.D., van Drongelen, W., Yuehtman, M., Suzuki, A., 1997.
Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44, 867 – 880.
Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C., Friston, K.J., Evans,
A.C., 1996. A unified statistical approach for determining significant
signals in images of cerebral activation. Hum. Brain Mapp. 4, 58 – 73.