www.elsevier.com/locate/ynimg
NeuroImage 23 (2004) 806 – 826

Hierarchical Bayesian estimation for MEG inverse problem
Masa-aki Sato,a,b,* Taku Yoshioka,a Shigeki Kajihara,c Keisuke Toyama,c Naokazu Goda,dKenji Doya,a,b and Mitsuo Kawatoaa

ATR Computational Neuroscience Laboratories, Seika, Soraku, Kyoto 619-0288, Japan
CREST, Japan Science and Technology Corporation, Japanc
Technology Research Laboratory, Shimadzu Co., Seika, Soraku, Kyoto 619-0237, Japand
National Institute for Physiological Sciences, Myodaiji, Okazaki 444-8585, Japanb

Received 24 November 2003; revised 1 March 2004; accepted 22 June 2004

Source current estimation from MEG measurement is an ill-posed
problem that requires prior assumptions about brain activity and an
efficient estimation algorithm. In this article, we propose a newhierarchical Bayesian method introducing a hierarchical prior thatcan effectively incorporate both structural and functional MRI data. In
our method, the variance of the source current at each source locationis considered an unknown parameter and estimated from the observed
MEG data and prior information by using the Variational Bayesian
method. The fMRI information can be imposed as prior information on
the variance distribution rather than the variance itself so that it givesa soft constraint on the variance. A spatial smoothness constraint, thatthe neural activity within a few millimeter radius tends to be similar
due to the neural connections, can also be implemented as ahierarchical prior. The proposed method provides a unified theory todeal with the following three situations: (1) MEG with no other data,(2) MEG with structural MRI data on cortical surfaces, and (3) MEG
with both structural MRI and fMRI data. We investigated theperformance of our method and conventional linear inverse methods
under these three conditions. Simulation results indicate that our
method has better accuracy and spatial resolution than the conventional linear inverse methods under all three conditions. It is also
shown that accuracy of our method improves as MRI and fMRI
information becomes available. Simulation results demonstrate that
our method appropriately resolves the inverse problem even if fMRI
data convey inaccurate information, while the Wiener filter method isseriously deteriorated by inaccurate fMRI information.
D 2004 Elsevier Inc. All rights reserved.
Keywords: MEG/EEG inverse problem; Hierarchical Bayesian method;
Variational Bayesian method; Variance estimation; MRI and fMRI
information; Spatial smoothness constraint* Corresponding author. ATR Computational Neuroscience Laboratories, 2-2-2 Hikaridai, Seika, Soraku, Kyoto 619-0288, Japan. Fax: +81 774
95 1259.
E-mail address: masa-aki@atr.jp (M. Sato).
Available online on ScienceDirect (www.sciencedirect.com.)
1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2004.06.037

Introduction
Functional magnetic resonance imaging (fMRI) and magnetoencephalography (MEG) are the major recording means of brain
activity. fMRI records brain activity with millimeter-order spatial
resolution, but its temporal resolution is on the order of several
seconds due to slow hemodynamic responses to neural activity(Bandettini, 2000; Belliveau et al., 1991; Churchland andSejnowski, 1988; Logothetis et al., 2001; Ogawa et al., 1990).
Conversely, MEG measures brain activity with millisecond-ordertemporal resolution, but its spatial resolution is poor due to theill-posed nature of the inverse problem for estimating source
currents from the electromagnetic measurement (Nunez, 1981). In
addition, the number of MEG sensors is generally insufficient toprovide a precise reconstruction of the source current (Grave dePeralta Menendez and Gonzalez Andino, 1998). Therefore, priorinformation on the source currents is essential to solve the inverse
problem. Inverse procedures are commonly classified as dipoleand distributed source methods. The dipole method (Aine et al.,
2000; Hari, 1991; Mosher et al., 1992) approximates brain
activity by a small number of current dipoles. Although thismethod gives good estimate when the number of active areas issmall, it is difficult to determine the appropriate number of dipole
sources for complicated spatio-temporal activity. In addition,neural current distribution over the cortical surface cannot be
estimated by the dipole method. The distributed source method
assumes distributed currents in the brain (Hamalainen et al.,1993). In a linear approach to resolve the inverse problem,
several prior assumptions have been used such as the (weighted)
minimum norm method (Hamalainen and Ilmoniemi, 1994;
Hamalainen et al., 1993; Wang et al., 1992), the maximum
smoothness method (Pascual-Marqui, 1999; Pascual-Marqui etal., 1994), the minimum L1-norm method (Uutela et al., 1999),
and others (Huang et al., 1997; Toyama et al., 1999). Bayesian
methods have also been proposed to incorporate nonlinear
smoothness constraints (Baillet and Garnero, 1997). Unfortunately, the prior assumptions are still insufficient to fully resolve

M. Sato et al. / NeuroImage 23 (2004) 806–826

the ill-posed nature of the inverse problem, and thus the spatial
resolution of these methods is still rather low. Attempts have been
made to overcome these limitations by supplementing theinformation of other imaging means (Ahlfors et al., 1999; Daleand Sereno, 1993; Dale et al., 2000; Fujimaki et al., 2002; Heinze
et al., 1994; Kajihara et al., 2004; Liu et al., 1998; Phillips et al.,
2002; Schmidt et al., 1999). The MRI image gives information on
the position and orientation of the cortical dipoles, while fMRI
provides topographical information on active dipoles. Although
fMRI has high spatial resolution, it has been pointed out that the
hemodynamic signals of fMRI may not precisely correspond toneural activity due to various factors such as the effects of noise
and artifacts. In particular, the fMRI activity corresponds to an
average of several thousands of MEG time series data because of
the time resolution difference between MEG and fMRI. Whenhigher brain functions are examined, several different processes
may occur within several seconds, and these activities areaveraged out in the fMRI data. Consequently, the fMRI active
areas may include inactive current areas in the MEG time slice
data, and the fMRI activity for the active current may be smearedby the temporal averaging. The recent approaches, such as the
Wiener filter or the Bayesian method (Dale et al., 2000; Kajihara
et al., 2004; Phillips et al., 2002; Schmidt et al., 1999), use thefMRI data as prior information on the source current variance byassuming that the current variances for fMRI active dipoles arelarge compared with those for fMRI inactive dipoles. Source
current estimation based on these methods may fail if the fMRI
data contain incorrect information.
In this article, we propose a new hierarchical Bayesian method
introducing a hierarchical prior that can effectively incorporate
both structural and functional MRI data. In our method, thevariance of the source current at each source location is considered
an unknown parameter and estimated from the observed MEG dataand prior information. The fMRI information can be imposed asprior information on the variance distribution rather than thevariance itself so that it gives a soft constraint on the variance.
Therefore, our method is capable of appropriately estimating thesource current variance from the MEG data supplemented with thefMRI data, even if fMRI data convey inaccurate information.
Accordingly, our method is robust against inaccurate fMRI
information. The spatial smoothness constraint that the neuralactivity within a few millimeter radius tends to be similar due to theneural connections can also be implemented as a hierarchical prior.
Because of the hierarchical prior, the estimation problem becomes
nonlinear and cannot be solved analytically. Therefore, theapproximate posterior distribution is calculated by using the
Variational Bayesian (VB) method (Attias, 1999; Sato, 2001).
The resulting algorithm is an iterative procedure that converges
quickly because the VB algorithm is a type of natural gradientmethod (Amari, 1998) that has an optimal local convergenceproperty. The VB method also provides a model selection criterion
and is used to select the most probable surface on which the source
current is located when there is no prior information on the source
position.
This paper presents principles and performance of the newhierarchical Bayesian method in comparison with the conventional
linear inverse methods for the three different situations: (1) MEG
with no other data, (2) MEG with structural MRI data on corticalsurfaces, and (3) MEG with both structural MRI and fMRI data.
We also examined the performance of our method and the Wiener
filter method for incorrect fMRI information.

807

Methods
Bayesian and linear filter approaches for MEG inverse problem
When neural current activity occurs in the brain, it produces amagnetic field observed by MEG. The relationship between themagnetic field B = {B m jm = 1:M}1 measured by M sensors and theprimary source current J = { J n jn = 1:N} in the brain is given by
B ¼ G d J;

ð1Þ

where G = { G m,n jm = 1:M, n = 1:N} is the lead field matrix. The
lead field G m,n represents the magnetic field B m produced by thenth unit dipole current. The MEG inverse problem is to estimate
the source current J from the observed magnetic field data B.
Conventional linear inverse filters can be derived by minimizing an error function defined by
2

E ðJ Þ ¼ btB  G d Jt þ JVd 0 d J:

ð2Þ

The first term in the right hand side (r.h.s) of Eq. (2) represents thereconstruction error, which measures the difference between theobserved MEG data and the reconstructed MEG signal calculated
from the estimated source current. To resolve the ill-posed natureof the MEG inverse problem, a regularization term is added to theerror function. Depending on the choice of the regularization
matrix 0, various linear inverse methods can be derived. For
example, the minimum norm solution (Hamalainen and Ilmoniemi,
1994) can be derived if 0 is the N-by-N identity matrix. The
parameter b controls the relative weight between the reconstruction
error and the regularization term. The minimum solution of Eq. (2)is given by

1
1
1
J ¼ L d B ¼ 1d B; ð3Þ
0 d GV d G d 0 d GV þ b I M
where I M denotes the M-by-M identity matrix.
These linear inverse filters can be also derived by using theBayesian method (Dale et al., 2000; Grave de Peralta Menendez
and Gonzalez Andino, 1998; Phillips et al., 2002). The probabilistic model for the source currents can be constructed assuming
Gaussian noise for the MEG sensors. Then, the probability
distribution, that the magnetic field B is observed for a given
current J, is given by
12PðBjJ Þ~exp  btB  G d J t ;ð4Þ2
where b corresponds to the inverse of the noise variance. In theBayesian method (Gelman et al., 1995), a prior assumption on thesource currents J is imposed as a prior probability distribution
P 0(J). The Bayesian method calculates the posterior probability
distribution for the source current, P(JjB), which represents theprobability that the source current is J under the observed MEG
data B and the prior information:
PðJ jBÞ ¼

PðBjJ ÞP0 ðJ Þ:
P ðB Þ

ð5Þ

1
In the following, the notation, m = 1:M, is used to represent m =1,. . ., M.

808

M. Sato et al. / NeuroImage 23 (2004) 806–826

The marginal likelihood P(B) appearing in the denominator isdefined byZPðBÞ ¼
dJ Pð BjJ ÞP0 ðJ Þ:
ð6Þ
Maximum a posteriori (MAP) estimates find the optimal parameter(source current value) that maximizes the posterior probability
distribution P(JjB). The linear inverse filters defined by Eq. (3)can be derived as the MAP solution by assuming the Normal
distribution prior for the source current with a fixed covariancematrix 1
0 ;


1
ð7Þ
P0 ðJ Þ~exp  J V d 0 d J :2
One can easily see that the maximization of the log-posteriorlog( P(JjB)) is equivalent to the minimization of the error
function (2).
The Wiener filter solution is derived by the minimization of theestimated current error, which is defined as the error between the true
current and the estimated current. The MAP solution (3) gives the
Wiener filter solution, if the prior current covariance matrix 1
0 is
the true current covariance matrix (Dale et al., 2000). However, thetrue current covariance is unknown. It has been proposed that thefMRI information can be used to determine the prior current
covariance matrix (Dale et al., 2000; Phillips et al., 2002).
Because of the underdetermined character of the MEG
inverse problem, that is, N N M, a class of the linear inverse
methods defined by Eq. (3) satisfies the following constraint (seeAppendix A):N
X

