• PNAS Physics Portal
  • Science Sessions: The PNAS Podcast Program

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

  1. Patrick L. Purdond,e,1,2
  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)

Significance

Subcortical structures play a critical role in brain functions such as sensory perception, memory, emotion, and consciousness. There are limited options for assessing neuronal dynamics within subcortical structures in humans. Magnetoencephalography and electroencephalography can measure electromagnetic fields generated by subcortical activity. But localizing the sources of these fields is very difficult, because the fields generated by subcortical structures are small and cannot be distinguished from distributed cortical activity. We show that cortical and subcortical fields can be distinguished if the cortical sources are sparse. We then describe an algorithm that uses sparsity in a hierarchical fashion to jointly localize cortical and subcortical sources. Our work offers alternative perspectives and tools for assessing subcortical brain dynamics in humans.

Abstract

Subcortical structures play a critical role in brain function. However, options for assessing electrophysiological activity in these structures are limited. Electromagnetic fields generated by neuronal activity in subcortical structures can be recorded noninvasively, using magnetoencephalography (MEG) and electroencephalography (EEG). However, these subcortical signals are much weaker than those generated by cortical activity. In addition, we show here that it is difficult to resolve subcortical sources because distributed cortical activity can explain the MEG and EEG patterns generated by deep sources. We then demonstrate that if the cortical activity is spatially sparse, both cortical and subcortical sources can be resolved with M/EEG. Building on this insight, we develop a hierarchical sparse inverse solution for M/EEG. We assess the performance of this algorithm on realistic simulations and auditory evoked response data, and show that thalamic and brainstem sources can be correctly estimated in the presence of cortical activity. Our work provides alternative perspectives and tools for characterizing electrophysiological activity in subcortical structures in the human brain.

Deep brain structures play important roles in brain function. For example, brainstem and thalamic relay nuclei have a central role in sensory processing (1, 2). Thalamocortical and hippocampal oscillations govern states of sleep, arousal, and anesthesia (3????8). Striatal regions are crucial for movement planning, while limbic structures like the hippocampus and amygdala drive memory, emotion, and learning (9???13). Altered signaling within the thalamus, striatum, hippocampus, and amygdala underlies pathologies such as autism, dementia, and depression (14). Much of our understanding of subcortical function comes from lesion studies and invasive electrophysiological recordings in animal models. Improved tools to characterize subcortical activity in humans would make it possible to analyze interactions between subcortical structures and other brain areas and could be used to better understand how subcortical activity relates to perception, cognition, behavior, and associated disorders.

At present, techniques for characterizing deep brain dynamics are limited. Invasive electrophysiological recordings (15) in humans are generally limited to regions that need to be monitored for clinical purposes. Functional magnetic resonance imaging (fMRI) can noninvasively measure deep brain activity with a spatial resolution of up to <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>2</mml:mn><mml:mtext>mm</mml:mtext></mml:mrow></mml:mrow></mml:mrow></mml:math>~1?2mm, but cannot record fast signals or oscillations. Magnetoencephalography (MEG) and electroencephalography (EEG) noninvasively measure fields generated by neural currents with millisecond-scale temporal resolution (16), and M/EEG source estimation (17) is widely used to spatially resolve neural dynamics to within <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>0.5</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>2</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:mrow></mml:mrow></mml:math>~0.5?2cm in the cerebral cortex (18??21). However, it remains an open question whether M/EEG can be used to estimate neural currents in deep brain structures.

The anatomy of deep brain structures poses significant challenges for source estimation with M/EEG. Subcortical structures produce smaller scalp M/EEG signals because subcortical structures are farther from the scalp than cortical structures. In addition, neurons in subcortical structures can have a closed-field geometry that further weakens the fields observed at a distance (22). A second, perhaps more fundamental problem is that subcortical structures are enclosed by the cortical mantle. Thus, measurements arising from activity within deep brain structures can potentially be explained by a surrogate distribution of currents on the cortical surface. This implies that it would be very difficult to estimate subcortical activity when cortical activity is occurring simultaneously.

Clearly, if the scalp M/EEG signals generated by subcortical structures are too small to measure, it is not possible to estimate subcortical sources. But in many cases, subcortical sources do generate measurable scalp signals, and in these cases it might be possible to estimate those sources, if only their fields could be distinguished from those generated by cortical sources. In many neuroscience studies, salient cortical activity can be restricted to a set of well-circumscribed areas (23, 24). We introduce the hypothesis that if the cortical generators are sparse in this sense, it might be possible to uniquely recover both the cortical and subcortical sources.

In this paper, we analyze M/EEG field patterns generated by cortical and subcortical sources and assess the extent to which sparse cortical and subcortical sources can be distinguished. We then introduce a hierarchical sparse estimation algorithm to illustrate how we can capitalize on these concepts to characterize both cortical and subcortical activity. We demonstrate the algorithm’s performance on simulated and experimental M/EEG data containing both cortical and subcortical activity.

Theory

Neural Sources and M/EEG Fields.

Primary neural currents (25), usually modeled with an ensemble of current dipoles, generate M/EEG fields <mml:math><mml:mi>??</mml:mi></mml:math>??,<mml:math display="block"><mml:mrow><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>M</mml:mi></mml:mrow></mml:msub><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>M</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub></mml:mrow><mml:mo>+</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>??N×T=??N×M??M×T+??N×T,[1]where <mml:math><mml:mi>??</mml:mi></mml:math>?? is the gain matrix determined by the quasistatic approximation of Maxwell’s equations, <mml:math><mml:mi>??</mml:mi></mml:math>?? contains the amplitudes of the current dipole sources, <mml:math><mml:mi>??</mml:mi></mml:math>?? is the noise, <mml:math><mml:mi>N</mml:mi></mml:math>N is the number of sensors, <mml:math><mml:mi>M</mml:mi></mml:math>M is the number of sources, and <mml:math><mml:mi>T</mml:mi></mml:math>T is the number of time points measured. To simplify notation, we assume that the data, the gain matrix, and the noise in Eq. 1 have been whitened to account for the spatial covariance <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>N</mml:mi></mml:mrow></mml:msub></mml:math>??N×N of the actual observation noise (26), so that <mml:math><mml:mi>??</mml:mi></mml:math>?? is Gaussian with zero mean and identity covariance matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>N</mml:mi></mml:mrow></mml:msub></mml:math>??N×N. When <mml:math><mml:mrow><mml:mi>T</mml:mi><mml:mo>=</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:math>T=1, we use notations <mml:math><mml:mi>??</mml:mi></mml:math>?? and <mml:math><mml:mi>??</mml:mi></mml:math>?? in lieu of <mml:math><mml:mi>??</mml:mi></mml:math>?? and <mml:math><mml:mi>??</mml:mi></mml:math>??.

We use high-resolution structural MRIs from healthy volunteers to delineate cortical surfaces and subcortical anatomic regions to define the locations and orientations of the elementary dipole sources (Materials and Methods). We place sources on cortical surfaces and in subcortical volumes and cluster proximal groups of dipoles into surface patches or volume subdivisions sized to homogenize signal strengths (Materials and Methods; Fig. 1 A and B; and SI Appendix, section SI I and Fig. S1). The resulting set <mml:math><mml:mi mathvariant="script">B</mml:mi></mml:math>B of <mml:math><mml:mi>K</mml:mi></mml:math>K patches and subdivisions, together called divisions, constitutes the distributed brain source space.

We can then group the columns of <mml:math><mml:mi>??</mml:mi></mml:math>?? and rows of <mml:math><mml:mi>??</mml:mi></mml:math>?? according to these divisions and rewrite Eq. 1 as<mml:math display="block"><mml:mrow><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:munderover><mml:mo largeop="true" movablelimits="false" symmetric="true">∑</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:mo>=</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mi>K</mml:mi></mml:munderover><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:mrow></mml:mrow><mml:mo>+</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>??N×T=∑k=1K??k??k+??N×T,[2]where <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k and <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k denote the gain matrix and source currents within the <mml:math><mml:mi>k</mml:mi></mml:math>kth division, respectively. We compute <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k for each division <mml:math><mml:mi>k</mml:mi></mml:math>k in <mml:math><mml:mi mathvariant="script">B</mml:mi></mml:math>B and decompose it into a set of eigenmodes (Materials and Methods). We denote gain matrices and source currents for a set of divisions <mml:math><mml:mrow><mml:mi mathvariant="script">F</mml:mi><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi></mml:mrow></mml:math>F?B by <mml:math><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow><mml:mtext>and</mml:mtext><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:mi>k</mml:mi></mml:mrow></mml:mrow><mml:mo>∈</mml:mo><mml:mi mathvariant="script">F</mml:mi></mml:mrow></mml:math>??F={??k}and??F={??k}?k∈F, respectively. Fig. 1 C and D illustrates field patterns for one cortical and one subcortical division.

