1. Introduction

While fMRI (functional Magnetic Resonance Imaging) yields high spatial resolution, brain dynamics are hardly resolved. This is due to two facts: First, there is a trade-off between signal-to-noise ratio (SNR), spatial, and temporal resolution of the imaging process. Secondly, there is an inherent low-pass filter in the event chain that leads from neuronal activity to hemodynamic reactions as they are measured. The combination with Electroencephalography (EEG) and Magnetoencephalography (MEG) source reconstruction techniques promises to add the desired temporal resolution. EEG and MEG sampling times are usually in the order of a millisecond, and the effects of the neuronal activity that are measured (electric potentials and magnetic fields) are of an instantaneous nature. They can be modeled using the quasi-static approximations of Maxwell’s equations. However, the localizing power of EEG or MEG alone is limited and typically in the order of one centimeter.

To achieve a combination of these modalities, the results of an fMRI analysis shall be used as constraints in the subsequent EEG/MEG source reconstruction. Doing this, one has to keep in mind that there are differences in the anatomical phenomena and in the necessary experimental setups (a single paradigm vs. statistically significant differences with respect to at least two paradigms). As a result, locations of significant hemodynamic differences (which we shall call “hotspots” subsequently) must not simply be equated with current sources. Instead, several problematic constellations may occur:

Current sources that do not show up in the fMRI images. This might happen, if neurons are not active long enough to cause increased blood flow.

fMRI hotspots without a corresponding current source. This situation could occur in the case of so-called silent sources, i.e. sources or source constellations that produce no measurable fields or potentials.

A slight displacement between fMRI hotspot and source location. This is to be expected due to the different locations of neurons and involved vessels, as well as due to the fact that an extended cortical patch, if curved, may not even contain the representative source location.

In this paper, we propose two strategies for the integration of fMRI with source reconstruction techniques that address these difficulties. One strategy has been developed for the dipole fit methods and another for distributed source models.

Furthermore, one has to deal with the co-registration problem between the imaging modalities involved [van den Elsen PA 1993]. Depending on the information available, several techniques could be used for finding the best match between sMRI (structural Magnetic Resonance Imaging) and fMRI image coordinates:

·      landmark-based registration. Here, corresponding landmarks have to be determined in both modalities. [Fuchs M et al 1995] (see Figure 1)

·      surface-based registration. Prerequisite is the segmentation of e.g. the skin or the brain surfaces in both modalities.

·      volume-based registration. No additional information is required. [Wells WM et al 1997]

·      scanner coordinate-based registration. Here, sMRI and fMRI must have been acquired during the same session, and the images must include position information in scanner coordinates.

For the surface-based registration and the scanner coordinate-based registration, we describe how a reliable match can be achieved.

Figure 1. Co-registration of MRI and CT data using anatomical landmarks (nasion and preauricular points).

2. Surface-Based Image Modality Co-Registration

The surface registration method is based on a contour fit algorithm [Huppertz HJ et al 1998,Willemsen S et al 1998] that matches a set of N digitized contour points r_i with a segmented surface such as the skin or the brain. The algorithm minimizes the L₁-norm of the distances of all contour points to the surface. These distances are determined efficiently by using a distance map D(r), which can be pre-computed from the segmented surface (see Figure 2).

with

(1)

The L₁-norm (sum of absolute values) is more robust with respect to outliers than the commonly used L₂-norm (sum of squares). In order to deal with local minima, three translation parameters (shifts) and three rotation parameters (Euler angles) are optimized independently, such that the best-fit translations are found for each set of rotations. This is done using two nested Nelder-Mead simplex algorithms. The translation parameters are initialized with the difference between the center of mass of the segmented surface and the center of a sphere fitted to the set of points. The rotation parameters are initialized by a coarse exhaustive search over the parameter space of the Euler angles with optimized translations.

with

Figure 2. A 3D distance map. Increasing intensity encodes increasing distance from the segmented surface.

For the surface-to-surface registration, one surface is converted into a set of contour points. These contour points must not include areas, which are prone to ambiguities or missing correspondences, such as the ears, the neck, or the chin. Subsampling (thinning) is applied in order to decrease numerical complexity, resulting in about 300 contour points. As one surface may be overall larger than the other, e.g. due to the segmentation thresholds used, the average distance of all contour points from the surface is not taken into account in this case, and the problem is now to

with

(2)

Results are shown in Figure 3.

Figure 3. Interlaced display of co-registration results obtained from segmented skin surfaces. The first modality is a 3D sMRI T₁ scan with 250 mm Field of View; the second modality is a gross average 3D T₁ data set normalized to Talairach coordinates with 230 mm Field of View. The two surfaces used for registration are shown as white outlines.

