Sylvain Baillet, John C. Mosher, and Richard M. Leahy: 14 November 2001 1053-5888/01/$10.00©2001IEEE

Sylvain Baillet, John C.
Mosher,
and Richard M. Leahy
T
he past 15 years have
seen tremendous ad-
vances in our ability to
produce images of hu-
man brain function. Applications
of functional brain imaging extend
from improving our understand-
ing of the basic mechanisms of
cognitive processes to better char-
acterization of pathologies that
impair normal function. Mag-
netoencephalography (MEG) and
electroencephalography (EEG)
(MEG/EEG) localize neural elec-
trical activity using noninvasive
measurements of external electro-
magnetic signals. Among the
available functional imaging tech-
©1997 PHOTODISC, INC. AND 1989-97 TECH POOL STUDIOS, INC.

niques, MEG and EEG uniquely
have temporal resolutions below
100 ms. This temporal precision
allows us to explore the timing of
basic neural processes at the level
of cell assemblies. MEG/EEG
source localization draws on a
wide range of signal processing
techniques including digital filter-
ing, three-dimensional image
analysis, array signal processing,
image modeling and reconstruc-
tion, and, more recently, blind
source separation and phase syn-
chrony estimation. In this article we describe the under- Introduction
lying models currently used in MEG/EEG source Functional brain imaging is a relatively new and
estimation and describe the various signal processing multidisciplinary research field that encompasses tech-
steps required to compute these sources. In particular niques devoted to a better understanding of the human
we describe methods for computing the forward fields brain through noninvasive imaging of the electrophysio-
for known source distributions and parametric and im- logical, hemodynamic, metabolic, and neurochemical
aging-based approaches to the inverse problem. processes that underlie normal and pathological brain
14 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

1053-5888/01/$10.00©2001IEEE
function. These imaging techniques are powerful tools
for studying neural processes in the normal working Functional brain imaging is a
brain. Clinical applications include improved under-
standing and treatment of serious neurological and multidisciplinary research field
neuropsychological disorders such as intractable epilepsy, that encompasses techniques
schizophrenia, depression, and Parkinson’s and Alzhei-
mer’s diseases. devoted to a better understanding
Brain metabolism and neurochemistry can be studied of processes that underlie normal
using radioactively labeled organic molecules, or probes,
that are involved in processes of interest such as glucose and pathological brain function.
metabolism or dopamine synthesis [1]. Images of dy-
namic changes in the spatial distribution of these probes,
the inherent ambiguity of the underlying static electro-
as they are transported and chemically modified within
magnetic inverse problem. Only by placing restrictive
the brain, can be formed using positron emission tomog-
raphy (PET). These images have spatial resolutions as models on the sources of MEG and EEG signals can we
high as 2 mm; however, temporal resolution is limited by achieve resolutions similar to those of fMRI and PET.
the dynamics of the processes being studied, and by pho- Reviews of the application of MEG and EEG to neurol-
ton-counting noise, to several minutes. For more direct ogy and neuropsychology can be found elsewhere [5]-[9].
studies of neural activity, one can investigate local We recommend [10] for a thorough review of MEG the-
hemodynamic changes. As neurons become active, they ory and instrumentation. This article provides a brief intro-
induce very localized changes in blood flow and oxygena- duction to the topic with an overview of the associated
tion levels that can be imaged as a correlate of neural activ- inverse problem from a signal processing perspective. In
ity. Hemodynamic changes can be detected using PET the next two sections we describe the sources of MEG and
[1], functional magnetic resonance imaging (fMRI) [2], EEG signals and how they are measured. Neural sources
and transcranial optical imaging [3] methods. Of these, and head models are then described, followed by the vari-
fMRI is currently the most widely used and can be readily ous approaches to the inverse problem in which the prop-
performed using a standard 1.5T clinical MRI magnet, al- erties of the neural current generators are estimated from
though an increasing fraction of studies are now per- the data. We conclude with a discussion of recent develop-
formed on higher field (3-4.5T) machines for improved ments and our perspective on emerging signal processing
SNR and resolution. fMRI studies are capable of produc- issues for EEG and MEG data analysis.
ing spatial resolutions as high as 1-3 mm; however, tem-
poral resolution is limited by the relatively slow
hemodynamic response, when compared to electrical Sources of EEG and MEG:
neural activity, to approximately 1 s. In addition to lim- Electrophysiological Basis
ited temporal resolution, interpretation of fMRI data is
hampered by the rather complex relationship between the MEG and EEG are two techniques that exploit what
blood oxygenation level dependent (BOLD) changes Galvani, at the end of the 18th century, called “animal
that are detected by fMRI and the underlying neural ac- electricity,” today better known as electrophysiology
tivity. Regions of BOLD changes in fMRI images do not [11]. Despite the apparent simplicity in the structure of
necessarily correspond one-to-one with regions of electri- the neural cell, the biophysics of neural current flow relies
cal neural activity. on complex models of ionic current generation and con-
MEG and EEG are two complementary techniques duction [12]. Roughly, when a neuron is excited by
that measure, respectively, the magnetic induction out- other—and possibly remotely located—neurons via an af-
side the head and the scalp electric potentials produced by ferent volley of action potentials, excitatory postsynaptic
electrical activity in neural cell assemblies. They directly potentials (EPSPs) are generated at its apical dendritic
measure electrical brain activity and offer the potential for tree. The apical dendritic membrane becomes transiently
superior temporal resolution when compared to PET or depolarized and consequently extracellularly electro-
fMRI, allowing studies of the dynamics of neural net- negative with respect to the cell soma and the basal den-
works or cell assemblies that occur at typical time scales drites. This potential difference causes a current to flow
on the order of tens of milliseconds [4]. Sampling of elec- through the volume conductor from the nonexcited
tromagnetic brain signals at millisecond intervals is membrane of the soma and basal dendrites to the apical
readily achieved and is limited only by the multichannel dendritic tree sustaining the EPSPs [13].
analog-to-digital conversion rate of the measurements. Some of the current takes the shortest route between
Unfortunately, the spatial resolving power of MEG and the source and the sink by traveling within the dendritic
EEG does not, in general, match that of PET and fMRI. trunk (see Fig. 1). Conservation of electric charges im-
Resolution is limited both by the relatively small number poses that the current loop be closed with extracellular
of spatial measurements—a few hundred in MEG or EEG currents flowing even through the most distant part of
versus tens of thousands or more in PET or fMRI—and the volume conductor. Intracellular currents are com-
NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 15

monly called primary currents, while extracellular cur- sume the cortex is about 4 mm thick, then a small patch 5
rents are known as secondary, return, or volume currents. mm × 5 mm would yield a net current of 10 nA-m, consis-
Both primary and secondary currents contribute to tent with empirical observations and invasive studies.
magnetic fields outside the head and to electric scalp po- At a larger scale, distributed networks of collaborating
tentials, but spatially structured arrangements of cells are and synchronously activated cortical areas are major con-
of crucial importance to the superposition of neural cur- tributors to MEG and EEG signals. These cortical areas are
rents such that they produce measurable fields. compatible with neuroscientific theories that model basic
Macrocolumns of tens of thousands of synchronously ac- cognitive processes in terms of dynamically interacting cell
tivated large pyramidal cortical neurons are thus believed assemblies [15]. Although cortical macrocolumns are as-
to be the main MEG and EEG generators because of the sumed to be the main contributors to MEG and EEG sig-
coherent distribution of their large dendritic trunks lo- nals [4], some authors have reported scalp recordings of
cally oriented in parallel, and pointing perpendicularly to deeper cortical structures including the hippocampus [16],
the cortical surface [14]. The currents associated with the cerebellum [17], and thalamus [18], [19].
EPSPs generated among their dendrites are believed to be
at the source of most of the signals detected in MEG and
EEG because they typically last longer than the rapidly Measuring EEG and MEG signals
firing action potentials traveling along the axons of ex- Electroencephalography
cited neurons [4]. Indeed, calculations such as those EEG was born in 1924 when the German physician Hans
shown in [10] suggest each synapse along a dendrite may Berger first measured traces of brain electrical activity in
contribute as little as a 20 fA-m current source, probably humans. Although today’s electronics and software for
too small to measure in MEG/EEG. Empirical observa- EEG analysis benefit from the most recent technological
tions instead suggest we are seeing sources on the order of developments, the basic principle remains unchanged
10 nA-m, and hence the cumulative summation of mil- from Berger’s time. EEG consists of measurements of a
lions of synaptic junctions in a relatively small region. set of electric potential differences between pairs of scalp
Nominal calculations of neuronal density and cortical electrodes. The sensors may be either directly glued to the
thickness suggest that the cortex has a macrocellular cur- skin ( for prolonged clinical observation) at selected loca-
rent density on the order of 100 nA/mm2 [10]. If we as- tions directly above cortical regions of interest or fitted in
Pyramidal Cell Assembly

Excitatory Post-Synaptic Potentials
Scalp Pyramidal Cell Assembly

Skull Neural Activation
CSF
Primary
Current
Secondary
Currents
Cortex
CNRS UPR640 - USC - LANL
▲ 1. Networks of cortical neural cell assemblies are the main generators of MEG/EEG signals. Left: Excitatory postsynaptic potentials
(EPSPs) are generated at the apical dendritic tree of a cortical pyramidal cell and trigger the generation of a current that flows
through the volume conductor from the non-excited membrane of the soma and basal dendrites to the apical dendritic tree sustain-
ing the EPSPs. Some of the current takes the shortest route between the source and the sink by travelling within the dendritic trunk
(primary current in blue), while conservation of electric charges imposes that the current loop be closed with extracellular currents
flowing even through the most distant part of the volume conductor (secondary currents in red). Center: Large cortical pyramidal
nerve cells are organized in macro-assemblies with their dendrites normally oriented to the local cortical surface. This spatial ar-
rangement and the simultaneous activation of a large population of these cells contribute to the spatio-temporal superposition of the
elemental activity of every cell, resulting in a current flow that generates detectable EEG and MEG signals. Right: Functional networks
made of these cortical cell assemblies and distributed at possibly mutliple brain locations are thus the putative main generators of
MEG and EEG signals.

