Transcript Bayesian filtering in MEG - Austrian Academy of Sciences

Particle Filtering in MEG: from single
dipole filtering to Random Finite Sets
A. Sorrentino
CNR-INFM LAMIA, Genova
methods for image and data analysis
[email protected]
www.dima.unige.it/~piana/mida/group.html
Co-workers
Genova group:
Cristina Campi
Annalisa Pascarella
Michele Piana
(Math Dep.)
(Comp. Sci. Dep.)
(Math. Dep.)
Long-time collaboration
Lauri Parkkonen
(Brain Research Unit, LTL, Helsinki)
Recent collaboration
Matti Hamalainen
(MEG Core Lab, Martinos Center, Boston)
Basics of MEG modeling
Biot-Savart
Poisson
bt ( r )

0
4

[ jt ( r ')  ( r ') Vt ( r ')]

j t ( r ' )  [ ( r ' ) Vt ( r ' )]
j t (r ' )
Accurate model
of brain
conductivity
 ( r ' ) Vt ( r ' )
( r  r ')
|( r  r ')|
3
Neural
current
dr '
Ohmic
term
Biot-Savart
bt (r )
2 approaches to MEG source modeling
Imaging approach
Parametric approach
Continuous current distribution
Model
N
jt (r )   qti (r  r i ) N large
i 1
Unknown
Method
Result
i
bt ( r )  G ({r }, r )  Qt
Regularization methods
Focal current
Mt
jt (r )   qti (r  rti )
M small
i 1
Mt
bt ( r ) 

i 1
qti 
r  rti
|r  rti |3
Non-linear optimization methods
Automatic current dipole estimate
Common approximations to solve this problem:
 Number of sources constant
 Source locations fixed
Mt
bt ( r ) 
r  rti
 qti |r r i |3
i 1
t
Common methods:
 Manual dipole modeling
 Automatic dipole modeling
Estimate the number of sources
Estimate the source locations
Least Squares for source strengths
Manual dipole modeling still the main reference method for comparisons
(Stenbacka et al. 2002, Liljestrom et al 2005)
Bayesian filtering allows overcoming these limitations
Bayesian filtering in MEG - assumptions
Two stochastic processes:
Markovian assumptions:
Our actual model
J1
J2
B1
B2
…
…
Jt …
Bt …
 ( jt 1 | jt ,..., j1 )   ( jt 1 | jt )
Markov process
Instantaneous propagation
 (bt | jt ,..., j1 )   (bt | jt )
 ( jt 1 | jt , b1 ,..., bt )   ( jt 1 | jt ) No feedback
J t 1  J t  J t
Bt  B( J t )  N t
The final aim:
 ( j1 | b1 ),  ( j2 | b2 ), ...,  ( jt | bt ),...
Bayesian filtering in MEG – key equations
Likelihood
function
“Observation”
 post ( jt | b1:t ) 
 (bt | jt )   prior ( jt | b1:t 1 )
 (bt )
b1
b2
Transition
kernel
“Evolution”
 prior ( j1 )
 post ( j1 | b1 )
 prior ( j2 | b1 )
 post ( j2 | b1:2 )
…
…
 prior ( jt 1 | b1:t )    ( jt 1 | jt )   post ( jt | b1:t ) djt
bt
Linear-Gaussian model  Kalman filter
Non-linear model  Particle filter
 prior ( jt | b1:t 1 )
 post ( jt | b1t: )
E
S
T
I
M
A
T
E
S
Particle filtering of current dipoles
The key idea: sequential Monte Carlo sampling.
jt  D t (single dipole space)
Draw random samples
(“particles”) from the prior
 prior ( j1 )
Update the
particle weights
 post ( j1 | b1 )
Resample and let
particles evolve
 prior ( j2 | b1 )
i
1 i 1,..., N
{j }
 prior ( j1 )   ( j1  j1i )
i
 post ( jt | b1 ,...,bt ) 
i
1
i
1 i 1,..., N
{j ,w }
 (bt | jt )   prior ( jt | b1 ,...,bt 1 )
 (bt )
 post ( j1 )   w1i ( j1  j1i )
i
 prior ( jt 1 | b1 ,...,bt )    ( jt 1 | jt )   post ( jt | b1 ,...,bt ) djt
{ j2i }i 1,..., N
 prior ( j2 | b1 )   ( j2  j2i )
i
A 2D example – the data
A 2D example – the particles
The full 3D case – auditory stimuli
S. et al., ICS 1300 (2007)
Comparison with beamformers and RAP-MUSIC
Pascarella et al., ICS 1300 (2007); S. et al. , J. Phys. Conf. Ser. 135 (2008)
Two quasi-correlated sources
Beamformers: suppression of correlated sources
Comparison with beamformers and RAP-MUSIC
Pascarella et al., ICS 1300 (2007); S. et al. , J. Phys. Conf. Ser. 135 (2008)
Two orthogonal sources
RAP-MUSIC: wrong source orientation,
wrong source waveform
Rao-Blackwellization
Campi et al. Inverse Problems (2008); S. et al. J. Phys. Conf. Ser. (2008)
  
 
Q
0 i  (r  rQi )  V  
B(r )  
  3  B (rQi )  Qi
4

|
r
 rQi |
i
Can we exploit the linear substructure?
 post ( jt | bt )   post (rt , qt | bt )   post (qt | rt , bt ) post (rt | bt )
Analytic solution
(Kalman filter)
Sampled
(particle filter)
Accurate results with much fewer particles
Statistical efficiency increased (reduced variance of importance
weights)
Increased computational cost
Bayesian filtering with multiple dipoles
A collection of spaces (single-dipole space D, double-dipole space,...)
A collection of posterior densities (one on each space)
Exploring with particles all spaces (up to...)
1 : D  

One particle = one dipole
2 : D D 

One particle = two dipoles
 3 : D  D D  
One particle = three dipoles
D
D D
Reversible Jumps (Green 1995) from one space
to another one
D D  D
Random Finite Sets – why
Non uniquess of vector representations of multi-dipole states:
(dipole_1,dipole_2) and (dipole_2,dipole_1)
same physical state, different points in D X D
Consequence: multi-modal posterior density
non-unique maximum
non-representative mean
Let (,,P) be a
probability space
 2 : D  D  
 2 (d1 , d2 )   2 (d2 , d1 )
A random finite set X of dipoles
is a measurable function
X :   P ( D)
Where P(D) is the set of all finite
subsets of D (single dipole space)
equipped with the
Mathéron topology
For some realizations,
X( )  {d }
X( ' )  {0}
X( ' ' )  {d1 , d 2 , d3}
Random Finite Sets - how
Probability measure of RFS: a conceptual definition PX ( A)  P({ | X()  A})
 Belief measure instead of probability measure
X (C)  P({ | X()  C})
A  B( F )
C  D
Multi-dipole belief measures can be derived from single-dipole probability measures
 Probability Hypothesis Density (PHD): the RFS-analogous of the conditional mean
The integral of the PHD in a volume = number of dipoles in that volume
Peaks of the PHD = estimates of dipole parameters
Model order selection: the number of sources estimated dynamically
RFS-based particle filter: Results
S. et al., Human Brain Mapping (2009)
Monte Carlo simulations:
1.000 data sets
Random locations (distance >2 cm)
Always same temporal waveforms
 2 time-correlated sources
 peak-SNR between 1 and 20
Results:
 75% sources recovered (<2 cm)
 Average error 6 mm, independent
on SNR
 Temporal correlation affects the
detectability very slightly
RFS-based particle filter: Results
S. et al., Human Brain Mapping (2009)
Comparison with manual dipole modeling
Data: 10 sources mimicking complex visual activation
The particle filter performed on average like manual dipole modeling performed by
uninformed users (on average 6 out of 10 sources correctly recovered)
In progress
Source space limited to the cortical surface
Two simulated sources
In progress
Two sources recovered with
orientation constraint
Only one source recovered
without orientation constraint
References
- Sorrentino A., Parkkonen L., Pascarella A., Campi C. and Piana M. Dynamical MEG source modeling with
multi-target Bayesian filtering Human Brain Mapping 30: 1911:1921 (2009)
-Sorrentino A., Pascarella A., Campi C. and Piana M. A comparative analysis of algorithms for the
magnetoencephalography inverse problem Journal of Physics: Conference Series 135 (2008) 012094.
-Sorrentino A., Pascarella A., Campi C. and Piana M. Particle filters for the magnetoencephalography
inverse problem: increasing the efficiency through a semi-analytic approach (Rao-Blackwellization)
Journal of Physics: Conference Series 124 (2008) 012046.
-Campi C., Pascarella A., Sorrentino A. and Piana M. A Rao-Blackwellized particle filter for
magnetoencephalography Inverse Problems 24 (2008) 025023
- Sorrentino A., Parkkonen L. and Piana M. Particle filters: a new method for reconstructing multiple
current dipoles from MEG data International Congress Series 1300 (2007) 173-176