University of Southern California Dissertations and Theses

Signal processing methods for interaction analysis of functional brain imaging data

(USC Thesis Other)

Signal processing methods for interaction analysis of functional brain imaging data

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SIGNAL PROCESSING METHODS FOR INTERACTION ANALYSIS OF
FUNCTIONAL BRAIN IMAGING DATA
by
Hua Brian Hui
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
May 2010
Copyright 2010 Hua Brian Hui
Dedication
to my family
ii
Acknowledgments
I am deeply indebted to my advisor, Prof. Richard Leahy, for his support, encourage-
ment, and invaluable advices. During the course of my Ph.D. studies, he guided me in
choosing interesting problems and in searching for solutions; he supported me finan-
cially and provided me a wonderful research environment; he has always been support-
ive to students and showed me what a great researcher and mentor should be, with his
enthusiasm and patience. I would like to take this opportunity to express my greatest
gratitude to him.
I also wish to thank the members of my guidance committee: Dr. Krishna Nayak,
Dr. Manbir Singh, Dr. Keith Jenkins, and Dr. Antonio Ortega for their suggestions and
valuable feedback.
Special thanks to Dr. Dimitrios Pantazis at USC, for invaluable help and support,
for being a considerate friend and a great collaborator. He has been like a mentor to me,
and I feel very lucky to know him.
I also want to Dr. Quanzheng Li at USC for being a wonderful friend and teaching
me how to focus and be effective.
I am grateful to Dr. John Mosher at the Cleveland Clinc for sharing his knowledge
and expertise on MEG. I also thank Dr. James Duncan at Yale, Dr. Jed Meltzer at
the NIH, and Dr. Karim Jerbi at Paris for providing data to conduct my research, and
offering valuable advice.
iii
I would like to thank my officemates Belma Dogdas, Juan Luis Poletti Soto and
Yu-teng Chang for being considerate and understandable, and also for our helpful late
afternoon discussions.
I also want to thank the other members of my research group: Sangtae Ahn, Francois
Tadel, Anand Joshi, Abhijit Chaudhari, Sanghee Cho, Anand Joshi, Zheng Li, Felix Dar-
vas, Sangeetha Somayajula, Syed Ashrafulla, Joyita Dutta, David Wheland, Yanguang
Lin, Ran Ren, Sergul Aydore, and Wentao Zhu for being part of big family. My graduate
studies would be different without them.
Most importantly, I want to thank my parents and my wife Dr. Janice Cheng for
providing unlimited support. Janice, thank you for all your help and sacrifices for me to
realize my dream!
iv
Contents
Dedication ii
Acknowledgments iii
List of Tables viii
List of Figures ix
Abstract xiii
Chapter 1: Introduction 1
Chapter 2: Overview of Interaction Analysis using Functional Imaging 7
2.1 Functional Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.1 Neural basis of interactions . . . . . . . . . . . . . . . . . . . . 8
2.1.2 Neural assemblies and large scale integration . . . . . . . . . . 8
2.1.3 Applications of functional interactions . . . . . . . . . . . . . . 9
2.1.4 Functional imaging for interaction analysis . . . . . . . . . . . 10
2.2 Noninvasive EEG and MEG . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Sources of EEG and MEG . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Instrumentation . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.3 Forward problem and inverse problem . . . . . . . . . . . . . . 12
2.2.4 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Intracranial Recordings . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.1 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Chapter 3: Mathematical Modeling of Functional Interaction 20
3.1 Different Categories of Interaction Measures . . . . . . . . . . . . . . . 20
3.2 Coherence and Correlation . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.1 Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
v
3.2.2 Coherence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.3 Imaginary part of the coherency . . . . . . . . . . . . . . . . . 25
3.3 Phase Synchrony . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 The MV AR Model Related Measures . . . . . . . . . . . . . . . . . . . 27
3.4.1 MV AR model . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4.2 Granger Causality . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.4.3 multivariate interactions measure . . . . . . . . . . . . . . . . . 29
3.4.4 Cross-talk effect on MV AR model based measures . . . . . . . 30
3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Chapter 4: Nulling Beamformer for Accurate Interaction Measurement 32
4.1 Standard Inverse Methods . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.1 Forward modeling . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.2 Regularized Minimum-norm Solutions . . . . . . . . . . . . . . 36
4.1.3 Linearly Constrained Minimum Variance (LCMV) Beamformer 37
4.2 A Nulling Beamformer to Avoid Signal Cancelation and cross-talk . . . 39
4.2.1 Nulling Beamformer . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.2 Eigenvector Constrained Nulling Beamformer . . . . . . . . . . 40
4.3 Tests for Significance for Interaction Networks . . . . . . . . . . . . . . 42
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4.1 Realistic simulation of interactions . . . . . . . . . . . . . . . . 45
4.4.2 MV AR analysis with focal sources . . . . . . . . . . . . . . . . 47
4.4.3 MV AR analysis with extended cortical sources . . . . . . . . . 49
4.4.4 Nulling beamformer with other interactions measures . . . . . . 52
4.4.5 Robustness to inaccurate patch size . . . . . . . . . . . . . . . . 53
4.4.6 Effect of noise . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
Chapter 5: Parametric phase synchrony and group analysis of working mem-
ory effect using intracranial EEG 62
5.1 Working Memory Experiment . . . . . . . . . . . . . . . . . . . . . . 63
5.1.1 Sternberg working-memory task . . . . . . . . . . . . . . . . . 64
5.2 Parametric Analysis of Phase Synchrony . . . . . . . . . . . . . . . . . 66
5.2.1 Circular data . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2.2 Representation of circular data . . . . . . . . . . . . . . . . . . 67
5.2.3 Summary statistics for circular data . . . . . . . . . . . . . . . 68
5.2.4 Circular data and phase synchrony . . . . . . . . . . . . . . . . 70
5.2.5 Probability model for circular data . . . . . . . . . . . . . . . . 71
5.2.6 Test for equality of concentration . . . . . . . . . . . . . . . . . 73
5.3 Image Registration for Group Analysis . . . . . . . . . . . . . . . . . . 78
5.3.1 Location of electrodes . . . . . . . . . . . . . . . . . . . . . . 78
vi
5.3.2 Surface parcellation . . . . . . . . . . . . . . . . . . . . . . . . 80
5.4 Group Analysis by Combiningp-values . . . . . . . . . . . . . . . . . 80
5.4.1 Methods for combiningp-values . . . . . . . . . . . . . . . . . 81
5.4.2 Combiningp-values for each region . . . . . . . . . . . . . . . 84
5.4.3 Controling false discovery rate for multiple comparisons . . . . 84
5.5 Network analysis using graph theory . . . . . . . . . . . . . . . . . . . 85
5.5.1 Representation of graph . . . . . . . . . . . . . . . . . . . . . . 86
5.5.2 Graph topology measure . . . . . . . . . . . . . . . . . . . . . 86
5.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.6.1 Frequencies with phase synchrony . . . . . . . . . . . . . . . . 88
5.6.2 Parametric test for equality of phase synchrony . . . . . . . . . 90
5.6.3 Image registration . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.6.4 Group analysis result . . . . . . . . . . . . . . . . . . . . . . . 93
5.6.5 Graph partition result . . . . . . . . . . . . . . . . . . . . . . . 95
5.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Chapter 6: Conclusion and Future Work 101
6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
References 105
vii
List of Tables
3.1 Different categories for interaction measure . . . . . . . . . . . . . . . 22
5.1 Different synchronized frequencies for all 13 subjects . . . . . . . . . . 89
viii
List of Figures
2.1 (Left) Whole-head Elekta Neuromag MEG system. (Right) MEG sen-
sors use low-temperature electronics cooled by liquid helium. . . . . . 13
4.1 Simulating cortical interactions with realistic oscillatory sources. (a) 5
dipole sources are selected on the cortical surface based on the locations
of surface-to-depth electrodes in a macaque monkey experiment. LFPs
from the implanted electrodes are assigned as time series to the sources.
(b) MEG data are simulated by projecting the LFP signals to the channel
topography of a 275-channel gradiometer MEG system. White Gaussian
noise is added to the data. (c) Different inverse methods are used to pro-
duce cortical activation maps. (d) The estimated source signals are eval-
uated for cortical interactions using MV AR models, coherence, phase
synchrony or other measures. (e) Significant interactions are found using
permutation tests while controlling for multiple comparisons. . . . . . 46
4.2 (Left) PDC interaction measures for five dipole sources using different
MEG inverse methods. Each subplot represents the PDC causal inter-
action between a pair of sources, where the subplot at i-th row and j-
th column represents the causal interaction from source j to source i.
The horizontal axis represents frequency from 0 to 50Hz, and the ver-
tical axis represents the value of PDC in the range of [0; 1]. (Right)
Significant PDC interactions found for five dipole sources with a per-
mutation test at an = 5% significance level corrected for multiple
comparisons. The arrows represent the directions of causal interactions,
with their radii indicating the maximal strengths of interactions over the
[0; 50] Hz frequency band. . . . . . . . . . . . . . . . . . . . . . . . . 48
ix
4.3 (Left) PDC interaction measures for cortical sources with extended sizes
of approximately 8 cm
2
. (Right) Significant PDC interactions found for
the five extended cortical sources with a permutation test at an = 5%
level corrected for multiple comparisons. . . . . . . . . . . . . . . . . . 50
4.4 (Top) PDC interactions computed using different inverse operators from
simulated MEG data (with noise). (Bottom) PDC interactions computed
using the same inverse operators but from simulated data without the
noise. The difference between them shows the levels of false interactions
caused by residual noise. . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.5 PDC interactions with confidence intervals. Interaction results are the
same as in figure 4.3a and 4.3d, but overlaid with [2; +2] confidence
lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.6 Coherence between all the possible pairs of extended cortical sources
estimated using different inverse methods. . . . . . . . . . . . . . . . . 53
4.7 Phase synchrony between sources 3 and 4, computed using different
inverse methods. The color scale shows the phase locking value between
the two sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.8 Simulation of cortical sources in the case where we have limited knowl-
edge of their true locations and extents. Left: true sources with approx-
imately 8 cm
2
spatial extents. Right: enlarged patches of around 29 to
45 cm
2
spatial extents. . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.9 Top: PDC interaction measures computed using (a) the true LPF sig-
nal, and reconstructed time series using the nulling beamformer with
(b) known patch sizes and (c) enlarged patch sizes. Bottom: significant
PDC interactions at an = 5% significance level after correcting for
multiple comparisons using a permutation test. . . . . . . . . . . . . . . 55
4.10 Comparison of the errors in estimated PDC interactions for different
inverse methods and different noise level. . . . . . . . . . . . . . . . . 57
x
4.11 PDC interaction computed from different noise types with SNR=3: (a)
the true interactions, (b) and (c) the interactions computed using nulling
beamformer from dipole sources simulation with white noise and brain
noise, (d), (e) and (f) from patch simulation with white noise using dif-
ferent inverse methods, (g), (h) and (i) from patch simulation with brain
noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.1 The time diagram of the Sternberg working-memory task . . . . . . . . 65
5.2 Circular data plot for an example of simulated von Mises distribution data 68
5.3 Linear histogram and circular histogram for circular data. The resultant
vector is denoted by red color in the right figure. . . . . . . . . . . . . . 69
5.4 The PDF of the von Mises distribution with different, where
0
= 0. . 72
5.5 The preoperative (left) and the postoperative (right) MRI of one subject.
Note the postoperative MRI is significantly deformed. . . . . . . . . . 79
5.6 The scheme for registration and surface parcellation . . . . . . . . . . . 79
5.7 Illustration of the narrow frequency (at 8Hz) phase locking at theta band
for subject 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.8 Histogram of phase difference between two electrodes for subject 1. The
red curve is the PDF of the fitted von Mises distribution . . . . . . . . . 90
5.9 Single-subject adjacency matrices of the memory load modulated inter-
action networks for two sided test (H
1
:
1
6=
4
) after FDR. White pix-
els indicate significant interactions between the electrodes correspond-
ing to the row and column; Black pixels indicate no significant interactions. 91
5.10 Adjacency matrix of single subject for the memory load modulated inter-
actions network for one sided test (H
1
:
1
<
4
) after FDR. . . . . . . 92
xi
5.11 Adjacency matrix of single subject for the memory load modulated inter-
actions network for one sided test (H
1
:
1
>
4
) after FDR. . . . . . . 93
5.12 All the electrodes (colored dots) from 12 subjects are mapped to com-
mon cortical surface using surface registration. Different color of the
electrodes indicate they are from different subjects. . . . . . . . . . . . 94
5.13 Adjacency matrix of parcellated regions of interest for the memory load
modulated interactions network for one sided test (H
1
:
1
<
4
) after
FDR at = 0:001. White pixels indicate significant interactions between
the electrodes corresponding to the row and column; Black pixels indi-
cate no significant interactions. . . . . . . . . . . . . . . . . . . . . . . 95
5.14 The memory task modulated significant network found on group level
using Fisher’s combinedp-value method. . . . . . . . . . . . . . . . . 96
5.15 Strength(degree) of each node is computed as the sum of all connection
weights for the node. . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.16 The strengths (degrees) of regions in a group level memory modulated
interaction network are indicated by the sizes of the spheres at the center
of each region of interest. . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.17 The connection of 8 individual nodes with the highest degree. Among
the 68 regions of interest, these 8 regions are those connected to other
regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
5.18 The community structure of the group level interaction network. 3 sub-
networks are found by maximizing the modularity. . . . . . . . . . . . 100
xii
Abstract
Modeling functional brain interaction networks using non-invasive EEG and MEG data
is more challenging than using intracranial recordings data. This is because most inter-
action measures are not robust to cross-talk (interference) between cortical regions,
which may arise due to the limited spatial resolution of EEG/MEG inverse procedures.
We describe a modified beamforming approach to accurately measure cortical interac-
tions from EEG/MEG data, designed to suppress cross-talk between cortical regions.
We estimate interaction measures from the output of the modified beamformer and
test for statistical significance using permutation tests. Since the underlying neuronal
sources and their interactions are unknown in real MEG data, we demonstrate the per-
formance of the proposed beamforming method in a novel simulation scheme, where
intracranial recordings from a macaque monkey are used as neural sources to simulate
realistic MEG signals. The advantage of this approach is that local field potentials are
more realistic representations of true neuronal sources than simulation models and there-
fore are more suitable to evaluate the performance of our nulling beamforming method.
Intracranial recordings have minimal cross-talk and therefore interactions can be
measured more reliably. However, performing group level studies is a challenging task
because of the sparsity and variable coverage of electrodes on each subjects’ brain.
We describe a set of group analysis procedures for intracranial EEG recordings, which
include registration of MRI volumes and cortical surfaces, and parcellation of anatom-
ical regions of interest. We use a parametric probability model to test for equality of
xiii
phase synchrony, and use Fisher’s combined p-value method to pool test results from
electrodes on individual subjects into the parcellated regions of interest. We apply our
group analysis procedure to intracranial EEG data recorded in a working memory exper-
iment and find an interaction network that is modulated by memory load.
xiv
Chapter 1
Introduction
Analyzing distributed functional interactions, or functional connectivity, between corti-
cal regions is a key issue in understanding how the human brain works. It will not only
improve our knowledge about the underlying mechanisms for normal brain function,
but also help us to better characterize the cause of pathological diseases. Functional
neuroimaging modalities such as positron emission tomography (PET), functional mag-
netic resonance imaging (fMRI), magnetoencephalography (MEG), electroencephalog-
raphy (EEG) and intracranial EEG (iEEG) provide us great opportunities to analyze the
cortical interactions. Among them, EEG, MEG and intractranial EEG have excellent
temporal resolution and offer great potential for better analysis of cortical interactions.
In this dissertation, we will present signal processing methods for functional interaction
analysis using these imaging modalities, and in particular using noninvasive EEG and
MEG to suppress linear cross-talk, and intracranial EEG recordings for group analysis
of interaction networks.
Noninvasive EEG and MEG externally measure the electromagnetic signals pro-
duced by electrical activity in the brain. Despite their relatively low spatial resolution,
their capability of monitoring neuronal activation at the millisecond time scale without
open-skull operation makes these modalities attractive for cortical interaction analysis.
However, until recently, most EEG and MEG literatures have focused on source local-
ization - the forward and inverse modeling of brain activations. Many inverse imaging
methods have been proposed to create cortical activation maps from the linearly mixed
MEG/EEG sensor measurements [Baillet et al., 2001, Lutkenhoner et al., 1996, Dar-
vas et al., 2004], which include dipole fitting [Scherg and V on Cramon, 1985, Mosher
1
et al., 1992], linearly constrained minimum variance (LCMV) beamforming [Van Veen
et al., 1997, Robinson and Vrba, 1999], and linearly regularized minimum-norm imag-
ing [Jeffs et al., 1987,Dale et al., 2000,Pascual-Marqui, 2002]. The research of EEG and
MEG inverse problem dated to 1980s, the research of functional interactions analysis,
however, has drawn little attention until recently.
As the research focus shifted to functional interaction analysis, many different meth-
ods have been proposed to model functional interactions, including spectral coher-
ence [Nunez et al., 1997], phase synchrony [Lachaux et al., 1999] and Granger causal-
ity [Brovelli et al., 2004]. They are generally used to characterize interactions from
depth electrode measurements [Brovelli et al., 2004, Ding et al., 2000, Kaminski et al.,
2001, Sehatpour et al., 2008], while in principle they extend to external EEG and MEG
measurements, the applications of most methods are limited by their sensitivity to the
cross-talk effect (or linear mixing). The broad sensitivity of MEG/EEG sensors intro-
duces a significant amount of linear mixing among the sensor measurements [Nunez
et al., 1997]. Since interaction measures generally are not able to distinguish true neu-
ronal interactions and the effect of instantaneous linear mixing, directly applying them
to raw MEG/EEG measurements would generate false positives.
Therefore, a more reasonable approach is to estimate the neural source signals using
existing inverse imaging methods and analyze the functional interactions from these
signal afterwards. Although this reduces the cross-talk effect in the raw sensor measure-
ments, most inverse methods primarily focus on creating accurate or unbiased source
localization and will not fully eliminating the cross-talk effect. For example, the lim-
ited resolution of minimum-norm imaging leaves substantial cross-talk between close
regions. Consequently, interaction methods do not generally perform well when used in
conjunction with linear imaging methods because of the linear cross-talk problem.
2
There are several possible approaches to solve the linear cross-talk problem. One is
to define new interaction measures which are less or not sensitive to linear cross-talk,
such as using the imaginary part of coherency [Nolte et al., 2004]. However, this method
still suffers from secondary cross-talk, which means that additional cortical sources can
leak into two cortical sites where we measure interaction and contribute to the imaginary
coherency, causing a mislocalization of cortical interactions. Another possible approach
is to use statistical test that can distinguish true interactions from linear cross-talk. It is
very difficult, however, to find such kind of statistical approach that can eliminate the
false interactions caused by cross-talk. The third approach is to use inverse methods
that are less sensitive to linear cross-talk, which is the approach we use to attack this
problem.
We propose the use of a nulling beamformer to accurately measure the cortical inter-
actions and reduce the cross-talk effect. The nulling beamformer is a modified version of
the LCMV beamformer approach [Van Veen et al., 1997], where additional nulling con-
straints are added to cancel signals from specific cortical locations. Furthermore, a set
of eigenvector constraints can be used to make the methods robust to specification of the
precise location and extent of the cortical sources of interest. To demonstrate the perfor-
mance of our method, we propose a novel simulation scheme: we simulate MEG cortical
sources with intracranially measured local field potentials (LFP) which are believed to
be better representations of true neuronal sources [Sutherling et al., 1988,Lachaux et al.,
2007] than strict simulation models. Because these signals are available before and after
linear mixing, we can compute interactions in both cases and evaluate the robustness of
our method to linear cross-talk. We measure the performance of nulling beamformer in
conjunction with different interaction measures and compare with those of the standard
LCMV beamformer and linear imaging methods.
3
Contrary to noninvasive EEG and MEG, intracranial EEG directly measures the elec-
trical field on the cortical surface and is much less sensitive to the linear cross-talk.
Therefore, it is more reliable for analysis of functional interactions. The second part of
the dissertation will be focused on interaction analysis using intracranial recordings.
We use a parametric paradigm to analyze the phase synchrony based on the circu-
lar statistics of the relative phases between signals [Mardia and Jupp, 2000, Stephens,
1969, Schelter et al., 2007]. The advantage of parametric analysis of phase syn-
chrony over the commonly used non-parametric methods [Jean-Philippe Lachaux,
1999,Le Van Quyen et al., 2001] is that we have the flexibility to conduct tests for signifi-
cance without computationally expensive resampling methods or surrogate methods. We
apply parametric tests for equality of phase synchrony between different experimental
conditions to intracranial EEG data of patients performing Sternberg working-memory
task.
The main challenge of applying intracranial EEG for interaction analysis is the
sparse spatial coverage of the implanted electrodes. Because intracranial recordings
generally are from patients undergoing monitoring for brain surgery, locations of elec-
trodes are chosen with clinical considerations, rather than basic research on functional
interaction networks. Therefore, intracranial recordings often do not measure signals
from interesting locations for interaction analysis. Moreover, because the locations of
electrodes are unique for each subject,it is not easy to pool information across subjects,
which makes group analysis using intracranial data very difficult.
We describe a new group analysis scheme for analyzing interactions from intracra-
nial recordings. We first map the electrodes from multiple subjects onto a common ref-
erence brain by registering the cortical surfaces extracted from MRI. We then parcellate
the cortical surface into a set of functionally congruent areas based on gyral and sulcal
cortical patterns. We combine the p-values of interaction for each pair of electrodes
4
between two cortical regions into a single summary p-value using Fisher’s combined
p-value approach. In this way, we reconstruct an interaction network for the regions of
interest on the group level.
We conduct network analysis on the group level network using graph theory based
methods. Graph theory has been introduced to structural and functional brain connec-
tivity studies recently [Gong et al., 2008, van den Heuvel et al., 2008, Iturria-Medina
et al., 2008, Reijneveld et al., 2007]. Graph theory can model the structure of complex
networks and allows us to see the organization of interactions on a large scale. Graph
theory also provides us the ability to partition the network into community structures,
which are groups of densely connected vertices with only sparse connections between
groups [Newman, 2006a, Newman, 2006b].
We apply the above analysis on the intracranial EEG data of a Sternberg working-
memory task and investigate memory load effects on the organization of the group net-
work.
This dissertation is organized as follows: Chapter 2 gives an overview of the electro-
physiological interaction networks, MEG/EEG and intracranial EEG, and their advan-
tages and limitations in analyzing the cortical interactions. Chapter 3 presents differ-
ent mathematical models for estimating functional interactions: coherence, phase syn-
chrony, Granger causality, and information theory based measures. We also discuss
their sensitivity on linear mixing (cross-talk). Chapter 4 begins with a description of
several inverse methods such as minimum-norm, spatial filtering, subspace correlation,
and matched filter techniques, to detect the neuronal current source configuration that
explains MEG measurements. It then introduces the use of nulling beamformer to accu-
rately estimate the interactions by suppressing the cross-talk. and compares the perfor-
mance of the nulling beamformer with other inverse methods by using a novel simula-
tion scheme. Chapter 5 first gives a brief introduction of working memory experiment,
5
followed by parametric analysis of phase synchrony. It then explains a novel group anal-
ysis scheme and applies it to the intracranial recordings of the aforementioned experi-
ment. Graph theory based network analysis is used on the group-level results and the
neuroscience aspect of the results is discussed in the end of the chapter. We conclude
the dissertation in Chapter 6 with several concluding remarks, and suggestions for future
work on functional interaction studies and their applications.
6
Chapter 2
Overview of Interaction Analysis using
Functional Imaging
Neurophysiological studies have shown evidences that the complicated brain activities,
such as sensorimotor, visual or cognitive functions, involve activation of multiple spe-
cialized brain regions. These regions not only are coactive, but also functionally con-
nected. The problem of interactions between these regions (also known as functional
connectivity or the binding problem) is essential for understanding how the brain works.
Functional neuroimaging modalities such as positron emission tomography (PET),
functional magnetic resonance imaging (fMRI), magnetoencephalography(MEG), elec-
troencephalography (EEG) and intracotical recording (intracranial EEG or local field
potential) have been used to assess functional interactions. These modalities have their
own advantages and disadvantages for analyzing interactions. Among them, MEG, EEG
and intracranial EEG have much higher temporal resolution and directly measure the
electrical brain activities, therefore have great potential in the functional interactions
studies.
In this chapter, we first describe the electrophysiological mechanism of neuronal
activity and functional interactions. We then present noninvasive functional imaging
modalities - MEG and EEG together with the electrophysiological origin, the instru-
mentation, the modeling and the clinic applications of MEG and EEG. We emphasize
their advantages in milliseconds temporal resolution and their limitation dues to the
intrinsic inverse problem. We also introduce the intracranial EEG recording and their
direct relationship with underlying functional interactions.
7
2.1 Functional Interactions
2.1.1 Neural basis of interactions
Neurons (also known as nerve cells) are the core components and basis of interactions.
Neurons are electrically excitable cells in the nervous system that process and transmit
information. A neuron is composed of the cell body (soma) containing the nucleus, a
dendritic tree receiving signals from other neurons and an axon conducting the nerve
signal. Neurons receive input on the cell body and dendritic tree, and transmit output
via the axon.
Neurons communicate with each other through chemical and electrical ”synapses”
with action potential, a propagating electrical signal that is generated by exploiting
the electrically excitable membrane of the neuron: When a pulse arrives at an axon
of a presynaptic cell, neurotransmitter molecules are released. These molecules bind
to receptors located on target cells, opening ion channels (mostly Na+, K+ and Cl-)
through the membrane. The resulting flow of charge causes an electrical current along
the interior of the postsynaptic cell, changing the postsynaptic potential (PSP). When
an excitatory PSP reaches the firing threshold at the axon hillock, it initiates an action
potential that travels along the axon with undiminished amplitude [Squire et al., 2008].
The human brain has a huge number of synapses. Each of these one hundred billion
neurons has on average 7,000 synaptic connections to other neurons. These synapses
forms the neural basis of function interactions on the brain.
2.1.2 Neural assemblies and large scale integration
Although communication between neurons is the basis of functional interactions, neuro-
physiological studies indicate that it is the large groups of neurons, or the neural assem-
blies that provide the functional elements of brain activity that execute the basic oper-
ations of informational processing [Fingelkurts et al., 2005]. These distributed neural
8
assemblies have emergent properties that do not exist at the level of individual neurons,
and their oscillations are more directly related to the specified aspects of neural tasks.
Furthermore, these assemblies become interdependent (functionally connected, or cou-
pled) when carrying out the complex brain operations. The functional interactions (also
known as the functional connectivity) is essential to understand the brain organization.
Studies also show the functional interactions not only exist among the nearby neu-
ral assemblies, they also appear among the farther apart assemblies, such as between
occipital and frontal lobes or across hemispheres. This is often referred as large scale
integration [Varela et al., 2001]. Since the farther apart cortical regions normally rep-
resent different functionalities, these large scale interactions are more interesting in a
sense that they represent the relationship between different neural functionalities.
2.1.3 Applications of functional interactions
Functional interaction has drawn great neuroscience attentions because it is closely
related to the perceptual binding problem and is also essential to the understanding of
cognitive functions.
The functional interaction analysis also has great pathological significance. Because
interactions represent the organization of the brain functions, it is natural that abnormal-
ity in interactions (lack of normal interactions or existence of abnormal interactions) has
been associated with neuropsychiatric disorders such as epilepsy, and Parkinson’s dis-
ease [Schnitzler and Gross, 2005]. For example, functional interaction analysis among
the cortical sources of epilepsy seizure will provide important information about the
prorogation of the epileptic activity across the brain, which will in turn help the neural
surgical planning.
9
2.1.4 Functional imaging for interaction analysis
To assess the functional interactions among the neural assemblies, different functional
imaging modalities have been used: positron emission tomography (PET), functional
magnetic resonance imaging (fMRI), magnetoencephalography (MEG), and electroen-
cephalography (EEG), intracranial recording. Each of them has different temporal and
spatial resolution and has its own advantage and disadvantage for analyzing the func-
tional interactions. In the following section, we will focus on the mostly used electrode
imaging modalities: intracranial recording, MEG and EEG.
2.2 Noninvasive EEG and MEG
MEG and EEG are complementary noninvasive techniques for measuring neuronal
activity in the human brain [Baillet et al., 2001]. MEG measures the magnetic induction
outside the head by electrical activity in neural cell assemblies, while EEG measures the
scalp electric potentials produced by electrical activity.
2.2.1 Sources of EEG and MEG
Both EEG and MEG signals come from the electrophysiological activation of neu-
ral assemblies we discussed in section 2.1.1. Although both the intracellular currents
(commonly called primary currents), and the extracellular currents (commonly called
secondary currents), contribute to magnetic fields outside the head and electric scalp
potentials, the currents associated with their large dendritic trunks are believed to be the
primary source of the MEG and EEG signal. In contrast, the temporal summation of
primary currents is not as effective in producing measurable external fields as the den-
dritic currents flowing in neighboring fibers, so that they are believed to contribute little
to MEG and EEG measurements [Baillet et al., 2001].
10
2.2.2 Instrumentation
EEG and MEG, as noninvasive ways to measure the neuronal activity in the brain, pro-
vide a safe and easy access to brain activation compared with the invasive intracranial
recordings. In addition, they require no radioactivity or injected contrast agents. There
are no known risks associated with MEG/EEG recordings so far.
EEG
The typical EEG scalp voltages are on the order of tens microvolts and can be mea-
sured using relatively low-cost scalp electrodes and high-impedance high-gain ampli-
fiers. Measurement of EEG comes from the electric potential differences between scalp
electrodes pairs measured by either the sensors glued to the skin or those in an elastic
cap for rapid attachment.
In the past, the electrodes transmitted the EEG signal through wires and the results
were typically recorded on a continuous roll of paper. Nowadays, they are more likely
to be stored in a computer’s memory or on a disk, sometimes even sent wirelessly to the
computer. Technological advances have make it possible for more expanded EEG mea-
surement and the combination of EEG measurement with video monitoring of seizure
activity in epileptic studies.
MEG
In contrast to the scalp potentials, the magnetic induction produced by the neural cur-
rent is extremely weak, typically 50-500 femtoTeslas, i.e. 10
9
or 10
8
times smaller
that the geomagnetic field of the earth [Hamalainen et al., 1993]. Therefore, the MEG
measurement is typical done in magnetically shield rooms and the Superconducting
QUantum Interference Device (SQUID) introduced in the late 1960s by James Zimmer-
man [Zimmerman et al., 1970]. The first measurement of human brain magnetic fields
11
using a SQUID magnetometer was carried out by David Cohen [Cohen, 1972] at the
Massachusetts Institute of Technology, and consisted of spontaneous alpha activity of a
healthy subject and abnormal brain activity of an epileptic patient.
The SQUID is a superconductive ring in which the magnetic flux is quantized in
very small units. When a constant biasing current is maintained in the SQUID device,
the measured voltage oscillates as the magnetic flux increases; one period of voltage
variation corresponds to an increase of one flux quantum. Counting the oscillations
allows one to evaluate the flux change that has occurred, and therefore detect magnetic
fields on the order of a few femtoTeslas. The sensitivity of the SQUID can be increased
to by attaching a coil of flux transformer placed in the vicinity of human head. These
flux transformers, according to its shape, can be configured as a first-order planar or axial
gradiometer, a second-order axial gradiometer, or a simple magnetometer [Hamalainen
et al., 1993]. The gradiometer configurations produce measurements proportional to the
spatial gradient of the magnetic field thus attenuate the noise sources distant from the
gradiometer and offering robustness to interference from distant magnetic field sources.
Modern MEG systems consist of whole-head sensor arrays with more than around
300 SQUID sensors placed in a liquid helium-filled dewar, with the flux-transformer
pickup coils surrounding a helmet structure (figure 2.1). They also include shielded
rooms made of successive layers of mu-metal, copper, and aluminum effectively atten-
uate high frequency disturbances.
2.2.3 Forward problem and inverse problem
Since the neuronal currents are not measured directly in EEG and MEG as in intracranial
recording, they have to been estimated based on the MEG or EEG measurement signals
from the sensor array. This estimation problem is referred as inverse problem. However,
before we could make such an estimation, we have to first understand how the scalp
12
Figure 2.1: (Left) Whole-head Elekta Neuromag MEG system. (Right) MEG sensors
use low-temperature electronics cooled by liquid helium.
potentials or external magnetic fields are generated for a specific set of neural current
sources, which is often called forward problem.
Forward modeling
Because the useful frequency spectrum for electrophysiological signals is typically
below 100Hz, the physics of EEG and MEG can be described with the quasi-static
approximation of Maxwell’s equations. The propagation of electromagnetic fields inside
the head is estimated based on the conductivities of the scalp, skull, gray and white mat-
ter, cerebrospinal fluid and other tissue types. Although head models that consist of
a set of nested concentric spheres with isotropic and homogeneous conductivities have
closed form solutions and work surprisingly well, we can get more accurate solutions by
using realistic head models based on anatomical information from high-resolution mag-
netic resonance (MR) or x-ray computed tomography (CT) volumetric images. These
13
more realistic models can be solved using numerical methods such as boundary element
methods (BEM) or finite element methods (FEM) [Darvas et al., 2004].
Inverse methods
To make inferences about the brain activity that gives rise to a set of MEG or EEG data,
we must solve the inverse problem, i.e. find a neuronal current source configuration that
explains the MEG or EEG measurements. Inverse methods for MEG/EEG typically fall
into two categories: dipole fitting methods and imaging methods.
The dipole fitting methods [Baillet et al., 2001], use a few equivalent current dipoles
to approximate the electrical current source and estimate the unknown location and
moment of the current dipoles based on MEG or EEG signal using a least squares
fitting or more recently using scanning methods. Scanning methods are based on the
dipole model, but involve scanning a source volume or surface and detecting sources
at those positions at which the scan metric produces a local peak [Baillet et al., 2001].
Examples of these methods include the MUSIC (MUltiple Signal Classification) algo-
rithm [Mosher et al., 1992] and the LCMV (Linearly Constrained Minimum Variance)
beamformer [Van Veen et al., 1997].
The imaging methods are based on the assumption that the primary sources are intra-
cellular currents in the dendritic trunks of cortical pyramidal neurons that are aligned
normally to the cortical surface. Consequently, a tessellated representation of the cere-
bral cortex can be extracted from a coregistered MR image and the inverse problem
is solved for a current dipole located at each vertex of the surface. Since the position
and orientation of the dipoles are fixed in this case, image reconstruction becomes a
linear problem and can be solved using standard techniques. The typical imaging recon-
struction methods includes regularized minimum-norm inverse methods with different
normalization factors.
14
We will discuss the inverse methods in more details in chapter 4
2.2.4 Resolution
MEG and EEG have temporal resolutions below 100 ms, which allow us to explore
the dynamics and functional interactions of the brain process at neural assembly level.
However, due to the underdetermined nature of the inverse problem in MEG and MEG
and the relatively small number of sensor measurement, their spatial resolution are lower
than intracranial recordings, PET and fMRI. It goes from a few millimeters to a couple
of centimeters, depending on the experiment. Focal cortical sources tend to spread over
multiple cortical sulci and gyri.
Because of this low spatial resolution, the estimated signal from various inverse
methods at any cortical location will be a linear mixing of the true neural signal at this
location with some signals from the rest of the brain. This is the common cross-talk
problem of the EEG and MEG inverse methods, we will discuss it more in the later
chapters.
Comparing MEG with EEG, the scalp electrical potentials are often affected by var-
ious inhomogeneities in the head where as the magnetic fields are mainly produced by
currents in intracranial space and less sensitive to the inhomogeneous conductivity of
skull and scalp. Therefore, MEG has better spatial accuracy than EEG in determining
the locations of source activity. However, because the magnetic field measured in MEG
is very weak compared with ambient noise, it is more sensitive to outside disturbances
include power-line fields, electronic devices, elevators, and radio-frequency waves.
2.2.5 Applications
Applications for MEG and EEG include both basic and clinical research. EEG and
MEG has been widely used for functional brain mapping. Evoked response fields have
15
been used to identify somatosensory, motor, and vision related activity. Several stud-
ies [Bowyer et al., 2005] have also used MEG and EEG to localize cortical activity
related to language, memory and other cognitive tasks using either equivalent current
dipoles or distributed cortical imaging. A wide range of signal processing techniques
including image modeling and reconstruction, blind source separation, and chaos the-
ory are under investigation to reveal complex cognitive processes such as attention and
working memory.
One of the most important clinical applications for MEG and EEG is the detection,
classification and localization of abnormal neuronal activities in epilepsy patients. EEG,
and recently MEG have been successfully used to localize three different spontaneous
interictal signal components: epileptic spikes, slow-wave activity, and fast-wave activ-
ity [Maudgil, 2003]. Neurosurgical planning of medically intractable epilepsy often
includes the identification of epileptogenic lesions with EEG and MEG [Stefan et al.,
2003, Ossadtchi et al., 2004].
Generally, EEG and MEG are regarded as two complementary noninvasive imaging
techniques. While spatial resolution is significantly lower than that of PET and fMRI,
the ability to monitor neuronal activation at the millisecond time scale using these two
modalities opens a unique window on human brain. Coupled with new functional inter-
action analysis tools, MEG will lead to new insights into the functions of the human
brain with applications in both clinical and cognitive neuroscience.
2.3 Intracranial Recordings
Intracranial recordings (also known as intracranial electroencephalography or iEEG as
opposed to scalp EEG) are recorded directly by electrodes implanted inside the brain.
Intracranial recordings usually measure the local field potentials (LFP) and occasionally
spikes [Lachaux et al., 2003].
16
Intracranial recordings by far provide the most robust measure of functional inter-
action since they are the direct summation (or spatial average) of currents from large
numbers of neurons in the local neighborhood. They are minimally affected by vol-
ume conduction (the interference from the rest of the brain), muscle and eye movements
artifacts. They are directly related to the underlying electrophysiological activation and
suitable for interactions.
Various studies have used LFP from cat and primate [Brovelli et al., 2004] [Wors-
ley et al., 2004] to assess the interactions. Human intracranial recordings are mostly
acquired from patients with medically intractable epilepsy, Parkinson’s disease or brain
tumors. Several studies have also used the human intracranial recordings for interaction
analysis [Kahana, 2006] [Zaveri et al., 1999].
2.3.1 Resolution
Intracranial recordings have both high spatial and high temporal resolution and have the
”millimeter-millisecond” resolution.
Similar with noninvasive EEG and MEG, the temporal resolution of intracranial
recordings theoretically is that of the electrophysiological phenomena it measures, and
it can even goes to sub-millisecond. And it is only limited by the sampling rate of the
acquiring device in practice. This is very important characteristics for analyzing the
functional interactions since we are more interested in the interdependence of the ”tem-
poral patterns” or the oscillations. Higher temporal resolution will provide us greater
opportunities to discover the interactions.
The spatial resolution of intracranial recordings can goes from the micro-scale res-
olution of single or multi-unit from the micro-electrodes, to the meso-scale resolution
of the local field potentials (LFP) from surface strips or grids of electrode or stereo-
tactically placed depth electrodes. Since we are interested in the signals from neural
17
assemblies in large-scale integrations, meso-scale LFP is the mostly used one. Stud-
ies has confirmed that the spatial resolution of the functional imaging is also important
for analyzing the interaction because it determines how accurate the signals we mea-
sure will represent the underlying electrophysiological activation. As we will discuss
in chapter 3, spatial resolution represents how much the computed interaction will be
affected by the interference (cross-talk or volume conduction): the worse the resolution,
the more cross-talk from other neighborhood neural assemblies will be contained in the
measured signals. Fortunately, in intracranial recording such cross-effect is very small
and can be neglected most time for functional interaction analysis.
2.3.2 Limitations
Despite the advantages of intracranial recordings for functional interaction analysis,
there are also some limitations. Human functional interaction analysis using intracra-
nial recordings always have ambiguity from pathology. Since the human intracranial
recordings were mostly recorded from patients with epilepsy or Parkinson’s disease. It
is possible the interactions in the brains of these patients are different from the normal
human subjects, either due to the pathology or the medication. Therefore, we have to
be careful when determining what kind of the task to be included in the experiments in
order to be independent of pathology and when interpreting the estimated interactions
from intracranial recordings.
Another limitation of intracranial recordings for interactions is the limitation of spa-
tial coverage in the brain. Because in all studies, the surface strips or grids of elec-
trodes are placed only on certain areas on the cortex. It is almost impossible to implant
electrodes to cover the whole brain surface simultaneously. Therefore, this limitation of
spatial coverage will prevent us from using intracranial recordings to analyze long-range
interactions such as between motor areas and frontal areas, or cross hemisphere. And
18
we also need to carefully select the area to implant the electrodes based on the specific
tasks to perform.
2.4 Summary
In this chapter, we have discussed the neural basis of functional interactions, and intro-
duced different functional imaging methods for assessing the functional interactions.
We have shown the advantage of intracranial recordings as direct measurement of the
neural activity and their limitations. We discussed the noninvasive measurement of the
neural activity - MEG and EEG, also the limit spatial resolution due to the inverse prob-
lem and the resulting cross-talk effects. We will discuss what these effects would mean
for functional interactions in the next chapter.
19
Chapter 3
Mathematical Modeling of Functional
Interaction
We have discussed the neural basis of functional interactions and functional imaging.
To quantitatively analyze the interactions, many different methods have been proposed
to measure their strengths and/or directions, including spectral coherence [Nunez et al.,
1997], phase synchronization [Lachaux et al., 1999], Granger causality [Brovelli et al.,
2004], etc. In this chapter, we will study some commonly used interaction measures and
their capability on handling linear cross-talk.
3.1 Different Categories of Interaction Measures
As the functional interaction analysis attracts more interests, many interactions mea-
sures have been proposed for different purposes and targeted at different applications.
For better understanding of the interactions, we categorize them into based on different
criteria:
Linear and nonlinear
Interaction measures such as correlation, coherence [Nunez et al., 1997] and Granger
causality [Brovelli et al., 2004] capture the linear relationship between cortical sig-
nals, but can not detect more complicated non-linear relationship. Whereas some
interaction measures capture both linear and nonlinear relationship such as phase syn-
chrony [Lachaux et al., 1999]. Although a small portion of nonlinear relationship can
be approximated by linear interactions, most of chaotic nonlinear relationship will be
missed by using linear measures.
20
Parametric and nonparametric
The measures discussed in this chapter are essentially different model for interac-
tions, they can be either complicated mathematical models involving different numbers
of parameters, or more general nonparametric approach. Parametric models such as
Granger causality and other MV AR model based measures have more detection power
when the assumptions are satisfied. However, they are also more sensitive to noise in
signals and produce erroneous results when assumptions are violated. Nonparametric
measures are simpler and more robust compared with parametric models. Example of
nonparametric interaction measures are coherence, standard form of phase synchrony,
information based measure, etc. In chapter 5, we use a parametric model for phase
synchrony under the assumption of the statistical distribution of relative phase follows
the von Mises distribution. We will see the parametric model has more flexibility to
compare phase synchrony between different experimental conditions.
Pairwise and multivariate
When analyzing functional interactions, we can either compute the interaction pair-
wisely between two signals or among multiple sites. The pairwise measures are easy to
define and compute, and most commonly used interaction measures, such as coherence,
phase synchrony and Granger causality, are pairwise. Recently, multivariate measures
which quantify the interactions among multiple sites become popular. Although they
are more complicated and computationally more expensive to compute, they have the
advantage of distinguishing true interaction between two signals apart from the inter-
actions caused by a third sources interacting with both. MV AR model based measures
such as Directed Transfer Function (DTF) and Partial Directed Coherence (PDC) are
multivariate measures. Because they are also parametric model, the number of the mul-
tivariate signals can not be very large due to the limited degrees of freedom.
21
Interaction Linear or Nonlinear Parametric or Nonparametric Pairwise or Multivariate
Coherence Linear Nonparametric Pairwise
Phase Synchrony Nonliear Nonparametric Pairwise
Granger Causality Linear Parametric Pairwise
DTF/PDC Linear Parametric Multivariate
Information Theory Nonlinear Nonparametric Pairwise
Table 3.1: Different categories for interaction measure
We summarize the interaction measures and their categories in table 3.1.
3.2 Coherence and Correlation
3.2.1 Correlation
The simplest way to measure relationship between two signals is the linear correlation.
If signal x(t) increases as signal y(t) increases, we call x(t) and y(t) are positively
correlated. If signal x(t) decreases as signal y(t) increase, we call x(t) and y(t) are
negatively correlated. We can compute instantaneous correlation by
corr
xy
(t) = Ef[x(t)
x
(t)][y(t)
y
(t)]