Gn V M :ð8Þ

n¼1

The estimation gain Gn represents the nth estimated current J n
by Eq. (3) for a single unit dipole source at the nth position andshould be one for perfect reconstruction. The inequality (8) implies
that the linear inverse method cannot perfectly retrieve morecurrent dipoles than the number of sensors M.
This constraint can be avoided by using the varianceinformation. To simplify the explanation, let us assume that there
are active currents in n = 1:K and the current
 covariance matrix
1
1
is
given
as
a
diagonal
matrixwith
00 n;n ¼ 1 (for n = 1:K)and 1
0 n;n ¼ 0 (for n = (K + 1):N). The estimated currents fornull variance points J n (n = (K + 1):N) vanish from Eq. (3) withthis covariance matrix. Therefore, the estimation gains for n = (K +1):NPK also vanish, Gn ¼ 0, and the constraint (8) becomesn ¼ 1 G n V M . Consequently, the inequality does not constrain
the estimation gains for active dipoles, Gn (n = 1:K), if the numberof active dipoles K is less than the number of sensors M.
Introduction of activity dependent variance information reduces theeffective degree of freedom by suppressing the estimation gains forsmall variance points, and this makes the estimation of active
dipoles more accurate.
However, the source current estimation based on the priorcurrent variance such as the Wiener filter method may fail if thevariance information is incorrect.
Hierarchical Bayesian method
We propose a new hierarchical Bayesian method to overcome
the above difficulties. In our method, the variances (diagonalpart of the covariance matrix) are considered unknown parameters and estimated from the observed MEG data by introducinga hierarchical prior on the current variance. The fMRI
information can be imposed as prior information on the variancedistribution rather than the variance itself so that it gives a soft
constraint on the variance. The spatial smoothness constraint thatneurons within a few millimeter radius tend to fire simultaneously due to the neural interactions can also be implemented
as a hierarchical prior.
Hierarchical prior
Let us suppose a time sequence of MEG data B 1:T u {B(t)jt =1:T} is observed. The MEG inverse problem in this case is toestimate the primary source current J 1:T u {J(t)jt = 1:T} from theobserved MEG data B 1:T . We assume a Normal prior for thecurrent:
"
#
T
1 X
JVt
ð Þ d A d Jðt Þ ;ð9Þ
P0 ðJ 1:T jAÞ~exp 
2 t¼1
where A is the diagonal matrix with diagonal elements A = {a n jn =1:N}. We also assume that the current variance A1 does notchange over period T. The current inverse variance parameter A isestimated by introducing an Automatic Relevance Determination(ARD) hierarchical prior (Neal, 1996):
NP0 ðAÞ ¼ P C an jā0n ; c0na ;
n¼1



c
C ajā;c u a1 ac=ā CðcÞ1 eac=ā ;ð10Þwhere C(ajā,c) represents the Gamma
R l distribution with mean āand degree of freedom c. CðcÞ u 0 dttg  1 et is the Gamma
function.
We briefly explain the difference between a conventional
Normal prior and the hierarchical prior. In the case of the Normal
prior, the value of the inverse variance parameter a n is given asprior information. For small variance (large a n value, a L ), theprior probability distribution for the current, P 0( J n ja L ), isconcentrated near zero (Fig. 1A), so large current is highly
penalized by the prior information. For large variance (small a nvalue, a S ), P 0( J n ja S ) spreads over large values (Fig. 1B), so theestimated current likely becomes large. In the case of thehierarchical prior, the inverse variance parameter a n is considereda random variable. Therefore, there are nonzero probabilities forany value of the inverse variance parameter a n , and the priorprobability distribution for a n is given by the Gamma distribution(10). The hyperparameters c 0na control the spread of the
distribution and represent the reliability (or confidence) of thehierarchical prior. The hyperparameters ā 0n represent the meanvalue of a n in the hierarchical prior. If the prior mean a 0n is small,the prior probability P 0(a n = a S ) that a n takes a small value, a S(large variance), increases (Fig. 1C). In contrast, the priorprobability P 0(a n = a L ) that a n takes a large value, a L (small
variance), increases for large ā 0n (Fig. 1D). Therefore, the priorinformation on the current variance is imposed as a soft constraint
in the variance estimation.
When the fMRI data are not available, we have no informationon the current variance. Lack of knowledge can be represented as a
non-informative prior, which assumes that any parameter value has

M. Sato et al. / NeuroImage 23 (2004) 806–826

809

Fig. 1. Hierarchical Bayesian method. (A) Probability distribution of a Normal prior P 0 ( J n ja L ) for the current J n with a large inverse variance (small variance)parameter a L is plotted. (B) That with a small inverse variance (large variance) parameter a S . (C) Hierarchical prior distribution P 0(a n ) for the inverse varianceparameter a n with a small prior mean ā 0n is plotted. The probability for small inverse variance a S is large compared with that for large inverse variance a L . (D)That with a large prior mean ā 0n . The probability for small inverse variance a S is small compared with that for large inverse variance a L . (E) Flow chart of the
Variational Bayesian (VB) method. In the J-step, the inverse filter L(A1) is recalculated according to (28) using the estimated covariance matrix A1 in theprevious iteration. The current J 1:T is then estimated from the observed MEG data B 1:T using the inverse filter according to Eq. (29). The current variance a 1
is estimated in the a-step by using the estimated time sequence of the current J 1:T and the fMRI information according to Eq. (30). The current covariancematrix A1 is then calculated by using the current variance and the smoothing filter W according to Eq. (31). The free energy is calculated at the end of thesteps. If the free energy has not converged, the above J- and a-steps are repeated.

the same probability.2 Therefore, we use a non-informative priorfor the current inverse variance parameter a n in this case, that is,
c 0na = 0 and P 0(a n ) = a n1. When the fMRI data is available, fMRI
information is imposed as the prior mean variance ā10n weighted
1with the confidence parameter c 0na . The prior mean variance ā 0n
isincreased for the fMRI active area, so the probability that a n takes
a small value (large variance) becomes large (Fig. 1C). On the1other hand, ā 0n
is decreased for the fMRI inactive area, so theprobability that a n takes a large value (small variance) becomeslarge (Fig. 1D).

correlation function is characterized by the smoothing filter W. If
we assume the Gaussian smoothing filter for W, the correlation incurrent activity between two points decreases exponentially as thesquared distance increases. Parameters A and L control the
diagonal and the off-diagonal part of the current covariancematrix, respectively. As the variance parameter a n1 increases,
the variance of the nth current increases. As the smoothness
parameter k n1 increases, the correlation between the nth and theneighboring currents increases. Parameters A and L are estimatedby introducing the ARD hierarchical prior for them.

Smoothness prior
The smoothness constraint can be imposed as the structural
constraint on the off-diagonal part of the covariance matrix, which
enforces high correlations between neighboring current activity:
"
#
T
1 X
P0 ðJ 1:T jA; LÞ~exp J VðtÞ d A d J ðt Þ ;
ð11Þ
2 t¼1

Variational Bayesian method
Because of the hierarchical prior, the estimation problem
becomes nonlinear and cannot be solved analytically. Therefore,
the approximate posterior distribution is calculated by using the
Variational Bayesian (VB) method (Attias, 1999; Sato, 2001). In
the VB method (see Appendix B.1), the calculation of theposterior distribution is reformulated as the maximization problemof the free energy for a trial distribution that approximates the true
posterior. The optimization problem is solved by imposing
constraints on the trial distribution. The resulting algorithm is an
iterative procedure. The method alternately estimates the current
and the variance from the observed MEG data B 1:T and the priorvariance information given by the fMRI data. It has been proved(Sato, 2001) that the VB method converges into a local maximumof the free energy and is equivalent to the natural gradient method(Amari, 1998), which has an optimal local convergence property.
The maximum free energy provides a model selection criterion andis used to select the most probable surface on which the source

where the current covariance matrix 1
A is given by1
þ W d K1 d W V:
1
A ¼A

ð12Þ

A and K are the diagonal matrices with the diagonal elements
A = {a n jn = 1:N} and L = {k n jn = 1:N}, respectively. The second
term in the r.h.s. of Eq. (12) represents the correlation in current
activity between the neighboring points. The spatial profile of the2

We assume the probability distribution is flat with respect to log(a n ).

810

M. Sato et al. / NeuroImage 23 (2004) 806–826

current is located when there is no prior information on the source
position (see Surface model selection). The detailed algorithm isderived in Appendix B.
Here, we explain some of the essential equations of the VB
method for the hierarchical prior Eqs. (9) and (10). The
algorithmic flow of the VB method is summarized in Fig. 1E(see also Appendix B.4). At the start of the algorithm, the initialvalue of the estimated variance a n1 is set to the prior meanvariance ā 0n1 and the covariance matrix is calculated as
1
1
A ¼ diagðA Þ. Namely, the prior variance information givenby the fMRI data is also used as the initial value of the current
variance. Consequently, the VB current estimate at the first
iteration is equal to the Wiener filter solution. When there is nofMRI information, all of the current variances are set to the samevalue, and thus the first VB current estimate is given by theminimum norm solution in this case.
In the J-step (Fig. 1E), the inverse filter L 1
is recalculatedA1using the estimated covariance matrix A . It depends on theobserved data B 1:T through the estimated variance A1. The
current is then estimated using the inverse filter (see Eqs. (28) and(29) in Appendix B.4):


1
1
L 1
¼ 1
A
A d G V d G d A d G V þ b I Md Bðt Þ:J ðtÞ ¼ L 1
A

1

;

ð13Þ

This equation has the same form as the linear inverse Eq. (3). The
crucial difference is that the variance of the current is estimated in
our method, while it is given as prior information in the linear
inverse method.
The variance is estimated in the a-step by using the estimatedtime sequence of the current J 1:T (Fig. 1E). The update equation
for the current variance a n1 is somewhat complicated (see Eq. (30)in Appendix B.4), but it becomes simple when the algorithm
converges:
 P
T
1
2
þ T2 T1
t ¼ 1 J n ðt Þ
c
ā
0n
0na
1


an ¼
:
d G n;nc0na þ T2 L 1
A

ð14Þ

The first term in the r.h.s. numerator of Eq. (14) represents
the contribution of the prior mean variance ā1
0n . It is weighted
by a prior confidence factor c 0na . The second term is the data
term and represents the average magnitude of the estimatedcurrent. It is weighted by half of the data number T/2. Therefore,
the estimated variance is a weighted average of the data term
and the prior contribution. The covariance matrix is then1calculated as 1
A ¼ diagðA Þ. When the smoothness constraint
is imposed by Eq. (11), the covariance matrix is calculated by
Eq. (12) using the smoothing filter W and the estimatedsmoothness parameter L1 as well as A1 (see Eq. (31) inAppendix B.4). The free energy is calculated at the end of thesteps. If the free energy is not converged, the above J- and asteps are repeated.
Since the VB method estimates the current variance according
to the observed MEG data as well as the fMRI information, the VB
method can appropriately estimate the current variance value, evenwhen the fMRI data is biased toward an incorrect variance value.
We evaluated this point by using anatomically realistic simulation
data, as described in the following sections.

Surface model selection
When structural MRI data is available, the cortical surface
obtained from the MRI image data can be used as a surface model
M for current estimates and a dipole source current J n is assumed
in each vertex point on the surface M. If there is no structural
information by MRI, the shape and position of the cortical surface
is unknown. In this case, we consider a set of surface models Mwith different radii and compare the probability that the source
current is on the surface M. Since there is no prior information on
the surface position, any surface may have the same priorprobability, that is, P0 ðMÞ = constant. Then, the probabilityPðMjB1:T Þ, that the source current is on the surface M under theobserved data B 1:T, is proportional to the marginal likelihood forthe surface model M; PðB1:T jMÞ:3PðMjB1:T Þ~PðB1:T jMÞP0 ðMÞ~PðB1:T jMÞ:

ð15Þ

Namely, we can find the most probable surface M by finding
the surface with the maximum marginal likelihood, which isequivalent to finding the maximum log-marginal likelihood logð PðB1:T jMÞÞ because the logarithm is a monotonically increasing
function.
The log-marginal likelihood is approximated by the maximum
free energy for the surface model M in the VB method (seeAppendix B.1). Therefore, the most probable surface can be
determined by finding the surface model with the maximum freeenergy. This procedure gives more accurate estimation of thesource current than using the distributed current model in threedimensional space.

Results
We evaluated the performance of our VB method to resolve
the MEG inverse problem in comparison with conventional linear
inverse filters. Basic estimation capabilities were evaluated byconducting simulations for a single dipole source and a single
distributed source randomly placed on a cortical surface under
three situations: (1) MEG with no other data, (2) MEG withstructural MRI data on cortical surfaces, and (3) MEG with both
structural MRI and fMRI data. Precision of the source localization was measured by a localization error, which was definedas the Euclidean distance between the true current and theestimated current peak. Accuracy of the current amplitude
estimation for a single dipole source was measured by theestimation gain, which was defined as the ratio between theestimated current and the source current at the source position.
The spatial profile of the current estimates for a single dipole
source was also measured by a point-spread curve, which was
defined as the estimated current profile around the source dipole.
These measures were averaged over 500 cases randomly
generated on the cortical surface.
More realistic estimation performances were also evaluatedusing simulated visual cortical responses under three othersituations: (1) MEG with structural MRI data, (2) MEG withstructural MRI and correct fMRI data, and (3) MEG with structural
MRI and incorrect fMRI data. To evaluate spatio-temporal patterns

3Here, the dependence on the surface model M is explicitly
expressed.

M. Sato et al. / NeuroImage 23 (2004) 806–826

of the estimated currents, local currents were defined as the sum of
the estimated currents within the areas corresponding to the active
source currents. The precision of the estimation was assessed as theestimation gain at the source current peak and the normalized root
mean squared error (RMSE) for the temporal sequences of the local
current.
MRI-based cortical surface modelA polygon model of cortical surfaces (number of vertex points,
24,917; intervertex distance, 1.5 F 0.5 mm) was constructed basedon MRI image slices of an experimental subject (YO) using
BrainVoyager 2000 (Goebel and Max Plank Society, 2000) (Fig.
2A). A single current dipole was assumed at each vertex point
perpendicular to the cortical surface. One hundred and ninety-nine
gradiometer sensors were assumed to surround the brain (blue dots
in Fig. 2A), which reproduced the geometry of the whole-head
MEG recording device used in our laboratory (Shimadzu SBI200).
For calculation of the lead field, the brain structures wereapproximated as a sphere with a 100-mm radius, and the Sarvas
equation (Sarvas, 1987) was used. In the following simulations, weused this MRI-based cortical surface model as a forward model togenerate simulated MEG signals.

811

VB and linear inverse filter estimates for single dipole anddistributed current sources
To investigate the basic capabilities of the VB estimate, we
conducted simulations for two types of sources, that is, a single
dipole source and a single circular distributed source composed of
uniformly distributed dipoles within a 7-mm radius from the source
center on the cortical surface. Noiseless and noisy conditions,where MEG signals contain Gaussian sensor noise with N/S of 0.1,
were tested under three situations: (1) MEG with no other data, (2)MEG with structural MRI data, and (3) MEG with both structural
MRI and fMRI data. Precision of the source localization wasmeasured by the localization error, which was defined as the
Euclidean distance between the true current (center) position X 0
and the maximum (absolute value) peak position X peak of theestimated current:
Localization Error ¼ tX 0  X peak t:
A single dipole or a single distributed source was located at oneof 500 randomly selected vertices among 24,917 vertices (redpoints in Figs. 2B and C), and the localization errors for the 500
cases were averaged out, since the estimation performance depends
on the depth and orientation of the current sources. The accuracy ofFig. 2. MRI-based cortical surface model. (A) MRI-based cortical surface model obtained from MRI image data. 199 MEG sensors are also plotted as blue dots.
(B) Five hundred points randomly selected from the MRI-based cortical surface model are shown by red dots (Top view). (C) Bottom view. (D) Variable-radiussurface model for current estimates is plotted together with 199 MEG sensors.

812

M. Sato et al. / NeuroImage 23 (2004) 806–826

the amplitude estimation and the spatial profile of the current
estimates for a single dipole source were also assessed by averagesof the estimation gain and the point-spread curve, respectively. The
results were compared with those for the conventional linear
inverse filters.
Before going into the detailed results, the estimation proceduresof the VB method are explained for the three situations.
MEG with no other information
Variable-radius surface model
In the situation where no other information is available than
MEG data, we assumed a variable-radius surface model for current
estimates. The surface consisting of 2282 polygons (number ofvertex points, 1185; intervertex distance, 5.7 mm at 70 mm radius)
with a variable radius R was constructed, and a tangential current
dipole with two degrees of freedom was assumed at each vertex
point (Fig. 2D). The surface radius was determined so as tomaximize the free energy (see Surface Model Selection). We used
the hierarchical prior without smoothness constraints Eq. (9) andassumed the non-informative prior for P 0(A).

Single current source
We first explain the VB method for a single sinusoidal current
dipole at X 1 = (radius r 1 = 74.5 mm, elevation angle h 1 = 69.48,
azimuthal angle / 1 = 110.28) of the MRI-based cortical surface
model (blue dot in Fig. 3A) under the noisy condition with
Gaussian sensor noise (N/S = 0.1).
The VB method was applied for each surface with a radius R byrepeating the current (J-step) and the variance (a-step) estimations
alternately according to Appendix B.4 until the free energy
converged to the maximum (Fig. 1E). The free energy obtainedfor variable R peaked at R = 75 mm (Fig. 3C), which roughly
agreed with the radius of the current source (74.5 mm). To see theeffect of noise, the free energy curve under noiseless condition was
also shown in a blue curve of Fig. 3B. The free energy was sharply
tuned at the peak position when there was no noise. As the noise
amplitude increased, the peak profile of the free energy curve
became broader as shown in Fig. 3C. However, the peak almost
stayed at the same position and there was no local peak around the
maximum peak position in the presence of noise. Therefore, thesource depth was safely determined by the free energy maximum.
On the other hand, the reconstruction error (green curve in Fig. 3B)

Fig. 3. Single and multiple source estimation. (A) A single dipole source at X 1 = (74.5 mm, 69.48, 110.28) in the MRI-based cortical surface model is shown by
a blue dot. The other dipole source at X 2 = (57.1 mm, 47.28, 152.38) used in multiple source estimation is also shown by a green dot. (B) Free energy curve fora single source estimation under noiseless condition is plotted by a blue line. The free energy was maximized for each surface with radius R. The free energy
peaked at 75 mm, which well agreed with the depth of the current source 74.5 mm. The reconstruction error under noiseless condition is also plotted by a green
line. (C) Free energy curve for a single source estimation under noise condition. The free energy maximized for each surface with radius R is plotted. The freeenergy peaked at 75 mm. (D) The reconstruction error of the dipole fit under noise condition was calculated for each radius. A single dipole was fitted to theobserved MEG data for each radius R by constraining the radius of the dipole to R. The error had the minimum at 75 mm. (E) Peak sinusoidal current estimatedby VB method under noise condition for the surface with R = 75 mm is plotted on a flattened map using pseudocolor representation of current intensity. The
source position is depicted by a blue dot. (F) Estimated currents on the maximum free energy surface for two dipole sources are plotted on the flattened map.
Two local current peaks were identified near the angular positions of the two sources and depicted by blue and green crosses. Two local surfaces around the
local current peaks are shown by blue and green circles. (G) The free energies maximized for a set of two local surfaces with radii R 1 and R 2 are plotted. The
free energy as a function of R 1 (R 2) is plotted by a blue (green) line, where R 2 (R 1) is fixed at 57 mm (70 mm). (H) Current estimates for a set of two local
surfaces with R 1 = 75 mm and R 2 = 57 mm are shown by pseudocolor representation. The two source positions are depicted by blue and green dots.

M. Sato et al. / NeuroImage 23 (2004) 806–826

of this distributed source model under the noiseless condition was
nearly flat for the radius larger than 50 mm, so it was not a goodmeasure to estimate the source depth. This is because the magnetic
field produced by a deep focal source can be reproduced by
distributed superficial sources. We also examined the standard
single dipole fit for this problem. The reconstruction error (Fig.
3D) of the dipole model under the same noise condition (N/S =0.1) showed the similar (reversed) tuning curve as the free energy(Fig. 3C).
We then estimated the source current for the maximum freeenergy surface. The sinusoidal current peak estimated by the VB
method for the surface with R = 75 mm was correctly focusedaround the true dipole position (blue dot in Fig. 3E) with alocalization error of 2.7 mm.
Multiple current sources
To explain the VB method for a multiple dipole case, we
assumed another sinusoidal current dipole of the same current
strength at X 2 = (57.1 mm, 47.28, 152.38) (green dot in Fig. 3A), in
addition to that for the single dipole case.
The maximum free energy surface was determined in the same
way as for the single current source case. The free energy peaked at
R = 70 mm. There were two current peaks at X̂ 1 and X̂ 2 on thissurface (blue and green crosses in Fig. 3F), and local surfaces (blueand green circles) were extracted as assemblies of vertex points X
within 308 angular coordinate differences from X̂ 1 and X̂ 2, that is,(XV X̂ i )/(jXtX̂ i j) z cos (308) for i = 1,2. Then, a new surface
model was constructed as the set of the two local surfaces with theradii R 1 and R 2, which were determined by alternately maximizing
the free energy with respect to R 1 and R 2. The free energy was firstmaximized by varying R 2 while R 1 was fixed at 70 mm. The freeenergy peaked at R 2 = 57 mm (green line in Fig. 3G). Next, the freeenergy was maximized by varying R 1 while R 2 was fixed at 57
mm. The free energy peaked at R 1 = 75 mm (blue line in Fig. 3G).
In principle, it is necessary to iterate this procedure until the freeenergy converges to a maximum value. However, no increase in
the free energy was found with further iteration in this case, and R 1
and R 2 were estimated at 75 and 57 mm, respectively. Plots ofsinusoidal current peaks estimated for the two local surfaces (R 1
and R 2, 75 and 57 mm) indicate that the VB current estimates werecorrectly focused to the two true dipole sources (blue and green
dots in Fig. 3H). Localization errors for X 1 and X 2 were 2.7 and2.6 mm, respectively.
The VB method can recover multiple current sources from
MEGs if the estimated current peaks in the maximum free energy
surface can be discriminated and correspond to all of the current
sources. To specify each source position, a new surface model isconstructed from a set of local surfaces for the individual current
peaks. The VB estimates can be completed by maximizing free813

energy by alternately changing each local surface radius. To
investigate whether the VB method can discriminate two sources
placed at nearby points, simulations were conducted by randomly
generating 350 sets of two dipole sources (maximum distance, 35mm) in the MRI-based cortical surface model. When the Euclidean
distance between the two sources was small, the VB method often
estimated only a single current peak (Table 1). The percentages that
only a single source was recovered were 81.8%, 40.4%, and 14.6%
for the distance ranges 0–5, 5–10, and 10–15 mm, respectively.
When only a single current peak was estimated, the localization
error for each source was defined as the Euclidean distance
between each source and the single estimated current peak.
Localization errors for the two sources were averaged over 350
cases according to the range of the distance, that is, 0–5,..., 30–35
mm (Table 1). Localization errors slightly increased arounddistances of 5–20 mm. However, the localization errors decreased
for distances less than 5 mm, since the single estimated current
peak could approximate both sources within 5-mm localization
error. Table 1 also shows that the mean localization errors for thetwo sources remained below 7 mm across the distances of 0–35
mm. These results indicate that the minimal distance for the VB
estimates to resolve two current sources is roughly 15 mm if noother information than MEG is available.
MEG with structural MRI information
When structural MRI data were available, the MRI-based
cortical surface model for current sources (Fig. 2A) could be used
as a surface model for current estimates. The VB estimate wasconducted in two steps: (1) low-resolution current estimation using
a coarse surface model and (2) high-resolution estimation using afine surface model. The coarse surface model (number of vertex
points, 2494; average intervertex distance, 6 mm) was constructedby reducing the number of vertex points down to one-tenth of thatfor the MRI-based cortical surface model. A locally distributed
current source with a Gaussian profile was assumed at each vertexof the coarse surface model. Then, the current in the coarse surface
model Z L induces the current J = W L Z L in the MRI-based corticalsurface model, where W L denotes the Gaussian profile of the
distributed currents. The Gaussian profile is given by (W L )nj = exp((d nj /r)2), where d nj denotes the shortest path length between thenth vertex point of the MRI-based cortical surface model and thejth vertex point of the coarse surface model. The magnetic filed B
produced by the coarse model current Z L is given by B = Gd J =Gd W L d Z L . Therefore, the effective lead field matrix for Z L isgiven by Gd W L . The coarse model current Z L was estimated by
the VB method with the hierarchical prior Eq. (9). The noninformative prior was assumed for P 0(AL ), where AL was aninverse variance parameter for Z L . The VB algorithm in Appendix

Table 1Discriminability of two sources
Distance (mm)% of failure
Localization error

0–5
81.8
3.62

5–10
40.4
6.95

10–15
14.6
6.43

15–20
2.0
5.50

20–25
2.4
3.89

25–30
0.0
3.78

30–35
0.0
3.04

Discriminability of two sources. Discriminability of two sources for the VB method with no other information was investigated by randomly generating 350
sets of two dipole sources (maximum distance, 35 mm) in the MRI-based cortical surface model. The percentages that only a single source is recovered arelisted according to the range of the Euclidean distance between the two sources, that is, 0–5,. . ., 30–35 mm. Localization errors for two sources were averaged
over 350 sets according to the range of the distance and listed in millimeters.

814

M. Sato et al. / NeuroImage 23 (2004) 806–826

B.4 was modified such that G was replaced by Gd W L , and W = 0.
The current in the MRI-based cortical surface model is given by

1
1
1
J ðtÞ ¼ 1d Bðt Þ;A d G V G d A d G V þ b I M
1
1
A ¼ W L d AL d W L V;

ð16Þ

where A L1 = diag(AL1) denotes the variance matrix for Z L . The
comparison between Eqs. (13) and (16) implies that this current
model is equivalent to assuming a hierarchical prior with thecurrent covariance matrix, (W L d A L1d W LV ).
High-resolution VB estimation was conducted for local peakcurrents of the coarse surface model. As in the case of MEG for
multiple dipoles with no other information, local surfaces wereextracted from the MRI-based cortical surface model as an assemblyof vertex points within 35 mm from the vertex point of each local
current peak. The VB current estimate was performed for a finesurface model constructed by a set of extracted local surfaces using
the hierarchical prior, Eqs. (11) and (12), with a non-informative
prior. The Gaussian filter for the smoothness constraints was
assumed as W nm = exp((d nm /r)2), where d nm denotes the shortest
path length between the nth and the mth vertex point of the MRIbased cortical surface model. The full width half maximum(FWHM) of the Gaussian filter was assumed to be 8 mm. We
tested several Gaussian filters with different FWHM. The biologically plausible value, FWHM = 3.5 mm according to the point
spread map of fMRI (Engel et al., 1997), was too small in thepresent MRI-based cortical surface model. The number of points
within a 3.5-mm diameter for each vertex point was very few(average number, 2 F 1), which caused inhomogeneity in thesmoothness constraints. The value FWHM = 8 mm (averageneighbor number, 8 F 2) was selected as the compromise between
spatial resolution and smoothness of the estimation results.
The performance of the VB estimates with structural MRI
information was compared with that of the conventional linear
inverse filters. The minimum norm estimate (MN) (Hamalainen andIlmoniemi, 1994) and the minimum norm estimate with smoothness
constraints (MNC) were tested using the same MRI-based corticalsurface model. To avoid the regularization problem of determining
the regularization coefficient, simulation was conducted only for thenoiseless case. The MN estimate is defined as the linear inverse
estimate (3) with 0 = I N and b 1 = 0. The MNC estimate is defined2
as that with 1
0 nm ¼ expð ðdnm =rÞ Þ (FWHM = 8 mm) and1b= 0. This smoothness constraint with Gaussian profile iscomparable with the smoothness constraint used in the highresolution VB estimation.

MEG with MRI and fMRI information
If fMRI data is available in addition to structural MRI data, that
information may be used as hierarchical priors. fMRI signals wereassumed to uniformly spread out in a circular area of a 15-mm
radius around the point source brain activity: the fMRI signal V n at
the nth vertex point was one in the 15-mm radius around the source
current (center) position, and V n = 0 for other positions. The VB
current estimate was performed in the same way as for the highresolution VB with MRI information, except for the introduction offMRI signals into the hierarchical priors defined by Eq. (32) (seeAppendix B.4).
The performance of the VB estimates was also compared withthat of the Wiener filter estimates (Dale et al., 2000; Kajihara et al.,
2004; Phillips et al., 2002), using the same MRI and fMRI
information. For the noiseless case, the Wiener filter estimate
2
is given by Eq. (3) with 1
0 nm ¼ Vn Vm expð ðdnm =rÞ Þ(FWHM = 8 mm) and b 1 = 0.
Performance improvement by combining MEG with MRI and fMRI
Table 2 compares localization errors of the VB estimates
averaged over 500 randomly placed dipoles and distributed sourcesin the MRI-based cortical surface model. Noiseless and noisy
conditions were tested under three situations: (1) MEG with noother data (NO), (2) MEG with MRI (MRI), and (3) MEG withboth MRI and fMRI (MRI + fMRI). One can see that localization
errors in the VB estimates reduce for all cases as more information
becomes available.
Interestingly, the performance of the VB estimate is slightly
better under the noisy condition than that of the noiselesscondition. In the noisy condition, dynamic sources with sinusoidal
waveform were assumed and the time sequence of MEG data wasused to estimate the current variance in the VB method. On theother hand, static sources were assumed in the noiseless case, andonly a single set of MEG data was used to estimate the current
variance. This result shows the usefulness of the time sequence
information of MEG data.
The contribution of MRI and fMRI information to improve
localization errors for single dipole sources is also shown by scatterplots of localization errors as a function of the lead field norm
2tG n t = M(mM= 1G m,n) (Figs. 4A–C). The lead field norm
represents the magnitude of MEG produced by a single dipole
source, and it becomes small as the source location becomes deep
and the current orientation approaches the radial direction. Local-

Table 2
Localization error for single source
MethodVB

Information
Dipole source
Distributed source

No Noise
N/S = 0.1
No Noise
N/S = 0.1

MN

MNC

WeinerNO

MRI

MRI + fMRI

MRI

MRI

MRI + fMRI

4.85 (8.19)
4.64 (7.72)
5.90 (7.71)
5.57 (7.44)

1.33 (3.49)
0.71 (2.08)
3.64 (3.30)
2.81 (2.70)

0.28 (0.98)
0.11 (0.57)
3.12 (1.61)
2.40 (1.94)

21.6 (31.1)–
20.1 (28.0)–

17.0 (24.6)–
16.4 (24.0)–

1.30 (1.90)–
1.65 (1.53)–

Summary of the localization error (mm). The localization errors averaged over 500 randomly placed dipole and distributed sources in the MRI-based corticalsurface model are listed in millimeters. The standard deviations are listed in parentheses. Noiseless and noisy conditions, where MEG signals contain Gaussian
sensor noise with N/S of 0.1, were tested under three different situations depending on the available information: no other data than MEG (NO), MEG with
MRI (MRI), and MEG with MRI and fMRI (MRI + fMRI). The estimation methods are the VB method (VB), the minimum norm estimate (MN), the minimum
norm estimate with smoothness constraints (MNC), and the Wiener filter method (Wiener). The noisy conditions were not tested for MN, MNC, and Wiener
filter methods as indicated by – in the table.

M. Sato et al. / NeuroImage 23 (2004) 806–826

815

Fig. 4. Localization error for a single dipole source. (A–C) Scatter plots of localization errors for single dipole source as a function of the lead field norm. (A)The VB method for MEG with no other data (NO). (B) That for MEG with MRI (MRI). (C) That for MEG with both MRI and fMRI (MRI + fMRI). (D) Pointspread curves of the current estimates around the source dipole were calculated by averaging over the 500 dipole sources as a function of the Euclidean distancefrom the source dipole position. Estimation methods are VB estimate with NO other information (solid red line), VB estimate with MRI (solid green line), VB
estimate with MRI and fMRI (solid blue line), minimum norm estimate (MN) with MRI (dashed cyan line), minimum norm estimate with smoothness
constraints (MNC) with MRI (dashed magenta line), and Wiener filter estimate with MRI and fMRI (dashed black line).

ization errors for dipole sources with small lead field norms werelarger than those with large lead field norms in all situations (Figs.
4A–C) because the estimation of the sources with small lead fieldnorms was difficult. The scatter plots also show that introduction of
MRI and fMRI information significantly decreases localization
errors, especially for large lead field norms. In addition, thepositions of dipole sources were perfectly recovered for large leadfield norms when both MRI and fMRI information were imposed.
Table 2 also shows that the precision of the VB estimates with
MRI information is far better than those of the MN and MNCestimates both for single and distributed sources. The precision of
the VB estimates with MRI and fMRI information is comparablewith that of the Wiener filter, that is, that of the VB estimates isslightly better than that of the Wiener filter for a single dipole
source, but the converse is true for a single distributed source(Table 2). It should be noted that the localization error of less than
7 mm for the distributed source implies that the source position iscorrectly estimated inside the distributed source, since the radius of
the distributed source is 7 mm. Therefore, the differences of
localization errors less than 7 mm do not have much meaning forthe distributed source, unless the detailed spatial profile of theestimated currents are examined.
The spatial profile of the estimated currents for single dipole
source was assessed as the point-spread curve. It was defined by
the estimated current profile around the source dipole and averaged

across the 500 dipole sources as a function of the Euclidean
distance from the source dipole position. If the current estimation isperfect, the value of the point-spread curve at the origin, which isequal to the estimation gain, should be one and that for positive
distance should be zero. Fig. 4D shows the point-spread curve of
the VB estimates for three situations: (1) MEG with no other data
(NO), (2) MEG with MRI (MRI), and (3) MEG with both MRI andfMRI (MRI + fMRI). The point-spread curve for NO case (red line
in Fig. 4D) was rather broadly spread with FWHM = 5 mm. The
estimation gain (height of the point-spread curve at the origin) was
0.68 in this case. The point-spread curves for MRI and MRI +fMRI cases (green and blue lines in Fig. 4D) were sharply tuned at
the source position with FWHM = 1 mm, and the estimation gains
were 0.48 and 0.81, respectively. These results show that the
introduction of MRI and fMRI information significantly increases
the spatial resolution of the VB method.
The point-spread curves of the MN and the MNC methods(dashed cyan and magenta lines in Fig. 4D) were almost flat andnearly overlapped with the x-axis, which indicated very poor
resolution of these methods. The point-spread curve of the Wiener
filter estimate (dashed black line in Fig. 4D) was worse than that of
the VB estimate for MRI + fMRI (FWHM = 4 mm; estimation
gain, 0.27).
These results show that the precision of the source localization
for the VB method is comparable with that for the Wiener filter

816

M. Sato et al. / NeuroImage 23 (2004) 806–826

method when structural MRI and correct fMRI data are available.
On the other hand, the estimated current profile of the Wiener filter
method have much larger spread than that of the VB method for asingle dipole source. Both measures are much better for the VB
method than for the MN and the MNC methods when only
structural MRI data is available.
VB and linear inverse filter estimates for simulated corticalresponses
Comparison of the VB performance for single and distributed
current sources among the three situations indicates that the MRI
and fMRI information effectively improves the spatial resolution as
well as the reliability of the VB estimates. However, mismatches
could occur between MEG and fMRI information for several
reasons, including noise artifacts and the slow hemodynamic
response of fMRI. We investigated how these mismatches affect
the VB and the Wiener filter estimates for MEGs by conducting
simulations reproducing realistic visual cortical responses.
We identified cortical areas in a subject’s (MA) brain including
V1v, V1d, V2v, V2d, V3, and VP (Fig. 5A). This was based on thecortical retinotopy determined by fMRI using a rotating pole andexpanding/shrinking annulus stimuli (Engel et al., 1994; Linden etal., 1999; Sereno et al., 1995). We simulated sequential activities
evoked in visual cortical areas by lower-quadrant visual fieldstimulation. Sinusoidal currents uniformly distributed in vertex
points within a circular area (radius, 5 mm) of V1v were assumed
and those in V2v and VP were activated with 50- and 100-ms
delays to those in V1v, respectively (Figs. 5C and D). White
Gaussian sensor noise with N/S = 0.1 was added to the generated
MEG signals.
To evaluate spatio-temporal patterns of the estimated currents,
we defined local V1v, V2v, and VP currents as the sum of theestimated currents within a 5-mm radius of the local peak current
positions in the V1v, V2v, and VP areas, respectively. The
precision of the estimation was assessed as the estimation gains
for the V1v, V2v, and VP local currents at 50, 100, and 150 ms,
respectively, where each source current became maximum. The
normalized root mean squared errors (RMSEs) were also calculated for the temporal sequences of the local V1v, V2v, and VP
currents. The estimation gain should be one and RMSE should bezero for perfect estimation.
MEG with no fMRI information
We first investigated the performance of the VB estimate for thesimulated visual cortical responses, in the situation where MEGs
were supplemented with only structural MRI information. The VB
estimate was conducted according to the method of high-resolution
estimation for MEG combined with MRI. Time-lapse current maps
using pseudocolor representation revealed that the VB estimate

Fig. 5. Simulated visual responses. (A) The cortex model used for the simulated visual responses. The visual area is depicted on the left hemisphere. (B) The
visual areas are shown by using an inflated map of the cortex. Visual areas V1v, V1d, V2v, V2d, VP, and V3 are shown by different colors. The generated
source currents with a 5-mm radius in the V1v, V2v, and VP areas are shown in yellow. (C) Time course plot of the V1v (red line), V2v (green line), and VP(blue line) currents. Total current strength of each source is shown. Abscissa denotes the time in ms. (D) Time-lapse map of source currents at t = 50, 100, and
150 ms are shown on the inflated map. (E) Simulated fMRI activity given as the sum of the source signal amplitude and 0.1 N/S white Gaussian noise
convolved with a Gaussian smoothing filter (FWHM = 6 mm). (F) False-positive fMRI signal including V1d, V2d, and V3 activities in addition to the true
fMRI signals in V1v, V2v, and VP.

M. Sato et al. / NeuroImage 23 (2004) 806–826

817

Fig. 6. Estimated current with no fMRI information. (A and C) Time-lapse current maps estimated by the VB method (A) and the MNC method (C) with nofMRI information are plotted using pseudocolor representation at t = 50, 100, and 150 ms. (B and D) Time course plot of the V1v (red line), V2v (green line),
and VP (blue line) currents estimated by the VB method (B) and the MNC method (D) with no fMRI information are plotted. Each current is defined as acurrent sum within a 5-mm radius of the local peak current positions in the V1v, V2v, and VP areas.

The performances of MNC were very poor. The time-lapse
maps revealed that the estimated currents broadly spread out fromthe areas in V1v, V2v, and VP for which current sources wereassumed (Fig. 6C). The time plots showed much smaller peaks
for V1v, V2v, and VP than those for the VB estimates (Fig. 6D).
In fact, the estimation gains for MNC were very small and
RMSEs were large as can be seen in Table 3. In addition, thetime course of V2v currents was more broadened than that of thetrue source currents (compare green trace in Fig. 6D with those in
Fig. 5C). These findings indicate that the MNC estimate is not
feasible for the MEG of realistic visual cortical responses. The
performance of the VB method is much better than that of thevirtually recovered V1v and VP responses but left V2v responses
undetectable (Fig. 6A). Correspondingly, the time plots for local
currents indicated a rather fine match between the estimated andtrue source currents assumed in V1v and VP (compare red and blue
traces in Fig. 6B with those in Fig. 5C). However, the estimatedV2v current was far smaller than the true source current assumed inV2v (compare green trace in Fig. 6B with those in Fig. 5C). Thesefindings can be confirmed by the estimation gains and RMSEs for
VB with the no fMRI situation in Table 3. The current estimates forV2v were difficult, since the lead field for V2v was much smaller(roughly 0.1) compared with those for V1v and VP currents, andV2v has a very close distance (less than 2 mm) to VP.
Table 3
Accuracy of spatio-temporal current estimates
MethodVB

MNC

WeinerInformation

No fMRI

Correct fMRI

Incorrect fMRI

No fMRI

Correct fMRI

Incorrect fMRI

0.85
0.19
0.70
0.14
0.82
0.34

0.95
0.86
0.87
0.06
0.21
0.12

0.96
0.86
0.88
0.05
0.20
0.12

0.27
0.15
0.24
0.72
0.84
0.77

0.83
0.58
0.95
0.20
0.51
0.34

0.93
0.24
1.06
0.35
0.59
0.55

Estimation gain

RMSE

V1v
V2v
VP
V1v
V2v
VP

Accuracy of spatio-temporal current estimates. The performance of the VB method was investigated using simulated visual cortical responses under threesituations: MEG with MRI data (No fMRI), MEG with MRI and correct fMRI data (Correct fMRI), and MEG with MRI and incorrect fMRI data (Incorrect
fMRI). The precision of the estimation was assessed as the estimation gains for the V1v, V2v, and VP local currents at 50, 100, and 150 ms, respectively. The
normalized root mean squared errors (RMSEs) for the temporal sequences of the local V1v, V2v, and VP currents were also listed in the table. The MNC
method with no fMRI information and the Wiener filter method with correct fMRI and incorrect fMRI information were also tested.

818

M. Sato et al. / NeuroImage 23 (2004) 806–826

MNC estimate, but the VB estimate with no fMRI information isstill insufficient to recover the fine dynamic structure of visual
responses.
MEG combined with correct fMRI information
Next, we examined whether the correct fMRI information
improves the VB and the linear filter performance for the same
simulated cortical responses. fMRI signals were assumed to be thesum of the source signal and 0.1 N/S Gaussian white noise
convolved with a Gaussian smoothing filter for fMRI imaging(FWHM = 6 mm) (Fig. 5E). The VB estimate was conducted in thesame way as for the situation of MEG combined with MRI and
fMRI. The time-lapse maps and time plots indicated that the VB
estimates nearly correctly recovered the current activities in V1v,
V2v, and VP (compare Figs. 7A,B with Figs. 5D,C). The estimatedcurrents were slightly smaller than the true source currents
assumed in these cortical areas. The estimation gains for V1v,
V2v, and VP currents were 0.95, 0.86, and 0.87, and RMSEs for
V1v, V2v, and VP currents were 0.06, 0.21, and 0.12, respectively(Table 3).
The performances of the VB estimates were compared withthose of the Wiener filter (Dale et al., 2000; Kajihara et al., 2004;
Phillips et al., 2002).
were conducted as
 The Wiener filter estimates2
Eq. (3) with 1
0 nm ¼ CVn Vm expððdnm =rÞ Þ (FWHM = 8 mm)
for the same MEGs and fMRI signals. The coefficient C wasdetermined so as to maximize the correlation between the fMRI
signals and the temporal average of the estimated currents(Kajihara et al., 2004). The Wiener filter estimates recovered
activities in V1v, V2v, and VP (compare Fig. 7C with Fig. 5D), but
V2v current was rather smaller than the true source current
assumed in V2v. The estimation gains for V1v, V2v, and VP
currents were 0.83, 0.58, and 0.95, and RMSEs for V1v, V2v, andVP currents were 0.20, 0.51, and 0.34, respectively (Table 3). In
addition, the time plots show a tendency of some oscillations intime courses of V1v, V2v, and VP currents (compare Fig. 7D with
Fig. 5C). Consequently, correct fMRI information almost perfectly
improved the performance of the VB estimate, but this improvement was still insufficient for the linear inverse filter.
MEG combined with incorrect fMRI information
We further examined how false-positive fMRI signals affect the
VB estimate. The VB estimate was conducted for the same
simulated cortical responses as those for the correct fMRI situation
using incorrect fMRI information. fMRI signals were assumed tocontain false-positive signals in V1d, V2d, and V3 as strong as thecorrect fMRI signals in V1v, V2v, and VP in addition to the correctones in these areas (Fig. 5F). The performance of the VB estimate
remained almost unaffected by false-positive fMRI signals as canbe seen from the estimation gains and RMSEs in Table 3. The
time-lapse map and time course plot of the VB estimates were
almost the same as those for the correct fMRI situation (compare
Figs. 8A,B with Figs. 7A,B). No false currents were estimated inV1d, V2d, and V3, where false-positive fMRI signals wereassumed (compare Fig. 8A with Fig. 5D).

Fig. 7. Estimated current with correct fMRI information. (A and C) Time-lapse current maps estimated by the VB method (A) and the Wiener filter method (C)with correct fMRI information are plotted at t = 50, 100, and 150 ms. (B and D) Time course plot of the V1v (red line), V2v (green line), and VP (blue line)
currents estimated by the VB method (B) and the Wiener filter method (D) with correct fMRI information.

M. Sato et al. / NeuroImage 23 (2004) 806–826

819

Fig. 8. Estimated current with false-positive fMRI information. (A and C) Time-lapse current maps estimated by the VB method (A) and the Wiener filter
method (C) with false-positive fMRI information are plotted at t = 50, 100, and 150 ms. (B and D) Time course plot of the V1v (red line), V2v (green line), andVP (blue line) currents estimated by the VB method (B) and the Wiener filter method (D) with false-positive fMRI information.

In contrast, the Wiener filter estimate was seriously affectedby false-positive fMRI signals. The time-lapse map indicated that
substantial currents were estimated around the area of V1d, V2d,
and V3, where false-positive signals were assumed (compare Fig.
8C with Fig. 5F). Correspondingly, the time course of V1v, V2v,
and VP currents became much more oscillatory than those for thecorrect fMRI situation (compare Fig. 8D with Fig. 7D) due toincreased ambiguity caused by false-positive signals. As a result,
RMSEs became larger than those for the correct fMRI case
(Table 3).
These findings show that the Wiener filter estimates are
sensitive to false-positive fMRI signals but not the VB estimates.
Separation of simultaneously active currents
The spatial resolution of the VB and the Wiener filter estimates
were also examined. The resolution matrix defined by Eq. (18) (seeAppendix A) was calculated using the optimized inverse filters for
MEGs of the simulated visual responses with the correct fMRI.
The nth row of the resolution matrix represents the point-spread
profiles of the current estimates for a single dipole source at the nthvertex point. The spatial profiles of the current estimates for
simultaneously active V1v, V2v, and VP currents (Fig. 9A) were
calculated by summing the rows of the resolution matrix
corresponding to the active vertex points of the V1v, V2v, and
VP. Profile maps indicate that both the VB and the Wiener filter
methods nearly correctly recover the source current profiles in V1v,
V2v, and VP (compare Figs. 9B and C with Fig. 9A).
Likewise, the spatial profiles of the current estimates for
simultaneously active V1v, V2v, and VP (Fig. 9A) were determinedaccording to the resolution matrix of the VB inverse filter optimized
for the MEGs combined with the false-positive fMRI (Fig. 9D).
Profile maps indicate that the VB estimates for simultaneously
active V1v, V2v, and VP correctly recover the source current
profiles (compare Figs. 9E with A) under the false-positive fMRI
information. Conversely, the current estimates based on the
resolution matrix of the Wiener filters optimized for the same setof MEG and incorrect fMRI as those for the VB estimates wererather seriously affected by the false-positive fMRI signals(compare Figs. 9F with D).
MEG combined with insufficient fMRI information
fMRI also may not provide correct information for travelingcortical responses. We studied the performance of the VB estimate
assuming that a uniform circular current source (radius, 5 mm)stayed in V1 for the 200-ms period of time and traveled (distance,
50 mm; velocity, 0.5 m/s) to the other region in V1 (Fig. 10A).
This may correspond to V1 responses, where a spot near the
fixation point moves toward the upper-right direction.
fMRI signals were given as a sum of the temporal average ofthe source currents and 0.1 N/S Gaussian white noise convolved
with a Gaussian smoothing filter (FWHM = 6 mm) (Fig. 10A).
Time-lapse maps indicate that the VB estimates for MEGs
generated by the traveling V1 responses almost recovered thespatio-temporal patterns of the traveling source currents (compare
Figs. 10C with B). However, the estimates were slightly deformed
and weaker than the true traveling sources. On the contrary, theWiener filter estimates almost failed to recover the spatio-temporal
patters of the traveling source currents (compare Figs. 10D with B).

820

M. Sato et al. / NeuroImage 23 (2004) 806–826

Fig. 9. Current profiles for simultaneously active V1v, V2v, and VP currents. Spatial profiles of the current estimates for simultaneously active V1v, V2v, andVP currents are shown on the inflated map of the visual area. (A) Simultaneously active currents in V1v, V2v, and VP. (B and C) Current profiles estimated by
the VB method (B) and the Wiener filter method (C) with correct fMRI information. (D) False-positive fMRI signals in V1d, V2d, and V3 in addition to thecorrect fMRI signals in V1v, V2v, and VP. (E and F) Current profiles estimated by the VB method (E) and the Wiener filter method (F) with false-positive fMRI
information.

A static source, which corresponded to the large fMRI signal, wasrecovered most of the time. These results show the advantage ofthe VB estimate over the Wiener filter estimate for spatiotemporally changeable brain activities.

Discussion
We proposed a novel hierarchical Bayesian method to resolve
the MEG inverse problem, which gives a unified framework forthree possible situations: (1) MEG with no other data (NO), (2)MEG with structural MRI data on cortical surfaces (MRI), and(3) MEG with both structural MRI and fMRI data (MRI + fMRI).
Our method estimates the current variance from the time
sequence of the MEG data by introducing a hierarchical prior.
The fMRI information can be imposed as prior information on
the variance distribution rather than the variance itself so that itgives a soft constraint on the variance. The spatial smoothness
constraint can also be incorporated as a hierarchical prior. The
smoothness constraint reduces the effective degree of freedom inthe estimation model, and thus it reduces the ill-posedness of theinverse problem.
Because of the hierarchical prior, the estimation problem
becomes nonlinear, so the approximate posterior distribution of
the source current is calculated by using the Variational Bayesian(VB) method. An important difference between the VB method
and the MAP estimate used for the linear inverse methods is that
the VB estimate calculates the posterior probability distribution of
the source current while the MAP estimate calculates the optimalcurrent value that maximizes the posterior distribution. The VB
method takes account of the spread of the posterior distribution inestimating the current variance, so robust estimation for thecurrent variance can be achieved. In the MEG inverse problem,
the number of unknown parameters is much larger than thenumber of MEG sensors M. The parameter estimation of such asystem often suffers from overfitting to the observed data. The
VB method can avoid such overfitting because the model
complexity penalty term in the free energy (see Eq. (22) inAppendix B.1) freezes unnecessary degrees of freedom bysuppressing the estimation gain for small variance points. The
advantages of the VB estimate over the MAP estimate were also
shown by Penny et al. (2003) for fMRI data analysis. It is also
noted that the current variances are estimated from T time samplesof the MEG data. Therefore, the effective number of data forvariance estimation becomes (T  M). This improves the accuracyof the variance estimation.
Since the VB method solves the nonlinear estimation problem,
the VB solution converges into a local maximum of the free energy.
Therefore, the converged solution depends on the initial value of theestimated current variance. Fortunately, we have reasonablecandidates for the initial value of the current variance. When there
is no fMRI information, the initial value of the current variance is
set to the same value for all dipole currents. This setting corresponds
to the minimum norm solution, that is, the VB current estimate at
the first iteration is given by the minimum norm solution. The VB
method, then, improves the accuracy of the current estimate by
estimating the current variance from the time sequence of the MEG
data. A single hyperparameter a 0 controls the initial variance as
1
2
(a 0n
= a 0 /G ave), whichby the average of the lead
 P is normalizedPMN
2
2
field norm Gave ¼
G
n¼1
m ¼ 1 m;n =N. In all the simulationsof this paper, a 0 = 10 was used. When there is fMRI information, it
is also used as the initial value of the current variance (seeM. Sato et al. / NeuroImage 23 (2004) 806–826

821

Fig. 10. Traveling wave responses in visual area. (A) Traveling wave responses schematically shown by an arrow, assuming that a uniform circular current
source (radius, 5 mm) stayed in V1 for a 200-ms period of time and traveled (distance, 50 mm; velocity, 0.5 m/s) to the other region in V1. Simulated fMRI
activity given as the sum of the temporal average of the source currents and 0.1 N/S Gaussian white noise convolved with a Gaussian smoothing filter(FWHM = 6 mm) is also shown using pseudocolor representation. (B) Time-lapse map of the source currents at t = 200, 220,..., 300 ms. (C and D) The timelapse map of the currents estimated by the VB method (C) and the Wiener filter method (D).

Variational Bayesian method). Therefore, the VB current estimate
at the first iteration is given by the Wiener filter solution in thiscase. The initial variance is controlled by two hyperparameters
a min and a max as defined by (32) in the Appendix B.4. They wereset as a min = 103 and a max = 10 in all the simulations of this paper.
The above hyperparameter values were determined by the
maximum free energy criterion applied to one of the simulated
MEG data (see the procedure explained in the Appendix B.4). The
range of the searched values were 1 V a max V 105 and 1 V a max /a min V 105. The obtained solutions were not sensitive to smallchanges of the hyperparameter values.
There are other factors, which affect the estimation accuracy.
In the simulations, we used the same cortical surface model bothfor the forward magnetic field calculation and for the inverse
source current estimation. However, in the analysis of realexperimental data, there may be inaccuracy in the cortical surface
estimation. Since the cortical surfaces have highly folded
structures, small errors in the cortical position may cause large
errors in the normal direction and then cause large lead fielderrors. Therefore, high accuracy of the cortical surface estimationis required when the dipole direction is imposed as the structuralprior information. The effective radius (FWHM) of the smoothness (Gaussian) filter also affects the estimation performance as
explained in the MEG with Structural MRI Information section.
Too small filter radius disables the smoothness constraint while
too large filter radius smears the current estimation. An
appropriate value may depend on the resolution and the accuracyof the extracted cortical surfaces.
As pointed out previously, fMRI activity may not precisely
correspond to the MEG source current activity. In particular, thefMRI activity corresponds to an average of several thousands of
MEG time series data because of the time resolution difference
between MEG and fMRI. When higher brain functions are
examined, several different processes may occur within several
seconds, and these activities are averaged out in the fMRI data.
Consequently, the fMRI active areas may include inactive current
areas in the MEG time slice data, and the fMRI activity for theactive current may be smeared by the temporal averaging. Whenthe fMRI activity is used to determine the source current variancesuch as for the Wiener filter method, the fMRI information may
mislead the source current estimation in these situations. To deal
with these situations, we divide the entire MEG time series into

822

M. Sato et al. / NeuroImage 23 (2004) 806–826

several short periods in which the brain activity changes
moderately, and the VB method is applied to the MEG data foreach period. Although the VB estimate starts from the Wiener filter
solution, the current variance estimation using the short-period
MEG data may improve the accuracy of the source current
estimation.
An advantage of the VB method is the model selection
capability. If there is no structural MRI information, we consider
a variable-radius surface model for current estimates and the mostprobable surface is determined by the surface with the maximum
free energy. When there was no noise, the free energy was
sharply tuned at the peak position (Fig. 3B). As the noise
amplitude increased, the peak profile of the free energy curve
became broader as shown in Fig. 3C. However, the peak almost
stayed at the same position and there was no local peak aroundthe maximum peak position in the presence of noise. Therefore,
the source depth was safely determined by the free energy
maximum. We also examined the standard single dipole fit for the
single dipole source problem shown in Fig. 3. The reconstruction
error (Fig. 3D) of the dipole model under the same noise
condition showed the similar (reversed) tuning curve as the freeenergy (Fig. 3C).
In the simulations, the VB method iterated the current
estimation until the free energy converged. Therefore, theprocessing time would be an important concern in practice.
The VB estimate with the hierarchical prior Eq. (9) (50 MEG
time samples and 2494 vertex points) required about 1 min on a
XEON (2.2 GHz) processor, while that with the smoothness prior
Eq. (11) (200 MEG time samples and 819 vertex points) required 6
min.

patterns of brain activities due to the slow hemodynamic response,
which would commonly be the case in many experimental
protocols.
In this article, we have concentrated on the MEG inverse
problem with fMRI information. It should be mentioned that our
method can also be applied to electroencephalography (EEG) andthat fMRI information can be replaced by positron emission
tomography (PET) data.

Acknowledgments
We would like to thank Dr. Daniel Callan for careful reading of
the manuscript. This research was supported in part by the National
Institute of Information and Communications Technology (NICT).

Appendix A. Resolution matrix and estimation gain
We define the resolution matrix and estimation gain in thisAppendix and derive the inequality on the estimation Eq. (8). In theabsence of noise, the magnetic field B produced by the current J isgiven by B = G d J. The estimated current Ĵ is written asJˆ ¼ L d B ¼ L d G d Jð17Þusing the linear inverse method (3). The resolution matrix R andthe estimation gain G (Grave de Peralta Menendez and Gonzalez
Andino, 1998) are defined by
R ¼ L d G; G ¼ Rn;n jn ¼ 1: Ng:
ð18Þ

The advantages of the VB method over the linear filters wereconfirmed in simulation studies using a single source and realistic
visual responses. The simulation results can be summarized asfollows.

The (n,k) component of the resolution matrix R n,k represents
the nth estimated current when the unit current dipole is applied forthe kth current, that is, J k = 1 and J l = 0(l p k). Therefore, theresolution matrix R should be the identity matrix in the current
space for perfect reconstruction. Specifically, the estimation gain
Gn (n = 1:N) should be one. However, the following simple
calculation shows the relation:

(1)

N
X

Summary(2)(3)

When the same information is given, the VB method achieves
better performance than the conventional linear inverse filter
methods.
When the correct information is given, the performance of the
VB method improves as more information becomes available.
When inaccurate fMRI information is given, performancedegradation of the current estimate is much less for the VB
method than for the Wiener filter method.

Consequently, the VB method provides a unified framework to
resolve the MEG inverse problem for three possible situations:
MEGs combined with no other information, with only MRI, andwith both MRI and fMRI. The first version may be useful for
analysis of MEGs whose current sources are in deep brain
structures such as the thalamus and basal ganglia because reliable
estimation of dipole locations and orientations based on MRI isdifficult. The second version may also be useful for situations inwhich fMRI measurement using the same experimental conditions
as those for MEG is difficult. Finally, the third version provides themost reliable solution of the MEG inverse problem. It is importantthat the VB method may correctly resolve the inverse problem
even if fMRI signals fail to convey the precise spatio-temporal1
Gn ¼ TrR ¼ M  b1 Tr G d 1IM
0 d GV þ b

1

V M:

n¼ 1

ð19Þ
This equation holds because both b1 I M and G d 1
0 d GV are
positive semidefinite matrices. The inequality (19) was derived forthe first time in this paper, as far as we know.

Appendix B. Variational Bayesian method for hierarchicalprior
In the following, we explain the VB method for the hierarchicalprior with the smoothness constraint, (11) and (12). The
hierarchical prior, (9) and (10), is a special case of the hierarchicalprior with the smoothness constraint, and it is obtained by setting
W u 0.
B.1. Variational Bayesian method
First, we explain the main ideas of the VB method using thehierarchical prior, (9) and (10), with a single MEG data set B. The

M. Sato et al. / NeuroImage 23 (2004) 806–826

objective of the Bayesian estimation is to calculate the posterior
probability distribution of J for the observed data B:
ZPðJ jBÞ ¼
dAPðJ ;AjBÞ;PðJ ;AjBÞ ¼

PðJ ;A; BÞ;PðBÞPðBÞ ¼

The calculation of the marginal likelihood P(B) cannot be done
analytically. In the VB method, the calculation of the joint posterior
P(J, AjB) is reformulated as the maximization problem of the freeenergy. The free energy for a trial distribution Q(J,A) is defined by"#ZPðJ ;A;BÞ
F ðQ Þ ¼
dJ dAQðJ ;AÞlog
Q ðJ ; A Þ
ð20Þ

Eq. (20) implies that the maximization of the free energy F(Q)is equivalent to the minimization of the Kullback–Leibler distance(KL-distance) defined byZ
KL½ QðJ;AÞtPðJ;AjBÞ ¼ dJ dAQðJ;AÞlogðQðJ;AÞ=PðJ;AjBÞÞ:This measures the difference between the true joint posterior
P(J,AjB) and the trial distribution QðJ ; AÞ. Since the KL-distancereaches its minimum at zero when the two distributions coincide,the joint posterior can be obtained by maximizing the free energy
F ðQÞ with respect to the trial distribution Q. In addition, themaximum free energy gives the log-marginal likelihood log( P(B)).
The optimization problem can be solved using a factorization
approximation restricting the solution space4 (Attias, 1999; Sato,
2001):
QðJ; AÞ ¼ QJ ðJÞQA ðAÞ:ð23ÞIn the second step (a-step), the free energy F(Q) is maximized
with respect to Qa while QJ is fixed. The solution is given bydJ dAPðJ ; A;BÞ:¼ logð PðBÞÞ  KL½QðJ ; AÞtPðJ ; AjBÞ :degree of freedom that can be well specified from the observed
data. Therefore, (negative sign of) the free energy can beconsidered a regularized error function with a model complexity
penalty term. The maximum free energy is obtained by alternately
maximizing the free energy with respect to QJ and Qa . In the firststep (J-step), the free energy F(Q) is maximized with respect to QJ
while Qa is fixed. The solution is given by
QJ ðJ Þ~exp½hlogPðJ ;A;BÞia :PðJ ;A;BÞ ¼ PðBjJ ÞP 0 ðJ jAÞP 0 ðAÞ;
Z

823

ð21Þ

Under the factorization assumption (21), the free energy can be
written as
F ðQÞ ¼ hhlogPðJ ; A; BÞiJ ia  hlogQJ ðJ ÞiJ  hlogQa ðAÞia
¼ hlogPðBjJ ÞiJ  KL½QJ ðJ ÞQa ðAÞtP0 ðJ jAÞP0 ðAÞ ;ð22Þ
where hd iJ and hd ia represent the expectation values with respect
to QJ ðJ Þ and Qa (A), respectively. The first term in the second Eq.
(22) corresponds to the negative sign of the expected reconstruction error. The second term (KL-distance) measures the difference
between the prior and the posterior and corresponds to the effective

4
The VB method with this factorization assumption gives moreaccurate approximation than the mean field approximation (MFA), which is
widely used in the field of statistical physics (Parisi, 1988), because MFA
has a stronger factorization assumption, that is, that all of the variable
components can be factorized.

Qa ðAÞ~exp ½ logPðJ ; A; BÞ J :ð24Þ

The above J- and a-steps are repeated until the free energy
converges.
B.2. Hierarchical prior with smoothness constraints
We assume the following probabilistic model for the source
currents:1PðBjJ ; bÞ~exp  b ðB  G d J ÞVd G d ðB  G d J Þ ;
2
ð25Þ
1

where ðbG Þ denotes the covariance matrix of the sensor noise.
1
1
¼ M,
G is the normalized covariance matrix satisfying Tr G
1and b is the average noise variance.
The direct application of the VB method to Eqs. (11) and (12)cannot be solved in a closed form because the inverse covariancematrix A is a nonlinear function of A and L. To get the closed
form solution, we rewrite the smoothness prior by introducing an
auxiliary internal variable Z(t) = {Z n (t)jn = 1:N}:
P0 ðJ 1:T jZ 1:T ; A; bÞ
"
#
T
1 X
~exp  bðJ ðtÞ  W d Z ðt ÞÞV d A d ðJ ðtÞ  W d Z ðtÞÞ ;
2 t¼1
"
#
T
1 X
ZVðtÞ d K d Z ðt Þ :ð26Þ
P0 ðZ 1:T jL; bÞ~exp  b
2 t¼1
The inverse noise variance b is also introduced in Eq. (26) as acoefficient of A and K to get the closed form solution. Therefore,
the definitions of A and L in this Appendix are different from thosein the text by a factor b. Integration over Z in the prior (26) can beaccomplished as
P0 ðJ 1:T jA; b; LÞ ¼

Z

dZ 1:T P0 ðJ 1:T jZ 1:T ; A; bÞP0 ðZ 1:T jL; bÞ
"
#
T
1 X
J Vðt Þ d A d J ðt Þ :~exp  b
2 t¼1

Namely, the prior (26) is equivalent to the smoothness prior(11) and (12). Parameters A and L are estimated by introducing the
ARD hierarchical prior,N

P0 ðAÞ ¼ P Cðan jā0n ; c0na Þ;
n¼1NP0 ðLÞ ¼ P Cðkn jk̄0n ; c0nk Þ:
n¼1

824

M. Sato et al. / NeuroImage 23 (2004) 806–826

The inverse noise variance parameter b is also estimated byassuming the non-informative prior:

The expectation values of the current J(t) and the noise
variance b 1 are estimated using the inverse filter (28).

P0 ðbÞ ¼ b1 :d Bðt Þ;J ðt Þ ¼ L 1
A

B.3. Free energy

Z ðt Þ ¼ K1 d W V d A d J ðt Þ;

To calculate the joint posterior P(J 1:T, Z 1:T, A, b, LjB 1:T ), thefree energy for the trial distribution QðJ 1:T ; Z 1:T ; A; b; LÞ isdefined asZF ðQ Þ ¼
dJ 1:T dZ 1:T dbdAdL QðJ 1:T ; Z 1:T ; A; b; LÞ

;b ¼


 log


PðJ 1:T ; Z 1:T ; A; b; L; B1:T Þ;QðJ 1:T ; Z 1:T ; A; b; LÞ

1
NT;
2
"
T
1 X
1
;b b ¼ðBðtÞ  G d J ðt ÞÞV d G d ðBðt Þ  G d J ðtÞÞ2 t¼1
þ ðJ ðt Þ  W d Z ðt ÞÞV d A d ðJ ðt Þ  W d Z ðt ÞÞ#þ ZVðtÞ d A d Z ðtÞ :

ð27Þ

ð29Þ

In the a-step (Fig. 1E), the expectation values of the varianceparameters are estimated as

PðJ 1:T ; Z 1:T ; A; b; L; B1:T Þ ¼ PðB1:T jJ 1:T ; bÞ
 P0 ðJ 1:T jZ 1:T ; A; bÞP0 ðZ 1:T jL; bÞcna ¼ c0na þT;2

cna ¼ c0na þ P0 ðbÞP0 ðAÞP0 ðLÞ:1
cna a1n ¼ c0na an þ

T
;
2

 i
T 1 han 1  A1 d GV d 1
B d G n;n2

In this case, we assume the factorizationT
1 X
ðJ ðtÞ  W d Z ðt ÞÞ2n ;
þ b
2 t¼1

QðJ 1:T ; Z 1:T ; A; b; LÞ ¼ QJ ðJ 1:T ; Z 1:T ; bÞQa ðA; LÞ:The maximum free energy can be obtained by alternately
maximizing the free energy F ðQÞ with respect to QJ and Qa . In thefirst step (J-step), the free energy F ðQÞ is maximized with respect
to QJ . The maximum solution of QJ is a Normal-Gamma
distribution:

1
cnk k1n ¼ c0nk k0 þh i
 1  K1 d W V d GV d 1
d
G
d
W
B
n;n
T
1 X
Z 2 ðt Þ:þ b
2 t¼1 n