an elastic cap for rapid attachment with near uniform cov- sensor at MIT. SQUIDs can be used to detect and quantify
erage of the entire scalp. Research protocols can use up to minute changes in the magnetic flux through magnetome-
256 electrodes. ter coils in a superconducting environment. D.S. Cohen,
EEG has had tremendous success as a clinical tool, es- also at MIT, made the first MEG recording a few years
pecially in studying epilepsy, where seizures are charac- later [22].
terized by highly abnormal electrical behavior in neurons Recent developments include whole-head sensor arrays
in epileptogenic regions. In many clinical and research ap- for the monitoring of brain magnetic fields at typically 100
plications, EEG data are analyzed using pattern analysis to 300 locations. Noise is a major concern for MEG. Instru-
methods to associate characteristic differences in the data mental noise is minimized by the use of superconducting
with differences in patient populations or experimental materials and immersing the sensing setup in a Dewar
paradigm. The methods described here for estimating the cooled with liquid helium. High-frequency perturbations
location, extent and dynamic behavior of the actual cur- such as radiofrequency waves are easily attenuated in
rent sources in the brain are currently less widely used in shielded rooms made of successive layers of mu-metal, cop-
clinical EEG. per, and aluminum (Fig. 2). Low frequency artifacts created
Though dramatic changes in the EEG, such as interictal by cars, elevators, and other moving objects near the MEG
spikes occurring between epileptic seizures, may be readily system are attenuated by the use of gradiometers as sensing
visible in raw measurements, event-related signals associ- units. A gradiometer is the hardware combination of multi-
ated with, for example, presentation of a specific sensory ple magnetometers to physically mimic the computation of
stimulus or cognitive challenge, are often lost in back- the spatial gradient of the magnetic induction in the vicinity
ground brain activity. Dawson demonstrated in 1937 that of the head. Noise sources distant from the gradiometer
by adding stimulus-locked EEG traces recorded during produce magnetic fields with small spatial gradients and
several instances of the same stimulus, one could reveal hence are effectively attenuated using this mechanism.
spatio-temporal components of the EEG signal related As with EEG, MEG has potentially important applica-
with that stimulus, and background noise would be mini- tions in clinical studies where disease or treatments affect
mized. This method of “stimulus-locked” averaging of spontaneous or event-related neural activity. Again, stim-
event-related EEG is now a standard technique for noise ulus-locked averaging is usually required to reduce back-
reduction in event-related studies. Liquid Helium
Averaging, however, relies on the strong hypoth-
esis that the brain is in a stationary state during the Shielded Room
experiment with insignificant adaptation or habitu-
ar
ation to experimental conditions during repeated
Dew
exposure to a stimulus or task. In general, this Sensors
stationarity does not hold true, especially as the
number of trials increases, which has motivated new MEG
research approaches that study the inter-trial varia- Gantry
tions by greatly reducing the number of trials in Time-Evolving Scalp Magnetic Field Topographies
each average, or by analyzing the raw unaveraged −9.76 −0.16 9.44
EEG data [20], [21].
Magnetoencephalography
28.65 38.26 47.87
Typical EEG scalp voltages are on the order of tens of
microvolts and thus readily measured using relatively
low-cost scalp electrodes and high-impedance
high-gain amplifiers. In contrast, characteristic mag-
netic induction produced by neural currents is ex- ▲ 2. MEG instrumentation and typical signals. Typical scalp magnetic
traordinarily weak, on the order of several tens of fields are on the order of a 10 billionth of the earth’s magnetic field.
femtoTeslas, thus necessitating sophisticated sensing MEG fields are measured inside a magnetically shielded room for pro-
technology. In contrast to EEG, MEG was devel- tection against higher-frequency electromagnetic perturbations (left).
oped in physics laboratories and especially in MEG sensors use low-temperature electronics cooled by liquid helium
(upper right) stored in a Dewar (left and upper right). Scalp magnetic
low-temperature and superconductivity research
fields are then recorded typically every millisecond. The resulting data
groups. In the late 1960s, J.E. Zimmerman co-in-
can be visualized as time-evolving scalp magnetic field topographies
vented the SQUID (Superconducting QUantum In- (lower right). These plots display the time series of the recorded mag-
terference Device)—a supremely sensitive amplifier netic fields interpolated between sensor locations on the subject’s scalp
that has since found applications ranging from air- surface. This MEG recording was acquired as the subject moved his fin-
borne submarine sensing to the detection of gravita- ger at time 0 (time relative to movement (t=0) is indicated in ms above
tional waves—and conducted the first human every topography). Data indicate early motor preparation prior to the
magnetocardiogram experiment using a SQUID movement onset before peaking at about 20 ms after movement onset.

neural activity and a volume (or passive) current flow
The excellent time resolution of J V (r′ ) that results from the effect of the electric field in
the volume on extracellular charge carriers:
MEG and EEG gives us a unique
window on the dynamics of J(r′ ) = J p (r′ ) + J V (r′ ) = J p (r′ ) + σ(r′ ) E(r′ )
= J p (r′ ) − σ(r′ )∇V(r′ )
human brain functions.
where σ(r ′) is the conductivity profile of the head tissues,
ground noise to a point where event-related signals can be
which we will assume, for simplicity, to be isotropic, and,
seen in the data. The major attraction of MEG, as com-
from the quasi-static assumption, the electric field E(r′ ) is
pared to EEG, is that while EEG is extremely sensitive to
the effects of the secondary or volume currents, MEG is the negative gradient of the electric potential, V(r′ ).
more sensitive to the primary current sources in which we If we assume that the head consists of a set of contiguous
are typically more interested—this will become clearer in regions each of constant isotropic conductivityσ i , i = 1,...,3,
our review of the forward model below. Contributors to representing the brain, skull and scalp for instance, we can
the initial developments of MEG put great emphasis over rewrite the Biot-Savart law above as a sum of contributions
the past two decades on the use of inverse methods to from the primary and volume currents [10]:
characterize the true sources of MEG signals within the
µ0 r − r′
brain. More recently, EEG and MEG have come to be B(r) = B0 (r) +
4π
∑ (σ i
− σ j )∫ V (r′ )
S ij
r − r′
3
× dS ′ij ,
viewed as complementary rather than competing modali- ij
(2)
ties and most MEG facilities are equipped for simulta-
neous acquisition of both EEG and MEG data. As we where B0 (r)is the magnetic field due to the primary current
shall see, inverse methods for the two are very closely re- only. The second term is the volume current contribution to
lated and can even be combined and optimized for joint the magnetic field formed as a sum of surface integrals over
source localization [23]. the brain-skull, skull-scalp, and scalp-air boundaries.
This general equation states that the magnetic field can
The Physics of MEG and EEG: be calculated if we know the primary current distribution
and the potential V(r′ ) on all surfaces. We can create a
Source and Head Models similar equation for the potential itself, although the deri-
Given a set of MEG or EEG signals from an array of exter- vation is somewhat tedious [10], [24], yielding
nal sensors, the inverse problem involves estimation of
(σ i + σ j )V (r) = 2σ 0 V 0 (r)
the properties of the current sources within the brain that
produced these signals. Before we can make such an esti- 1 r − r′
mate, we must first understand and solve the forward − ∑ (σ − σ j )∫S ij V (r′ )
2 π ij i 3
⋅ dS ij′
r − r′ (3)
problem, in which we compute the scalp potentials and
external fields for a specific set of neural current sources.
for the potential on surface S ij where V 0 (r) is the potential
at r due to the primary current distribution.
Quasi-Static Approximation These two equations therefore represent the integral
of Maxwell Equations solutions to the forward problem. If we specify a primary
The useful frequency spectrum for electrophysiological current distribution J P (r′ ), we can calculate a primary
signals in MEG and EEG is typically below 1 kHz, and potential and a primary magnetic field,
most studies deal with frequencies between 0.1 and 100
Hz. Consequently, the physics of MEG and EEG can be 1 r − r′
described by the quasi-static approximation of Maxwell’s
V 0 (r) =
4 πσ 0 ∫ J P (r′ ) ⋅
r − r′
3
dr ′ ,
equations. The quasi-static current flow J(r′ ) at location
r′ is therefore divergence free and can be related rather µ0 r − r′
4π ∫
B0 (r) = J P (r′ ) × 3
dr ′ ,
simply to the magnetic field B(r)at location r through the r − r′ (4)
well-known Biot-Savart law
µ0 r − r′ the primary potential V 0 (r) is then used to solve (3) for
4π ∫
B(r) = J(r′ ) × 3
dv′ the potentials on all surfaces, and therefore solves the for-
r − r′ (1) ward problem for EEG. These surface potentials V(r) and
the primary magnetic field B0 (r) are then used to solve
where µ 0 is the permitivity of free space. We can partition (2) for the external magnetic fields. Unfortunately, the
the total current density in the head volume into two cur- solution of (3) has analytic solutions only for special
rent flows of distinct physiological significance: a primary shapes and must otherwise be solved numerically. We re-
(or driving) current flow J P (r′ ) related to the original turn to specific solutions of the forward problem below

but first we discuss the types of models used to describe Consider the special case of a current dipole of mo-
the primary current distributions. ment q located at rq in a multishell spherical head, and an
MEG system in which we measure only the radial compo-
nent of the magnetic field, i.e., the coil surface of the mag-
Source Models: Dipoles and Multipoles netometer is oriented orthogonally to a radial line from
Let us assume a small patch of activated cortex is centered the center of the sphere through the center of the coil. It is
at location rq and that the observation point r is some dis- relatively straightforward to show that the contributions
tance away from this patch. The primary current distribu- of the volume currents vanish in this case, and we are left
tion in this case can be well approximated by an with only the primary term B0 (r). Taking the radial com-
equivalent current dipole represented as a point source ponent of this field for the current dipole reduces to the
J P (r′ ) = qδ(r′− rq ), where δ(r) is the Dirac delta function, remarkably simple form:
with moment q ≡ ∫ J P (r′ )dr′. The current dipole is a
straightforward extension of the better-known model of r r µ r × rq
Br (r) ≡ ⋅ B(r) = ⋅ B0 (r) = 0 ⋅ q.
the paired-charges dipole in electrostatics. It is important r r 4π r r − r 3
to note that brain activity does not actually consist of dis-

