Deep brain source imaging with M/EEG and anatomical MRI

Abstract

A method for non-invasively resolving electrophysiological activity in sub-cortical structures located deep in the brain by comparing amplitude-insensitive M/EEG field patterns arising from activity in subcortical and cortical sources under physiologically relevant sparse constraints is disclosed. The method includes a sparse inverse solution for M/EEG subcortical source modeling. Specifically, the method employs a subspace-pursuit algorithm rooted in compressive sampling theory, performs a hierarchical search for sparse subcortical and cortical sources underlying the measurement, and estimates millisecond-scale currents in these sources to explain the data. The method can be used to recover thalamic and brainstem contributions to non-invasive M/EEG data, and to enable non-invasive study of fast timescale dynamical and network phenomena involving widespread regions across the human brain.

Claims

1. A method for non-invasively characterizing electrophysiological activity within deep brain regions based on data measurements acquired during a functional task of interest, the method comprising: acquiring measured data representative of brain activity of a subject during a functional task of interest, acquiring magnetic resonance images of brain anatomy of the subject, constructing superficial (cortical) source spaces and deep (subcortical) source spaces for regions of the brain from the magnetic resonance images, the cortical source spaces including a first set of cortical source spaces having patches with a first size, a second set of cortical source spaces having patches with a second size, the second size being smaller than the first size; determining a first relationship between the source spaces and the measured data; determining a minimum norm source current estimate (MNE) based on the source spaces, the first relationship, and the measured data, the MNE including a first subset of source spaces; refining the MNE a first time by using and determining a second relationship between the measured data, and cortical source spaces from the second set of cortical source spaces that at least partially overlap with the first subset of source spaces, the refined MNE including a second subset of source spaces, the second subset of source spaces including a cortical source space; further refining the MNE by using and determining a third relationship between the second subset of source spaces, the subcortical source spaces, and the measured data, the further refined MNE including a superficial source space and a deep source space that underlies the measured data.

2. The method of claim 1, wherein the measured data is obtained by at least one of an electroencephalography (EEG) and an magnetoencephalography (MEG).

3. The method of claim 1, wherein constructing source spaces from the magnetic resonance images of brain anatomy of the subject includes positioning dipoles within a given source space, and wherein at least one source space corresponds to at least one of: a neocortex, a hippocampus, a gray-white matter interface, and a subcortical volume.

4. The method of claim 1, wherein each cortical source space is a surface patch having dipoles, and wherein the subcortical source spaces include volume subdivisions and surface patches having dipoles.

5. The method of claim 4, wherein the surface patches and volume subdivisions are sized to homogenize current strengths.

6. The method of claim 4, wherein the surface patches and the volume subdivisions having high current densities are assigned finer subdivisions, and wherein the surface patches and the volume subdivisions having low current densities are assigned larger subdivisions to homogenize current strengths thereby enabling comparable consideration of the surface patches and volume subdivisions in an inverse solution.

7. The method of claim 1, wherein determining the second relationship includes searching for cortical source spaces within the second set of cortical source spaces that at least partially overlap with the first subset of source spaces that correlate to the measured data with least residual error to determine the second subset of source spaces.

8. The method of claim 1, wherein constructing cortical source spaces includes constructing a third set of cortical source spaces having patches with a third size, the third size being smaller than the first size and the second size, and further comprising refining the MNE a second time by using and determining a relationship between the measured data, and cortical source spaces from the third set of cortical source spaces that at least partially overlap with the second subset of source spaces, the refined MNE from the second time including a third subset of source spaces, the third subset of source spaces including a cortical source space, and wherein further refining the MNE includes using and determining a fourth relationship between the third subset of source spaces, the subcortical source spaces, and the measured data, to identify the superficial source space and the deep source space that underlies the measured data.

9. The method of claim 8, wherein the cortical source spaces from the second set of cortical source spaces that at least partially overlap with the first subset of source spaces include cortical source spaces from the second set of cortical source spaces that are closely located to the first subset of source spaces.

10. A method for imaging electromagnetic sources within deep brain structures based on M/EEG data identified while a subject is performing a functional task of interest, the method comprising: constructing source spaces in superficial (cortical) and deep (subcortical) regions of the brain from images of the brain, the source spaces including a first set of patches of the superficial and deep regions of a first preselected level of coarseness, and a second set of patches at a second preselected level of coarseness, each of the patches in each set being sized to homogenize current strengths, and wherein the first set of patches and the second set of patches span across the same surface; relating the source spaces with the M/EEG data; computing a minimum norm source current estimate (MNE) based on the relationship between the source spaces and the M/EEG data, the source spaces, and the M/EEG data, refining the MNE by iteratively reducing and expanding a subset of source spaces that are defined within the MNE to identify at least one superficial and at least one deep source whose activity underlies the M/EEG data.

11. The method of claim 10, wherein the homogenization of the patches enables weaker and stronger gain regions to have comparable consideration in an inverse solution.

12. The method of claim 10, wherein refining the MNE includes determining a correlation between a neighboring source space and a given source space within the subset of source spaces.

13. The method of claim 12, further comprising comparing the correlation to a mutual coherence threshold, wherein the neighboring source space is included in the MNE, based on the correlation being less than the mutual coherence threshold, and wherein the neighboring source space is excluded from the MNE, based on the correlation being greater than the mutual coherence threshold.

14. A method for electromagnetic source imaging based on non-invasive M/EEG recordings and MM-based anatomic measures to employ information within M/EEG field patterns for estimating source currents across superficial and deep brain regions, the method comprising: constructing source spaces, the source spaces including a first set of cortical source spaces having patches with a first size, a second set of cortical source spaces having patches with a second size, and a third set of subcortical source spaces having a third size, wherein the second size being smaller than the first size, and wherein the first set of source spaces and the second set of source spaces span across the same surface; computing a minimum norm source current estimate (MNE) using the source spaces and the non-invasive M/EEG recordings; iteratively refining the MNE until the MNE includes a subset of source spaces having at least one superficial and at least one deep source whose activity underlies the non-invasive M/EEG recordings.

15. The method of claim 14, wherein the iterative refinement of the MNE includes: comparing the first set of cortical source spaces and the M/EEG recordings to determine a first subset of cortical source spaces; determining a second set of cortical source spaces, the second set of cortical source spaces being from the second set of cortical source spaces, the second set of cortical source spaces at least partially overlapping with the first subset of cortical source spaces; and determining a second subset of cortical sources spaces within the second set of cortical source spaces using the M/EEG recordings, the second subset of cortical sources determining the subset of source spaces.

