IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 11, NOVEMBER 2011

3121

Removal of Spurious Coherence in MEG
Source-Space Coherence Analysis
Kensuke Sekihara*, Fellow, IEEE, Julia P. Owen, Stephan Trisno, and Srikantan S. Nagarajan, Senior Member, IEEE

Abstract—Source-space coherence analysis has become a popular method to estimate functional connectivity based on
MEG/EEG. Source-space analysis involves solving the inverse
problem, estimating the time courses of specific brain regions,and then examining the coherence between activities at different
brain regions. However, source-space coherence analysis can beconfounded by spurious coherence caused due to the leakage properties of the inverse algorithm employed. Such spurious coherence
is typically manifested as an artifactual large peak around the seedvoxel, called seed blur, in the resulting coherence images. This seedblur often obscures important details of brain interactions. This
paper proposes the use of the imaginary part of the coherence toremove the spurious coherence caused by the leakage of an imaging
algorithm. We present a theoretical analysis that explains how theuse of imaginary part can remove this spurious coherence. We thenpresent results from both computer simulations and experiments
using resting-state MEG data which demonstrate the validity of
our analysis.
Index Terms—Bioelectromagnetism, coherence analysis,
functional connectivity, magnetoencephalography (MEG), MEG
source imaging, source coherence.

I. INTRODUCTION
HERE has been tremendous interest in estimating the functional connectivity of neuronal oscillations across brain regions based on electromagnetic measurements such as magnetoencephalography (MEG) and electroencephalography (EEG).
However, prior studies using EEG or MEG have largely employed sensor-space analysis, in which brain interactions havebeen analyzed using raw sensor recordings [1]–[3]. In sensorspace analysis, the field spread across many sensors from asingle brain region leads to uncertainties in interpreting the estimation results of brain interactions, as pointed out in [4].
Recently, a number of studies have begun to use sourcespace analysis. In source-space analysis, voxel time courses
are first estimated by solving the inverse problem and brain

T

Manuscript received April 25, 2011; revised June 8, 2011; accepted July 3,
2011. Date of publication August 4, 2011; date of current version October 19,
2011. The work of K. Sekihara was supported by the Ministry of Education,
Science, Culture and Sports, Japan, under Grant 20500394 and Grant 22659210.
The work of S. S. Nagarajan was supported by the National Institute of Healthunder Grant R01-DC004855-01A1 and Grant R01-NS006435. Asterisk indicates corresponding author.
*K. Sekihara is with the Department of Systems Design and Engineering, Tokyo Metropolitan University, Tokyo 191-0065, Japan (e-mail:k-sekihara@nifty.com).
J. P. Owen, S. Trisno, and S. S. Nagarajan are with the Department of Radiology and Biomedical Imaging, Biomagnetic Imaging Laboratory, Universityof California, San Francisco, CA 94143 USA (e-mail: juliapowen@gmail.com;
strisno@berkeley.edu; srikantan.nagarajan@ucsf.edu).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TBME.2011.2162514

interactions are then analyzed using those estimated voxel timecourses [4]–[8]. Although a certain degree of inaccuracy exists
in the source inversion process, the source-space analysis has
the potential of providing more accurate information regarding
which brain regions are functionally coupled.
Either in sensor- or source-space analysis, some kind of ameasure for the interaction, called a connectivity metric, must
be computed. Among existing measures a widely used representative measure is coherence [5]–[7], [9]–[11]. In source-spacecoherence analysis, a typical procedure involves first setting areference point, called the seed point, and computing the coherence between the time courses from the seed point and anothervoxel’s location, referred to as the target location in this paper.
By scanning through all target locations in a brain, a 3-D mapping of source coherence, namely a source-coherence image,with respect to the seed location can be obtained.
A serious problem in source-coherence imaging arises fromthe spurious coherence caused by the leakage of an inverse
algorithm, and such leakages are more or less inevitable inany inverse algorithm [12]. One representative ramification ofthis spurious coherence is an artifactual large peak around theseed voxel, called seed blur, in the resulting coherence image.
Quite often, the seed blur dominates the resultant coherence
images, and obscures important details of the brain interactions.
Examples of the seed blur are shown in our computer simulations
and experiments using resting-state MEG data.
To remove such spurious coherence, this paper proposes theuse of the imaginary part of coherence, which is called the imaginary coherence. The use of imaginary coherence was originallyproposed by Nolte et al. [13] to remove the spurious coherence
caused by the volume conduction in EEG sensor-space coherence analysis. This paper extends their sensor-space analysis ofimaginary coherence to source-space coherence analysis, andpresents a detailed theoretical analysis that explains how theuse of imaginary coherence leads to the removal of the spurious
coherence caused by the leakage of source imaging algorithms.
The validity of our analysis is then demonstrated by our computer simulation and by experiments using resting-state MEG
data.
II. METHOD
A. Source-Space Coherence Analysis
Source-space coherence analysis (or source-coherence imaging) requires a two-step procedure. In the first step, voxel timecourses are estimated using an inverse algorithm. In the second
step, the coherence measure is computed using the voxel timecourses estimated in the first step. In this paper, the coherence

0018-9294/$26.00 © 2011 IEEE

3122

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 11, NOVEMBER 2011is computed between the time courses from the seed voxel andfrom each of the target voxel. Let us define the time course fromthe seed voxel as uS (t) and the time course from a target voxel
as uT (t), and let the spectra of the seed and the target voxels be
σS (f ) and σT (f ), respectively.
The coherence η(f ) is obtained by computing the correlationof these spectraη(f ) = σT (f )σS∗ (f )

