A comparison of minimum norm and MUSIC for a combined MEG / EEG sensor array

Many different algorithms for imaging neuronal activity with magnetoencephalography (MEG) or electroencephalography (EEG) have been developed so far. We validate the result of other authors that a combined MEG/EEG sensor array provides smaller source localisation errors than a single MEG or single EEG sensor array for the same total number of sensors. We show that Multiple Signal Classification (MUSIC) provides smaller localisation errors than an unweighted minimum norm method for activity located in the cortical sulcus regions. This is important for many medical applications, e.g. the localisation of the origin of epileptic seizures (focal epilepsy) that can be located very deep in the cortical sulcus.


Introduction
Imaging of neuronal activity can be very helpful in medical applications, e.g.before a resection surgery of brain areas including the origin of epileptic seizures (Pataraia et al., 2002).Functional Magnetic Resonance Imaging (fMRI) measures the oxygen consumption due to synaptic activity and provides a high spatial resolution of the order of a few millimeters but due to the slow haemodynamic response a relatively low temporal resolution of the order of one second (Baillet et al., 2001).In contrast MEG or EEG directly measure the magnetic field or electrostatic scalp potential produced by neuronal currents outside the brain and provide a much higher temporal resolution of a few milliseconds but a lower spatial resolution of a few centimeters.Today in MEG the magnetic fields of the order of a few femtotesla are measured with Superconducting Quantum Interference Devices (SQUIDs) and in EEG the potential is measured with electrodes directly at the scalp.With future uncooled and un-shielded magnetoelectric sensors developed by the research project SFB 855 of the Deutsche Forschungsgemeinschaft maybe also the magnetic field can be measured with sensors directly at the scalp.This could provide on one hand a higher SNR, on the other hand a better alignment of the head to the sensor array for applications in MEG/EEG video monitoring of epileptic seizures.It has been shown before (Sharon et al., 2007;Liu et al., 2002) that a combined SQUID-MEG/EEG sensor array provides smaller source localisation errors than single SQUID-MEG or single EEG sensor arrays for the same total number of sensors.In Sect. 5 we affirm these results with a simulation for a combined MEG/EEG sensor array with MEG magnetometers located directly at the scalp.EEG electrodes and MEG sensors contain different information of the primary currents and for a spherical conductor model even complementary information (Dassios et al., 2007;Dassios, 2008).Due to the non-uniqueness of the solution of the inverse problem (Fokas et al., 2004) lots of different algorithms for the localisation of neuronal activity have been developed in the last decades (Baillet et al., 2001).In Sect.6 we compare the localisation accuracy of an unweighted minimum norm method and the MUSIC (Mosher et al., 1992) beamforming technique for a focal activity located at the cortical gyrus and in the cortical sulcus region of the brain.We show that the localisation error for the cortical sulcus activity is smaller for the MUSIC algorithm than for the unweighted minimum norm method.

Realistic head model
In the simulations of this work we use a realistic head model with the realistic head layers cortex, innerskull, outerskull and scalp extracted from an MRI-Scan of the patients head.The layers are tesselated into a triangular mesh with a vertex Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V. to vertex distance of the order of 1 mm. Figure 1 shows the used layers.The cortical mesh contains 15 010, the innerskull surface 642, the outerskull surface 642 and the scalp surface 1082 vertices.The MRI geometries were provided by Sabine Meunier (Hôpital de la Salpêtrière, Paris) as part of the tutorial dataset distributed with the BrainStorm program (Tadel et al., 2011).

Forward model
The total current density j present in the brain can be separated into a primary current density j p that corresponds to the synaptic activity and a volume current density j V that corresponds to the passive response currents to fulfil the charge conservation ∇ • j = 0: To calculate the magnetic field and the electrostatic potential at the sensor positions we use current dipoles q j fixed at the vertices r j of the cortical mesh to model the primary current density (distributed model): We choose dipole orientations normal to the cortical surface (Baillet et al., 2001, p. 16) and consider n dipoles on the cortical mesh and p sensors.We use a piecewise homogeneous and isotropic conductor model with the 4 conductivity domains brain, skull, scalp and air separated by the 3 layers innerskull, outerskull and scalp as shown in Fig. 1.We use the default conductivity parameters of the BrainStorm program 1 (brain), 0.0125 (skull), 1 (scalp) and 0 (air).The magnetic field and electrostatic potential at the sensor positions is calculated with the Boundary Element Method (BEM) taking into account the contributions of both primary and volume currents (Mosher et al., 1999, p. 246).For the numerical solution of the BEM we use the forward modeling program OpenMEEG (Gramfort et al., 2010;Kybic et al., 2005).The signal of the MEG sensor is proportional to the magnetic flux through the sensor that can be calculated with the superposition principle: For a data time series t 1 ,...,t T the time dependent dipole amplitude matrix Q is defined by and the time dependent data matrix φ B is defined by Inserting Eqs. ( 4) and (5) into Eq.(3) yields a system of linear equations

Forward model for combined MEG/EEG
Using the superposition principle we can derive an analogue system of linear equations for the EEG signal: where V denotes the electrostatic scalp potential and an analogue system of linear equations for a combined MEG/EEG data matrix: where G B denotes the gain matrix for the MEG signal and G V denotes the gain matrix for the EEG signal.

Sensor arrays
Figure 2 shows the combined MEG/EEG sensor array that is used in the simulations.The MEG magnetometers (red) measure just the component of the magnetic field normal to the scalp, the EEG electrodes (green) measure the electrostatic potential at the scalp.Both MEG and EEG sensors are considered as point sensors without noise.

