Transcript Source Estimation in EEG - University College London

Source Estimation in EEG
The forward and inverse
problems
Christophe Phillips, Ir, Dr
Cyclotron Research Centre, University of Liège, Belgium
Agenda
• Introduction
Why , what for, where, how...
• Part I : the forward problem
From sources to electrodes
• Part II : the inverse problem
From electrodes to sources
EEG Recordings
Electroencephalography (EEG) is ‘simply’ about
recording electromagnetic signals produced by
neuronal activity :
EEG signal is spread in space and time.
EEG signal: origin
Neurone
Cell body
Dendrites
Head anatomy:
•gray matter,
•white matter,
•CSF,
•bone, air, skin, mucle, etc.
Axon
Synaptic
terminals
EEG signal: origin
• From a distance, postsynaptic potential (PSP) looks
like a current dipole oriented along the dentrite.
• Pyramidal cells have parallel dendrites, oriented
perpendicular to the cortical surface.
• Typical dipole strength : ~20 fA m (1fA = 10 A)
• Current-dipole moments required to be measured
outside the head : ~10 nA m (1nA = 10 A)
-15
-9
About 106 synapses
must be simultaneously
and effectively active to
produce an evoked
response.
About 40 to 200 mm2
of active cortex.
Cell body
Dendrites
Synaptic
terminals
Axon
EEG signal: origin
What about action potentials ?
• large but brief potential compared to PSP
PSP
AP
100 mV
10 mV
10 ms
1 ms
• modelled by a current-quadrupole
 the field decreases with distance as 1/r3,
compared to 1/r2 for the PSP dipole.
EEG signals are produced in large
part by synaptic current flow,
which is approximately dipolar.
EEG vs fMRI Recordings
•fMRI
•EEG
Signal change
Difference between haemodynamic and
electromagnetic signals produced by neuronal
activity as recorded by :
Fitted haemodynamic
response function
Preconditions for signal
detection
EEG signal
•Activation of a neural
population must be
synchronous in time.
•Active neural population
must be spatially
organised (parallel
fibers).
•The sources must be in
an ‘open-field’
configuration.
fMRI signal
•Neural activity needs not
be synchronous in time.
•Geometrical orientation
of the sources is totally
irrelevant.
Signal Sensitivity
EEG signal
fMRI signal
•Measured signal is
sensitive in relative
timing and amplitude
of activity.
•Measured signal
amplitude influenced
by both duration and
amplitude of activity.
•Critical neural activity
needs not be extended
in time.
•Signal detected only if
net haemodynamic
changes
Source Localisation in EEG
Forward Problem
Inverse Problem
Forward Problem
Head anatomy
Neurone
Cell body
Axon
Dendrites
Synaptic
terminals
Head model : conductivity layout
Source model : current dipole
Solution by Maxwell’s equations
Inverse Problem
Solving
 
v  f (r , j )  
v, potential at the electrodes, (Nel x 1)
, source(s) location,
r
where
[rx , ry , rz]’

j , source(s) orientation & amplitude
[jx , jy , jz]’
f, solution of the forward problem
, additive noise, (Nel x 1)
is an “ill-posed” problem, i.e. the inverse
problem does NOT have a unique solution,
(von Helmholtz, 1853).
Source Estimation in EEG
Part I : The forward problem
Agenda
Part I : Solving the forward problem
• Maxwell’s equations
• Analytical solution
• Numerical solution
• Pseudo-sphere approach
• Conclusion
Maxwell’s equations (1873)
 
Ohm’s law :
E 


j  E
 
B  0 
 
B
Continuity equation :
E  
t






 
E

j

  B  j  

t
t
Quasistatic approximation
Maxwell’s equations can be simplified
because :
• EEG frequencies are genrally below 100hz.
• Cellular electrical phenomena contain
mostly frequencies below 1kHz.
No propagation phenomenon,
i.e. changes in the field/sources are
«instantaneous».
Derivative terms can be discarded.
Simplified Maxwell’s equation
With
• quasi-static approximation of
 Maxwell’s equations
• and dipolar current sources j


 

  V   j
 conductivity
where V
 electric potential
j current source
i.e. a “simple” mathematical relationship
linking current sources and electric
potential, depending on the conductivity
of the media.
Solving the forward Problem
From


 

  V   j
 
find V  f ( j ,r )
where f(.) will depend on the conductivity 
of the media (and the boundary conditions).
The head
model is the
conductivity
layout
adopted !
Solving the forward Problem
From


 

  V   j
find
 
V  f ( j ,r )
f(.) can be estimated
• analytically, i.e. an exact solution
‘formula’ exists,
• numerically, i.e. numerical methods
provide an approximation of the
solution.
Agenda
Part I : Solving the forward problem
• Maxwell’s equations
• Analytical solution
• Numerical solution
• Pseudo-sphere approach
• Conclusion
Analytical solution
From


 

  V   j
find
 
V  f ( j ,r )
f(.) can be estimated analytically only for
particular cases:
• higly symmetrical geometry, e.g.
spheres, cubes, concentric spheres, etc.
• homogeneous isotropic conductivity
These are very restricted head models!
Analytical solution, cont’d
Human head is not spherical neither
is its conductivity homogeneous and
isotropic...
Agenda
Part I : Solving the forward problem
• Maxwell’s equations
• Analytical solution
• Numerical solution
• Pseudo-sphere approach
• Conclusion
Numerical solution, 1
From


 

  V   j
find
 
V  f ( j ,r )
f(.) can be estimated numerically for ANY
conductivity layout!
The most general method is the “Finite
Element Method”, or FEM:
• conductivity can be arbitrary, i.e.
anisotropic and inhomogeneous,
• potential is estimated throughout the
volume.
Numerical solution, FEM
Principles of the FEM:
• The (head) volume is tesselated into
small volume elements on which
Maxwell’s equation is solved locally.
• The conductivity is defined for every
volume elements individually.
Drawbacks:
• How to build up the model and define
the conductivity at each element ?
• Computation time is huge !
Numerical solution, 2
From


 

  V   j
find
 
V  f ( j ,r )
f(.) can be more easily estimated numerically
with some assumptions:
• volume divided into sub-volumes of
homogeneous and isotropic conductivity,
• potential is only estimated on the
surfaces seperating those sub-volumes.
This is the “Boundary Element Method”, or
BEM.
Numerical solution, BEM
The surfaces are tessellated into flat triangles and
the potential is approximated on each triangle as
a constant or linear function.
Example of BEM
head model :
Brain surface
The sources, current
dipoles, are placed
in the brain volume.
Skull surface
Scalp surface
BEM, Potential approximation
On each triangle of the surfaces, the potential
function can be approximated by :
a linear function between the
potential at the vertices (LPV)

V ( s3 )
V

V ( s1 )

s1

s3

V ( s2 )

s2
BEM application
Homogeneous volumes definition
BEM head
model
Scalp surface
with electrode
locations
BEM limitations
Superficial dipoles have
sharper potential
distribution.
BEM fails when the size of
surface elements is ‘large’
compared to the sharpness
of potential distribution.
Agenda
Part I : Solving the forward problem
• Maxwell’s equations
• Analytical solution
• Numerical solution
• Pseudo-sphere approach
• Conclusion
Solutions: Analytical vs. numerical
The head is NOT spherical:
 cannot use the exact analytical solution
because of model/anatomical errors.
Realistic model needs BEM solution:
 surfaces extraction
 computationnaly heavy
 errors for superficial sources
Could we combine the advantages of both
solutions ?
Anatomically constrained spherical head
models, or pseudo-spherical model.
Pseudo-spherical model
Scalp (or brain) surface
Best fitting sphere:
centre and radii (scalp, skull, brain)
Spherical transformation of
source locations
Leadfield for the spherical model
Applications, scalp surface
Fitted sphere and
scalp surface
Applications, cortical surface
Agenda
Part I : Solving the forward problem
• Maxwell’s equations
• Analytical solution
• Numerical solution
• Pseudo-sphere approach
• Conclusion
Conclusion
Solving the « Forward Problem » is not very
exciting neither easy but it is crucial for any
attempt at source estimation.
The key elements are the head model and, from it,
the solution used.
Model :
simple
 define 3 sphere shell
realistic
 extract volume surfaces
most realistic  extract volume conductivity
Solution :
analytic  easy and quick but anatomical errors
numeric  slower, more anatomically accurate but
numerical erros
So far we still have NOT
localised anything…
Source Estimation in EEG
Part II : The inverse problem
Agenda
Part II : Solving the inverse problem
• Introduction
• Equivalent current dipole solution
• Distributed linear solution
• Other solutions
• Conclusion
Inverse Problem
Solving
 
v  f (r , j )  
v, potential at the electrodes, (Nel x 1)
, source location,
r
where
[rx , ry , rz]’

j , source orientation & amplitude
[jx , jy , jz]’
f, solution of the forward problem
, additive noise, (Nel x 1)
Function f is
• linear w.r.t. the source orientation & amplitude
• non-linear w.r.t. the source location
Parameters
When Ns sources are present the problem
to solve is
Ns


v   f (rl , j l )  
l 1
For each source, there are 6 parameters :
or
• 3 for the location, [x y z] coordinates,
• 3 for the orientation and amplitude, [jx jy jz]
components
• 3 for the location, [x y z] coordinates
• 2 for the orientation, [q j] angles
• 1 for the amplitude, j intensity/strength
Parameters, cont’d
The inverse problem is « ill posed », i.e. in
general there is no unique solution:
• Number and location of active sources are
unknown!
• Measurements from just Ne electrodes.
To uniquely solve the inverse problem
assumptions/constraints on the solution
MUST be adopted.
Those constraints define
the form of the solution !
Inverse Problem, example
Example :
4 different networks but with the same
measurable output: 2V and 4Ω.
2V
12 Ω 6 Ω
12 Ω 6 Ω
12 Ω
2V
6Ω
4Ω
6V
S1
3V
S2
If we constrain the solution to have
• the smallest source
• the smallest but deeper source
• source along the 12Ω resistor
S3
:
 solution S3
 solution S2
 solution S1
Equivalent Current Dipole
With the ECD solution :
A priori fixed number of sources
considered, usually less than 10
 over-determined but nonlinear problem
 iterative fitting of the 6 parameters of
each source.
Problem : How many ECDs a priori ?
The number of sources limited : 6xNs < Ne
Advantage : Simple focused solution.
But is a single (or 2 or 3 or…) dipole(s)
representative of the cortical activity ?
Distributed Linear solution
With the DL solution :
“All” possible (fixed) source locations (>103)
considered simultaneously
 largely under-determined but linear problem :
v  Lj 
 external constraints required to calculate
a unique solution
Problem :
What should be the constraints ?
Advantage : 3D voxelwise results.
Agenda
Part II : Solving the inverse problem
• Introduction
• Equivalent current dipole solution
• Distributed linear solution
• Other solutions
• Conclusion
ECD solution
For Ns sources, the problem can be rewritten as
Ns
 
 
v   f (rl , j l )   L (rl ) j l  L j
Ns
l 1
l 1
which is an over-determined but non-linear problem.
Once the location is fixed

ˆ
j  Lv
The cost function to be optimised is :
2
2
cost  v  vˆ  v  Lˆ
j  v  LL v

To be iteratively minimised only w.r.t the
i.e. 3 parameters per source.
2

rl,
ECD solution, cont’d
For Ns sources,
the cost function can be minimised

for the rl using any nonlinear procedure,
e.g. gradient descent, simplex algorithm, etc.

cost  v  LL v
2
At each iteration the leadfield L must be recalculated
(many times) as the source locations are modifed
Once the location of the sources is determined,
their intensity is obtained by

ˆ
j Lv
ECD solution, cont’d
The iterative optimisation procedure can only find
a local minimum
the starting location(s) used can influence
the solution found !
1D example of
optimisation
problem:
Cost function
For an ECD solution, initialise the dipoles
• at multiple random locations and repeat the
fitting procedure  focal cluster of solutions ?
• at a «guessed» solution spot.
Value of parameter
Local
minimum
Local
minimum
Global minimum
ECD applications
Epileptic patient:
• EEG recorded from 21 electrodes
• FDG-PET with electrode markers
EEG power
EEG data
ECD applications, cont’d
PET scan used to:
- generate a
pseudo-sphere
model
- define the
electrodes
location
ECD applications, cont’d
First peak, above F4
Second and third peak
ECD applications, cont’d
Second
source
First source
Third
source
Hypometabolic region
ECD limitations
• How many dipoles ? The more sources, the
better the fit… in a mathematical sense !!!
• Is a dipole, i.e. a punctual source, the right
model for a patch of activated cortex ?
• What about the influence of the noise ? Find
the confidence interval.
• (Is the sECD a good approach ? Given that you
find what you put in.)
Agenda
Part II : Solving the inverse problem
• Introduction
• Equivalent current dipole solution
• Distributed linear solution
• Other solutions
• Conclusion
Distributed Linear solution
With the DL solution :
“All” possible (fixed) source locations (>103)
considered simultaneously
 largely under-determined but linear problem :
v  Lj 
 external constraints required to calculate a
unique solution
Problem : What should be the constraints ?
Advantage : three-dimensional voxelwise results.
Constraints and priors
“Hard” (anatomical) :
•Orientation ( cortical
sheet)
•Location (in gray matter)
•Local spatial coherence
included in lead field
matrix L : v = L j
spatially informed basis
functions (sIBF) Bs:
j=Bs ks
“Soft” (functional or other probabilistic) :
•Other source priors
weighting matrix H :
reduced weight at priors
Example of “Soft” priors:
- location prior as found by other modalities (PET/fMRI), or
defined “by hand”
- long/short distance spatial coherence between/inside areas
- depth constraints
Weighted Minimum Norm
solution
The linear instantaneous problem:
v=Lj+
with var( ) =
C
can be solved using a Weighted Minimum Norm
(WMN) approach:
ˆj  arg min
C
1 / 2

Lj  v  
2

2
2
2
 H 1 j   H 2 j   H 3 j  ...
