1. Introduction

An ideal technology for the study of the human brain activity would have high spatial-temporal resolution to reveal the rapid evolution of patterns of event-related cortical activation. At the moment, such a technology does not exist. Thus, there is increasing interest in the combination of neuroimaging techniques having high temporal but low spatial resolution (i.e. electroencephalography, EEG and magnetoencephalography MEG) with those having high spatial but low temporal resolution (positron emission tomography, PET and functional magnetic resonance, fMR). Combined EEG and MEG data increase stability and accuracy in the solution of the linear inverse source estimate of human event-related electromagnetic cortical activity. To date, the linear inverse source estimate from combined EEG and MEG data has been based on (i) preliminary normalization by the covariance of background electromagnetic noise; (ii) realistic MR-constructed head model; (iii) multi-dipole source model; (iv) regularization procedure; and (v) boundary or finite element mathematics (Fuchs et al., 1998). Commonly used regularization schemes assume that between-sensor covariance of background electromagnetic noise is zero, although this assumption rarely holds in the real case. The rationale for the combined use of EEG and MEG data in the linear inverse source estimate is that the MEG activity is poorly contaminated by head volume conduction effects and subcortical fields, while the EEG (but not MEG) is sensitive to the activation of radial cortical sources (Kristeva et al., 1991). The use of fMR as a constraint in the linear inverse source estimate of EEG or MEG data is a promising approach to improve the spatial information content of the source EEG-MEG solutions. In the present study, linear inverse source estimates from combined EEG and MEG data were computed using realistic head and source models, boundary element mathematics, and a regularization scheme taking into account that between-sensor covariance of background electromagnetic noise can be nonzero. Furthermore, fMR data were used as a constraint in the linear inverse source solution of the EEG data. These technologies were used to re-evaluate the controversial hypothesis that bilateral primary sensorimotor areas (M1-S1) are involved in both the preparation and execution of voluntary unilateral finger movements.

2. Methods

2.1 Subjects and task

The present study was carried out on two healthy, right-handed (Edinburgh Inventory) male volunteers. The experiments were undertaken with the understanding and written consent of each participant. General procedures were approved by the local institutional ethics committee. During the experiments, subjects lay on a non-magnetic wooden bed placed in a dimly-lit, sound-damped, and magnetically-shielded room (Vacuumschemelze, Hanau). Head was stabilized by a vacuum cast. Motor task consisted of brisk, internally triggered unilateral right middle finger extensions followed by passive return to the original resting position (inter-movement interval: 2-5 sec). Subjects were asked to avoid blinking, eye movements, and respiration immediately before and during the movement. A brief training was performed to render stable and reproducible the motor performance. During these sessions, surface electromyographic (EMG) activity of extensor digitorum muscle of both sides was monitored (Ag-AgCl cup electrodes; 1-100 Hz bandpass; 400 Hz sampling rate) to control operating muscle response and involuntary mirror movements. Furthermore, surface EMG activity of bilateral axial and proximal muscles was also recorded from the two participants to monitor co-activation of these muscles in concomitance with the finger movement. No notable co-activation of axial and proximal muscles was observed.

2.2 EEG recordings

EEG activity was recorded (0.1-100 Hz bandpass) with 128 electrodes (linked earlobe electric reference) during the movement. Electrode positions and reference landmarks were digitized for subsequent integration between the EEG, MEG, and MR data. Electrooculogram (0.1-100 Hz passband) and electromyogram (EMG, 1-100 Hz passband) from m. extensor digitorum of both sides were also recorded. All data were acquired (400 Hz sampling rate) from 3 sec before to 1 sec after the onset (zerotime) of the EMG response in the operating muscle. About 200 single trials were collected. Artifact-free (eye and/or mirror movements) single trials were averaged with respect to zerotime.

2.3 MEG data acquisition