Fields Generated by Subcortical Sources Can Be Explained by Currents on the Cortex.

To analyze distinctions between subcortical and cortical fields, we investigated the extent to which MEG field patterns arising from subcortical currents can be explained by cortical surrogates.

We simulated field pattern <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">V</mml:mi></mml:msub></mml:math>??V corresponding to unit current in the VPL thalamus (Fig. 2 A and B) and assessed whether <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">V</mml:mi></mml:msub></mml:math>??V could be explained by some distribution of cortical source currents. Specifically, we fitted the subcortical field pattern with cortical sources; i.e., we computed the cortical minimum-norm estimate to explain <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">V</mml:mi></mml:msub></mml:math>??V. We found that the resulting currents are small and broadly distributed across several cortical patches (Fig. 2C). Further, the goodness of fit between the field pattern for the cortical estimate (Fig. 2D) and <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">V</mml:mi></mml:msub></mml:math>??V was <mml:math><mml:mrow><mml:mn>98.5</mml:mn><mml:mo>%</mml:mo></mml:mrow></mml:math>98.5%, showing that <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">V</mml:mi></mml:msub></mml:math>??V can be explained by some combination of sources in the full dense cortical space.

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.

Analysis of Forward Solutions.

To generalize the above result, we used principal angles (27, 28) to characterize the relationship between the cortical and subcortical field patterns. Principal angles quantify the correlation between linear subspaces, in this case the space spanned by MEG fields arising from sources in different cortical and subcortical regions. For the example in Fig. 2, the maximum principal angle between subspaces spanned by the subcortical and cortical gain matrices was <mml:math><mml:msup><mml:mn>0</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>. Further, the maximum principal angle between any subcortical gain matrix and the cortical gain matrix was <mml:math><mml:msup><mml:mn>0</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>. We conclude that the presence of the full cortical source space makes it impossible to unambiguously estimate currents in simultaneously active subcortical sources.

Sparsity Makes It Possible to Distinguish Fields from Subcortical and Cortical Sources.

We next studied the extent of subspace correlation when subcortical sources are active together with only a small subset of cortex. As an example, we examined the scenario of median-nerve somatosensory stimulation, which elicits activity in ventroposterolateral (VPL) thalamus, primary and secondary sensory cortices (S1, S2), and posterior parietal cortex (PPC) (29), and analyzed forward solutions for a source space encompassing these brain areas (Fig. 3A).

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.

Specifically, we considered all possible configurations of subcortical and cortical currents in these divisions, computed principal angles between the subcortical and cortical gain matrices corresponding to each possible configuration (Materials and Methods), and plotted the distribution of angles (Fig. 3B). A large proportion of the principal angles are high (median <mml:math><mml:msup><mml:mn>43.9</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>43.9°), indicating that many different configurations of activity within the sparse cortical and subcortical divisions can be distinguished from one another. We also computed principal angles for all mutually exclusive configurations of activity within the cortical divisions in Fig. 3A and found comparable principal angles (median <mml:math><mml:msup><mml:mn>63.5</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>63.5°). Therefore, in principle, the problem of distinguishing subcortical sources from sparse cortical sources is similar in difficulty to that of distinguishing sparse cortical sources from one another. We also illustrate typical subcortical and cortical field patterns (Materials and Methods) for source current configurations corresponding to the various angles in this distribution (Fig. 3 C and D). Subcortical and cortical field patterns with principal angles as low as <mml:math><mml:msup><mml:mn>45</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>45° are clearly distinguishable.

This example represents a conservative scenario, but illustrates an approach for characterizing the extent to which subcortical sources can be resolved for any given cortical source distribution. We also found similar trends for other more general cases (SI Appendix, section SI II and Figs. S2 and S3). These results lead us to conclude that spatial sparsity constraints can enable distinctions between cortical and subcortical field patterns. Based on this analysis, we introduce and test an inverse algorithm that uses a sparse cortical representation to achieve localization of simultaneous subcortical and cortical activity.

Inverse Algorithm

Electromagnetic Inverse Problem.

The electromagnetic inverse problem is to estimate source currents <mml:math><mml:mi>??</mml:mi></mml:math>?? underlying M/EEG data <mml:math><mml:mi>??</mml:mi></mml:math>??, given forward gain matrix <mml:math><mml:mi>??</mml:mi></mml:math>?? for divisions distributed across the brain. This inverse problem is commonly solved using the linear <mml:math><mml:msub><mml:mi>l</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>l2 minimum-norm estimator (MNE),<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:mi>??</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mi>??</mml:mi></mml:mrow></mml:mrow><mml:mo>,</mml:mo></mml:mrow></mml:math>??^(??,??)=??MNE(??)??,[3]where <mml:math><mml:mrow><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:mi>??</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>??^(??,??) is an estimate of <mml:math><mml:mi>??</mml:mi></mml:math>??, and the MNE estimator <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup></mml:math>??MNE is a function of <mml:math><mml:mi>??</mml:mi></mml:math>??. The performance of <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup></mml:math>??MNE can be assessed using the resolution matrix<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:mi>??</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mi>??</mml:mi></mml:mrow></mml:mrow><mml:mo>.</mml:mo></mml:mrow></mml:math>??(??)=??MNE(??)??.[4]The ideal <mml:math><mml:mrow><mml:mi>??</mml:mi><mml:mo>=</mml:mo><mml:mi>??</mml:mi></mml:mrow></mml:math>??=?? corresponds to a <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup></mml:math>??MNE which exactly recovers the current locations and amplitudes, in the absence of noise (30). In practice, source estimates are biased and more extended than the true sources. This nonideal behavior can be analyzed using the spatial dispersion (SD) and dipole localization error (DLE) metrics (31) (Materials and Methods), which indicate how far the inverse solution for a given source spreads from the actual source location. We analyze the performance of the MNE on the subcortical and cortical source estimation problems and describe our proposed inverse algorithm.

Estimation Performance—Distributed Cortical and Subcortical Sources.

We first considered the problem of estimating neural currents in a set of divisions <mml:math><mml:mi mathvariant="script">B</mml:mi></mml:math>B distributed across both cortical and subcortical structures in the brain (Fig. 4A). Fig. 4B shows the MNE resolution matrix <mml:math><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>??MNE(??B). We found that estimates for the cortical sources concentrate around the diagonal (Fig. 4B, upper left), implying good resolution for cortical sources. On the other hand, estimates for subcortical sources have low amplitudes on the diagonal and instead spread to cortical sources (Fig. 4B, upper right and lower left). Fig. 4C shows the distribution of SD and DLE across all divisions. The median SD is <mml:math><mml:mrow><mml:mn>4.23</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>4.23cm (close to the radius of the human brain), and the median DLE is <mml:math><mml:mrow><mml:mn>1.89</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>1.89cm. The cortical DLE is ~<mml:math><mml:mrow><mml:mn>0.5</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>0.5cm and the subcortical DLE is in the 2- to 3-cm range. These findings are consistent with our principal angle analyses (Fig. 2).

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.

Estimation Performance—Sparse Cortical and Distributed Subcortical Sources.

Earlier we found the principal angles improve when only sparse subsets of cortical sources are active alongside deep sources (Fig. 3). Therefore, we assessed whether the resolution of the MNE inverse solution improves similarly. We considered the somatosensory stimulation example from Fig. 3 and constructed a composite source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br comprising the sparse cortical divisions in Fig. 3A alongside all subcortical divisions (Fig. 4D). Fig. 4E shows the MNE resolution matrix <mml:math><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>??MNE(??Br). We see that the estimates for subcortical sources do not spread significantly to the cortex (Fig. 4E, low values for upper right and lower left). However, the estimates for subcortical sources tend to spread across the subcortical source space (Fig. 4E, lower right, off-diagonal portions). Fig. 4F shows the resolution error metrics across all divisions in <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br. The median SD is <mml:math><mml:mrow><mml:mn>2.10</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>2.10cm, and the median DLE is <mml:math><mml:mrow><mml:mn>0.805</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>0.805cm. This is an improvement from the previous case, but still not as accurate as needed to resolve many subcortical sources. Given a sparse cortical source space as a starting point, we anticipate that it might be possible to use a subsequent sparse estimation step to reduce spread among subcortical sources.

Estimation Performance—Sparse Cortical and Sparse Subcortical Sources.

Sparse estimation procedures based on <mml:math><mml:msub><mml:mi>l</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>l1-norm minimization, projection pursuit, and/or Bayesian theory are effective at pruning out spurious features, while identifying relevant sparse features in noisy high-dimensional problems (32??35). We have recently shown that subspace pursuit can accurately estimate multiple sparse cortical sources underlying MEG data (36). Thus, we assessed whether a similar subspace pursuit algorithm could reduce the spread among spurious subcortical sources and enable improved resolution for subcortical source estimates.

We continue with our analysis of the resolution matrix in the composite source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br (Fig. 4G, faded background), this time with subspace pursuit. Since the subspace pursuit algorithm is nonlinear, a closed-form resolution matrix in the sense of Eq. 4 does not exist; thus the performance of subspace pursuit must be characterized empirically instead. To this end, we simulated unit currents <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:math>??i in each <mml:math><mml:mrow><mml:mi>i</mml:mi><mml:mtext>th</mml:mtext></mml:mrow></mml:math>ith division <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>i</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>Br(i) within <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br, one at a time, to generate corresponding noiseless field patterns <mml:math><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>i</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:msub><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>??i=??Br(i)??i. Then, for each <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:math>??i, we performed subspace pursuit and constructed an empirical resolution matrix <mml:math><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>SP</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>??SP(??Br) (Materials and Methods). The resulting matrix, shown in Fig. 4H, has a near-diagonal structure for the majority of cortical and subcortical sources. Moreover, Fig. 4I shows SD and DLE, across all divisions in <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br. The median SD and DLE are the same and equal to <mml:math><mml:mrow><mml:mn>0.737</mml:mn><mml:mtext>cm</mml:mtext></mml:mrow></mml:math>0.737cm. This is a substantial improvement over previous solutions that do not use sparsity constraints (Fig. 4H vs. Fig. 4 C and E).

Hierarchical Subspace-Pursuit Inverse Algorithm.

The above results suggest that it is possible to resolve both cortical and subcortical sources by applying sparsity in both domains. In previous work, we developed a sparse estimation algorithm for cortical divisions (36), in which sets of cortical divisions were nested in successively finer resolutions (i.e., smaller patches or divisions), and subspace pursuit was applied to derive sparse estimates in successively finer resolutions, which formed a hierarchy from coarse to fine resolution. We therefore intuited that subcortical sources could be viewed as an additional, ultimate step in this hierarchical refinement process, achieved by adding a set of subcortical divisions to the final set of sparse cortical sources and applying subspace pursuit.

Inverse Algorithm.

Inputs: Data <mml:math><mml:mi>??</mml:mi></mml:math>??, distributed gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub></mml:math>??B, and target sparsity level <mml:math><mml:mi>L</mml:mi></mml:math>L.