16. The method of claim 14, further comprising acquiring magnetic resonance images of the brain and using the images to construct the source spaces in for the cortical and subcortical regions of the brain.

17. The method of claim 14, wherein the third set of subcortical source spaces include a plurality of subcortical volume subdivisions.

18. The method of claim 17, wherein the patches of the first set of cortical source spaces and the plurality of subcortical volume subdivisions of the third set of subcortical source spaces are sized to have homogenous current strength.

19. A subspace pursuit process for characterizing electrophysiological activity within deep brain regions based on non-invasive M/EEG measurements, the pursuit process being performed on a plurality of source spaces, the source spaces including a plurality of cortical surface patches and subcortical volume subdivisions for regions of the brain, the process comprising: acquiring measured data representative of brain activity of a subject, the measured data being indicative of electrophysiology activity within deep brain regions of the subject; computing a minimum norm source current estimate (MNE) using the cortical surface patches and the measured data; refining the MNE by performing a subspace pursuit to identify a first cortical surface patch that correlates to the measured data, the first cortical surface patch having a first size; further refining the MNE by performing a second subspace pursuit on a subset of cortical surface patches each having second size to determine a second cortical surface patch that correlates to the measured data, the second cortical surface patch being smaller than the first cortical surface patch and the subset of cortical surface patches, the second size being smaller than the first size, the subset of cortical surface patches at least partially overlapping with the first cortical surface patch; and performing a third subspace pursuit between the second cortical surface patch and the subcortical volume subdivisions to determine cortical and subcortical brain regions whose activity correlate best to the measured data.

20. The process of claim 19, wherein the patches are sized, such that a given patch has a first current strength and a second patch has a second current strength, and wherein the first current strength and the second current strength are homogenous.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a functional diagram illustrating electromagnetic source imaging with M/EEG and anatomical MRI.

(2) FIG. 2 is a flow chart illustrating the general series of steps for deep brain electromagnetic source imaging using M/EEG recordings and anatomical MRI.

(3) FIG. 3A is a flow chart illustrating the steps for defining cortical and subcortical source spaces.

(4) FIG. 3B illustrates a series of cortical surfaces divided into patches of varying sizes. Each patch is a different color. The leftmost frame of FIG. 3B includes a cortical surface divided into patches of a first size, the central frame of FIG. 3B includes a cortical surface divided into patches of a second size, and the leftmost frame of FIG. 3B includes a cortical surface divided into patches of a third size.

(5) FIG. 3C illustrates subcortical volumes comprising subdivided anatomic segmentations. Each subdivision is a different color. The leftmost frame of FIG. 3C is a coronal plane view, the central frame of FIG. 3C is a horizontal plane view, and the rightmost frame of FIG. 3C is a sagittal plane view.

(6) FIGS. 4A-D illustrate that sparsity constraints enable distinctions between deep subcortical and superficial cortical contributions to M/EEG field patterns. Each panel contains histograms of principal angles between MEG field patterns generated by activity in pairs of brain regions. The principal angles quantify distinctions between subspaces defining field patterns from the two regions (0 degree indicates identical, and 90 degree indicates fully orthogonal). (A vs. B) Full cortical source space can generate most any field pattern that a subcortical source space can, while sparse cortical sources cannot. (B vs. C) Sparse subcortical and cortical source spaces generate field patterns that are just as distinct as those generated by pairs of sparse cortical source spaces. (C vs. D) Pairs of subcortical regions also generate distinct field patterns, like pairs of cortical regions, but exhibit a bit more overlap due to their greater distance from sensors. (B-D) A large proportion of the principal angles are high and close to orthogonal when sparse regions are compared.

(7) FIGS. 5A-F illustrate an approach to exploit sparsity constraints for the deep brain source imaging problem. In particular, the diagram shows a hierarchical search to identify relevant sparse (and distinguishable) subcortical and cortical sources underlying the measurement. The hierarchical reduction of the full brain source space to relevant sparse distinguishable subsets is achieved as below: (A) the space of all possible sources across the brain. (B) A sparse cortical region most relevant to the measurement under study. (C) Refined sparse cortical regions within and around the region in Panel B. Different colors are used for neighboring patches to highlight patch boundaries. (D) Refined sparse cortical region most relevant to the measurement. (E) Joint source space comprising the sparse cortical region and subcortical volumes. The sparse subset of the cortical source space allows subcortical areas to show through in the inverse solution. It is possible to search across these sparse and therefore distinguishable sets of sources to localize appropriate subcortical and cortical sources. (F) Final sparse cortical and subcortical regions most relevant to the measurement under study.

(8) FIG. 6 illustrates a hierarchical subspace pursuit algorithm for deep brain source imagingin accordance with the approach disclosed abovein simplified flow chart form.

(9) FIGS. 7A-B illustrate auditory evoked response recordings obtained to validate the efficacy of the deep brain imaging approach presented above. Each panel presents stimulus-locked average auditory evoked responses obtained on a healthy volunteer presented with a broadband click train stimulus. Data shown are averages across 11200 epochs. (A) Auditory brainstem responsesthe red bar marks periods when early peaks are seen in the data, the timing of these peaks is consistent with expected location of wave V. Time courses are raw recordings filtered 500-1625 Hz, rectified and averaged across channels. (B) Middle latency responses. Common peaks are marked and latencies follow known trends. Time courses are raw recordings filtered 0-300 Hz, displayed across multiple channels.

(10) FIGS. 8A-F illustrate cortical and subcortical estimates for the auditory brainstem recordings of FIG. 7. Estimates are obtained with the proposed sparse hierarchical deep brain imaging algorithm. (A) Reduced cortical source space obtained as a result of initial source space reduction performed upon the middle latency recordings. Estimates shown on the inflated surface are snapshots at 25 msec and thus comprise patches in primary auditory areas. Subspace pursuit (SP) and minimum norm estimation (MNE) algorithms were then applied on the reduced source space, to localize sources of the auditory brainstem recordings. The leftmost frame of FIG. 8A is a first view of the reduced cortical source space, the central frame of FIG. 8A is a second view of the reduced cortical source space, and the rightmost frame of FIG. 8A is a third view of the reduced cortical source space. (B-C) Representative MRI maps showing snapshots of subspace pursuit and MNE estimates at 5 msec (around appearance of wave V). MRI slice locations in the surface right-anterior-superior (RAS) coordinate system are 115 mm (coronal for the leftmost frame of FIGS. 8B and 8C), 126 mm (horizontal for the central frame of FIGS. 8B and 8C) and 142 mm (sagittal for the leftmost frame of FIGS. 8B and 8C). (D) Summary of mean activity by anatomic region for each estimation algorithm. Legend: cortical (co), putamen (r/lp), caudate (r/lc), thalamus (r/lt), hippocampus (r/lh), lateral geniculate (r/ll), medial geniculate (r/lm), and inferior colliculus (ic), where r/l refer to right/left. (E-F) Time courses of source current estimates obtained with SP and MNE algorithms. (D-F) SP estimates indicate wave V arises from focal sources in the inferior colliculusas in established literature. SP is more specific to the inferior colliculus than MNE. All estimates displayed are resultants across dipoles within their respective patch or volume subdivisions. Each time course shown corresponds to a selected subdivision in the respective anatomic regions indicated in legends.

DETAILED DESCRIPTION OF THE DISCLOSURE

(11) The present disclosure describes methods of electromagnetic source imaging with non-invasive M/EEG recordings and MRI-based anatomic measures that can be used to acquire gain-insensitive information in M/EEG arising from subcortical and cortical structures, and provides a hierarchical subspace pursuit algorithm to estimate neural currents in the subcortical structures. Distinctions between field patterns can be used to localize subcortical sources and distinguish subcortical versus cortical contributions using the hierarchical subspace pursuit process or algorithm.

(12) Referring now to FIG. 2, a flowchart illustrates a series of steps for estimating source currents from deep brain regions with a hierarchical subspace pursuit algorithm. Initially, M/EEG data are acquired during a neurophysiologic paradigm of interest. These recordings are complemented with magnetic resonance images of the subject's brain anatomy (step 10). Next, distributed source spaces are constructed from the images (step 12) and electromagnetic forward solutions G are computed using Maxwell's equations (step 14). Analysis of these forward solutions reveals that a sparse inverse solution would enable distinctions between cortical and subcortical sources of electromagnetic fields recorded with M/EEG. Therefore, we disclose a sparse hierarchical inverse solution for deep brain electromagnetic source imaging. In sequence, this involves computation of the minimum norm estimates (MNE, step 16), refinement of these estimates with a subspace pursuit scheme (step 18) and iterative evaluation of the MNE+subspace pursuit procedure across hierarchies of cortical and subcortical source spaces. This iterative process enables identification of both superficial and deep locations of brain activity, and allows estimation of the source currents therein (step 20). Each of these steps is described more fully below.

Acquire Images

(13) Referring still to FIG. 2, and now also to FIG. 3A, after the magnetic resonance images are acquired, distributed neural source spaces are constructed (step 12). Image analysis software can be used to reconstruct and triangulate neocortical and hippocampal surfaces (step 22), as well as to segment the subcortical volumes by anatomic region (step 24). Subsequently, dipoles are positioned at vertices in the triangulated surface mesh for neocortex and hippocampus, at the gray-white matter interface, and assigned orientations fixed normal to the neocortical and hippocampal surfaces respectively (step 26). The neocortical and hippocampal dipole placements are performed consistently as both structures have a similar layered cyto-architecture, and as the hippocampal inner structure is not accessible with anatomical MRI. Dipoles are also placed in subcortical volumes (thalamus, caudate, putamen, amygdala, medial and lateral geniculate nuclei, brainstem) extracted from the segmentations, as well as in specific smaller nuclei (medial and lateral geniculate) defined based on standard anatomic landmarks and segmentations. (step 28) Although other configurations are possible, the dipoles can be positioned at a 1 mm voxel spacing, and oriented randomly. Although a number of different segmentation software packages could be used, the FreeSurfer Software Suite available at freesurfer.net has been used in exemplary applications. Other suitable software packages could include the brainVISA/Anatomist package (brainvisa.info) and the FMRIB Software Library (fsl.fmrib.ox.ac.uk/fsl/fslwiki).

Construct Source Spaces

(14) Referring still to FIG. 2, and now also to FIG. 3A, after the magnetic resonance images are acquired, distributed neural source spaces are constructed (step 12). Image analysis software can be used to reconstruct and triangulate neocortical and hippocampal surfaces (step 22), as well as to segment the subcortical volumes by anatomic region (step 24). Subsequently, dipoles are positioned at vertices in the triangulated surface mesh for neocortex and hippocampus, at the gray-white matter interface, and assigned orientations fixed normal to the neocortical and hippocampal surfaces respectively (step 26). The neocortical and hippocampal dipole placements are performed consistently as both structures have a similar layered cyto-architecture, and as the hippocampal inner structure is not accessible with anatomical MRI. Dipoles are also placed in subcortical volumes (thalamus, caudate, putamen, amygdala, medial and lateral geniculate nuclei, brainstem) extracted from the segmentations, as well as in specific smaller nuclei (medial and lateral geniculate) defined based on standard anatomic landmarks and segmentations. (step 28) Although other configurations are possible, the dipoles can be positioned at a 1 mm voxel spacing, and oriented randomly. Although a number of different segmentation software packages could be used, the FreeSurfer Software Suite available at http://freesurfer.net/ has been used in exemplary applications. Other suitable software packages could include the brainVISA/Anatomist package (http://brainvisa.info/) and the FMRIB Software Library (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/).

(15) To reduce the dimensionality of the inverse problem, and make it better posed for the number of sensors typically available in M/EEG recordings, the dipoles are grouped into surface patches for neocortex and hippocampus, and volume subdivisions for subcortical volumes (steps 30, 32, 34, 36). This process does not cause loss of information as small groups of nearby dipole clusters generate similar or correlated electromagnetic fields at the external sensors. The process for grouping dipoles into patches and subdivisions is described below.

(16) First, cortical and hippocampal dipoles are grouped into patches. This can be accomplished within the MNE software package available at martinos.org/mne/stable/index. In one example, the cortical and hippocampal surfaces were each approximated with the topology of recursively subdivided icosahedra, i.e., with a solid figure with twenty plane faces. For each icosahedron subdivision, detailed surface geometry information were used to derive surface patches.

(17) Second, subcortical dipoles are grouped into volume subdivisions obtained by subdividing the anatomic subcortical segmentations derived from the FreeSurfer Software Suite. In particular, the volume subdivisions are sized to homogenize current strengths within dipole groups across the brain. In one example, the volume subdivisions were sized to homogenize current strengths across cortical and subcortical regions as follows. The surface current densities , and volume current densities , were specified across brain regions in a structure-specific manner (steps 30, 32). Using neocortical surface current density and average cortical patch area , the average current strength for cortical patches was computed as . Next, total volumes V for anatomically distinct subcortical regions were computed. Then, the number of subdivisions for a given subcortical volume segmentation was derived as =V/. Finally, using a k-means based clustering algorithm, the dipoles within the anatomical region are grouped into n distinct equally sized clusters. The typical volume of a subcortical subdivision can range 150-800 mm.sup.3, depending on structure.

(18) The process of homogenizing current strengths enables weaker and stronger gain regions to have comparable consideration in an inverse solution, while implicitly assigning higher resolution (finer subdivisions) to regions with high current densities and lower resolution (larger subdivisions) to regions with low current densities. While the above-mentioned methods detail one way for sizing the surface patches and volume subdivisions to homogenize current strengths across sources, several other decomposition methods can be used.

(19) Referring again to FIG. 3A (step 38), and now also to FIG. 3B-3C, example cortical and subcortical source spaces resulting from the above source space construction techniques are illustrated. In FIG. 3B, cortical patches with three different sizes (varying coarseness) for use in a hierarchical inverse solution are shown. The set of coarse patches has an average area 2500 mm.sup.2 (on left) and is denoted .sup.(1), the set of finer patches with average area 650 mm.sup.3 (middle) is denoted .sup.(2), and the set of even finer patches with average area 175 mm.sup.2 (right) is denoted .sup.(3). In FIG. 3C, a subcortical source space S comprising volume subdivisions across caudate, putamen, hippocampus, amygdala, thalamus, lateral and medial geniculate nuclei, and brainstem is shown. These subdivisions have volume ranging 150-800 mm.sup.3 and are illustrated in different planes of view (coronal, horizontal and sagittal).

Compute Forward Solutions

(20) Referring again to FIGS. 1 and 2, in step 14, it is possible to calculate a series of electromagnetic forward solutions for each source in the above source space. For a given source (surface patch or volume subdivision) of known location, Maxwell's equations determine the mapping between source activity and electromagnetic fields at the sensors i.e., the electromagnetic forward solution G. Then, activity in dipoles located within a given source region, denoted by source currents X, generates a vector y of field patterns at the M/EEG sensors as below:
y.sub.N1=G.sub.NMX.sub.M1
where N is the number of sensors and M is the number of dipole sources in the region of interest. The mapping specified in G is linear as it arises from a quasi-static approximation to Maxwell's equations.

(21) In an exemplary M/EEG application, a Neuromag Vectorview MEG system (Elekta Neuromag, Helsinki, Finland) comprising 306 sensors across 102 sensor locations (2 planar gradiometers and 1 magnetometer per location) and a standard 70-lead EEG electrode configuration were used to calculate the forward solutions. To ensure alignment between the brain source coordinates and sensor coordinates used for forward solutions, the positions of M/EEG sensors and corresponding fiduciary points were digitally co-registered with the MRI data acquired for source space construction (step 10). For each regional subdivision in the source space, MEG and EEG forward solutions were generated within the MNE software package, using a three-compartment boundary element model derived from the MRI data. Dipoles within 5 mm of the inner skull bounding surface were excluded to avoid numerical errors. To account for non-homogeneity between the current densities across regions in the brain, the G for each region was multiplied by regional current strength denoted as c. For each subcortical region, current strength was derived as c=V where and V are the regional current density and volume respectively. For each neocortical and hippocampal patch, current strength was derived as c=A where and A are the surface current density and area respectively.

(22) For any subdivision in the source space, the columns of G obtained above represent fields arising from small collections of neighboring dipoles, and thus are highly correlated. Thus, information in the forward solutions for individual regions can be condensed into a low-dimensional approximation i.e., into a low-rank basis derived using a reduced singular value decomposition
G=U.sub.NBS.sub.BBW.sub.BM
=U.sub.NPS.sub.PP
where P<B is the number of eigen-modes retained in the low-rank approximation. For each subdivision, P is chosen such that the reduced order represents at least 95% of the total spectral energy in G:

(23) $\frac{\underset{b=1}{\overset{P}{.Math.}}{S}_{\mathrm{bb}}}{\underset{b=1}{\overset{B}{.Math.}}{S}_{\mathrm{bb}}}0.95.$
Typically P=26 modes are sufficient to capture 95% of the spectral energy in G. This approach (a) enables reduction of the dimensionality of G, (b) minimizes representation error on average, and (c) implicitly defines reduced subspaces for regional forward solutions for use in inverse solutions. All following analyses use these low-rank approximations of the forward models.

(24) Analyze Forward Solutions to Motivate Inverse Approach

(25) The eigenmodes in represent normalized field patterns that can arise from some distribution of source activity in a region of interest. To understand how to distinguish subcortical and cortical contributions to M/EEG data within an inverse solution, it is useful to compare field patterns arising from subcortical vs. cortical activity. These comparisons can be quantified using metrics such as principal angles which are established metrics for subspace comparisons.

(26) Principal angles between subspaces of possible field patterns from regions i and j, defined by the reduced order electromagnetic forward models .sub.i and .sub.j respectively, are denoted .sub.i,j. The angles can be computed using a formulation that uses the singular value decomposition:

(27) ${\hat{G}}_i^{}{\hat{G}}_j={\mathrm{Up}}_{i,j}\underset{i,j}{.Math.}{\mathrm{Wp}}_{i,j},\underset{i,j}{.Math.}=\mathrm{diag}\left(s_1,s_2,\mathrm{.Math.}{s}_n\right)=\mathrm{diag}\left(\cos {}_1,\cos {}_2,\mathrm{.Math.}\cos {}_n\right)$${}_{i,j}=\left[{}_1,{}_2\mathrm{.Math.}{}_n\right].$

(28) FIG. 4A shows the principal angles between field patterns arising from activity in the subcortical regional subdivisions and the full cortical source space respectively. The low angles illustrate that the full cortical source space can explain most field patterns arising from a subcortical subdivision. This is because the full cortical case has several eigen-modes spanning most possible field patterns including those generated by activity within a subcortical volume subdivision. FIG. 4B illustrates angles between field patterns from sparse subcortical and cortical regions, which are significantly higher than those in FIG. 4A, illustrating near orthogonality between field patterns generated by sparse pairs of regions. FIG. 4C illustrates that pairs of cortical regions also give rise to near orthogonal field patterns. The angle trends in FIG. 4B-4C are similar, illustrating that subcortical and cortical sources can be distinguished using differences between their respective field patterns. FIG. 4D illustrates field patterns arising from pairs of subcortical regions, and shows similar trends to FIG. 4C. Therefore, in every case where sparse regions are compared, the field patterns are more distinguishable than the full cortical space case of FIG. 4A.

A Sparse Inverse Solution for Deep Brain Electromagnetic Source Imaging

(29) The above results demonstrate the key idea that sparsity constraints within an inverse solution can enable robust distinctions between subcortical and cortical contributions to M/EEG data (field patterns). Furthermore, sparsity constraints are suitable in a variety of neurophysiologic scenarios as only a small proportion of the brain has salient activity in many cognitive experiments and clinical conditions. These observations motivate a sparse inverse solution for the deep brain electromagnetic source imaging problem.

(30) One such sparse inverse solution is outlined in FIG. 2, steps 16, 18 and 20. The approach is to hierarchically reduce the space of all possible sources to sparse subsets most relevant to a given M/EEG recording. The electromagnetic inverse problem and steps in the hierarchical iterative sparse solution are described formally below.

(31) The electromagnetic inverse problem is: to estimate the source currents X underlying M/EEG measurements Y, given the electromagnetic forward solution G linking the two:
Y.sub.NT=G.sub.NMX.sub.MT+V.sub.NT
Here, the observation noise V is assumed to be Gaussian, temporally uncorrelated, and specified completely by spatial covariance Q.sub.NN independent of X, N is the number of sensors, M is the number of dipole sources across the brain, and T is the number of time points in the measurement.
For a distributed source space comprising K regions across the brain, columns of G and rows of X can be grouped by region and the measurement equation can be rewritten as:

(32) $Y_{NT}=\left[G_1,G_2,\mathrm{.Math.},G_K\right]\left[\begin{array}{c}{x}_1\\ x_2\\ \mathrm{.Math.}\\ x_k\end{array}\right]+V_{NT}$
where the set of regions (patches or subdivisions) is [1, 2, . . . , K], G.sub.k is the electromagnetic forward solution for the k.sup.th region, and x.sub.k denotes the source current of dipoles in the k.sup.th region across time. Thus, the norm of x.sub.k can be thought of as the quantum of activity from the k.sup.th region, and indicate the least squares match between the possible field patterns from a region and field patterns within the M/EEG data.

(33) Compute MNE Estimates:

(34) The classic solution for this problem is given by the minimum l.sub.2-norm estimate (MNE) which projects measurement Y onto the range of G while accounting for noise characteristics and prior information. The MNE solution for source currents X given full forward solution G, M/EEG measurements Y, noise statistics Q and prior covariance R on source currents is given by:
x.sup.MNE(Y,G,Q,R)=RG(GRG+Q).sup.1Y

(35) However, the above-mentioned principal angle results suggest that a simple application of MNE on the full brain source space, with G specifying the composite forward solution for all subcortical and cortical sources, will likely localize to the cortical source spaces even if the activity arises from a deeper subcortical source. Instead, as distinctions in subcortical and cortical field patterns arise when sparse source spaces are considered, localizing subcortical contributions requires refining the MNE estimates. An ideal inverse algorithm would have the ability to refine the classical MNE solution in two ways: (a) reduce the full distributed source space to accurate sparse cortical subsets alongside deep sources, and (b) concentrate the source current estimates into the appropriate subcortical sources.

(36) Refine MNE Estimates with Subspace Pursuit:

(37) Both of the above objectives are served well by projection pursuit methods, which identify sparse projections to best explain multivariate data. First, pursuit methods can find sparse cortical regions most relevant to the measured fields, amidst the high-dimensional dictionary of all regions and their corresponding field patterns. This is because they deal uniquely well with the curse of dimensionality. Second, upon identifying the sparse cortical contributors, the pursuit search can be repeated in the reduced source space to search for weak but distinct subcortical contributors. This is because pursuit methods are robust to noise and outliers, and can recover small amplitude features in large dictionaries. These methods can be implemented with subspace pursuit algorithms, which find a sparse projection matching a given dictionary (e.g. matching in the least squares or MNE sense), remove the component along that projection and iterate to find new projections, till all matching projections have been found.

(38) Formally, the subspace pursuit (SP) estimate is written as:
[H,{circumflex over (X)}.sub.H]=SP(Y,G,Q,R)
where H is a set of indices denoting the subset of regions deemed to be most relevant to the measurement. Implicitly H gives a subset of regions whose normalized field patterns match those in the measurement after compensating for noise, and {circumflex over (X)}.sub.H specifies source current estimates that can together generate these different field patterns to quantitatively explain measurement amplitudes.
Intuitively, this approach starts with the MNE estimate, iteratively picks out the relevant sparse components of the MNE estimate, and refines the MNE estimate until all relevant projections have been found.

(39) The above MNE+pursuit process is repeated and adapted across hierarchies of cortical and subcortical source spaces, so as to focus the estimate iteratively into the relevant regions (iteratively reduce the source space to the solution) for greater spatial resolution and accuracy. FIGS. 5A-F illustrates one instance of such a hierarchical source space reduction procedure. Starting from a nave space of all possible sources (FIG. 5A), it is possible to find a sparse cortical solution (FIG. 5B), perform a finer search in the neighborhood of this sparse cortical solution (FIG. 5C), obtain a refined sparse cortical solution (FIG. 5D). The sparse cortical solution is distinguishable from subcortical sources in general. Then, the search can be repeated with a combination of the reduced/sparse cortical space and the full subcortical space (FIG. 5E) to identify relevant subcortical and cortical sources (FIG. 5F). This iterative process enables identification of both superficial and deep locations of brain activity, and allows estimation of the source currents therein (FIG. 2, step 20).

(40) Formally, if the algorithm starts with a source space C.sup.(1), with forward model denoted as G.sup.(1),, MNE+subspace pursuit is performed to pick L regions that can best explain measured fields:

(41) $\left[H_1,{\hat{X}}_{H_1}\right]=\mathrm{SP}\left(Y,G^{\left(1\right)},Q,R,L\right)$
where H.sub.1 denotes the chosen set of L regions in C.sup.(1) (FIG. 5B) and {circumflex over (X)}.sub.H.sub.1 denotes the source current estimates in these chosen cortical regions. Next, we repeat the estimation in a finer cortical source space C.sup.(2)with patches overlapping H.sub.2 and their nearest neighbors (FIG. 5C). The forward model for these patches is denoted as G.sup.(2) and we perform subspace pursuit to pick L regions that can best explain measured fields:

(42) $\left[H_2,{\hat{X}}_{H_2}\right]=\mathrm{SP}\left(Y,G^{\left(2\right)},Q,R,L\right)$
where H.sub.2 denotes the chosen set of L regions in C.sup.(2) (FIG. 5D) and {circumflex over (X)}.sub.H.sub.2 denotes the source current estimates in these chosen regions. This systematic reduction of possible cortical sources can be repeated to the level of fineness relevant to cortical activity in the task under study. This hierarchical reduction process results in a highly spatially resolved subset of relevant sparse cortical sources, denoted as .sub.(sp) (e.g. as in FIG. 5D). Then, in the final stage search, the subcortical source space S, containing volume subdivisions of caudate, putamen, hippocampus, amygdala, thalamus, lateral and medial geniculate nuclei, brainstem, is considered alongside the sparse cortical source space .sub.(sp) (FIG. 5E). We denote the composite or joint source space as J=[.sub.(sp), S]. The forward model for J denoted as G.sub.J=[.sub.sp, G.sub.S], implicitly contains nearly decorrelated columns, and we perform subspace pursuit to pick L regions that can best explain measured fields:

(43) $\left[H_J,{\hat{X}}_{H_J}\right]=\mathrm{SP}\left(Y,G_J,Q,R,L\right)$
where H.sub.J specifies the set of chosen regions whose dipole activity best explains the data (e.g. as in FIG. 5F), while {circumflex over (X)}.sub.H.sub.J gives the source current estimates in these subdivisions. We denote these dipole currents as {circumflex over (x)} for notational convenience.

Hierarchical Subspace Pursuit-Based Algorithm for Deep Brain Electromagnetic Source Imaging

(44) Referring now to the flowchart in FIG. 6, the algorithm is illustrated in a stepwise manner. Steps 50, 52 and 54 describe the preprocessing of the M/EEG data and estimation of the noise covariance Q and source covariance R. Then in step 56, the target sparsity level L and the starting source space are initialized. Following this, an iterative process is entered. In step 58, the correlation threshold is calculated. Then, the MNE estimates are obtained (step 60) and reduced via subspace pursuit to the L sources that best explain the data while satisfying the correlation threshold (step 62). If there are additional source spaces to search within the hierarchy (step 64), the search space is updated (step 66) and steps 58-62 are repeated. The repetitions are performed until all of the hierarchical source spaces have been searched. The cortical and subcortical source locations most relevant to the measurement are then identified (step 68), and the neural source currents in these locations that explain the M/EEG data are estimated. We further detail the steps within the iterations below.

(45) The mutual coherence thresholds , described above, are used in the subspace pursuit algorithm to enforce incoherence during search. To estimate these thresholds, the maximum correlations between the modes of forward solutions from pairs of neighboring subdivisions (e.g., neighboring patches or nearby volumes within an anatomic region) can be computed. These can be averaged within a neighborhood to obtain the average neighborhood maximum correlation, and then across neighborhoods to obtain the threshold . In essence, these thresholds ensure that at each iterative subspace pursuit stage, correlation between the current support and new regions entering the solution is <, preventing new regions from clustering around neighbors of the current support. Rather, it enforces that new regions entering the solution explain substantially different field patterns than the current support and its neighbors can. For a given subject, one threshold was set for each cortical patch decomposition .sup.(1), .sup.(2), and .sup.(3), and used for the successive hierarchies during cortical source space reduction. The final joint cortical and subcortical estimation stage used the minimum of the threshold from the last cortical stage and the group of subcortical volumes. This allows any of the cortical sources in reduced source space .sub.(sp) to enter the final solution. Beyond this, the exact value of the threshold does not matter at the final stage because the sparse subcortical and cortical field patterns are implicitly nearly orthogonal.

(46) The MNE Estimates:

(47) As MNE is an established source imaging technique, it can be computed with known procedures. Specifically, the estimated noise covariance Q is eigen-decomposed to Q=U.sub.Q.sub.Q.sup.2U.sub.Q. Then, the measurement equation is whitened by premultiplying with Q.sup.1/2=.sub.Q.sup.1U.sub.Q. This gives the whitened measurement Y=Q.sup.1/2Y and whitened forward model {tilde over (G)}=Q.sup.1/2G. Then, a computationally stable rewrite of the MNE estimate is:

(48) ${\hat{X}}^{\mathrm{MNE}}\left(\tilde{Y},\tilde{G},R\right)=R{{\tilde{G}}^{}\left(\tilde{G}R{\tilde{G}}^{}+I_N\right)}^{-1}\tilde{Y}$
where N is the number of sensors. Computing these current estimates requires specification of the prior covariance matrix R. In practice, R is unknown, and is usually written in terms of a regularization parameter as {tilde over (R)}=R/.sup.2, giving:

(49) ${\stackrel{.Math.}{X}}^{\mathrm{MNE}}\left(\tilde{Y},\tilde{G},\tilde{R}\right)=\tilde{R}{{\tilde{G}}^{}\left(\tilde{G}\tilde{R}{\tilde{G}}^{}+{}^2{I}_N\right)}^{-1}\tilde{Y}$
Given the SNR of the data, we set .sup.2=1/SNR and assume {tilde over (R)} is diagonal. This gives a convenient choice for the elements of {tilde over (R)}:

(50) ${\tilde{R}}_{\mathrm{ii}}=\frac{\mathrm{TR}\left(I_N\right)}{\mathrm{TR}\left(\tilde{G}{\tilde{G}}^{}\right)}$
where TR refers to trace of the matrix in brackets. Intuitively, high SNR values enforce source current amplitudes to explain peaks in the data better than the noise. By contrast, low SNR values weight the source covariance more and allow larger discrepancy between measured and predicted data. Standard settings of SNR in the range of 9-25 were used for MEG evoked response recordings where the goal is to capture peaks better than background noise. For EEG recordings, as forward solutions have greater blurring due to the relatively low conductivity of the scalp and skull, greater information is contained in background below peaks, and thus SNR was set to 1. The same SNR value was used through the source space hierarchies.

(51) The subspace pursuit algorithm used to refine and identify a sparse subset of MNE estimates for a given source space is denoted as SP (Y, G, Q, R, L)SP ({tilde over (Y)}, {tilde over (G)}, {tilde over (R)}, L), and is performed as follows:

(52) Initialization:

(53) 1. Correlation threshold: Compute a correlation threshold mutual coherence) specifying the degree of orthogonality required amongst the selected columns of {tilde over (G)}. This is set as the average worst case correlation between the forward solutions of neighboring cortical patches in the source space under consideration.
2. Support: Compute the standard MNE estimates {circumflex over (X)}.sup.MNE ({tilde over (Y)}, {tilde over (G)}, {tilde over (R)},) as specified, and the l.sub.2-norms across time for each row of {circumflex over (X)}.sup.MNE. Specify initial support H.sup.(0) as the set of L rows in {circumflex over (X)}.sup.MNE with largest l.sub.2-norms satisfying the degree of orthogonality specified by .
3. Residual: F.sup.(0)={tilde over (Y)}{tilde over (G)}(H.sup.(0)) {circumflex over (X)}.sup.MNE({tilde over (Y)},{tilde over (G)}(H.sup.(0)), {tilde over (R)})

(54) Iteration starting at l=1:

(55) 1. Support Expansion: H.sup.(l)=H.sup.(l1){L rows in {circumflex over (X)}.sup.MNE (F.sup.(l1), {tilde over (G)}, {tilde over (R)}) with largest l.sub.2-norms satisfying the degree of orthogonality specified by }.

(56) 2. Estimation on Expanded Support: Z.sup.(l)={circumflex over (X)}.sup.MNE(, (H.sup.l)),{tilde over (R)}).

(57) 3. Support Trimming: Update H.sup.(l) to L rows in Z.sup.(l) with largest l.sub.2-norms satisfying the degree of orthogonality specified by .

(58) 4. Residual Update: F.sup.(l)={tilde over (Y)}{tilde over (G)}(H.sup.(l)){circumflex over (X)}.sup.MNE({tilde over (Y)},{tilde over (G)}(H.sup.(l)),{tilde over (R)}).

(59) 5. Stopping Criterion: If H.sup.(l)=H.sup.(l1), set outputs H=H.sup.(l) and {circumflex over (X)}.sub.H.sup.(end)={circumflex over (X)}.sup.MNE(Y,{tilde over (G)}(H.sup.(l)),{tilde over (R)}) and end iterations.

(60) The final H specifies the subset of regions in the source space that are most relevant to the measurement, while {circumflex over (X)}.sub.H.sup.(end) specifies the currents in modes of the forward solutions of these regions. Note that the {circumflex over (X)}.sub.H.sup.(end) is unit-less as the forward solutions are scaled by current strengths. To interpret the regional currents on the same units scale and in physical dipole coordinates, all {circumflex over (X)}.sub.H.sup.(end) were scaled by the current strength c and projected back into dipole space to obtain {circumflex over (X)}.sub.H=cW{circumflex over (X)}.sub.H.sup.(end). Overall, this subspace pursuit algorithm offers an efficient means to iteratively concentrate the MNE into the few regions most relevant (in the minimum l.sub.2-norm sense) to the measurement.

(61) Data Analysis Across Hierarchies

(62) The above MNE+subspace pursuit algorithm is repeated across hierarchies to yield estimates {circumflex over (X)}.sub.H.sub.J at the final hierarchy. In practice, forward solutions can be computed across the hierarchies of source spaces and grouped by region. Time windows of interest can be selected based on latencies of expected features. The number of relevant regions or the target sparsity level L can be specified to reflect the spatial extent/number of regions expected to be active during a given paradigm. The same L can be maintained across hierarchies. The number of relevant regions can be pre-specified for convenience, or the algorithm could be rewritten to iteratively add one projection at a time until all significant field patterns in the measurement are explained and the residual is white across sensors. The procedure is very similar to employing a formal model selection and goodness of fit procedure. Both MEG and EEG recordings can be used for localization if simultaneous recordings are available. The two sets of recordings can be processed separately (filtering, averaging and covariance estimation), their forward solutions can be computed separately, and the measurement vectors and forward solutions can be concatenated after accounting for respective noise covariances separately for consistency of units. The final source current estimates {circumflex over (X)}.sub.H.sub.J specifies the currents in modes of the forward solutions, and can be projected back into dipole space for reporting. Then, the estimated source currents can be summarized as vector resultant currents by region. Spatial brain activity maps, regional source current time courses and summary root mean square values of dipole source current estimates can be computed and displayed for understanding the regional distribution of current estimates.

Experimental Results

(63) To validate the algorithm, auditory responses evoked with a train of click stimuli were analyzed during resting eyes open condition. Auditory responses were chosen because these responses comprise distinct M/EEG peaks and established latencies corresponding to a stereotypical progression of activity from the cochlea, through inferior colliculus to medial geniculate thalamus and auditory cortex, and thus serve as a good testcase for validating a subcortical source localization algorithm. Simultaneous M/EEG auditory evoked recordings (AEPs) were obtained during binaural stimulation and corresponding anatomic MRI using the paradigm described below.