Movement-related MEG activity was recorded (0.16-250 Hz passband), in separate blocks, from the left and the right hemisphere by a dewar (diameter: 16 cm) including an array of 25 sensors. The sensor array was centered on C3 or C4 site of the 10-20 international system, which are roughly overlying the hand representation of the left and right M1-S1, respectively. This array comprised 9 magnetometers with a 80 mm2 integrated pick-up coil (plus 3 reference channels to be used for noise cancellation) and 16 axial gradiometers (250 mm2 area, 8 cm baseline). The noise spectral density of each sensor channel was 5-7 fT/vHz at 1 Hz. Vertical and horizontal eye movements were controlled by recording electrooculographic activity (0.16-250 Hz bandpass) from a pair of Ag-AgCl cup electrodes placed laterally and medially just above and below the right eye, respectively. The EMG activity of the extensor digitorum muscle of both sides was recorded as in the training session (see above). Right finger movements were performed in blocks lasting about 10 min (10 min inter-block interval). All data were gathered in continuous mode (1000 Hz sampling rate). Positions of the sensor array with respect to subjects anatomical landmarks (nasion, inion, and preauricular points) were detected after 1 or 2 recording blocks, in order to register subtle head movements across the experimental session. Positions of these fiducial landmarks were also digitized off-line for the subsequent integration between MEG and MR data.

2.4 Realistic head and source models

Sixty-four T1-weighted sagittal MR images were acquired (30 msec repetition time, 5 msec echo time, and 3 mm slice thickness without gap). These images were processed with contouring and triangulation algorithms for the construction of subject’s head model. The MR-constructed head model reproduced scalp, skull, and dura mater surfaces with 910 triangles for each surface. Source model was built with the following procedure: (i) the points belonging to the MR images of the cortex were selected with a semiautomatic procedure (thresholding algorithm); (ii) these points were subsampled from 12,000-14,000 to 2,400-4,100; however, the general features of the neocortical envelope were well preserved especially in correspondence of pre- and postcentral gyri and frontal mesial area (mesh refinement), (iii) the subsampled points were triangulated (5,000-7,000 triangles); and (iv) an orthogonal unitary equivalent current dipole was placed in the center of each triangle forming the cortex compartment. Regions of interest were drawn on the reconstructed cortical surface to model supplementary motor area (SMA) and left and right M1-S1.

2.5 FMRI acquisition and analysis

FMR images were acquired from the two volunteers by means of a Siemens Magnetom Vision 1.5 T scanner equipped with 25 mT/m gradients. The motor task consisted of repetitive right middle finger extensions at variable frequency rate (0.5-2 Hz) to make similar the act to that performed during the EEG and MEG recordings. The experimental design included the acquisition of 64 volumes of fMR for the rest-movement-rest cycle. The temporal resolution of each volume acquisition was of 5 sec. Influences of signal intensity drift and head motion on fMR were corrected by normalization and automated image registration algorithms, respectively. The fMR images were co-registered with anatomic images of the subjects cortex, which were obtained by the acquisition of T1-weighted conventional spin-echo-axial-oblique sequence. Percent change of fMR signal intensity due to the movement was computed for each voxel according to the procedure of Kim and co-workers (PC map; Kim et al., 1993 a,b). The difference between the mean rest and movement-related signal intensity was calculated voxel-by-voxel. The rest-related fMR signal intensity was obtained by averaging the pre-movement and recovery fMR. Bonferroni-corrected Students t-test was used to minimize alpha inflation effects due to multiple statistical voxel-by-voxel comparisons (Type I error; p<0.05). Only voxels with a statistically significant PC activation were considered. The fMR values of these voxels were assigned and summed by a minimum distance projection algorithm to the nearest dipoles placed onto the modeled cortical surface.

2.6 Electrical forward solution

Let a head model be constituted by electrically homogeneous and isotropic compartments simulating scalp, skull, and dura mater. The forward solution specifying the potential distribution (V) on these compartments Sk (k = 1,...3) due to a dipole is given by the Fredholm integral equation of the second kind

(1)

with

(2)

where (i) is the potential due to a dipole located in an infinite homogeneous medium; (ii) is the conductivity inside the surface S of the multicompartment head model; (iii) is the conductivity outside the surface Sj; (iv) m is the total number of compartments within the head model; and (v) is the solid angle subtended by the surface element dS located in (point of observation ). A numerical solution of the Fredholm integral equation can be obtained by decomposing the surfaces Sk (k = 1,...3) into triangle panels and by using boundary element techniques (Meijs et al 1993). With the boundary-element techniques, a discrete version of the Fredholm integral equation is given by

v = g + Ω v

(3)