3. Image Modality Co-Registration using Scanner Coordinates

The scanner coordinate-based registration method makes use of the Image Position and Image Orientation tags, which are part of the ACR-NEMA standard and can be found in DICOM images. The Image Position is a coordinate vector in [mm]-space that describes the position of the upper left corner of the image with respect to the scanner. The Image Orientation are two normal vectors in the same coordinate system. They describe the orientations of the x and y image axes, respectively. The coordinate system thus defined is tracked down the pipeline of image transformations (shift, zoom, rotation, mirroring, stacking), which are applied when scanner images are read in and can finally serve for co-registration between modalities. The scanner coordinate-based co-registration can be computed on the fly when ACR-NEMA or DICOM images are read in. Results are shown in Figure 4. If scanner coordinates are available, their use for co-registration purposes is practical and straightforward.

Figure 4. Co-registration using scanner coordinates. The first modality is a 3D sMRI T₁ scan with 256 mm Field of View, the second modality are 5 T₂ slices with 250 mm Field of View and a slice-to-slice distance of 4 mm that have been acquired during the same session. fMRI acquisition has been performed for the same slices.

4. Volume Conductor Modeling

The accuracy of EEG and MEG source localization with respect to source locations has to be in the order of few millimeters, if fMRI priors shall successfully be used. This can – for sources in most parts of the head – only be achieved by using realistic volume conductor models. Two popular modeling techniques are the BEM (boundary element method) and the FEM (finite element method). While the BEM [Fuchs, M et al 1998] uses surface triangulations of the borders between compartments of equal isotropic conductivities as a geometric model, most FEM [Buchner H et al 1997] implementations work on tetrahedral volume tessellations, where each tetrahedron may be assigned an individual, anisotropic conductivity.

The generation of these geometric models is a nontrivial task, because the segmentation of the compartment borders is not straightforward, and certain minimal mesh element sizes, smoothness constraints, and inter-compartment distances have to be observed in order to avoid instabilities of the algorithms. An automatic procedure can be used to generate BEM or FEM meshes containing the brain, skull, and skin compartment borders (Figure 5) from sMRI. It speeds up model generation and ensures model quality. The whole model definition takes some three to four minutes on a PIII/600 PC and comprises the following steps:

1.    Cortex segmentation. The cortex forms the basis of the model. Segmentation is based upon a 3D region-growing algorithm. The segmented white matter, which is smoothed and enlarged, serves as a maximum volume containing the cortical structures. The cortex may also be used for constrained source reconstructions.

2.    Inside of the skull. As the inside of the skull does not give reliable contrast in T₁-weighted images, it is derived as the smoothed and slightly enlarged shape of the cortical envelope. This guarantees a minimum distance between cortical sources and innermost compartment border.

3.    Outside of the skull. In order to ensure a minimum inter-compartment distance, the smoothed union of the segmented skull-skin transition and the enlarged inside of the skull is used as the outside of the skull.

4.    Skin. In a similar way, the segmented skin and the enlarged outside of the skull define the skin compartment. The outside of the skull, strongly enlarged, also serves as a maximum volume for the skin. This avoids the inclusion of the difficult to model and less important basal regions into the skin compartment.

5.    Mesh generation (Figure 6). For the BEM, the three compartment borders are triangulated, using individual triangle sizes. For the FEM, a tetrahedral tessellation with regularly shaped tetrahedra is performed.

 

Figure 5. Compartment borders (white) overlaid onto MR slices.

a)b)

Figure 6. a) BEM meshes; b) FEM mesh.

5. fMRI-Constrained Dipole Models

Dipole models are parameterizations of the brain activity that use few current dipoles (point sources), whose number, locations, components, and temporal characteristics have to be determined. One dipole is interpreted as an abstraction for an extended patch of cortical tissue that has a predominant orientation. The rationale for this abstraction is, that a single neuron can be modeled as a microscopic current dipole. If location, orientation, and temporal activation of neurons are correlated, they produce a measurable far field, which resembles the one of a current dipole at the same location. A class of neurons for which all of these necessary correlations are certainly given are the pyramidal cells in the cortical gray matter with their aligned orientations perpendicular to the cortical sheet.

Usually, number and temporal characteristics are assumed to be known a priori, while the remaining parameters are fitted. Common temporal characteristics include the moving dipole model, where independent locations and components are determined for each time point, the rotating dipole model, where the location stays constant over a latency range, and the fixed dipole model with constant location and orientation and a time-varying strength only.

This leads to an optimization problem such as

(3)

Here, m are the measured data, j are the unknown dipole components for given locations r, and L is the (linear) lead field matrix for given locations r. The lead field matrix describes the relation between dipole components and measured data as defined by the conducting properties of the head. m and j may comprise several time points. For given locations r, D and thus j can be computed non-iteratively using linear algebra, while an iterative solver has to be used for finding the optimal r.

For fMRI-constrained dipole models, rotating or fixed dipoles are associated with each fMRI hotspot. To this end, the hotspots are used as seed points for the dipole fit, while a maximum distance constraint applies. Dipole locations, components, and time courses are fitted. In order to account for unexplained data, an additional dipole is fitted, which is not related to any hotspot (Figure 7). A regularized solution is obtained that is an extension to Eq. 3 in the sense that it suppresses source components with small data overlap [Fuchs M et al 1998] and that a penalty term P is added, which assures that the maximum distance criterion is met. The resulting optimization problem is

(4)

Here, C is the sensor weighting resp. SNR normalization matrix. W_overlap is the non-diagonal overlap weighting matrix for given C, m, and L. SNR normalization is an approach that allows for the combination of EEG and MEG, the suppression of noisy channels, and a straightforward assessment of the goodness-of-fit in a c²-sense [Pflieger ME et al]. Due to the noise normalization performed by C, the regularization parameter l can be pre-calibrated with respect to the SNR. The vicinity constraint for N sources and seed points s_i, and maximum distances d_i is imposed using the penalty term

with

(5)

Again, the optimization problem can be solved analytically for given locations r. In the scope of this paper, d_i = 5 mm was used.

Figure 7. Schematic drawing for three hotspots and a total of four dipoles: dipoles are loosely fixed to fMRI hotspots, using the hotspot locations as seeds and constraining the dipoles to stay within a maximum distance from their seed points. One additional dipole is fitted, which is not tied to a hotspot.

6. fMRI-Constrained Distributed Source Models

For fMRI-constrained distributed source models, the cortical sheet segmented from sMRI is used as the source space [Wagner M et al 1995,Dale AM et al 1999,Liu AK et al 1998], resulting in some ten thousand fixed locations and normals with distances of 2 to 3 mm. At each of these cortical locations, a dipolar source is assumed that accounts for the neuronal activity of the small surrounding patch. The orientation of the source is fixed to the pyramidal cell orientation, i.e. perpendicular to the cortical surface.

a)	b)

Figure 8. a) A middle layer of the gray matter sheet used as the source space for current density reconstructions (detail). Current dipoles are computed at each vertex location. Their orientations are normal with respect to the triangulated surface; b) Spatial relation between cortical sheet, fMRI hotspot, and the somewhat larger area of increased location weights.

A depth-normalized, regularized inverse solution with a minimal weighted L₂-norm of the source strengths is computed:

(6)

The first term is the data term. It measures the discrepancy between measured and forward calculated data. Here, m are the measured data, j are the unknown currents, L is the lead field matrix, and C is the sensor weighting resp. SNR normalization matrix.

The second term is the model term. It measures the discrepancy between the computed currents and an implicitly assumed model. Here, W_fMRI is the diagonal fMRI-induced location-weighting matrix, and W_depth is the diagonal depth normalization matrix common to weighted least squares methods. The location weights W_fMRI depend on the fMRI-induced significances. Different weighting functions can be used, including a total suppression of non-hotspot activity, a damping of non-hotspot activity, and arbitrary functions of the fMRI correlations (Figure 9). In the scope of this paper, different weight factors are used, while locations are either identified as hotspots (significance 1.0) or not (significance 0.0). Unconstrained solutions are obtained if W_fMRI = 1.

Figure 9. Relation between fMRI significance and weighting factor for a given location. a) fMRI hotspots are enhanced by a factor of 1.4; b) fMRI hotspots are enhanced by a factor of 2.0; c) fMRI hotspots are enhanced as a function of their significance.

7. A Simulation Study

In order to test the constrained source reconstruction methods, forward calculated 61 electrode EEG data source (Figure 10) were generated for three dipolar sources. Each source showed a depolarization / repolarization sequence as activation pattern. Noise was added to achieve an SNR of 15 for the dipole fits and SNR’s of 30 to 8 for the current density reconstructions. Functional MRI data were simulated with three hotspots. Two of these hotspots coincided with two of the three source locations, while one hotspot had no corresponding source (Figure 11).

 

Figure 10. Simulation setup with 61 electrodes and source plane.

a)	b)
c)

Figure 11. a) Three simulated fMRI hotspots (black); b) Three simulated sources (black). From left to right: sources 1, 2, and 3; c) Source waveforms.

7.1 Dipole Fit Results

In all 4-dipole fits, three dominating sources at the correct locations were found. Sources 1 and 2 interfere slightly in the rotating dipole fits but not in the fixed dipole fits. Results are shown in Figure 12 and Figure 13.

a)	b)

Figure 12. a) Fit results of 4 rotating dipole fit (black) and fMRI hotspots (seeds, gray). A small displacement between sources 1 and 2 and the hotspots used as seeds can be observed, which is due to the maximum distance constraint; b) Rectified time courses of the 4 rotating dipoles. Time courses are all positive; source orientation reverses instead. The fourth source is the one at the third hotspot; it should have no activity at all.

a)	b)

Figure 13. a) Results of 4 fixed dipole fit. A small displacement between sources 1 and 2 and the hotspots used as seeds can be observed, which is due to the maximum distance constraint; b) Time courses of the 4 fixed dipoles. The fourth source is the one at the third hotspot; it should have no activity at all.

a)	b)
c)

Figure 14. a) Results of fMRI-constrained MNLS current density reconstruction (SNR = 15, clipped) and simulated source locations (black); b) Results of unconstrained MNLS current density reconstruction (SNR = 15, clipped) and simulated source locations (black); c) Scale.

7.2 Current Density Results

In the fMRI-constrained Minimum Norm Least Squares (MNLS) current density solutions, all source locations were retrieved when using a weight factor of 1.4, even in the case of a source without corresponding hotspot. However, the time courses of sources 1 and 2 interfere. Results are shown in Figure 14 a). For comparison, the results of an unconstrained solution with W_fMRI = 1 are shown in Figure 14 b).

Over the range of weight factors between 1.0 (unweighted) and 10.0 (strongly weighted), a transition from unconstrained results (data driven) to a mere reflection of hotspot locations (constraint driven) can be observed. For all methods and SNR’s, however, a weight factor of 1.4 turned out to be the best choice: the solution would benefit from the fMRI constraints in such a way that sources at hotspot locations were focalized, while the additional source was still retrieved. In Table 1 to Table 5, entries with a weight factor of 1.4 are emphasized.

For SNR’s between 30 and 8, the fMRI-constrained MNLS method (Table 1 to Table 3) was able to correctly localize all three sources without showing false maxima.

The Minimum L₁ Norm method [Fuchs M et al 1999, Wagner M et al 1998] (Table 4) generally exposes less smearing artifacts than MNLS. When used together with fMRI constraints, activity is focalized to hotspot locations, while the third source is still being localized.

LORETA [15,16] (Table 5) incorporates lead field weighting differently than minimum norm methods, as weights are applied to local smoothness measures rather than to the source strengths themselves. While focalizing to hotspot locations, LORETA was not able to retrieve the third source, which is only revealed in the unconstrained solutions.

Table 1: Minimum Norm Least Squares, SNR = 30. For all MNLS reconstructions, fMRI constraints with weights 1.4 focalize sources with corresponding hotspot, while the source without corresponding hotspot is still retrieved.

Table 2: Minimum Norm Least Squares, SNR = 15

Table 3: Minimum Norm Least Squares, SNR = 8. Even for this relatively low SNR, all three sources are found. The unconstrained solution is extremely smeared due to the high level of noise in the data.

Table 4: Minimum L₁ Norm, SNR = 15. For L₁ norm reconstructions, fMRI constraints with weight 1.4 focalize sources with corresponding hotspot, while the source without corresponding hotspot is still retrieved. The unconstrained solution shows artifacts in the cerebellum, which is typical for very sensitive measures like the L₁-norm.

Table 5: Cortical LORETA, SNR = 15. fMRI constraints focalize sources with corresponding hotspot, but the source without corresponding hotspot is not reconstructed.

 

8. Discussion

Depending on the available information, different co-registration methods for merging sMRI and fMRI coordinate systems may be useful. The usage of scanner coordinates seems to be the most straightforward method, if both modalities have been acquired during the same session and the necessary information is part of the image file format. If this is not the case or not possible, other methods have to be applied. Landmark- and surface-based registrations both leave the user with the responsibility of choosing matching landmarks resp. surfaces. Volume-based registration is independent of such influences, but about an order of magnitude slower than the surface-based approach.

Dipole fits can benefit from fMRI constraints: Meaningful seed points for the location fit are obtained, a crucial issue in all multi-dipole fits. A reconstructed dipole in the vicinity of each fMRI hotspot yields the corresponding source time course. Spatially unconstrained dipoles are then necessary to account for remaining activity. However, the temporal discrimination of close sources can only be achieved by an additional fixed dipole constraint.

Current density reconstructions react upon fMRI constraints in two ways: Activity in the vicinity of fMRI hotspots is bundled. Remaining activity can be localized correctly, if its field distribution cannot be generated from sources within the hotspots, and if the weight factor is not too large. The temporal discrimination of close sources is enhanced with respect to unconstrained reconstructions, if the correct source orientations are dominating within the hotspots.