q
(5)
crete sets of physical current dipoles, but rather that the
Note that this magnetic field measurement is linear in
dipole is a convenient representation for coherent activa-
the dipole moment q but highly nonlinear with respect to
tion of a large number of pyramidal cells, possibly extend-
its location rq . Although we do not reproduce the results
ing over a few square centimeters of gray matter. here, the magnetic fields for arbitrary sensor orientation
The current dipole model is the workhorse of and the scalp potentials for the spherical head models can
MEG/EEG processing since a primary current source of both be written in a form similar to (5) as the inner product
arbitrary extent can always be broken down into small re- of a linear dipole moment with a term that is nonlinear in
gions, each region represented by an equivalent current location [28]. While it may not be immediately obvious,
dipole. This is the basis of the imaging methods described this property also applies to numerical solution of (2) and
later on. However, an identifiability problem can arise (3), i.e., to arbitrary geometries of the volume conductor,
when too many small regions and their dipoles are re- and the measured fields remain linear in the dipole mo-
quired to represent a single large region of coherent acti- ment and nonlinear in the dipole location [10], [28].
vation. These sources may be more simply represented by From (5) we can also see that due to the triple scalar
a multipolar model. The multipolar models can be gener- product, Br (r) is zero everywhere outside the head if q
ated by performing a Taylor series expansion of the points towards the radial direction rq . A more general re-
3
Green’s function G(r , r′ ) = r − r′ / r − r′ about the cen-
sult is that radially oriented dipoles do not produce any
troid of the source. Successive terms in the expansion give
external magnetic field outside a spherically symmetric
rise to the multipolar components: dipole, quadrupole,
volume conductor, regardless of the sensor orientation
octupole, and so on. The first multipolar definitions for
[31]. Importantly, this is not the case for EEG, which is
electrophysiological signals were established in
sensitive to radial sources, constituting one of the major
magnetocardiography [25]. In MEG, the contributions
differences between MEG and EEG data.
to the magnetic field from octupolar and higher order
terms drop off rapidly with distance, so that restricting
sources to dipolar and quadrupolar fields is probably suf- Realistic Head Models
ficient to represent most plausible cortical sources [26]. We have described how the forward models have
An alternative approach to multipolar models of brain closed-form solution for heads with conductivity profiles
sources can be found in [27]. that can be modeled as a set of nested concentric homoge-
neous and isotropic spheres. In reality, of course, we do
Head Models not have heads like this—our heads are anisotropic,
inhomogeneous, and not spherical. Rather surprisingly,
Spherical Head Models
the spherical models work reasonably well, particularly
Computation of the scalp potentials and induced mag- for MEG measurements, which are less sensitive than
netic fields requires solution of the forward equations (3) EEG to the effects of volume currents, which, in turn, are
and (2), respectively, for a particular source model. When affected more than primary currents by deviations from
the surface integrals are computed over realistic head the idealized model [32].
shapes, these equations must be solved numerically. Ana- More accurate solutions to the forward problem use
lytic solutions exist, however, for simplified geometries, anatomical information obtained from high-resolution
such as when the head is assumed to consist of a set of volumetric brain images obtained with MR or X-ray
nested concentric homogeneous spherical shells repre- computed tomography (CT) imaging. Since MR scans
senting brain, skull, and scalp [28]-[30]. These models are now routinely performed as part of most MEG proto-
are routinely used in most clinical and research applica- cols, this data is readily available. To solve (2) and (3) we
tions to MEG/EEG source localization. must extract surface boundaries for brain, skull, and scalp

