1.Introduction

New insights into the normal functions and pathologies of the human brain have been provided by the emergence in the last two decades of novel non-invasive neuroimaging techniques ([Belliveau 1991;Kwong 1992;Ogawa 1992;Pascual-Marqui 1994;Valdés-Sosa 2000]). Current scanning strategies in functional Magnetic Resonance Imaging (fMRI) and newer-generation of Positron Emission Tomography cameras (PET) now can dynamically scan the entire brain with high sensitivity, allowing for mapping of hemodynamic measures as well as those of neurotransmitter and receptor uptake and regulation. Current hemodynamic measurements, particularly functional MRI (fMRI), provide excellent spatial resolution (millimeters). Unfortunately, fMRI is limited in its temporal resolution by the latency of the hemodynamic response (in the order of seconds) ([Kim 1997;Menon 1998]), curtailing the ability of this technique to following the high frequency components of temporal evolution of the neural activity underlying brain functions. Conversely, ElectroEncephaloGraphy/MagnetoEncephaloGraphy (EEG/MEG), provide excellent temporal resolution (in the order of milliseconds), but uncertain spatial localization of brain generators ([Hamalainen 1993]). The stage is clearly set for the development of new techniques with both high spatial and high temporal resolution obtained by combining information from electrophysiological and hemodynamic imaging modalities ([Ahlfors 1999;Heinze 1998;Heinze 1994;Ives 1993;Korvenoja 1999;Luck 1999;Opitz 1999;Snyder 1995;Warach 1996]).

The basic problem for electrophysiological imaging is the ill-posedness of the EEG/MEG inverse problem, i.e. for sufficiently general configurations of primary current there is no unique solution. An approach to solve this non-uniqueness, adopted in previous work, has consisted in constraining the possible sites for EEG/MEG sources to be limited to those in which significant activation is found in comparable fMRI experiments ([Dale 2000;George 1995;Liu 1998]). We shall denote this type of procedure as fMRI “constrained EEG/MEG inverse solution” (fMRI EEG/MEG), to emphasize the asymmetrical nature of the inference, i.e. spatial information is provided essentially by the fMRI and temporal information by the EEG/MEG. While straightforward conceptually, the fMRI EEG/MEG approach suffers from a number of drawbacks:

• The fMRI is considered as the “truth” for spatial information notwithstanding the extremely low signal to noise ratio of this imaging modality. In fact in ([Dale2000]) it is recognized that serious bias may occur when an actual EEG/MEG source is not detected by the hemodynamic response.
• As a consequence spatial electrophysiological information is not used in inference.
• A more serious criticism is that there is no principled way of estimating the various “tuning parameters” that are involved in the model.

This paper proposes an alternative procedure, based on a fully Bayesian Model that entails formulating a “Forward Problem” for both the electrophysiological and hemodynamic measurements. It is a specification of a more general formulation stated for the first time in ([Valdés-Sosa 1998]). Since both types of measurements contribute to the localization of “active” brain regions, we shall denote this type of procedure as “Symmetrical fMRI-EEG/MEG Fusion” (fMRI EEG/MEG). The relative weight given to either type of imaging modality is determined by the data itself. An additional advantage of this approach is that all “tuning Parameters” are estimated from the data.

This paper is organized in as follows. In Section 2 a general formulation of fMRI-EEG/MEG image fusion is set out, following ([Valdés-Sosa1998]) closely. Different possible approaches to image fusion are then discussed and the reasons for choosing a Symmetrical model is explained. In Section 3 the Symmetrical fMRI-EEG/MEG model is fully specificied and procedures for the estimation of parameters are described. In Section 4 simulation results are presented that indicate that the new proposed here has substantial gains in performance according to a number of measures of quality of Tomographic methods based on the Point Spread Function (PSF). Section 5 describes the application of the fusion method to actual data, a joint MEG-fMRI Somatosensory experiment. Finally in Section 6 limitations of the present model are discussed and future directions of work outlined.