2
1
2
2
2
3
which can also be expressed as, where Vk  H kT H k
ˆj  arg min
Constrained
solution
Lj  v 
T




C 1 Lj  v    j T 21V1  22V2  23V3  ... j
Sources likelihood
Sources priors
WMN solution, cont’d
Simplified form of equation:
• only one source constraints H is employed, and
• the corresponding hyperparameter  is simply used to take
into account the noise contained in the data
C
ˆ
j  arg min
noise variancecovariance matrix
1 / 2

Lj  v  
2
 Hj
2
balancing between
model and constraint fit
2

Noise
regularisation
With known  and C , and spatial basis function Bs
such that j=Bs ks , the solution of the WMN problem is
ˆ
j
 B s TB s v 
TB s
T
Bs

L

C  L Bs   H H Bs
1
2
T
Bs
]
1 T
Bs
L C 1
WMN & Bayesian estimate
Under Gaussian assumptions for the distributions
of j and  , the WMN solution


ˆ
j  L C L   H H
T
1
2
T
]
1 T
L C 1 v
is connected to the Bayesian estimate.
The conditional expectation (or posterior mean) of
the sources is j given by
E ( j ,v ) 
L C
T
1


]
1 1 T
1
j

1
T
L C
LC v
 C j LT LC j L  C 
]
v
where C j is the prior covariance of the sources.
C j1  2 (H T H ) or
Then j  N (0,C j ) where  1
2
2
2
C


V


V


 j
1 1
2 2
3V 3  ...
Hierarchical “parametric
empirical Bayes” approach
In the context of a 2-level hierarchical parametric
empirical Bayes (PEB) approach, the source
localisation problem
v  Lj 
can be expressed as:
C j
By defining 
C 
v  Lj  1 , 1  N (0,C  )

 j  0   2 ,  2  N (0,C j )
 µ1C1  µ2C 2  µ3C 3  ..
 1Ce1  2Ce2  3Ce3  ...
and using an EM algorithm, j , i and µi can be
jointly estimated.
Hierachical PEB, cont’d
Hierarchical PEB, cont’d
As Cj (and C) is defined as a linear combination of
users defined covariance basis functions, the true
source (and noise) covariance matrix can be more
precisely approximated.
Interpretation:
With the 2-level approach, the unknown parameters j
are assumed to have zero mean (due to the
shrinkage priors at the 2nd level) and some
variance Cj.
To render some location more likely to be active, the
local variance can be increased: the activity at
some location is less constrained (larger variance),
thus more likely to be different from zero.
Restricted Maximum
Likelihood solution
The 2-level hierachical model can be collapsed into:
v  1  L 2  L0  1  L2
Because of the 2nd level shrinkage prior, there are
no “fixed effects” and only “random effects”.
The covariance partitioning implied is
 
E vv t  C v
 C   LC j Lt

t


C


LC
L
 i ei  i i 
i
i
The only unkowns are the hyperparameters i and i,
which can be estimated through a ReML procedure.
Somatosensory ERP
Simple right wrist electrical stimulation
P45
59 electrodes, ERP data
Stim
N20
ERP at electrode CP5
Somatosensory ERP
ERP,, cont’d
cont ’d
Somatosensory
Scalp mapping of the average N20 peak
nose
right ear
left ear
1
0
-1
Somatosensory ERP, cont’d
Reconstruction of N20
without (left) and with
(below) location prior.
Location prior
Advantages of DL solutions
Hierachical PEB and ReML (EM) seem a
good approach to EEG source localisation
as:
• LE and RMSE greatly reduced by the
introduction of accurate location priors
• accurate and inaccurate location priors used
simultaneously  no effect of inaccurate
location priors on the source reconstruction
• accurate noise variance estimate, on average
Interestingly, the seeded ECD approach is simply
an over constrained DL solution.
Limitations of DL solutions
The solution calculated totally relies on the
constraints employed, therefore
• if only inaccurate priors are used
 wrong/biased solution
• in general, without location priors
 distributed/blurred solution
 spurious sources, depending on the
threshold used
Using a single constraint is not advisable...
Agenda
Part II : Solving the inverse problem
• Introduction
• Equivalent current dipole solution
• Distributed linear solution
• Other solutions
• Conclusions
Conclusions
Depending on
• the data available,
• your knowledge about the activation and
• the assumptions you can/want to make
you should decide wich approach or
solution to use.
There does NOT exist any ideal
and general solution so far…
Special thanks to:
• Prof. Karl J. Firston (a)
• Prof. Mick D. Rugg (b)
• Dr. Pierre Maquet (c)
• Dr. Jeremie Mattout (a)
• Dr. Olivier David (a)
(a) Wellcome Department of Imaging Neuroscience,
University College London, UK.
(b) University of California Irvine, USA.
(c) Cyclotron Research Centre, University of Liège,
Belgium.