geneous, but it is unlikely that the EIT approach will be
One of the most exciting current able to produce high-resolution images of spatially vary-
ing anisotropic conductivity.
challenges in functional brain A second approach to imaging conductivity is to use
mapping is the question of how magnetic resonance. One technique uses the shielding ef-
fects of induced eddy currents on spin precession and
to best integrate data from could in principle help determine the conductivity profile
different modalities. at any frequency [39]. The second technique uses diffu-
sion-tensor imaging with MRI (DT-MRI), which probes
from these images. Many automated and semiautomated the microscopic diffusion properties of water molecules
methods exist for surface extraction from MR images, within the tissues of the brain [40]. The diffusion values
e.g., [33], [34]. The surfaces can then be included in a can then be tentatively related to the conductivity of these
boundary element method (BEM) calculation of the for- tissues [41]. Both of these MR-based techniques are still
ward fields. under investigation, but given the poor signal-to-noise
While this is an improvement on the spherical model, ratio (SNR) of the MR in bone regions, which is of criti-
the BEM methods still assume homogeneity and isotropy cal importance for the forward EEG problem, the poten-
within each region of the head. This ignores, for example, tial for fully three-dimensional impedance tomography
anisotropy in white matter tracts in the brain in which with MR remains speculative.
conduction is preferentially along the axonal fibers com-
pared to the transverse direction. Similarly, the sinuses Algebraic Formulation
and diploic spaces in the skull make it very inhomo- With the introduction of the source and head models for
geneous, a factor that is typically ignored in BEM calcula- solution of the forward problem, we can now provide a
tions. The finite element method (FEM) can deal with all few key definitions and linear algebraic models that will
of these factors and therefore represents a very powerful clarify the different approaches taken in the inverse meth-
approach to solving the forward problem. Typically BEM ods described in the next section. As we saw above, the
and FEM calculations are very time consuming and use of magnetic field and scalp potential measurements are linear
these methods may appear impractical when incorpo- with respect to the dipole moment q and nonlinear with re-
rated as part of an iterative inverse solution. In fact, spect to the location r q . For reasons of exposition it is con-
through use of fast numerical methods, precalculation, venient to separate the dipole magnitude q ≡ q from its
and look-up tables and interpolation of precalculated orientationΘ = q / q which we represent in spherical coor-
fields, both FEM and BEM can be made quite practical dinates by Θ = {θ,ϕ}. Let m(r) denote either the scalp elec-
for applications in MEG and EEG [35]. tric potential or magnetic field generated by a dipole:
One problem remains: these methods need to know
the conductivity of the head. Most of the head models m(r) = a(r , rq ,Θ)q, (6)
used in the bioelectromagnetism community consider
typical values for the conductivity of the brain, skull, and where a(r , rq ,Θ) is formed as the solution to either the
skin. The skull is typically assumed to be 40 to 90 times magnetic or electric forward problem for a dipole with
more resistive than the brain and scalp, which are as- unit amplitude and orientation Θ.
sumed to have similar conductive properties. These val- For the simultaneous activation of multiple dipoles
ues are measured in vitro from postmortem tissue, where located at rqi , and by linear superposition, we can sim-
conductivity can be significantly altered compared to in ply sum the individual contributions to obtain
vivo values [36]. Consequently, recent research efforts m(r) = ∑ i a(r , rqi ,Θ i )q i . For the simultaneous EEG or
have focused on in vivo measures. MEG measurements made at N sensors we obtain
Electrical impedance tomography (EIT) proceeds by
injecting a small current (1-10 microA) between pairs of  m(r1 )   a(r1 , rq1 ,Θ1 ⋅⋅⋅ a(r1 , rqp ,Θ p )   q1 
   
EEG electrodes and measuring the resulting potentials at m=   =   
 
all electrodes. Given a model for the head geometry, EIT  m(rN )  a(rN , rq1 ,Θ1 ) ⋅⋅⋅ a(rN , rqp ,Θ p ) q p 
 
solves an inverse problem by minimizing the error be-
tween the measured potentials of the rest of the EEG = A({rqi ,Θ i })S T ,
(7)
leads and the model-based computed potentials, with re-
spect to the parameters of the conductivity profile. Re- where A(rqi ,Θ i ) is the gain matrix relating the set of p di-
cent simulation results with three or four-shell spherical poles to the set of N discrete locations (now implicitly a
head models have demonstrated the feasibility of this ap- function of the set of N sensor locations), m is a generic
proach though the associated inverse problem is also fun- set of N MEG or EEG measurements, and the matrix S is
damentally ill-posed [37], [38]. These methods are a generalized matrix of source amplitudes, defined below.
readily extendible to realistic surface models as used in Each column of A relates a dipole to the array of sensor
BEM calculations in which each region is assumed homo- measurements and is called the forward field, gain vector,

or scalp topography, of the current dipole source sampled between the data and the fields computed from the
by the N discrete locations of the sensors. This model can estimated sources using a forward model. Each dipole
be readily extended to include a time component t, when represented in the matrix A({rqi ,Θ i }) comprises three
considering time evolving activities at every dipole loca- nonlinear location parameters r qi , a set of two nonlinear
tion. For p sources and T discrete time samples, the orientation parameters Θ i = (θ i , φ i ), and the T linear di-
spatio-temporal model can therefore be represented as pole amplitude time series parameters in the vector s i .
For p dipoles, we define the measure of fit in the
 m(r1 ,1) ⋅⋅⋅ m(r1 ,T )   s 1T  least-square (LS) sense as the square of the Frobenius norm
 
M=  = A r ,Θ
({ i i })  
  2
 m(rS ,1) ⋅⋅⋅ m(rS ,T )  s Tp  J LS ({rqi ,Θ i },S) = M − A({rqi ,Θ i })S T .
  F (9)
T
= A({ri ,Θ i })S . (8) A brute force approach is to use a nonlinear search pro-
gram to minimize J LS over all of parameters ({rqi ,Θ i },S)
The corresponding time series for each dipole are the col- simultaneously; however, the following simple optimal
umns of the time series matrix S, where S T indicates the modification greatly reduces the computational burden.
matrix is transposed. Because the orientation of the dipole For any selection of {rqi ,Θ i }, the matrix S that will mini-
is not a function of time, this type of model is often referred mize J LS is
to as a “fixed” dipole model. Alternative models that allow
these dipoles to “rotate” as a function of time were intro- S T = A + M, (10)
duced in [42] and are extensively reviewed in [43].
where A + is the pseudoinverse of A = A({rqi ,Θ i }). If A
is of full column rank, then the pseudoinverse may be ex-
Imaging Electrical Activity in the Brain: plicitly written as A + = ( A T A) −1 A T [44], [45]. We can
The Inverse Problem then solve (9) in {r qi ,Θ i } by minimizing the adjusted
cost function:
Parametric and imaging methods are the two general ap-
proaches to estimation of EEG and MEG sources. The
parametric methods typically assume that the sources can
J LS ({r
qi
,Θ i }) = M − A( A + M )
2
be represented by a few equivalent current dipoles of un- = ( I − AA + ) M

2
= PA⊥ M
2
,
known location and moment to be estimated with a non- F F (11)
linear numerical method. The imaging methods are based
on the assumption that primary sources are intracellular where PA⊥ is the orthogonal projection matrix onto the
currents in the dendritic trunks of the cortical pyramidal left null space of A. Thus, the LS problem can be opti-
neurons, which are aligned normally to the cortical sur- mally solved in the limited set of nonlinear parameters
face. Thus a current dipole is assigned to each of many {rqi ,Θ i } with an iterative minimization procedure. The
tens of thousands of tessellation elements on the cortical linear parameters in S are then optimally estimated from
surface with the dipole orientation constrained to equal (10) [43], [45]. Minimization methods range from
the local surface normal. The inverse problem in this case Levenberg-Marquardt and Nelder-Meade downhill
is linear, since the only unknowns are the amplitudes of simplex searches to global optimization schemes using
the dipoles in each tessellation element. Given that the multistart methods, genetic algorithms and simulated
number of sensors is on the order of 100 and the number annealing [46].
of unknowns is on the order of 10,000, the problem is se- This least-squares model can either be applied to a sin-
verely underdetermined, and regularization methods are gle snapshot or a block of time samples. When applied se-
required to restrict the range of allowable solutions. In quentially to a set of individual time slices, the result is
this section we will describe parametric and imaging ap- called a “moving dipole” model, since the location is not
proaches, contrasting the underlying assumptions and constrained [47]. Alternatively, by using the entire block
the limitations inherent in each. of data in the least-squares fit, the dipole locations can be
fixed over the entire interval [42]. The fixed and moving
dipole models have each proven useful in both EEG and
Parametric Modeling MEG and remain the most widely used approach to pro-
Least-Squares Source Estimation cessing experimental and clinical data.
In the presence of measurement errors, the forward A key problem with the LS method is that the number
model may be represented as M = A({rqi ,Θ i })S T + ε, of sources to be used must be decided a priori. Estimates
where ε is a spatio-temporal noise matrix. Our goal is to can be obtained by looking at the effective rank of the data
determine the set {rqi ,Θ i } and the time series S that best using SVD or through information-theoretic criteria, but
describe our data. The earliest and most straightforward in practice expert data analysts often run several model or-
strategy is to fix the number of sources p and use a nonlin- ders and select results based on physiological plausibility.
ear estimation algorithm to minimize the squared error Caution is obviously required since a sufficiently large

number of sources can be made to fit any data set, regard- min
T
tr{C y } subject to W T A(rq ) = I , (13)
W
less of its quality. Furthermore, as the number of sources
increases, the nonconvexity of the cost function results in wher e C y = E[ yy T ] = W T C m W a nd C m = E[mm T ].
increased chances of trapping in undesirable local min- Solving (13) using the method of Lagrange multipliers
ima. This latter problem can be approached using sto- yields the solution [49]:
chastic or multistart search strategies [46], [48].
[ ]
−1
The alternatives described below avoid the noncon- W = A(rq ) T C m−1 A(rq ) A(rq ) T C m−1 .
(14)
vexity issue by scanning a region of interest that can range
from a single location to the whole brain volume for pos-
Applying this filter to each of the snapshot vectors
sible sources. An estimator of the contribution of each
m(t ), t = 1,...,T , in the data matrix M produces an esti-
putative source location to the data can be derived either
mate of the dipole moment of the source at rq [51]-[53].
via spatial filtering techniques or signal classification indi-
ces. An attractive feature of these methods is that they do By simply changing the location rq in (13), we can pro-
not require any prior estimate of the number of underly- duce an estimate of the neural activity at any location.
ing sources. Unfortunately, the transient and often correlated na-
ture of neural activation in different parts of the brain will
often limit performance of the LCMV as correlations be-
Beamforming Approaches tween different sources will result in partial signal cancel-
A beamformer performs spatial filtering on data from a lation. However, simulation results [51], [54] and recent
sensor array to discriminate between signals arriving evaluations on real data [55] seem to indicate LCMV-
from a location of interest and those originating else- based beamforming methods are robust to moderate lev-
where. Beamforming originated in radar and sonar signal els of source/interference correlation. Similarly, model-
processing but has since found applications in diverse ing errors in the constraint matrix A(rq ) or imprecise
fields ranging from astronomy to biomedical signal pro- dipole locations can result in signal attenuation or even
cessing [49], [50]. cancellation. More elaborate constraints may be designed
Let us consider a beamformer that monitors signals by using eigenvectors that span a desired region to be ei-
from a dipole at location rq , while blocking contributions ther monitored (gain=1) or nulled (gain=0) [49], but as
from all other brain locations. If we do not know the ori- the number of constraints increases, the degrees of free-
entation of the dipole, we need a vector beamformer con- dom are reduced and the beamformer becomes less adap-
sisting of three spatial filters, one for each of the Cartesian
tive to other unknown sources.
axes, which we denote as the set {Θ1 ,Θ 2 ,Θ 3 }. The out-
In the absence of signal, the LCMV beamformer will
put of the beamformer is computed as the three element
produce output simply due to noise. Because of the vari-
vector y(t ) formed as the product of a 3 × N spatial filter-
able sensitivity of EEG and MEG as a function of source
ing matrix W T with m(t ), the signal at the array at time t,
i.e., y(t ) = W T m(t ). location, the noise gain of the filter will vary as a function
of location rq in the constraint A(rq ). A strategy to ac-
The spatial filter would ideally be defined to pass signals
within a small distance δ of the location of interest r q with a count for this effect when using the LCMV beamformer
gain of unity while nulling signals from elsewhere [51]. in a scanning mode is to compute the ratio of the output
Thus the spatial filter should obey the following constraints: variance of the beamformer to that which would have
been obtained in the presence of noise only. It is straight-
 I forward to show that this ratio is given by
r − rq ≤ δ: passband constraint
W T A(r) = 
 0 r − rq > δ: stopband constraint,
(12)
var(r ) =
tr{[ A(r ) C
q
T −1
m ]
A(rq ) }
−1
tr{[ A(r ) C A(r )] }

q −1
T −1
where A(r) = [a(r ,Θ1 ), a(r ,Θ 2 ), a(r ,Θ 3 )] is the N × 3 for- q ε q
(15)
ward matrix for three orthogonal dipoles at location r.
There are insufficient degrees of freedom to enforce a strong
where C ε is an estimate of the noise-only covariance [53].
stop-band constraint over the entire brain volume, so that a
This neural activity index can be extended to statistical
fixed spatial filter is impractical for this application.
parametric mapping (SPM) as in the synthetic aperture
Linearly constrained minimum variance (LCMV)
beamforming provides an adaptive alternative in which magnetometry (SAM) technique [53]; the recent para-
the limited degrees of freedom are used to place nulls in metric mapping method in [56] uses a similar idea to this,
the response at positions corresponding to interfering except that the linear operator applied to the data is a min-
sources, i.e., neural sources at locations other than rq . imum-norm imaging, rather than spatial filtering, matrix.
This nulling is achieved by simply minimizing the output In low noise situations, the signal covariance can be
power of the beamformer subject to a unity gain con- ill-conditioned, and therefore the inverse may be regular-
straint at the desired location rq . The LCMV problem can ized by replacing C −1
m with [C m + λI ]
−1
where λ is a small
be written as positive constant in (14) and (15) [53].

From Classical to RAP-MUSIC clude the use of prewhitening to account for spatial corre-
The multiple signal classification approach (MUSIC) was lations in background brain activity [60] and use of
developed in the array signal processing community [57] time-frequency methods to better select the signal
before being adapted to MEG/EEG source localization subspace of interest [61].
[43]. We will restrict our brief description of the MUSIC One distinct advantage of MUSIC over LCMV meth-
approach here to dipole sources with fixed orientation, al- ods is the relaxation of the requirement of orthogonality
though it can be extended to rotating dipoles and between distinct sources. MUSIC requires the weaker as-
multipoles. As before, let M = A({rqi ,Θ i })S T + ε be an sumption that different sources have linearly independent
N × T spatio-temporal matrix containing the data set un- time series. In noiseless data, partially correlated sources
der consideration for analysis, and let the data be a mix- will still result in a cost function equal to zero at each true
ture of p sources. Let M = UΣV T be the singular value dipole location. In the presence of noise, MUSIC will fail
decomposition (SVD) of M [44]. In the absence of noise, when two sources are strongly or perfectly correlated.
the set of left singular vectors is an orthonormal basis for This problem can be corrected by adjusting the concept of
the subspace spanned by the data. Provided that N > p, single dipole models to specifically allow sets of synchro-
the SNR is sufficiently large, and noise is i.i.d. at the sen- nous sources [58].
sors, one can define a basis for the signal and noise
subspaces from the column vectors of U. The signal Imaging Approaches
subspace is spanned by the p first left singular vectors in
Cortically Distributed Source Models
U, denoted U S , while the noise subspace is spanned by
the remaining left singular vectors. The best rank p ap- Imaging approaches to the MEG/EEG inverse problem
proximation of M is given by M S = (U S U ST ) M and consist of methods for estimating the amplitudes of a dense
PS⊥ = I − (U S U ST ) is the orthogonal projector onto the set of dipoles distributed at fixed locations within the head
noise subspace. volume. In this case, since the locations are fixed, only the
linear parameters need to be estimated, and the inverse
We can define the MUSIC cost function as:
problem reduces to a linear one with strong similarities to
2 those encountered in image restoration and reconstruc-
PS⊥ a(r ,Θ)
J (r ,Θ) = 2
, tion, i.e., the imaging problem involves solution of the lin-
a(r ,Θ)
2
ear system M = AS T for the dipole amplitudes, S.
2 (16) The most basic approach consists of distributing di-
poles over a predefined volumetric grid similar to the
which is zero when a(r ,Θ) corresponds to one of the true
ones used in the scanning approaches. However, since
source locations and orientations, r = rqi and Θ = Θ i ,
primary sources are widely believed to be restricted to the
i = 1,..., p [43]. As in the beamforming approaches, an ad-
cortex, the image can be plausibly constrained to sources
vantage over least-squares is that each source can be
lying on the cortical surface that has been extracted from
found by scanning through the possible set of locations
an anatomical MR image of the subject [62]. Following
and orientations, finding each source in turn, rather than
segmentation of the MR volume, dipolar sources are
searching simultaneously for all sources. By evaluating
placed at each node of a triangular tessellation of the sur-
J(r ,Θ) on a predefined set of grid points and then plot-
face of the cortical mantle. Since the apical dendrites that
ting its reciprocal, a “MUSIC” map is readily obtained
produce the measured fields are oriented normal to the
with p peaks at or near the true locations of the p sources.
surface, we can further constrain each of these elemental
Although we do not show the details here, (16) can be dipolar sources to be normal to the surface. The highly
modified to factor the dipole orientation out of the cost convoluted nature of the human cortex requires that a
function. In this way, at each location we can test for the high-resolution representation contains on the order of
presence of a source without explicitly considering orien- ten to one hundred thousand dipole “pixels.” The inverse
tation. If a source is present, a simple generalized problem is therefore hugely underdetermined and imag-
eigenanalysis of a 3 × 3 matrix is sufficient to compute the ing requires the use of either explicit or implicit con-
dipole orientation [43], [58]. Once all of the sources are straints on the allowed current source distributions.
found, their time series can be found, as in the Typically, this has been accomplished through the use of
least-squares approach, as S T = A + M where is the regularization or Bayesian image estimation methods.
pseudoinverse of the gain matrix corresponding to the
sources found in the MUSIC search.
Recursively applied and projected MUSIC Bayesian Formulation of the Inverse Problem
(RAP-MUSIC) is a recent improvement to the original For purposes of exposition, we will describe imaging
MUSIC scanning method, which refines the MUSIC cost methods from a Bayesian perspective. Consider the prob-
function after each source is found by projecting the sig- lem of estimating the matrix S of dipole amplitudes at
nal subspace away from the gain vectors a(ri ,Θ i ) corre- each tessellation element from the spatio-temporal data
sponding to the sources already found [59]. Other matrix M, which are related in the noiseless case by
extensions of MUSIC for MEG and EEG applications in- M = AS T . The ith row of S contains the amplitude image

across the cortex at time i. From Bayes theorem, the pos- where C S−1 is the inverse spatial covariance of the image;
terior probability for the amplitude matrix S conditioned this model assumes that the image is independent from one
on the data M is given by time sample to the next. The corresponding energy func-
tion U(S) is quadratic in S and the minimum is given by
p( M / S) p(S)
p(S / M ) =
p( M ) (17) S T = WW T A T ( AWW T A T + λI) −1 M = Fλ M (23)
where p( M / S) is the conditional probability for the data where we have factored C S−1 = WW T . We note that for
given the image and p(S) is a prior distribution reflecting this case, the posterior is Gaussian and the MAP estima-
our knowledge of the statistical properties of the un- tor is equivalent to the minimum mean squared error esti-
known image. While Bayesian inference offers the poten- mator or Wiener solution.
tial for a full statistical characterization of the sources We can also interpret (21) as a Tikhonov regularized
through the posterior probability [63], in practice images form of the inverse problem [65], [66], where the first
are typically estimated by maximization of the posterior term measures the fit to the data and the last is a regulariz-
or log-posterior probability: ing function that measures smoothness of the image. The
scalar λ is the regularization parameter that can be chosen
S = arg max p( M|S) p(S) ≡ arg max ln p( M|S) + ln p(S). using cross-validation methods or the L-curve. Within this
S S (18) regularized interpretation of (21), several forms of W have
been proposed for MEG/EEG imaging applications:
The term p(M / S) is the log likelihood for the data ▲ i) the identity matrix which produces a regularized
that depends on the forward model and the true source minimum norm solution [67];
distribution. Typically, MEG and EEG data are assumed ▲ ii) the column normalized minimum norm in which W
to be corrupted with additive Gaussian noise that we as- is a diagonal matrix with elements equal to the norm of
sume here is spatially and temporally white (generaliza- the corresponding column of A [68];
tions for colored noise are straightforward). The log ▲ iii) W computes a spatial derivative of the image of first
likelihood is then simply given by, within a constant, order [69] or Laplacian form [70];
1 2 ▲ iv) W is diagonal with elements equal to some estimate
ln p( M|S) = M − AS T . of the source power at that location, which may be com-
2σ 2 F
(19) puted from the output of a beamformer or MUSIC scan
evaluated for each dipole pixel in turn [62], [71].
The prior is a probabilistic model that describes our ex-
The underdetermined nature of the inverse problem in
pectations concerning the statistical properties of the
MEG/EEG is such that these linear methods produce
source for which we will assume an exponential density
very low-resolution solutions. Focal cortical sources tend
1 to spread over multiple cortical sulci and gyri. In some ap-
p(S) = exp{−β f (S)} plications, this may be sufficient to draw useful inferences
z (20)
from the resulting images. However, the images formed
where β and z are scalar constants and f (S) is a function of do not reflect the generally sparse focal nature of event-re-
lated cortical activation that is visualized using the other
the image S. This form encompasses both multivariate
functional imaging modalities of PET and fMRI. In an at-
Gaussian models and the powerful class of Gibbs distri-
tempt to produce more focal sources, the FOCUSS
butions or Markov random field models [64]. Com-
method [72] uses an iterative reweighting scheme in
bining the log likelihood and log prior gives the general
which the diagonal weight matrix W is updated at each it-
form of the negative log posterior whose minimization
eration to equal the magnitude of the current image esti-
yields the maximum a posteriori or MAP estimate:
mate. This approach does indeed produce sparse sources,
2 but can be highly unstable with noisy data.
U (S) = M − AS T + λf (S), An interesting approach to the interpretation of mini-
(21)
mum norm images formed using (23) was proposed by
where λ = 2βσ 2 . We can now give a brief overview of the Dale et al. [56] in which an image of SNR is computed by
imaging methods as special cases of minimization of the normalizing each pixel value computed using (23) with an
energy function in (21). estimate of the noise sensitivity of that pixel, i.e., for the
case of white Gaussian noise, each value in S T is normal-
ized by the noise sensitivity given by the corresponding di-
Linear Imaging Methods agonal elements of Fλ FλT . This has the interesting
In the case of a zero mean Gaussian image, the log prior property of generally reducing the amount by which activ-
has the form: ity spreads across multiple sulci and gyri when compared
to the standard minimum norm image; these images can
{ }
f (S) = tr SC S−1 S T , (22) also be used to make statistical inferences about the proba-