In all that follows the following conventions will hold. Scalars, vectors and matrices will be denoted by lower case, bold lower case and bold uppercase typefaces respectively. When defining a vector or matrix the dimensions will be expressed in a subscript, e.g.  indicates that the matrix  has r rows and c columns.  will denote the identity matrix of order N. For the vector ,  will denote its transpose. For a matrix , , will denote its determinant and trace, respectively. The matrix  is the diagonal matrix with the elements of  on the main diagonal.  will denote the column vector formed by stacking all the smaller column vectors of . Probability densities will be expressed according to the notation in Table 3. For example the r-dimensional multivariate normal density with mean  and covariance matrix will be denoted as . We shall denote an independent identically distributed (i.i.d.) random variable from this distribution is denoted as .

2. The General fMRI-EEG/MEG Fusion Problem.

In general, data fusion is the process by which data from a multitude of sensors are used to yield an optimal estimate of a specified state vector pertaining to the observed system ([Mohamad 1997;Richardson 1991]). A natural framework for combining all this information about the statistical properties and interdependence between the variables included in the model is Bayesian Inference Theory. In a Bayesian estimation paradigm, the a priori knowledge about the parameters of interest is included in the form of prior density distributions. The measurement process modifies this a priori information through the specification of the likelihood terms. As a result, all the necessary information for the estimation of our model is included in the posterior distribution of the model parameters, when the data is given. Formally this is stated as Bayes’ Theorem :

Where D indicates the data for the problem and H the model to be estimated. The term  is known as the “likelihood” and  is the “a priori” probability density. In Hierarchical Bayesian Modeling the likelihood and a priori are decomposed into the product of simpler, known densities. This is the procedure to be followed below.

For the particular case of fMRI-EEG/MEG fusion this general formulation is now interpreted as follows.

2.1 Forward Problems for the EEG, MEG, and fMRI: Likelihoods

The likelihoods for the observed data are obtained by consideration of the corresponding forward problems. Without loss of generality we will assume that there is no time correlation. This is an unrealistic assumption in practice. However this complication will be omitted for the sake of clarity of exposition.

For electrophysiological measurements the voltage/magnetic field recorded in the sensor array distributed over the scalp (EEG/MEG) is related to the primary current density (PCD) inside the head by a Fredholm integral equation of the first kind ([Riera J. 1998]). Once a grid of points is selected inside the brain for the solution of the inverse problem the resulting discretization of the problem yields a system of linear algebraic equations that is represented as:

   

where,  are the vectors of measurements in the  EEG electrodes and  MEG sensors respectively for each time instant ; and are the vectors of measurement errors for the EEG and MEG sensor respectively. ,  are the electric and magnetic lead field matrices multiplied by the known vectors of orientations of the current density for each point in the gray matter. In other words it has been assumed that the orientation of the current density vector for each possible source is constant in time and directed perpendicular to the cortex. Thus, only the strength and sign of neuronal PCD  at the  points in the solution grid will be estimated. A further assumption is that,

 i.i.d

where the hyperparameters ,  are the sensor noise variances. From these expressions it follows that the likelihoods for the electrophysiological measurements read:

The other source of information is the fMRI time series. For simplicity we assume that the spatial sampling of this neuroimaging modality is matched to that of the solution grid for the electrophysiological inverse solution. Here it is assumed that, an unobserved fast process, the “Net Metabolic Demand” (NMD) for all voxels at time t, , is convolved with a known hemodynamic matrix response function [Friston 1995a;Friston 1995b]:

   

or, in a more compact notation:

with

.

The actually observed fMRI time series,, is a sub-sampled version of  (since the acquisition rate of the MRI scanner is much lower than that of electrophysiological measurements). This is formalized by use of the sub-sampling matrix

where   is the number of fMRI time samples. We shall assume that the fMRI measurement errors are distributed as

 i.i.d.

which leads to the following likelihood for the fMRI measurements:

 

2.2 Mechanisms of Data Generation: Posterior Density for Fusion