g (3.1)
Note that this value depends on the power of the signal. To make it independent of the
power, we can define the correlation coefficients,

xy
(t) =
Ef[x(t)
x
(t)][y(t)
y
(t)]

g
q
E

jx(t)
x
(t)j
2

E

jy(t)
y
(t)j
2

(3.2)
22
If the signals are stationary, the ergodicity allows us to compute the correlation coef-
ficent using the average over time:

xy
=
T
P
t=0
f[x(t)
x
(t)][y(t)
y
(t)]

g
s
T
P
t=0
jx(t)
x
(t)j
2
T
P
u=0
jy(u)
y
(u)j
2
(3.3)
This instantaneous correlation (coefficient) is used in resting state network of fMRI
to measure the relationship between signals from different voxels at 0.1Hz [Biswal et al.,
1995] [Fox and Raichle, 2007] [Scott J. Peltier, 2003] [Shehzad et al., 2009]. Because
of the low frequency, any causal effect will be ignored, and the relationship will appear
as instantaneous.
However, when dealing with signal at 1-50Hz in electrophysiogical signal such as
MEG and EEG, we have to acknowledge the fact that the brain activities do not happen
simultaneously and there are delays between signals at different areas of the brain. A
more general definition is to compute the correlation with different delays using the
cross-correlation:
R
xy
(t;) = Ef[x(t)
x
(t)][y(t)
y
(t)]