QJ ðJ 1:T ; Z 1:T ; bÞ ¼ QJ jb ðJ 1:T ; Z 1:T jbÞQb ðbÞ:The conditional posterior distribution QJ jb ðJ 1:T ; Z 1:T jbÞ is theNormal distribution for (J 1:T, Z 1:T ) with the noise varianceparameter b 1. The posterior distribution for b, Qb (b), is theGamma distribution. The expectation values of these parameters are calculated in the J-step. In the second step (a-step),
the free energy F ðQÞ is maximized with respect to Qa . The
maximum solution of Qa is the Gamma distribution for A and
L. The expectation values of A and L are calculated in thea-step.
B.4. VB algorithm
The VB algorithm for the hierarchical prior with the smoothness constraint, (11) and (12), is summarized here (Fig. 1E). All ofthe quantities in the following equations represent the expectation
value with respect to the posterior distribution. In the J-step (Fig.
1E), the inverse filter L(A1) is calculated using the estimatedcovariance matrix A1 in the previous iteration:
J ¼ GV d G d G þ A ;

ð30ÞThe covariance matrix 1
A is then recalculated as
1
þ W d K1 d W V:
1
A ¼ A

ð31Þ

The above J- and a-steps are repeated until the free energy
converges. In the simulations, we further imposed a constraint L =jA, which implies that the correlation in the current activities is
proportional to the current variance, to reduce the number of
unknown parameters. Under this constraint, the expectation valuesof A, j, and L are calculated in the a-step by applying thefollowing update rule after the update of (30):

1
1
a1
n ¼ cna an þ cnk kn j =ðcna þ cnk Þ;

j1 ¼

N
1 X
an =kn ;N n¼1

kn ¼ jan :1
B ¼ G d 1
A d GV þ G ;


1
1
L 1
¼ 1
A
J d GV d G ¼ A d GV d B :

T 1
k
2 n

fMRI information was imposed as
ð28Þ

 2
2
a1
0n ¼ amax Vn þ amin =Gave

ð32Þ

M. Sato et al. / NeuroImage 23 (2004) 806–826

where V n denotes the fMRI signal for the nth current. It is
normalized P byPthe average of the lead field normNM
2G2ave ¼
n¼1
m ¼ 1 Gm;n =N. The hyperparameters a max
and a min control the maximum and minimum of the priormean variance, respectively. The optimal value of the hyperparameters (a max, a min) can be determined by maximizing thefree energy with respect to (a max, a min). (1) The free energy ismaximized for a given MEG signal and a sampled value of(a max, a min). (2) Then, the optimal value of (a max, a min) thatmaximizes the free energies obtained for the sampled values in
the first step is determined.
B.5. Posterior distribution for the source current
After the free energy converges, the posterior distribution forthe current J 1:T is obtained by integrating the internal variable Z 1:T
and the inverse noise variance parameter b: "jJ jT =2 C cb þ NT =2
b¯

 1þ
QðJ 1:T Þ ¼
N T=2
2cb
C cb
2pc
b

#ðcb þ N T =2Þ
TX¯ðJðt Þ  J̄ðt ÞÞV d J ðJ ðtÞ  J ðt ÞÞ;
t¼1

where the expectation value of the current and the inverse noise
variance obtained in the above algorithm are denoted by J̄ and b̄ todistinguish them from the variables J and b, respectively. The
distribution is the student t distribution. By using this posterior
distribution, it is possible to calculate the confidence interval and tvalue statistics.
References
Ahlfors, S.P., Simpson, G.V., Dale, A.M., Belliveau, J.W., Liu, A.K.,
Korvenoja, A., Virtanen, J., Huotilainen, M., Tootell, R.B.H., Aronen,
H.J., Ilmoniemi, R.J., 1999. Spatiotemporal activity of a cortical
network for processing visual motion revealed by MEG and fMRI.
J. Neurophysiol. 82, 2545 – 2555.
Aine, C., Huang, M., Stephen, J., Christner, R., 2000. Multistart algorithms
for MEG empirical data analysis reliably characterize locations and timecourses of multiple sources. NeuroImage 12, 159 – 172.
Amari, S., 1998. Natural gradient works efficiently in learning. NeuralComput. 10, 251 – 276.
Attias, H., 1999. Inferring parameters and structure of latent variable
models by variational Bayes. Proc. 15th Conference on Uncertainty in
Artificial Intelligence. Morgan-Kauffman, pp. 21 – 30.
Baillet, S., Garnero, L., 1997. A Bayesian approach to introducing
anatomo-functional priors in the EEG/MEG inverse problem. IEEE
Trans. Biomed. Eng. 44, 374 – 385.
Bandettini, P.A., 2000. The temporal resolution of functional MRI. In:
Moonen, C.T.W., Bandettini, P.A. (Eds.), Functional MRI. Springer,
Berlin, pp. 205 – 220.
Belliveau, J.W., Kennedy, D.N., McKinstry, R.C., Buchbinder, B.R.,
Weisskoff, R.M., Cohen, M.S., Vevea, J.M., Brady, T.J., Rosen, B.R.,
1991. Functional mapping of the human visual cortex by magnetic
resonance imaging. Science 254, 716 – 719.
Churchland, P.S., Sejnowski, T.J., 1988. Perspectives on cognitive neuroscience. Science 242, 741 – 745.
Dale, A.M., Sereno, M.I., 1993. Improved localization of cortical activity
by combining EEG and MEG with MRI cortical surface reconstruction:a linear approach. J. Cogn. Neurosci. 5, 162 – 176.
Dale, A.M., Liu, A.K., Fischl, B.R., Buchner, R.L., Belliveau, J.W.,
Lewine, J.D., Halgren, E., 2000. Dynamic statistical parametric825

mapping: combining fMRI and MEG for high-resolution imaging ofcortical activity. Neuron 26, 55 – 67.
Engel, S.A., Rumelhart, D.E., Wandell, B.A., Lee, A.T., Glover, G.H.,
Chichilnisky, E.J., Shadlen, M.N., 1994. fMRI measurements of humanvisual cortex. Nature 369, 525.
Engel, S.A., Glover, G.H., Wandell, B.A., 1997. Retinotopic organizationin human visual cortex and the spatial precision of functional MRI.
Cereb. Cortex 7, 181 – 192.
Fujimaki, N., Hayakawa, T., Nielsen, M., Kosche, T.R., Miyauchi, S., 2002.
An fMRI-constrained MEG source analysis with procedures for
dividing and grouping. NeuroImage 17, 324 – 343.
Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 1995. Bayesian Data
Analysis. Chapman & Hall, London.
Goebel, R., Max Planck Society, 2000. Brain Voyager 2000, Version 4.2.
http://www.brainvoyager.de.
Grave de Peralta Menendez, R., Gonzalez Andino, S., 1998. A criticalanalysis of linear inverse solutions to the neuroelectromagnetic inverse
problem. IEEE Trans. Biomed. Eng. 45, 440 – 448.
Hamalainen, M.S., Ilmoniemi, R.J., 1994. Interpreting magnetic fields of thebrain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35 – 42.
Hamalainen, M.S., Hari, R., Ilmoniemi, R.J., Knuutila, J., Lounasmaa,
O.V., 1993. Magnetoencephalography—Theory, instrumentation, andapplications to noninvasive studies of the working human brain. Rev.
Modern Phys. 65, 413 – 497.
Hari, R., 1991. On brain’s magnetic responses to sensory stimuli. J. Clin.
Neurophysiol. 8, 157 – 169.
Heinze, H.J., Mangun, G.R., Burchert, W., Hinrichs, H., Scholz, M., Munte,
T.F., Gos, A., Scherg, M., Johannes, S., Hundeshagen, H., Gazzaniga,
M.S., Hillyard, S.A., 1994. Combined spatial and temporal imaging of
brain activity during visual selective attention in humans. Nature 372,
543 – 546.
Huang, M., Aaron, R., Shiffman, C.A., 1997. Maximum entropy method
for magnetoencephalography. IEEE Trans. Biomed. Eng. 44, 98 – 102.
Kajihara, S., Ohtani, Y., Goda, N., Tanigawa, M., Toyama, K., 2004.
Wiener filter magnetoencephalography of visual cortical activities.
Brain Topogr. (in press).
Linden, D.E.J., Kallenbach, U., Heineckeo, A., Singer, W., Goebel, R., 1999.
The myth of upright vision. A psychophysical and functional imaging
study of adaptation to inverting spectacles. Perception 28 (4), 469 – 481.
Liu, A.K., Belliveau, J.W., Dale, A.M., 1998. Spatiotemporal imaging of
human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proc. Natl. Acad.
Sci. U. S. A. 95, 8945 – 8950.
Logothetis, N.K., Pauls, J., Augth, M., Trinth, T., Oeltermann, A., 2001.
Neurophysiological investigation of the basis of the fMRI signal.
Nature 412, 150 – 157.
Mosher, J.C., Lewis, P.S., Leahy, R.M., 1992. Multiple dipole modelling
and localization from spatio-temporal MEG data. IEEE Trans. Biomed.
Eng. 39, 541 – 557.
Neal, R.M., 1996. Bayesian Learning for Neural Networks. SpringerVerlag, New York.
Nunez, P.L., 1981. Electric Fields of the Brain: The Neurophysics of EEG.
Oxford Univ. Press, New York.
Ogawa, S., Lee, T.-M., Kay, A.R., Tank, D.W., 1990. Brain magnetic
resonance imaging with contrast-dependent oxygenation. Proc. Natl.
Acad. Sci. U. S. A. 87, 9868 – 9872.
Parisi, G., 1988. Statistical Field Theory. Perseus Books.
Pascual-Marqui, R.D., 1999. Review of methods for solving the EEG
inverse problem. Int. J. Bioelectromagn. 1, 75 – 86.
Pascual-Marqui, R.D., Michel, C.M., Lehmann, D., 1994. Low resolution
electromagnetic tomography: a new method for localizing electricalactivity in the brain. Int. J. Psychophysiol. 18, 49 – 65.
Penny, W., Kiebela, S., Friston, K., 2003. Variational Bayesian inference for
fMRI time series. NeuroImage 19, 727 – 741.
Phillips, C., Rugg, M.D., Friston, K.J., 2002. Anatomically informed basis
functions for EEG source localization: combining functional andanatomical constraints. NeuroImage 16, 678 – 695.

826

M. Sato et al. / NeuroImage 23 (2004) 806–826

Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem. Phys. Med. Biol. 32, 11 – 22.
Sato, M., 2001. On-line model selection based on the variational Bayes.
Neural Comput. 13, 1649 – 1681.
Schmidt, D.M., George, J.S., Wood, C.C., 1999. Bayesian inference applied
to the electromagnetic inverse problem. Hum. Brain Mapp. 7, 195 – 212.
Sereno, M.I., Dale, A.M., Reppas, J.B., Kwong, K.K., Belliveau, J.W.,
Brady, T.J., Rosen, B.R., Tootell, R.B.H., 1995. Borders of multiple
visual areas in human revealed by functional magnetic resonanceimaging. Science 268, 889 – 893.

Toyama, K., Yoshikawa, K., Yoshida, Y., Kondo, S., Tomita, Y., Takanashi,Y., Ejima, Y., Yoshizawa, S., 1999. A new method for magnetoencephalography (MEG): three dimensional magnetometer-spatial filter
system. Neuroscience 91, 405 – 415.
Uutela, K., Hamalainen, M.S., Somersalo, E., 1999. Visualization ofmagnetoencephalographic data using minimum current estimates.
NeuroImage 10, 173 – 180.
Wang, J.-Z., Williamson, S., Kaufman, J., 1992. Magnetic source images
determined by a lead-field analysis: the unique minimum-norm leastsquares estimation. IEEE Trans. Biomed. Eng. 39, 665 – 675.