The basic assumption underlying our image fusion approach is that all observable quantities are a consequence of neural activation. Let us denote the net amount of neural activation for all voxels at time t as . We shall assume that both the PCD and the NMD at each voxel are generated by some deterministic or stochastic mechanisms  and  such that:

where the symbol “” denotes operations still to be defined. Based on these arguments, and taking into consideration the expressions for the likelihoods for the EEG, MEG and fMRI defined in the previous section, we now set down a general model for Hierarchical Bayesian fMRI-EEG/MEG Fusion which reads as follows:

Note that this model requires in addition to the likelihoods, a priori distributions for the PCD, for the activation as specified in as well as for all the hyperparameters by means of the densities , and . Estimation of all these parameters may be carried out by any of the usual methods of Bayesian computation.

2.3 Possible Solutions Based on Detailed Physiological Knowledge.

Knowledge is accumulating about the physiological changes of the brain that accompany neural activation, for example, elicited by a given stimulus, a modest increase in the cerebral metabolic rate of oxygen (CMRO₂) is accompanied by a much larger increase in local blood flow ([Frostig 1990;Grinvald 1991;Malonek 1996;Malonek 1997]). Because of this imbalance, local capillary and venous blood are more oxygenated during activation. The Primary Current Density generated in the extra cellular medium, by the neural activation is the basis for the EEG measurements (millisecond range), and the decrease in local deoxy-hemoglobin concentration that take place in the venous system, is the basis for functional Magnetic Resonance Imaging (second range), exploiting the blood oxygenation level dependent (BOLD) effect. If detailed mathematical models were available for  and as well as , then the specific formulation of model could be relatively straightforward. Some advances have been made in terms of formulating the pair of equations (in our notation)  and  [Friston 2000]. While of great promise this line of modeling is not yet complete for practical application.

2.4 fMRI Constrained EEG/MEG Inverse Solutions (fMRI EEG.).

This type of fusion may be reanalyzed in the framework of model . In fact, this model is equivalent to simplifying, for the EEG case, the full model to:

It should be noted that the presence of an activation factor is ignored and therefore inference is not symmetric.

3. Symmetrical fMRI-EEG/MEG Image Fusion (EEG ↔ fMRI).

3.1 Complete Model Specification

For initial exploration we shall assume the simplest possible model for the fMRI forward problem matrix  that is equivalent to assuming that:

i.e. that the hemodynamic measurement is just a summarization of NMD at certain critical time instants. This might be more appropriate for PET measurements than fMRI.

Additionally, we shall assume the following model for the relation between PCD, NMD and neural activation:

with  a constant introduced to match the units of PCD and NMD,  and  denoting i.i.d.  variables and   the Hadamaard product. As a consequence of this specification the a priori densities of the PCD and NMD are:

defining

and

Continuing the hierarchical specification, the a prior for  is taken as an Inverse Wishart:

where  is a scaling factor, and  is are the degrees of freedom In order to completely specify our model, we include non-informative inverse Gamma distributions for the hyperparameters:

, , , ,

3.2 Estimation Procedure

The Log-posterior probability distribution of all the parameters when the data is given, takes the form:

We are looking for the set of parameters  that maximize . These are found by the Iterated Conditional Maximization method for finding the Maximum A Posteriori (MAP) estimators for the parameters. We iterate between the following expressions:

with

Where, ,  and  are the elements of , the diagonal of ,  and , respectively, with i = 1..Ng.

An intuitive interpretation of these estimators is that activity is estimated as a linear combination of PCD and NMD neuroimages, and in turn, influences their estimation. Note that the estimation of all hyperparameters comes from the data itself.

4. Simulation Results

In this section we report the results of comparing measures of quality of the Tomographic methods derived from the Pont Spread Function (PSF) of Low Resolution Electric Tomography (LORETA, [Pascual-Marqui1994]), EEG fMRI, and fMRIEEG. This was carried out for the grid of points inside the head that represents the possible generators (i.e. gray matter) for which the PCD is computed. In the simulations the same spatial resolution for EEG/MEG Inverse Problem, and fMRI is assumed, and a grid of 3433 points equally spaced inside the gray matter is constructed, matching the spatial locations of the voxels in the fMRI. The simulation consisted in generating the EEG/MEG and fMRI data, using the forward problem equations (1) and (6). The PCD and NMD used correspond to the activation of a single unit magnitude dipole located in each point of the grid with a Gaussian variation in time (Figure 1). This curve was discretized at five time points were simulated. A measurement noise level of 10% of the peak-to-peak voltage/fMRI signal is assumed in the computations. The PSF in each point, was computed at the time instant where the signal to noise ratio is maximum in both EEG/MEG and fMRI. All measurement units have been omitted for the sake of simplicity.

Based on the PSF for both EEG fMRI and LORETA, some quality measures for characterizing electrophysiological Tomography techniques was computed, namely, the Localization Error (LE), Full Width at Half Maximum (FWHM) of the PSF, and visibility. The definition used for LE was that described in [Pascual-Marqui1994]: the distance between the location of the simulated dipole, and the maximum of the estimated PCD, We define the FWHM directly from the definition by using a nonparametric estimation of the PSF. The measure of visibility was defined as the maximum amplitude of the PSF. The density estimator used was the Naradaya-Watson estimator [Stone 1977].

The measures of quality are summarized in form of nonparametric regression curves in Figures 2-4. Additionally nonparametric standard errors are also included in the figures. Nonparametric regressions were carried out also using the NW estimator.

Figure 2 shows the LE. For LORETA this error decreases monotonically with the eccentricity of the source, from 25 to 14 mm, while for EEG fMRI the LE is constant and almost zero for all points. It is clear that the spatial information offered by the fMRI measurements plays an important role in the improvement of the LE.

Figure 3 shows the FWHM. There is a significant reduction of the FWHM for EEG fMRI with respect to LORETA for all points, indicating an increased spatial sensitivity.

One difficulty of linear inverse methods is the poor visibility of deep sources, which are masked by those sources closer to the sensors. Figure 4. Illustrates that the visibility of LORETA increases exponentially with the eccentricity with a maximum value of approximately 10^-2. In contrast EEG fMRI shows an almost constant behavior, near the optimal value of 1.

In addition to the results of regression equations for all grid points, it is illustrative to compare the actual tomographic solutions for selected grid points. In this case the results of LORETA and EEG fMRI are compared with a fMRI constrained minimum norm type electrophysiological solution (fMRIEEG) previously described in the literature ([Dale 1993;Dale2000]).

In Figure 5-11 shows the results for a source simulated to be in the temporal region. Figure 5 shows the actual source. In Figure 6-8 the temporal evolution of the activation, the modulus of the NMD and the PCD are shown as estimated with EEG fMRI are shown. The modulus of the PCD is reconstructed with a small decrease in the amplitude of the simulated dipole (see fig. 1), which is consistent with the small decay of the visibility with the distance to the sensors. In the case of the hemodynamic variable, the activation is slightly overestimated. Figure 9 shows the map of the PCD obtained with LORETA showing a wide region of activation, and a significant bias in the spatial location of the maximum activity that is due to the smoothing properties of the technique. When using fMRIEEG method (fig. 10) a peak at the actual location of the simulated dipole is observed probably an inheritance of the minimum norm part of the method. The estimation by EEG ↔ fMRI (Fig. 11) is almost indistinguishable from the ideal PCD used for the simulations.

It was also thought convenient to explore the results of a simulation using a distributed source. Therefore a Gaussian function was placed in the occipital region as the simulated source (Fig. 12). The PCD estimated by LORETA is shown in Fig. 13. As expected, LORETA performs very well in this case though there is a small error in the localization error and a reduced visibility. The solution for fMRIEEG method is shown in Fig 14. The reconstruction is quite poor, with maximal activation at the points nearest to the sensors. A sharp frontier for the estimated source is also shown showing that this method is excessively biased for those points where is fMRI activity. For the EEG/MEG fMRI model, there is also a good reconstruction of the spatial distribution of the activity (fig 15), and, unlike LORETA, there is no localization error of the maximum. Nevertheless, there is a significant bias in the estimation of the amplitude of the source.

A quantitative comparison of the simulation results is summarized in tables 1-2.

5. Analysis of a Somatosensory Experiment

The description of the somatosensory experiment analyzed here is described fully in [Hoechstetter 2000].

Anatomical MRI: A T1 anatomical image was image was obtained and used for registration of MEG and fMRI procedures. The scalp and cortex were extracted using software developed at the Montreal Neurological Institute by MacDonald [MacDonald et al. 2000] resulting in triangulated surfaces with normal vectors.

MEG data: Somatosensory evoked fields were recorded from a healthy adult using a Neuromag-122 whole head MEG system with 118 valid sensors used. Prior to recording head position was determined by four coils attached to the scalp. These served to transform the sensor positions and orientations into the MRI space. A 3 sphere model was fitted to the head using the scalp surface.

Brief tactile pressure pulses were delivered by finger clips (BTI) at a constant inter-stimulus interval of 1.03 s. to the tip of the left index finder. Data were sampled at 769 Hz. Event Related Magnetic Fields (ERF) to the SS stimuli were obtained from 739 individual trials, yielding a very good signal to noise ratio. The time course of the ERF was assessed by means of a statistical analogue of the Global Field Power recently described by [F.Carbonell 2001] (consisting of a Hotelling’s T2 statistic over the MEG sensors evaluated at each time point) in order to test the presence of a mean vector. The composite alpha level for the whole time epoch is derived from the theory of random fields [Worsley 1994]. The results fo this computation are shown in Figure 16. As can be seen 3 major peaks are present and these were selected for further processing. In other words the magnetic data to be processed consisted of 3 vectors .

Sources were assumed to be restricted to a 2048-point grid on the cortical surface as defined by the cortical triangulation. The magnetic lead field was obtained using the sensor position and source positions and orientations. Figure 17-a shows the SPM of the LORETA solution for the first vector of the MEG data.

fMRI data: A total of 190 frames where gathered, one every 2 sec. Two conditions were gathered in a blocked design: Baseline vs. Activated state, lasting 20 seconds each. The stimulus was the same as for the MEG recording. Ten slices were selected to pass through the region of interest. The activity in cortex was estimated by an interpolation procedure, sampling restricted to the same points used for the cortical grid described in the previous section. Figure 17-b shows the surface SPM map for this data.

Results of applying the EEG fMRI fusion method are shown in the SPM map in Figure 17-c. As can be seen the localization of activity is much more focal than for LORETA.

6. Final Considerations

The results presented seem to indicate significant advantages for electrophysiological and fMRI image fusion.

The simulation results are very promising. It should be taken into consideration that the cortical PSF based on optical measurements has been estimated in primary visual cortex of the macaque to be 1.5 mm (antero-posterior) by 2.7 mm (medio lateral) ([Grinvald 1994]), while the fMRI resolution in human visual cortex has been estimated to be 3.5 mm ([Engel 1997]). EEG fMRI preserves the good spatial resolution of the imaging techniques mentioned before, with an average FWHM of the order of 1.4 mm. This is much better performance than that of other techniques. What is more important there are considerable improvements in visibility which seems to be the weakness of linear inverse solutions.

The preliminary results presented for the somatosensory data show that the methods proposed may be used for the analysis of actual data, a concern with many current Bayesian image analysis methods.

There is a number of directions in which the present work is being extended. One is concerned with the introduction of temporal correlations in the a priori information. More importantly, the models expressed above are valid for a single experimental condition. This is not the most frequent situation either in fMRI or electrophysiology in which either event related data or block experiments are gathered under some experimental design. The extension of Statistical Parametric Mapping methods for the fusion model is a natural development of the Bayesian approach presented here which will be subject of subsequent presentations.

7. Acknowledgements

We wish to thank Prof. Michael Scherg for providing the experimental data analyzed in Section

References