g (3.4)
For stationary signal, it will be independent of timet,
R
xy
() =
T
X
t=0
f[x(t)
x
(t)][y(t)
y
(t)]

g (3.5)
We can see the instantaneous correlation defined above is actuallyr() when = 0.
If we are interested in the frequency band at which the signals are correlated, we can
take Fourier transform onr
xy
() and compute the cross-spectral densityS
xy
:
S
xy
(f) =FfR
xy
()g =
Z
+1
1
R
xy
()e
j2f
d (3.6)
23
Similarily, we can define the auto-correlation and auto-spectral density:
S
xx
(f) =
Z
+1
1
R
xx
()e
j2f
d
where
R
xx
() =
T
X
t=0
x(t)x

(t)
3.2.2 Coherence
Coherence is a measure of the linear relationship between two signals at a specific
frequency [Carter, 1987, Carter et al., 1973, Gardner, 1992]. The complex coherence
(sometimes referred to as coherency) is defined as the cross-spectrum normalized by the
square root of the power spectra of both signals.
C
xy
(f) =
S
xy
(f)
p
S
xx
(f)
p
S
yy
(f)
(3.7)
where S
xy
(f) is the cross-spectral density function between signal x(t) and y(t), and
S
xx
(f) and S
yy
(f) are their power spectral densities. In practice, the magnitude of
the complex coherencejC
xy
(f)j is generally used for coherence. Some studies use the
square of magnitudejC
xy
(f)j
2
, which is referred to as Magnitude Squared Coherence
(MSC).
The magnitude of coherence takes values between 0 and 1, with 0 indicating that the
two signals are perfectly uncorrelated at frequencyf. A value of unity indicates perfect
correlation. Since coherence measure the linear dependence between two signals, it is
less sensitive to the nonlinear interdependence. [Gardner, 1992] has shown that nonlin-
ear interdependence such asy(t) =x(t)
2
could produce low degree of coherence. As a
measure for linear interactions, it is not surprising that coherence could not distinguish
24
the true interdependence from the linear mixing cross-talk. Suppose the estimate of
signalx(t) is contaminated with the cross-talk fromy(t), i.e.
^
x(t) =x(t) +y(t)
^
y(t) =y(t) (3.8)
Then the resulting coherence computed from
^
x(t) and
^
y(t) will contain an additional
term due to the cross-talk term in
^
x(t).
^
C
xy
(f) =C
xy
(f) +
p
S
yy
(f)
p
S
xx
(f)
(3.9)
Therefore, coherence might show false interactions due to the cross-talk effect even
when the real signals are totally independent from each other. One has to be very care-
ful when applying coherence measure to assess the interactions using the functional
imaging techniques contains cross-talk. For example,although event related coherence
in the sensor domain of EEG and MEG has been widely reported [], it is possible some
of the results are simply effect of cross-talk.
3.2.3 Imaginary part of the coherency
As indicated in Eq.(3.7), phase shift between the two signals can be represented by
phase of the complex coherence. Linear mixing is instantaneous and can produce no
phase shift, therefore in the presence of only linear cross-talk, we would expect the
coherence to be entirely real, just as the
p
Syy (f)
p
Sxx(f)
term in Eq.(3.9) due to the cross-
talk. Conversely, true interactions that are distinguishable from cross-talk produce an
imaginary component in coherence. For this reason, [Nolte et al., 2004] proposed to use
the imaginary part of coherency as a robust measure of cortical interaction to avoid the
false interactions found by coherence due to the cross-talk effect.
25
Imaginary of complex coherence reduces the cross-talk effect in coherence measure.
However, it has some limitations: First, it only shows partial interactions when the real
interactions is also instantaneous. Such instantaneous true interactions will be ignored
as the real part of coherency. Second, this method still suffers from secondary cross-
talk when there are additional correlated sources correlated with both signals of interest.
Therefore, in practice, even this measure can be affected by cross-talk because of the
effects of additional sources, noise and finite sample size.
3.3 Phase Synchrony
Phase synchrony, also known as phase locking, measures the frequency-specific tran-
sient phase-locking between two signals with a ratio of n : m [Lachaux et al., 1999]
[Le Van Quyen et al., 2001]:
n
1
(t)m
2
(t) = constant (3.10)
wheren,m are integers indicating the ratios of possible frequency locking. In practice,
we are mostly interested in the case thatn =m = 1. In this case we check the constancy
of phase difference within a limited window T. However, because of the background
noise existing in the brain signals, we would never be able to get perfectly constant
phase difference. We have to check the phase difference for an approximate constant:
=
1
(t)
2
(t) constant (3.11)
In another word, the degree of phase locking is estimated as peaks in the distribution
of phase difference in the unit circle over a window T. Therefore, the estimation of
phase locking generally involves two steps: (1) compute the instantaneous phase of
each signal, (2) estimate the degree of the constancy of the phase difference.
26
(1) Computation of instantaneous phase
The instantaneous phase (f;t) of a signal can be computed either using the
Hilbert transform or the complex Morlet wavelet transform [Lachaux et al., 1999]
[Le Van Quyen et al., 2001].
(2) Estimation of peakiness of the phase differences
To compute the peakiness of the phase differences, we average the complex number
e
j(
1

2
)
over the multiple trials.
PLV (f;t) =
1
N
trial

N
trial
X
trial=1
e
j(
1;trial
(f;t)
2;trial
(f;t))

(3.12)
This trial-average of the phase difference at each point (f;t) is referred to as phase
locking value (PLV). As shown in the appendix of [Lachaux et al., 1999], phase syn-
chrony is also sensitive to the cross-talk effect: two sources with no relationship between
them can still show non-zero phase synchrony if the signal from one source leaks to the
other.
3.4 The MV AR Model Related Measures
3.4.1 MV AR model
A new group of interaction measures have been proposed based on Multivariate Autore-
gressive (MV AR) model with ability to estimate the direction as well as the strength
of the interactions. The MV AR model is a straightforward extension of the univari-
ate autoregressive time series model. It incorporates interactions among multiple time
series. With x(t) = [x
1
(t);:::;x
p
(t)]
T
representing a p-dimensional time series, the
MV AR model can be expressed as:
x(t) + A(1)x(t 1) + + A(m)x(tm) = n(t) (3.13)
27
where n(t) = [n
1
(t);n
2
(t);:::;n
p
(t)]
T
is a zero mean uncorrelated noise vector with
covariance matrix . A(k)’s are the pp coefficient matrices. The coefficients for
a particular data set are estimated by solving a multivariate version of the Yule-Walker
equations [Kaminski et al., 2001].
After fitting time series into the MV AR model, we can performe multivariate spectral
analysis by transforming (3.13) into the frequency domain [Kus et al., 2004]:
x(f) = A
1
(f)n(f) = H(f)n(f) (3.14)
Where vector x(f) represents the Fourier transform of the time series and H(f) is a
pp transfer function matrix,
H(f) [A(f)]
1

"
m
X
k=0
A(k)e
2ikf
#
1
(3.15)
Using this general form, a number of different functional connectivity metrics can
be computed: Granger Causality [Brovelli et al., 2004], the Directed Transfer Function
[Kaminski et al., 2001] , and Partial Directed Coherence [Kus et al., 2004].
3.4.2 Granger Causality
Granger causality defines the pairwise causal influence from time series X to time series
Y if the variance of the autoregressive prediction error of Y is reduced by inclusion of
past measurements from X. [Granger, 1969]. Geweke [Geweke, 1982] extended this
idea by formulating Granger causality as frequency decomposition. Unlike coherence
and phase synchrony, Granger casualty can provide us the information about the direc-
tion of the causal influence.
28
3.4.3 multivariate interactions measure
The pairwise interaction metrics we mentioned above, although effective for cortical
interactions of different nature, all have a problem of ambiguity for indirect interactions
caused by additional sources. The estimated significant interaction between X and Y
can simply caused by the simultaneous interaction between X and Z, and between Y
and Z. Therefore, we need to use multivariate interaction metrics to avoid this ambiguity
problem.
The Directed transfer function (DTF) [Kaminski et al., 2001] and Partial Directed
Coherence (PDC) [Kus et al., 2004] are extensions of pairwise Granger causality for
multivariate time series using MV AR model. Both measures are frequency dependent
matrix with the (i;j)-th entry of the matrix representing a measure of the interaction
from sourcej to sourcei. They are defined using the elements of H(f) or its inverse
A(f) as the following:
DTF
2
ij
(f) =
jH
ij
(f)j
2
p
X
k=1
jH
ik
(f)j
2
PDC
2
ij
(f) =
A
2
ij
(f)
p
X
k=1
A
2
kj
(f)
(3.16)
Both of DTF and PDC are normalized to [0; 1] range with different normalized fac-
tors (the denominator) [Kus et al., 2004]: DTF is normalized to all the inflows of the
activities to the destination channeli; PDC is normalized to all the outflows of the activ-
ities from the source channelj. Both of them give similar results except that DTF could
not distinguish the indirect interactions from the direct ones [Kus et al., 2004].
29
3.4.4 Cross-talk effect on MV AR model based measures
Because linear interactions are concerned in MV AR model based analysis, the linear
cross-talk present in the signals are regarded as cortical interactions and produces erro-
neous results. Assuming the estimated sources ~ x(t) are the blurred version of the true
sources x(t) by the cross-talk through U:
~ x(t) = Ux(t) (3.17)
U is the resolution matrix defined in [de Peralta-Menendez and Gonzalez-Andino,
1998]. The MV AR model for the estimated sources ~ x(t) becomes:
~ x(t) +
~
A(1)~ x(t 1) + +
~
A(m)~ x(tm) = ~ n(t) (3.18)
The coefficient matrices
~
A(k) relate to those of the MV AR model for the true sources
by
~
A(k) = UA(k)U
1
(3.19)
Substitute (3.19) for the transfer function matrix
~
H(f)
~
H(f) =

m
X
k=0
UA(k)U
1
e
2ikf
!
1
= UH(f)U
1
(3.20)
Where H(f) is the transfer function matrix computed from the true sources time series
x(t).
The result in (3.20) shows that the transfer function matrix is contaminated by the
resolution matrix U. Therefore the DTF defined using
~
H(f) and the PDC using
~
A(f)
in (3.16) will be blurred because of the cross-talk.
30
3.5 Summary
In this chapter, we have studied several commonly used interaction measures, and cat-
egorized them according to different criteria. We noted that all of the interaction mea-
sures listed here are sensitive to cross-talk caused by EEG/MEG inverse operation, and
will produce significant false interactions if applied directly. In next chapter, we will
describe a novel inverse method to suppress the cross-talk effect in order to estimate
true interactions.
31
Chapter 4
Nulling Beamformer for Accurate
Interaction Measurement
In chapter 3, we have discussed several different types of mathematical models for mea-
suring the functional interactions. These measures have been widely used to character-
ize interactions from depth electrode measurements [Brovelli et al., 2004, Ding et al.,
2000, Kaminski et al., 2001, Sehatpour et al., 2008]. While in principle they extend
to EEG and MEG measurements, their utility is limited by their sensitivity to the cross-
talk effect (or linear mixing). The broad spatial sensitivity of MEG/EEG sensors [Nunez
et al., 1997] introduces a considerable amount of linear mixing among the sensor mea-
surements. Since interaction measures are limited in their ability to distinguish between
true neuronal interactions and the effect of instantaneous linear mixing, applying them
directly to raw MEG/EEG measurements would generate false positives. To better ana-
lyze the cortical interaction, the neuronal sources have to be identified and the current
densities at the sources have to be estimated [Hadjipapas et al., 2005] from the inverse
procedures.
Many inverse imaging methods have been proposed to create cortical activa-
tion maps from the linearly mixed MEG/EEG sensor measurements [Baillet et al.,
2001,Lutkenhoner et al., 1996,Darvas et al., 2004]. These include dipole fitting [Scherg
and V on Cramon, 1985, Mosher et al., 1992], linearly constrained minimum variance
(LCMV) beamforming [Spencer et al., 1992, Van Veen et al., 1997, Robinson and Vrba,
1999], and minimum norm imaging [Jeffs et al., 1987,Dale et al., 2000,Pascual-Marqui,
2002]. Although the estimation of neural source signals using inverse imaging reduces
32
the cross-talk effect present in the raw sensor measurements, most inverse methods pri-
marily focus on creating accurate or unbiased source localization instead of fully elim-
inating the cross-talk effect. For example, the limited resolution of minimum-norm
imaging leaves substantial cross-talk between nearby regions, as we demonstrate later
in the paper.
An LCMV beamformer has higher resolution than minimum-norm imaging when
the cortical sources are focal, which makes them more suitable for assessing interactions
[Hadjipapas et al., 2005]. However, the underlying assumption of adaptive beamformer
is that the neural sources are incoherent. When signals exhibit coherent behavior, the
beamformer will fail to form deep nulls at the locations of other coherent sources due to
partial cancellation [Reddy et al., 1987]. Therefore, the output of adaptive beamformer
can suffer from signal cancellation and cross-talk effect which can confound subsequent
interaction analysis.
There are several possible approaches to deal with the linear cross-talk problem.
One is to define new interaction measures that are less sensitive to this effect. An exam-
ple is to use the imaginary part of complex coherence which is equal to zero when
signals are just linearly mixed [Nolte et al., 2004]. However, this method still suffers
from secondary cross-talk; additional pairs of coherent sources can both leak to the two
cortical sites between which we measure interactions, causing a mislocalization of cor-
tical interactions. An alternative approach is to use statistical testing to distinguish true
interactions from cross-talk. This is often attempted using permutation tests or other
non-parametric methods in which surrogate data are used to establish the distribution
of the interaction measure under the case when no interactions are present. Methods
based on phase randomization [Prichard and Theiler, 1994, Andrzejak et al., 2003] and
permutation of trial indices [Brovelli et al., 2004] have been used for this purpose. How-
ever in practice it is very difficult, or in some cases impossible, to find a permutation
33
scheme in which cross-talk is retained while true interactions are removed, as is required
to establish an appropriate threshold. The third approach is to use inverse methods that
are less sensitive to linear cross-talk so that it has little effect on the interaction mea-
sures. We follow this approach by modifying the standard LCMV beamformer. We also
use permutation tests, but rather than using them to differentiate true interactions from
cross-talk as described, we use them to test for the statistical significance of interactions.
We propose the use of a nulling beamformer to address the problem of cross-talk.
The nulling beamformer is a modified version of the LCMV beamformer [Spencer et al.,
1992, Van Veen et al., 1997], where additional nulling constraints are added to cancel
signals from specific cortical locations. Furthermore, a set of eigenvector constraints
can be used to make the methods robust to misspecification of the precise locations and
extents of the cortical sources of interest. We also use the nulling beamformer in a novel
simulation scheme to demonstrate the performance of our method, we simulate MEG
cortical sources using intracranially measured local field potentials (LFPs), which are
believed to be more realistic representations of true neuronal sources [Sutherling et al.,
1988, Lachaux et al., 2007] than simulation models. Since these signals are available
before and after linear mixing, we can compute interactions in both cases and evaluate
the robustness of our method to linear cross-talk. We measure the performance of the
nulling beamformer in conjunction with different interaction measures and compare it
with that of the standard LCMV beamformer and linear imaging methods.
We make the assumption that the approximate locations of neuronal sources of inter-
est have been identified. This can be achieved either by using a standard beamformer
or inverse imaging method, where peaks of reconstructed activation maps are potential
source locations, or with reference to published neuroscience studies that have identified
cortical regions of interest for a particular network. Use of eigenvector constraints, as
we describe below, provides some robustness to misspecification of these locations.
34
The contents of this chapter are organized as following: We first describe several
standard inverse imaging methods: regularized minimum-norm, linearly constrained
minimum variance (LCMV) beamformer. Then we introduce the nulling beamformer
and its strength in assessing cortical interaction analysis by canceling the cross-talk
effects. We extend this methodology to sources covering extended cortical areas by
using eigenvector constraints. After that, we discuss the non-parametric permutation
approach to establish statistical significance and propose a new simulation scheme to
evaluate the performance of nulling beamformer. We end this chapter with the interac-
tion results from the nulling beamformer compared with other inverse methods.
4.1 Standard Inverse Methods
Due to the characteristics of noninvasive EEG and MEG, various inverse methods have
been widely used to localize the cortical activationHowever, since the goal of this work
to accurately measure the cortical interactions and suppress false interaction caused by
the cross-talk effect. We will focus more on the aspect about how to estimate sources
with less cross-talk effect and accurately compute interaction measures based on them.
4.1.1 Forward modeling
As we discussed in chapter 1, MEG forward model relates the sources to the sensor
measurement. The forward model can generally written as:
M = GS + N (4.1)
where M (n
channels
xn
timepoints
) represents the measured magnetic field at each sensor
as a function of time, which is linearly related with the brain activation S (3n
sources
x
n
timepoints
) via forward matrix (lead field matrix) G (n
channels
x 3n
sources
) and N repre-
sents additive noise for the sensor measurements. The lead field fori
th
spatial location
35
is denoted as G
i
(n
channels
x 3). Alternatively, we may constrain the solution to sources
aligned normally to the cortical surface. Then, the forward operator G has dimensions
(n
channels
xn
sources
), and the lead field for spatial locationi is g
i
(n
channels
x 1).
The forward matrix G depends on the shape and conductivity of the selected head
model. It can be computed using a simplified spherical head model (one sphere model
or overlapping sphere model), or more accurately using either boundary element model
(BEM) or finite element model (FEM) that account for the true shape and conductivity
within the head [Ermer et al., 2001].
4.1.2 Regularized Minimum-norm Solutions
From the forward model, a cortical map of the source time series X can be computed
by applying a linear inverse method. However, due to the characteristics of noninvasive
EEG and MEG, the forward matrix G is ill-conditioned and the inverse problem is highly
ill-posed. Therefore there are no unique solution for this inverse problem. Various
approaches have been proposed to search for the “best” estimate of distributed current
density according to different criteria by introducing implicit or explicit constraints on
the allowed current source distributions. In minimum-norm methods, we look for the
solutin with the smallest norm from all possible solutions. Different norm will achieve
current density map with different characteristic. For example, L-2 norm generally gives
smooth solution and L-1norm produces sparse solution [Komssi et al., 2004]. The most
popular one, Tikhonov regularization [Tikhonov et al., 1977] minimum-norm solution,
solves the problem:
^
X = arg min
X

kM-GXk
2
2
+kXk
2
2

(4.2)
36
where is the tikhonov regularization parameter. We can easily compute analytical
solution of this quadratic problem:
^
X = (G
T
G +I)
1
G
T
M (4.3)
= HM
where H = (G
T
G +I)
1
G
T
is the inverse operator.
Minimum-norm method can also be interpreted in Bayesian formulation as many
other similar methods [Baillet et al., 2001].
4.1.3 Linearly Constrained Minimum Variance (LCMV) Beam-
former
The LCMV beamforming method is a spatial filtering technique first applied in radar
and sonar signal processing [Van Veen and Buckley, 1988], that has been widely used in
the analysis of EEG and MEG data [Spencer et al., 1992,Van Veen et al., 1997,Robinson
and Vrba, 1999, Sekihara et al., 2001, Gross et al., 2001]. It is based on the assumption
that the measured signal m at the EEG/MEG sensors is generated by a small numberN
of focal neural sources s(q
i
) at locationsq
i
(i = 1:::N):
m =
N
X
i=1
g(q
i
)s(q
i
); (4.4)
where the neural signal s(q
i
) is a scalar if the source is modeled as orientation con-
strained with respect to the cortical surface, or a 31 vector withx,y andz components
if the source is modeled as orientation free. The sensitivity or lead field g(q
i
) is a vector
in the orientation constrained case, or a 3-column matrix in the orientation free case.
An LCMV beamformer constructs a spatial filter whose output ^ s(q
i
) at locationq
i
is
represented as
^ s(q
i
) = w
T
(q
i
)m: (4.5)
37
The weights of the spatial filter w(q
i
) are selected to minimize the variance, or power
at the filter output subject to passing signals from a cortical region of interest with unit
gain:
min
w(q
i
)
tr[w
T
(q
i
)C
m
w(q
i
)] subject to w
T
(q
i
)g(q
i
) = I; (4.6)
where C
m
denotes the spatial covariance matrix of the measurement data, and I is a
3x3 identity matrix in the orientation free case or simply the scalar 1 in the orienta-
tion constrained case. This optimization problem can be readily solved using Lagrange
multipliers [Van Veen et al., 1997]:
w(q
i
) = C
1
m
g(q
i
)[g(q
i
)
T
C
1
m
g(q
i
)]
1
: (4.7)
These weights allow the beamformer to adaptively reduce noise and interference, while
passing the desired signal through the filter without attenuation. Therefore, the LCMV
beamformer can suppress cross-talk from other sources when the sources are non-
coherent. However, when signals exhibit coherent behavior, the minimum in (4.6)
is achieved by allowing non-zero gain with respect to the interfering sources so they
can partially cancel the signal from the source of interest, which in turn reduces the
total output power of the beamformer [Van Veen et al., 1997, Van Veen and Buckley,
1988, Brookes et al., 2007]. In that case, the output signals from the beamformer are
reduced and not accurate enough for estimating the interactions. Consequently, the
beamformer in its standard form is of limited use in investigating interactions where we
want to observe correlations between cortical regions.
38
4.2 A Nulling Beamformer to Avoid Signal Cancelation
and cross-talk
4.2.1 Nulling Beamformer
To suppress signal cancellation and cross-talk effects in the LCMV beamformer, we
have to make sure that the filter output at each source locationq
i
will not be affected by
signals from the other locationsq
j
(j6= i). This can be achieved by forcing additional
nulling constraints, i.e. zero gain conditions at interfering source locations [Hui and
Leahy, 2006, Dalal et al., 2006, Hui et al., 2010]:
w
T
(q
i
)g(q
j
) = 0 for everyj2f1N;j6=ig: (4.8)
We combine these nulling constraints with the unit gain condition of the traditional
LCMV beamformer in (4.6). By combining all the gain vectors for the N sources of
interest into a matrix G = [g(q
1
) g(q
N
)], the beamformer design problem can then
be written as follows:
min
w(q
i
)
tr[w
T
(q
i
)C
m
w(q
i
)] subject to w(q
i
)
T
G = f
T
i
; (4.9)
where f
i
= e
i