Notation: Denote the distributed cortical and subcortical source spaces as <mml:math><mml:mrow><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>C</mml:mi></mml:msub><mml:mtext>and</mml:mtext><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>S</mml:mi></mml:msub></mml:mrow><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi></mml:mrow></mml:math>BCandBS?B, respectively. Denote <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mo>?</mml:mo><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:mrow></mml:math>Hr?Br as the set of <mml:math><mml:mi>L</mml:mi></mml:math>L brain divisions whose estimated neural currents <mml:math><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub></mml:math>??^Hr best explain data <mml:math><mml:mi>??</mml:mi></mml:math>??.

  • 1. Do subspace pursuit on the distributed cortical source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>C</mml:mi></mml:msub></mml:math>BC: <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>C</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>C</mml:mi></mml:msub></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mtext>SP</mml:mtext><mml:mrow><mml:mo>(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>C</mml:mi></mml:msub></mml:msub><mml:mo>,</mml:mo><mml:mi>L</mml:mi><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>[HC,??^HC]=SP(??,??BC,L).

  • 2. Construct <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mrow><mml:mi>C</mml:mi><mml:mo>,</mml:mo><mml:mtext>refined</mml:mtext></mml:mrow></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>C</mml:mi></mml:msub><mml:mo>∪</mml:mo><mml:mrow><mml:mtext>neighbors of</mml:mtext><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi>C</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:mrow></mml:math>BC,refined=HC∪neighbors ofHC in a finer subdivision of cortical patches.

  • 3. Repeat subspace pursuit on the coarse-to-fine hierarchy of cortical source spaces <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mrow><mml:mi>C</mml:mi><mml:mo>,</mml:mo><mml:mtext>refined</mml:mtext></mml:mrow></mml:msub></mml:math>BC,refined: <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:msub><mml:mi>C</mml:mi><mml:mtext>sp</mml:mtext></mml:msub></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:msub><mml:mi>C</mml:mi><mml:mtext>sp</mml:mtext></mml:msub></mml:msub></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mtext>SP</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mrow><mml:mi>C</mml:mi><mml:mo>,</mml:mo><mml:mtext>refined</mml:mtext></mml:mrow></mml:msub></mml:msub><mml:mo>,</mml:mo><mml:mi>L</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>[HCsp,??^HCsp]=SP(??,??BC,refined,L).

  • 4. Construct the composite space of sparse cortical sources and distributed subcortical sources: <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:msub><mml:mi>C</mml:mi><mml:mtext>sp</mml:mtext></mml:msub></mml:msub><mml:mo>∪</mml:mo><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi mathvariant="script">S</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>Br=[HCsp∪BS].

  • 5. Repeat subspace pursuit on the composite sparse space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub></mml:math>Br: <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mtext>SP</mml:mtext><mml:mrow><mml:mo>(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub></mml:msub><mml:mo>,</mml:mo><mml:mrow><mml:mi>α</mml:mi><mml:mi>L</mml:mi></mml:mrow><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>[Hr,??^Hr]=SP(??,??Br,αL), where <mml:math><mml:mrow><mml:mi>α</mml:mi><mml:mo>></mml:mo></mml:mrow></mml:math>α>1.

Outputs: Cortical and subcortical source locations <mml:math><mml:mrow><mml:mi mathvariant="script">H</mml:mi><mml:mo>=</mml:mo><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub><mml:mo>?</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mn>1,2</mml:mn><mml:mo>,</mml:mo><mml:mi mathvariant="normal">…</mml:mi><mml:mo>,</mml:mo><mml:mi>K</mml:mi><mml:mo stretchy="false">]</mml:mo></mml:mrow></mml:mrow></mml:math>H=Hr?[1,2,…,K]; and the estimated time courses of neural currents at these locations <mml:math><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">H</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:msub><mml:mi mathvariant="script">H</mml:mi><mml:mi mathvariant="normal">r</mml:mi></mml:msub></mml:msub></mml:mrow></mml:math>??^H=??^Hr.

Subspace pursuit (steps 1, 3, and 5).

For a source space comprising a subset of brain divisions <mml:math><mml:mrow><mml:mi mathvariant="script">F</mml:mi><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi></mml:mrow></mml:math>F?B with gain matrices <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub></mml:math>??F, subspace pursuit (SP) estimates the locations and time courses of neural currents to explain data series <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>N</mml:mi><mml:mo>×</mml:mo><mml:mi>T</mml:mi></mml:mrow></mml:msub></mml:math>??N×T:<mml:math display="block"><mml:mrow><mml:mrow><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mi mathvariant="script">H</mml:mi><mml:mo>,</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">H</mml:mi></mml:msub><mml:mo stretchy="false">]</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mtext>SP</mml:mtext><mml:mrow><mml:mo>(</mml:mo><mml:mi>??</mml:mi><mml:mo>,</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mi>L</mml:mi><mml:mo>)</mml:mo></mml:mrow></mml:mrow></mml:mrow><mml:mo>.</mml:mo></mml:mrow></mml:math>[H,??^H]=SP(??,??F,L).[5]The pursuit procedure finds a sparse subset of dictionary <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub></mml:math>??F that best explains data <mml:math><mml:mi>??</mml:mi></mml:math>??, computes residuals remaining to be explained, and iterates to find new relatively uncorrelated subsets of <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub></mml:math>??F that best explain these residuals, until all matching subsets of <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">F</mml:mi></mml:msub></mml:math>??F are found (34?36) (SI Appendix, section SI III).

Hierarchical construction (steps 2 and 4).

We first apply SP on a distributed cortical source space to identify the subset of cortical divisions that best explain the measured fields. Then, to improve accuracy of this sparse subset of cortical sources, we perform SP across a coarse-to-fine hierarchy of cortical source spaces (36) (SI Appendix, Fig. S4). Subsequently, we construct a composite space of the sparse cortical source estimates and distributed subcortical sources and use SP on this composite sparse space to identify the subset of subcortical and cortical sources that best explain the data. This process of using sparse estimation across increasingly refined hierarchies enables systematic reduction of the distributed source space and the gain matrix: pruning out sources not important for explaining the data, implicitly decorrelating the columns of the gain matrix, and concentrating estimates into subsets of the brain whose neural currents best explain the data. Overall, given data <mml:math><mml:mi>??</mml:mi></mml:math>??, distributed gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub></mml:math>??B, and target sparsity level <mml:math><mml:mi>L</mml:mi></mml:math>L, the hierarchical SP algorithm identifies the sparse subset <mml:math><mml:mi mathvariant="script">H</mml:mi></mml:math>H that specifies locations for both cortical and subcortical sources and estimates the time courses <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">H</mml:mi></mml:msub></mml:math>??H of neural currents at those locations.

Data Examples

We illustrate the performance of the algorithm by analyzing noisy evoked response simulations and experimental data. First, we preprocess the measurements <mml:math><mml:mi>??</mml:mi></mml:math>?? and estimate the noise covariance matrix <mml:math><mml:mi>??</mml:mi></mml:math>??. Next, we use the MRI data to construct the distributed source space, i.e., brain divisions <mml:math><mml:mi mathvariant="script">B</mml:mi></mml:math>B, and compute the forward solutions <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub></mml:math>??B. Then, we use the hierarchical SP inverse algorithm to estimate the locations and time courses of sources that best explain the M/EEG data. We perform all data analyses at the individual subject level. We specify the target sparsity level <mml:math><mml:mi>L</mml:mi></mml:math>L for SP empirically, based on a conservative approximation of the expected number of active divisions. We maintain <mml:math><mml:mi>L</mml:mi></mml:math>L across the cortical hierarchies and increase it by a factor of <mml:math><mml:mrow><mml:mi>α</mml:mi><mml:mo>~</mml:mo><mml:mrow><mml:mn>1.5</mml:mn><mml:mo>?</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:mrow></mml:math>α~1.5?2 for the final hierarchy, i.e., the composite hybrid source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br comprising sparse cortical sources and distributed subcortical sources. Estimates displayed in this section are obtained at the final hierarchy level <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br. We note that <mml:math><mml:mi>L</mml:mi></mml:math>L specifies the number of eigenmodes included in the sparse solution and that currents from a given patch or subdivision may be represented by more than one eigenmode. As a result, the number of subdivisions seen on the topographical maps may be less than <mml:math><mml:mi>L</mml:mi></mml:math>L.

Somatosensory Evoked Simulations.

We simulated MEG evoked responses mimicking those elicited by electrical stimulation of the right median nerve at the wrist. Fig. 5 A and B illustrates the spatial and temporal features of simulated fields in the sensor space. The simulated fields include additive Gaussian noise, with a signal-to-noise ratio (SNR) of <mml:math><mml:mrow><mml:mn>7</mml:mn><mml:mtext>dB</mml:mtext></mml:mrow></mml:math>7dB, similar to resting eyes-open recordings. Fig. 5 C and D display the spatial and temporal patterns of simulated currents in the source space. Specifically, the evoked responses comprise early currents in the left somatosensory region of thalamus (L Som Th), followed by currents in the left primary somatosensory cortex (L S1, near the postcentral gyrus), and later currents in the left PPC (L PPC) and bilateral secondary somatosensory cortices (L/R S2). The simulated thalamic current time course has a periodic on/off pattern up to <mml:math><mml:mrow><mml:mn>250</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>250ms poststimulus. This pattern was chosen to assess source estimation performance for the challenging case of phasic, temporally overlapping subcortical and cortical activity. We note that the MEG fields arising from thalamic source currents (e.g., <mml:math><mml:mrow><mml:mn>0</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>15</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:mrow></mml:math>0?15ms in Fig. 5B) lie below the observation noise and are significantly smaller than those arising from cortical currents (e.g., <mml:math><mml:mrow><mml:mn>70</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>100</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:mrow></mml:math>70?100ms in Fig. 5B).

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.

For source estimation, we used a source space different from that used for the simulation, to recreate a scenario closer to what might occur in practice, in which the true generating sources and the source space parcellation might not correspond. We set the sparsity level <mml:math><mml:mi>L</mml:mi></mml:math>L to <mml:math><mml:mn>8</mml:mn></mml:math>8 and <mml:math><mml:mn>12</mml:mn></mml:math>12 for the cortical and hybrid hierarchies, respectively. The procedure refines source current estimates across cortical hierarchies (SI Appendix, Fig. S5) and culminates in the final hybrid hierarchy <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br. The final spatial distributions and time courses of estimated source currents (Fig. 5 E and F) closely resemble those of the simulated ground truth (Fig. 5 C and D) for both cortical and subcortical sources. The estimated left thalamic time course in Fig. 5F matches the simulation in shape and phase and further is not contaminated by leakage from cortical sources.

Further, we found that our algorithm offers significant gains in performance compared with other methods that do not use principles of sparsity and hierarchy (SI Appendix, Figs. S6 and S7). We conclude that both sparsity and hierarchy are necessary to accurately resolve locations and time courses of the thalamic source currents (SI Appendix, Fig. S6). Further, using sparsity in a hierarchical fashion helps recover the true distribution of mean source activity across anatomic regions (SI Appendix, Fig. S7).

Auditory Evoked Response Experiments.

We recorded auditory evoked responses elicited by binaural stimulation with a train of clicks (37, 38) during the resting eyes-open condition simultaneously with MEG and EEG. Auditory responses comprise distinct M/EEG peaks at established latencies corresponding to a progression of activity from the cochlea, through the inferior colliculus in the brainstem, to the auditory cortex (37), and thus serve as a suitable test case for validating a subcortical source estimation algorithm.

The M/EEG evoked responses are shown in Fig. 6 A and B. We see early auditory brainstem response (ABR) peaks in both EEG and MEG at <mml:math><mml:mrow><mml:mn>5.8</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>5.8ms and <mml:math><mml:mrow><mml:mn>6.2</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>6.2ms, a low-amplitude Po feature in the EEG at <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>10</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:mrow></mml:math>~10ms, and prominent midlatency response (MLR) peaks Na-Pa in MEG and EEG channels at <mml:math><mml:mrow><mml:mn>18</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>25</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:mrow></mml:math>18?25ms poststimulus. The ABR peaks are consistent with the brainstem wave V known to arise from the inferior colliculus (IC), the Po feature marks the end of brainstem components, and the Na-Pa peaks correspond to cortical responses known to arise in the auditory cortex (37).

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.

We performed hierarchical SP and set sparsity levels <mml:math><mml:mi>L</mml:mi></mml:math>L to <mml:math><mml:mn>4</mml:mn></mml:math>4 and <mml:math><mml:mn>8</mml:mn></mml:math>8 for the cortical and hybrid hierarchies, respectively. Fig. 6 C and E shows localization of the 18-?to?25-ms Na-Pa MLR peaks to bilateral auditory cortices and the associated time courses. The auditory areas compose the reduced cortical source space, which, along with the distributed subcortical sources, forms the hybrid source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br for estimation of deep sources underlying the ABR data. Fig. 6 D and F illustrates the localization of the wave V ABR peaks to bilateral ICs. Although the recorded wave V peaks do not have very high SNR, the source time course at the IC peaks at <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mn>5</mml:mn><mml:mo>?</mml:mo><mml:mn>6</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>~5?6ms and drops off after <mml:math><mml:mrow><mml:mn>10</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>10ms, as expected. Consistent results in another subject are shown in (SI Appendix, Fig. S9). We compared performance to algorithms that do not use sparsity and hierarchy (SI Appendix, Fig. S10) and found that hierarchical SP is necessary to estimate specific subcortical sources even for filtered recordings containing temporally separated early latency responses.

Discussion

The extent to which subcortical activity can be estimated from M/EEG measurements has been a controversial subject. Our key finding is that the M/EEG fields from subcortical sources can be distinguished from those generated by the cortex when the underlying cortical activity is sparse, i.e., limited to a subset of cortical regions. In this scenario, the problem of distinguishing subcortical from cortical sources has a similar level of ambiguity to that of resolving different cortical sources. To demonstrate how this insight might be practically applied, we developed a sparse hierarchical algorithm to estimate both sparse cortical and subcortical sources.

It is known that deep and superficial sources exhibit different M/EEG field patterns (39), but the degree to which this information could be used to resolve multiple distributed subcortical and cortical sources has remained unclear. Analyses of cortical and subcortical field patterns assuming that entire structures can be simultaneously active have provided evidence for substantial correlation (40), consistent with our Fig. 2. However, we reasoned that it would be unlikely to observe synchronous activity, the major determinant of MEG/EEG (41), simultaneously within the entirety of cortex and any given subcortical structure. Thus, we analyzed sparse subdivisions of cortical and subcortical structures and found clear distinctions in the ensuing field patterns (Fig. 3). Although previous work and the data presented here show that cortical and subcortical sources cannot, in general, be unambiguously resolved, we found that if the distribution of cortical and subcortical sources is sparse, the problem becomes tractable.

These observations motivated us to create a hierarchical SP inverse algorithm to find the set of sparse cortical and subcortical sources that best explain the M/EEG data. Our analyses of various source estimators (Fig. 4) showed that our algorithm has a performance superior to alternatives for the subcortical structures and similar to existing approaches for the cortical structures (36, 42??45).

If the locations of activity in cortical and subcortical structures are known, and each active area can be modeled with an equivalent current dipole, linear least squares can be used to estimate source current time courses (46, 47). However, fitting the locations of multiple dipolar sources usually requires tailored and often interactively guided fitting strategies. Our approach, on the other hand, automatically finds the constellation of sources in a variety of conditions, including those where the source activities may overlap in time. Further, instead of isolated current dipoles we use concise distributed dipole representations within relevant subcortical anatomical subdivisions. Other methods, such as magnetic field tomography (48) and the linearly constrained minimum-variance beamformer (49, 50), that have been applied to locate deep sources implicitly, use some of the principles we formalize here. Our analyses on simulated and experimental data show that algorithms using both sparsity and hierarchy can resolve simultaneously active cortical and subcortical sources under realistic low SNR conditions (Figs. 5 and 6).

Our algorithm relies on being able to specify, at least approximately, the sparsity level <mml:math><mml:mi>L</mml:mi></mml:math>L for the cortical and subcortical sources. This parameter can be set based on domain knowledge of the number of expected sources. In this work, we chose conservative values for <mml:math><mml:mi>L</mml:mi></mml:math>L exceeding the true or expected level of sparsity and found that the algorithm performed well. In cases where prior domain knowledge does not reasonably exist, formal model selection methods such as cross-validation could be used to choose <mml:math><mml:mi>L</mml:mi></mml:math>L. Several recent publications have proposed sparsity-based algorithms for M/EEG cortical source estimation (20, 21, 36, 51?????????61). These and other methods (62) might also be applied in a hierarchical fashion to estimate both cortical and subcortical sources. Further, techniques using distributed sparse representations and dynamical sparsity constraints could be used to estimate subcortical sources in conditions involving more extended areas of cortex (63, 64).

Our analysis implies that sparsity is not merely an assumption, but rather a crucial requirement for cortical and subcortical sources to be jointly estimated. The key question to answer is the extent to which the cortical fields can be distinguished from those of potential subcortical sources. Since the cortical areas that participate will vary according to the task, fulfillment of this requirement must be evaluated on case-by-case basis. Notably, even in complex cognitive tasks such as picture naming (23) and viewing or mimicking of lip forms (24), salient activity can be limited to a small number of regions. In each of these regions the source may be focal (well represented by a dipole) or extended (a cortical patch). The methods we have presented in this paper illustrate one approach for characterizing the extent to which subcortical activity can be recovered using M/EEG.

Materials and Methods

MRI and M/EEG Acquisition.

We acquired MRI and M/EEG data from four healthy subjects aged <mml:math><mml:mrow><mml:mn>25</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:mn>45</mml:mn><mml:mi mathvariant="normal">y</mml:mi></mml:mrow></mml:mrow></mml:math>25?45y. All subjects provided written informed consent. All studies were approved by the Human Research Committee at Massachusetts General Hospital. Data from one subject were used for the principal angles analysis and simulations. We recorded auditory evoked fields in three other subjects. In one subject, technical problems made the data unusable. For each subject, we obtained T1-weighted structural MRI, a <mml:math><mml:mn>306</mml:mn></mml:math>306-channel MEG, and a <mml:math><mml:mn>70</mml:mn></mml:math>70-electrode EEG recording (SI Appendix, section SI I).

Source Space Construction.

We used FreeSurfer to reconstruct neocortical and hippocampal surfaces and segment subcortical volumes from the MRI (65?67). We placed dipoles with orientations normal to the triangulated surface mesh for neocortex and hippocampus at the gray–white matter interface, with <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mtext>mm</mml:mtext></mml:mrow></mml:mrow></mml:math>~1mm spacing. We placed triplets of orthogonal dipoles in subcortical volumes covering the thalamus, putamen, caudate, brainstem (midbrain), and amygdala, at <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>1</mml:mn><mml:mtext>mm</mml:mtext></mml:mrow></mml:mrow></mml:math>~1mm voxel spacing. To reduce the dimensionality of the source space, we grouped neighboring cortical dipoles into “patches” (68?70). We grouped neighboring subcortical dipoles into “subdivisions” sized to produce signals with similar amplitudes to the cortical patches (SI Appendix, section SI I). For cortical patches with an average area of <mml:math><mml:mrow><mml:mn>175</mml:mn><mml:msup><mml:mtext>mm</mml:mtext><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:math>175mm2, this sizing procedure resulted in <mml:math><mml:mn>209</mml:mn></mml:math>209 subcortical subdivisions with volumes ranging from 175- to 1,800-?<mml:math><mml:msup><mml:mtext>mm</mml:mtext><mml:mn>3</mml:mn></mml:msup></mml:math>mm3 (Fig. 1B). Regions with higher current density have finer divisions (higher resolution) than those with low current density (e.g., <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mn>200</mml:mn><mml:msup><mml:mtext>-mm</mml:mtext><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:math>~200-mm3 divisions in striatum vs. <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mn>1</mml:mn><mml:mo>,</mml:mo><mml:mn>795</mml:mn><mml:msup><mml:mtext>-mm</mml:mtext><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:math>~1,795-mm3 divisions in thalamus).

Forward Solutions.

We derived a three-compartment boundary-element model from the MRI data and numerically computed forward solutions using the MNE software package (68, 69). To account for the different sensor types, units, and noise levels in the M/EEG measurements, we whitened the gain matrices, using an estimate of the observation noise covariance matrix (69). For MEG simulation studies, we set the noise covariance matrix to be similar to typical resting eyes-open recordings, <mml:math><mml:mrow><mml:mi>??</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mtext>diag</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mo stretchy="false">[</mml:mo><mml:mrow><mml:msup><mml:mi mathvariant="normal">g</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>,</mml:mo><mml:msup><mml:mi mathvariant="normal">g</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>,</mml:mo><mml:msup><mml:mi mathvariant="normal">m</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>,</mml:mo><mml:mi mathvariant="normal">…</mml:mi><mml:mo>,</mml:mo><mml:msup><mml:mi mathvariant="normal">g</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>,</mml:mo><mml:msup><mml:mi mathvariant="normal">g</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>,</mml:mo><mml:msup><mml:mi mathvariant="normal">m</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:mrow><mml:mo stretchy="false">]</mml:mo></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:math>??=diag([g2,g2,m2,…,g2,g2,m2]), where <mml:math><mml:mrow><mml:mi mathvariant="normal">g</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mrow><mml:mn>?2.5</mml:mn><mml:mtext>fT</mml:mtext></mml:mrow><mml:mo>/</mml:mo><mml:mtext>cm</mml:mtext></mml:mrow></mml:mrow></mml:math>g=?2.5fT/cm and <mml:math><mml:mrow><mml:mi>m</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mn>?10</mml:mn><mml:mtext>fT</mml:mtext></mml:mrow></mml:mrow></mml:math>m=?10fT. For M/EEG experimental studies, we estimated noise covariance matrices from the resting eyes-open data. To account for differences in current strength across brain divisions, we scaled gain matrices for each division by the regional current strengths (71) (SI Appendix, section SI I). We constructed the reduced-dimensionality M/EEG gain matrices <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k for each division <mml:math><mml:mrow><mml:mi>k</mml:mi><mml:mi>?</mml:mi><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:mi mathvariant="script">B</mml:mi><mml:mo>:</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>,</mml:mo><mml:mn>?2</mml:mn><mml:mo>,</mml:mo><mml:mi mathvariant="normal">…</mml:mi><mml:mo>,</mml:mo><mml:mi>K</mml:mi></mml:mrow><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:mrow></mml:math>k?{B:?1,?2,…,K}, using a singular value decomposition, retaining components to capture <mml:math><mml:mrow><mml:mo>></mml:mo><mml:mrow><mml:mn>95</mml:mn><mml:mo>%</mml:mo></mml:mrow></mml:mrow></mml:math>>95% of the total spectral energy (36, 70).

Analysis of Forward Solutions.

We simulated the field pattern <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k for a division <mml:math><mml:mi>k</mml:mi></mml:math>k by activating the most significant eigenmode of <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k with a unit current. We assessed the degree to which the field pattern <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k could be explained by some distribution of currents in a region <mml:math><mml:mrow><mml:mi mathvariant="script">R</mml:mi><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>k</mml:mi><mml:mo>?</mml:mo><mml:mi>R</mml:mi></mml:mrow></mml:mrow></mml:math>R?B(k?R) by fitting <mml:math><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi></mml:msub><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi></mml:msub></mml:mrow></mml:math>??R??R to <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k, using the regularized <mml:math><mml:msub><mml:mi>l</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>l2 minimum-norm criterion. We computed <mml:math><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">R</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:msubsup><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi><mml:mo>′</mml:mo></mml:msubsup><mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi></mml:msub><mml:msubsup><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi><mml:mo>′</mml:mo></mml:msubsup></mml:mrow><mml:mo>+</mml:mo><mml:mrow><mml:msup><mml:mi>λ</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mi>??</mml:mi></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:msup><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>??^R=??R′(??R??R′+λ2??)?1??k, setting the regularization parameter <mml:math><mml:msup><mml:mi>λ</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:math>λ2 to <mml:math><mml:mrow><mml:mn>1</mml:mn><mml:mo>/</mml:mo><mml:mn>9</mml:mn></mml:mrow></mml:math>1/9. We then quantified the goodness of fit between the best-fitting field patterns, <mml:math><mml:mrow><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">R</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">R</mml:mi></mml:msub><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">R</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>??^R=??R??^R and the original pattern <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub></mml:math>??k, as <mml:math><mml:mrow><mml:mn>1</mml:mn><mml:mo>?</mml:mo><mml:mrow><mml:msub><mml:mrow><mml:mo fence="true">||</mml:mo><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:mo>?</mml:mo><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">R</mml:mi></mml:msub></mml:mrow><mml:mo fence="true">||</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msub><mml:mo>/</mml:mo><mml:msub><mml:mrow><mml:mo fence="true">||</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:mo fence="true">||</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msub></mml:mrow></mml:mrow></mml:math>1?||??k???^R||2/||??k||2. For visualization of the fields, we dewhitened <mml:math><mml:msub><mml:mrow><mml:mover accent="true"><mml:mi>??</mml:mi><mml:mo stretchy="false">^</mml:mo></mml:mover></mml:mrow><mml:mi mathvariant="script">R</mml:mi></mml:msub></mml:math>??^R and mapped it to a virtual grid of <mml:math><mml:mn>304</mml:mn></mml:math>304 magnetometers distributed evenly across the Elekta helmet (72).

We used principal angles (27, 28) to quantify the correlation between putative field patterns arising from currents in two nonoverlapping regions: <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi></mml:mrow></mml:math>R1?B and <mml:math><mml:mrow><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub><mml:mo>?</mml:mo><mml:mi mathvariant="script">B</mml:mi></mml:mrow></mml:math>R2?B. We define a configuration of currents in a collection of patches or subdivisions as a combination of eigenmodes from those patches or subdivisions. Any putative field pattern arising from a configuration of currents in <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>R1 is defined by some subset <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub></mml:math>??R1 of the eigenmodes of the gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub></mml:math>??R1. If <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub></mml:math>??R1 comprises <mml:math><mml:mrow><mml:mi>n</mml:mi><mml:mn>1</mml:mn></mml:mrow></mml:math>n1 eigenmodes, there are <mml:math><mml:mrow><mml:msup><mml:mn>2</mml:mn><mml:mrow><mml:mi>n</mml:mi><mml:mn>1</mml:mn></mml:mrow></mml:msup><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:math>2n1?1 subsets of eigenmodes <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R1} collectively corresponding to all possible current configurations within <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>R1. Similarly, for region <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>R2, the gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub></mml:math>??R2 comprises <mml:math><mml:mrow><mml:mi>n</mml:mi><mml:mn>2</mml:mn></mml:mrow></mml:math>n2 eigenmodes, and there are <mml:math><mml:mrow><mml:msup><mml:mn>2</mml:mn><mml:mrow><mml:mi>n</mml:mi><mml:mn>2</mml:mn></mml:mrow></mml:msup><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:math>2n2?1 subsets of eigenmodes <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R2} collectively corresponding to all possible current configurations within <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>R2. The degree to which a field pattern arising from any current configuration in <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>R1 can be explained by a field pattern arising from any current configuration in <mml:math><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>R2 is specified by the set of principal angles <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi mathvariant="normal">Θ</mml:mi><mml:mrow><mml:mn>1</mml:mn><mml:mo>,</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{Θ1,2} between each subset of <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R1} and each subset in <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R2}. We computed the set of principal angles <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi mathvariant="normal">Θ</mml:mi><mml:mrow><mml:mn>1</mml:mn><mml:mo>,</mml:mo><mml:mn>2</mml:mn></mml:mrow></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{Θ1,2} across all <mml:math><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msup><mml:mn>2</mml:mn><mml:mrow><mml:mi>n</mml:mi><mml:mn>1</mml:mn></mml:mrow></mml:msup><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msup><mml:mn>2</mml:mn><mml:mrow><mml:mi>n</mml:mi><mml:mn>2</mml:mn></mml:mrow></mml:msup><mml:mo>?</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>(2n1?1)(2n2?1) combinations of subsets of <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R1} and <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R2}.

The principal angles analysis does not account for the gains of the eigenmodes. To ensure that the principal angles histogram is not biased toward smaller eigenmodes, we confined the histogram to include only those subsets of <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R1} and <mml:math><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:math>{??R2} that contained the largest eigenmodes from each patch or subdivision. To illustrate how the principal angles correspond to field patterns with varying levels of similarity, we selected a representative combination of eigenmodes <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub></mml:math>??R1 and <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub></mml:math>??R2 having angles <mml:math><mml:msup><mml:mn>35</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>35°, <mml:math><mml:msup><mml:mn>45</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>45°, <mml:math><mml:msup><mml:mn>60</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>60°, and <mml:math><mml:msup><mml:mn>86</mml:mn><mml:mo>°</mml:mo></mml:msup></mml:math>86°; generated the field <mml:math><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mn>1</mml:mn></mml:msub><mml:mo>=</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:msub></mml:mrow></mml:math>??1=??R1; projected (minimum <mml:math><mml:msub><mml:mi>l</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>l2-norm) <mml:math><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mn>2</mml:mn></mml:msub><mml:mo>=</mml:mo><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">R</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:msub></mml:mrow></mml:mrow></mml:math>??2=??R2??R2 to <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>??1; and plotted the <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>??2 most similar to <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mn>1</mml:mn></mml:msub></mml:math>??1.

Analysis of Resolution Matrices for Inverse Solutions.

We used the resolution matrix to assess performance of the minimum <mml:math><mml:msub><mml:mi>l</mml:mi><mml:mn>2</mml:mn></mml:msub></mml:math>l2-norm (MNE) and SP inverse solutions. The MNE solution for the distributed source space <mml:math><mml:mi mathvariant="script">B</mml:mi></mml:math>B, with prewhitened gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub></mml:math>??B, is <mml:math><mml:mrow><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:mo>)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:msubsup><mml:mi>????</mml:mi><mml:mi mathvariant="script">B</mml:mi><mml:mo>′</mml:mo></mml:msubsup><mml:msup><mml:mrow><mml:mo>(</mml:mo><mml:mrow><mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:msubsup><mml:mi>????</mml:mi><mml:mi mathvariant="script">B</mml:mi><mml:mo>′</mml:mo></mml:msubsup></mml:mrow><mml:mo>+</mml:mo><mml:mrow><mml:msup><mml:mi>λ</mml:mi><mml:mn>??</mml:mn></mml:msup><mml:mi>??</mml:mi></mml:mrow></mml:mrow><mml:mo>)</mml:mo></mml:mrow><mml:mrow><mml:mo>?</mml:mo><mml:mn>??</mml:mn></mml:mrow></mml:msup></mml:mrow></mml:mrow></mml:math>??MNE(??B)=????B′(??B????B′+λ????)???. We specified the prior source covariance <mml:math><mml:mi>??</mml:mi></mml:math>?? so that <mml:math><mml:mrow><mml:mrow><mml:mtext>Tr</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mtext>??</mml:mtext><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:msubsup><mml:mrow><mml:mi mathvariant="normal">??</mml:mi><mml:mtext>??</mml:mtext></mml:mrow><mml:mi mathvariant="script">B</mml:mi><mml:mo>′</mml:mo></mml:msubsup></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:mtext>Tr</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mtext>??</mml:mtext><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mi>N</mml:mi></mml:mrow></mml:math>Tr(??B????B′)=Tr(??)=N, set <mml:math><mml:mrow><mml:msup><mml:mi>λ</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mo>=</mml:mo><mml:mrow><mml:mn>?1</mml:mn><mml:mo>/</mml:mo><mml:mn>9</mml:mn></mml:mrow></mml:mrow></mml:math>λ2=?1/9 (72), and computed the MNE resolution matrix <mml:math><mml:mrow><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow><mml:mo>=</mml:mo><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:msub><mml:mi>??</mml:mi><mml:mi mathvariant="script">B</mml:mi></mml:msub></mml:mrow></mml:mrow></mml:math>??MNE(??B)=??MNE(??B)??B (30). We used a similar process to compute the MNE resolution matrix <mml:math><mml:mrow><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>??MNE(??Br). The SP estimates for the source space <mml:math><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:math>Br, with prewhitened gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub></mml:math>??Br, were obtained using the procedure in SI Appendix, section SI III. We characterized the SP resolution matrix <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>SP</mml:mtext></mml:msup></mml:math>??SP empirically. We simulated unit currents in the most significant eigenmode of each brain division <mml:math><mml:mi>i</mml:mi></mml:math>i and used the gain matrix <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mi>i</mml:mi><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:msub></mml:math>??Br(i) to generate the noiseless (prewhitened) MEG fields <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:math>??i. We then used SP to estimate the source location <mml:math><mml:mrow><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:mtext>SP</mml:mtext><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi mathvariant="normal">y</mml:mi><mml:mi mathvariant="normal">i</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub></mml:msub><mml:mo>,</mml:mo><mml:mi>L</mml:mi><mml:mo>=</mml:mo><mml:mn>?1</mml:mn><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>qi=SP(yi,??Br,L=?1) and computed the source current <mml:math><mml:msub><mml:mi>β</mml:mi><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:mrow><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mi>i</mml:mi></mml:mrow><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:msub></mml:math>β{qi,i}, using a least-squares fit of <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:msub><mml:mi mathvariant="script">B</mml:mi><mml:mi>r</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:msub></mml:math>??Br(qi) to <mml:math><mml:msub><mml:mi>??</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:math>??i. The estimated <mml:math><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub></mml:math>qi and <mml:math><mml:msub><mml:mi>β</mml:mi><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:mrow><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mo>,</mml:mo><mml:mi>i</mml:mi></mml:mrow><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:msub></mml:math>β{qi,i} together specify the locations and magnitudes of the nonzero elements of the empirical resolution matrix <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>SP</mml:mtext></mml:msup></mml:math>??SP. As <mml:math><mml:mrow><mml:mi>L</mml:mi><mml:mo>=</mml:mo><mml:mn>1</mml:mn></mml:mrow></mml:math>L=1, only the <mml:math><mml:mrow><mml:msub><mml:mi>q</mml:mi><mml:mi>i</mml:mi></mml:msub><mml:mtext>th</mml:mtext></mml:mrow></mml:math>qith element of column <mml:math><mml:mi>i</mml:mi></mml:math>i of <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>SP</mml:mtext></mml:msup></mml:math>??SP is nonzero. In the resolution matrices, we ordered the sources according to their physical proximity to each other. Therefore, when sources are estimated accurately, the resolution matrix has a diagonal appearance. We used the resolution matrices <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>MNE</mml:mtext></mml:msup></mml:math>??MNE and <mml:math><mml:msup><mml:mi>??</mml:mi><mml:mtext>SP</mml:mtext></mml:msup></mml:math>??SP to compute the spatial dispersion <mml:math><mml:mrow><mml:msub><mml:mtext>SD</mml:mtext><mml:mi>i</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:msqrt><mml:mrow><mml:msubsup><mml:mo largeop="true" symmetric="true">∑</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:mo>=</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mi>N</mml:mi></mml:msubsup><mml:mrow><mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:msub><mml:mi>d</mml:mi><mml:mrow><mml:mi>k</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub><mml:mo>?</mml:mo><mml:mrow><mml:mo fence="true">||</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>k</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub><mml:mo fence="true">||</mml:mo></mml:mrow></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup><mml:mo>/</mml:mo><mml:mrow><mml:msubsup><mml:mo largeop="true" symmetric="true">∑</mml:mo><mml:mrow><mml:mi>k</mml:mi><mml:mo>=</mml:mo><mml:mn>1</mml:mn></mml:mrow><mml:mi>N</mml:mi></mml:msubsup><mml:msup><mml:mrow><mml:mo fence="true">||</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>k</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub><mml:mo fence="true">||</mml:mo></mml:mrow><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:mrow></mml:msqrt></mml:mrow></mml:math>SDi=∑k=1N(dki?||??ki||)2/∑k=1N||??ki||2 and the dipole localization error <mml:math><mml:mrow><mml:msub><mml:mtext>DLE</mml:mtext><mml:mi>i</mml:mi></mml:msub><mml:mo>=</mml:mo><mml:msub><mml:mi>d</mml:mi><mml:mrow><mml:mi>j</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub></mml:mrow></mml:math>DLEi=dji, where <mml:math><mml:mrow><mml:mi>j</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mtext>arg</mml:mtext><mml:mrow><mml:msub><mml:mi>max</mml:mi><mml:mi>k</mml:mi></mml:msub><mml:mrow><mml:mo stretchy="false">{</mml:mo><mml:mrow><mml:mo fence="true">||</mml:mo><mml:msub><mml:mi>??</mml:mi><mml:mrow><mml:mi>k</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub><mml:mo fence="true">||</mml:mo></mml:mrow><mml:mo stretchy="false">}</mml:mo></mml:mrow></mml:mrow></mml:mrow></mml:mrow></mml:math>j=argmaxk{||??ki||} and <mml:math><mml:msub><mml:mi>d</mml:mi><mml:mrow><mml:mi>j</mml:mi><mml:mi>i</mml:mi></mml:mrow></mml:msub></mml:math>dji is the distance between centroids of divisions <mml:math><mml:mi>j</mml:mi></mml:math>j and <mml:math><mml:mi>i</mml:mi></mml:math>i (31).

Source Estimation Algorithm.

At each hierarchy level, we performed source estimation using the SP algorithm in ref. 36, modified to use relaxed mutual coherence thresholds. The mutual coherence thresholds were set equal to the mean-max correlation among gain matrices from neighborhoods of brain divisions under consideration. Our relaxed mutual coherence thresholds adapt to changing levels of gain matrix correlation across hierarchy levels, giving correlation thresholds in the <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mrow><mml:mn>0.75</mml:mn><mml:mo>?</mml:mo><mml:mn>0.90</mml:mn></mml:mrow></mml:mrow></mml:math>~0.75?0.90 range. Additional details are in SI Appendix, section SI III, Table S2, and Fig. S4.

Somatosensory Evoked Response Simulations.

We simulated activity in five regions of interest: a <mml:math><mml:mrow><mml:mn>1</mml:mn><mml:msup><mml:mtext>-cm</mml:mtext><mml:mn>3</mml:mn></mml:msup></mml:mrow></mml:math>1-cm3 volume in the left somatosensory thalamus including the ventral posterior area and four <mml:math><mml:mrow><mml:mo>~</mml:mo><mml:mn>600</mml:mn></mml:mrow></mml:math>~600- to <mml:math><mml:mrow><mml:mn>800</mml:mn><mml:msup><mml:mtext>-mm</mml:mtext><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:math>800-mm2 surface patches in primary and secondary somatosensory cortices and the posterior parietal area. The simulated cortical current time courses were Gabor atoms of the form <mml:math><mml:mrow><mml:mi>A</mml:mi><mml:msup><mml:mi>e</mml:mi><mml:mrow><mml:mo>?</mml:mo><mml:mrow><mml:mrow><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mi>t</mml:mi><mml:mo>?</mml:mo><mml:msub><mml:mi>t</mml:mi><mml:mi>o</mml:mi></mml:msub></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow><mml:mo>/</mml:mo><mml:mn>2</mml:mn></mml:mrow><mml:msup><mml:mi>σ</mml:mi><mml:mn>2</mml:mn></mml:msup></mml:mrow></mml:mrow></mml:msup></mml:mrow></mml:math>Ae?(t?to)/2σ2, where <mml:math><mml:mi>A</mml:mi></mml:math>A, <mml:math><mml:msub><mml:mi>t</mml:mi><mml:mi>o</mml:mi></mml:msub></mml:math>to, and <mml:math><mml:mi>σ</mml:mi></mml:math>σ denote the amplitude, delay, and width of the evoked component and were set based on previous studies (15, 36, 58). The simulated thalamic current time course consisted of <mml:math><mml:mn>10</mml:mn></mml:math>10 repetitions of <mml:math><mml:mrow><mml:msup><mml:mi>cos</mml:mi><mml:mn>2</mml:mn></mml:msup><mml:mrow><mml:mo stretchy="false">(</mml:mo><mml:mrow><mml:mrow><mml:mn>2</mml:mn><mml:mo>?</mml:mo><mml:mi>π</mml:mi><mml:mo>?</mml:mo><mml:mi>f</mml:mi><mml:mo>?</mml:mo><mml:mi>t</mml:mi></mml:mrow><mml:mo>+</mml:mo><mml:mi>?</mml:mi></mml:mrow><mml:mo stretchy="false">)</mml:mo></mml:mrow></mml:mrow></mml:math>cos2(2?π?f?t+?) with <mml:math><mml:mrow><mml:mi>f</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mn>?100</mml:mn><mml:mtext>Hz</mml:mtext></mml:mrow></mml:mrow></mml:math>f=?100Hz, <mml:math><mml:mrow><mml:mi>?</mml:mi><mml:mo>=</mml:mo><mml:mrow><mml:mi>π</mml:mi><mml:mo>/</mml:mo><mml:mn>3</mml:mn></mml:mrow></mml:mrow></mml:math>?=π/3, duration <mml:math><mml:mrow><mml:mn>15</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>15ms, and repetition period <mml:math><mml:mrow><mml:mn>25</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>25ms. We simulated currents over a time range of <mml:math><mml:mrow><mml:mn>0</mml:mn><mml:mo>–</mml:mo><mml:mn>250</mml:mn><mml:mo>?</mml:mo><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>0–250?ms, with <mml:math><mml:mrow><mml:mn>3</mml:mn><mml:mtext>-kHz</mml:mtext></mml:mrow></mml:math>3-kHz sampling rate. We calculated the MEG signals using the MRI-based forward model and added white Gaussian observation noise with mean zero and variance specified to achieve a SNR of <mml:math><mml:mrow><mml:mn>7</mml:mn><mml:mtext>dB</mml:mtext></mml:mrow></mml:math>7dB.

Auditory Evoked Response Recordings.

We recorded M/EEG simultaneously while presenting binaural broadband click trains (<mml:math><mml:mrow><mml:mn>0.1</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>0.1ms duration, <mml:math><mml:mrow><mml:mn>65</mml:mn><mml:mo>–</mml:mo><mml:mrow><mml:mrow><mml:mn>80</mml:mn><mml:mo>?</mml:mo><mml:mtext>dB</mml:mtext></mml:mrow><mml:mo>/</mml:mo><mml:mtext>nHL</mml:mtext></mml:mrow></mml:mrow></mml:math>65–80?dB/nHL intensity, <mml:math><mml:mrow><mml:mn>110</mml:mn><mml:mtext>ms</mml:mtext></mml:mrow></mml:math>110ms interstimulus interval, <mml:math><mml:mrow><mml:mn>9.09</mml:mn><mml:mtext>Hz</mml:mtext></mml:mrow></mml:math>9.09Hz click rate) as subjects rested with eyes open (38, 73). We preprocessed the data, performed stimulus-locked averaging to compute the ABR and the MLR, estimated observation noise covariances using the filtered baseline recordings (38) (SI Appendix, Fig. S8), whitened the filtered MEG and EEG data using their respective noise covariances, and used both MEG and EEG for source estimation. Further details are in SI Appendix, section SI IV.

Materials Availability.

All code and data from this paper are shared online on our MNE website: martinos.org/mne.

Acknowledgments

We acknowledge data collection assistance from Samantha Huang, Stephanie Rossi, and Tommi Raij; helpful discussions on source space construction with Koen Van Leemput, Imam Aganj, Doug Greve, and Bruce Fischl; and helpful discussions on canonical correlations and sparsity with Demba Ba and Emery N. Brown. This work was supported by NIH Grants P41-EB015896, 5R01EB022889, and 5R01-EB009048, NIH Grant 1S10RR031599-01 (to M.S.H.), and NIH Grant DP2-OD006454 (to P.L.P.).

Footnotes

  • ?1M.S.H. and P.L.P. contributed equally to this work.

  • ?2To whom correspondence may be addressed. Email: patrickp{at}nmr.mgh.harvard.edu or msh{at}nmr.mgh.harvard.edu.
  • Author contributions: P.K., B.B., M.S.H., and P.L.P. designed research; P.K., G.O.-H., J.A., S.K., and J.E.I. performed research; P.K., G.O.-H., J.A., S.K., B.B., M.S.H., and P.L.P. analyzed data; and P.K., M.S.H, and P.L.P. wrote the paper.

  • The authors declare no conflict of interest.

  • This article is a PNAS Direct Submission.

  • This article contains supporting information online at www.danielhellerman.com/lookup/suppl/doi:10.1073/pnas.1705414114/-/DCSupplemental.

This is an open access article distributed under the PNAS license.

References

  1. ?
    .
  2. ?
    .
  3. ?
    .
  4. ?
    .
  5. ?
    .
  6. ?
    .
  7. ?
    .
  8. ?
    .
  9. ?
    .
  10. ?
    .
  11. ?
    .
  12. ?
    .
  13. ?
    .
  14. ?
    .
  15. ?
    .
  16. ?
    .
  17. ?
    .
  18. ?
    .
  19. ?
    .
  20. ?
    .
  21. ?
    .
  22. ?
    .
  23. ?
    .
  24. ?
    .
  25. ?
    .
  26. ?
    .
  27. ?
    .
  28. ?
    .
  29. ?
    .
  30. ?
    .
  31. ?
    .
  32. ?
    .
  33. ?
    .
  34. ?
    .
  35. ?
    .
  36. ?
    .
  37. ?
    .
  38. ?
    .
  39. ?
    .
  40. ?
    .
  41. ?
    .
  42. ?
    .
  43. ?
    .
  44. ?
    .
  45. ?
    .
  46. ?
    .
  47. ?
    .
  48. ?
    .
  49. ?
    .
  50. ?
    .
  51. ?
    .
  52. ?
    .
  53. ?
    .
  54. ?
    .
  55. ?
    .
  56. ?
    .
  57. ?
    .
  58. ?
    .
  59. ?
    .
  60. ?
    .
  61. ?
    .
  62. ?
    .
  63. ?
    .
  64. ?
    .
  65. ?
    .
  66. ?
    .
  67. ?
    .
  68. ?
    .
  69. ?
    .
  70. ?
    .
  71. ?
    .
  72. ?
    .
  73. ?
    .

Online Impact

    <acronym id="UPyyYwe"></acronym>
    <rt id="UPyyYwe"></rt>
    <acronym id="UPyyYwe"></acronym>
    <acronym id="UPyyYwe"><optgroup id="UPyyYwe"></optgroup></acronym><acronym id="UPyyYwe"><small id="UPyyYwe"></small></acronym>
    <tr id="UPyyYwe"><optgroup id="UPyyYwe"></optgroup></tr>
    <tr id="UPyyYwe"><optgroup id="UPyyYwe"></optgroup></tr>
    <acronym id="UPyyYwe"></acronym>
  • 8189251275 2018-02-18
  • 6298941274 2018-02-18
  • 8345181273 2018-02-18
  • 207841272 2018-02-18
  • 2683681271 2018-02-18
  • 5067491270 2018-02-18
  • 2051721269 2018-02-18
  • 2999231268 2018-02-18
  • 183621267 2018-02-18
  • 5236401266 2018-02-18
  • 2592991265 2018-02-18
  • 9896941264 2018-02-18
  • 1171081263 2018-02-18
  • 983551262 2018-02-18
  • 3896031261 2018-02-18
  • 4643431260 2018-02-18
  • 4122621259 2018-02-18
  • 336531258 2018-02-17
  • 6455421257 2018-02-17
  • 5128821256 2018-02-17