bility of a source being present at each location. This is an the MRF model can capture local interaction properties
alternative to the statistical parametric mapping described between image pixels and their neighbors. The total num-
in our discussion of beamforming methods. ber, J, of these functions depends on the number of pixels
and the number of different ways in which they are al-
Non-Gaussian Priors lowed to interact with their neighbors. Among the sim-
In an attempt to produce more physiologically plausible plest MRF image models are those in which each of the
images than can be obtained using linear methods, a large potential functions involves a pair of neighboring pixel
number of researchers have investigated alternative meth- values. To model smoothness in images an appropriate
ods that can collectively be viewed as selecting alternative choice of potential function might be the squared differ-
(nonquadratic) energy functions f (S) in (21). From a ence between these neighboring pixels.
regularization perspective, these have included entropy In the case of MEG/EEG, the model should reflect the
metrics and L p norms with values of p < 2, i.e., f (S) = S p observation that cortical activation appears to exhibit a
[73]. For the latter case, solutions will become increas- sparse focal structure, i.e., during an event related
ingly sparse as p is reduced. For the special case of p = 1, MEG/EEG study, most of the cortex is not involved in the
the problem can be modified slightly to be recast as a lin- response, and those areas correspond to focal regions of ac-
ear program. This is achieved by replacing the quadratic tive cell assemblies. To capture these properties, a highly
log-likelihood term with a set of underdetermined linear nonconvex potential function defined on the difference be-
inequality constraints, where the inequalities reflect ex- tween each pair of neighboring pixel values was used in
pected mismatches in the fit to the data due to noise. The [75]. This prior has the effect of favoring the formation of
L1 cost can then be minimized over these constraints us- small discrete regions of active cortex surrounded by re-
ing a linear simplex algorithm. The attraction of this gions of near-zero activity. An alternative model was pro-
approach is that the properties of linear programming posed in [76] where a binary random field, x, was used to
problems guarantee that there exists an optimal solution indicate whether each dipole pixel was either active ( x = 1)
for which the number of nonzero pixels does not exceed or inactive ( x = 0). A MRF was defined on this binary field
the number of constraints, or equivalently the number of to capture the two desired properties of sparseness and spa-
measurements. Since the number of pixels far outweighs tial clustering of active pixels; the parameters of this prior
the number of measurements, the solutions are therefore could then be adjusted to achieve differing degrees of
guaranteed to be sparse. This idea can be taken even fur- sparseness and clustering [76]. Interactions between
ther by using the L p quasi-norm for values of p < 1. In this neighboring pixels can be described in both space and
case, it is possible to show that there exists a value 0 < p < 1 time. In [75] for instance, temporal smoothing is included
for which the resulting solution is maximally sparse [68]. in the prior in addition to the spatial sparseness term.
An alternative to the use of simple algebraic forms for
the energy function f (S) is to explic-
itly define a prior distribution that
captures the desired statistical prop-
erties of the images. This can be done
using the class of Markov random
field (MRF) models [64], [74].
MRFs are a powerful framework,
which have been extensively investi- (a) (c) (e) (f)
gated in image restoration and recon-
struction for statistical modeling of a
range of image properties. A key
property of MRFs is that their joint
statistical distribution can be con-
structed from a set of potential func-
tions d efined on a loc a l (b) (d) (g) (h)
neighborhood system. Thus, the en- ▲ 3. Examples of surfaces extracted from high-resolution MR images. The following sur-
ergy function f (S) for the prior can faces were all extracted using an automated method described in [34]: (a) high-resolu-
be expressed as tion brain surface, (b) smoothed brain surface, (c) skull surface, (d) scalp surface. The
surfaces in (b), (c), and (d) are used as input to a boundary element code for comput-
J
ing forward MEG and EEG fields. The remaining figures are high-resolution cortical sur-
f (S) = ∑ Φ j (S), faces extracted from the same MR image and corrected to have the topology of a
j =1 (24) sphere. These can be used for cortically-constrained MEG or EEG imaging: (e) high res-
olution cortical surface, (f) smoothed representation obtained using relaxation methods
where Φ j (S) is a function of a set of similar to that described by [62] to allow improved visualization of deep sulcal features,
dipole pixel sites on the cortex that (g) high resolution and (h) smoothed representations of the cortex with approximate
are all mutual neighbors. In this way, measures of curvature overlaid. Figure courtesy of David W. Shattuck [90].

The MRF-based image priors lead to nonconvex Limitations of Imaging Approaches and Hybrid Alternatives
[75] and integer [76] programming problems in com- The imaging approaches are fundamentally limited by the
puting the MAP estimate. Computational costs can be huge imbalance between the numbers of spatial measure-
very high for these methods since although the priors ments and dipole-pixels. As we have seen, methods to
have computationally attractive neighborhood struc- overcome the resultant ambiguity range from mini-
tures, the posteriors become fully coupled through the mum-norm based regularization to the use of physiologi-
likelihood term. Furthermore, to deal with noncon- cally based statistical priors. Nonetheless, we should
vexity and integer programming issues some form of emphasize that the class of images that provide reasonable
deterministic or stochastic annealing algorithms must fits to the data is very broad, and selection of the “best”
be used [77]. image within the class is effectively done without regard
▲ 4. MEG, from modeling to imaging. Upper frame: The head

modeling step consists in designing spherical (3-shell spherical
volume conductor model; left) or realistic head models from the
subject’s anatomy using the individual MRI volume (piecewise
homogenous BEM model, right). Right frame: Three representa-
tive imaging approaches were applied to identify the MEG gen-
erators associated to the right-finger movement data set
introduced in Fig. 2. Top: The LS-fit approach produced a sin-
gle-dipole model located inside the contralateral (left) central
sulcus. Location is adequate but this model is too limited to
make any assumption regarding the true spatial extension of
the associated neural activation. The dipole time series (not
shown here) indicated little premotor activation and much stron-
ger somatosensory activity about 20 ms after the movement on-
set. Center: Minimum-Norm imaging of the cortical current
density map has much lower spatial resolution. The estimated
neural activation is spread all over the central sulcus region. Its
extension to the more remote gyral crowns is certainly
artifactual. The source time series in the central sulcus area (not
shown here) revealed similar behavior as in the LS-fit case. Bot-
tom: RAP-MUSIC modeling followed by cortical remapping: the
RAP-MUSIC approach generated a 3-source model: one in the
somatosensory regions of each hemispheres and one close to
the post-supplementary motor area (SMA). Cortical remapping
of the contralateral source revealed activation in the
omega-shaped region of the primary sensory and motor hand
areas (contralateral sensori-motor cortex, cSM). The cortical
patch equivalent to the ipsilateral source was located in the
ispsilateral somato-sensory region (iSS). Time series of the corti-
cal activations were extracted in the [-400, 100] ms time win-
dow (bottom left). Sustained pre-motor activation occurred in all
the above-mentioned areas; but only the SMA and cSM time
series had clear peaks at about 20 ms following the movement
onset, revealing motor activation of the contralateral finger and
its associated somatosensory feedback. Premotor activation in
the iSS could be related to active control of the immobility of the
ipsilateral fingers.

to the data. In contrast, the dipolar and multipolar meth- Signal Denoising and Blind Source Separation
ods control this ambiguity through a more explicit speci- An area of intense interest at the moment is the use of
fication of the source model. This may lead to improved blind source separation and independent component
confidence in the estimated sources, but at the potential analysis (ICA) methods for analysis of EEG and MEG
cost of missing sources that do not conform to the chosen data. Electrophysiological data is often corrupted by ad-
model, and to the added complexity of interpreting the ditive noise that includes background brain activity, elec-
resulting solutions. trical activity in the heart, eye-blink and other electrical
Recently we have been exploring the idea of remap- muscle activity, and environmental noise. In general,
ping estimated dipolar and multipolar solutions onto cor- these signals occur independently of either a stimulus, or
tex as a hybrid combination of the parametric and the resultant event-related responses. Removal of these
imaging approaches [78]. In this way, we can first rapidly interfering signals is therefore an ideal candidate for ICA
find a solution to the inverse problem using, for example, methods that are based on just such an independence as-
the MUSIC scanning method. We then fit each source in sumption. Successful demonstrations of denoising have
turn to the cortex by solving a local imaging problem to been published using mutual information [84], entropy
compute an equivalent patch of activated cortex whose [85], and fourth-order cumulant [86] based approaches.
magnetic fields or scalp potentials match those of the esti- These methods perform best when applied to raw
mated dipole or multipole (see Fig. 4 for an illustration (unaveraged) data; one enticing aspect of this approach is
on the data presented in Fig. 2). that after noise removal, it may be possible to see event-
A second hybrid approach to source estimation draws related activity in the unaveraged denoised signals [87].
elements from the imaging and parametric approaches by This is important since much of the information in
specifying a prior distribution consisting of a set of acti- MEG/EEG data, such as signals reflecting non time-
vated cortical regions of unknown location, size and ori- locked synchronization between different cell assemblies
entation. By constructing and sampling from a posterior [88] may be lost during the averaging process.
distribution using Markov chain Monte Carlo methods, In addition to denoising, ICA has also been used to de-
Schmidt et al. [79] are able to investigate the parameter compose MEG/EEG data into separate components,
space for this model and provide estimates, together with each representing physiologically distinct processes or
confidence values, of the true source distribution. As with sources. In principle, localization or imaging methods
the other physiologically based Bayesian models, this ap- could then be applied to each of these components in
proach has high computational costs. turn. This decomposition is based on the underlying as-
sumption of statistical independence between the activa-
tions of the different cell assemblies involved, which still
Emerging Signal Processing Issues remains to be validated experimentally. This approach
Combining fMRI and MEG/EEG could lead to interesting new ways of investigating data
One of the most exciting current challenges in functional and developing new hypotheses for methods of neural
brain mapping is the question of how to best integrate communications. This is currently a very active and po-
data from different modalities. Since fMRI gives excellent tentially fruitful research area.
spatial resolution with poor temporal resolution, while
MEG/EEG gives excellent temporal resolution with poor
Conclusion and Perspectives
spatial resolution, the data could be combined to provide
insight that could not be achieved with either modality As we have attempted to show, MEG/EEG source imag-
alone. One manner in which this has been done is to find ing encompasses a great variety of signal modeling and
regions of activation in fMRI images and use these to in- processing methods. We hope that this article serves as an
fluence the formation of activated areas on the introduction that will help to attract signal-processing re-
MEG/EEG images. This can be done by modifying the searchers to explore this fascinating topic in more depth.
covariance matrix C S−1 in (22) so that activated pixels in We should emphasize that this article is not intended to
the fMRI images are more likely to be active in the MEG be a comprehensive review, and for the purposes of pro-
images [56]. This approach works exceedingly well when viding a coherent introduction, we have chosen to pres-
the areas of activation in the two studies actually corre- ent the field from the perspective of the work that we have
spond, but can lead to erroneous results if areas actively done over the last several years.
contributing to the MEG/EEG signal do not also pro- The excellent time resolution of MEG/EEG gives us a
duce activation in fMRI, or if hemodynamic response im- unique window on the dynamics of human brain func-
aged with fMRI occurs at some distance from the tions. Though spatial resolution is the Achilles’ heel of
electrical response measured with MEG/EEG [80]-[82]. this modality, future progress in modeling and applying
The non-Gaussian Bayesian methods could be similarly modern signal processing methods may prove to make
modified to include fMRI information but would be sub- MEG/EEG a dependable functional imaging modality.
ject to the same kind of errors. This issue remains an open Potential advances in forward modeling include better
research problem [83]. characterization of the skull, scalp and brain tissues from

MRI and in vivo estimation of the inhomogeneous and Richard M. Leahy received the B.Sc. and Ph.D. degrees in
anisotropic conductivity of the head. Progress in inverse electrical engineering from the University of Newcastle
methods will include methods for combining MEG/EEG upon Tyne, U.K., in 1981 and 1985, respectively. In
with other functional modalities and exploiting signal 1985 he joined the University of Southern California
analysis methodologies to better localize and separate the (USC), where he is currently a Professor in the Depart-
various components of the brain’s electrical responses. Of ment of Electrical Engineering-Systems and Director of
particular importance are methods for understanding the the Signal and Image Processing Institute. He holds joint
complex interactions between brain regions using sin- appointments with the Departments of Radiology and
gle-trial signals to investigate transient phase synchroni- Biomedical Engineering at USC and is an Associate
zation between sensors [88] or directly within the Member of the Crump Institute for Molecular Imaging
MEG/EEG source map [89]. and the Laboratory for Neuro Imaging at UCLA. He is
an Associate Editor of the IEEE Transactions on Medical
Imaging and and was Co-Chair of the 2001 Conference
Acknowledgments on Information Processing in Medical Imaging. His re-
The authors are grateful to Marie Chupin and David W. search interests lie in the application of signal and image
Shattuck for their help in preparing the illustrations. This processing theory to biomedical inverse problems. His
work was supported in part by the National Institute of current research involves the development of methods
Mental Health under Grant R01-MH53213, by the Na- for anatomical and functional imaging with applications
tional Foundation for Functional Brain Imaging, Albu- in neuroimaging, oncology, and gene expression.
querque, New Mexico, and by Los Alamos National
Laboratory, operated by the University of California for References
the United States Department of Energy, under Contract
[1] S.R. Cherry and M.E. Phelps, “Imaging brain function with positron emis-
W-7405-ENG-35.
sion tomography,” in Brain Mapping: The Methods, A.W. Toga and J.C.
Mazziotta, Eds. New York: Academic, 1996.
Sylvain Baillet graduated in applied physics from the [2] A. Villringer and B. Chance, “Noninvasive optical spectroscopy and imag-
Ecole Normale Supérieure, Cachan and in signal pro- ing of human brain function,” Trends Neurosci., vol. 20, pp. 435-442, 1997.
cessing from the University of Paris-Sud, in 1994. In [3] A.M. Howseman and R.W. Bowtell, “Functional magnetic resonance imag-
1998, he completed the Ph.D. program in electrical en- ing: Imaging techniques and contrast mechanisms,” Philos. Trans. R. Soc.
gineering from the University of Paris-Sud at the Insti- London B, Biol. Sci., vol. 354, pp. 1179-1194, 1999.
tute of Optics, Orsay, and at the Cognitive [4] P.L. Nunez, Electric Fields of the Brain. New York: Oxford, 1981.
Neuroscience and Brain Imaging Laboratory at La
[5] C. Baumgartner, E. Pataraia, G. Lindinger, and L. Deecke,
Salpêtrière Hospital, Paris. From 1998 to 2000 he was a “Magnetoencephalography in focal epilepsy,” Epilepsia, vol. 41, suppl. 3,
Post-Doctoral Research Associate with the NeuroIm- pp. S39-47, 2000.
aging group at the Signal and Image Processing Insti- [6] C. Del Gratta, V. Pizzella, K. Torquati, and G.L. Romani, “New trends in
tute, University of Southern California, Los Angeles. magnetoencephalography,” Electroencephalogr. Clin. Neurophysiol. Suppl., vol.
He is now a Research Scientist with the National Center 50, pp. 59-73, 1999.
for Scientific Research (CNRS) and the Cognitive Neu- [7] R. Hari and N. Forss, “Magnetoencephalography in the study of human
roscience and Brain Imaging Laboratory, La Salpêtrière somatosensory cortical processing,” Philos. Trans. R. Soc. London B, Biol. Sci.,
Hospital, Paris, France. His research interests involve vol. 354, pp. 1145-1154, 1999.
methodological and modeling issues in brain functional [8] V. Jousmaki, “Tracking functions of cortical networks on a millisecond
imaging. timescale,” Neural Netw., vol. 13, pp. 883-889, 2000.
[9] M. Reite, P. Teale, and D.C. Rojas, “Magnetoencephalography: Applica-

tions in psychiatry,” Biol. Psychiatry, vol. 45, pp. 1553-1563, 1999.
John C. Mosher received his bachelor’s degree in electrical
enginering with highest honors from the Georgia Insti- [10] M. Hämäläinen, R. Hari, R. Ilmoniemi, J. Knuutila, and O. Lounasmaa,
tute of Technology in 1983. From 1979-1983 he was “Magnetoencephalography. Theory, instrumentation and applications to
the noninvasive study of human brain function,” Rev. Mod. Phys., vol. 65,
also a cooperative education student with Hughes Air- pp. 413-497, 1993.
craft Company in Fullerton, CA. From 1983-1993, he
[11] M. Piccolino, “Luigi Galvani and animal electricity: Two centuries after
worked at TRW in Los Angeles. While at TRW, he re- the foundation of electrophysiology,” Trends Neurosci., vol. 20, pp.
ceived his M.S. and Ph.D. degrees in electrical engineer- 443-448, 1997.
ing from the Signal and Image Processing Institute of the
[12] A. Hodgkin, The Conduction of the Nervous Impulse. London, U.K.: Liver-
University of Southern California. Upon graduation, he pool Univ. Press, 1964.
joined the Los Alamos National Laboratory, New Mex-
[13] P. Gloor, “Neuronal generators and the problem of localization in electro-
ico, where he researches the forward and inverse model- encephalography: Application of volume conductor theory to electroen-
ing problems of electrophysiological recordings. His cephalography,” J. Clin. Neurophysiol., vol. 2, pp. 327-354, 1985.
interests also include the general source localization and [14] P.L. Nunez and R.B. Silberstein, “On the relationship of synaptic activity
imaging problems, both in neuroscience work and in to macroscopic measurements: Does co-registration of EEG with fMRI
other novel applications of sensor technology. make sense?,” Brain Topogr., vol. 13, pp. 79-96, 2000.

[15] F.J. Varela, “Resonant cell assemblies: A new approach to cognitive func- [35] J.J. Ermer, J.C. Mosher, S. Baillet, and R.M. Leahy, “Rapidly re-comput-
tions and neuronal synchrony,” Biol. Res., vol. 28, pp. 81-95, 1995. able EEG forward models for realistic head shapes,” Phys. Med. Biol., vol.
46, no. 4, pp. 1265-1281, Apr. 2001.
[16] C.D. Tesche and J. Karhu, “THETA oscillations index human
hippocampal activation during a working memory task,” Proc. Natl. Acad. [36] L. Geddes and L. Baker, “The specific resistance of biological materi-
Sci. USA, vol. 97, pp. 919-924, 2000. als—A compendium of data for the biomedical engineer and physiologist,”
[17] C.D. Tesche and J. Karhu, “Somatosensory evoked magnetic fields arising Med. Biol. Eng., vol. 5, pp. 271-293, 1967.
from sources in the human cerebellum,” Brain Res., vol. 744, pp. 23-31, [37] T.C. Ferree, K.J. Eriksen, and D.M. Tucker, “Regional head tissue con-
1997. ductivity estimation for improved EEG analysis,” IEEE Trans. Biomed. Eng.,
[18] C.E. Tenke, C.E. Schroeder, J.C. Arezzo, and H.G. Vaughan, “Interpreta- vol. 47, pp. 1584-1592, 2000.
tion of high-resolution current source density profiles: A simulation of [38] S. Goncalves, J.C. de Munck, R.M. Heethaar, F.H. Lopes da Silva, and
sublaminar contributions to the visual evoked potential,” Exp. Brain Res., van B.W. Dijk, “The application of electrical impedance tomography to re-
vol. 94, pp. 183-192, 1993. duce systematic errors in the EEG inverse problem—A simulation study,”
[19] R.R. Llinas, U. Ribary, D. Jeanmonod, E. Kronberg, and P.P. Mitra, Physiol. Meas., vol. 21, pp. 379-393, 2000.
“Thalamocortical dysrhythmia: A neurological and neuropsychiatric syn- [39] Y. Yukawa, N. Iriguchi, and S. Ueno, “Impedance magnetic resonance
drome characterized by magnetoencephalography,” Proc. Natl. Acad. Sci. imaging with external AC field added to main static field,” IEEE Trans.
USA, vol. 96, pp. 15222-15227, 1999. Magn., vol. 35, pp. 4121-4123, 1999.
[20] R. Vigario and E. Oja, “Independence: A new criterion for the analysis of [40] H.A. Rowley, P.E. Grant, and T.P. Roberts, “Diffusion MR imaging.
the electromagnetic fields in the global brain?,” Neural Netw., vol. 13, pp. Theory and applications,” Neuroimag. Clin. N. Am., vol. 9, pp. 343-361,
891-907, 2000. 1999.
[21] P.A. Karjalainen, J.P. Kaipio, A. Koistinen, and M. Vauhkonen, [41] D.S. Tuch, V.J. Wedeen, A.M. Dale, J.S. George, and J.W. Belliveau,
“Subspace regularization method for the single trial analysis of evoked po- “Conductivity mapping of biological tissue using diffusion MRI,” Ann. N.
tentials,” IEEE Trans. Biomed. Eng., vol. 46, pp. 849-860, 1999. Y. Acad. Sci., vol. 888, pp. 314-316, 1999.
[22] D. Cohen, “Magnetoencephalography: Evidence of magnetic fields pro-
[42] M. Scherg and D. von Cramon, “Two bilateral sources of the late AEP as
duced by alpha rhythm currents,” Science, vol. 161, pp. 664-666, 1972.
identified by a spatio-temporal dipole model,” Electroencephalogr. Clin.
[23] S. Baillet, L. Garnero, G. Marin, and J.P. Hugonin, “Combined MEG and Neurophysiol., vol. 62, pp. 32-44, 1985.
EEG source imaging by minimization of mutual information,” IEEE Trans.
[43] J. Mosher, P. Lewis, and R. Leahy, “Multiple dipole modeling and local-
Biomed. Eng., vol. 46, pp. 522-534, 1999.
ization from spatio-temporal MEG data,” IEEE Trans. Biomed. Eng., vol.
[24] D.B. Geselowitz, “On the magnetic field generated outside an 39, pp. 541-557, 1992.
inhomogeneous volume conductor by internal current sources,” IEEE
[44] G.H. Golub and C.F. Van Loan, Matrix Computation. Baltimore, MD:
Trans. Magn., pp. 346-347, 1970.
The John Hopkins Univ. Press, 1983.
[25] T. Katila and P. Karp, “Magnetocardiography: Morphology and multipole
[45] G.H. Golub and V. Pereyra, “The differentiation of pseudo-inverses and
presentations,” Biomagnetism, an Interdisciplinary Approach, S.J. Williamson,
nonlinear least squares problems whose variables separate,” SIAM J. Numer.
G.L. Romani, L. Kaufman, and I. Modena, Eds. New York: Plenum, pp.
Anal., vol. 10, pp. 413-432, Apr. 1973.
237-263, 1983.
[46] K. Uutela, M. Hamalainen, and R. Salmelin, “Global optimization in the
[26] J.C. Mosher, R.M. Leahy, D.W. Shattuck, and S. Baillet, “MEG source
localization of neuromagnetic sources,” IEEE Trans. Biomed. Eng., vol. 45,
imaging using multipolar expansions,” in Proc. 16th Conf. Information Pro-
pp. 716-723, 1998.
cessing in Medical Imaging, IPMI’99, Visegrád, Hungary, June/July 1999
(Lecture Notes in Computer Science), A. Kuba, M. Sámal, and A. [47] C.C. Wood, “Application of dipole localization methods to source identi-
Todd-Pokropek, Eds. New York: Springer-Verlag, June/July 1999, pp. fication in human evoked potentials,” Ann. N.Y. Acad. Sci., vol. 388, pp.
15-28. 139-155, 1982.
[27] G. Nolte and G. Curio, “Current multipole expansion to estimate lateral [48] M. Huang, C. Aine, S. Supek, E. Best, D. Ranken, and E. Flynn,
extent of neuronal activity: A theoretical analysis,” IEEE Trans. Biomed. “Multi-start downhill simplex method for spatio-temporal source localiza-
Eng., vol. 47, pp. 1347-1355, 2000. tion in magnetoencephalography,” Electroencephalogr. Clin. Neurophysiol.,
[28] J.C. Mosher, R.M. Leahy, and P.S. Lewis, “EEG and MEG: Forward so- vol. 108, pp. 32-44, 1998.
lutions for inverse methods,” IEEE Trans. Biomed. Eng., vol. 46, pp. [49] B.D. van Veen and K.M. Buckley, “Beamforming: A versatile approach to
245-259, 1999. spatial filtering,” IEEE ASSP Mag. (see also IEEE Signal Processing Mag.),
[29] J. de Munck, “The potential distribution in a layered spheroidal volume vol. 5, pp. 4-23, Apr. 1988.
conductor,” J. Appl. Phys., vol. 64, pp. 464-470, 1988. [50] H. Krim and M. Viberg, IEEE Signal Processing Mag., vol. 13, pp. 67 -94,
[30] Z. Zhang, “A fast method to compute surface potentials generated by di- July 1996.
poles within multilayer anisotropic spheres,” Phys. Med. Biol., vol. 40, pp. [51] B.D. van Veen, W. van Dronglen, M. Yuchtman, and A. Suzuki, “Local-
335-349, 1995. ization of brain electrical activity via linearly constrained minimum variance
[31] J. Sarvas, “Basic mathematical and electromagnetic concepts of the spatial filtering,” IEEE Trans. Biomed. Eng., vol. 44, pp. 867-880, 1997.
biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11-22, 1987. [52] M. Spencer, R. Leahy, J. Mosher, and P. Lewis, “Adaptive filters for mon-
[32] R.M. Leahy, J.C. Mosher, M.E. Spencer, M.X. Huang, and J.D. Lewine, itoring localized brain activity from surface potential time series,” in Conf.
“A study of dipole localization accuracy for MEG and EEG using a human Rec. Twenty-Sixth Asilomar Conf. Signals, Systems and Computers, vol. 1. Pa-
skull phantom,” Electroencephalogr. Clin. Neurophysiol., vol. 107, pp. cific Grove, CA: pp. 156-161, 1992.
159-173, 1998. [53] S.E. Robinson and J. Vrba, “Functional neuroimaging by Synthetic Aper-
[33] B. Fischl and A.M. Dale, “SURFER: A software package for cortical surface- ture Magnetometry (SAM),” in Recent Advances in Biomagnetism, T.
based analysis and visualization,” NeuroImage, vol. 9, no. 6, p. S243, 1999. Yoshimoto, M. Kotani, S. Kuriki, H. Karibe, and N. Nakasato, Eds.
Sendai, Japan: Tohoku Univ. Press, 1999, pp. 302-305.
[34] D.W. Shattuck, S.R. Sandor-Leahy, K.A. Schaper, D.A. Rottenberg, and
R.M. Leahy, “Magnetic resonance image tissue classification using a partial [54] K.S. Sekihara, S.S. Nagarajan, D. Poeppel, and Y. Miyashita, “Spatio-tem-
volume model,” NeuroImage, vol. 13, no. 5, pp. 856-876, May 2001. poral activities of neural sources from magnetoencephalographic data using