I. The operator
represents Kronecker product,
e
i
= [0 0 1 0 0]
is an indicator vector whosei-th element is one and the rest are zero, and I is defined as
in (4.6).
Similarly to the traditional LCMV beamformer, this minimization problem can be
readily solved using Lagrange multipliers [Van Veen et al., 1997]:
w(q
i
) = C
1
m
G[G
T
C
1
m
G]
1
f
i
: (4.10)
39
The resulting weight w(q
i
) has a unit gain at the location of interestq
i
and zero gains at
the other N-1 locations, thus nulling possible interference between them.
4.2.2 Eigenvector Constrained Nulling Beamformer
The standard adaptive beamformer assume focal sources. In reality, however, broad
cortical regions represented by multiple dipole sources can be active in response to sen-
sory, motor or cognitive tasks. When the total number of sources representing a cortical
region is small, we can apply the above nulling constraints on all sources individually,
to suppress the cross-talk from the entire cortical patch. But when that number is large,
applying a nulling constraint on every point is impossible, because our beamformer,
which has the same dimension as our MEG/EEG measurements, provides us with only
limited degrees of freedom. To overcome this problem, we use eigenvector constraints,
as described in [Van Veen and Buckley, 1988].
Consider a patch represented withk cortical sources at locations (vertices)q
1
q
k
and their lead field matrix G
p
= [g(q
1
) g(q
k
)]. Forcing point constraints at all loca-
tions in the patch can be enforced with the constraint:
w
T
(q
i
)G
p
= r
d
; (4.11)
where r
d
is the desired response vector for the whole region, for example r
d
= 0 for
nulling constraints with the dimension of the zero matrix being 1 x k for the orientation
constrained case, and 3 x 3k for the orientation free case. Instead of using constraints on
every location within the region, we minimize the degrees of freedom used by allowing
a least squares approximation to the desired response on the whole region:
min
w(q
i
)
jjw
T
(q
i
)G
p
r
d
jj
2
: (4.12)
40
To limit the number of constraints to a small numberL for each region, we use a rankL
approximation of G
p
from itsL largest singular values:
G
p
U
L

L
V
T
L
; (4.13)
where
L
is anLL diagonal matrix containing theL largest singular values of G
p
,
and U
L
and V
L
are matrices containing theL corresponding singular vectors.
By substituting (4.13) in (4.12), we find that the weights minimizing the squared
error should satisfy the following equation:
w
T
(q
i
)U
L
= r
d
V
L

1
L
: (4.14)
This has the same form as the constraint in (4.9), and therefore the solution to the eigen-
vector constrained nulling beamformer derives directly from the solution to the point
constrained one. The only difference between the two beamformers is that instead of
forcing constraints on individual locations, we now force constraints on the principal
eigenvectors of cortical patches.
The advantage of this eigenvector constraint method is that it only requires a small
number of constraints for an extended cortical source. Therefore, when the extents
and true positions of the sources are only approximately known, we can still apply the
nulling constraints by controlling a relatively larger patch around the true source.
In addition to applying extended nulling constraints on interfering patch sources, we
can also impose the constraint of approximately unit gain using eigenvector constraint
on an extended cortical patch of interest where r
d
= 1 of dimension 1k for orientation
constrained case and r
d
= 1
I, a stack of identity matrices with total dimension of
3 3k, for the orientation free case.
41
Combining the eigenvector nulling constraints and eigenvector unit constraints, we
can formulate the eigenvector constrained nulling beamformer similarly to the mini-
mization problem in (4.9) and obtain the weights of the spatial filter using Lagrange
multipliers.
min
w(q
i
)
tr[w
T
(q
i
)C
m
w(q
i
)] subject to
w
T
(q
i
)[U
1
L
U
N
L
] = [r
1
d
V
1
L

1
L
1
r
N
d
V
N
L

N
L
1
] (4.15)
where the superscript indicates the patch number. For estimating the source signals
at the ith patch, we have r
i
d
= 1 for orientation constrained case or r
i
d
= 1
I for
orientation free case, and r
j
d
= 0 for everyj6=i as in (4.11). In (4.15) we assume that
the number of eigenvector constraints is the same for all the patches, however the result
generalizes directly to the case where different dimensions are used for each patch.
4.3 Tests for Significance for Interaction Networks
We assume that the MEG/EEG study is repeated multiple times, and the data are col-
lected as a set of stimulus-locked event-related trials. We estimate the data covariance
matrix C
m
over all trials and compute the weights w(q
i
) for the nulling beamformer
with point sources (4.10), or extended sources (4.15). Time series at each cortical region
of interest for each trial are estimated as in (4.5).
Since we are only interested in the induced response of the experiment, we remove
the evoked response by subtracting from each source its average time series. The
source time series are then used to estimate any of the interaction measures described
above. For example, we can estimate coherence or phase synchrony between any pair
of sources, or fit an MV AR model simultaneously to all sources. Provided our modeling
assumptions are correct and the sources have been localized correctly, the interaction
measures should be unaffected from cross-talk effects between the identified sources.
42
To test for the statistical significance of interactions, we rely on permutation tests
[Brovelli et al., 2004, Nichols and Holmes, 2002, Pantazis et al., 2005]. We assume
that the induced response signals from different trials are independent and identically
distributed. Under the null hypothesis that sources have no interactions between each
other, we can exchange the trial labels for each source separately. These permutation
samples would be statistically equivalent to the original data under the null hypothe-
sis, since they would share no interactions among the sources. Under the alternative
hypothesis that some sources are interacting, randomizing trial indices would destroy
this interaction, and therefore the computed non-zero values of interaction measures
(coherence, phase synchrony, or MV AR causality) in the permutation resamples would
merely appear by chance, which allows us to establish statistical significance for the
original data. We use the permutation samples to create the empirical distribution of a
given interaction measure, and estimate a threshold that leaves (typically 5%) of the
distribution on the right side. We then reject the null hypothesis of no interaction if the
statistic of the original data falls above this threshold.
The MV AR interaction measures, DTF and PDC produce causality estimates
between each pair of sources for each frequency. In this case, we have a multiple com-
parisons problem corresponding to one hypothesis test per frequency and source pair. To
control false positives, we follow the maximum statistic approach [Nichols and Holmes,
2002, Pantazis et al., 2005]: for each permutation sample we compute the maximum of
interactions over all frequencies in a range of interest and over all possible pairs. We use
these values to estimate the maximum distribution and define a level threshold that
controls the family-wise error rate (FWER) over all source pairs and frequencies. This
threshold is applied to the original MV AR measures, identifying significant interactions
among all sources in the network.
43
We note that this procedure can be applied only because we process the data to pre-
clude cross-talk between estimated signals. In the presence of cross-talk this method
would not adequately control false positives, as explained in more detail in the Discus-
sion section. Similar methods have been applied from depth electrode data [Brovelli
et al., 2004] where cross-talk is also not a major issue.
We also use another resampling scheme, the bootstrap, to measure the stability of our
interaction measures. The variability of the estimated interaction between sources can
be attributed to two factors: the trial-to-trial variability in the true interaction between
the sources, and the variability introduced from the linear mixing of sources due to
inaccurate inverse modeling. To measure these factors, we apply a bootstrap resam-
pling scheme. In a simulation setting, the variability of the true interaction is estimated
by bootstrapping the true source signals (resampling with replacement over trials) and
recomputing the interaction measures with the MV AR model. We then measure the vari-
ability of the estimated source interactions by using the nulling beamformer model with
the MV AR model, and test their stability when both models are combined.
4.4 Results
We investigate the ability of the nulling beamformer in detecting source interactions
when used in conjunction with coherence, phase synchrony, or MV AR models. To sim-
ulate realistic sources, we use intracranial EEG recordings experimentally measured
from macaque monkeys [Bressler et al., 1993]. We compare the nulling beamformer
against other inverse methods, and test its ability to measure true interactions when the
exact locations and extents of the sources are not known.
44
4.4.1 Realistic simulation of interactions
To evaluate the performance of our method to detect cortical interactions, we need to
know the ground truth of the interaction network. Since this information is unavailable
in real MEG data, we resort to a simulation with known source configurations and inter-
action profiles. In [Hui and Leahy, 2006] we used time series generated from an MV AR
model as cortical sources signals. Although this simulation is helpful to illustrate the
effectiveness of the nulling beamformer, the time series do not resemble the oscillatory
nature and noise properties of real EEG/MEG data, and therefore may not fully indi-
cate the performance of the nulling beamformer. In this study, we used experimentally
measured intracranial LFPs [Bressler et al., 1993] as time series for the simulated MEG
sources. Intracranial LFP recordings should reasonably resemble the electrophysiologi-
cal neuronal sources that generate the electrical field measured in extracranial EEG and
the magnetic field measured in MEG [Sutherling et al., 1988, Lachaux et al., 2007].
The LFP data were existing and well studied recordings that used surface-to-depth
electrodes from an experiment involving a highly trained macaque monkey performing a
GO-NOGO visuomotor pattern discrimination task [Brovelli et al., 2004]. We selected a
segment of 300 milliseconds after stimulus onset for the GO condition from each of 350
trials recorded on 5 electrodes, and assigned them to sources in a human cortical surface
at locations approximately corresponding to the original placement of the electrodes
on the monkey cortical surface (Figure 4.1). Both focal dipole sources and extended
patch sources were simulated. We modeled the geometry of a CTF 275-channel axial
gradiometer MEG system (VSM MedTech Ltd., CTF Systems, Coquitlam, BC, Canada)
and the lead field matrix was estimated using BrainStorm [Mosher et al., 2005] based on
an overlapping spheres head model [Huang et al., 1999]. Gaussian white noise was sim-
ulated and added to the MEG sensors to model instrumentation noise with SNR=6. The
45
Figure 4.1: Simulating cortical interactions with realistic oscillatory sources. (a) 5
dipole sources are selected on the cortical surface based on the locations of surface-to-
depth electrodes in a macaque monkey experiment. LFPs from the implanted electrodes
are assigned as time series to the sources. (b) MEG data are simulated by projecting
the LFP signals to the channel topography of a 275-channel gradiometer MEG system.
White Gaussian noise is added to the data. (c) Different inverse methods are used to pro-
duce cortical activation maps. (d) The estimated source signals are evaluated for cortical
interactions using MV AR models, coherence, phase synchrony or other measures. (e)
Significant interactions are found using permutation tests while controlling for multiple
comparisons.
SNR was defined as the ratio of the Frobenius norm of the signal-magnetic-field spatio-
temporal matrix to that of the noise matrix for each trial as in [Sekihara et al., 1997].
The source time series at the 5 locations of interest were estimated using the standard
LCMV beamformer, nulling beamformer, and cortically constrained linear minimum
L2 norm imaging. The different interaction measures (coherence, phase synchrony and
MV AR model) were then applied to each of the sets of estimated time series.
46
4.4.2 MV AR analysis with focal sources
Figure 4.2 shows the estimated interaction network of the 5 dipole cortical sources
with MV AR analysis. We used the Partial Directed Coherence (PDC) as our interaction
measure, but we observed that the Directed Transfer Function (DTF) produced quali-
tatively similar results (see [Baccala and Sameshima, 2001], [Baccala and Sameshima,
2001] for the differences between PDC and DTF). The MV AR model was fitted to the
LFP time series representing the true sources, as well as the reconstructed time series
with Tikhonov regularized minimum-norm [Tikhonov et al., 1977], LCMV beamformer
and the nulling beamformer. For minimum-norm imaging, we chose the regularization
parameter as 5% of the maximum singular value of the lead field matrix. We subtracted
the ensemble mean from the true and estimated source signals to obtain the induced
responses. The PDC results appear as a 5x5 matrix of subplots where the subplot at the
i-th row and j-th column represents the interaction from cortical region j to region i.
Because PDC between each region pair was computed as a frequency spectrum, we use
the height of the shaded area to indicate the value of the computed PDC spectrum.
The permutation method described in Section 4.3 was used to threshold the inter-
action network at error level = 5% with correction for multiple comparisons. For
each permuted sample, we computed the maximum value of PDC over the frequency
band of [0; 50] Hz for each pair of sources, then get the maximum among all the pairs.
We generated the empirical distribution of the maximum values from 2000 permutation
samples. The PDC values above the threshold were regarded as significant interactions
and shown in figure 4.2 (right).
The MV AR results from the nulling beamformer are in close agreement with those
from the true source time series. This indicates that the absence of cross-talk prevented
the inclusion of spurious interactions in the estimated network. In contrast, the MV AR
47
PDC Significant interactions
(a) True source (b) Minimum-norm (c) LCMV (d) Nulling
Figure 4.2: (Left) PDC interaction measures for five dipole sources using different MEG
inverse methods. Each subplot represents the PDC causal interaction between a pair of
sources, where the subplot ati-th row andj-th column represents the causal interaction
from sourcej to sourcei. The horizontal axis represents frequency from 0 to 50Hz, and
the vertical axis represents the value of PDC in the range of [0; 1]. (Right) Significant
PDC interactions found for five dipole sources with a permutation test at an = 5% sig-
nificance level corrected for multiple comparisons. The arrows represent the directions
of causal interactions, with their radii indicating the maximal strengths of interactions
over the [0; 50] Hz frequency band.
48
results from the minimum-norm method and the LCMV beamformer are noticeably
affected by linear mixing, producing both type I and II errors (false positives and true
negatives). Even for source pairs where a true interaction is identified, the PDC measure
is distorted as a result of linear mixing and has a different frequency profile.
4.4.3 MV AR analysis with extended cortical sources
As described before, the nulling beamformer successfully cancels linear mixing in the
case of point sources where only one degree of freedom is needed per source. A more
realistic simulation involves extended sources, where we assigned the same time series
to neighboring surface elements representing a cortical patch. We assumed the activ-
ity is uniform within each patch. We used the same source locations as before, but
simulated extended sources with patch sizes of approximately 8 cm
2
each as shown in
figure 4.3a. We then applied the nulling beamformer with eigenvector constraints as in
(4.15), where we assumed the patch locations and sizes are exactly known (we will relax
this assumption later).
Since we relied on the principle eigenvectors to approximate the response for
extended sources in (4.14), we first determined the number of eigenvector constraints
L to be used for each source patch. Our criterion was to keep all eigenvectors with
eigenvalues at least 10% of the maximal eigenvalue for each patch, which resulted in
L = 4 for all the patches.
We applied the MV AR analysis to time series produced by the same inverse methods
as before, and the results are shown in figure 4.3. For the minimum-norm method we
computed the time series at the center of each cortical patch of interest. For LCMV , we
used the same eigenvector constraints for each source patch in turn that were used in
the nulling beamformer, but without the nulling constraints. The nulling beamformer
with eigenvalue constraints reproduced the interactions present in the true source time
49
PDC Significant interactions
(a) True source (b) Minimum-norm (c) LCMV (d) Nulling
Figure 4.3: (Left) PDC interaction measures for cortical sources with extended sizes of
approximately 8 cm
2
. (Right) Significant PDC interactions found for the five extended
cortical sources with a permutation test at an = 5% level corrected for multiple com-
parisons.
series. This is not true for the minimum-norm and LCMV beamformer methods, where
the linear mixing severely affected their performance.
50
minimum-norm LCMV Beamformer Nulling beamformer
with noise
without noise
Figure 4.4: (Top) PDC interactions computed using different inverse operators from
simulated MEG data (with noise). (Bottom) PDC interactions computed using the same
inverse operators but from simulated data without the noise. The difference between
them shows the levels of false interactions caused by residual noise.
To test whether the inaccurate estimation of interactions in figure 4.3 is caused by
cross-talk or the difference in residual noise level for each method, we applied the same
inverse operators to the simulated MEG data without adding any noise. Figure 4.4 shows
the estimated PDC with and without noise; the difference between them for each inverse
operator is minimal, therefore most of the false interactions from minimum-norm and
LCMV beamformer were not caused by the residual noise, but by cross-talk.
To test the reliability of the combination of nulling beamformer and MV AR analy-
sis to detect interactions, we also performed bootstrap analysis as described in Section
4.3. Figure 4.5 shows the same MV AR interaction results as in figure 4.3 (a) and (d) ,
but also with confidence interval lines plotted at2 standard errors. The nulling beam-
former produces interactions with variability close to the variability of the true source
interactions.
51
(a) True source interactions (b) Nulling beamformer
Figure 4.5: PDC interactions with confidence intervals. Interaction results are the same
as in figure 4.3a and 4.3d, but overlaid with [2; +2] confidence lines.
4.4.4 Nulling beamformer with other interactions measures
Even though we mainly used the nulling beamformer in conjunction with MV AR anal-
ysis, it has similar benefits when used with other interaction measures. Here we show
the performance of the nulling beamformer with coherence and phase synchrony. We
simulated the same extended sources as Section 4.4.3 and estimated the source time
series at the center of each region of interest using each of the inverse methods. We
then computed the coherence and phase synchrony between each pair of the estimated
signals.
In figure 4.6, we contrast the accuracy in measuring coherence using the nulling
beamformer versus other inverse methods. The coherence between all pairs of extended
cortical sources has been estimated using all inverse methods. The nulling beamformer
produces results very close to the true coherence across all frequencies. This is not
the case with the LCMV beamformer and the minimum-norm method, which typically
overestimate the coherence for all frequencies.
In addition to coherence, we also estimated the phase synchrony between sources
3 and 4 using the wavelet method [Lachaux et al., 1999, Le Van Quyen et al., 2001].
The same extended sources were used as before, and the results from different inverse
methods are shown in figure 4.7. The true source interaction has a phase synchrony peak
52
Figure 4.6: Coherence between all the possible pairs of extended cortical sources esti-
mated using different inverse methods.
at 10-15Hz and 0.1-0.3sec, which was accurately detected with the nulling beamformer.
Both Minimum-norm and LCMV beamformer method overestimated the magnitude of
the synchrony and identified additional false interactions.
4.4.5 Robustness to inaccurate patch size
Here we relax the requirement that the nulling beamformer knows the exact locations
and extents of the cortical sources. This is important, because we can only approxi-
mately specify the locations of cortical sources in practice. We used the same extended
cortical source simulation data as in Section 4.4.3, but we now assume that the exact
locations and sizes of the sources are unknown when calculating the nulling beamformer
53
(a) True sources synchrony (b) Minimum-norm method
(c) LCMV beamformer (d) Nulling beamformer
Figure 4.7: Phase synchrony between sources 3 and 4, computed using different inverse
methods. The color scale shows the phase locking value between the two sources.
with (4.15). In particular, to design the nulling beamformer, we use larger patches to
estimate the eigenvalue constraints. Those patches, shown in figure 4.8 (right), have
sizes from 29 to 45 cm
2
as compared to the 8 cm
2
approximate sizes of the true cortical
sources. We usedL = 6 for the number of the eigenvector constraints for all the patches,
based on the same criterion used in Section 4.4.3.
Figure 4.9 compares the PDC results for the true cortical sources, the nulling beam-
former with known patch sizes and locations, and the nulling beamformer with larger
patches that approximate the true sources. The PDC results from the larger regions are
similar to those from the true interactions with the exception of small amount of false
interactions from source 2 to source 5 and source 2 to source 4, which are detected as
significant after permutation. Despite this false positive, the results demonstrate that the
54
Figure 4.8: Simulation of cortical sources in the case where we have limited knowledge
of their true locations and extents. Left: true sources with approximately 8 cm
2
spatial
extents. Right: enlarged patches of around 29 to 45 cm
2
spatial extents.
nulling beamformer has a certain degree of robustness to approximation of the cortical
locations with substantially larger cortical patches.
True Interactions
Patch size known Larger patch size
(a) (b) (c)
Figure 4.9: Top: PDC interaction measures computed using (a) the true LPF signal, and
reconstructed time series using the nulling beamformer with (b) known patch sizes and
(c) enlarged patch sizes. Bottom: significant PDC interactions at an = 5% significance
level after correcting for multiple comparisons using a permutation test.
4.4.6 Effect of noise
In this section we explore how noise affects the performance of the nulling beamformer,
LCMV beamformer and minimum-norm in detecting cortical interactions. We consider
not only white noise in the sensor domain, which models instrumentation noise, but
55
also background brain activity. The latter was created using experimental data from the
same MEG sensor configuration as our synthesized data, and it was recorded during a
baseline period when the subject was resting with eyes open. We considered several
levels of SNR for both sensor noise and background brain activity.
We define the PDC estimation error as the Frobenius norm of the difference between
the estimated PDC and the ground truth. In figure 4.10 we show how the above error
depends on the SNR level, and in figure 4.11 we show the estimated PDC coefficients for
SNR=3 for all methods. The nulling beamformer has lower error than the LCMV beam-
former and minimum-norm for all cases with SNR>0.75. Also, the nulling beamformer
with point sources is more robust to noise than the one with patch sources.
Ongoing brain activity affects the performance of all the methods more severely than
channel white noise. This is expected, as ongoing brain activity lies in a similar subspace
as our sources of interest, so cross-talk increases in all inverse methods. For example,
in the nulling beamformer a large part of the brain noise lies in areas not explicitly
canceled, therefore these signals can potentially leak to the true source locations.
Figure 4.10 reveals interesting behavior: PDC estimation error decreases for the
nulling beamformer as SNR increases, but the opposite is true for the LCMV beam-
former and minimum-norm. This potentially unexpected behavior of the LCMV beam-
former and minimum-norm can be explained by considering how noise affects the esti-
mated interactions. Without noise, both methods tend to overestimate interactions (fig-
ure 4.4 (top), versus figure 4.3a where the true interactions are shown). When noise
is increased, both true and false interactions are suppressed (figure 4.11), and there-
fore approach the true interactions. At very low SNR levels, most methods produce an
error close to 6, which is actually the Frobenius norm of the true interactions. Since we
defined the error as the Frobenius norm of the difference between the estimated PDC
56
and the ground truth PDC, this indicates that at very low SNR most methods estimate
close to zero PDC in all frequencies.
Figure 4.10: Comparison of the errors in estimated PDC interactions for different inverse
methods and different noise level.
4.5 Discussion
The nulling beamformer accurately detected interactions in all cases tested, namely focal
cortical sources, and patches with known or unknown exact locations and sizes. This
is because it effectively suppressed interference from the nulled sources and therefore
linear mixing was greatly reduced from these locations. Any interaction measure can
be used in conjunction with the nulling beamformer. In contrast, when other inverse
methods were used, interaction analysis was prone to significant errors because of linear
mixing and/or signal cancellation.
57
(a) True Interactions (b) Nulling/white noise (Point) (c) Nulling/brain noise (Point)
(d) Minimum-norm/white noise (e) LCMV/white noise (f) Nulling/white noise
(g) Minimum-norm/brain noise (h) LCMV/brain noise (i) Nulling/brain noise
Figure 4.11: PDC interaction computed from different noise types with SNR=3: (a) the
true interactions, (b) and (c) the interactions computed using nulling beamformer from
dipole sources simulation with white noise and brain noise, (d), (e) and (f) from patch
simulation with white noise using different inverse methods, (g), (h) and (i) from patch
simulation with brain noise.
Linear mixing can affect estimation of interactions in two ways. It may cause type I
errors, in which case false interactions are detected when there is no true interaction, or
type II errors, in which case we fail to detect the existing true interactions. For exam-
ple, the minimum-norm and LCMV beamformer resulted in several type I errors in the
MV AR analysis with focal and extended sources, including false interactions between
sources 1 and 3, 2 and 4, 2 and 5, and many others in figure 4.2 and figure 4.3. Both
the minimum-norm and LCMV beamformer produced type I errors in the estimation of
58
coherence, and minimum-norm even produced type I errors in phase synchrony. Exam-
ples of type II errors include the failure of minimum-norm method to detect the PDC
interaction from source 2 to 1 in the patch configuration (Figure 4.3).
We should note that the performance of the standard LCMV beamformer is highly
dependent on the choice of the data used in estimating the data covariance in (4.7). Here
we have selected a time window corresponding to the data for which we are interested
in analyzing coherence as the basis for estimating the beamformer weights. In practice,
using a different time window will produce a different response in the beamformer and
hence different degrees of cross-talk. The important issue of the dependence of the
beamformer on the time window used has recently been explored in [Brookes et al.,
2008].
In the presence of coherent source activity, the LCMV beamformer can cause signal
cancellation, as mentioned in Section 2.1. Therefore, inaccurate estimation of source
interactions with the LCMV beamformer can be caused either by linear mixing of inco-
herent third sources, or self cancellation of coherent ones.
An important limitation of the nulling beamformer is the requirement that we know
the locations and extents of the sources to be canceled. However, we demonstrated that
there is some tolerance in specifying the exact spatial configuration of the sources when
eigenvalue constraints are used. In particular, when the exact locations of the sources
are not known, we can use evidence from other inverse imaging methods or previously
published studies to approximate the locations of activated cortical areas and force a null
in the principal subspace of the corresponding patches. For cortical patches spanning a
few cm
2
, a few eigenvalues are enough to null the majority of the subspace.
Our permutation approach relies on the fact that channel noise and ongoing brain
activity on cortical areas outside the nulled sites have minimal impact on the estimated
source time series with the nulling beamformer. In this case, the linear mixing primarily
59
occurs between the identified cortical sites, and the nulling beamformer can effectively
cancel it by zeroing the interfering sources. Therefore, the original data in the permu-
tation test contain only the true interactions, which are completely destroyed by ran-
domly exchanging the trial indices of each source in the permuted data. However, since
we cannot completely ignore linear mixing from noise and other cortical sources, our
permutation test is potentially not exact, but rather liberal; the original data have inter-
actions due to linear mixing which are not canceled with the nulling beamformer, but
are completely eliminated in the permutation samples. Our permutation test results in
the MV AR analysis did not suffer from false positives when the locations of the sources
were known, indicating that the method is relatively robust, but we cannot guarantee an
exact level control of the error rate in real data situations.
The bootstrapping analysis demonstrated that the MV AR models, when combined
with the nulling beamformer, are reliable estimators of true cortical interactions. Since
linear mixing is suppressed with the nulling beamformer, the only major variability in
the interaction measurements is the one introduced by the MV AR model estimation.
The nulling beamformer is more sensitive to noise than the LCMV beamformer
for two reasons. First, because the nulling beamformer uses additional constraints to
suppress the cross-talk, it has less degrees available to suppress noise than the LCMV
beamformer. In our simulations, we observed that the white noise gain increased from
the LCMV beamformer to the nulling beamformer. We also observed that the point
constrained nulling beamformer had greater noise gain than the eigenvector constrained
one. Second, when several spatial constraints are enforced, there is increased chance
that the subspace of the nulling constraints has small subspace angles relative to the sub-
space for the unit constraints. This can cause the norm of the weight vectors to become
larger, and therefore the noise gain increases. In our simulations, sources with small
angles between the subspaces of the nulling and unit constraints had greater noise gain.
60
Therefore, when designing a nulling beamformer, the overlap between the aforemen-
tioned subspaces should be considered. The estimation of the angle between subspaces
is discussed in [Bjorck and Golub, 1973].
In comparing our network results to those in [Brovelli et al., 2004], who used the
same LFP data, we see some differences in the locations exhibiting significant interac-
tions. This is due in part to differences in the time-window used for analysis. However,
in exploring these differences further we found that PDC and Granger causality can
produce interactions with different amplitudes from the same data. After testing for sig-
nificance these differing amplitudes can result in differences in the final network model.
While the methodological implications of these differences are clearly important, a com-
parison of different interaction measures is beyond the scope of this paper.
4.6 Summary
We have described a method to detect cortical interaction networks while at the same
time controlling for linear mixing among several cortical sources. The method relies on
the nulling beamformer and its extension with eigenvector constraints. We demonstrated
that the method is superior to other inverse imaging methods, such as minimum norm
imaging and conventional LCMV beamforming, in estimating cortical interactions.
61
Chapter 5
Parametric phase synchrony and group
analysis of working memory effect
using intracranial EEG
In the previous chapter, we have shown that non-invasive brain imaging techniques, such
as EEG and MEG, have limitations for assessing brain interactions. In particular, their
limited spatial resolution causes cross-talks effects and limits their ability to detect true
interactions. Invasive recordings from intracranial EEG provide us a unique opportunity
to directly measure cortical activations and analyze functional interactions with minimal
cross-talk effects. However, intracranial EEG has its own limitations: the implanted
electrodes only sample a limited set of cortical locations. Moreover, because most
intracranial EEG is recorded during patients under monitoring before brain surgery, most
electrodes implantation schemes focus on monitoring pathological conditions instead of
reflecting brain regions related to the experiments. In addition, because pathological
conditions are different for each patient, the coverage of brain regions among multiple
subjects varies significantly. Therefore, group analysis of interactions using intracranial
recording is very challenging.
In this chapter, we describe a novel group analysis scheme for investigating func-
tional interactions using intracranial recordings, which involves brain image registra-
tion, combining p-values, and parametric tests of phase synchrony. We apply the
scheme to experimental intracranial EEG recordings from subjects performing a Stern-
berg working-memory task. The chapter is organized as follows: in section 5.1, we
introduce the working-memory experiments during which the intracranial EEG data
62
was recorded. In section 5.2, we describe a parametric statistical test for equality of
phase synchrony. The test is used to establish statistical significance of interactions for
each subject. In section 5.3, we discuss the image registration method used to map the
anatomical information from multiple subjects. Then in section 5.4, we describe the
methods used to find the significant interactions at group level by combining thep-value
computed for interactions of each electrode pair. In section 5.6 we present results of our
analysis and discuss their neuroscience significance.
5.1 Working Memory Experiment
Memory is an important and fundamental cognitive function of the brain, which still
is not fully understood. Working memory, as a cognitive psychological concept, refers
to the memory system for temporarily storing and manipulating the information neces-
sary for complex cognitive tasks [Baddeley, 1992]. Working memory definition evolves
from the short-term memory concept, and has drawn a great amount of attentions in
neuroimaging studies.
Many studies have used functional Magnetic Resonance Imaging (fMRI) [Badde-
ley, 1992] [O’Sullivan et al., 1996] or positron emission tomography(PET) [O’Sullivan
et al., 1996] to study the organization of working memory. However, field potential
recordings from noninvasive image modalities such as EEG and MEG, or from invasive
intracranial EEG recording, offer us a unique opportunity to study oscillatory activities
associated with working memory. By using EEG/MEG or intracranial recordings, we
are able to obtain the information of active cortical areas involved in different working
memory tasks, and analyze the frequency dependent nature of the oscillatory activities
as well.
The majority of the previous work was focused on analyzing the oscillatory power
associated with the working memory task on different frequency bands. [Klimesch,
63
1999] gave an in-depth review of EEG alpha and theta oscillatory power in working
memory tasks. [Jensen and Tesche, 2002] used MEG to study the frontal theta activ-
ity modulation associated with different memory tasks. [Mitchell et al., 2008] studied
the involvement of frontal-midline theta in processing of short-term memory and atten-
tion. [Npflin et al., 2008] showed the stability of working memory related EEG spectra
for individual subjects over time. Intracranial EEG has also been used for the same pur-
pose. [Raghavachari et al., 2001] observed gating effects of the theta oscillations during
the work memory task using intracranial EEG. [Meltzer et al., 2008] used intracranial
recording to find the changes of oscillation power modulated by the load of memory
task .
Compared with oscillatory power, functional interactions during the working mem-
ory task are less studied. [Wu et al., 2007] used non-invasive EEG to study large-scale
synchronization at the theta frequency, and [Sarnthein et al., 1998] [Kumar et al., 2009]
studied coherence using EEG. [Jensen and Tesche, 2002]. [Axmacher et al., 2008] [Sia-
pas et al., 2005] [Siegel et al., 2009] analyzed the functional connectivity using intracra-
nial EEG or local field potentials.
5.1.1 Sternberg working-memory task
The intracranial EEG data used in this work were recorded from 12 epilepsy patients
performing the Sternberg working memory task [Meltzer et al., 2008], while they were
under invasive monitoring of seizure activities at the Yale New-Haven Hospital. The
original data were recorded from 17 patients, but two were excluded because of large-
scale epileptiform activity, one because of low-quality anatomical scans, and two more
because of recording artifacts.
The Sternberg work-memory task is a classic test in cognitive neuroscience, which
allows the researcher to parametrically control the memory load of the working memory
64
Figure 5.1: The time diagram of the Sternberg working-memory task
experiment by keeping other components of the task constant [Sternberg, 1966]. In
the Sternberg task, the subjects are asked to maintain a list of items, such as digits and
letters in their working memory during a memory retention period. The number of items
is parametrically controlled and referred to as memory load.
The Sternberg working-memory task used in the experiment is a slightly modified
version of the one used in a related EEG-fMRI study [Meltzer et al., 2007] to make it
easy for the patients under monitoring condition to perform. Three memory load levels
were used, representing 1, 2 and 4 digits that the patients needed to remember. At the
beginning of each trial, a “start“ cue was presented to indicate the beginning of each
trial and inform the patients to disregard previous digits. A horizontal list of 1, 2, or 4
digits was presented for 4 seconds, and replaced by a matching list of horizontal lines
during the 6 seconds of the memory retention period (or delay period). The patients
were then asked to press a button within 3 seconds to respond to a probe digit appeared
between two question marks. The experiment was organized as 12 sessions, with each
session lasting 4 minutes and containing 15 trials. Within each session, the memory
load level was kept constant. Each memory load level was conducted in 4 sessions. The
time sequence of each trial is shown in figure 5.1.
The intracranial EEG data were recorded using a 128-channel Ceegraph clinical
monitoring system (BioLogic Systems Corp, Mundelein, IL). The implanted electrodes
were either square grids or rectangular strips of stainless steel subdural electrodes on
frontal, temporal, parietal, occipital and midline surfaces, except for some multi-contact
65
depth electrodes in the medial temporal lobe. Each of the electrodes had an exposed
surface of 2.3 mm in diameter, and spaced at 1 cm intervals. Although the placement
of electrodes within each patient was determined solely by clinical criteria, the overall
coverage of the brain among all the subjects are broad enough for us to conduct group
analysis using the data.
Intracranial electrodes were monopolar and the reference electrodes were placed
between the inner and outer laminae of the skull, as far as possible from other electrodes.
Signals were sampled at 256 Hz and electronic marker signals were recorded from the
stimulus presentation computer if available.
5.2 Parametric Analysis of Phase Synchrony
From these unique intracranial recordings, we are interested in investigating the func-
tional interactions involved in the memory retention period at the group level. Among
many interaction metrics discussed in chapter 3, we choose to use phase synchrony
after preliminary tests of several measures on the data. We notice that phase synchrony
detects strong synchrony at a narrow frequency in theta band during the memory reten-
tion period as opposed to coherence, which produces spectrally smooth interactions. It
is also consistent with the previous results reported in studies of non-invasive EEG and
MEG.
As discussed in chapter 3, phase synchrony is estimated as the degree of peakiness
in the distribution of the relative phases on the unit circle over a time window, and it
is quantified as the magnitude of the average of all complex phasese
j
. In this way,
phase synchrony is represented as a single real number, the phase locking value. It is
convenient to use but it does not preserve the information of the statistical distribution
of the relative phase. For example, it is difficult if we want to compare phase synchrony
between conditions using this definition.
66
In this chapter, we treat the phase differences as circular data, and use circular statis-
tics/directional statistics to analyze the phase synchrony.
5.2.1 Circular data
Circular data refers to the data measured in the form of angles (phases) or two-
dimensional orientations, or other cyclic data which can be transformed into the form of
angles, such as times on a clock, months in a year, etc. Circular data are widely seen in
various scientific and engineering fields such as biology, geography, social sciences and
engineering, and their analysis needs special treatment.
5.2.2 Representation of circular data
Because each observation for circular data can be represented as a data point on the
unit circle, and each point on the unit circle can be uniquely determined by the angle
, we can use the unit vector x in Cartesian coordinate system to represent a circular
observation:
x = (cos; sin)
T
(5.1)
Alternatively, we can represent it as a unit complex number
z =e
j
= cos +j sin (5.2)
This is the same representation as the phase difference used in chapter 3 for phase syn-
chrony. We will later show that the standard definition of the phase locking value is one
of the summary statistics for circular data.
Another way to represent circular data is to use diagrams. We can use a circle data
plot to represent circular data, as shown in figure 5.2. Each observation is plotted as one
point on the unit circle.
67
Figure 5.2: Circular data plot for an example of simulated von Mises distribution data
5.2.3 Summary statistics for circular data
Circular data have very unique properties compared to the common linear data, which
can not be fully captured by the classical statistical analysis. An entire field of circular
statistics has been developed for the special treatment of circular data. In this subsection,
we will show some typical summary statistics for circular data.
Circular data can be graphically summarized using a linear histogram or a circular
histogram. In figure 5.3, we show the histogram plots for some simulated circular data.
If we assume that circular observations are given by the angles,
1
; ;
i
; ;
n
,
or by the vectors, x
1
; ; x
i
; ; x
n
, the center of mass x in the Cartesian coordinates
is (

C;

S) [Mardia and Jupp, 2000]
R, x =
1
n
n
X
i=1
x
i
= (

C;

S)
T
(5.3)
68
Figure 5.3: Linear histogram and circular histogram for circular data. The resultant
vector is denoted by red color in the right figure.
where the first trigonometric moments of the circular data is defined as

C =
1
n
n
X
i=1
cos
i
(5.4)

S =
1
n
n
X
i=1
sin
i
(5.5)
Here R is referred as the resultant vector of the circular data. For example, the
resultant vector for the simulated data in 5.3 is shown in red color in the right figure.
The mean direction is given by

= arctan

S

C
(5.6)
This mean direction is the counterpart of the zero order moment (mean) for linear (non-
circular) data.
The magnitude of the resultant vector is called the mean resultant length, [Mardia
and Jupp, 2000]

R =
p

C
2
+

S
2
(5.7)
69
Because all of the circular data are on the unit circle, the resultant will alway fall
into the unit circle. Therefore, the length of the resultant is constrained in the range:
0

R 1 (5.8)

R is a very useful summary statistic because it tells us how concentrated the data is
towards the center of mass.

R = 0 means the circular data is evenly distributed on the
unit circle without any concentration; While

R = 1 means the circular data fall into
the same direction all the time and they have perfect concentration. We can see that the
concentration measure is opposite to the variance in non-circular data. For circular data,
a statistic similar to variance is called the circular variance or dispersion. It is defined as
V = 1

R (5.9)
5.2.4 Circular data and phase synchrony
With circular data and the summary statistics, we can now treat the phase difference as
circular data and investigate the connection between the standard phase locking value
and the circular statistics.
Phase synchrony between two signals is measured as the degree of concentration
of their phase differences (relative phases). After being wrapped into [0; 2] range, the
relative phases become a set of circular data. The phase locking value is just the resultant
length of these circular data by definition, which measures the exact concentration of the
relative phases. A large resultant length indicates the concentrated phase differences, or
strong phase synchrony, and vice versa. Therefore, we can compute the phase locking
value using different ways of evaluating the concentration of the relative phases.
One of the advantages of treating phase synchrony as circular statistics is that it pro-
vides us an effective way to test for significant phase locking values. Without using the
circular statistics concept, one will need to test for the significance of the phase locking
70
values by using surrogating methods or resampling methods. In resampling methods,
one can resample the original signals using bootstrapping or permutation test to get an
empirical distribution of the phase locking values on the resampled data, and compute
the threshold from it. Both methods are computationally expensive. By treating the
relative phases as circular data, we have more flexibility. To find the significant degree
of concentration, we can resample the wrapped phase differences, or use parametric sta-
tistical test to find out whether the circular data is significantly different from uniform
distribution.
Another advantage of treating phase synchrony using circular statistics is that we
can perform other types of test on the phase difference. For example, we can test for
equality of the phase synchrony between different experimental conditions by testing
the equality of concentration between two sets of phase differences. It will be further
discussed in the following sections.
5.2.5 Probability model for circular data
In order to perform parametric tests for phase synchrony, the basic probability model
used for circular data is shown in this section. Because of the 2 ambiguity of the circu-
lar data, not all probability models used for linear data are valid for circular data. Several
distribution models have been used for circular data such as the von Mises distribution,
wrapped normal distribution, wrapped Cauchy distribution and so on [Jammalamadaka
and Sengupta, 2001,Fisher, 1996,Mardia and Jupp, 2000] . Among them, the von Mises
distribution is mostly used and has most relevant statistical properties.
71
Figure 5.4: The PDF of the von Mises distribution with different, where
0
= 0.
Von Mises distribution
The von Mises distribution is also known as the circular normal distribution. It is a cir-
cular analogue of the normal distribution in linear data. The probability density function
(PDF) of the von Mises distribution for is
f() =
1
2I
0
()
expf cos(
0
)g(0< 2) (5.10)
where
0
is the mean angle or mean direction, and is the concentration parameter.
Here
0
and are analogous to the mean and the variance in the normal. I
0
()is the
modified Bessel function of zeroth order.
I
0
() =
1
2
Z
2
0
e
cos
d (5.11)
The shape of the von Mises distribution is unimodal and symmetric around
0
shown
in figure 5.2.5. Different values of determine the shape of PDF. The larger the ,
the sharper the shape of the PDF is. The von Mises distribution has an interesting
72
limiting behavior. As increases, the von Mises distribution converges to a normal
distribution N(0; 1=). As reaches the smallest value, the distribution becomes the
uniform distribution over an interval of 2.
Phase Locking and the von Mises distribution
The von mises distribution is the maximum entropy probability distribution for a circular
variable given the mean angle and the concentration parameter, in the same way that the
normal distribution is the maximum entropy distribution for all real-valued distributions
with given mean and standard deviation. In practice, if only the mean and the standard
deviation are known for a distribution defined on (1; +1), it is reasonable to assume
the distribution is normal. For distributions defined for circular data, if only the mean
angle and the concentration are known, it is good practice to assume the distribution is
von Mises.
For this reason, we choose the von Mises distribution as the probability model for
testing of phase synchrony. We assume that the phase differences between 2 signals
at certain frequency and over some time window follow the von Mises distribution.
To validate this underlying assumption, we fit the phase difference at the frequency of
interest into the von Mises distribution and test the goodness of fit.
5.2.6 Test for equality of concentration
The phase synchrony is reflected by the concentration of the relative phases. If we
assume the underlying distribution of the relative phases is the von Mises distribution,
any test for phase locking value is equivalent to the test on the concentration parameters
in Eq. (5.10). We can test for phase synchrony against the hypothesis of no phase
synchrony. By using a parametric probability model, we only need to test for uniformity
of relative phases. A lot of such tests exist in the circular statistics literature [Mardia
and Jupp, 2000]. However, it is more interesting in neuroscience to test for the changes
73
of phase synchrony rather than to test against the lack of synchrony. In a neuroscience
study, such as the work memory experiment we described at the beginning of this chap-
ter, we are interested in testing the changes of the phase synchrony between different
experimental conditions. In this section, we will discuss the test for equality of phase
synchrony by using the test for equality of concentration of relative phase.
If we assume the underlying distributions for the two conditions are the von Mises
distributionsM(
1
;
1
) andM(
2
;
2
), the goal is to test
H
0
:
1
=
2
against H
1
:
1
6=
2
(5.12)
where both the mean angles
1
and
2
, and the concentration parameters
1
and
2
are
unknown.
Because is the concentration parameter for von Mises distribution, and the resul-
tant length

R is a circular statistic which also measure the concentration, we will use the
approach in [Mardia and Jupp, 2000] to derive the test of equality on from the limit-
ing distribution of

R. Similar results have also been reported in [Allefeld and Kurths,
2004b].
Limiting distribution of resultant length
If the underline distribution of the circular data is von Mises distributionM(0;), where
> 0, the distribution of

R is asymptotically normal [Mardia and Jupp, 2000], with
mean and variance given by,
E

R

' A() +
1
2n
(5.13)
nV AR(

R) =

1A()
2

A()

1
4n
2
+O(n
2
) (5.14)
where
A() =
I
1
()
I
0
()
(5.15)
74
I
0
() andI
1
() are modified Bessel function
I
p
() =
1
2
Z
2
0
cospe
cos
d (5.16)
For large sample,

R is asymptotically normal with
E

R

'A()
nV AR(

R) =

1A()
2

A()

=g()
Note that both the mean and the variance of

R will depend on. We can use variance-
stabilizing transformations to find a functionh to transform

R such thath(

R) has con-
stant variance [Mardia and Jupp, 2000].
We can writeh(

R) using a first-order Taylor series expansion:
h(

R) h

E

R

+

R E

R

h
0

E

R

(5.17)
h (A()) +

RA()

h
0
(A()) (5.18)
and
h(

R)h (A())

RA()

h
0
(A()) (5.19)

h(

R)h (A())

2

RA()

2
[h
0
(A())]
2
(5.20)
By taking the expectation of both sides, we have
V AR

h(

R)

V AR

R

[h
0
(A())]
2
(5.21)

1
n
g() [h
0
(A())]
2
(5.22)
(5.23)
75
To get the constant variance forh(

R), we need
nV AR

h(

R)

1 (5.24)
and
g() [h
0
(A())]
2
1 (5.25)
h(x) =
Z
x
0
1
p
g(A
1
(u))=n
du (5.26)
Now, we can compute the transformationh(x) directly from the exact expression for
A() using Eq. (5.15), or using some approximation ofA() for different ranges of
as in [Mardia and Jupp, 2000].
Test for equality of concentration
From section 5.2.6, we can get the asymptotic distribution of h(

R) using a variance-
stabilizing transformation as a normal distribution with a mean depending on and a
constant variance
E

h

R

' h (A()) (5.27)
nV AR

h

R

' 1 (5.28)
Then the test for equality of the concentration parameters in Eq.(5.12) can be for-
mulated as a test for equality ofh (A(
1
)) andh (A(
2
)) after the variance-stabilizing
transformation. This is equivalent to a two-samplez-test forh(

R
1
) andh(

R
2
). Under
the null hypothesisH
0
, the statistic
z =
h

R
1

h

R
2

q
V AR

h

R
1

+ V AR

h

R
2

(5.29)
=
h

R
1

h

R
2

q
1
n
1
+
1
n
2
(5.30)
76
is distributed approximately asN(0; 1). We can then use the normal distribution to find
out thep-value of the estimated statisticz.
[Mardia and Jupp, 2000] used an approximation forA() at different ranges of
to get the variance-stabilizing transformationh(x), and got different analytical forms of
the test as follows:
For

R< 0:45,
z =
2
p
3
g
1
(2

R
1
)g
1
(2

R
2
)
q
1
n
1
4
+
1
n
2
4
N(0; 1) (5.31)
whereg
1
(x) = sin
1
(ax) anda =
q
3
8
.
For 0:45<

R< 0:70,
z =
g
2
(

R
1
)g
2
(

R
2
)
c
3
q
1
n
1
3
+
1
n23
N(0; 1) (5.32)
whereg
2
(x) = sinh
1

Rc
1
c
2
(ax) andc
3
= 0:893.
For

R> 0:70, [Mardia and Jupp, 2000] actually used the high concentration approx-
imation instead. UnderH
0
,
(n
1
R
1
)=(n
1
1)
(n
2
R
2
)=(n
2
1)
F
n
1
1;n
2
1
: (5.33)
Monte Carlo simulations on test for equality of concentration
We use a Monte Carlo simulation to investigate the detection power of the test for equal-
ity of concentration in the last section. We first randomly generated two sets of circular
data from the von Mises distribution with different
1
and
2
, and then used the above
method to test for equality of concentration. By repeating the simulation and test for
1000 times, we computed the probability of rejecting the null hypothesis H
0
:
1
=
2
.
We repeated such procedure by varying the value of
2
in the simulation, and generated
the change of probability of rejectionH
0
as a function of
2
.
77
5.3 Image Registration for Group Analysis
Because electrodes were implanted at different locations in each patient, to perform
group level analysis we register the subjects’ brains into a common reference brain, and
estimate the electrode positions on the common brain.
5.3.1 Location of electrodes
The positions of electrodes were initially marked on the postoperative computed tomog-
raphy (CT) at Yale University. CT scans were then registered to postoperative MRI
scans, and the electrode positions were transformed into the postoperative MRI coor-
dinates. Since we are interested in the electrodes implanted the cortical surface, it is
important to align the cortical structure among all the subjects by surface registration.
Therefore, we need to extract the brain surface from the 3D MRI volume images. Unfor-
tunately, the subjects’ brain can be deformed due to the process of implanting electrodes
and change of pressure when the skull is open. As shown in figure 5.5, the postoperative
MRI can be very different from a normal brain, and cause failure of surface extraction
software. We instead use additional registration steps to address this problem. We first
use non-rigid 3D volume registration to register the postoperative MRI with the preoper-
ative MRI for the same subject, then perform the surface extraction on the preoperative
MRI. The procedure of the registration is in figure 5.3.1.
We use Automated Image Registration (AIR) software (University of California,
Los Angeles) for 3D volume registration between preoperative MRI and postoperative
MRI. Although the registration is done for intra-subject alignment, we use the 10th
order polynomials based 3D nonlinear registration provided in AIR because of the brain
deformation in postoperative condition.
We then use the software FreeSurfer [Dale et al., 1999,Fischl et al., 1999] to extract
the cortical surface from preoperative MRI, and register the 3D cortical surfaces from
78
Figure 5.5: The preoperative (left) and the postoperative (right) MRI of one subject.
Note the postoperative MRI is significantly deformed.
Figure 5.6: The scheme for registration and surface parcellation
all the subjects to a common brain. In FreeSurfer, the cortical surface of each subject
is first mapped to a sphere with a parameterization for the subsequent cortical align-
ment, and the registration of the brain of each subject is then computed with respect to
the intermediate spherical representation [Fischl et al., 1999]. Because FreeSurfer uses
variance of the convexity, we are able to get consistent folding structure such as the
central sulcus, the sylvian fissure, and apply the the same registration to extrapolate the
electrode positions on the cortical surface.
79
5.3.2 Surface parcellation
Although the fact that electrodes from all the subjects have a broad coverage of the
cortical surface, indicates the possibility of doing group analysis, we note that there are
almost no overlapping electrodes after the mapping. It is, therefore, impractical to do
the group analysis on the each electrode across subjects. In order to find the group effect
of the working memory task, we have to conduct analysis on a higher level. We divide
the cortical surface into a set of gyral based regions of interest, and perform the analysis
on the regional level by combining information from electrodes on the same region of
interest. We use the automated parcellation and labellings provided in FreeSurfer to
generate the regions of interest [Desikan et al., 2006, Fischl et al., 2004]. Each cortical
hemisphere was subdivided into 34 regions based on a dataset of manually labeled 40
brain using the intermediate spherical representation. Each of the parcellated cortical
region will contain reasonable amount of electrodes; therefore, we can use these regions
to conduct further group analysis.
5.4 Group Analysis by Combiningp-values
We discussed the parametric test for equality of phase synchrony in section 5.2. For
a single subject study, we conduct a test for phase locking between different memory
loads for each pair of electrodes at certain frequency band, and find the p-value for
this test. We can then find the statistically significant p-value over the whole brain for
that subject. But for group analysis, we need to be able to combine these p-values
from multiple subjects to find meaningful group level statistical significance. Here we
propose to use a method of combiningp-values and the parcellated regions of interest
from section 5.3.
80
5.4.1 Methods for combiningp-values
Summarizing evidence from multiple hypothesis tests has always been an important
topic in statistics. When it is impossible to generate a single combined hypothesis,
combination of p-values is used to summarize all the single hypothesis test. These
methods of combiningp-values are very popular for data fusion or meta-analysis.
Suppose we have k independent tests (k 2), and the null hypothesis of each is
H
01
; ;H
0k
. The combined null hypothesisH
0
will be that all of the individual null
hypothesis is true, and the alternative hypothesis is that at least one of them is not true.
Many methods have been proposed for combining p-values into a single hypothesis
test [Loughin, 2004]. Most methods are based on the fact that under the null hypothesis
a p-value from an continuous test statistic has a uniform distribution on the interval
[0,1]. Different combining methods will then give statistics with different distributions
and reject regions. The most popular methods are the minimum p-value method, the
maximum p-value method, Fisher’s combined p-value method [Fisher, 1996], inverse
normal method and Lipt´ ak’s combinedp-value method [Lipt´ ak, 1959].
Assume thep-values for the individual hypothesis testsH
01
; ;H
0k
arep
1
; ;p
k
.
The different combinedp-value are defined as follows:
(1) Maximalp-value:
p = max (p
1
; ;p
k
) (5.34)
The combined p-value is the maximum of all p-values of the individual tests. In this
case, we will need all of thep-values to be smaller than the designed level in order to
reject the null hypothesisH
0
. The rejection region forH
0
is max (p
1
; ;p
k
)<.
(2) Minimalp-value:
p = min (p
1
; ;p
k
) (5.35)
81
The minimum combinedp-value allows us to reject the null hypothesis as long as one
of the individualp-values is smaller than the designed level. The rejection region for
H
0
is min (p
1
; ;p
k
)< in this case.
(3) Fisher’s combinedp-value:
Fisher’s combined p-value is based on the fact that under the null hypothesis,
2 ln(p
i
) follows the
2
distribution with 2 degrees of freedom. If the tests are indepen-
dent, the following statistics
X
2
=2
k
X
i=1
ln(p
i
) (5.36)
will have a chi-square distribution with 2k degrees of freedom under the null hypothesis.
We can then find thep-value for the combined hypothesis by using the PDF of chi-square
distribution. When the individual p-value is small,X
2
will be large, we will reject the
null hypothesis that each of the individual null hypothesis is true. The rejection region
forH
0
in Fisher’s combinedp-value method isp
1
p
k
< c. Herec is determined by
the desired value of.
(4) Inverse normal method for combinedp-value:
Because the p-value follows uniform distribution under null hypothesis, we can
transform thep-value into a Gaussian distributed random variable,i.e.

1
(p
i
)N(0; 1); (5.37)
where (x) is the inverse of the cumulative distribution function (CDF) of standard
normal distribution. Therefore,
Z =
1
p
k
k
X
i=1

1
(p
i
); (5.38)
82
will also follow a standard normal distribution N(0; 1). We can then use the normal
distribution to find thep-value of the combined hypothesis. The rejection region forH
0
is then
1
p
k
k
X
i=1
Z
1p
i
>Z

: (5.39)
(5) Lipt´ ak’s combinedp-value:
Lipt´ ak [Lipt´ ak, 1959] extended the inverse normal method above by allowing each
test to have different weightsw
i
, where
k
X
i=1
w
2
i
= 1: (5.40)
Eq.5.38 becomes
Z =
k
X
i=1
w
i

1
(p
i
): (5.41)
It also follows the standard normal distribution under the null hypothesis. We can com-
pute the combinedp-value in the same way as the inverse normal method. Similarly, the
rejection region forH
0
is then
k
X
i=1
w
i
Z
1p
i
>Z

: (5.42)
We can see that inverse normal method is actually a special case for Lipt´ ak’s com-
binedp-value method where
w
i
=
1
p
k
(5.43)
Several studies have shown that there is no best combination ofp-values [Birnbaum,
1954, Won et al., 2009]. [Owen, 2009] suggested that the inverse normal combining
function does well in problems where evidence against the combined null is spread
among more than a small fraction of the individual tests, or when the total evidence
is weak. Fisher’s combined p-value method works best when the evidence is at least
83
moderately strong and is concentrated in a relatively small fraction of the individual
tests.
5.4.2 Combiningp-values for each region
In the case of analyzing brain interactions from intracranial EEG, instead of having
the full data, the estimated current density map throughout the cortex in noninvasive
EEG/MEG, or 3D volume image over the entire brain volume in fMRI, we only have
a very sparse measurement scattering on the cortex. Therefore, we cannot perform the
group analysis solely based on the inter-subject registration as in MEG or fMRI. We
have to base the analysis on a set of large regions of interest.
For each parcellated regions, we are interested in how the regions interact with oth-
ers over all the subjects. We combine the p-values from test for the phase synchrony
between the electrodes across the regions to get the interactions between regions. For
example, to compute connectivity between frontal gyrus and central gyrus, we find all
the subjects which have electrodes in both regions. For each cross-region electrodes
pair, we test for significance of difference in phase synchrony between memory condi-
tions and compute thep-value for the test. We then combine all thep-values from all
these electrodes pairs across the subjects, and generate ap-value for interaction between
these two regions.
5.4.3 Controling false discovery rate for multiple comparisons
After combiningp-values for the group level regions of interest, we get onep-value for
each pair of interactions. Finding which of these group levelp-values are significant is
a multiple comparison problem [Nichols and Hayasaka, 2003] [Pantazis et al., 2003].
To correct for multiple comparison, we use a false discovery rate (FDR) procedure to
control the expected proportion of errors among the rejected hypothesis [Benjamini and
84
Hochberg, 1995] [Benjamini and Yekutieli, 2001] [Genovese et al., 2002]. For example,
FDR at a level = 0:05 means that 5 of these will be false positives if we reject 100
null hypotheses. We choose FDR because when many of our hypotheses are rejected,
the error from a single erroneous rejection is not crucial for drawing conclusions about
the group level effect.
5.5 Network analysis using graph theory
Once the interaction network is constructed from functional data, we can use graph
theory to characterize the organization of the resulting network. A graph G=(V , E) is
defined as a set of nodes (V) and the edges (E) between them. A graph can be defined
as a directed graph where the edges start from one node and ends at the other, or an
undirected graph when the edges are just links between two nodes without directionality.
A graph can also be a binary (unweighted) graph or weighted graph. Graph theory based
network analysis has been used in many fields of science and engineering. It has also
been used in cortical network analysis [Bullmore and Sporns, 2009] for both structural
[Gong et al., 2008,Iturria-Medina et al., 2008] and functional network analysis [Deuker
et al., 2009,van den Heuvel et al., 2008,Honey et al., 2007,Reijneveld et al., 2007,Sporns
et al., 2007, Stam et al., 2009]. In this work we focus on the application of graph theory
in functional network analysis.
For functional interaction networks, we define the nodes as the regions of interest
on the cortical surface, and treat the functional interactions between them as the edges
between nodes. The graph topology can be quantitatively characterized using different
representations and measures.
85
5.5.1 Representation of graph
A graph can be represented by its adjacency matrix. The adjacency matrix A of a graph
is a matrix with rows and columns labeled by graph nodes, with the entryA
ij
at position
(i;j) representing the edge between theith node andjth node. For binary graphs, the
entry is 1 if there is an edge from nodei to nodej and is 0 if there is no edge from node
i to nodej; for weighted graph, the entry is assigned as the weight of the edge fromi to
j. The adjacency matrix is always symmetric for undirected graph.
5.5.2 Graph topology measure
Node degrees of the graph
The degree of a node is the most fundamental measure in graph theory. It is defined as
the number of connections that link it to the rest of the network. For the weighted graph,
the degree is the sum of the weights of all the edges connected to the node.
D(i) =
X
j2V
j6=i
A
ij
(5.44)
Clustering coefficient
The clustering coefficient is a measure of the likelihood that other nodesj connected to
the nodei will also be connected to each other [Watts and Strogatz, 1998]. It quanti-
fies the number of connections that exist between the nearest neighbors of a node as a
proportion of the maximum number of possible connections.
C
i
=
P
k6=i
P
j6=i
j6=k
A
ik
A
ij
A
kj
P
k6=i
P
j6=i
j6=k
A
ik
A
ij
(5.45)
86
The clustering coefficient for the whole graph G is defined as the average of clustering
coefficients around each node.
C =
1
n
X
i2V
C
i
(5.46)
Random networks generally have low average clustering whereas complex networks
have high clustering coefficients.
Path length
Path length indicates the minimum number of edges to traverse from nodei to nodej.
For the weighted graph, it is defined as the inverse of the edge weight, i.e.
L
ij
=
1
A
ij
(5.47)
The length of a weighted path between any two nodes is defined as the sum of the
lengths of the edges of this path. The path of the whole graph is then defined as the
average weighted path length.
Centrality
The centrality (or betweenness centrality) of a node measures how many of the shortest
paths between all other node pairs in the graph pass through it. A high centrality node
is the one crucial to efficient communication, therefore it is the most ”‘central”’ node of
the graph. Mathematically, it is defined as
b
i
=
X
k;j2V
k6=j

kj
(i)

kj
(5.48)
where
kj
is the total number of shortest paths fromk toj and
kj
(i) is the number of
the paths passing through nodei [DallAsta et al., 2006] .
87
Hubs of graph
Hubs are the nodes with high degree or high centrality.
Community Structure
Some networks, such as the functional interaction networks, exhibit a community struc-
ture. Such networks (graphs) are made up of clusters of strongly inter-connected nodes
and can be partitioned into subnetworks. Several methods have been used for graph
partitioning. One approach is based on the optimization of a modularity parameter Q
defined as follows,
Q =
1
n
X
i;j

A
ij

k
i
k
j
2m

group
ij
(5.49)
wherem =
1
2
P
ij
A
ij
is the total number of the edges in the graph,k
i
is the degree of
nodei and
group
ij
equals 1 if nodesi andj are in the same community and 0 otherwise.
Q measures the difference between the fraction of the edges connecting nodes within
communities and the same fraction for a random network with the same partition. By
maximizing the modularity of the graph (subgraph), we can iteratively partition the net-
work into subnetworks until no partitions can increase the modularity of the subgraphs.
5.6 Results
In this section, we apply the methods described above to the intracranial EEG data and
analyze the functional interaction network modulated by the memory load in working
memory experiments.
5.6.1 Frequencies with phase synchrony
For each subject, we compute the phase synchrony using the nonparametric method. We
first estimate the instantaneous phase for the recordings during memory retention period
88
Figure 5.7: Illustration of the narrow frequency (at 8Hz) phase locking at theta band for
subject 1
at each electrode using complex Morlet wavelet, then compute the phase locking value
by averaging the complex phase over all the trials. From the computed phase locking
values, we observe strong phase synchrony at a narrow frequency in theta band as shown
in figure 5.7. We also observe that the narrow frequency bands in which the signals are
synchronized vary among the subjects, although they are all located in theta band (48
Hz). We list the central frequency of these synchronized frequency band in table 5.1.
subject 1 2 3 4 5 6 7 8 10 11 12 13
frequency(Hz) 7 4 5 5 7 7 7 7 5 8 6 5
Table 5.1: Different synchronized frequencies for all 13 subjects
The synchronized frequencies are consistent with the theta band oscillatory power
and phase synchrony reported in [Raghavachari et al., 2006, Haenschel et al., 2007,
89
Jensen and Tesche, 2002,Wu et al., 2007,Siapas et al., 2005]. We will focus our analysis
on these frequency bands for each subject respectively.
5.6.2 Parametric test for equality of phase synchrony
For each subject, we will test the change of phase synchrony between 2 different mem-
ory load conditions, load 1 and load 4, at the frequency found above in table 5.1. Instead
of using complex Morlet wavelet again to find the instantaneous phase, we use Hilbert
transform to compute the phase in order to avoid the leakage from nearby frequencies
caused by frequency smoothing in complex Mortlet wavelet. We then get the phase dif-
ferences for each memory load, and use them as the circular data in testing for equality
of concentration. We plot relative phase data for one pair of electrodes in figure 5.8 as
an example and show the fitting of data using von Mises distribution.
Figure 5.8: Histogram of phase difference between two electrodes for subject 1. The
red curve is the PDF of the fitted von Mises distribution
We apply the parametric test for equality of phase synchrony between memory load
1 and 4 for each subject, and show the result of the adjacency matrix for each subject
90
after thresholding based on FDR in figure 5.9, 5.10 and 5.11 for either two sided and
one sided test.
Figure 5.9: Single-subject adjacency matrices of the memory load modulated interaction
networks for two sided test (H
1
:
1
6=
4
) after FDR. White pixels indicate significant
interactions between the electrodes corresponding to the row and column; Black pixels
indicate no significant interactions.
Comparing figure 5.9, 5.10 and 5.11, we notice that the one-sided test
1
<
4
has more significant interactions (the white dots in the figure) than the one-sided test

1
>
4
for most subjects. This means that the interactions are more likely to increase
as the working memory task becomes harder.
91
Figure 5.10: Adjacency matrix of single subject for the memory load modulated inter-
actions network for one sided test (H
1
:
1
<
4
) after FDR.
5.6.3 Image registration
The result of surface registration with the positions of the electrodes on the common
brain surface are shown in figure 5.12. The electrodes are represented by the colored
dots and the different colors indicate the subjects that the electrodes are from.
The result of the parcellation is represented by the coloring of the brain regions in
the same figure. Each of 68 parcellated regions is represented by a unique color. The
result shows that most of the parcellated cortical region contains a reasonable amount of
electrodes, which makes the group analysis on the regions possible.
92
Figure 5.11: Adjacency matrix of single subject for the memory load modulated inter-
actions network for one sided test (H
1
:
1
>
4
) after FDR.
5.6.4 Group analysis result
We use the Fisher’s combinedp-value method to find thep-value of interaction among
the parcellated regions of interest and establish the significance by controlling the FDR.
The significant interactions are shown in figure 5.13
The interaction network corresponding to this adjacency matrix is shown in figure
5.14. Because it is hard to see the structure of the network when we plot all the sig-
nificant interactions, we use the graph theory based network analysis to extract more
detailed information from this network.
We first compute the strength (degree) of each node (region of interest) for the net-
work in figure 5.6.4. The strength of each region is indicated in figure 5.16 by the size
93
Figure 5.12: All the electrodes (colored dots) from 12 subjects are mapped to common
cortical surface using surface registration. Different color of the electrodes indicate they
are from different subjects.
of sphere at the center of the each region. The nodes with the most connections with
others are called the hubs defined by degrees. We observe the following hubs regions
in the network : left superior frontal, right superior frontal, right precentral, left middle
temporal, right caudal middle frontal and left lateral occipital.
In figure 5.17, we show how these important regions of interest are connected to
other nodes in the graph. From the results of the individual regions, we observe that
most of those regions are active for both intra-hemisphere and in inter-hemisphere. In
addition, these important regions are also more likely to connect each other than these
other regions.
94
Figure 5.13: Adjacency matrix of parcellated regions of interest for the memory
load modulated interactions network for one sided test (H
1
:
1
<
4
) after FDR at
= 0:001. White pixels indicate significant interactions between the electrodes corre-
sponding to the row and column; Black pixels indicate no significant interactions.
5.6.5 Graph partition result
We applied the modularity based graph partition method to the group-level interaction
network. We found 3 subnetworks through the partition procedure by maximizing the
modularity of the graph. The subnetworks are shown in figure 5.18. The first subnetwork
covers the left frontal region, the second subnetwork is mostly on the right frontal region
and the bilateral sensory motor area, and the third subnetwork on the bilateral parietal
and occipital areas.
95
Top view Bottom view
Left view
Right view
Figure 5.14: The memory task modulated significant network found on group level using
Fisher’s combinedp-value method.
5.7 Discussion
We used a parcellation of the cortical surfaces based on gyral anatomy , and grouped
the interactions between electrodes from two different regions together by combining
the p-value computed using the parametric test of phase synchrony. However, we want
to note that our group analysis procedures are not limited to the particular parcellation
scheme of gyri-based regions of interest. We can use other parcellation schemes with
smaller regions as long as we have enough electrodes within each region of interest. For
example, we can use the finer parcellation described in [Fischl et al., 2004], or use other
manually defined parcellation scheme by defining the regions of interest.
96
Figure 5.15: Strength(degree) of each node is computed as the sum of all connection
weights for the node.
We used the parametric test for finding the change of phase synchrony between dif-
ferent memory conditions. Parametric test has its computational advantages, especially
for large multi-subject data as we used in this work. However, we note that nonpara-
metric tests, such as bootstrapping or permutation test, can also be used to test the dif-
ferences of phase synchrony. The nonparametric tests will be more accurate when the
underlying probability assumption of von Mises distribution can not be guaranteed.
In addition, our methods are not limited to the specific interaction measure of phase
synchrony used in this chapter. we can use other interaction measures described in
chapter 3. Because different interaction measures reflect different aspects of functional
interactions, we will not get exactly the same network as we get for phase synchrony.
In this chapter, we limited our investigation to the narrow frequency bands within
theta band. Our group analysis can be also be used for interactions over other frequency
97
Top view
Bottom view
Left view
Right view
Figure 5.16: The strengths (degrees) of regions in a group level memory modulated
interaction network are indicated by the sizes of the spheres at the center of each region
of interest.
bands as well. Our choice of frequency is basically based on the following two consid-
erations: First, theta activities have been widely accepted as the active frequency band
for working memory task. Second, our preliminary investigation of the phase synchrony
over all frequency bands shows the dominated narrow frequency within theta band. We
picked the narrow frequency individually for each subject. If one observes a similar
narrow frequency synchrony in other frequency, the same procedure can be used for that
frequency band as well.
98
Figure 5.17: The connection of 8 individual nodes with the highest degree. Among the
68 regions of interest, these 8 regions are those connected to other regions.
99
Subnetwork 1 Subnetwork 2 Subnetwork 3
Figure 5.18: The community structure of the group level interaction network. 3 subnet-
works are found by maximizing the modularity.
5.8 Summary
We have described a novel set of group analysis procedures for analyzing interactions
using intracranial recording based brain image registration, combining statistical sig-
nificance and parametric test of phase synchrony. We applied them on real intracranial
EEG recordings from subjects performing a Sternberg working-memory task. The group
analysis results we got are partially consistent with the working memory related inter-
actions.
100
Chapter 6
Conclusion and Future Work
6.1 Conclusion
We have described the nulling beamformer inverse method for accurate estimation of
cortical interaction networks by suppressing cross-talk effects. We reviewed the com-
monly used functional interaction measures and showed their sensitivity to linear mixing
among cortical sources. We presented the nulling beamformer methods as a modifica-
tion of the standard LCMV beamformer by explicitly forcing nulling constraints on other
sources when estimating the source signals. In this way, we control the linear mixing
among several cortical sources. We first described the nulling beamformer for dipole
sources, then extended it for patch sources using eigenvector constraints.
We have used a novel simulation scheme to show the effectiveness of the nulling
beamformer. We simulated MEG cortical sources with intracranially measured local
field potentials (LFP) as better representations of true neuronal sources. We have com-
puted interactions before and after linear mixing and evaluate the robustness of our
method to linear cross-talk. We measured the performance of nulling beamformer in
conjunction with different interaction measures and demonstrated that the method is
superior to other inverse imaging methods, such as minimum-norm imaging and con-
ventional LCMV beamforming, in estimating cortical interactions.
Even though the nulling beamformer method is effective in addressing the cross-
talk problem in cases where we approximately know the locations of the interfering
sources, the information is not always known in real data experiments. Therefore, the
nulling beamformer should be better viewed as a step towards ameliorating the linear
101
mixing problem, rather than a complete solution. Careful interpretation of results and
fine tuning of the eigenvector constraints are important for a reliable estimation of the
interacting cortical sources.
Interaction measures are insensitive to cross-talk, such as the imaginary part of the
complex coherence [Nolte et al., 2004], providing another way of exploring cortical
interaction networks. However, they suffer from considerable limitations. For exam-
ple, imaginary coherence is sensitive to time delay, and therefore completely blind to
perfectly synchronous signals with zero time lag. Such a limitation does not exist with
the nulling beamformer. For this reason, methods such as the nulling beamformer and
imaginary coherence are complimentary, because they provide a different view of ongo-
ing cortical networks, and they should both be considered when exploring experimental
MEG and EEG data.
On the other hand, intracranial EEG recordings are free of the linear cross-talk, but
the challenge of using intracranial data is the sparse and distinct locations of implanted
electrodes, which makes the group analysis from multiple subjects difficult. We have
presented a set of novel analysis procedures for assessing functional interactions among
multiple subjects using intracranial EEG data. We used a parametric model of phase
synchrony to test for the equality of interactions under different memory load condi-
tions in a working memory experiment. We have shown by using a von Mises distribu-
tion based probability model, we can use the circular statistics to efficiently compute the
significant of phase synchrony without computationally expensive resampling methods.
To prepare for multi-subject analysis, we have used 3D volumetric registration and sur-
face registration to map the implanted electrodes onto a commonly cortical surface. We
have proposed a group analysis for these intracranial recordings based on a gyri-based
regions of interest parcellation of the cortical surfaces. We grouped the interactions
102
between electrodes from two different regions together by combining thep-value com-
puted using the parametric test. We found the significant interactions network among
all the regions of interest under the control of false discovery rate (FDR).
We have also used graph theory based network analysis on the significant interaction
network found using group analysis. We have shown that our results of the most con-
nected regions and sub-networks are consistent with the regions and interactions found
for working memory experiment in the literatures.
6.2 Future Work
Over the past 6 years, I have investigated different aspects of signal processing in func-
tional interaction analysis. I have addressed the importance of false interactions in
MEG/EEG caused by cross-talk effect, explored different possible approaches and pro-
posed the nulling beamformer as a solution to this problem. I have investigated func-
tional interactions using intracranial recordings and proposed group analysis procedures.
I have also studied other aspects of functional analysis and identified several interesting
problems for further investigation in the field of interactions analysis:
Multivariate nonlinear interaction measure
I have investigated a wide range of mathematical models for interaction analysis and
believe there is a clear need for multivariate interaction analysis. Although MV AR
model based multivariate measures has been proposed, they are limited to modeling lin-
ear interactions. A multivariate nonlinear interaction model will not only eliminate the
ambiguity caused by additional sources interacting with both signals in a bivariate inter-
action analysis, it will also allow us to capture the chaotic nature of brain activity. There
has been several attempts at this kind of measures [Allefeld and Kurths, 2004a] [Pereda
103
et al., 2005], but these approaches are not yet well established. Therefore, I hope further
investigation would help us find an effective multivariate nonlinear measure.
Analysis of phase-amplitude interaction and phase synchrony
between different frequencies using circular statistics
I have noticed the possibility of analyzing phase-amplitude interactions and phase syn-
chrony between different frequencies using circular statistics during the course of apply-
ing parametric model for phase synchrony. The phase-amplitude interaction has been
studied by [Canolty et al., 2006] using averaging methods similar to those used in non-
parametric phase synchrony detection. I believe it is possible to formulate the problem
as a circular data problem and use circular statistics in a similar way as the circular
statistic based phase synchrony analysis described in Chapter 5. For the same reason,
I believe we could extend the circular statistic to phase synchrony estimation between
two different frequencies. During my preliminary research, I have come across several
“correlation” defined for circular data, such as circular-circular correlation and circular-
linear correlation, which would possibly facilitate such analysis.
Resting state network using electrophysiological data
Resting state network analysis using fMRI has been well studied during the past years
[Cordes et al., 2002, Garrity et al., 2007, Shehzad et al., 2009, Weissenbacher et al.,
, Wu et al., 2006]. However, such resting network has not been reliably established
for electrophysiological data such EEG, MEG or intracranial EEG despite many efforts
[Cathy Nangini, 2008, Laufs, 2008, Leopold et al., , Mantini et al., 2007, Stam et al.,
2006]. I believe finding such networks using electrophyisological data is still possible,
with the help of new signals processing techniques. I also believe that the simultane-
ous EEG-fMRI recordings devices will play an important role in helping us find such
networks.
104
References
[Allefeld and Kurths, 2004a] Allefeld, C. and Kurths, J. (2004a). An approach to multi-
variate phase synchronization analysis and its application to event-related potentials.
International journal of bifurcation and chaos, 14(2):417–426.
[Allefeld and Kurths, 2004b] Allefeld, C. and Kurths, J. (2004b). Testing for phase
synchronization. International Journal of Bifurcation and Chaos in Applied Sciences
and Engineering, 14(2):405–416.
[Andrzejak et al., 2003] Andrzejak, R., Kraskov, A., Stgbauer, H., Mormann, F., and
Kreuz, T. (2003). Bivariate surrogate techniques: Necessity, strengths, and caveats.
Phys. Rev. E J1 - PRE, 68(6):066202.
[Astolfi et al., 2008] Astolfi, L., Cincotti, F., Mattia, D., De Vico Fallani, F., Tocci, A.,
Colosimo, A., Salinari, S., Marciani, M., Hesse, W., Witte, H., Ursino, M., Zavaglia,
M., and Babiloni, F. (2008). Tracking the time-varying cortical connectivity patterns
by adaptive multivariate estimators. Biomedical Engineering, IEEE Transactions on,
55(3):902–913.
[Axmacher et al., 2008] Axmacher, N., Schmitz, D. P., Wagner, T., Elger, C. E., and
Fell, J. (2008). Interactions between medial temporal lobe, prefrontal cortex, and
inferior temporal regions during visual working memory: A combined intracranial
eeg and functional magnetic resonance imaging study. J. Neurosci., 28(29):7304–
7312.
[Baccala and Sameshima, 2001] Baccala, L. and Sameshima, K. (2001). Partial
directed coherence: a new concept in neural structure determination. Biological
Cybernetics, 84(6):463–474.
[Baddeley, 1992] Baddeley, A. (1992). Working memory. Science, 255(5044):556–
559.
[Baillet et al., 2001] Baillet, S., Mosher, J., and Leahy, R. (2001). Electromagnetic
brain mapping. Signal Processing Magazine, IEEE, 18(6TY - JOUR):14–30.
105
[Benjamini and Hochberg, 1995] Benjamini, Y . and Hochberg, Y . (1995). Controlling
the false discovery rate: a practical and powerful approach to multiple testing. Journal
of the Royal Statistical Society, 57:289–300.
[Benjamini and Yekutieli, 2001] Benjamini, Y . and Yekutieli, D. (2001). The control
of the false discovery rate in multiple testing under dependency. The Annals of
Statistics, 29(4):1165–1188.
[Birnbaum, 1954] Birnbaum, A. (1954). Combining independent tests of significance.
Journal of the American Statistical Association, 49(267):559–574.
[Biswal et al., 1995] Biswal, B., Yetkin, F., Haughton, V ., and Hyde, J. (1995). Func-
tional connectivity in the motor cortex of resting human brain using echo-planar mri.
Magnetic Resonance in Medicine, 34(4):537–541.
[Bjorck and Golub, 1973] Bjorck and Golub (1973). Numerical methods for computing
angles between linear subspaces. Mathematics of computation, pages 579–594.
[Bowyer et al., 2005] Bowyer, S., Moran, J., Weiland, B., Mason, K., Greenwald, M.,
Smith, B., Barkley, G., and Tepley, N. (2005). Language laterality determined by
MEG mapping with MR-FOCUSS. Epilepsy and Behavior, 6(2):235–241.
[Bressler et al., 1993] Bressler, S., Coppola, R., and Nakamura, R. (1993). Episodic
multiregional cortical coherence at multiple frequencies during visual task perfor-
mance. Nature, 366(6451):153–156.
[Brookes et al., 2007] Brookes, M., Stevenson, C., Barnes, G., Hillebrand, A., Simp-
son, M., Francis, S., and Morris, P. (2007). Beamformer reconstruction of correlated
sources using a modified source model. NeuroImage, 34(4):1454–1465.
[Brookes et al., 2008] Brookes, M. J., Vrba, J., Robinson, S. E., Stevenson, C. M.,
Peters, A. M., Barnes, G. R., Hillebrand, A., and Morris, P. G. (2008). Optimising
experimental design for meg beamformer imaging. NeuroImage, 39(4):1788–1802.
[Brovelli et al., 2004] Brovelli, A., Ding, M., Ledberg, A., Chen, Y ., Nakamura, R., and
Bressler, S. (2004). Beta oscillations in a large-scale sensorimotor cortical network:
Directional influences revealed by granger causality. PNAS, 101(26):9849–9854.
[Bullmore and Sporns, 2009] Bullmore, E. and Sporns, O. (2009). Complex brain net-
works: graph theoretical analysis of structural and functional systems. Nat Rev
Neurosci, 10(3):186–198.
[Canolty et al., 2006] Canolty, R., Edwards, E., Dalal, S., Soltani, M., Nagarajan, S.,
Kirsch, H., Berger, M., Barbaro, N., and Knight, R. (2006). High gamma power
is phase-locked to theta oscillations in human neocortex 10.1126/science.1128115.
Science, 313(5793):1626–1628.
106
[Carter, 1987] Carter, G. (1987). Coherence and time delay estimation. Proceedings of
the IEEE, 75(2):236–255.
[Carter et al., 1973] Carter, G., Knapp, C., and Nuttall, A. (1973). Estimation of the
magnitude-squared coherence function via overlapped fast fourier transform process-
ing. Audio and Electroacoustics, IEEE Transactions on, 21(4):337–344.
[Cathy Nangini, 2008] Cathy Nangini, F.T., S. J. G. (2008). A novel method for inte-
grating meg and bold fmri signals with the linear convolution model in human pri-
mary somatosensory cortex. Human Brain Mapping, 29(1):97–106.
[Cohen, 1972] Cohen, D. (1972). Magnetoencephalography: detection of the brain’s
electrical activity with a superconducting magnetometer. Science, 175(4022):664.
[Cordes et al., 2002] Cordes, D., Haughton, V ., Carew, J., Arfanakis, K., and Maravilla,
K. (2002). Hierarchical clustering to measure connectivity in fmri resting-state data.
Magnetic Resonance Imaging, 20(4):305–317.
[Dalal et al., 2006] Dalal, S., Sekihara, K., and Nagarajan, S. (2006). Modified beam-
formers for coherent source region suppression. Biomedical Engineering, IEEE
Transactions on, 53(7TY - JOUR):1357–1363.
[Dale et al., 1999] Dale, A., Fischl, B., and Sereno, M. (1999). Cortical surface-based
analysis I. Segmentation and surface reconstruction. Neuroimage, 9(2):179–194.
[Dale et al., 2000] Dale, A., Liu, A., Fischl, B., Buckner, R., Belliveau, J., Lewine, J.,
and Halgren, E. (2000). Dynamic statistical parametric mapping: combining fmri
and meg for high-resolution imaging of cortical activity. Neuron, 26(1):55–67.
[DallAsta et al., 2006] DallAsta, L., Barrat, A., Barthelemy, M., and Vespignani, A.
(2006). Vulnerability of weighted networks. Journal of Statistical Mechanics: Theory
and Experiment, 2006:P04006.
[Darvas et al., 2004] Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., and Leahy, R.
(2004). Mapping human brain function with meg and eeg: methods and validation.
NeuroImage, 23(Supplement 1Darvas04):S289–S299.
[de Peralta-Menendez and Gonzalez-Andino, 1998] de Peralta-Menendez, R. and
Gonzalez-Andino, S. (1998). A critical analysis of linear inverse solutions to
the neuroelectromagnetic inverse problem. IEEE Transactions on Biomedical
Engineering, 45(4):440–448.
[Desikan et al., 2006] Desikan, R., S´ egonne, F., Fischl, B., Quinn, B., Dickerson, B.,
Blacker, D., Buckner, R., Dale, A., Maguire, R., Hyman, B., et al. (2006). An auto-
mated labeling system for subdividing the human cerebral cortex on MRI scans into
gyral based regions of interest. Neuroimage, 31(3):968–980.
107
[Deuker et al., 2009] Deuker, L., Bullmore, E. T., Smith, M., Christensen, S., Nathan,
P. J., Rockstroh, B., and Bassett, D. S. (2009). Reproducibility of graph metrics of
human brain functional networks. NeuroImage, 47(4):1460–1468.
[Ding et al., 2000] Ding, M., Bressler, S., Yang, W., and Liang, H. (2000). Short-
window spectral analysis of cortical event-related potentials by adaptive multivari-
ate autoregressive modeling: data preprocessing, model validation, and variability
assessment. Biol Cybern, 83(1):35–45.
[Ermer et al., 2001] Ermer, J. J., Mosher, J. C., Baillet, S., and Leahy, R. M. (2001).
Rapidly recomputable eeg forward models for realistic head shapes. Physics in
Medicine and Biology, 46(4):1265–1281.
[Fingelkurts et al., 2005] Fingelkurts, A., Fingelkurts, A., and Kahkonen, S. (2005).
Functional connectivity in the brain–is it an elusive concept? Neuroscience &
Biobehavioral Reviews, 28(8):827–836.
[Fischl et al., 1999] Fischl, B., Sereno, M., and Dale, A. (1999). Cortical surface-based
analysis II: Inflation, flattening, and a surface-based coordinate system. Neuroimage,
9(2):195–207.
[Fischl et al., 2004] Fischl, B., van der Kouwe, A., Destrieux, C., Halgren, E., Segonne,
F., Salat, D., Busa, E., Seidman, L., Goldstein, J., Kennedy, D., et al. (2004). Auto-
matically parcellating the human cerebral cortex. Cerebral Cortex, 14(1):11.
[Fisher, 1996] Fisher, N. (1996). Statistical analysis of circular data. Cambridge Univ
Pr.
[Fox and Raichle, 2007] Fox, M. and Raichle, M. (2007). Spontaneous fluctuations
in brain activity observed with functional magnetic resonance imaging. Nat Rev
Neurosci, 8(9):700–711.
[Gardner, 1992] Gardner, W. A. (1992). A unifying view of coherence in signal pro-
cessing. Signal Process., 29(2):113–140.
[Garrity et al., 2007] Garrity, A., Pearlson, G., McKiernan, K., Lloyd, D., Kiehl,
K., and Calhoun, V . (2007). Aberrant ”default mode” functional connectivity in
schizophrenia. American Journal of Psychiatry Am J Psychiatry, 164(3):450–457.
[Genovese et al., 2002] Genovese, C., Lazar, N., and Nichols, T. (2002). Threshold-
ing of statistical maps in functional neuroimaging using the false discovery rate.
NeuroImage, 15:870–878.
[Geweke, 1982] Geweke, J. (1982). Measurement of linear dependence and feed-
back between multiple time series. Journal of the American Statistical Association,,
77(378):304–313.
108
[Gong et al., 2008] Gong, G., He, Y ., Concha, L., Lebel, C., Gross, D. W., Evans, A. C.,
and Beaulieu, C. (2008). Mapping anatomical connectivity patterns of human cere-
bral cortex using in vivo diffusion tensor imaging tractography. Cereb. Cortex, page
102.
[Granger, 1969] Granger, C. (1969). Investigating causal relations by econometric
models and cross-spectral methods. Econometrica, 37(3):424–438.
[Gross et al., 2001] Gross, J., Kujala, J., Hamalainen, M., Timmermann, L., Schnitzler,
A., and Salmelin, R. (2001). Dynamic imaging of coherent sources: Studying neural
interactions in the human brain 10.1073/pnas.98.2.694. Proceedings of the National
Academy of Sciences, 98(2):694–699.
[Hadjipapas et al., 2005] Hadjipapas, A., Hillebrand, A., Holliday, I., Singh, K.,
and Barnes, G. (2005). Assessing interactions of linear and nonlinear neuronal
sources using meg beamformers: a proof of concept. Clinical Neurophysiology,
116(6):1300–1313.
[Haenschel et al., 2007] Haenschel, C., Uhlhaas, P. J., and Singer, W. (2007).
Synchronous oscillatory activity and working memory in schizophrenia.
Pharmacopsychiatry, 40:S54–S61.
[Hamalainen et al., 1993] Hamalainen, M., Hari, R., Ilmoniemi, R., Knuutila, J., and
Lounasmaa, O. (1993). Magnetoencephalographytheory, instrumentation, and appli-
cations to noninvasive studies of the working human brain. Reviews of modern
Physics, 65(2):413–497.
[Honey et al., 2007] Honey, C. J., Ktter, R., Breakspear, M., and Sporns, O. (2007).
Network structure of cerebral cortex shapes functional connectivity on multiple time
scales. Proceedings of the National Academy of Sciences, 104(24):10240–10245.
[Huang et al., 1999] Huang, M. X., Mosher, J. C., and Leahy, R. M. (1999). A sensor-
weighted overlapping-sphere head model and exhaustive head model comparison for
meg. Physics in Medicine and Biology, 44(2):423–440.
[Hui and Leahy, 2006] Hui, H. and Leahy, R. (2006). Linearly constrained meg beam-
formers for mvar modeling of cortical interactions. In Biomedical Imaging: Nano to
Macro, 2006. 3rd IEEE International Symposium on, pages 237–240–.
[Hui et al., 2010] Hui, H. B., Pantazis, D., Bressler, S. L., and Leahy, R. M.
(2010). Identifying true cortical interactions in meg using the nulling beamformer.
NeuroImage, 49(4):3161–3174.
[Iturria-Medina et al., 2008] Iturria-Medina, Y ., Sotero, R. C., Canales-Rodrguez, E. J.,
Alemn-Gmez, Y ., and Melie-Garca, L. (2008). Studying the human brain anatomical
network via diffusion-weighted mri and graph theory. NeuroImage, 40(3):1064.
109
[Jammalamadaka and Sengupta, 2001] Jammalamadaka, S. and Sengupta, A. (2001).
Topics in circular statistics. World Scientific Pub Co Inc.
[Jean-Philippe Lachaux, 1999] Jean-Philippe Lachaux, E.R., J. M. F. J. V . (1999). Mea-
suring phase synchrony in brain signals. Human Brain Mapping, 8(4):194–208–.
[Jeffs et al., 1987] Jeffs, B., Leahy, R., and Singh, M. (1987). An evaluation of
methods for neuromagnetic image reconstruction. Biomedical Engineering, IEEE
Transactions on, BME-34(9Jeffs87):713–723.
[Jensen and Tesche, 2002] Jensen, O. and Tesche, C. D. (2002). Frontal theta activity
in humans increases with memory load in a working memory task. European Journal
of Neuroscience, 15(8):1395–1399.
[Kahana, 2006] Kahana, M. J. (2006). The cognitive correlates of human brain oscilla-
tions. J. Neurosci., 26(6):1669–1672.
[Kaminski et al., 2001] Kaminski, M., Ding, M., Truccolo, W., and Bressler, S. (2001).
Evaluating causal relations in neural systems: Granger causality, directed trans-
fer function and statistical assessment of significance. Biological Cybernetics,
85(2):145–157.
[Klimesch, 1999] Klimesch, W. (1999). Eeg alpha and theta oscillations reflect cog-
nitive and memory performance: a review and analysis. Brain Research Reviews,
29(2-3):169–195.
[Komssi et al., 2004] Komssi, S., Huttunen, J., Aronen, H. J., and Ilmoniemi, R. J.
(2004). Eeg minimum-norm estimation compared with meg dipole fitting in the local-
ization of somatosensory sources at s1. Clinical Neurophysiology, 115(3):534–542.
[Kumar et al., 2009] Kumar, S., Rao, S. L., Chandramouli, B. A., and Pillai, S. V .
(2009). Reduction of functional brain connectivity in mild traumatic brain injury
during working memory. Journal Of Neurotrauma, 26(5):665–675.
[Kus et al., 2004] Kus, R., Kaminski, M., and Blinowska, K. (2004). Determination
of eeg activity propagation: pair-wise versus multichannel estimate. Biomedical
Engineering, IEEE Transactions on, 51(9Kus04):1501–1510.
[Lachaux et al., 1999] Lachaux, J., Rodriguez, E., Martinerie, J., and Varela, F. (1999).
Measuring phase synchrony in brain signals. Human Brain Mapping, 8(4):194–208.
[Lachaux et al., 2003] Lachaux, J., Rudrauf, D., and Kahane, P. (2003). Intracranial
EEG and human brain mapping. Journal of Physiology-Paris, 97(4-6):613–628.
110
[Lachaux et al., 2007] Lachaux, J.-P., Baillet, S., Adam, C., Ducorps, A., Jerbi, K.,
Bertrand, O., Garnero, L., and Martinerie, J. (2007). A simultaneous MEG and
intracranial EEG study of task-related brain oscillations. International Congress
Series New Frontiers in Biomagnetism. Proceedings of the 15th International
Conference on Biomagnetism, Vancouver, BC, Canada, August 21-25, 2006,
1300:421–424.
[Laufs, 2008] Laufs, H. (2008). Endogenous brain oscillations and related networks
detected by surface eeg-combined fmri. Human Brain Mapping, 29(7):762–769.
[Le Van Quyen et al., 2001] Le Van Quyen, M., Foucher, J., Lachaux, J.-P., Rodriguez,
E., Lutz, A., Martinerie, J., and Varela, F. (2001). Comparison of hilbert transform
and wavelet methods for the analysis of neuronal synchrony. Journal of Neuroscience
Methods, 111(2):83–98.
[Leopold et al., ] Leopold, D., Murayama, Y ., and Logothetis, N. Very slow activ-
ity fluctuations in monkey visual cortex: Implications for functional brain imaging.
Cerebral Cortex, 13:422–433.
[Lipt´ ak, 1959] Lipt´ ak, T. (1959). On the combination of independent tests.
[Loughin, 2004] Loughin, T. M. (2004). A systematic comparison of methods for com-
bining p-values from independent tests. Computational Statistics & Data Analysis,
47(3):467–485.
[Lutkenhoner et al., 1996] Lutkenhoner, B., Greenblatt, R., M.Hamalainen, Mosher,
J., Scherg, M., Tesche, C., and Sosa, P. V . (1996). Comparison between different
approaches to the biomagnetic inverse problem - workshop report.
[Mantini et al., 2007] Mantini, D., Perrucci, M., Del Gratta, C., Romani, G., and Cor-
betta, M. (2007). Electrophysiological signatures of resting state networks in the
human brain. Proceedings of the National Academy of Sciences, 104(32):13170–
13175.
[Mardia and Jupp, 2000] Mardia, K. and Jupp, P. (2000). Directional statistics. Wiley
New York.
[Maudgil, 2003] Maudgil, D. (2003). Brain imaging in epilepsy. Remedica Pub Ltd.
[Meltzer et al., 2007] Meltzer, J. A., Negishi, M., Mayes, L. C., and Constable, R. T.
(2007). Individual differences in eeg theta and alpha dynamics during working
memory correlate with fmri responses across subjects. Clinical Neurophysiology,
118(11):2419–2436.
111
[Meltzer et al., 2008] Meltzer, J. A., Zaveri, H. P., Goncharova, I. I., Distasio, M. M.,
Papademetris, X., Spencer, S. S., Spencer, D. D., and Constable, R. T. (2008). Effects
of working memory load on oscillatory power in human intracranial eeg. Cereb.
Cortex, 18(8):1843–1855.
[Mitchell et al., 2008] Mitchell, D. J., McNaughton, N., Flanagan, D., and Kirk, I. J.
(2008). Frontal-midline theta from the perspective of hippocampal ”theta”. Progress
in Neurobiology, 86(3):156–185.
[Mosher et al., 2005] Mosher, J., Baillet, S., Darvas, F., Pantazis, D., Kucukaltun-
Yildirim, E., and Leahy, R. (2005). BrainStorm electromagnetic imaging software.
2005. International Journal of Bioelectromagnetism.
[Mosher et al., 1992] Mosher, J., Lewis, P., and Leahy, R. (1992). Multiple dipole mod-
eling and localization from spatio-temporal meg data. Biomedical Engineering, IEEE
Transactions on, 39(6Mosher92):541–557.
[Newman, 2006a] Newman, M. E. J. (2006a). Finding community structure in networks
using the eigenvectors of matrices. Phys. Rev. E, 74(3):036104–19.
[Newman, 2006b] Newman, M. E. J. (2006b). Modularity and community structure in
networks. Proceedings of the National Academy of Sciences, 103(23):8577–8582.
[Nichols and Holmes, 2002] Nichols, T. and Holmes, A. (2002). Nonparametric per-
mutation tests for functional neuroimaging: A primer with examples. Human Brain
Mapping, 15(1):1–25.
[Nichols and Hayasaka, 2003] Nichols, T. E. and Hayasaka, S. (2003). Controlling the
familywise error rate in functional neuroimaging: A comparative review. Statistical
Methods in Medical Research, 12(5):419–446.
[Nolte et al., 2004] Nolte, G., Bai, O., Wheaton, L., Mari, Z., V orbach, S., and Hallett,
M. (2004). Identifying true brain interaction from eeg data using the imaginary part
of coherency. Clinical Neurophysiology, 115(10):2292–2307.
[Npflin et al., 2008] Npflin, M., Wildi, M., and Sarnthein, J. (2008). Test-retest relia-
bility of eeg spectra during a working memory task. NeuroImage, 43(4):687–693.
[Nunez et al., 1997] Nunez, P., Srinivasan, R., Westdorp, A., Wijesinghe, R., Tucker,
D., Silberstein, R., and Cadusch, P. (1997). Eeg coherency: I: statistics, refer-
ence electrode, volume conduction, laplacians, cortical imaging, and interpretation at
multiple scales. Electroencephalography and Clinical Neurophysiology, 103(5):499–
515.
112
[Ossadtchi et al., 2004] Ossadtchi, A., Baillet, S., Mosher, J., Thyerlei, D., Sutherling,
W., and Leahy, R. (2004). Automated interictal spike detection and source localiza-
tion in magnetoencephalography using independent components analysis and spatio-
temporal clustering. Clinical Neurophysiology, 115(3):508–522.
[O’Sullivan et al., 1996] O’Sullivan, B., Large, M., Woodham, B., Barrett, N., Smith,
G., Karayanidis, F., Michie, P., Watson, J., and Kavanagh, D. (1996). The role of the
human prefrontal cortex in working memory and attention – a study using positron
emission tomography. NeuroImage, 3(1, Supplement 1):S503–S503.
[Owen, 2009] Owen, A. (2009). Karl Pearsons meta-analysis revisited. The Annals of
Statistics, 37(6B):3867–3892.
[Pantazis et al., 2003] Pantazis, D., Nichols, T., Baillet, S., and Leahy, R. (2003). Spa-
tiotemporal localization of significant activation in meg using permutation tests. Inf
Process Med Imaging, 18:512–23.
[Pantazis et al., 2005] Pantazis, D., Nichols, T. E., Baillet, S., and Leahy, R. M. (2005).
A comparison of random field theory and permutation methods for the statistical
analysis of MEG data. NeuroImage, 25(2):383–394.
[Pascual-Marqui, 2002] Pascual-Marqui, R. (2002). Standardized low-resolution brain
electromagnetic tomography (sloreta): technical details. Methods Find Exp Clin
Pharmacol, 24 Suppl D:5–12.
[Pereda et al., 2005] Pereda, E., Quiroga, R., and Bhattacharya, J. (2005). Nonlinear
multivariate analysis of neurophysiological signals. Progress in Neurobiology, 77(1-
2):1–37.
[Prichard and Theiler, 1994] Prichard, D. and Theiler, J. (1994). Generating surrogate
data for time series with several simultaneously measured variables. Phys. Rev. Lett.
J1 - PRL, 73(7):951 LP – 954.
[Raghavachari et al., 2001] Raghavachari, S., Kahana, M. J., Rizzuto, D. S., Caplan,
J. B., Kirschen, M. P., Bourgeois, B., Madsen, J. R., and Lisman, J. E. (2001). Gating
of human theta oscillations by a working memory task. J. Neurosci., 21(9):3175–
3183.
[Raghavachari et al., 2006] Raghavachari, S., Lisman, J. E., Tully, M., Madsen, J. R.,
Bromfield, E. B., and Kahana, M. J. (2006). Theta oscillations in human cortex during
a working-memory task: Evidence for local generators. J Neurophysiol, 95(3):1630.
[Reddy et al., 1987] Reddy, V ., Paulraj, A., and Kailath, T. (1987). Performance analy-
sis of the optimum beamformer in the presence of correlated sources and its behavior
under spatial smoothing. Acoustics, Speech, and Signal Processing [see also IEEE
Transactions on Signal Processing], IEEE Transactions on, 35(7Reddy87):927–936.
113
[Reijneveld et al., 2007] Reijneveld, J. C., Ponten, S. C., Berendse, H. W., and Stam,
C. J. (2007). The application of graph theoretical analysis to complex networks in
the brain. Clinical Neurophysiology, 118(11):2317–2331.
[Robinson and Vrba, 1999] Robinson, S. and Vrba, J. (1999). Functional neuroimaging
by synthetic aperture magnetometry (SAM). Recent Advances in Biomagnetisml,
Japan: Tohoku Univ. Press, pages 302–305.
[Sarnthein et al., 1998] Sarnthein, J., Petsche, H., Rappelsberger, P., Shaw, G. L., and
von Stein, A. (1998). Synchronization between prefrontal and posterior association
cortex during human working memory. Proceedings of the National Academy of
Sciences of the United States of America, 95(12):7092–7096.
[Schelter et al., 2007] Schelter, B., Winterhalder, M., Timmer, J., and Peifer, M. (2007).
Testing for phase synchronization. Physics Letters A, 366(4-5):382–390.
[Scherg and V on Cramon, 1985] Scherg, M. V on Cramon, D. and V on Cramon, D.
(1985). Two bilateral sources of the late aep as identified by a spatio-temporal dipole
model. Electroencephalogr Clin Neurophysiol, 62(1):32–44.
[Schnitzler and Gross, 2005] Schnitzler, A. and Gross, J. (2005). Normal and patholog-
ical oscillatory communication in the brain. Nature Reviews Neuroscience Nat Rev
Neurosci, 6(4):285–296.
[Scott J. Peltier, 2003] Scott J. Peltier, T.A.P., D. C. N. (2003). Detecting low-frequency
functional connectivity in fmri using a self-organizing map (som) algorithm. Human
Brain Mapping, 20(4):220–226.
[Sehatpour et al., 2008] Sehatpour, P., Molholm, S., Schwartz, T., Mahoney, J., Mehta,
A., Javitt, D., Stanton, P., and Foxe, J. (2008). A human intracranial study of
long-range oscillatory coherence across a frontal-occipital-hippocampal brain net-
work during visual object processing 10.1073/pnas.0708418105. Proceedings of the
National Academy of Sciences, 105(11):4399–4404.
[Sekihara et al., 2001] Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., and
Miyashita, Y . (2001). Reconstructing spatio-temporal activities of neural sources
using an meg vector beamformer technique. Biomedical Engineering, IEEE
Transactions on, 48(73-d beamformer: x, y and z direction):760–771.
[Sekihara et al., 1997] Sekihara, K., Poeppel, D., Marantz, A., Koizumi, H., and
Miyashita, Y . (1997). Noise covariance incorporated meg-music algorithm: a method
for multiple-dipole estimation tolerant of the influence of background brain activity.
Biomedical Engineering, IEEE Transactions on, 44(9):839–847.
114
[Shehzad et al., 2009] Shehzad, Z., Kelly, A. M. C., Reiss, P. T., Gee, D. G., Gotimer,
K., Uddin, L. Q., Lee, S. H., Margulies, D. S., Roy, A. K., Biswal, B. B., Petkova, E.,
Castellanos, F. X., and Milham, M. P. (2009). The Resting Brain: Unconstrained yet
Reliable. Cereb. Cortex, page 256.
[Siapas et al., 2005] Siapas, A. G., Lubenov, E. V ., and Wilson, M. A. (2005). Pre-
frontal phase locking to hippocampal theta oscillations.
[Siegel et al., 2009] Siegel, M., Warden, M. R., and Miller, E. K. (2009). Phase-
dependent neuronal coding of objects in short-term memory. Proceedings of the
National Academy of Sciences.
[Spencer et al., 1992] Spencer, M., Leahy, R., Mosher, J., and Lewis, P. (1992). Adap-
tive filters for monitoring localized brain activity from surface potential time series.
pages 156–161 vol.1.
[Sporns et al., 2007] Sporns, O., Honey, C. J., and Ktter, R. (2007). Identification and
classification of hubs in brain networks. PLoS ONE, 2(10):e1049.
[Squire et al., 2008] Squire, L., Bloom, F., and Spitzer, N. (2008). Fundamental neuro-
science.
[Stam et al., 2006] Stam, C., Jones, B., Manshanden, I., van Cappellen van Walsum, A.,
Montez, T., Verbunt, J., de Munck, J., van Dijk, B., Berendse, H., and Scheltens, P.
(2006). Magnetoencephalographic evaluation of resting-state functional connectivity
in alzheimer’s disease. NeuroImage, 32(3):1335–1344.
[Stam et al., 2009] Stam, C. J., de Haan, W., Daffertshofer, A., Jones, B. F., Manshan-
den, I., van Cappellen van Walsum, A. M., Montez, T., Verbunt, J. P. A., de Munck,
J. C., van Dijk, B. W., Berendse, H. W., and Scheltens, P. (2009). Graph theoretical
analysis of magnetoencephalographic functional connectivity in alzheimer’s disease.
Brain, 132(1):213–224.
[Stefan et al., 2003] Stefan, H., Hummel, C., Scheler, G., Genow, A., Druschky, K.,
Tilz, C., Kaltenhauser, M., Hopfengartner, R., Buchfelder, M., and Romstock, J.
(2003). Magnetic brain source imaging of focal epileptic activity: a synopsis of 455
cases. Brain, 126(11):2396.
[Stephens, 1969] Stephens, M. A. (Mar., 1969). Tests for the von mises distribution.
Biometrika, 56(1):149–160.
[Sternberg, 1966] Sternberg, S. (1966). High-speed scanning in human memory.
Science, 153(736):652–4.
115
[Sutherling et al., 1988] Sutherling, W., Crandall, P., Darcey, T., Becker, D., Levesque,
M., and Barth, D. (1988). The magnetic and electric fields agree with intracranial
localizations of somatosensory cortex. Neurology, 38(11):1705.
[Tikhonov et al., 1977] Tikhonov, A., Arsenin, V ., and John, F. (1977). Solutions of
ill-posed problems.
[van den Heuvel et al., 2008] van den Heuvel, M., Stam, C., Boersma, M., and Hul-
shoff Pol, H. (2008). Small-world and scale-free organization of voxel-based resting-
state functional connectivity in the human brain. NeuroImage, 43(3):528–539.
[Van Veen and Buckley, 1988] Van Veen, B. and Buckley, K. (1988). Beamforming: a
versatile approach to spatial filtering. ASSP Magazine, IEEE, 5(2):4–24.
[Van Veen et al., 1997] Van Veen, B., Van Drongelen, W., Yuchtman, M., and Suzuki,
A. (1997). Localization of brain electrical activity via linearly constrained minimum
variance spatial filtering. Biomedical Engineering, IEEE Transactions on, 44(9):867–
880.
[Varela et al., 2001] Varela, F., Lachaux, J.-P., Rodriguez, E., and Martinerie, J. (2001).
The brainweb: Phase synchronization and large-scale integration. Nature Reviews
Neuroscience Nat Rev Neurosci, 2(4):229–239.
[Watts and Strogatz, 1998] Watts, D. J. and Strogatz, S. H. (1998). Collective dynamics
of /‘small-world/’ networks. Nature, 393(6684):440–442.
[Weissenbacher et al., ] Weissenbacher, A., Kasess, C., Gerstl, F., Lanzenberger, R.,
Moser, E., and Windischberger, C. Correlations and anticorrelations in resting-state
functional connectivity mri: A quantitative comparison of preprocessing strategies.
NeuroImage, In Press, Uncorrected Proof.
[Won et al., 2009] Won, S., Morris, N., Lu, Q., and Elston, R. C. (2009). Choosing
an optimal method to combine ¡i¿p¡/i¿-values. Statistics in Medicine, 28(11):1537–
1553.
[Worsley et al., 2004] Worsley, K., Taylor, J., Tomaiuolo, F., and Lerch, J. (2004). Uni-
fied univariate and multivariate random field theory. NeuroImage, 23(Supplement
1TY - JOUR):S189–S195.
[Wu et al., 2007] Wu, X., Chen, X., Li, Z., Han, S., and Zhang, D. (2007). Binding of
verbal and spatial information in human working memory involves large-scale neural
synchronization at theta frequency. NeuroImage, 35(4):1654–1662.
[Wu et al., 2006] Wu, X., Yao, L., Long, Z.-y., Lu, J., and Li, K.-c. (2006). Functional
connectivity in the resting brain: An analysis based on ica. In Neural Information
Processing, pages 175–182–.
116
[Zaveri et al., 1999] Zaveri, H. P., Williams, W. J., Sackellares, J. C., Beydoun, A.,
Duckrow, R. B., and Spencer, S. S. (1999). Measuring the coherence of intracranial
electroencephalograms. Clinical Neurophysiology, 110(10):1717–1725.
[Zimmerman et al., 1970] Zimmerman, J., Thiene, P., and Harding, J. (1970). Design
and Operation of Stable rf-Biased Superconducting Point-Contact Quantum Devices,
and a Note on the Properties of Perfectly Clean Metal Contacts. Journal of Applied
Physics, 41:1572.
117

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Measuring functional connectivity of the brain

PDF

Statistical signal processing of magnetoencephalography data

PDF

Signal processing methods for brain connectivity

PDF

Causality and consistency in electrophysiological signals

PDF

Functional connectivity analysis and network identification in the human brain

PDF

Brain connectivity in epilepsy

PDF

Diffusion MRI of the human brain: signal modeling and quantitative analysis

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

PDF

On the electrophysiology of multielectrode recordings of the basal ganglia and thalamus to improve DBS therapy for children with secondary dystonia

Asset Metadata

Creator Hui, Hua Brian (author)

Core Title Signal processing methods for interaction analysis of functional brain imaging data

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 03/28/2010

Defense Date 12/21/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag beamformer,cross-talk,functional interaction,intracranial EEG,MEG,OAI-PMH Harvest,phase synchrony

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Leahy, Richard M. (committee chair), Nayak, Krishna (committee member), Singh, Manbir (committee member)

Creator Email brianhuahui@gmail.com,hhui@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2879

Unique identifier UC1416528

Identifier etd-Hui-3497 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-292402 (legacy record id),usctheses-m2879 (legacy record id)

Legacy Identifier etd-Hui-3497.pdf

Dmrecord 292402

Document Type Dissertation

Rights Hui, Hua Brian

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu