# Sparsity enables estimation of both subcortical and cortical activity from MEG and EEG

1. aAthinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Charlestown, MA 02129;
2. bHarvard–Massachusetts Institute of Technology Division of Health Sciences and Technology, Cambridge, MA 02139;
3. cInstitute for Infocomm Research, Agency for Science, Technology and Research (A*STAR), Singapore 138632, Singapore;
4. dDepartment of Anesthesia, Critical Care, and Pain Medicine, Massachusetts General Hospital, Boston, MA 02114;
5. eHarvard Medical School, Boston, MA 02115;
6. fDepartment of Neurology, Massachusetts General Hospital, Charlestown, MA 02129;
7. gDepartment of Electrical & Computer Engineering, University of Maryland, College Park, MD 20742;
8. hDepartment of Neuroscience and Biomedical Engineering, Aalto University School of Science, Espoo 02150, Finland;
9. iThe Swedish National Facility for Magnetoencephalography (NatMEG), Department of Clinical Neuroscience, Karolinska Institute, Stockholm 17177, Sweden
1. Edited by Robert Desimone, Massachusetts Institute of Technology, Cambridge, MA, and approved September 18, 2017 (received for review March 31, 2017)

1. View larger version:
Fig. 2.

Fields generated by subcortical sources can be explained by currents on the cortex. (A and B) Example of unit source current in left ventroposterolateral (VPL) thalamus (sized <mml:math><mml:mrow><mml:mn>1.5</mml:mn><mml:msup><mml:mtext>cm</mml:mtext><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:math>1.5cm3) and the corresponding noiseless MEG field pattern. (C and D) Distribution of currents on cortical surface patches (sized <mml:math><mml:mrow><mml:mn>650</mml:mn><mml:msup><mml:mtext>mm</mml:mtext><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:math>650mm2) that reproduce the MEG field pattern generated by the subcortical source. C and D show the fitted cortical currents and MEG field pattern. The source current plots in A vs. C are the resultant currents from dipoles within a subdivision. The field maps in B and D are normalized. The fitted cortical field pattern is indistinguishable from the simulated subcortical pattern. This analysis illustrates how a subcortically generated field can be explained by some distribution of cortical currents.

2. View larger version:
Fig. 3.

Sparsity makes it possible to distinguish fields from subcortical and cortical sources. (A) Sources of activity derived from stimulation of the right median nerve: left VPL thalamus, primary and secondary somatosensory areas (S1, bilateral S2), and posterior parietal cortex (PPC). (B) Normalized histogram of principal angles, which quantify the correlation between fields arising from all representative combinations of activity within this neurophysiological source space. The orange histogram shows the distribution of subcortical vs. cortical angles, while the green histogram shows the distribution of cortico-cortical angles. (C) MEG field pattern resulting from activity in VPL thalamus. (D) Representative fields from example cortical source sets, whose gain matrices have the indicated principal angles with the subcortical gain matrix, that best fit the subcortical field pattern in C. All field map color scales are normalized to emphasize the spatial patterns. The spatial profiles of the cortical and subcortical MEG field patterns are distinct, even for a principal angle of <mml:math><mml:msup><mml:mn>30</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>30°, and substantially so for principal angles >45°. These distinctions suggest the feasibility of resolving simultaneous subcortical and cortical activity.

3. View larger version:
Fig. 4.

An analysis of how sparsity and hierarchy influence subcortical source estimation. (A) Illustration of all brain divisions considered. (B) Minimum-norm estimator (MNE) resolution matrix for the source space in A. (C) Summary dispersion and error metrics for the resolution matrix in B. Cortical estimates concentrate around the diagonal (low localization error), whereas subcortical estimates spread significantly to the cortex (high spatial dispersion). (D) A reduced space composed of sparse cortical regions that generate somatosensory evoked potentials combined with all subcortical volumes. (E and F) MNE resolution matrix and associated performance metrics for the reduced source space in D. The sparse subset of the cortical source space allows subcortical activity to be estimated, albeit with significant spread to nondiagonal regions. (G) Final sparse cortical and subcortical source regions identified using an inverse solution employing sparsity constraints. The faded subcortical regions show the hierarchically reduced subcortical source space, while the foreground subcortical regions show estimated sources in the thalamus. (H and I) Empirical resolution matrix (one active source per column) and associated performance metrics for the sparse solution. Estimates mostly concentrate on and around the diagonal for both cortical and subcortical sources. B, E, and H show left/right (l/r) cortex (l/rco), hippocampus (r/lh), amygdala (r/la), putamen (r/lp), caudate (r/lc), thalamus (r/lt), and brainstem (bs). All resolution matrices order sources based on physical proximity. Therefore, when sources are estimated accurately, the resolution matrix has a diagonal appearance. The blue boxes are used to delineate the position of the cortical, left thalamic, and right caudate sources in the resolution matrices. The changes in the color-scale range highlight the 3?10× increase in recovered source amplitude when sparse estimation is applied across progressively refined hierarchies. Overall, hierarchical sparsity enables focal spatial resolution with minimal dispersion (or point spread) for inverse solutions incorporating both cortical and subcortical sources.