a vector beamformer,” in Proc. ICASSP-01, Salt Lake City, UT, 2001., pp. [73] K. Matsuura and Y. Okabe, “Selective minimum-norm solution of the
2021-2026. biomagnetic inverse problem,” IEEE Trans. Biomed. Eng., vol. 8, pp.
608-615, 1995.
[55] J. Gross, J. Kujala, M. Hamalainen, L. Timmermann, A. Schnitzler, and
R. Salmelin, “Dynamic imaging of coherent sources: Studying neural inter- [74] S.Z. Li, Markov Random Field Modeling in Computer Vision. New York:
actions in the human brain,” Proc. Nat. Acad. Sci. USA, vol. 98, pp. Springer-Verlag, 1995.
694-699, 2001. [75] S. Baillet and L. Garnero, “A Bayesian approach to introducing
[56] A.M. Dale, A.K. Liu, B.R. Fischl, R.L. Buckner, J.W. Belliveau, J.D. anatomo-functional priors in the EEG /MEG inverse problem,” IEEE
Lewine, and E. Halgren, “Dynamic statistical parametric mapping: Com- Trans. Biomed. Eng., vol. 44, pp. 374-385, 1997.
bining fMRI and MEG for high-resolution imaging of cortical activity,” [76] J.W. Phillips, R.M. Leahy, and J.C. Mosher, “MEG-based imaging of fo-
Neuron, vol. 26, pp. 55-67, 2000. cal neuronal current sources,” IEEE Trans. Med. Imag., vol. 16, pp.
[57] R.O. Schmidt, “Multiple emitter location and signal parameter estima- 338-348, 1997.
tion,” IEEE Trans. Antennas Propagat., vol. 34, pp. 276-280, 1986. [77] S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions and
[58] J.C. Mosher and R.M. Leahy, “Recursive MUSIC: A framework for EEG the Bayesian restoration of images,” IEEE Trans. Pattern Anal. Machine
and MEG source localization,” IEEE Trans. Biomed. Eng., vol. 45, pp. Intell., PAMI-6, pp. 721-741, 1984.
1342-1354, 1998. [78] J.C. Mosher, S. Baillet, and R.M. Leahy, “EEG source localization and
imaging using multiple signal classification approaches,” J. Clin.
[59] J.C. Mosher and R.M. Leahy, “Source localization using recursively ap-
Neurophysiol., vol. 16, pp. 225-238, 1999.
plied and projected (RAP) MUSIC,” IEEE Trans. Signal Processing, vol. 47,
pp. 332-340, 1999. [79] D. Schmidt, J. George, and C. Wood, “Bayesian inference applied to the elec-
tromagnetic inverse problem,” Hum. Brain Map., vol. 7, pp. 195-212, 1999.
[60] K. Sekihara, S. Miyauchi, and H. Koizumu, “Covariance incorporated
MEG-MUSIC algorithm and its application to detect SI and SII when large [80] A. Puce, T. Allison, S.S. Spence, D.D. Spencer, and G. McCarthy, “Com-
background brain activity exists,” NeuroImage, vol. 3, no. 3, p. S92, 1996. parison of cortical activation evoked by face measured by intracranial field
potentials and functional MRI: Two case studies,” Hum. Brain Map., vol.
[61] K. Sekihara, S.S. Nagarajan, D. Poeppel, S. Miyauchi, N. Fujimaki, H.
5, pp. 298-305, 1997.
Koizumi, and Y. Miyashita, “Estimating neural sources from each time-fre-
quency component of magnetoencephalographic data,” IEEE Trans. Biomed. [81] A.J. Blood and A.W. Toga, “Optical intrinsic signal imaging responses are
Eng., vol. 47, pp. 642-653, 2000. modulated in rodent somatosensory cortex during simultaneous whisker
and forelimb stimulation,” J. Cereb. Blood Flow Metab., vol. 18, pp.
[62] A. Dale and M. Sereno, “Improved localization of cortical activity by com- 968-977, 1998.
bining EEG and MEG with MRI surface reconstruction: A linear ap-
proach,” J. Cogn. Neurosci., vol. 5, pp. 162-176, 1993. [82] S.P. Ahlfors, G.V. Simpson, A.M. Dale, J.W. Belliveau, A.K. Liu, A.
Korvenoja, J. Virtanen, M. Huotilainen, R.B. Tootell, H.J. Aronen, and
[63] J.O. Berger, Statistical Decision Theory and Bayesian Analysis, 2nd ed. New R.J. Ilmoniemi , “Spatiotemporal activity of a cortical network for process-
York: Springer-Verlag, 1985. ing visual motion revealed by MEG and fMRI,” J. Neurophysiol., vol. 82, pp.
[64] Markov Random Fields: Theory and Applications, R. Chellappa and A. Jain, 2545-2555, 1999.
Eds. San Diego, CA: Academic, 1991. [83] P.L. Nunez and R.B. Silberstein, “On the relationship of synaptic activity
to macroscopic measurements: Does co-registration of EEG with fMRI
[65] G. Demoment, “Image reconstruction and restoration: Overview of com-
make sense?,” Brain Topogr., vol. 13, pp. 79-96, 2000.
mon estimation structures and problems,” IEEE Trans. Acoust. Speech, Sig-
nal Processing, vol. 37, pp. 2024-2036, 1989. [84] T.W. Lee, M. Girolami, A.J. Bell, and T.J. Sejnowski, “A unifying infor-
mation-theoretic framework for independent component analysis,” Comput.
[66] A. Tikhonov and V. Arsenin, Solutions to Ill-Posed Problems. Washington
Math. Appl., vol. 39, pp. 1-21, 2000.
D.C.: Winston, 1977.
[85] A. Hyvarinen, “Fast ICA for noisy data using Gaussian moments,” in Proc.
[67] Y. Okada, “Discrimination of localized and distributed current dipole
1999 IEEE Int. Symp. Circuits and Systems VLSI, vol. 5. Piscataway, NJ,
sources and localized single and multiple sources,” Biomagnetism, an Inter-
1999, pp. 57-61.
disciplinary Approach, W. Weinberg, G. Stroink, and T. Katila, Eds. New
York: Pergamon 1983, pp. 266-272. [86] J.F. Cardoso, “Blind signal separation: Statistical principles,” Proc. IEEE,
vol. 86, pp. 2009-2025, 1998.
[68] B. Jeffs, R. Leahy, and M. Singh, “An evaluation of methods for
neuromagnetic image reconstruction,” IEEE Trans. Biomed. Eng., vol. 34, [87] S. Makeig, M. Westerfield, T.P. Jung, J. Covington, J. Townsend, T.J.
pp. 713-723, 1987. Sejnowski, and E. Courchesne, “Functionally independent components of
the late positive event-related potential during visual spatial attention,” J.
[69] J.Z. Wang, S.J. Williamson, and L. Kaufman, “Magnetic source imaging Neurosci., vol. 19, pp. 2665-2680, 1999.
based on the minimum-norm least-squares inverse,” Brain Topogr., vol. 5,
1993, pp. 365-371. [88] E. Rodriguez, N. George, J.P. Lachaux, J. Martinerie, B. Renault, and
F.J. Varela, “Perception’s shadow: Long-distance synchronization of human
[70] R.M. Pascual-Marqui, C. Michel, and D. Lehman, “Low resolution elec- brain activity,” Nature, vol. 397, pp. 430-433, 1999.
tromagnetic tomography: A new method for localizing electrical activity in
[89] O. David, L. Garnero, and F.J. Varela, “A new approach to the
the brain,” Int. J. Psychophysiol., vol. 18, pp. 49-65, 1994.
MEG/EEG inverse problem for the recovery of cortical phase-synchrony,”
[71] R. Greenblatt, “Probabilistic reconstruction of multiple sources in the bio- in Proc. 17th Conf. Information Processing in Medical Imaging (Lecture Notes
electromagnetic inverse problem,” Inverse Problems, vol. 9, pp. 271-284, 1993. in Computer Science, vol. LNCS2082), M. Insana and R Leahy, Eds.
Springer, 2001, pp. 272-285.
[72] I.F. Gorodnitsky, J.S. Georg, and B.D. Rao, “Neuromagnetic source im-
aging with FOCUSS: A recursive weighted minimum norm algorithm,” [90] D.W. Shattuck, “Automated segmentation and analysis of human cerebral
Electroencephalogr. Clin. Neurophysiol., vol. 95, no. 4, pp. 231-251, Oct. cortex imagery,” Ph.D. dissertation, Univ. of Southern California, Aug.
1995. 2001.

Sylvain Baillet, John C. Mosher, and Richard M. Leahy: 14 November 2001 1053-5888/01/$10.00©2001IEEE

Uploaded by

Copyright:

Available Formats

Sylvain Baillet, John C. Mosher, and Richard M. Leahy: 14 November 2001 1053-5888/01/$10.00©2001IEEE

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Sylvain Baillet, John C. Mosher, and Richard M. Leahy: 14 November 2001 1053-5888/01/$10.00©2001IEEE

Uploaded by

Copyright:

Available Formats

Sylvain Baillet, John C.

©1997 PHOTODISC, INC. AND 1989-97 TECH POOL STUDIOS, INC.

14 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 15

Pyramidal Cell Assembly

Scalp Pyramidal Cell Assembly

CNRS UPR640 - USC - LANL

16 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

each average, or by analyzing the raw unaveraged −9.76 −0.16 9.44

EEG data [20], [21].

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 17

18 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

to note that brain activity does not actually consist of dis-

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 19

20 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

be represented by a few equivalent current dipoles of un- = ( I − AA + ) M

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 21

tr{[ A(r ) C A(r )] }

22 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 23

24 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 25

▲ 4. MEG, from modeling to imaging. Upper frame: The head

26 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 27

[9] M. Reite, P. Teale, and D.C. Rojas, “Magnetoencephalography: Applica-

28 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

NOVEMBER 2001 IEEE SIGNAL PROCESSING MAGAZINE 29

30 IEEE SIGNAL PROCESSING MAGAZINE NOVEMBER 2001

You might also like