where the elements of matrix Ω , vector v, and vector g are defined as follows: (i) vi is the potential value in the center of mass of the i-th triangle; (ii) gi is the potential value generated by a source in the center of mass of the i-th triangle; and (iii) Ωij is the matrix element proportional to the solid angle subtended by the j-th triangle at the center of mass of the i-th triangle. The numerical solution of the Fredholm integral equation can be improved using the correction proposed by Hämäläinen and Sarvas (1989). The linear system of equation 3 is singular since the potential distribution generated on the scalp compartment by an equivalent dipole is determined up to a constant. This singularity can be removed by using a deflation procedure that yields the potential distribution (V) on the compartment surfaces Sk (k = 1,...3).

2.7 Magnetic forward solution

The magnetic field of a current dipole in a piecewise homogenous conducting medium is given by

(4)

where (i) B₀ is the field of the current dipole in the free space with each surface integral extending to an interface between homogeneously conducting media; (ii) n is a unit vector orthogonal to the surface S_j; and (iii) the summation index j runs over the surfaces. The surface integrals take into account the effects of volume currents. Notice that these currents are equivalent to layers of current dipoles orthogonal to the interfaces, with dipole moment per unit area V(σ⁺- σ^-). The magnetic lead field matrix was calculated by adding, for each dipole, the free space term and the corresponding volume current term. Since a boundary element model was used, the volume currents were modeled by an array of fictitious dipoles, located at the centroids of each triangle of the reconstructed surfaces, with moment V(σ⁺- σ^-) times the triangle area

2.8 Integrated electric and magnetic solution

The forward solution specifying the potential scalp field due to an arbitrary dipole source configuration was computed on the basis of the linear system

Ax = b

(5)

where A is the matrix composed of the electric and magnetic lead fields, b is the measurement array of EEG and MEG data and x is the array of the unknown cortical dipole strengths.

2.9 EEG-MEG regularization

Since the number of dipoles was much higher than the spatial measurements, the linear system of equation 5 has infinite solutions. However, a reliable solution of this system can be obtained with a priori assumption of a minimum energy constraint. The general formulation of the linear inverse problem based on this assumption is

X = argmin (||D(Ax-b)|| 2+ λ2 ||C x||2)

(6)

where D is the sensor equalization matrix obtained by the Cholesky decomposition of the inverted covariance matrix S of the EEG and MEG sensors (DTD = S-1) and C is the matrix that regulates how each EEG or MEG sensor is influenced by dipoles located at different depths into the source model. The covariance matrix S was derived from the EEG and MEG data by maximum likelihood estimation on a set of EEG and MEG data with maximum background electromagnetic noise (i.e., no event-related electromagnetic signal). The matrix C was given by:

C = Wii VT

(7)

where W_ii is a diagonal matrix in which the i-th element is equal to the norm of the i-th column of the normalized lead field matrix DA, and V is the matrix of the principal dipole components obtained by the singular value decomposition of the normalized lead field matrix DA (Fuchs et al 1998). An optimal regularization of the linear system Ax = b was obtained by the L-curve approach (Hansen 1992). The L-curve plots the residual norm versus the solution norm at different lambda (λ) values. Computation of the L-curves and optimal λ correction values were performed by original Hansen’s routines. With equation 7, an optimal regularization value of 1 was always found (Fuchs et al 1998).

2.10 EEG-fMR regularization

For the source space metric computation, we normalized the electric lead field matrix by the column norm, to balance the much more visibility of the superficial than deep cortical sources of the EEG data. Furthermore, the column normalization for each modeled cortical source accounted for the movement-related percentage intensity values (α) of the integrated fMR data (Kim et al., 1993 a,b). The metric for the source space was associated with a diagonal matrix C, whose the i-th term was

Cii = ||Ai||2 * g (αi)-2

(8)

where ||Ai|| is the norm of the i-th column of the lead field matrix A and g(αi) is a function of the statistically significant percentage increase of the fMR signal, assigned to the i-th dipole of the modeled source space. The g(αi) function was expressed by

g(αi) = 1 + K αi; αi >0;

(9)