4. View larger version:
Fig. 5.

Sparse hierarchical estimates recover simulated somatosensory responses. (A and B) Spatial distribution and time courses (one color per channel) of simulated MEG fields in sensor space. (C and D) Spatial distribution and time courses of simulated source currents in source space. Inflated views show sources located in somatosensory (S1, S2) and parietal (PPC) cortices. MRI views show thalamic source locations. The sagittal section passes through the left thalamus. The somatosensory thalamus (Th) is activated in a periodic on/off pattern. (D vs. B) While the cortical source currents contribute large-amplitude MEG signals, fields due to thalamic sources are not visible above the simulated noise. (E and F) Spatial distribution and time courses of estimated source currents in the source space. All topographical snapshots are at <mml:math><mml:mrow><mml:mn>84</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>84ms (top gray arrows in time-course plots), the color scales and slice locations are the same in C–E, and all source currents are plotted in terms of the resultant magnitudes across dipoles within each region in the key. (E and F vs. C and D) Estimated source locations and time courses closely match the simulated ground truth. The thalamic source estimate follows the true phasic on/off pattern. Although there is a stray source estimate in right thalamus, it is weak and relatively constant over time. In SI Appendix, Figs. S6 and S7 compare the performance of our algorithm to alternatives that do not use sparsity and hierarchy.

5. View larger version:
Fig. 6.

Cortical and subcortical source estimates for evoked auditory responses. Stimulus-locked average auditory evoked responses were recorded from a healthy volunteer presented with a broadband click train stimulus. Time courses are averages across <mml:math><mml:mrow><mml:mn>11,170</mml:mn></mml:mrow></mml:math>11,170 epochs filtered between 500 Hz and 1,625 Hz for the auditory brainstem response (ABR), and between 30 Hz and 300 Hz for the middle latency response (MLR). (A) MLR time courses displayed across channels, one color per channel. Red labels denote common peaks occurring at the expected poststimulus latencies. The Na and Pa peaks (shaded gray section) are particularly prominent. (B) ABR time courses rectified and averaged across channels. The shaded gray section marks the 5.0- to 6.5-ms period poststimulus, when peaks consistent with ABR wave V appear in the recordings. (C and E) Sparse cortical estimates for middle-latency recordings (<mml:math><mml:mrow><mml:mn>30</mml:mn><mml:mo>?</mml:mo><mml:mn>300</mml:mn><mml:mtext>Hz</mml:mtext><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>30?300Hz): snapshots at <mml:math><mml:mrow><mml:mn>25</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>25ms (top black arrow, E). The activity is localized to the Heschl’s gyrus and superior temporal gyrus, consistent with auditory cortical processing. The source time courses from these areas have peaks consistent with the Na and Pa peaks in the scalp recordings in A. (D and F) Sparse hierarchical estimates for early-latency recordings (<mml:math><mml:mrow><mml:mn>500</mml:mn><mml:mo>?</mml:mo><mml:mn>1,625</mml:mn><mml:mtext>Hz</mml:mtext><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:math>500?1,625Hz), obtained using a source space comprised of sparse subsets of cortex in C and the distributed subcortical space. The spatial plots display the source activity at <mml:math><mml:mrow><mml:mn>5</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>5ms (top black arrow, F). The activity is localized primarily to the inferior colliculus. A weak stray source is also seen in right amygdala. The brainstem source time courses show peaks consistent with the ABR wave V peaks in B. The color scales and slice locations are maintained for topographical snapshots across C and D, and all source currents are resultant magnitudes across dipoles within each region in the key. The color scales in C and D have units nAm and <mml:math><mml:mrow><mml:mn>0.1</mml:mn><mml:mtext>nAm</mml:mtext></mml:mrow></mml:math>0.1nAm, respectively. Overall, our sparse hierarchical algorithm recovers cortical and subcortical sources consistent with the auditory stimuli presented. SI Appendix, Fig. S10 compares the performance of our algorithm to alternatives that do not use sparsity and hierarchy.

#### Online Impact

• 8686301327 2018-02-22
• 1879481326 2018-02-22
• 9332351325 2018-02-22
• 7384141324 2018-02-22
• 8918371323 2018-02-22
• 7638311322 2018-02-22
• 9654151321 2018-02-22
• 1588961320 2018-02-22
• 5712971319 2018-02-22
• 5536211318 2018-02-22
• 4417061317 2018-02-22
• 3024201316 2018-02-21
• 4658931315 2018-02-21
• 3216561314 2018-02-21
• 1965251313 2018-02-21
• 970811312 2018-02-21
• 609011311 2018-02-21
• 3219131310 2018-02-21
• 613261309 2018-02-21
• 6972481308 2018-02-21