EEG Source Estimation via Cortical Graph Smoothing

Introduction

The EEG source estimation problem consists of inferring brain activity from electrical potentials measured on the surface of the scalp. This is a challenging inverse problem as the number of degrees of freedom of activity inside the brain greatly exceeds the number of available measurements (e.g. scalp electrodes). Within the past ten years, the development of diffusion weighted imaging and associated modalities (e.g. DTI, DSI, HARDI, etc) has enabled non-invasive measurement of the anatomical connections between brain regions of individual subjects.  In this work, we seek to exploit the availability of connectivity information in order to improve EEG source localization. We describe the cortical connectivity using a symmetric weighted graph, and introduce a prior for source activity by penalizing the sums of squares of differences of cortical activity across the edges. This encourages solutions which have consistent activation across connected cortical patches.

Methods

This work is done in the context of ongoing work on building accurate numerical models of human head electromagnetics, for EEG and ERP analysis [1]. Subject specific models are based on head tissue geometry segmented from a T1 MRI image, using the BrainK software package developed at the NeuroInformatics Center. The cortical surface is computed as the outer surface of segmented gray matter, and partitioned into N_d=2400 roughly equal area patches. Cortical sources are modeled by placing current dipoles at the center of each of the cortical patches, oriented perpendicular to the cortical surface. The relation between source currents J ∈ R^N_d and measured potentials Φ ∈ R^N_e  at the electrodes, in the noise-free case, is given by Φ = KJ, where the columns of the N_e×N_d  lead-field matrix K are determined by solving the inhomogeneous Poisson equation describing the electromagnetics of volume conduction.

We compute the connectome matrix by aligning the partitioned cortical surface with tract streamlines computed by whole-brain tractography, using software developed at the Computational Radiology Laboratory. The raw data for tractography consists of 10 unweighted (b=0) and 60 directional diffusion weighted images (b=700 s/mm^2), which are used to compute the diffusion tensor image. The tractography algorithm employed differs from simple streamlining through the use of both tensor deflection and directional inertia, which are employed to encourage tract streamlines to pass through regions of reduced anisotropy caused by crossing fibers. Additionally, sub-voxel resolution of the tract trajectory is enabled by use of log-Euclidean tensor interpolation. Tracts are initialized by seeding 30 random points in all image voxels with fractional anisotropy exceeding 0.4, yielding approximately 15 million streamlines which were then thresholded by retaining only tracts with both endpoints within 1 cm of the cortical surface. Start and end regions for each tract in this reduced set were determined from the cortical patches closest to the tract endpoints. The i,j entry of the N_d×N_d tractography-based connectome matrix A^tr is then given by summing the reciprocals of tract lengths, for all tracts connecting regions i and j.

In addition to tractography, which captures the long-range non-local connections between cortical regions, we also include a local component of the connectome based purely on spatial adjacency of cortical patches. This local component is denoted A^loc, given by a^loc_i,j=1 if cortical patches i and j are adjacent on the cortical surface, or 0 otherwise. 

The overall weighted graph connectome has adjacency matrix A=τ^tr A^tr + τ^loc A^loc, where τ^tr and τ^loc are regularization constants. The graph smoothing penalty is given by ∑_i,ja_i,j(J_i −J_j)².

This expression can be written as J^T L J, using the graph Laplacian operator L defined by L=D-A, where D is the diagonal degree matrix  with nonzero entries d_i,i=∑_j a_i,j

The cortical graph smoothing (CGS) estimate for the source currents, given observed potentials Phi at some timepoint, is defined as the solution to