where the factor K tuned fMR solutions in the source space for the time varying electromagnetic component b. With eq. 9, high K value (i.e., Kαi about 10) produced terms Cii, which were roughly one order of magnitude lower than those used by taking into account the only column normalization. On the other hand, a low K value (i.e., Kαi <<1) resulted in a Cii value roughly proportional to the squared column norm, which completely disregarded the fMR solution (i.e., 0% fMR solution or 100% EEG solution). Optimal value of the K factor was obtained by using the indexes proposed to evaluate the quality of cortical re-constructions in the linear inverse problem (Grave de Peralta et al., 1997).

3. Results

Fig. 1 plots mean MEG and EEG wave forms computed in a subject from the lateral-frontal and medial-parietal areas of both hemispheres in association with the right finger movement. During the movement preparation, the MEG wave forms were represented by a bilateral, contralaterally preponderant slow magnetic shift, starting at about -0.5 sec and culminating at about the zerotime (readiness field-motor field peak; RF-MFp). This complex showed contralateral lateral-frontal positivity (outward current flow)/medial-parietal negativity (inward current flow) and ipsilateral lateral-frontal negativity/medial-parietal positivity. The MEG wave forms presented also a transient magnetic shifts peaking at about +110 msec (motor evoked field 1 peak, MEF1p). With respect to the RF-MFp, the MEF1p was higher in amplitude and had reversed polarity over the contralateral hemisphere. Furthermore, it had low amplitude and non-reversed polarity over the ipsilateral hemisphere (i.e. negativity). On the other hand, the EEG wave forms were characterized by a slow negative shift starting at about -1 sec and peaking close to the movement onset in the frontal, central, and parietal leads (readiness potential peak, RPp). The analogous of the MEF1p was the component peaking at about +110 msec (movement-related response 1 peak, MRR1p).

The mean amplitude maps of Fig. 2 show the topography of the subject's cortically projected electromagnetic activity illustrated in Fig. 1 (recorded data). An intense, dipolar contralateral lateral-frontal and medial-parietal MEG field was observed at RF-MFp and MEF1p, the MEF1p being reversed in polarity compared to the RF-MFp. In addition, RF-MFp presented a low-amplitude MEG field reversed in polarity over the ipsilateral lateral-frontal and medial-parietal areas. The corresponding raw EEG potential distributions were characterized by large and distant negative and positive maxima, preponderant in the side contralateral to the movement. The electric field was tilted of 90° with respect to the magnetic field. Amplitude 3-D maps of Fig. 3 illustrate linear inverse source estimates from EEG, MEG, and combined EEG-MEG data shown in previous figures. These estimates for the movement preparation and execution periods mapped the strength of the dipoles used as a cortical source model. There were circumscribed zone of negativity and positivity in cortical regions roughly corresponding to the M1-S1 of both sides and SMA. The linear inverse source estimate of the combined EEG-MEG data could integrate in a unique solution features of linear inverse source estimate of these data considered separately. Similar results were observed in the other subject.

Figure 1. Wave forms of MEG (left) and EEG (right) activity peaking (p) in the contralateral and ipsilateral hemispheres of a subject during the pre-movement (readiness field-motor field peak, RF-MFp; readiness potential peak RPp) and performance (movement-evoked field 1 peak, MEF1p; movement-related response 1 peak, MRR1p) periods of voluntary unilateral right middle finger extension.

Fig. 4 shows the time evolution (wave forms) of the cortical activation as indicated by linear inverse source estimates in the subject's regions of interest (M1, S1, SMA). FMR data were used as a constraint in the linear inverse source estimation. A map of the linear inverse source estimate at MRR1p (+110 msec) is also shown. This map represents the linear inverse source estimate forwarded over the subject's dura mater model. Modeled wave forms of the cortical activity peaked at +115 msec in the S1 contralateral to the movement. Such an activity peaked few msec later in the SMA and ipsilateral M1 and S1.

Figure 2: Movement-related EEG and MEG activity recorded (128 and 50 channels respectively) about 110 msec after the onset of electromyographic response accompanying the movement. Color maps of the recorded electric potentials and magnetic fields are projected over the subject's MR-constructed cortical model for illustrative purposes. Percent color scale (256 values) is normalized with reference to the maximum amplitude calculated for each map. Maximum negativity (-100%) is coded in violet and maximum positivity (+100%) is coded in red.

Figure 3: Linear inverse source estimates from EEG, MEG, and combined EEG-MEG data shown in figure 2. These estimates for the movement preparation (RF-MFp, RPp) and execution (MEF1p, MRR1p) periods mapped the strength of the dipoles used as a cortical source model. Percent color scale (256 values) in which maximum negativity (-100%) is coded in red and maximum positivity (+100%) is coded in violet.

Figure 4. Wave forms of the linear inverse source estimates of MRR1p computed in the subject's cortical regions of interest (M1, S1, SMA) by using fMR data as a constraint. The illustrated color map represents this estimate forwarded over the subject's MR-constructed dura mater model. Percent color scale (256 values) is used.

4. Discussion

This study presented performances of advanced neuroimaging techniques for the modeling of human movement-related cortical activity from combined EEG, MEG, and fMR data. Linear inverse source estimates of the combined EEG and MEG data were regularized not assuming that between-sensors noise covariance was zero. The regularization scheme extended to the linear inverse source problem for the underestimated cases (more dipoles than sensors) the mathematical approach recently used for the overestimated case (more sensors than dipoles) (Fuchs et al., 1998). Linear inverse source estimates of the combined EEG and MEG data enhanced spatial information content when compared to those obtained from the EEG and MEG data considered separately. In fact, combined EEG and MEG linear inverse source estimates represented in a unique solution spatial features presented in the EEG and MEG linear inverse source estimates considered separately. The results of the present study modeled a maximum activation in the contralateral M1-S1 and in the mesial-frontal cortical area (including SMA). In addition, a minor activation was modeled in the ipsilateral M1-S1 supporting the working hypothesis of a bilateral M1 S1 activation during the preparation and execution of unilateral finger movements. Plausibly, the activation of the ipsilateral M1-S1 during the movement execution is mainly related to the processing of movement-evoked somatosensory information supplied by double crossed and uncrossed pathways (Urbano et al., 1996). Putative double crossed pathway would include dorsal-column lemniscal system and transcallosal M1-S1 connections (Rouiller al., 1994; Wiesendanger et al., 1996). Whereas, putative uncrossed pathways to ipsilateral M1-S1 may comprise spinoreticular, spinomesencephalic, and spinocerebellar pathways (Scheibel, 1984). The ipsilateral M1-S1 activation accompanying the movement execution would subserve transcallosal inhibition of the small uncrossed pyramidal pathway originating from the contralateral M1-S1 (Sadato et al., 1997). This M1-S1 activation might also provide some additional control of the operating hand by means of ipsilateral uncrossed pyramidal pathway (Kristeva et al., 1991). The view of an ipsilateral finger motor representation in M1-S1 is in line with transcallosal M1-S1 connections, uncrossed (about 10%) pyramidal fibers (Wiesendanger et al, 1996), single neuron recordings in monkeys (Tanji et al., 1988), transcranial magnetic stimulation in normal subjects and patients with callosal dysfunction (Meyer et al., 1995; Wassermann et al., 1994). Furthermore, this view is in line with clinical (Haaland et al, 1981, 1989; Jones et al, 1989) and PET (Chollet et al., 1991; Weiller et al., 1996) studies in unilateral hemisphere stroke patients. In contrast, some PET (Colebacht et al., 1991; Roland et al., 1980; Shibasaki et al., 1993) and functional MR (Boecker et al., 1994; Rao et al., 1993, 1995) scans in normal subjects showed no involvement of the ipsilateral M1-S1 in the motor control. A possible explanation for this is that these neuroimaging techniques may have an insufficient sensitivity for a stable detection of slight M1-S1 signal as that related to stereotyped ipsilateral movements. In fact, PET and fMR scans showed ipsilateral M1 and S1 activation only when the motor task consisted of controlled thumb-ulnar finger opposition sequences (Boecker et al., 1994; Kim et al., 1993 a,b; Schroeder et al., 1995).

In conclusion, the proposed neuroimaging technologies used highly sampled EEG and MEG data, realistic MR-constructed subject's head models, fine cortical source model, and fMR data as a physiological constraint of the linear inverse source estimate. These technologies provided evidence in favor of the hypothesis that in right-handed subjects not only the contralateral M1-S1 but also the ipsilateral M1-S1 would subserve the preparation and execution of volitional unilateral one-digit movements. It can be speculated that EEG, MEG, and fMR integration could be largely useful for a study with maximum spatial-temporal resolution of dynamic cortical responses in normal and diseased subjects.

References