The magnitude coherence between the seed and the target
voxels is expressed as





σT σ
S∗ 


|
η| =  
(5)
.
2
2
 |
σT | |
σS |  
Using (4), we haveS∗  = σT σS∗ 

σT σ

(1)

+ d1 |σT |2  + d2 |σS |2  + d1 d2 σS σT∗ (6)where the superscript ∗ indicates the complex conjugate, and thebracket · indicates the ensemble average. In practical applications, this ensemble average is computed by averaging acrossmultiple trials. When only a single continuous dataset is measured, the single dataset is divided into many trials and coherence is obtained by averaging across these trials.
In conventional source-coherence imaging, the magnitude coherence |η(f )| is computed and mapped by scanning throughthe target voxels. It is, however, apparent from (1) that if theseed and target spectra contain common components that do
not result from true brain interactions, then the magnitude coherence may contain spurious components. When coherence iscomputed between estimated voxel time courses, the leakage ofthe imaging algorithm used is a major source of such spurious
coherence.

|
σT |2  = |σS |2  + d21 |σT |2  + 2d1  (σT σS∗ )

(7)

|
σS |2  = |σT |2  + d22 |σS |2  + 2d2  (σT σS∗ )(8)

|σT (f )|2 |σS (f )|2 

andwhere  ( · ) indicates the real part. Therefore, even when there
is no true source interaction, i.e., even when σT σS∗  = 0 andη | has a nonzero value, which is equal toσS σT∗  = 0, ||
η| = 

.

(9)

σT σS∗  + σS σT∗  = 2 (σS σT∗ )1) Leakage Effects in the Magnitude Coherence: The estimated source time courses at the seed voxel u
S (t) and at thetarget voxel u
T (t) can be expressed asuS (t) = uS (t) + d1 uT (t) + cS (t)(2)

uT (t) = uT (t) + d2 uS (t) + cT (t)(3)

where uS (t) and uT (t) are the true source time courses of theseed and the target locations, respectively.1 In the aforementioned equations, d1 uT (t) indicates the leakage of the target
signal in the estimated seed signal, and d2 uS (t) indicates theleakage of the seed signal in the estimated target signal. The
real-valued constants d1 and d2 express the relative amount ofthese leakages.
In the aforementioned equations, cS (t) and cT (t) express theinterference terms, which may include the leakage from othersources, contributions from external disturbances, and that of
the sensor noise. The influence of these interference terms canbe considered separately from the influence of the leakage terms
in Section V. In the arguments here, we ignore cS (t) and cT (t)and the estimated spectra at the seed and target voxels, σ
S (f )and σ
T (f ), are, respectively, expressed asσT = σT + d2 σS

the cross spectrum σT σ
S∗  can be expressed asS∗  = (1 − d1 d2 )σT σS∗  + d1 |σT |2 σT σ
+ d2 |σS |2  + 2d1 d2  (σT σS∗ ) .

(10)

By taking the imaginary part of (10), we can derive

and

and

(|σS |2  + d21 |σT |2 ) (|σT |2  + d22 |σS |2 )

Such spurious coherence is typically manifested as an artifactual large peak around the seed voxel, called seed blur, in the
resulting coherence images.
2) Leakage Effects in the Imaginary Coherence: Using (6)and the relationship

B. Leakage Effects of Imaging AlgorithmsσS = σS + d1 σT|d1 |σT |2  + d2 |σS |2 |

 (
σT σ
S∗ ) = (1 − d1 d2 ) (σT σS∗ )(11)where  ( · ) indicates the imaginary part. The aforementioned
equation indicates that, aside from the multiplicative constant,the imaginary part of the cross spectrum between the seed andtarget voxels is equal to  (σT σS∗ ), which is the imaginary
part of the true cross spectrum.
S∗  is normalTo computecoherence, the cross spectrum 
σT σ

2
2
σT | |
σS | . Thus, we obtain the imaginary partized by |of the estimated coherence  (η ) as (η) = 
=

 (
σT σ
S∗ )
σS |2 |σT |2 |(1 − d1 d2 ) (σT σS∗ )= Ω  (η)σS |2 |σT |2 |(12)

where  (η) indicates the true value of the imaginary coherence.
Using (7) and (8), Ω is obtained as

(4)

Ω=

where we omit the explicit notation of (f ) for simplicity.

(1 − d1 d2 )
√
φ1 φ2

(13)

where
1 In this paper, the estimated values are indicated bytrue values.

 to distinguish their

φ1 = 1 + d21

|σT |2  (σT σS∗ )+ 2d1
2
|σS | 
|σS |2 (14)

SEKIHARA et al.: REMOVAL OF SPURIOUS COHERENCE IN MEG SOURCE-SPACE COHERENCE ANALYSIS

and|σS |  (σT σS∗ )
φ2 = 1 + d22+
2d.
2
|σT |2 
|σT |2 
2(15)

Equation (12) shows that when  (η) = 0, we have  (η ) = 0,indicating that no spurious coherence has been generated.
However, (12) also indicates that the value of  (η ) differsfrom the true value  (η), i.e., the amplitude of the estimatedimaginary coherence is biased and the bias is represented by
Ω, which is called the amplitude bias factor in this paper. Let
us calculate the amplitude bias factor assuming a simple scenario in which d1 = d2 = d and |σS |2  = |σT |2 . Under theseassumptions, the bias factor is simplified to
1 − d2
1 − |d|1 + |d|≥Ω=≥.
21 − |d|1 + d + 2(η)d
1 + |d|(16)

The aforementioned equation shows how the amount of leakage|d| affects the amplitude bias Ω, and it shows that if |d| issmall such as |d| < 0.1, the amplitude bias is less than 10%.
The narrow-band adaptive spatial filter used in our investigation
generally has small leakages [12], and thus, this amplitude bias is
generally small. We check this fact in our computer simulation.
C. Statistical Thresholding of an Imaginary-Coherence Image
In practical applications, we need to assess the statisticalsignificance of the imaginary coherence. In this paper, the surrogate data method [14], [15] is used for this assessment. In thismethod, the surrogate voxel spectra are created by multiplying
random phases with the original voxel spectra. The surrogate
spectra for the seed and the target voxels σ
S (f ) and σ
T (f ) are
expressed asσS = σS e2π ε SandσT = σT e2π ε T(17)where εS and εT are the uniform random numbers between 0 and
1, and the explicit notation of (f ) is again omitted for simplicity.
Note that the surrogate spectra have the same original power
spectra but the phase relationship is destroyed by multiplying
the random phases to the original spectra. The value of theimaginary coherence calculated using the surrogate spectra σ
Sand σ
T is expressed as


 
σT σ
S∗ e2π Δ ε S∗ )
 (
σT σ
= 
(18)
ξ=
|
σT |2 |σS |2 |σT |2 |σS |2 where Δε = εT − εS . The generation of ξ is repeated with
B times, and a total B values of ξ, which are denotedξ 1 , ξ 2 , . . . , ξ B , are obtained. These ξ 1 , . . . , ξ B can form an empirical null distribution at each target voxel. We then could
derive a voxel-by-voxel statistical threshold using this empirical null distribution. However, the statistical threshold derived inthis manner does not take the multiple comparisons into accountand it generally leads to a situation in which many false-positive
voxels arise, i.e., many voxels that contain no brain interaction
are found to be interacting. To avoid this problem, the statisticalsignificance is determined using a procedure that takes multiple
comparisons into account. For this purpose, we use the maximal
statistics [16], [17].

3123

To utilize maximum statistics, the values ξ β (β = 1, . . . , B)
are first standardized and converted into pseudo-t values, suchthat
ξ β − ξ¯
(19)
Tβ =
σξ
where ξ¯ and σξ2 are the average and the variance of ξ β (β =1, . . . , B), respectively. Since these ξ¯ and σξ2 are obtained ateach target voxel, the values at the jth target voxel are denoted¯ and σ 2 (j). The maximum value of T β obtained at theas ξ(j)ξ
jth voxel is denoted as T m ax (j). Denoting a total number ofvoxels NV , we have T m ax (1), . . . , T m ax (NV ) to form the null
distribution. We then sort these values in an increasing orderT m ax (1̃) ≤ T m ax (2̃) ≤ · · · ≤ T m ax (ÑV )where T m ax (k̃) is the kth minimum value. We set the levelof the statistical significance to α, and choose T m ax (p̃) wherep̃ = αNV where αNV indicates the maximum integer not
greater than αNV . The threshold value for the jth voxel ξ th (j)is derived as¯ξ th (j) = T m ax (p̃)σξ (j) + ξ(j).

(20)

We evaluate, at the jth voxel, the statistical significance of theimaginary coherence by comparing | (η ) | with |ξ th (j)|. Whenth| (η ) | > |ξ (j)|, the imaginary coherence at the jth voxel
 (η ) is considered to be statistically significant. Conversely,η ) |,  (η ) is not considered to be statisticallyif |ξ th (j)| > | (significant, and it is set to 0.
III. COMPUTER SIMULATION
A. Data Generation
Computer simulations were performed to verify the arguments in Section II-B. In our simulations, we used a sensor alignment of the 275-sensor array from the Omega (VMS Medtech
Ltd., Coquitlam, BC, Canada) neuromagnetometer. A single
vertical plane (x = 0 cm) was assumed at the middle of thewhole-head sensor array, and three sources were assumed toexist on this plane. The source–sensor configuration and the coordinate system are depicted in Fig. 1. The (y, z) coordinatesof the three sources were (−1.0, 9.5), (1.5, 10.0), and (1.0, 7.5)cm, respectively.
Multiple-trial measurements were simulated, in which a totalof 120 trial recordings were generated. Each trial consists of 600
time points where the sampling interval was assumed to be 2 ms.
The time courses of the three sources had the same shape withtrial-to-trial time jitters. The time jitter between the first andsecond sources and the jitter between the second and the thirdsources were both generated using Gaussian random numberswith a standard deviation equal to 15 time points. There is a
50-time-point mean delay between the time courses of the first
and the second sources, and also between those of the second
and the third sources. Here, the amount of jitter controls thevalue of coherence, and setting the jitters in the aforementioned
manner results in the following values: |(η1,2 )| = 0.27 and|(η2,3 )| = 0.23 where η1,2 indicates the coherence between

3124IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 11, NOVEMBER 2011

Fig. 3. Results of reconstructing source power distribution on the plane x = 0
cm. (a) SIR was set equal to 4.0. (b) SIR was set equal to 0.25. Asterisks show
the location of the seed that was set for computing the coherence images shownin Fig. 4.

Defining the jth-source time course given as uj (t), the signal
magnetic recordings bs (t) were computed usingbs (t) =3


Aj uj (t)lj

(21)

j =1

Fig. 1. Coordinate system and source–sensor configuration used in the computer simulation. The coordinate origin was set at the center of the sensor coillocated at the center of the array. The plane at x = 0 cm is shown. The large
circle shows the cross section of the sphere used for the forward calculation,and the small circles show the locations of the three sources. The square shows
the region for which reconstruction results are shown in Figs. 3–5.

where lj is the lead field vector of the jth source and Aj is its
amplitude. Here, the powers of the three sources were set equalin the sensor domain, i.e., the relationshipA21 l1 2 = A22 l2 2 = A23 l3 2
held. The lead field was calculated using the spherical homogeneous conductor model [18] with the sphere origin set to(0, 0, 4) cm.
The simulated sensor recordings b(t) were generated byadding spontaneous MEG signal to the computed signal recordings bs (t), such thatb(t) = bs (t) + γbI (t)(22)

where bI (t) is the spontaneous MEG measured using the same
275 whole-head sensor array, and γ is a constant that controls
the signal-to-interference ratio (SIR) of the simulated sensor
recordings. In our computer simulations, an SIR was set equal
to 4.0 and 0.25, where SIR is defined as the ratio of the average signal power  bs (t) 2  to the average interference power
 bI (t) 2 . Note that the coherence analysis is usually performed using non-averaged trial data. Therefore, an SIR of 4.0
is considered a significantly high SIR and the SIR of 0.25 canbe considered a typical SIR for non-averaged trials.
B. Source-Coherence Imaging
Fig. 2. Time courses assigned to the three sources in the computer simulation.
The time courses in the top, middle, and bottom panels are those for the first,
second, and the third sources, respectively. These time courses have trial-to-trial
time jitters, and the time courses in the first three trials are shown. The solid,
broken, and dashed-dotted lines, respectively, show the time courses for the first,
second, and third trials.

the first and second sources and η2,3 that between the second andthird sources. Also, (η1,2 ) and (η2,3 ) indicate the imaginaryparts of these coherence values. The source time courses in thefirst three trials are shown in Fig. 2.

The voxel time-course estimation was performed using thenarrow-band adaptive spatial filter [12], [19] with a data covariance matrix obtained with a frequency band between 8 and
12 Hz. This is because the spectra of source time courses arebandlimited to this bandwidth. Reconstructed source power images on the plane x = 0 cm are shown in Fig. 3 for both SIR
cases. The reconstructed images show that the three sources areresolved for the both SIR cases. We then set a seed point at
the location of the second source to compute coherence images.
The seed point is marked by asterisks in Fig. 3.
The resultant coherence images are shown in Fig. 4. To obtainthese results, statistical thresholding using the surrogate data

SEKIHARA et al.: REMOVAL OF SPURIOUS COHERENCE IN MEG SOURCE-SPACE COHERENCE ANALYSIS

Fig. 4. Results of imaging the source coherence on the plane x = 0 cm.
The seed was set at the second source location indicated by the asterisk.
(a) Magnitude-coherence image when SIR was set equal to 4.0. (b) Imaginarycoherence image when SIR was set equal to 4.0. (c) Magnitude-coherence image
when SIR was set equal to 0.25. (d) Imaginary-coherence image when SIR was
set equal to 0.25. Statistical thresholding using the surrogate data method described in Section II-C was applied with the number of surrogate datasets to
200 and the level of the statistical significance set to 0.99.

method described in Section II-C was applied. Here, the numberof surrogate datasets B was set to 200, and the level of thestatistical significance α was set to 0.99. The upper panels in
Fig. 4(a) and (b) show the results for the case of the SIR equal
to 4.0, and the lower panels in Fig. 4(c) and (d) show the results
for the case of the SIR equal to 0.25. Also, the left panels in
Fig. 4(a) and (c) show the magnitude-coherence images and theright panels in Fig. 4(b) and (d) show the imaginary-coherence
images.
When the SIR is equal to 4.0, i.e., when the SIR is very high,the magnitude-coherence image [see Fig. 4(a)] manages to show
the first and the third sources, although a large seed blur exists atthe second source location. That is, the two sources interactingwith the second source can be observed in spite of the existenceof the seed blur when the SIR is very high. In the imaginarycoherence image [see Fig. 4(b)], the intensity of the seed blur ismuch reduced, and the first and the third sources, which interact
with the second source, are clearly shown.
When the SIR is equal to 0.25, i.e., when the SIR is typical of coherence analysis, the seed blur dominates and it obscures the other sources in the magnitude-coherence image [seeFig. 4(c)]. Thus, it is difficult to obtain information on the interacting sources from the magnitude-coherence image. On thecontrary, in the imaginary-coherence image [see Fig. 4(d)], theintensity of the seed blur is much reduced and the two sourcesthat interact with the second source can clearly be observed.
The narrow-band adaptive spatial filter generally has a small
leakage [12]. In this computer simulation, the leakage constants
d1 and d2 are approximately equal to 0.01 between the first
and the second sources, and also between the second and thethird sources. Therefore, the bias factor Ω in (16) is nearly
equal to 1, and the amplitude bias should be small. Actually, in

3125

Fig. 5. Results of imaging the source coherence on the plane x = 0 cm when
the seed was set at a location where no active source existed. The first source
was made inactive, and the seed was set at the first source location indicatedby the asterisks. (a) Reconstructed source power distribution. (b) Magnitudecoherence image. (c) Imaginary-coherence image. Statistical thresholding using
the surrogate data method described was applied with the number of surrogate
datasets to 200 and the level of the statistical significance set to 0.99.

our computer simulation, the estimated imaginary coherence,η3,2 )| are equal to 0.275 and 0.23 when the|(η1,2 )| and |(SIR was set to 4, and 0.26 and 0.21 when the SIR was set
to 0.25. These values of the estimated imaginary coherence arevery close to their true values, |(η1,2 )| = 0.27 and |(η2,3 )| =0.23, and thus, the amplitude bias is not evident in this computer
simulation.
Next, we performed computer simulation in which no source
activity existed at the seed location. Here, the first source wasmade inactive and the seed was set at the first source location.
The results of the source power reconstruction in this case areshown in Fig. 5(a). The results contain the second and the thirdsources but do not contain the first source. The seed location
is indicated by the asterisk at the first source location. The
thresholded magnitude-coherence and the imaginary-coherence
images are shown in Fig. 5(b) and (c). The magnitude-coherence
image [see Fig. 5(b)] shows a large blurred peak at the seedlocation, as if the large region surrounding the seed location is
synchronously active. This large region represents the spurious
coherence expressed in (9). The imaginary-coherence image[see Fig. 5(c)], however, shows no activity.
IV. EXPERIMENTS WITH RESTING-STATE MEG DATA
We performed source-coherence imaging using resting-state
MEG data. Three healthy subjects participated in the study,
which was approved by the ethics committee of University ofCalifornia, San Francisco. The data were acquired using a 275channel CTF Omega 2000 whole-head MEG system (VSM
MedTech Ltd.). When collecting the data, subjects were laid
down with their eyes closed, and 60-s-long single continuous datasets were acquired with a 1200 Hz sampling rate. The

3126

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 11, NOVEMBER 2011

Fig. 6. Power spectra of the sensor data for three subjects. The power spectrum
is the mean power spectrum across all sensors and all trials. The two red lines
denote the beta range used for the coherence analysis.

single trial data were divided into 2-s-long trials with half overlap, resulting in a total of 60 trials. The power spectra of thesensor data for three-subject cases are shown in Fig. 6. Eachpower spectrum was computed using the mean power spectrum
across all sensors and all trials.
The voxel time courses were estimated using the narrow-band
adaptive spatial filter [12], [19] with the frequency band tuned
to the individual subjects’ beta band, approximately between 14
and 27 Hz. The two red lines in Fig. 6 denote the beta range foreach individual subject used for the narrow-band spatial filter.
The source reconstruction was performed using the NUTMEG
toolbox [20], and the results are shown in Fig. 7 for the threesubject cases. In these results, the reconstructed voxel power,
which is the squared mean of the voxel time course, was mapped,
and three orthogonal MR slices at the voxel with the maximum
intensity are shown. In all three-subject cases, a large blurredpeak exists in a middle, posterior region, and the peak maxima
are not necessarily located near primary motor areas.
We then performed source-coherence imaging by setting seedvoxels at the left motor area. To obtain seed voxels, we selected
all the voxels whose location was determined to be within theprecentral gyrus. These location labels were obtained by normalizing the subjects MRI to the Montreal Neurological Institutetemplate brain. We then computed the centroid for these voxels
and selected the nine nearest voxels to the centroid, resulting inten seed voxels in the left precentral gyrus for each subject. The
seed voxels for each individual subject are shown in the twoleftmost columns of Fig. 8.
We computed the magnitude and imaginary coherence between the seed voxels and the remaining voxels, resulting in ten
coherence images corresponding to each of the ten seed voxels.
The final coherence image was obtained by averaging those ten

Fig. 7. Results of source power reconstruction for the three-subject cases. The
results from the first, second, and third subjects are, respectively, shown in theupper, middle, and bottom rows. The reconstructed voxel power is computed
using the squared mean of the voxel time course. The normalized voxel power
is color-coded according to the color bar, and overlaid onto each subject’s MRI.
Three orthogonal MR slices at the voxel with the maximum intensity are shown.

coherence images. The statistical thresholding using the surrogate data method was applied to select voxels with statisticallysignificant values of coherence. Here, the number of surrogate
datasets was set to 200, and the level of statistical significance α
was set to 0.99. The results are shown in Fig. 8 in which the twocentral columns show the magnitude-coherence images, and thetwo rightmost columns show the imaginary-coherence images.
The results here show that the magnitude-coherence images
are dominated by the seed blur. That is, for all three subjects,
large blurred peaks are located near the seed voxels and noother peaks are seen in the magnitude-coherence images. On
the contrary, the imaginary-coherence images are able to detect
the right motor cortex and show the connectivity between theleft and right motor cortices in the beta band. The seed blur that
dominates the magnitude-coherence images is not seen in theimaginary-coherence images. These results effectively demonstrate that the artifactual seed blur does not exist when using
imaginary coherence as the connectivity metric.
V. DISCUSSION
In the analysis so far, we have ignored the interference termscS (t) and cT (t). Here, we discuss the influence of these interference terms, and we use the expressionsσS = σ̄S + CSandσ
T = σ̄T + CT .

(23)

Here, CS and CT are the spectra of the interference terms,
and σ̄S and σ̄S are the source spectra including only the leakage terms, which are defined such that σ̄S = σS + d1 σT andσ̄T = σT + d2 σS . We assume that the interference spectra are
independent from the source spectra, i.e., the relationships,
σ̄T CT∗  = 0, σ̄T CS∗  = 0, σ̄S CT∗  = 0, and σ̄S CS∗  = 0,
hold.

SEKIHARA et al.: REMOVAL OF SPURIOUS COHERENCE IN MEG SOURCE-SPACE COHERENCE ANALYSIS

3127

Fig. 8. Results of imaging beta-band coherence between the left and right primary motor cortices. The coherence images were obtained using the resting-state
MEG data, and results for three subjects are shown. In each of the three pair of columns, a coronal view is presented on the left, and a sagittal view on the right.
The two leftmost columns show the ten seed voxels in the precentral gyrus. The two middle columns show the magnitude-coherence images. The two rightmost
columns show the imaginary-coherence images. The value of magnitude coherence is color-coded according to the upper color bar and the value of imaginary
coherence to the lower color bar. The results from the first, second, and third subjects are, respectively, shown in the upper, middle, and bottom rows.

Then, using (23), the cross spectrum σT σ
S∗  is expressed asS∗  = σ̄T σ̄S∗  + CT CS∗ 

σT σ

(24)

where CT CS∗  indicates the spurious component due to the interference terms. The aforementioned equation indicates that the
spurious component caused by the interference terms is additive
to that caused from the leakage terms. Taking the imaginary partof both sides of (24), we deriveS∗ ) =  (σ̄T σ̄S∗ ) +  (CT CS∗ ) .
 (
σT σ

(25)

σT σ
S∗ ) =  (σ̄T σ̄S∗ ),
If CT CS∗  is real valued, we have  (and taking the imaginary part removes CT CS∗ .
To further investigate what scenarios make CT CS∗  realvalued, we must take the property of inverse algorithms into
consideration. We assume that the voxel time courses are
estimated using a spatial filter with real-valued weight vectors.2 Under this assumption, the interference time coursescS (t) and cT (t) are expressed as cS (t) = wT (r S )bI (t) andcT (t) = wT (r T )bI (t), where w(r S ) and w(r T ) are realvalued weight vectors of the spatial filter at the seed and thetarget locations, respectively, and bI (t) is the sensor recordings
due to the interference and sensor noise.
By Fourier transforming cS (t) and cT (t), we obtain CS (f ) =wT (r S )β I (f ) and CT (f ) = wT (r T )β I (f ), where β I (f ) is
the Fourier transform of bI (t). The cross spectrum CT CS∗  is,thus, expressed as
CT CS∗  = wT (r T )β I (f )β H
I (f )w(r S )(26)

where the superscript H indicates the Hermitian transpose (thecomplex conjugation plus matrix transpose).
2 Most of existing spatial filters including the narrow-band adaptive spatial
filter used in our experiments employ real-valued weight vectors.

Since the weight w(r) is assumed to be a real-valued vector,
in order for CT CS∗  to be real valued, the matrix β I (f )β H
I (f )
should be real valued. The diagonal elements of β I (f )β H
I (f )
are real valued. The (j, k)th off-diagonal element of this matrix
is given by
(j )

(k )∗

βI (f )βI∞(f ) =−∞(j )(j,k )RI(τ )e−2π if τ dτ

(27)

(k )

where βI (f ) and βI (f ) is the jth and kth elements of β I (f ),(j,k )
respectively. Here, RI (τ ) is the cross correlation defined as(j,k )RI∞(τ ) =
−∞

(j )

(j )

(k )

bI (t)bI (t + τ )dt

(28)

(k )

where bI (t) and bI (t) are the jth and the kth elements of(j,k )bI (t), respectively. Note that since bI (t) is real valued, RI (τ )is real valued.
(j,k )
In general, we can decompose RI (τ ) into(j,k )

RI(j,k )(τ ) = RE(j,k )(τ ) + RO

(τ )

(29)

(j,k )(j,k )
where RE (τ ) is an even function and RO (τ ) is an oddfunction. Since the Fourier transform of a real even function is
(j )
(k )∗
real valued, the real part of βI (f )βI (f ) is caused from the(j,k )
even-function component RE (τ ), and the imaginary part of
(j )
(k )∗
βI (f )βI (f ) is caused from the odd-function component(j,k )RO (τ ).
Assuming that the interference bI (t) is a stationary process,(j,k )(k ,j )
(j,k )
we have the relationship RI (−τ ) = RI (τ ) = RI (τ ).
(j,k )
Thus, RI (τ ) is an even function and CT CS∗  is real valued.
As a result, taking the imaginary part of CT CS∗  can effectively remove the spurious cross spectrum CT CS∗ . However,3128

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 11, NOVEMBER 2011

(j,k )

when bI (t) deviates from a stationary process, RO (τ ) may
become significant, and using the imaginary coherence may beless effective in removing CT CS∗ .
Imaginary coherence has already been applied to MEG
source-space connectivity analysis, and resting-state functional
connectivity mappings have been performed using imaginary
coherence for schizophrenia patients [21] and for brain tumor
patients [7], [22]. In these publications, however, the primary
scope has been on the clinical implications of the connectivity
results and no rigorous arguments have been presented regarding
the choice of imaginary coherence.
This paper has presented detailed theoretical analysis on theproperties of imaginary coherence in source-space coherence
imaging. Our analysis and experiments show that magnitudecoherence images are confounded by the spurious coherence
caused by the leakage properties of the inverse algorithm. The
spurious coherence is manifested as the seed blur, which is an
artifactual large peak around the seed voxel in the resulting coherence images. This paper proposes the use of the imaginary
coherence to remove this seed blur, and presents a detailed theoretical analysis that explains how the use of imaginary coherence
leads to the removal of the seed blur. The effectiveness of using imaginary coherence has been validated by our computer
simulation and experiments using resting-state MEG data.
It should be pointed out that the effectiveness of using imaginary coherence holds, regardless of the source imaging algorithm used to estimate voxel time courses. The use of theadaptive spatial filter in our computer simulation and the experiments was mainly due to the convenience in our softwaredevelopments. It should, however, be emphasized that the use ofimaginary coherence cannot retrieve source activities the source
imaging algorithm fails to reconstruct.
Our analysis shows that the estimated imaginary coherence
has the amplitude bias due to the leakage of an imaging algorithm employed. Recently, Brookes et al. [23] reported thatthe length of a time window for coherence computation also
affects the amplitude of estimated coherence. Further studies
are needed to examine whether other factors exist that also produce the amplitude bias, and whether such amplitude biases inimaginary coherence are really problematic for the estimationof functional connectivity with this metric.

ACKNOWLEDGMENT
We would like to thank S. Honma and A. Findlay for their
help in measuring the resting-state MEG data.

REFERENCES[1] Z. Liu, M. Fukunaga, J. A. de Zwart, and J. H. Duyn, “Large-scale spontaneous fluctuations and correlations in brain electrical activity observedwith magnetoencephalography,” NeuroImage, vol. 51, pp. 102–111, 2010.
[2] T. Montez, S. Poil, B. F. Jones, I. Manshanden, J. P. A. Verbunt, B. W. van
Dijk, A. B. Brussaard, A. van Ooyen, C. J. Stam, P. Scheltens, and
K. Linkenkaer-Hansen, “Altered temporal correlations in parietal alphaand prefrontal theta oscillations in early stage Alzheimer disease,” Proc.
Nat. Acad. Sci. U.S.A., vol. 106, pp. 1614–1619, 2009.

[3] A. P. Georgopoulos, E. Karageorgiou, A. C. Leuthold, S. M. Lewis, J.
K. Lynch, A. A. Alonso, Z. Aslam, A. F. Carpenter, A. Georgopoulos, L.
S. Hemmy, I. G. Koutlas, F. J. Langheim, J. R. McCarten, S. E. McPherson,
J. V. Pardo, P. J. Pardo, G. J. Parry, S. J. Rottunda, B. M. Segal, S. R. Sponheim, J. J. Stanwyck, M. Stephane, and J. J. Westermeyer, “Synchronousneural interactions assessed by magnetoencephalography: A functional
biomarker for brain disorders,” J. Neural Eng., vol. 4, pp. 349–355,
2007.
[4] J.-M. Schoffelen and J. Gross, “Source connectivity analysis with MEG
and EEG,” Human Brain Mapping, vol. 30, pp. 1857–1865, 2009.
[5] J. Gross, J. Kujara, M. Hämäläinen, L. Timmermann, A. Schnitzler, and
R. Salmelin, “Dynamic imaging of coherent sources: Studying neuralinteractions in the human brain,” Proc. Nat. Acad. Sci. U.S.A., vol. 98,
pp. 694–699, 2001.
[6] P. Belardinelli, L. Ciancetta, M. Staudt, V. Pizzella, A. Londei, N. B.
G. L. Romani, and C. Braun, “Cerebro-muscular and cerebro-cerebral
coherence in patients with pre- and perinatally acquired unilateral brain
lesions,” NeuroImage, vol. 37, pp. 1301–1314, 2007.
[7] A. G. Guggisberg, S. M. Honma, A. M. Findlay, S. S. Dalal, H. E. Kirsch,
M. S. Berger, and S. S. Nagarajan, “Mapping functional connectivity
in patients with brain lesions,” Ann. Neurol., vol. 63, pp. 193–203,
2007.
[8] D. Mantini, M. G. Perucci, C. D. Gratta, G. L. Romani, and M. Corbetta,“Electrophysiological signatures of resting state networks in the humanbrain,” Proc. Nat. Acad. Sci. U.S.A., vol. 104, pp. 13170–13175, 2007.
[9] L. Astolfi and F. Babiloni, Estimation of Cortical Connectivity in Humans.
New York: Morgan & Claypool Publishers, 2008.
[10] W. H. R. Miltner, C. Braun, M. Arnold, H. Witte, and E. Taub, “Coherenceof gamma-band EEG activity as a basis for associative learning,” Nature,
vol. 397, pp. 434–436, 1999.
[11] P. L. Nunez, R. Srinivasan, A. F. Westdorf, R. S. Wijesinghe, D. M. Tucker,
R. B. Silberstein, and P. J. Cadusch, “EEG coherency I: Statistics, reference electrode, volume conduction, Laplacians, cortical imaging, and interpretation at multiple scale,” Electroencephalogr. Clin. Neurophysiol.,vol. 103, pp. 499–515, 1997.
[12] K. Sekihara and S. S. Nagarajan, Adaptive Spatial Filters for Electromagnetic Brain Imaging. Berlin, Germany: Springer-Verlag, 2008.
[13] G. Nolte, O. B. L. Wheaton, Z. Mari, S. Vorbach, and M. Hallett, “Identifying true brain interaction from EEG data using the imaginary part ofcoherency,” Clin. Neurophysiol., vol. 115, pp. 2292–2307, 2004.
[14] J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. D. Farmer, “Testingfor nonlinearlity in time series: The method of surrogate data,” Physica
D, vol. 58, pp. 77–94, 1992.
[15] L. Faes, G. D. Pinna, A. Porta, R. Maestri, and G. Nollo, “Surrogate data
analysis for assessing the significance of the coherence function,” IEEE
Trans. Biomed. Eng., vol. 51, no. 7, pp. 1156–1166, Jul. 2004.
[16] Y. Benjamini and Y. Hochberg, “Controlling the false discovery rate: A
practical and powerful approach to multiple testing,” J. Roy. Statist. Soc.
B, vol. 57, pp. 289–300, 1995.
[17] T. E. Nichols and S. Hayasaka, “Controlling the familywise error rate infunctional neuroimaging: A comparative review,” Statist. Methods Med.
Res., vol. 12, pp. 419–446, 2003.
[18] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[19] S. S. Dalal, A. G. Guggisberg, E. Edwards, K. Sekihara, A. M. Findlay,
R. T. Canolty, M. S. Berger, R. T. Knight, N. M. Barbaro, H. E. Kirsch,
and S. S. Nagarajan, “Five-dimensional neuroimaging: Localization of
the time-frequency dynamics of cortical activity,” NeuroImage, vol. 40,
pp. 1686–1700, 2008.
[20] S. S. Dalal, J. M. Zumer, A. G. Guggisberg, M. Trumpis, D. D. E. Wong,
K. Sekihara, and S. S. Nagarajan, “MEG/EEG source reconstruction, statistical evaluation, and visualization with NUTMEG,” Comput. Intell.
Neurosci., art. 75893, 2011.
[21] L. B. Hinkley, J. P. Owen, M. Fisher, A. M. Findlay, S. Vinogradov,and S. S. Nagarajan, “Cognitive impairments in schizophrenia as assessed
through activation and connectivity measures of magnetoencephalography(MEG) data,” Frontiers Human Neurosci., vol. 3, art. 73, 2009.
[22] J. Martino, S. M. Honma, A. F. Findlay, A. G. Guggisberg, M. S. Berger,
H. E. Kirsch, and S. S. Nagarajan, “Resting functional connectivity inpatients with brain tumors in eloquent areas,” Ann. Neurol., vol. 69,
pp. 521–532, 2011.
[23] M. J. Brookes, J. R. Hale, J. M. Zumer, C. M. Stevenson, S. T. Francis,
G. R. Barnes, J. P. Owen, P. G. Morris, and S. S. Nagarajan, “Measuringfunctional connectivity using MEG: Methodology and comparison with
fcMRI,” NeuroImage, vol. 56, pp. 1082–1104, 2011.

SEKIHARA et al.: REMOVAL OF SPURIOUS COHERENCE IN MEG SOURCE-SPACE COHERENCE ANALYSIS

Kensuke Sekihara (F’09) received the M.S. and
Ph.D. degrees from the Tokyo Institute of Technology, Tokyo, Japan, in 1976 and 1987, respectively.
From 1976 to 2000, he was with Central ResearchLaboratory, Hitachi, Ltd., Tokyo, Japan. From 1985 to
1986, he was a Visiting Research Scientist at Stanford
University, Stanford, CA, and from 1991 to 1992, he
was at Basic Development, Siemens Medical Engineering, Erlangen, Germany. From 1996 to 2000, he
was with “Mind Articulation” research project sponsored by Japan Science and Technology Corporation.
He is currently a Professor in the Department of Systems Design and Engineering, Tokyo Metropolitan University, Tokyo, Japan. He is the author of the bookentitled AdaptiveSpatial Filters for Electromagnetic Brain Imaging (1st ed. New
York: Springer, 2008). His research interests include the electromagnetic source
imaging, reconstruction and inverse methods, multidimensional signal processing, and statistical signal processing, especially their application to functionalsource imaging.
Dr. Sekihara is a member of the IEEE Medicine and Biology Society and theIEEE Signal Processing Society. He is a Fellow of the International Society ofFunctional Source Imaging.

Julia P. Owen received the B.S.E. degree from The
Cooper Union, New York, NY, in 2003, and the Ph.D.
degree in bioengineering from the Joint Program in
Bioengineering, University of California, San Francisco (UCSF), and University of California, Berkeley,in 2011.
She did her graduate research with Dr. Nagarajan
at the Biomagnetic Imaging Laboratory, UCSF. Her
research interests include the development of source
localization algorithms for magnetoencephalography
and their use with functional connectivity metrics.

3129

Stephen Trisno received the B.S. degree in bioengineering from the University of California, Berkeley.
In the summer of 2010, he started his research
as an Undergraduate Researcher at the Biomagnetic
Imaging Laboratory, University of California, San
Francisco. He aspires to go to medical school after
completing his degree.

Srikantan S. Nagarajan (SM’09) received the M.S
and Ph.D. degrees in biomedical engineering from
Case Western Reserve University, Cleveland, OH.
He was a Postdoctoral Fellow at the Keck Center
for Integrative Neuroscience, University of California, San Francisco (UCSF). He is currently a Professor in the Department of Radiology and Biomedical Imaging and the Department of Bioengineering and Therapeutic Sciences, UCSF, and a member of the Joint Graduate Program in Bioengineering, UCSF/University of California, Berkeley. His
research interests include the areas of neural engineering and machine learning,
better understanding neural mechanisms of sensorimotor learning and speech
motor control, and developing algorithms for improved functional brain imaging, biomedical signal processing, and brain–computer interfaces.