Crosstalk error
The most simple solution of the inverse problem is the unweighted minimum norm solution giving the solution with the smallest norm || Q||.In a forward-inverse-simulation we can insert φ = GQ into Eq.( 9) and get Q = G † GQ.If q j denotes the given amplitude of the j -th dipole, qi = j (G † G) ij q j is the reconstructed amplitude of the i-th dipole.The crosstalk error (Liu et al., 2002) measures how the reconstructed amplitude of the dipole at position i depends on the given dipole amplitudes at other positions j .So the crosstalk error quantifies the mislocalisation of activity.The cortical average of the crosstalk error is given by

Results
Figures 3, 4 and 5 show the crosstalk localisation error plotted on a flat representation of the right cortex for a single EEG, single MEG and a combined MEG/EEG sensor array, respectively.The localisation errors are smallest for the combined MEG/EEG sensor array, especially in the cortical sulcus regions (black arrow).Also the averaged crosstalk C is smallest in the case of the combined MEG/EEG sensor array.

MUSIC
The minimum norm solution of Eq. ( 9) for a distributed model reconstructs the amplitudes of a huge number of dipoles (n = 15 010).In contrast the MUSIC-algorithm tries to represent the measured fields with a smaller number of dipoles r n, namely the number r of active areas in the brain.Assuming that the time series and the field topologies of the r active areas are linearly independent the number of  active areas r can be estimated from the rank of the data matrix r ≈ rank(φ) (Mosher and Leahy, 1998).The positions l i of the r active areas can be found as the r zero points of a cost function u: where g(l i ) are the column vectors of G and P ⊥ is a rank (p − r) orthogonal projector that can be calculated with a Singular Value Decomposition (SVD) of the data matrix φ (Mosher et al., 1992).In praxis the reciprocal of the cost function u is plotted on the cortex and the positions of the r active areas can be found as r sharp peaks in the MUSICmap.In our simulations we consider the number of active areas to be known as one.

Activation in cortical sulcus
Figures 6 and 7 show the results of the simulation in the combined MEG/EEG sensor array with 177 MEG and 177 EEG sensors for the minimum norm and MUSIC, respectively.The given activity is represented by a cortical patch of 8 dipoles activated in the cortical sulcus normally to the cortex with maximum activity at the white dot in the middle.
The localisation error is shown with a white double-arrow and is defined as the distance between the maximum of the given activity and the maximum of the reconstructed activity.The localisation error of the minimum norm solution is 1.8 cm (Fig. 6), the localisation error of the MUSIC-map is 0.4 cm (Fig. 7).The MUSIC-map provides a smaller localisation error than the minimum norm solution.

Activation in cortical gyrus
Figures 8 and 9 show the results of the simulation in the combined MEG/EEG sensor array with 177 MEG and 177 EEG sensors for the minimum norm and MUSIC, respectively.The given activity is located at the cortical gyrus.The localisation error of the minimum norm solution is 0.5 cm (Fig. 8), the localisation error of the MUSIC-map is 0.7 cm (Fig. 9).The difference between the localisation errors of the minimum norm and the MUSIC solution for a cortical gyrus activation is smaller than for a cortical sulcus activation.For a cortical gyrus activation minimum norm and MUSIC provide comparable localisation errors.

Conclusions
In Sect. 5 we showed that a combined MEG/EEG sensor array provides smaller crosstalk localisation errors than a single EEG and a single MEG sensor array for the same total number of sensors.In Sect.6 we compared the unweighted minimum norm solution and the MUSIC solution of the inverse problem for a given one source activation in the cortical sulcus and at the cortical gyrus.For the cortical sulcus activation the localisation error of the minimum norm solution is 1.4 cm bigger than for the MUSIC solution.In the case of the cortical gyrus activation the minimum norm solution and the MUSIC solution are more comparable: the localisation error of the MUSIC solution is 0.2 cm bigger than the localisation error of the minimum norm solution.This motivates the application of MUSIC beamforming techniques for the localisation of deep one source activation (e.g.focal epilepsy) in the brain.

Fig. 1 .
Fig. 1.Realistic head geometry extracted from MRI. Number of dipoles q j on the cortical mesh n = 15 010.

Fig. 2 .
Fig. 2. Combined MEG/EEG sensor array with 177 normal MEG magnetometers and 177 EEG electrodes.In the single MEG simulation the EEG electrodes are substituted by MEG magnetometers, in the single EEG simulation the MEG magnetometers are substituted by EEG electrodes.

Fig. 6 .
Fig. 6.Minimum norm solution (colored) for a given activity in the cortical sulcus and the combined MEG/EEG sensor array from Fig. 2. Left: reconstructed dipole amplitudes plotted in MRI-Scan.Right: activity plotted on a flat representation of the cortex.

Fig. 7 .
Fig. 7. MUSIC-map (colored) for a given activity in the cortical sulcus and the combined MEG/EEG sensor array from Fig. 2. Left: MUSIC-map plotted in MRI-Scan.Right: activity plotted on a flat representation of the cortex.

Fig. 8 .
Fig. 8. Minimum norm solution (colored) for a given activity at the cortical gyrus and the combined MEG/EEG sensor array from Fig. 2. Left: reconstructed dipole amplitudes plotted in MRI-Scan.Right: activity plotted on a flat representation of the cortex.

Fig. 9 .
Fig. 9. MUSIC-map (colored) for a given activity at the cortical gyrus and the combined MEG/EEG sensor array from Fig. 2. Left: MUSICmap plotted in MRI-Scan.Right: activity plotted on a flat representation of the cortex.