In the biomagnetic inverse problem the main interest is the activation of a region of interest, i.e. the power dissipated in that region. The Bayesian power imaging method (BPI) provides a quantified probability that the activation of a region of interest is above a given threshold. This paper introduces the method and derives the equations used. The method is illustrated in this paper using both experimental and simulated data

Hasson, R.

Swithenby, S. J.

English

Open Research Online (The Open University)

Open Research OnlineThe Open University’s repository of research publicationsand other research outputsThe Bayesian Power Imaging (BPI) method formagnetic source imagingConference or Workshop ItemHow to cite:Hasson, R. and Swithenby, S. J. (2001). The Bayesian Power Imaging (BPI) method for magnetic sourceimaging. In: Biomag2000, Proc. 12th Int. Conf. on Biomagnetism, Aug 2000, Helsinki, Finland, Helsinki Univ. ofTechnology, Espoo, Finland, pp. 697–700.For guidance on citations see FAQs.c© 2001 The AuthorVersion: Version of RecordLink(s) to article on publisher’s website:http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.5.2153&rep=rep1&type=pdfCopyright and Moral Rights for the articles on this site are retained by the individual authors and/or other copyrightowners. For more information on Open Research Online’s data policy on reuse of materials please consult the policiespage.oro.open.ac.ukThe Bayesian Power Imaging (BPI) methodfor magnetic source imagingR. Hasson1 and S.J. Swithenby21Department of Applied Mathematics, The Open University, Milton Keynes, United Kingdom;2Department of Physics, The Open University, Milton Keynes, United Kingdom1 IntroductionIn the biomagnetic inverse problem the main inter-est is the activation of a region of interest, i.e. thepower dissipated in that region. The Bayesian powerimaging method (BPI) provides a quantified probabil-ity that the activation of a region of interest is abovea given threshold. This paper introduces the methodand derives the equations used. The method is illus-trated in this paper using both experimental and sim-ulated data. Another paper [1] in this volume extendsthe method to compare task and control experiments.2 MethodsBefore we can begin to describe the method we needa few standard definitions, as used in [2]. Let m 2RN be the vector of measurements at a single timeinstant, i.e. miis the measurement of the ith channel,i = 1; : : : N . Let ~Li(~r) be the lead field of the ithdetector.For simplicity we will use the same Hilbert space ofcurrents as in [2], namely L2(Q) the square integrablefield defined on the brain volume Q with the follow-ing inner product between currentsD~j1;~j2E=ZQ~j1(~r) ~j2(~r)!(~r) d~rwhere !(~r) is a weighting distribution defined on thesource space Q. Define the Gramm-Schmidt matrixP , using this inner product, by Pik=D~Li;~LjE.This section continues by first deriving expressionsfor the a posteriori probability distribution on themeasurements using the results from [2]. In turn thisis used to derive an a posteriori probability distribu-tion on the activation of brain regions. It is this latterprobability distribution that is mapped out in sourcespace.2.1 Distribution on measurementsThe Bayesian method derived in [2] provides the aposteriori probability distribution of current distribu-tions~j(~r). The results of [2] are stated in terms of theprojection of this a posteriori probability distributiononto one dimensional subspaces of L2(Q) by usinga test current ~t(~r). A special case (namely when theassumed prior current distribution ~jprior is identicallyzero) of the results derived in [2] are that the expectedmean is uT(P + D) 1m and the expected varianceis:2uT PD 1P + P 1u+1 ~t ;~t  uTP 1uwhere ui=D~t(~r) ; !(~r)~Li(~r)E. 2D is the covari-ance matrix of the assumed Gaussian noise and  is aregularization parameter. In this Bayesian setting theregularization parameter  is equal to 2=2 where2 is the variance of the assumed prior probabilitydistribution on currents~j(~r).In [2] the test current ~t was scanned across a set ofvoxels to produce an image. Here the approach is dif-ferent. By setting ~t(~r) = !(~r)~Li(~r) we can projectthe a posteriori distribution onto the expected mea-surement in the ith channel. So, the probability dis-tribution represents our best guess of the true value ofthe measurements given our prior knowledge of themeasurement geometry and the data vector m. Theexpected mean and variance can be calculated usingthe equations above. The formulae can be simplifiedby noting that uj=D!(~r)~Li(~r) ; !(~r)~Lj(~r)E= Pij,i.e. the vector u is the ith column of P . It follows thatthe mean and variance for all channels can be calcu-lated in a single matrix equation. In this case the ex-pected mean is P (P + D) 1m and the covariancematrix is given by2P PD 1P + P 1P +1 P   PP 1Pwhich simplifies to 2D(P + D) 1P .A few checks on the reasonableness of the aboveequations are possible by looking at extreme cases ofthe regularization parameter  .As  ! 0, the mean tends to m and the covariancetends to 2D. This is reasonable because, in thisBayesian setting,  ! 0 corresponds to having ei-ther perfect data (i.e.  ! 0) or no prior knowledge(i.e.  ! 1). In this case, the formulae show thatthe best estimate of the measurements is given by thedata with covariance equal to the assumed covarianceof the noise.As  !1, both the mean and the covariance matrixtend to zero. This is also reasonable since  ! 1corresponds to either worthless data (i.e.  ! 1) orperfect prior knowledge (i.e.  ! 0). The formulaeabove indicate that the best estimate of the measure-ments is equal to our prior prejudices (in this case zerocurrent density) with complete certainty.2.2 Distribution on activation estimatesThe above formulae can be used in the direct powerestimation algorithm derived previously [3]. In thismethod an optimal estimation matrix Yfor the thvoxel is derived and the expression mTYm is used toestimate the activation of the th voxel. The expres-sion derived in [3] for Yis(P2+D)Y(P2+D) = PXP (1)where Xis a matrix which is related to the choice ofregion of interest. In this paper the region of interest isa single voxel in a source space so the expression forXis simply LkLkTwhere Lkis a column of the gainmatrix for the problem, i.e. Lkis the measurementvector that would result from a dipole at position ~roriented in direction bek.The matrix Y in Equation 1 is then given by;Y=2Xk=0	k	kT (2)where 	k= (P2+D) 1PLk. For this special caseit is more efficient not to form the matrix Y but tocalculate the activation directly using 	k;A(m) = mTYm =2Xk=0	kTm2(3)Here the method is put in a Bayesian context be-cause we now know the a posteriori distribution ofmeasurement channels. Unfortunately the probabil-ity distribution of the quadratic form mTYm is notknown in a closed form. Analytic approximationshave been derived to calculate the probability distri-bution of quadratic forms in Gaussian random vari-ables (e.g. [4]) but these methods are computationallyexpensive. In practice the quickest way of achieving agiven level of accuracy for the probability is by MonteCarlo integration. First the measurements are normal-ized by factoring out the known covariance matrix:m?=1(P + D) 1=2P1=2m (4)Then the probability can be estimated by picking MGaussian random vectors of length N with zero meanand unit covariance matrix, say zi, i = 1; : : : M , andcomputing the probability:P(x) =1MMXi=1(x A(zi)) (5)where P(x) is the one-sided probability of the esti-mate x and  is the Heaviside step function.Maps of probability can then be computed by scan-ning  across a defined source space.3 ResultsThe effectiveness of the algorithms has been tested intwo simulated experiments and on real data.The first simulated experiment has a simple geometry(Figure 1). The head is modelled as a homogeneousconducting sphere of radius 8:9 cm with its centre at(0; 0; 0:07 cm). The source space is a 6 cm6 cmsquare thin lamina consisting of 3333 voxels in theplane z =  0:01 cm with centre (0,0,-0.01 cm). Themeasurement instrument is a hexagonal array of 37second order axial gradiometers with baseline 5 cmwith the lowest ’sensing’ coils in the plane z = 4 cm.The sources for the first numerical study weretwo current dipoles with positions (0 cm, 0 cm) and( 2:25 cm, 2:25 cm) and oriented along the x andy axes respectively. Uncorrelated Gaussian noise isadded to the pure signal so that the resulting signal tonoise ratio was 9:7.In analyzing the experimental data there are manypossible methods of determining noise levels etc. Inthis paper we have tried to keep the analysis choicessimple in order to concentrate on the fundamentalsof the method. With this in mind, all the examples-0.06 0.0 0.06x/m-0.060.00.06y/mFigure 1: A plan view of the experiment geometry.The grid depicts the source space voxels and the cir-cles represent the position of the measurement coilsof the second order gradiometers.have the same choices for the parameters: the noiseis assumed to be uncorrelated Gaussian noise of vari-ance 2, the regularization parameter is fixed to be0:1  trace(P )=N and the number of trials, M , forthe Monte-Carlo probability estimation method was1000.If the Monte-Carlo probability estimation method isapplied directly to the raw power estimates, large ar-eas of the source space are found to be significant.This is because it is very unlikely for there to be zeropower dissipated in each voxel. An easy calculationshows that the expected power for the th voxel isgiven by trace(Y). The sensitivity to the (unknown)variance of the noise in the measurements can be min-imized by using the mean and standard deviation ofthe power estimates across the source space as a mea-sure of expected power. This is equivalent to the as-sumption that the desired signal produces significantpower estimates on only a small proportion of thesource space. In this paper the examples show prob-ability maps of the probability that the activation ismore than three standard deviations above the meanactivation. The resulting thresholded image is shownin Figure 2. The two sources are located accuratelyand with a very high significance.In the second example, a current dipole source wasembedded in the innermost shell of a 3-shell modelconducting sphere. The shell radii are 73mm, 80mmand 86mm and the corresponding conductivities are0:33 Sm 1, 0:0042 Sm 1 and 0:33 Sm 1. A detectorsystem consisting of 65 EEG electrodes was used tomeasure the scalp potential. The electrode positions-0.03 0.0 0.03x/m-0.030.00.03y/mFigure 2: A probability map across the source spaceshown in Figure 1. The scale goes from white (repre-senting a probability of less than 0.1) to black (prob-ability greater than 0.9).and a representation of the potentials measured by theelectrodes are shown in Figure 3. Also shown in Fig-ure 3 is the conducting shell used as the source spacefor the reconstruction.Figure 3: An EEG simulated dipolar dataset, togetherwith the corresponding BPI map. The white blob thatappears on the shell is the region where the activationis greater than the threshold (at the 99.9% significancelevel).The third and last example involves real MEG datafrom an evoked response study of face-processing [6].The experimental instrument used is the Neuromag-122TM [5]. Human subjects were presented brieflywith photographs of human faces. It is known thatthe early response to face images involves widespreadactivity in the posterior brain but there is limited ev-idence for the precise distribution and timecourse ofthe neuronal sources. One suggestion is that there arethree major areas of activity; in occipital cortex andboth right and left ventral occipito-temporal cortex[7]. Strong occipital activity (starting about 100 msafter the stimulus) is expected to lead to concurrentactivity in the two other regions with a stronger re-sponse in the right hemisphere [6].The result shown in Figure 4 is derived from the re-sponse in an individual subject at a latency of 160msafter the stimulus. At this time, there is a globalpeak in the total field power for this experiment. Thesource space was restricted to the surface of the cor-tex, as revealed by MRI and also shown in Figure 4.The findings of the BPI algorithm are consistent withearlier studies in that there is a significant response inboth right and left cortices, with a stronger responsein right hemisphere. Further work is needed to refinethese predictions as the source space has a complexfolded structure and we have not yet disentangled thecontributions of superficial and deeper source regionsFigure 4: The BPI map on the cortical source spaceused. The 122-channel helmet detector system is alsoshown.4 DiscussionThe main impetus behind the development ofquadratic algorithms to estimate power rather thancurrent density is that other modalities (e.g. fMRI,PET) naturally produce images of power dissipated.So if multi-modal comparisons are to be made then itis sensible to produce power images from MEG data.The BPI method has the additional advantage of be-ing able to produce spatial images of the variation ofprobability that are easy to interpret.In this paper the BPI method was shown to work ef-fectively for simulated MEG and EEG data. It wasalso demonstrated for experimental MEG data andproduced results consistent with earlier analyses.AcknowledgementsWe would like to thank Dr. Sven Bräutigam for manyhelpful conversations and for the preparation of theface processing data.References[1] R. Hasson, J.A. James, and S.J. Swith-enby. The Bayesian Power Imaging (BPI)test for task/control experiments. Submitted toBiomag2000, Helsinki, August 2000.[2] R. Hasson and S.J. Swithenby. A Bayesian testfor the appropriateness of a model in the bio-magnetic inverse problem. Inverse Problems,15(6):1439–1454, December 1999.[3] R. Hasson and S.J. Swithenby. Activationcurves from optimally shaped regions. InT. Yoshimoto, M. Kotani, S. Kuriki, H. Karibe,and N. Nakasato, editors, Recent Advances inBiomagnetism, pages 205–208. Tohoku Univer-sity Press, Sendai, 1999. ISBN 4-925085-19-0C3047.[4] J.P. Imhof. Computing the distributionof quadratic forms in normal variables.Biometrika, 48(3/4):419–426, December 1961.[5] J.E.T. Knuutila, A.I. Ahonen, M.S. Hämäläinen,M.J. Kajola, P.P. Laine, O.V. Lounasmaa, L.T.Parkkonen, J.T.A Simola, and C.D. Tesche. A122-channel whole cortex squid system for mea-suring the brains magnetic-fields. IEEE Trans.Mag., 29 (6):3315–3320, 1993.[6] S.J. Swithenby, A.J. Bailey, S. Bräutigam, O.E.Josephs, V. Jousmäki, and C.D.Tesche. Neu-ral processing of human faces: a magnetoen-cephalographic study. Expt. Brain Res., 118(4):501–510, 1998.[7] S.T. Lu, M.S. Hämäläinen, R. Hari, R.J. Il-moniemi, O.V. Lounasmaa, M. Sams, andV. Vilkman. Seeing faces activates 3 separateareas outside the occipital visual cortex in man.Neuroscience, 43 (2/3):287–290, 1991.

The Bayesian Power Imaging (BPI) method for magnetic source imaging

Abstract

Similar works

Full text

Available Versions

Open Research Online (The Open University)