International Journal of Bioelectromagnetism Vol. 8, No. 1, pp. 44-55, 2006.

www.ijbem.org

Simple techniques for EEG source imaging

Rolando Grave de Peralta Menendez^a, Patrice Morier^a, Fabienne Picard^b, Theodor Landis^b,
Sara L. Gonzalez Andino^a

^aElectrical Neuroimaging Group, Neurology Department, Geneva University Hospital. Switzerland
^bNeurology Department. Geneva University Hospital. Switzerland.

Correspondence: Rolando Grave ^aElectrical Neurimaging Group, Neurology Dept.,
Geneva University Hospital, 24 Rue Micheli du Crest, 1211 Geneva 14, Switzerland.
E-mail: Rolando.Grave@hcuge.ch, phone +41 22 3728323, fax +41 22 3728333

Abstract. This manuscript describes and illustrates techniques that can be used to simplify and increase the reliability of the solution to the bioelectromagnetic inverse problem. These  steps simplify the head model, the source model and the scalp maps to be localized. The motivation of these techniques is physiological rather than mathematical and they are meant to convert electrical neuroimaging in a reliable, practical tool for brain researchers. The  simplified head model (SYSMAC) considers the individual subject MRI while keeping the computational simplicity of the spherical approximation. The inverse solution can be displayed onto the original anatomical image possibilitating the analysis of patients. We then describe ELECTRA source model based on the fact that the EEG is generated exclusively by irrotational sources. This constraint leads to a better posed inverse problem reducing threefold the unknowns. Finally, a multivariate space time-frequency decomposition (MaSTiF) of the EEG is described to isolate scalp maps of interest in the time-frequency domain. By using a time-frequency decompositition, simultaneous processes occurring at different frequencies are localized separately, decreasing the number of simultaneously active sources. Consequently, the reliability of the inverse solution result is enhanced and the number of ghost sources diminished. These techniques are illustrated in the analysis of sleep spindles, a typical example of brain oscillatory phenomena.

Keywords: Inverse problem, ELECTRA, source localization, time-frequency, neural oscillations, brain rhythms.

1. Introduction

The localization of the generators of electric or magnetic activity recorded at/near the scalp, have been for decades one of the main goals of brain researchers. Not even the recent development of functional imaging techniques (fMRI, PET) has reduced the importance of this problem. This is due not only to the indisputable high temporal resolution of bioelectromagnetic techniques but also because of clinical reasons. The localization of some brain disorders such as the epilepsy is not always reflected in anatomical or functional images. Functional imaging techniques measure physical magnitudes that are weakly related to electrical activity within the brain. Furthermore, almost all normal mental processes are known to occur within 50-500 milliseconds after the presentation of the stimuli. In this fast processing of information many different brain structures are activated serially or in parallel. Some of these areas remain active only a few milliseconds which certainly impede their identification by means of hemodynamic or metabolic image modalities. In addition the existence of a coherent oscillatory rhythmic activity between different areas seems to be the basic principle of transmission of information in the brain. Such oscillatory activity cannot be detected by functional imaging techniques.  This explains why in spite of the mathematical difficulties associated with the solution of the bioelectromagnetic inverse problem; it constitutes still a basic goal for the researchers in the field.

The obstacles to solve the bioelectromagnetic inverse problem have been known for more than 150 years. Helmholtz showed that the bioelectromagnetic inverse problem has no unique solution even for an infinite number of noise-free measurements. Since in practice, the available number of EEG or MEG sensors is technically bounded, it appears that any attempt to deduce the distribution of sources within the brain is doomed to failure. However, it is feasible to obtain estimates for the electrical activity within the brain from electromagnetic signals, provided sufficient a priori information is incorporated in the analysis and the level of detail demanded in the reconstruction is in agreement with the sensitivity of the measurements, i.e., we cannot expect to reconstruct details at the neuronal level from scalp recorded data.

The goal of this paper is to describe and illustrate different techniques that can be used to simplify the solution to the bioelectromagnetic inverse problem while increasing the reliability of the obtained images. The approaches described here comprise steps to simplify the head model, the source model and the scalp maps to be localized. The motivation of these simplified techniques is physiological rather than mathematical and they are meant to convert electrical neuroimaging in a practical and reliable tool for brain researchers.

The manuscript is divided into a theoretical section and the application section. The theoretical section discusses the basis of the techniques, afterward illustrated with the analysis of sleep spindle data. In particular, we include within the theoretical section detailed description of several steps that can be used to simplify the solution, namely:

The second part of the paper illustrates the use of these techniques in the analysis of sleep spindles. Sleep spindles are strong oscillatory phenomena consisting of waxing-and-waning field potentials of 7-14 Hz, grouped in sequences that are short-lasting with a mean duration of around one second. They are considered as sleep maintaining events that block the transfer of sensory information into the cerebral cortex during sleep. Spindles occur with highest density during the stage 2 of sleep and are generated as the result of synaptic interactions of neurons of the reticular thalamic (RE) neurons, thalamocortical relay cells and cortical pyramidal neurons (Steriade 1995).

2. Material and Methods

One fascinating field of application of inverse problem theory is the exploration of the human brain. Almost, if not all the available techniques used to image the brain, e.g., PET, SPECT, MRI, CAT, involve the solution of an inverse problem. Probably the most mathematically challenging of these problems is the one related to the construction of a tomography of the neural activity of the brain based on the electric and magnetic fields measured at/near the scalp surface.

This inverse problem, usually referred as the neuroelectromagnetic inverse problem or the source localization problem or the EEG/MEG source imaging, can be represented by a (first kind) Fredholm linear integral equation. This equation denotes the relationship between the data measured at the external point, d(s), and the superposition of the contribution of the unknown current source density distribution at locations r inside the brain (Fuchs, et al. 1999; Hamalainen 1993; Hamalainen 1992).

The (vector) lead field function L(s,r) contains all the information about the boundary conditions, as well as the media conductivities or permitivities for the electric and magnetic cases, respectively.

Under real conditions, neither the measurements nor the lead field function are known for arbitrary surface/brain locations. However, assuming that the integral equation can be approximated by a discrete sum, Eq. 1 can be represented by an underdetermined system of linear equations:

d = Lj

(2)

Vectors d and j and matrix L represent the discretization of the continuous functions, i.e., d_k = d(s_k), j_m = j(r_m), and L_km = w_kmL(s_k,r_m) and w_km are the quadrature weights for m=1 to Number of solution points and for k=1 to Number of sensors.

All linear solutions of Eq. 2 can be obtained solving a variational problem (Grave de Peralta-Menendez 1998). This yields the inverse matrix G that, when applied to the measured data, produces the estimated current density vector , i.e.:

= Gd

(3)

Substitution of the measured data, as described in Eq.2, into Eq. 3 yields the following fundamental equation for underdetermined linear systems:

= Gd = GLj = Rj

(4)

Here, R=GL denotes the resolution matrix describing the relationship between the estimates and the original magnitudes. The rows of R are called resolution kernels and represent the way current source estimates are distorted by the reconstruction procedure G and the media properties L.

In principle, to solve the inverse problem one has to model two aspects: the sources (generators) and the volume conductor. The volume conductor model is the model used to describe the geometry and electrical parameters (conductivity, magnetic permeability) of the head. Although the most accurate volume conductor model is the realistic one extracted from MRI or CAT of the individual subjects, such images are hardly available in practice for the study of normal subjects. For these reasons and their mathematical simplicity, multicompartimental spherical models are still widely used to model the head.

Because of the shortcomings of spherical head models, efforts have been made over recent years to combine the computational efficiency of spherical head models with more realistic and accurate descriptions of the head shape. Huang (Huang, et al. 1999) proposed a sensor-weighted overlapping-sphere (OS) head model for rapid calculation of more realistic head shapes. The volume currents associated with primary neural activity were used to fit spherical head models for each individual MEG sensor such that the head is more realistically modeled as a set of overlapping spheres, rather than a single sphere. Comparisons of the OS model, with more realistic models and the multishell spherical model showed that the OS model has accuracy similar to the realistic methods but is faster to compute. This model has been extended to EEG by (Ermer, et al. 2001) with improvements in localization accuracy and speed similar to those obtained for MEG. Following the same direction we present in next section a simple alternative to compute realistic head models based on subject magnetic resonance image (MRI).

2.1 Symmetric spherical head model with anatomical constrains (SYSMAC)

Based on the transformation of the MRI into a sphere, SYSMAC shares the simplicity of the spherical model computations and the detailed description of the geometry (i.e. scalp, brain, etc.) of boundary element methods. The main difference with SMAC (Spinelli, et al. 2000), resides in the improvement of the alignment of the MRI image and the electrode configuration (projected on the spherical model). This guarantees that the symmetry plane of both systems perfectly coincides, thereby avoiding that the estimated source activity might change from one hemisphere to the other (Michel, et al. 2004).

The main characteristic of SYSMAC is that anatomical landmarks (Inion, Nasion and Vertex) are selected at the central plane separating both hemispheres (interhemispheric fissure or brain symmetry plane). This information is used as a hard constraint for the computation of the center of the best fitting sphere. Importantly, the fitting is only applied to a subset of points of the smooth part of the MRI, i.e., scalp points with a positive projection on the direction orthogonal to the nasion-inion line and pointing towards the vertex. Once the center is correctly located at the symmetry plane of the MRI, two orthogonal rotations are enough to make the nasion-inion line of both systems parallel. An additional orthogonal rotation will make parallel the pre-auricular lines. Typical errors in the alignment of the theoretical (spherical) model and the MRI are represented in Figure 1. The misalignments results in asymmetric lead field matrices and thus in location biased source reconstructions.

Figure 1. Alignment between the spherical model and MRI. A) Correct: plane of symmetry of both models coincides. B) and C) Incorrect: Symmetry planes of both systems are not aligned due to incorrect position of the center (C) or lack of parallelism (B).

In summary the main advantages of SYSMAC can be itemized as follows:

2.2 Multivariate Space Time Frequency (MaSTiF) decomposition of the EEG

The starting point for the analysis is the recorded EEG or MEG data D (with columns d as in Eq.2 and 3). This data matrix consists of N_s rows determined by the sensors and N_t  columns associated to the temporal evolution of these signals. Since the sensors are associated with locations over/near the scalp, the columns of this matrix can be interpreted as spatial patterns while its rows as temporal patterns. The following decomposition emphasizes this interpretation in terms of the superposition of multiple components:

where K defines the number of components yielding the noise free data. The i-th column of S and the i-th row of T are the spatial and temporal patterns associated to the i-th component.

The decomposition in Eq 5 is not unique and some other constraints are needed. A very general framework comprising several standard and emerging techniques is the Blind Signal Processing approach (Cichocki 2004), that includes, among others, the Principal Components Analysis (PCA) method and the Independent Component Analysis (ICA) method. Whereas PCA searches for orthogonal patterns in both the space and the time domain, ICA searches for spatial patterns with statistically independent time courses (temporal patterns). Although the assumptions about the orthogonality or the independence of the temporal patterns are debatable, these methods (PCA and ICA) have two clear limitations namely, 1) they have no information about the spectral components of the signal, i.e., insensitive to non stationary data, and 2) they provide no information about the reliability of the estimated spatial patterns.

An alternative decomposition of the data that takes into account the spectral components of the signal and the multivariate nature of the measurements was proposed in (Gonzalez Andino, et al. 2000). This multivariate space-time frequency decomposition produces for each time t and frequency value w, a (space) data map d(t,w) representing the (approximate) contribution of that time and that frequency to the observed (measured) map d(t). As for the PCA and ICA, these maps can be submitted to source localization procedures to obtain the generators associated to it.

MaSTiF decomposition method is based on the time frequency representation of all sensors considered together. Lets denote by C(t,w) the column vector composed by the complex time frequency coefficient of all sensors at (t,w). If all sensors have  approximately the same phase at (t,w), then C(t,w) is approximately equal to the real vector d(t,w) times a complex number. If this is true, the ratio between the lowest eigenvalue and the highest eigenvalue of the 2D covariance matrix computed from the complex numbers in vector C(t,w) should be small (close zero). Otherwise it can increase up to one. This measure called simplicity test can be plotted for all time frequency values to produce an image m(t,w) indicating the time frequency points where the approximation can be granted. This is extremely useful to compute the brain sources associated to spikes or to study epilepsy propagation during the first milliseconds of seizure onset. For an example of its application to seizure propagation analysis see (Gonzalez Andino, et al. 2000).

2.3 The  Irrotational source model: ELECTRA

Several (theoretical) source models have been used to solve Eq.1 and 2 and thus to describe the sources of the electromagnetic activity of the brain, e.g., dipoles, monopoles, current density vector. However, none of these theoretical source models actually exists within the brain nor is any physically measurable. Instead, real measurements are the result of quantifiable potentials at different “measurable” levels.  At the microscopic (neuron) level, this is the membrane potential. At the macroscopic (region) level, this is the local field potential (LFP). Through volume conduction, these potentials arrive at the scalp where they are measured as the Electroencephalogram (EEG). It is then natural to question whether potentials inside the brain can be related to and thus computed from potentials measured at the scalp.

A positive answer to this question can be given if we notice that macroscopic primary sources, i.e. the generators of the EEG, are dominated by microscopic secondary (volume) currents or in Plonsey words (Plonsey 1982) that “the fields measured do not even arise from J [the current source density vector field] but rather from secondary sources only. These secondary sources, in turn, depend on both the electrical field and the interfaces, and hence are related to divergence of J and the geometry”. We would note that this kind of source corresponds to a potential distribution inside the brain.

A definitive theoretical argument can be obtained if we notice that, according to the Helmholtz theorem, the current density vector field can be written as the sum of a solenoidal vector field plus an irrotational vector field plus the gradient of a harmonic function. Based on Green identities, it follows that only the irrotational current contributes to the measured potentials (EEG). In mathematical parlance, it means that the EEG generators fulfill:

× J = 0 <=> J = φ

(6)

where φ is a potential field within the brain. Assuming piece-wise constant conductivities it is very easy to show (using Poisson equation) that φ has the same sources and sinks as the EEG potential. This argument explains why we use the term Local Field Potential to denote these non-invasive estimates. Still these estimates are far from having the spatial resolution of LFP recorded with micro-electrodes in animals. Our non-invasive estimates own, in the best case, a spatial resolution comparable to that of intracranial recordings in epileptic patients.

Moreover, plotting the modulus of the estimated primary current obtained by Eq.1 or Eg.2, which we would note has thus far been the common procedure used to depict inverse solutions results, does not reflect the actual generators. Instead, the actual generators are determined by the sources and the sinks obtained from the Laplacian of potential field φ or the divergence of the primary current density vector (•J ).

The irrotational source model (ELECTRA) corresponds to the solution of one of the following problems (Grave de Peralta et al., 2000; Grave de Peralta et al., 2004): 1) The estimation of an irrotational current density vector J = φ with the standard vector lead field. 2) The estimation of a scalar field, the current source density (CSD), •J = I with a scalar lead field. 3) The estimation of a scalar field, the potential distribution φ inside the brain, with a transformation of the standard vector lead field.

In summary, the main advantages of the irrotational source model of ELECTRA are:

A final theoretical point to discuss is the possibility to use the irrotational source model with magnetic measurements (MEG). While the general electromagnetic formulation cannot exclude the existence of rotational sources (rotor or curl different from zero), the results of Plonsey about the sources of the bioelectric and biomagnetic fields are rather conclusive: “Even if the divergence and curl of the primary source were independent (and hence were both needed to define the primary source), because the secondary sources all arise from the divergence of the primary source the magnetic field reflects the same source component as the electric field”. In our opinion, this argument speaks in favor of applying the ELECTRA source model also to magnetic measurements.

3. A practical example: Studying sleep spindles.

We selected sleep spindles as the example to illustrate this technique because they are one of the oldest oscillatory phenomena known but little information is available about its functional role. The study of spindles, a pure oscillatory phenomenon, requires the use of non-invasive high temporal resolution techniques neuroimaging techniques as the one proposed here. Techniques based on the hemodynamic responses are of little help in their study.

Some recent studies link, sleep spindles with episodic memory consolidation (Stickgold 2004). Regarding declarative memory consolidation, sleep may provide a state during which newly acquired memory contents, which are temporarily stored in the hippocampus, can be easily transferred to the neocortex for integration into long-term memories (Buzsaki 1998). Siapas and Wilson (Siapas and Wilson 1998) reported that hipocampal ripples, known to play a role in memory formation, occur in a temporal correlation to sleep spindles in single cell derivations.

3.1 Data Collection and Preprocessing

Nocturnal EEG was collected from 64 electrodes in an epileptic patient using a sampling rate of 200 Hz. Electrodes were positioned according to the extended 10/20 system as indicated in Figure 2. Preprocessing consisted in a nearest neighbour based interpolation of electrode P5 and transformation of the data to the average reference. Sleep spindles were automatically selected from the stage 2 of sleep using a semiautomatic procedure described below.

To select sleep spindles we used the concept of spectral envelope. A set of two frontal, two central and two parietal EEG channels were a priori selected. The EEG at these channels was narrowly band filtered within two frequency ranges (11-13 Hz and 13-15 Hz). The band-pass filter was based on a simple direct and inverse Fourier transform with pruning of undesired frequencies. The envelope was computed for the filtered signals in both ranges by transforming the filtered data to z-scores followed by a Hilbert transform. The Hilbert transform returns the analytic signal whose magnitude is the complex envelope of the original signal. Spindles were marked if their envelope showed the characteristic waxing and waning properties characteristic of sleep spindles in frontal, central or parietal derivations. Typical examples of the signals and their envelopes are shown in Figure 3.

Figure 2. Flat view of electrode positions for the recorded nocturnal EEG.

Figure 3. EEG signal at a frontal derivation (blue) filtered between 11 and 13 Hz and the spectral envelope.

3.2 (MaSTiF) decomposition of the sleep spindles.

EEG segments containing two or more consecutive spindles were selected for analysis. The goal was to study the whole spatio-temporal distribution of the spindle generators not only at the time of spindles but also between spindles. The rational to analyze the inter spindles interval derives from the experimental observation that in rats cortical spindles tend to follow hipocampal ripples (Siapas and Wilson 1998).  It has been suggested that the hippocampal circuitry may be biasing the initiation of the spindle–ripple episodes either through direct anatomical projections (Jay and Witter, 1991) or indirectly through neural systems such as the entorhinal cortex or the basal ganglia. We were therefore interested in observing if temporal lobe structures were activated during the inter-spindles interval.

Figure 4 shows some original EEG traces for a segment where three spindles were identified. The vertical markers indicate the timing of the spindles as detected from their spectral envelope. Figure 5 (lower panel) shows the modulus of the MaSTiF decomposition, d(t,w), for the EEG segment shown in Figure 4 and its associated simplicity test map (Figure 5, upper panel). The simplicity map indicates the time frequency points where the map seems to be produced by a simple generator configuration, i.e., all sensors have approximately the same phase for the given time-frequency pair.

In this particular example the MaSTiF approximation was computed using the Fourier transform based on overlapping windows of duration 2 seconds. This leads to a frequency resolution of 0.5 Hz. The overlapping of consecutive windows was 1.9 seconds. Thus, the analysis reflects an effective temporal resolution of 10 ms and a frequency resolution of 0.5 Hz. While the frequency resolution of this method is similar to the one achievable by other spectral methods, increasing the overlap can increase the temporal resolution up to the original resolution of the data. Although such excellent temporal resolution can be achieved by time-frequency (TF) based analysis procedures, we would note that TF’s alone do not directly lead to interpretable potential maps as done by MaSTIF.

Figure 4. Original EEG traces with markers (green vertical lines) indicating the spindles identified from the spectral envelope.

Figure 5. MaSTiF approximation for the spindle data shown in Figure 4. The lower panel depicts the global time frequency energy. Darker colors indicate strong energy concentrations around the specified time-frequency region. The upper panel shows the results of the simplicity test. Darker colors indicate simple maps likely to be generated by simple configurations of the generators.

 The results of Figure 5 indicate the presence of the several time-frequency energy spots (see lower panel). In particular, the three spindles shown in Figure 4 are clearly delineated as strong time frequency energy peaks at 3.6, 6.8 and 11 secs. Interestingly, the plot indicates the simultaneity of 12 and 14 Hz spindles in nearly all cases. This is suggestive of the coexistence of different neural processes within a single spindle. The evidence that these are in fact different neural processes that coexist in time is given by the differences in topography of the simultaneous spindles at 12 and 14 Hz shown in Figure 6. This plot presents an expanded view of the middle spindle occurring between 6 and 8 secs. This spindle is selected because the simplicity plot suggests the existence of simple maps for both frequencies, i.e., 12 and 14 Hz. Importantly, while the existence of different scalp maps is a certain evidence of difference in neural generators, the opposite does not hold. Similar scalp maps can be produced by different configuration of generators due to the existence of silent sources, i.e., sources that produce no scalp EEG (e.g. closed fields). In fact this is the origin of the non-uniqueness of the neuroelectromagnetic inverse problem.

Figure 6. Topographic distribution of sleep spindles for the 12 Hz (left) and 14 Hz (right) time frequency spot. Top plots show the spatial distribution of waveshapes obtained after MaSTiF for the 12 Hz (left) and 14 Hz (right). The sequence of maps is shown for each range to demonstrate the stability of the map sequence during the spindle. Notice that except for polarity inversions, the same scalp spatial pattern remains during the whole spindle cycle. Slower spindles show frontal distribution.

To localize the generators of the sleep spindles we used the simplified SYSMAC head model and the variant of ELECTRA source model based on the estimation of intracranial potentials. The SYSMAC model was constructed on the basis of a standard head model since the subject individual MRI was not available. A total of 4024 knots, homogeneously distributed over the gray matter, composed the solution space for an effective resolution of 6 millimeters. Conductivities for the lead field construction were assumed as in Stok.

The localization results for the spindles at 12 and 14 Hz shown in Figure 6 are displayed in Figure 7 and 8 respectively. While the localization results cover the same period than in Figure 6, we have reduced the number of images shown over the interval to facilitate visualization of the results. The localization results are very stable over the spindle period as expectable from the stability observed for the scalp maps. The slower spindles at 12 Hz have maximum (minimum) that localize to the middle frontal lobe (Brodmann area 9). The frontal localization was consistently observed over several explored spindles although some spindles showed maxima at the inferior frontal gyrus rather than the middle frontal gyrus shown here. Hemispheric lateralization was not clear though most spindles start at the left or right hemisphere and progress towards a more bilateral involvement around the spindle peak.