(64) Simultaneous MEG and EEG AEP recordings, and structural MRIs were obtained on two healthy volunteers aged 25-45 years screened for standard MRI contraindications and normal audiometry (no evidence of hearing loss in 0-4 kHz range). A train of broadband clicks having 0.1 msec duration, intensity 65-80 dB/nHL and inter-stimulus interval 110 msec (click rate 9.09 Hz, corresponding to highest AEP SNR), generated within the Presentation software (Version 17.1, Neurobehavioral Systems, Inc., Albany, Calif., USA), was delivered binaurally during eyes open resting condition. Subjects were asked to sit still, not pay any attention to the sounds and imagine a dot at the center of the screen. M/EEG data were recorded at 5 kHz sampling with filter cutoffs set to 0.03-1660 Hz. After acquiring 2 min of pre-stimulus baseline eyes open data, AEPs were recorded in 5 runs of 5.5 min each, yielding 10000-16000 epochs for averaging. At the start of each run, standard checks were performed to ensure no stimulus artifacts, subjects were allowed a blink break if needed, and HPI recordings were obtained. Non-stimulus baseline recordings were repeated at the end of the study and MEG empty-room recordings were also obtained. The study was approved by the Partners Human Research Committee at Massachusetts General Hospital, and written informed consent was obtained from all subjects.

(65) The raw data were preprocessed to remove power line noise using a comb notch filter (MATLAB) with 30 notches at harmonics of 60 Hz, each having bandwidth 1 Hz. Artifactual channels (marked by inspection) and eye-blink epochs (peak to peak EOG<150 uV in 1-40 Hz band) were both excluded. The early auditory brainstem response (ABR) and later middle latency response (MLR) components of the auditory evoked potential were processed separately by band-pass filtering the preprocessed data to 500-1625 Hz, and 30-300 Hz respectively. In each case, the data were processed to (a) compensate for interference, (b) obtain stimulus locked average evoked responses, and (c) estimate noise covariances. First, standard signal space projections (SSP), computed using principal component analysis on the empty room MEG recordings, were used to compensate for environmental interference. Second, stimulus timings to denote epochs for averaging were compensated for a sound tube delay of 9.775 msec (as measured by recording sounds delivered at the ear piece with an optical microphone). Grand averages across runs were performed without HPI correction if HPI head coordinates across runs were within 2%. Overall, a total of 11200 epochs were used for the averages considered in our analysis. Third, observation noise covariances were estimated using the baseline eyes open recording. These data give a good measure of background noise arising from ongoing brain activity and systematic instrumental disturbances. As the MEG noise covariance is typically ill-conditioned, external disturbances to the covariance estimate were suppressed by applying the SSP operator, and then diagonal loading was applied. All SSP, averaging and covariance calculations were performed in bands relevant to components of interest.

(66) FIGS. 7A-B summarize the M/EEG evoked responses. The early auditory brainstem responses (ABR) are shown in FIG. 7A, while later middle latency responses (MLR) are in FIG. 7B.

(67) Examining FIG. 7A, we note that the EEG ABRs have higher SNR than MEG ABRs. Thus marking the early peaks on EEG and looking for coincidental peaks in the MEG recordings, we see that the MEG magnetometers also show ABRs at the expected latencies while the gradiometers do not exhibit significant features at these latencies. While it is difficult to resolve individual early waves corresponding to cochlear nucleus and lower brainstem structures like the superior olivary complex, there is clear evidence of evoked peaks in the 5-7 msec time range, which is consistent with the expected timing (5.6-5.8 msec) of the brainstem wave V known to arise from inferior colliculus. Proceeding to the MLR in FIG. 3-8B, we see that the EEG has a low amplitude Po feature, which is thought to mark the end of brainstem components and potentially the start of thalamic components. However, Po is not seen clearly on the MEG and the thalamic involvement at these latencies is controversial and has not been resolved to date. Peaks Na-Pa reflect classic middle latency activity in the (supratemporal) primary auditory cortex (A1). Peaks beyond 40 msec reflect activity in secondary auditory cortices and higher association regions. Thus, as the most prominent peaks with well-known locations are the early brainstem wave V and the later Na-Pa peaks, we asked if source localization can reveal activity in (a) inferior colliculus during 5-7 msec, and (b) A1 during 20-30 msec.

(68) FIG. 8 illustrates the spatial maps and time courses of source estimates obtained on the auditory evoked responses. FIG. 8A shows localization of MLR peaks (Na-Pa) to bilateral primary auditory cortices (A1). These areas comprise the reduced cortical source space Csp, which along with the subcortical volume sources S, form the combined source space J for final estimation. We focus on describing estimates for the early ABR components where subcortical structures are activated. FIG. 8B-D illustrates spatial maps and time courses of the estimates within source space J. Subspace pursuit estimates (FIG. 8B) localize specifically to inferior colliculus, while MNE estimates (FIG. 8C) localize both to inferior colliculus and sources in caudate and thalamus (with the latter estimates having higher amplitude than the former). FIG. 8D quantitatively summarizes the regional distributions of dipole current amplitudes by regionshowing that SP refines the non-specific MNE estimate to concentrate activity into the most relevant regions. Time courses of the SP and MNE estimates are in FIG. 8E-F. SP estimates of activity in inferior colliculus peak at the time interval where wave 5 is expected. As the number of relevant regions L was set to 1 (only 1 key region is expected to be active during 0-10 msec), the sparse estimate shows only one time course. The MNE estimates show activity in multiple regions peaking at the time intervals of interest, making it hard to discern specific activity expected in inferior colliculus.

(69) Finally, all the above estimates were obtained using both MEG and EEG data (where EEG has higher SNR wave 5 peaks). Accuracy was also effectively validated on unimodal MEG data, which has lower SNR, showing that the algorithm does not rely on having high amplitude peaks or high SNR in the data.

(70) Although specific embodiments are described above, it will be apparent to those of ordinary skill that a number of variations can be made within the scope of the disclosure. For example, it will be apparent that processing equipment and hardware can be used to perform the data acquisitions, and data analysis methods, processes, and algorithms described above. Computing equipment for performing these tasks can include a processor; various types of memory components; user input devices, such as keyboards, mice, and data input drives; and output devices such as displays, printers, and data transmission devices. These components can be provided as part of a medical imaging system, in a network connected to a medical imaging system, as a stand-alone system, or as part of a network. Although flow charts illustrating steps are provided and described above, these charts are by way of illustration. The steps shown in any given flow chart are not required in all cases, and the order of the steps can, in some cases, be varied. It should be understood, therefore, that the methods and apparatuses described above are only exemplary and do not limit the scope of the invention, and that various modifications could be made by those skilled in the art that would fall within the scope of the invention. To apprise the public of the scope of this invention, the following claims are made: