V. Volkov 1 , A. Zherdetsky 1 , S. Turovets 2 , Allen D. Malony 3

1. A 3D Vector-Additive Iterative Solver for the Anisotropic Inhomogeneous Poisson Equation in the Forward EEG problem V. Volkov1, A. Zherdetsky1, S. Turovets2, Allen D. Malony3 Department of Mathematics and Mechanics1Belarusian State University, Minsk, Belarus Neuroinformatics Center2 Department of Computer & Information Science3 University of Oregon ICCS 2009

2. Background: Observing Dynamic Brain Function • Brain activity occurs in cortex • Observing brain activity requires • high temporal and spatial resolution • Cortex activity generates scalp EEG • EEG data (dense-array, 256 channels) • High temporal (1msec) / poor spatial resolution (2D) • MR imaging (fMRI, PET) • Good spatial (3D) / poor temporal resolution (~1.0 sec) • Want both high temporal and spatial resolution • Need to solve source localization problem!!! • Find cortical sources for measured EEG signals

3. Computational Head Models • Source localization requires modeling • Goal: Full physics modeling of human head electromagnetics • Step 1: Head tissue segmentation • Obtain accurate tissue geometries • Step 2: Numerical forward solution • 3D numerical head model • Map current sources to scalp potential • Step 3: Conductivity modeling • Inject currents and measure response • Find accurate tissue conductivities • Step 4: Source localization • Applies to optical transport modeling optical (NIR) electrical

4. Source Localization • Mapping of scalp potentials to cortical generators • Signal decomposition (addressing superposition) • Anatomical source modeling (localization) • Source modeling • Anatomical Constraints • Accurate head model and physics • Computational head model formulation • Mathematical Constraints • Criteria (e.g., “smoothness”) to constrain solution • Current solutions limited by • Simplistic geometry • Assumptions of conductivities, homogeneity, isotropy

5. Theoretical and Computational Modeling image Governing Equations ICS/BCS Continuous Solutions  Finite-DifferenceFinite-ElementBoundary-ElementFinite-VolumeSpectral Discretization System of Algebraic Equations Discrete Nodal Values mesh TridiagonalADISORGauss-SeidelGaussian elimination  Equation (Matrix) Solver Approximate Solution  (x,y,z,t)J (x,y,z,t)B (x,y,z,t) solution

6. Governing Equations (Forward Problem) • Given positions and magnitudes of current sources • Given geometry and head volume  • Calculate distribution of electrical potential on scalp  • Solve linear Poisson equation on  in  withwith no-flux Neumann boundary condition  (U)=SJ, in  (U)  n = 0, on  (x,y,z) = head tissues conductivity (known) SJ = is the current source (known) U =U( x,y,z,t) is the electrical potential (to find)

7. Governing Equations (Inverse Problem) • Inverse problem uses general tomographic structure • Given distribution of the head tissue parameters  • Predict measurements values Up given a forward model F, as nonlinear functional Up=F() • Then an appropriate cost function is defined and minimized against the measurement set V: • Find global minimum using non-linear optimization

8. Modeling Anisotropy • According to Ohm’s law the current density, J, and electrical field, E, are related by • J is in the same direction as E when σ is a scalar • If σ is a tensor (anisotropic), the direction of the current density, J, is different from the direction of the applied electrical field, E:

9. Inhomogeneity and Anisotropy in Human Head • Inhomogeneous • Conductivity depends on location • Anisotropy • Conductivity depends on orientation • Human head tissues are inhomogeneous • White matter (WM): includes fiber tracts • Gray matter (GM): cortex mainly • Cerebrospinal fluid (CSF): clear conductive fluid • Skull: highly resistive, different components • Scalp • Image segmentation used to identify head tissues • Conductivities can not be directly measured accurately

10. t r Skull and Brain Anisotropy Parameterization • Skull is more conductive tangentially than radially • Diffusion is preferential along white matter tracts MRI DT brain map (Tuch et al, 2001) • Linear relation between conductivity and diffusion tensor eigenvalues •   = K (d - d0), λ= 1, 2 , 3

11. Modeling Head Electromagnetics • Current forward problem (isotropic Poisson equation) • Used in Salman et al., ICCS 2005 / 2007 • Anisotropic Forward Problem • If we model anisotropy with existing principal axes then the tensor is symmetrical - 6 independent terms: ij = ji • Numerical implementation so far deals with the orthotropic case: ii are different , but all other components of ij = 0, ij  (U)=SJ in  ,  - scalar function of (x,y.z) (U)  n = 0 on   (ijU)=SJin  , ij - tensor function of (x,y.z) (U)  n = 0 on 

12. Transformation to the Global XYZ • We assume that the conductivity tensor is diagonal in the local coordinate system (for every voxel) • The transformation from the local to the global Cartesian system for any voxel j: σjglobal=RTj σjlocal Rj , where rotation matrix Rj is defined by the local Euler angles αβγ, (sine: s and cosine: c) :

13. Anisotropic Poisson Equation • In the global Cartesian coordinate system, the anisotropic Poisson equation is expressed as • If we model anisotropy with the existing principal axes the tensor is symmetrical - 6 independent terms: ij = ji

14. Finite Difference Modeling of Poisson Equation • Finite difference approximation of the second order accuracy with mixed derivatives can be made with a minimal stencil of 7 points in • 2D (7 point stencil) [Volkov, Diff. Equations, 1997] • Generalization to 3D leads to a13-point stencil[Volkov, ICCS 2009] • The whole problemcomputational domain isrepresented by a 3Dcheckerboard lending itselffor domain decomposition(partitioning)

15. Numerical Equations • Consider even (or only odd) mesh cells, each of them having eight neighboring computational cells • Every internal node of this checkerboard grid belongs simultaneously to two neighboring cells • Natural to introduce two components of an approximate numerical solution, ( , ), where m=1, …, 8 • In these notations, the finite difference approximation, L, of the differential operator in the Poisson equation in an arbitrary node of the grid can be represented as • Am are vectors with components given by coefficients of the finite difference approximation

16. Expression for Operators A1 and A8

17. Iterative Solution • Numerical scheme is equivalent to a system of finite difference equations with a 13 diagonal system matrix and dimension N3 , where N is a total number of voxels • Elementary per-voxel step of the iterative process solves a system of linear algebraic equations • Computational complexity per iteration is Q=NQ0 /8, where Q0 is the computational cost for solving the linear system with a matrix 8 X 8 • N/8 is a number of computational cells in the checkerboard discretization • Assuming Gaussian elimination algorithm, Q0 ~ (2/3)83 ≈341 floating operations per–cell at one iteration

18. 3D Anisotropic Simulations (88x128x128 voxels) All tissues isotropic Anisotropic skull σr/σt= 1/10

19. Matlab Prototype Performance (1) Quasi-Minimal Residual (QMR) (2) BiConjugate Gradient (BiCG) (3) Vector-additive method (a,b,c) preconditioners (none, Jacobi, incomplete Cholesky) Vector-additive method not optimized Heterogeneouscoefficients 1e-04 accuracy

20. Summary • 3D finite volume algorithm for solving the anisotropic heterogenious Poisson equation based on the vector-additive implicit methods with a 13-points stencil • Variable iterative parameters to improve the convergence rate in the heterogeneous case • First attempt to implement the vector-additive numerical scheme for a 3D anisotropic problem • Believe the 3D vector additive method has better parallelism potential due to its cell-level data decomposition, especially as head volumes scale to 2563 voxels • Currently developing GPGPU implementation