Figure 7. ELECTRA localization for the slower spindles at 12 Hz. The sequence of maps represents the temporal evolution of localization over the spindle cycle. Maxima (positive potential) are represented in red an marked with a red cross and minima in blue with a blue cross indicating the exact position. Note that the maxima and minima remain stable over the window except for the poalrity inversions.

The maxima for faster spindles at 14 Hz were systematically localized at parietal or parieto-occipital areas as shown in Figure 8. Curiously, the frontal spindles progressively involved anterior cingulate areas toward the spindle maximum while parietal spindles recruited motor cortical areas during their development. Apparently, the increase in amplitude of the spindles is not only due to the increase in activity at frontal or parietal areas but also to the progressive recruitment of additional structures anatomically connected to the site where the spindle initially developed, i.e., frontal cortex to cingulate and parietal cortex to motor areas.

Figure 8. ELECTRA localization for the faster spindles at 14 Hz. Note that the maxima and minima remain stable at parietal locations over the spindle cycle except for the polarity inversions. Note the progressive recruitment of sensory-motor cortex over the cycle.

One interesting finding of this analysis was that we often observe strong activity localized at brain voxels located at near the putative position of the thalamus in our head model. Such deep activity corresponded even with the maximum or the minimum at some punctual times that preceded the spindles. Such strong thalamic activity was never observed at the spindle during the times of maximal spindle waning or waxing. Later observation reveals that deep brain activity can be observed at the scalp and thus localized with inverse solutions when the cortical areas are relatively silent and the deep activity is strong. This is likely to be the case for sleep spindles for which thalamus is known to be a basic generator. One example of such activation is presented in Figure 9.

Figure 9. Voxels at the putative localization of the thalamus can appear as the maxima or the minima of the activity recovered by ELECTRA inverse solution. Here the maximum is indicated by the red cross. Thalamic activity was often observed preceding the start of the spindle, i.e. when cortical areas were practically silent.

Figure 10. Stable temporal lobe activity is observed during the inter-spindles intervals.

We were also interested in analyzing the stability and generators of the activity during the inter-spindles intervals. Experimental evidence (Sirota, et al. 2003) indicates that sleep spindles promote hipoccampas ripples after a certain temporal delay. It is therefore interesting to explore the areas activated after spindles. We did in fact found stable activity after all spindles explored. The stability occurred preferentially after the slower 12 Hz spindles that showed maxima at temporal brain areas. In general the frontal maximum observed during the spindle shifted toward temporal areas of the same hemisphere where the frontal maximum was found. This effect is illustrated in Figure 10. The sequence of maps show the end of the spindle with the right frontal maximum and minimum followed by the increase in activation at the right middle temporal lobe that soon becomes the maximum. Activity is progressively shifted towards the right superior temporal gyrus. The axial slide shown in Figure 11, at the level of the blue line shown in the rightmost inset, illustrates that the temporal activity observed after spindles is not restricted to the cortex but also involves deep temporal lobe structures such as the hippocampal/amygdala formation.

Figure 11. Simultaneous cortical and subcortical activity (amygdala/hippocampus formation) is observed during the interspindles interval.

3. Discussion and Conclusions

We have described in this manuscript some techniques that help to simplify the solution of the bioelectromgnetic inverse problem. The use of the proposed head model, the SYSMAC, permits a realistic selection of individual anatomical details for lead field design and consequently realistic inverse solution displays whereas keeping the computational simplicity of spherical head models.  The SYSMAC model constitutes a refinement of the originally proposed SMAC that guarantees a correct hemispheric alignment. This is a non-trivial detail for multiple applications of the inverse problem such as epileptic focus localization or neurocognitive studies aimed to assess interhemipheric functional differences in tasks such as language. The use of a solution space that conforms to the individual subject brain is essential for functional localization given the large interindividual variability observed in cortical anatomy and the variability in functions. The images presented in this study are an example of the level of realism that is attainable with this easy to compute head model.

The simplification in the source model introduced by the variant of ELECTRA used here, based on the estimation of intracranial potentials, is both conceptual and computational. This variant reduces threefold the number of unknowns to be estimated with the same amount of data when compared with source models based on the whole current density vector. This reduction complies with the idea that the part of the current density vector that is silent to the measurements should not be searched for. Silent sources cannot be recovered unless complete a priori information about them is provided by alternate simultaneous measurement modalities. Furthermore, ELECTRA estimates the temporal course of intracranial potentials amenable to direct comparisons with actual intracortical recordings in patients. This is not the case for current density estimates or dipolar moment estimates that are often incorrectly compared with direct recordings. Because the current density is the first spatial derivative of the potentials, a direct comparison between current estimates and intracranial recordings makes as little sense as comparing apples to oranges.

The MaSTiF decomposition tackles a very important problem in modern Neuroscience by providing a mathematically founded tool for the analysis of oscillatory brain phenomena with a temporal resolution that can be identical to the one in the original data. Importantly, the assumption of stationarity of the time series, hardly plausible for most electrophysiological processes is avoided in this approach based on linear time-frequency decomposition.  Notice that the principle behind the MaSTiF decomposition differs from the pure mathematical rationale behind decompositions as ICA or PCA. MaSTiF, is compatible with the experimentally supported idea (Andersen and Buneo 2002) that the same neural population codes within different frequency bands information about task parameters and intentions. These concurrent oscillatory phenomena difficult the analysis and interpretation of the recorded data since simultaneously occurring oscillations cannot be separated in the temporal domain.

The MaSTiF decomposition is not a data reduction technique as ICA or PCA. It rather increases the amount of data by splitting concurrent processes. The plots described in the paper, i.e., the time frequency global energy plot and the simplicity plot compensate such data splitting. They serve to provide a fast view of long spontaneous multivariate data resuming in a single image the whole information contained in the spatial, spectral and temporal domains. The plots help to immediately track oscillatory events of interest, to isolate them and to decide based on the simplicity, the scalp maps where reliable localization is to be expected. Moreover, it is precisely the data splitting that allows for better and more reliable localization results. First, the MaSTiF method filters out noise that tends to be confined to frequency bands different from the ones reflecting actual electrophysiological processes. Second, an instantaneous scalp map is approximately, the superposition of all scalp maps generated by MaSTiF for all frequencies at the corresponding time. This trivially implies that a single time-frequency map generated by MaSTiF is simpler in terms of the amount of active sources than the corresponding instantaneous map. The existence of simultaneously active sources complicates the task of linear inverse solutions leading to an increase of ghost and lost sources (Grave de Peralta-Menendez 1998). Thus, as any method that minimizes the likelihood of multiple active sources results, the localization results obtained for the MaSTiF transformed data are more credible than the localization results obtained for the original raw EEG.

The analysis of sleep spindles described here serves to illustrate the techniques. Interpreting the results in neurophysiological terms would be too precipitated since a single subject was considered in the analysis. Still, our results are concordant with previous studies in terms of the separation (Durka, et al. 2005) of fast and slow spindles and their differential cerebral localization. In agreement with Anderer (Anderer, et al. 2001) we found that slow spindles had maxima at the frontal lobe while faster spindles at 14 Hz preferentially activated parietal areas. We did also observe a progressive recruitment during the spindle cycle of Anterior Cingulate Areas by frontal spindles and of sensory-motor areas by parietal spindles. This result requires to be confirmed over a larger amount of subjects. Interestingly, during the interspindles interval we observe systematic activation of temporal lobe structures at both the cortical surface and subcortical structures systematically involved in memory. Such relationship, if proved to be consistent over a larger number of subjects could constitute an additional evidence to support the proposed role of sleep spindles in memory consolidation.

Our results show that deep brain activity, in particular thalamic activity, can be recovered by inverse solution methods. This might surprise those who believe that scalp measured fields (EEG and MEG) are not sensitive to deep sources. In fact, there are some experimental studies showing that thalamic activity can be detected by external MEG recordings (Ikeda, et al. 2002) even though the MEG is no more sensitive  than the EEG to deep sources (Malmivuo, et al. 1997). Sources are silent, i.e., not measurable by the scalp EEG or MEG, not because they are deep or superficial but because they have specific field configurations (e.g., closed fields) or physical properties (e.g., irrotational sources produce no EEG). This means that even if spatially smeared, deep sources are reflected in the measurements and can be thus localized by inverse solutions.

Acknowledgements

The Swiss National Science Foundation (grants 3152A0-100745, the IM2.MI on Brain Machines Interfaces) and the European IST Programme FET Project FP6-003758 provided financial support. This paper only reflects the authors' views and funding agencies are not liable for any use that may be made of the information contained herein. We thank Dennis Brunet as main author of Cartool software used to build the images.

References

Anderer P, Klosch G, Gruber G, Trenker E, Pascual-Marqui RD, Zeitlhofer J, Barbanoj MJ, Rappelsberger P, Saletu B. (2001): Low-resolution brain electromagnetic tomography revealed simultaneously active frontal and parietal sleep spindle sources in the human cortex. Neuroscience 103(3):581-92.

Andersen RA, Buneo CA. (2002): Intentional maps in posterior parietal cortex. Annu. Rev. Neurosci. 25:189-220.

Buzsaki G. (1998): Memory consolidation during sleep: a neurophysiological perspective. J. Sleep. Res. 7 Suppl 1:17-23.

Cichocki A. (2004): Blind signal processing methods for analyzing multichannel brain signals. International Journal of Bioelectromagtism 6(1).

Durka PJ, Matysiak A, Montes EM, Sosa PV, Blinowska KJ. (2005): Multichannel matching pursuit and EEG inverse solutions. J. Neurosci. Methods 148(1):49-59.

Ermer JJ, Mosher JC, Baillet S, Leah RM. (2001): Rapidly recomputable EEG forward models for realistic head shapes. Phys Med Biol 46(4):1265-81.

Fuchs M, Wagner M, Köhler T, Wischmann H. (1999): Linear and non-linear current density reconstructions. J. Clin. Neurophysiol. 16:267-295.

Gonzalez Andino SL, Michel CM, Lantz G, Grave de Peralta Menendez R. (2000): Non stationary source approximation: an alternative to improve localization procedures. NeuroImage 11(5, Supplement 1):S502.

Grave de Peralta Menendez R, Gonzalez Andino SL, Morand S, Michel CM, Landis T. (2000): Imaging the electrical activity of the brain: ELECTRA. Hum. Brain Mapp. 9(1):1-12.

Grave de Peralta Menendez R, Murray MM, Michel CM, Martuzzi R, Gonzalez Andino SL. (2004): Electrical neuroimaging based on biophysical constraints. NeuroImage 21(2):527-539.

Grave de Peralta-Menendez R, Gonzalez-Andino, Sara L. (1998): A critical analysis of linear inverse solutions to the neuroelectromagnetic inverse problem. IEEE Transactions on Biomedical Engineering 45(4):440-48.

Hämäläinen M, ; Hari, R; Ilmoniemi RJ; Knuutila J; Lounasmaa OV. (1993): Magnetoencephalography-Theory, Instrumenstation and Applications to noninvasive studies of the working human brain. Reviews of Modern Physics 65(2):413-497.

Hämäläinen MS. (1992): Magnetoencephalography: a tool for functional brain imaging. Brain Topogr. 5(2):95-102.

Huang MX, Mosher JC, Leahy RM. (1999): A sensor-weighted overlapping-sphere head model and exhaustive head model comparison for MEG. Physics in Medicine and Biology (2):423-440.

Ikeda H, Leyba L, Bartolo A, Wang Y, Okada YC. (2002): Synchronized Spikes of Thalamocortical Axonal Terminals and Cortical Neurons Are Detectable Outside the Pig Brain With MEG. J. Neurophysiol. 87(1):626-630.

Malmivuo J, Suihko V, Eskola H. (1997): Sensitivity distributions of EEG and MEG measurements. IEEE Trans. Biomed. Eng. 44(3):196-208.

Michel CM, Murray MM, Lantz G, Gonzalez S, Spinelli L, Grave de Peralta R. (2004): EEG source imaging. Clin Neurophysiol 115(10):2195-222.

Plonsey R. (1982): The nature of sources of bioelectric and biomagnetic fields. Biophys. J. 39(3):309-12.

Siapas AG, Wilson MA. (1998): Coordinated interactions between hippocampal ripples and cortical spindles during slow-wave sleep. Neuron 21(5):1123-8.

Sirota A, Csicsvari J, Buhl D, Buzsaki G. (2003): Communication between neocortex and hippocampus during sleep in rodents. PNAS 100(4):2065-2069.

Spinelli L, Andino SG, Lantz G, Seeck M, Michel CM. (2000): Electromagnetic inverse solutions in anatomically constrained spherical head models. Brain. Topogr. 13(2):115-25.

Steriade M. (1995): Thalamic origin of sleep spindles: Morison and Bassett (1945). J. Neurophysiol. 73(3):921-2.

Stickgold R. (2004): Dissecting sleep-dependent learning and memory consolidation. Comment on Schabus M et al. Sleep spindles and their significance for declarative memory consolidation. Sleep 2004;27(8):1479-85. Sleep 27(8):1443-5.

Stok CJ. (1987): The influence of model parameters on EEG/MEG single dipole source estimation. IEEE Trans. Biomed. Eng. 34(4):289-96.