Mutual information spectrum for selection of event-related spatial components. Application to eloquent motor cortex mapping

Spatial component analysis is often used to explore multidimensional time series data whose sources cannot be measured directly. Several methods may be used to decompose the data into a set of spatial components with temporal loadings. Component selection is of crucial importance, and should be supported by objective criteria. In some applications, the use of a well defined component selection criterion may provide for automation of the analysis. In this paper we describe a novel approach for ranking of spatial components calculated from the EEG or MEG data recorded within evoked response paradigm. Our method is called Mutual Information (MI) Spectrum and is based on gauging the amount of MI of spatial component temporal loadings with a synthetically created reference signal. We also describe the appropriate randomization based statistical assessment scheme that can be used for selection of components with statistically significant amount of MI. Using simulated data with realistic trial to trial variations and SNR corresponding to the real recordings we demonstrate the superior performance characteristics of the described MI based measure as compared to a more conventionally used power driven gauge. We also demonstrate the application of the MI Spectrum for the selection of task-related independent components from real MEG data. We show that the MI spectrum allows to identify task-related components reliably in a consistent fashion, yielding stable results even from a small number of trials. We conclude that the proposed method fits naturally the information driven nature of ICA and can be used for routine and automatic ranking of independent components calculated from the functional neuroimaging data collected within event-related paradigms.

1. Introduction

Spatial decomposition is one of the key techniques applied to exploratory analysis of multichannel data in general, and to spatial-temporal electro- and magnetoencephalographic (EMEG) signals in particular. The most commonly used methods to obtain both spatial components and the corresponding temporal loadings are independent component analysis (ICA) (Comon, 1994), principal component analysis (PCA) (Golub and Van Loan, 1996) and factor analysis (FA) (Child, 2006).

The most frequently used approach for analysis of stimulus-locked averaged EMEG data is PCA, which can be performed using the singular value decomposition (SVD) of the stimulus-locked averaged data matrix (Lagerlund et al., 1997). This analysis is followed by thresholding the singular values (SV) spectrum to identify the subspace capturing the largest amount of data variance for a given approximation rank (Vandewalle, 1988). This technique is inherently power driven. Its application to the identification of the repetitive task-related signal subspace from the averaged ERP/F data relies on the assumption that the individual evoked responses are sufficiently well phase-locked to the stimulus. In that case, the stimulus-locked summation results in an enhanced relative power of the phase-locked component (Misulis, 1994).

SVD is the optimal method for signal subspace detection measured by subspace correlation for a given approximation rank (Vandewalle, 1988). However, the actual value of signal subspace rank, R, is, in general, unknown. Finding an estimate of R is not a trivial task. It is often done by visual inspection of the SV spectrum. The method is based on identifying the target index, R, of a singular component just preceding a sharp drop in power, followed by a slow decaying plateau in the SV spectrum. However, a large disparity of activation amplitudes, spatial proximity of the neuronal sources and powerful noise sources may result in the absence of a clear cut division between task-related and noise components. In addition, in realistic conditions, the recordings are often contaminated by spatially colored brain activity and/or spatially coherent artifacts. Under these circumstances, component selection based on the SV spectrum may be misleading. As a motivating example, consider the top panel of Figure 1 that shows the SV spectrum calculated for the averaged data obtained from the simulated MEG timeseries containing the contribution of two non-synchronous dipolar sources. Although the subspace spanned by the first $\widehat{R}$ = 2 singular topographies almost exactly matches the true subspace (subcorr([a₁ a₂], [u₁ u₂]) = [1, 0.987]) the spectrum of SVs fails to provide evidence that the second component contains task-related signal.

FIGURE 1

Figure 1. In some practical cases, the first singular topographies remain reasonably good estimators of signal subspace but inspection of the SV spectrum fails to reveal this as illustrated by a model example here. While the subspace spanned by the first two singular topographies and the actual simulated subspace practically coincide the SV spectrum (top panel) fails to reveal the fact that the second singular component also belongs to the signal subspace. On the contrary the MI spectrum (bottom panel) computed using raw data projected onto the left singular vectors demonstrates a very clearly cut separation of the task-related and task unrelated subspaces. This figure also introduces the measure of task-related subspace identifiability used in the paper. Since the correct signal subspace rank value is R = 2, we use discriminating indicator $q=log\left(\frac{S(2)}{S(3)}\right)-log\left(\frac{S(3)}{S(4)}\right)$ that formalizes the strategy employed by the human observers and estimates the amount of drop between the second and the third components referenced to the ratio of the two largest components of the noise range spectrum (with indices 3 and 4) immediately following the two signal components (with indices 1 and 2).

ICA is one of the most widely used approaches to blind source separation, popular for exploratory analysis of multidimensional data. In the analysis of EMEG data from evoked response experiments, this decomposition may be used both to isolate task-related components (Makeig et al., 1996; Vigario et al., 2000) and to remove artifacts (Jung et al., 2000). It can also be used for estimation of source timeseries when proper forward modeling is unavailable and for the estimation of the signal subspace in cases when the experimental paradigm can not guarantee sufficiently accurate stimulus locking, e.g., in voluntary movement paradigm (Ossadtchi et al., 2000; Delorme, 2010).

The application of ICA to routine analysis of EMEG datasets is limited by the absence of standard approaches for ordering the independent components (Hyvärinen et al., 2001). In the most typical scenario a human observer visually identifies the desired components by exploring their timecourses and topographies. Since the raw EMEG data are often very noisy, it can be difficult to determine which components should be selected. Additionally, such manual selection is often daunting, and selection based on power (Delorme, 2010) or stim-locked averaged power (AP) (Hyvärinen et al., 2001) does not ensure that the components are event-related for reasons essentially similar to those just described for the SV spectrum. It should be also noted that such power-driven ordering methods may be inappropriate, since ICA as an information (rather than power) driven technique. An alternative method, using the correlation metrics between each estimated component and the event trigger, is critically dependent on signal shapes, and is therefore highly unreliable. For completeness we will mention that for some methods of blind source separation, such as AMUSE (Tong et al., 1991), components may have an intrinsic order but such an ordering is not very useful in the context of analysis of EMEG data from ERP studies. These problems hinder efficient utilization of ICA for batch-mode processing of EMEG datasets, and affect the objectivity of the results obtained with manual analysis.

The independent components sorting problem has received considerable attention in the fMRI data analysis literature. Gu et al. (2001) and Esposito et al. (2002) have introduced methods for component ordering based on spatial characteristics. Lu and Rajapakase (2003) suggested ranking based on component timecourse kurtosis. Himberg et al. (2004) used clustering of a succession of ICA realizations to select relevant components. Yang et al. (2008) describes a method for components selection based on the reproducibility principle. In the application of ICA to EEG and MEG, a technique based on measuring the amount of spatial component variance explained by the electromagnetic model was proposed by Grosse-Wentrup and Buss (2008). However, an accurate forward model is required to fully benefit from this approach.

In the current paper, we present a novel mutual information (MI) based approach for ICA components sorting. Moritz et al. (2003) has described a somewhat related method for component ranking, based on the spectral manifestation of stimulus periodicity. However, the periodicity assumption is not always fulfilled, especially in voluntary movement paradigms. In addition, the spectral measure uses only first and second order statistical moments, while our MI-based method implicitly employs higher order moments for estimating the amount of task-related signal present in a component. We report the performance of our new MI based approach and compare it against more conventionally used AP driven technique. Additionally, we demonstrate an application of the MI Spectrum to sorting InfoMax ICA components obtained from real MEG recordings obtained from an experiment designed to non-invasively map primary motor cortex (M1 zone).

2. Materials and Methods

2.1. EMEG Signal Model and Preliminaries

EMEG data recorded by a K − sensor array during the i-th repetition of a neuromotor or cognitive task can be written as the following linear combination

$\begin{array}{cc}{x}^i(t)=\left[a_1,\dots ,a_R\right]\text{}\left[\begin{array}{c}{f}_1^i(t)\\ ⋮\\ f_R^i(t)\end{array}\right]\text{}+\left[b_1,\dots ,b_L\right]\text{}\left[\begin{array}{c}{p}_1^i(t)\\ ⋮\\ p_L^i(t)\end{array}\right]\text{}+n(t) & (1)\end{array}$

For the i − th epoch, a multichannel signal at each instance of time, xⁱ(t) is a noisy additive mixture of source topographies [a₁, …, a_R] weighted by the corresponding stimulus-locked activation timeseries [fⁱ₁(t), …, fⁱ_R(t)], along with a task-unrelated contribution from sources with topographies [b₁, …, b_L] activated with task-unrelated timeseries [pⁱ₁(t), …, pⁱ_L(t)], and a random noise vector n(t). Topographies of task related sources form an R-dimensional signal subspace and topographies of task unrelated sources form an L-dimensional coherent interference subspace. ERP experiments are usually accompanied by a binary stimulus signal s(t) marking the task onset. In neuro-motor experiments this binary signal may be derived from the myographic activity record or from the accelerometer signal, using a thresholding procedure. Usually the goal of data analysis is to identify the task related signals and extract the task-related signal subspace to be used subsequently for neuronal source localization. For completeness, we may include induced sources whose activation power is locked to the task-onset moment with random phase. However, since we are interested in analysis of ERP's (which are phase-locked by definition), we do not include the induced component in (1).

Under the ideal and largely unrealistic conditions when activations fⁱ_r(t) are exactly reproducible across trials, time locked to the stimulus, and spatially coherent task-unrelated components are absent, the identification of task related components can be done optimally using the SVD of the stimulus locked average data matrix X = [x(0), …, x(T)], where $\overline{x}(t)={\displaystyle {\sum }_{i=1}^M{x}_i(t)}$, M is the task repetitions count, and T is the interval of interest duration (Vandewalle, 1988). The SVD yields the averaged data matrix decomposition X = USV^T. Columns of the orthonormal matrix U are the singular topographies, S is a diagonal matrix of SVs, and columns of the right singular matrix V are the singular activations. Task related components are chosen to be the first $\widehat{R}$ components ranked by power. $\widehat{R}$ is determined typically by visual analysis of the SV spectrum. Optionally, the SV spectrum of random matrix may be used as a reference in this task (Golub and Van Loan, 1996).

ICA is usually applied to the raw (unaveraged) spatial-temporal matrix X(t) and yields a spatial unmixing matrix B and a collection of independent components z_i(t) obtained as z(t) = BX(t), z(t) = [z₁(t), …, z_K(t)]^T. We assume that some of these components contain task-related signal and the others do not. In correspondence to B we can put matrix F = B⁻¹ so that X(t) = Fz(t). Columns of F are called independent topographies and describe the profiles formed by the corresponding independent “sources” on the sensors.

2.2. Mutual Information Spectrum

We propose to assess the degree to which the i-th component z_i(t) is related to the task using the normalized MI spectrum, computed as I_i = I(z_i(t), e(t)). e(t) is the expanded stimulus line signal, computed by convolution of the original binary stimulus signal s(t) with expansion kernel k(t) as e(t) = s(t) ∗ k(t) to produce monotonic variations over the interval of interest around each event onset moment. In this work we used a centered (i.e., symmetric around the x-axis) ramp function as the expansion kernel k(t).

We used a simple scaled histogram method to compute the MI as the difference between the entropy of an independent component z_i(t) and its entropy conditioned on the expanded stimulus line signal e(t), i.e.,

$\begin{array}{cc}{I}_0\left(z_i(t),e(t)\right)=H\left(z_i(t)\right)-H\left(z_i(t)|e(t)\right)& (2)\end{array}$

where H(u) denotes the entropy of u.

As suggested by Strehl (2002), we use the geometric mean of the two marginal entropy values to obtain the normalized MI quantities as

$\begin{array}{cc}{I}_i=I\left(z_i(t),e(t)\right)=\frac{I_0\left(z_i(t),e(t)\right)}{\sqrt{H\left(z_i(t)\right)H\left(e(t)\right)}}& (3)\end{array}$

We then define the MI spectrum as the rank ordered elements I_i : {I_i ≥ I_{i + 1}}.

As with the more conventional SV spectrum, visual analysis of the MI spectrum can be used to estimate the signal subspace rank. Originally suggested by us in Ossadtchi et al., (2000), this measure of MI with the expanded stimulus signal is a power-invariant way to assess the degree of task-relatedness of raw (unaveraged) timeseries.

We have introduced the notion of MI spectrum in the context of ICA component selection. This method can also be used for ranking singular components obtained via SVD of the stimulus locked average data matrix. To compute the MI spectrum for such singular components, first project the raw unaveraged data matrix X(t) onto the left singular column vectors of U as z(t) = U^TX(t) and then apply the MI spectrum calculation procedure as described above.

In neuro-motor tasks, the EMG signal, m(t), can be recorded and used instead of the expanded stimulus line. Then the MI spectrum is calculated as I_i = I(z_i(t), m(t)). Since the EMG signal usually occupies a broader spectrum than that of EEG or MEG signals, it is beneficial to perform a zero-phase-shift band-pass filter to remove excessively low and high frequency components prior to computing the MI.

2.3. Statistical Testing

In automated applications, and for a more informed decisions during the visual analysis of the MI spectrum, we suggest the following randomization testing scheme to estimate the p-values to reject the null-hypothesis that a component is not task related.

The suggested scheme is based on the observation that for signal components that contain a statistically significant evoked response, the value of the MI is directly related to the consistent correspondence (not similarity!) between the shape of this component and the expanded stimulus signal. Therefore, when the actual task onset moments are randomized, this correspondence will be destroyed. The MI values for the task-related components will experience a significant drop, while those that pertain to the task-unrelated components will remain in the original range.

We suggest the following simple steps to generate surrogate data and assess statistical significance of the observed MI values. In what follows M is number of independent components, J-number of randomization iterations, N_t is the number of samples in the stimulus signal s(t) and $N_r={\displaystyle {\sum }_{t=1}^{N_t}s(t)}$ is the number of task repetitions.

1. for j = 1:J

1. Create a new, surrogate stimulus signal s*(t) by randomizing task onset moments:

s*(t) = 0, ∀t ∈ [1, N_t];

for k = 1 : N_r, t ← U (0, N_t), s*(t) = 1, end_k

2. Calculate surrogate expanded stimulus signal: e*(t) = s*(t) ∗ k(t)

3. Calculate the amount of normalized MI of all the components with this surrogate expanded stimulus:

for i = 1 : M, I*_ij = I (e*(t), z_i(t)), end_i

2. end_j

The i-th row of I*_ij measures the MI for the i-th component and the j-th randomization of the stimulus onset signal. In order to calculate the p-values of the null-hypothesis that the i-th component does not contain task-related signal, we compute the fraction of surrogate values I*_ij, j ∈ [1, J] that exceed the actual observed value I_i. If we define a logical function L(a, b) = 1 if(a > b) and L(a, b) = 0 otherwise, then p-values for the null-hypothesis that the i-th component contains no task related signal can be expressed as $p_i=\frac{1}{M}{\displaystyle {\sum }_{j=1}^{j<J}L(I_{ij}^{\ast },I_i)}$.

2.4. Multiple Comparison Correction

Since we are testing several hypotheses, we need to correct the calculated p-values for multiple comparisons. Simple Bonferroni correction is appropriate, since the components are independent or at least orthogonal (SVD case). We therefore conclude that a component can be considered as task-related at the significance level of α if $p_i<\frac{\alpha }{N_c}$, where N_c is the total number of components tested. An example of applying the suggested statistical testing scheme is illustrated in Figure 2, where the proposed procedure allows for the correct identification of task related components in a simulated scenario with R = 2 task-related sources with powerful, spatially coherent, interference. For this and subsequent simulations, we used the procedure described in the following section.

FIGURE 2

Figure 2. Raw components' log(p)-values computed using the suggested randomization scheme. The horizontal line corresponds to Bonferroni corrected threshold determined for α = 0.05. The data were simulated with two dipolar sources and powerful spatially coherent interference. SVD was applied to the stimulus-locked data matrix that yielded first singular topographies. Subspace correlation of the first two singular topographies with the true signal subspace was [1, 0.987]. We can see that the suggested randomization scheme is able to correctly detect the first two components that span the signal subspace.

2.5. Simulation Procedure

In order to illustrate the performance of the MI spectrum and compare it with a more standard power driven approach, we performed realistic simulations with the following procedure. To simulate the observed sensor signals, we used Equation (1). We obtained a realistic source configuration from analysis of a somatosensory MEG dataset recorded with a 67 channel CTF MEG system. We applied the RAP-MUSIC localization algorithm to [0–200 ms] range of the stimulus-locked average data and obtained R = 2 dipoles with topographies a₁ and a₂ and their corresponding activations f₁(t) and f₂(t) [see Equation (1)]. We used these dipoles as sources of task related activity in our simulations. To simulate task related activation timeseries we adapted a kernel-based model of evoked potentials described by Lange et al. (1997). This model includes random trial-to-trial variation in the latency and amplitudes of signal components. It is based on the decomposition of activation timeseries into a superposition of Gaussian kernels with varying amplitudes and delays. The model is justified by the fact that, given relatively poor spatial resolution of MEG, the dipole timeseries may be viewed as the sum of activations of several neuronal assemblies, each with different intensity and activation latency values. A graphical example of such a decomposition is shown in Figure 3. The simulated activation of the r-th dipole during the i − th epoch can be expressed formally as $f_r^i(t)={\displaystyle {\sum }_{k=1}^K{\beta }_k^r}{v}_k^r(t-{\theta }_k^r)$ with K kernels defined as

$\begin{array}{cc}{v}_k^r\left(t\right)=f_r(t)\frac{e^{-\frac{{\left(t-{\tau }_k^r\right)}^2}{2{\sigma }_k^{r2}}}}{{\displaystyle {\sum }_{l=1}^K{e}^{-\frac{{\left(t-{\tau }_l^r\right)}^2}{2{\sigma }_l^{r2}}}}}, r\in \left[1,R\right], k\in \left[1,K\right]& (4)\end{array}$

The model incorporates random variables β_k, θ_k, k = {1, …, K} representing amplitude and latency variations. The latency jitter values were independent for all components and were generated using a Gaussian random variable with mean of 50 ms and standard deviation 10 ms.The k-th kernel amplitude variation β_k was modeled as normally distributed random variable with mean of unity and standard deviation equal to 0.2.

FIGURE 3

Figure 3. To simulate trial-to-trial variation of the responses we used overlapping Gaussian kernel based model. At each trial we varied relative positions of the kernel centers, kernel amplitudes, and global response latency with respect to the binary stimulus signal. A typical response and its representation with a set of overlapping modulated kernels is shown.

We modeled brain noise with L = 1000 spatially coherent, task-unrelated cerebral sources whose locations and time series varied with each realization. The corresponding topographies b_l were calculated using locally fitted concentric spheres MEG forward model as implemented in EMSE Software Suite, Source Signal Imaging Inc., San Diego, CA, USA. The activation time series were narrow-band signals obtained via zero-phase filtering of realizations of Gaussian (pseudo)random process by the fifth order band-pass IIR filters in the bands corresponding to theta (4–7 Hz), alpha (8–12 Hz), beta (15–30 Hz) and gamma (30–50 and 50–70 Hz) activity. Their relative contributions were scaled in accordance with $\frac{1}{f}$ characteristic of the realistic EMEG spectrum. An additional narrow-band alpha-component (9–11 Hz) of occipital origin ([−0.05, 0.01, 0.06] in EMSE coordinate system) was also included. We scaled the brain noise components to match typical signal-to-noise ratio of real-life recordings.

3. Results

3.1. Discriminative Power of the MI Spectrum

Consider the situation when a task-related signal is generated by a pair of dipolar sources. When a pair of sources has highly correlated topographies or, in case of a large imbalance in source magnitude, the second singular component may be obscured and may not produce a pronounced SV distinguishable from the baseline. In this case, analysis of the SV spectrum will fail to provide the correct estimate of the signal subspace dimension. An example is illustrated in Figure 1, where the SV spectrum of the averaged data matrix obtained from a simulated dataset with R = 2 task-related sources does not have a significant drop between R = 2 and 3. On the contrary, the MI spectrum exhibits a very clear separation between task-related and task-unrelated parts, as it can be seen in the bottom panel of Figure 1.

In order to perform a more systematic evaluation of using MI to measure the extent to which a component is task-related, we performed a set of simulations with two dipolar sources in the presence of realistic brain noise. We varied the ratio of activation amplitudes of the two dipoles, performed SVD of the averaged data matrix and calculated the MI spectrum for the projections of continuous data onto the left singular vectors. We then compared the discriminating power of the MI and the SV spectra. To do so we introduced the discriminating indicator q. Since the correct rank value is R = 2 we used $q=log\left(\frac{S(2)}{S(3)}\right)-log\left(\frac{S(3)}{S(4)}\right)$, see Figure 1. q is sensitive to the drop between the second and the third components, referenced to the ratio of the two largest noise range spectrum values (with indices 3 and 4) immediately following the two signal components (with indices 1 and 2). Results are shown in Figure 4, illustrating the discriminating indicator q as a function of source amplitudes ratio. We performed this numerical experiment for a varying number of trials in a simulated dataset. We found that the proposed MI spectrum outperforms the SV spectrum for all trial counts, and also provides for a clearer seperation between the task-related and task-unrelated components. Also note that in most cases the correlation of the subspace spanned by the first two singular topographies and the true signal subspace spanned by a₁ and a₂ was sufficiently high to be considered as a correct estimate of the simulated signal subspace.

FIGURE 4

Figure 4. Discriminating power indicator q as a function of source amplitudes ratio for different number of trials calculated for MI and SV spectra. Each curve corresponds to a fixed datapoints count. We can see that the proposed MI based measure outperforms power based technique and produces a clearer cut between the task-related and task-unrelated components.

3.2. Receiver Operating Characteristics of MI

In this section we describe our experiments on exploring receiver operating characteristics (ROC) of the MI metrics. We consider the task of discriminating between the components that contain task-related signal and those that do not. Spatial components z_i(t) obtained from signals that can be represented using Equation (1) can be viewed as the superposition of signal and noise, expressed as

$\begin{array}{cc}z(t)=af(t)+{\sigma }_{p}p(t)+{\sigma }_{n}n(t),& (5)\end{array}$

where f(t) is task-related source activity of amplitude a. p(t) and n(t) are the contributions from spatially coherent and spatially white noise sources, with standard deviation values σ_p and σ_n respectively.

We simulated repetitions of the task-related signal, including the jitter and variations characteristic of realistic brain signals. For each of N_mc = 1000 Monte-Carlo realizations, we simulated N = 100 signals according to Equation (5). N₁ = 10 out of 100 signals had a = 1 and the remaining N₀ = 90 signals had a = 0, i.e., no task-related signal present. We simulated brain noise as described in the Simulations Procedure section. The goal was to detect the components that contain task related signals. We compared the MI values against the more traditionally used stimulus-locked averaged signal power (AP), calculated as

$\begin{array}{cc}{P}_i={\displaystyle \sum _{t\in W}{\overline{z}}_i^2(t)},& (6)\end{array}$

where z_i(t) is the stimulus-lock averaged i-th signal. The summation was performed over a 200 ms window centered ont the stimulus. We used the same windows to calculate both MI and AP measures.

To calculate the ROC curves, we applied thresholding to the MI and AP spectra separately, and marked as detected only those components whose corresponding MI or AP values exceeded the threshold. The threshold was originally chosen to be 0.05 of the largest value in the spectrum (AP or MI). In order to obtain the ROC curve we calculated the sensitivity $p_{sens}(\theta )=\frac{\text{NTP}}{N_1}$ and specificity $p_{spec}(\theta )=1-\frac{\text{NFP}}{N_0}$ for a succession of evenly spaced threshold values θ = 0.05 k max_i(I_i) or θ = 0.05 k max_i(P_i) for k = 1, …, 19.

The result is shown in Figures 5A,B. For all epoch counts, the MI based measure significantly outperforms the power-based characteristic, and also provides better sensitivity for any selected specificity. We have observed similar behavior when dealing with real MEG data, as described below in the “MI versus AP for small epoch counts in real data” section. The improved ROC may be explained by the fact that the MI based measure implicitly takes into account higher order statistical information as compared to the power based approach, where only the first and second order statistical moments are used.

FIGURE 5

Figure 5. (A) The family of ROC curves for Averaged Response Power and Mutual Information based detection for various number of epochs on the same plot. MI based detection clearly and significantly outperforms the conventional method. Even with 50 trials the MI based criterion (“Nrm. MI, 50 trials” curve) allows to achieve 70 percent of sensitivity with ideal specificity. We can see that for all counts of epochs the MI based measure significantly outperforms the power based characteristic. This can be explained by the fact that in calculation of MI we implicitly take into account higher order information as compared to the power based approach where only the first two statistical moments are used. (B) Area under ROC curve performance characteristic for Averaged Response Power and Mutual Information.

3.3. Application to M1 Mapping

Reliable mapping of the primary motor cortex (M1) based on functional neuroimaging provides an important complement to the use of structural data alone. However, since various zones forming the somatosensor complex appear to be in a coupled interaction even in the motor planning stage, the localization of M1 zone from the functional EEG and MEG data via standard approaches is problematic and often does not yield reliable results (Sanders et al., 1996; Gerloff et al., 1998).

Inspired by the work of Riehle (2005) we explored the possibility of using the information from activation timecourse morphology and looked for spatial properties of activations with sharp non-linear increase just preceding the movement onset. To do so we studied MEG-recorded brain responses during a voluntary index finger movement task performed by 18 healthy right-handed volunteers.

For computational feasibility we used a subset of 50 sensors located over the left sensory-motor region. These were selected based on the grand-average responses, as shown in Figure 6. Recordings from all experimental sessions in all subjects were concatenated into a single sequence and decomposed using the InfoMax-ICA approach. We obtained 50 independent components and ranked them according to the amount of MI with the expanded stimulus signal e(t) = s(t) ∗ k(t) in the 200 ms window centered around the movement onset moment marked by s(t). The choice of the time window is motivated by our interest in the early components reflecting the activity M1 zone.

FIGURE 6

Figure 6. We selected for ICA analysis 50 electrodes over the left sensory-motor cortex corresponding to the cortical representation of the right hand. Evoked response fields are shown. The selected sensors are framed.

We then focused on the first two components with the largest mutual information, as shown in Figure 7. The temporal dynamics of the first two components showed a slow activation increase starting as early as 400 ms before the actual movement onset. However, as illustrated in Figures 7A,B, the two components differed in their behavior during the interval directly preceding the movement onset. The component with the larger value of MI exhibited a sharp quasi-exponential growth starting at around 50–70 ms before the movement onset. The onset dynamics of the second component was smooth, corresponding to a quasi-linear growth. After movement onset both components show a pronounced negative deflection reaching a minimum at around 100 ms after the movement onset.

FIGURE 7

Figure 7. Relative contribution of independent components of the evoked magnetic activity accompanying movement onset. All derived components (abscissa) are ranked according to the amount of mutual information with the expanded stimulus signal. The MI values were normalized to the overall amount of information (ordinate). (A,B) The stimulus-locked average time courses of the two components (4 and 1) with the greatest share of mutual information with expanded stimulus signal. Abscissa: time (ms) in relation to the movement onset. Ordinate: component amplitude (arbitrary units). Low-pass filter with 30 Hz cutoff frequency was applied. Movement onset is shown by the vertical dotted line. Horizontal dotted line represents the background signal level. Dark-gray area shows time interval of exponential growth preceding movement onset of component 4.

We used MNE distributed source imaging to localize the neuronal generators underlying the topographies of the first two components with the highest values of MI. In both cases we observed activations in the area surrounding the central fissure of the hemisphere contralateral to movement, shown in Figure 8. After thresholding at p < 0.001 we observed that cortical areas subserving these two components do not overlap, as shown in Figure 8. Cortical sources of the first component localized primarily on the anterior slope of the central sulcus superior to the omega zone, shown in Figure 8C. This source most likely lies in M1, based on the anatomy. The cortical sources for the second component were located in the post-central sulcus and in the depth of the central sulcus inferior to the omega shaped zone, shown in Figure 8D.

FIGURE 8

Figure 8. Source localization of the two independent components with the largest amount of MI, (A) component 4, (B) component 2. Inflated surface of the left hemisphere obtained by averaging of the co-registered individual cortices of 18 subjects is shown. Light-gray areas represent gyri and dark-gray areas represent sulci. The central sulcus is indicated by arrows. Localization of the sources of both components was done via MNE approach followed by Bonferroni correction. (C,D) Illustrated patches correspond to current distribution thresholded at significance level of p < 0.001. (C) We can see that in agreement with our hypothesis the topography of the independent component with the largest amount of MI predominantly localizes anterior to the central sulcus and may represent the M1 zone. (D) The second component's sources localize to the post-central fissure and in the depth of the central fissure inferior to the omega shaped zone.

3.4. MI Versus AP for Small Epoch Counts in Real Data

We also compared the performance of the MI spectrum with the more conventional AP metric (6) when the number of epochs is limited. We took every 30th event and analyzed the data according to the scheme described above. Independently sorted MI and AP spectra for the first 15 components are shown in Figures 9A,B respectively. The MI spectrum clearly shows the presence of task-related signal in the first two components with original indices 4 and 1. The AP spectrum shows five seemingly task-related components (indices 8, 3, 4, 15, 1) standing out from the baseline. Components 1 and 4 are identified as task-related by the both measures of task-relatedness.

FIGURE 9

Figure 9. Analysis of a subset of trials. Normalized by the first element value MI (A) and AP (B) spectra calculated using a subset of events for the first 15 independent components. On the top of each stem actual (unsorted) component number. As prescribed by the MI spectrum we should select components 4 and 1 as task related. AP spectrum suggests components 8, 3, 4, 15, 1 for this role. In order to check which of the two measures provide the right answer we inspected stimulus-locked averages of the components suggested by the two measures, see Figure 10.

In order to check which of the two methods provided the correct answer, we performed stimulus-locked averaging of the first five components obtained by sorting in decreasing order the MI and AP spectra, shown in Figure 9. The results are shown in Figure 10. The first two components (4 and 1) identified by the MI spectrum (see Figure 9A) show a clear task related deflection. The remaining components do not have significant amount of stimulus-locked activity and therefore are most likely unrelated to the task. Three out of five components identified by the AP spectrum (first two and the fourth) do not exhibit any deflection resulting from coherent summation. Note also that components 1 and 4 are among the five components selected by the AP spectrum (the third and the fifth). Visual analysis of the averages obtained for the other three components suggested by the AP spectrum does not reveal the presence of significant amount of stimulus-locked activity in three out of five components.

FIGURE 10

Figure 10. In order to check which of the two methods provide the right answer we performed stimulus-locked averaging of the first five components in the order prescribed by sorting of MI and AP spectra (Figure 9). Left panel corresponds to MI prescribed components and the right panel - to AP. As it can be seen the first two components (4 and 1) emphasized by the MI spectrum (Figure 9A) show a clear task related deflection, the subsequent components do not have significant amount of stimulus-locked activity and therefore are most likely unrelated to the task as correctly indicated by the characteristic drop in the MI spectrum (Figure 9A). Three out of five (8, 3, 15) component averages prescribed by the AP spectrum do not exhibit the expected deflection.

Based on these observations, we conclude that both MI and AP spectra demonstrate identical sensitivity, as both were able to detect two clearly task-related components (1 and 4) that were also found using the full dataset and characterized by inverse modeling (see Figure 8). However, the MI spectrum exhibits optimal specificity, identifying two components (Figure 9A). Both of these components appear to have a task-related deflection in their stimulus-locked averaged profiles (Figure 10A). The specificity of AP based measure r is poor by comparison with the MI measure, since AP it identified 5 components (Figure 9B), including 3 false positives and 2 correct hits (Figure 10B). The observed behavior is consistent with our simulation studies, illustrated in Figures 5A,B, where the MI spectrum demonstrated significantly higher ROC characteristics, and provided higher specificity for any fixed sensitivity value compared to the AP spectrum.

4. Discussion

We describe a novel information-theoretic approach for spatial components ranking. Our method is based on the MI Spectrum which serves as a power-invariant measure of repetitive task-related signal in the temporal loadings of spatial components. Using realistic simulations we demonstrated that the task-relatedness measure, based on estimating the MI between a component and the expanded binary stimulus signal, allows for significantly higher detector characteristics when compared with conventional alternatives. It also provides a means for more clear-cut separation of task-related and task-unrelated components when compared with the standard power driven approach that is used in SVD, and sometimes used for ranking ICA components as well. The MI measure can be used for sorting the components obtained from any sort of spatial decomposition, as long as it is possible to calculate the quasi-continuous timeseries underlying the components of interest. The demonstrated advantage in performance over the power-driven measure makes the MI spectrum method a candidate for the routine use in ranking both SVD and ICA components in the analysis of ERP data. Since the MI method is insensitive to powerful non-task-related noise sources, it should also facilitate automatic unsupervised analysis of ERP data using ICA.

The method can be easily extended to extract not only the evoked (phase-locked to the stimulus) activity but also band-specific task induced activity that is characterized by random phase but stimulus-locked power fluctuations. Such an extension would require that the band-pass filtered components envelope should be calculated before MI spectrum estimation.

We have also investigated the MI method performance applied to an MEG dataset in a voluntary finger movement task. Such paradigms present special challenges, since they include large amount of random latency jitter when compared with an external stimulus driven paradigms. This increased jitter comes from the inevitable errors in the estimates of motion onset obtained from the accelerometer signal. Nevertheless, the MI measure supported a clear cut separation of four task related components (Figure 7). The component with the largest MI (index 4) demonstrated a non-linear increase of activation just prior to the motion onset. In agreement with the previous experiments on primates (Riehle, 2005), this component localized primarily to the anterior slope of the central sulcus superior to the omega zone (Figure 8C), and most likely originates in M1. Thus, we demonstrated the potential to localize M1 non-invasively on a group level, using a functional probe.

In order to compare the performance of the MI and power based measures on the experimental dataset, we used a reduced number of trials. As shown by simulations (Figures 4, 5), this reduction should increase the contrast between the performance characteristics of the two methods. It should also mimic more realistic scenarios, when only a single subject dataset is used for ICA analysis. Under these conditions, we demonstrated that the proposed MI based measure for a fixed sensitivity value yields significantly higher specificity than the more conventional power based measure. This result is in agreement with our simulation studies.

In the current work we used a simple histogram-based approach for calculation of MI, omitting any bias correction. For realizations of independent random processes Treves and Panzeri (1995) have shown that MI estimate bias is quadratically proportional to the number of histogram bins used to approximate the pdf of continuous random processes and is inversely proportional to the number of datapoints. Since we used a large number of datapoints (N ≈ 6.5 × 10⁵) in our experimental data analysis compared to the K = 10 bins used for approximation of the probability density functions, we do not expect a bias correction procedure to appreciably alter the observed MI. However, it has also been shown (Chrisman, 2013) that the bias decreases as the true MI between the timeseries pairs grow. This means that bias may result in MI values of task-unrelated components being overestimated, yielding a decreased contrast between the task-related and task-unrelated components in the MI spectrum. In our simulation studies we used a relatively small number of datapoints compared to a standard EEG/MEG data recording per single patient. Therefore, we expect that the observed performance (Figures 4, 5) may be further improved with a proper bias correction procedure. The use of a biased estimator in the statistical testing approach we implemented results in less sensitive tests, since the null-hypothesis distribution estimate appears to be “shifted to the right.” Selection of an appropriate bias correction method, however, requires a significant amount of additional numerical experiments and goes beyond the scope of this paper.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Dr. Richard Greenblatt, Dr. Alexandr Valjamae and the reviewers for the careful reading of the manuscript and their critical comments that helped to improve the readability and added to the clarity of presentation. This study was partly supported by the Russian Targeted Federal Program “Scientific and scientific-pedagogical personnel of innovative Russia” (Agreement 8488), President Grants for Government Support of Young Russian Scientists (5137.2013.4), Russian foundation for basic research (grants 11-06-12037-ofi-m-2011 and 11-06-00449-a), Basic Research Program of the National Research University Higher School of Economics, SPBU grant for priority research 0.37.522.2013 “Neuroeconomics and managerial decision making: interdisciplinary study,” by the funds provided to Alexei Ossadtchi by Source Signal Imaging Inc. as well as Alexei Ossadtchi's personal funds.

References

Child, D. (2006). The Essentials of Factor Analysis. 3rd Edn. London; New York, NY: Continuum.

Chrisman, L. (2013). Estimation of mutual information. Available online at: http://blog.lumina.com/2013/estimation-of-mutual-information/

Comon, P. (1994). Independent component analysis: a new concept? Signal Process. 36, 287–314. doi: 10.1016/0165-1684(94)90029-9

CrossRef Full Text

Delorme, A. (2010). Headit and eeglab: open resources for electrophysiological discovery. Psychophysiology 47:S6. [15th Annual Meeting Portland Marriott Downtown Waterfront (Portland, OR)].

Esposito, F., Formisano, E., Seifritz, E., Goebel, R., Morrone, R., Tedeschi, G., et al. (2002). Spatial independent component analysis of functional MRI time-series: to what extent do results depend on the algorithm used? Hum. Brain Mapp. 16, 146–157. doi: 10.1002/hbm.10034

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gerloff, C., Richard, J., Hadley, J., Schulman, A. E., Honda, M., and Hallett, M. (1998). Functional coupling and regional activation of human cortical motor areas during simple, internally paced and externally paced finger movements. Brain 121(Pt 8), 1513–1531. doi: 10.1093/brain/121.8.1513

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Golub, G., and Van Loan, C. (1996). Matrix Computations. 3rd Edn. Baltiomore; London: Johns Hopkins University Press.

Grosse-Wentrup, M., and Buss, M. (2008). Multiclass common spatial patterns and information theoretic feature extraction. IEEE Trans. Biomed. Eng. 55, 1991–2000. doi: 10.1109/TBME.2008.921154

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gu, H., Engelien, W., Feng, H., Silbersweig, D. A., Stern, E., and Yang, Y. (2001). Mapping transient, randomly occurring neuropsychological events using independent component analysis. Neuroimage 14, 1432–1443. doi: 10.1006/nimg.2001.0914

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Himberg, J., Hyvarinen, A., and Esposito, F. (2004). Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage 22, 1214–1222. doi: 10.1016/j.neuroimage.2004.03.027

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hyvärinen, A., Hoyer, P., and Inki, M. (2001). Topographic independent component analysis. Neural Comput. 13, 1527–1558. doi: 10.1162/089976601750264992

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jung, T., Makeig, S., Humphries, C., Lee, T., McKeown, M., Iragui, V., et al. (2000). Removing electroencephalographic artifacts by blind source separation. Psychophysiology 37, 163–178. doi: 10.1111/1469-8986.3720163

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lagerlund, T., Sharbrough, F., and Busacker, N. (1997). Spatial filtering of multichannel electroencephalographic recordings through principal component analysis by singular value decomposition. J. Clin. Neurophysiol. 14, 73–82.

Pubmed Abstract | Pubmed Full Text

Lange, D., Pratt, H., and Inbar, G. (1997). Modelling and estimation of single evoked brain potential components. IEEE Trans. Biomed. Eng. 44, 791–799. doi: 10.1109/10.623048

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lu, W., and Rajapakase, J. (2003). Eliminating indeterminacy in ICA. Neurocomputing 50, 271–290. doi: 10.1016/S0925-2312(01)00710-X

CrossRef Full Text

Makeig, S., Bell, A., Jung, T., and Sejnowski, T. (1996). “Independent component analysis of electroencephalographic data,” in Advances in Neural Information Processing Systems 8: Proceedings of the 1995 Conference, Vol. 8, eds D. Touretzky, M. Mozer, and M. Hasselmo (Cambridge; London: MIT Press), 145–151.

Misulis, K. (1994). Spehlmann's Evoked Potential Primer: Visual, Auditory, and Somatosensory Evoked Potentials in Clinical Diagnosis. 2nd Edn. Oxford: Butterworth-Heinemann Medical.

Moritz, C. H., Rogers, B. P., and Meyerand, M. E. (2003). Power spectrum ranked independent component analysis of a periodic fMRI complex motor paradigm. Hum. Brain Mapp. 18, 111–122. doi: 10.1002/hbm.10081

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ossadtchi, A., Baillet, S., Mosher, J., and Leahy, R. (2000). “Using mutual information to select event-related components in ICA,” in Biomag 2000: Proceedings of the 12-th International Conference on Biomagnetism (Espoo: Helsinki University of Technology), 199.

Riehle, A. (2005). “Preparation for action: one of the key functions of motor cortex,” in Motor Cortex in Voluntary Movements: A Distributed System for Distributed Functions, eds A. Riehle and E. Vaadia (Boca Raton, FL: CRC Press), 213–240.

Sanders, J. A., Lewine, J. D., and Orrison, W. W. (1996). Comparison of primary motor cortex localization using functional magnetic resonance imaging and magnetoencephalography. Hum. Brain Mapp. 4, 47–57. doi: 10.1002/(SICI)1097-0193(1996)4:1<47::AID-HBM3>3.0.CO;2-P

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Strehl, A. (2002). Cluster ensembles a knowledge reuse framework for combining multiple partitions. J. Mach. Learn. Res. 3, 583–617. doi: 10.1162/153244303321897735

CrossRef Full Text

Tong, L., Liu, R.-W., V, C. S., and Huang, Y.-F. (1991). Indeterminacy and identifiability of blind identification. IEEE Trans. Circuits Syst. 38, 1–13. doi: 10.1109/31.76486

CrossRef Full Text

Treves, A., and Panzeri, S. (1995). The upward bias in measures of information derived from limited data samples. Neural Comput. 7, 399–407. doi: 10.1162/neco.1995.7.2.399

CrossRef Full Text

Vandewalle, J., and De Moor, B. (1988). “A variety of applications of sigular value decomposition in identification and signal procesing,” in SVD and Signal Processing, Algorithms, Applications and Architectures, ed E. Deprettere (Amsterdam: Elsevier), 43–91.

Vigario, R., Sarela, J., Jousmaki, V., Hamalainen, M., and Oja, E. (2000). Independent component approach to the analysis of eeg and meg recordings. IEEE Trans. Biomed. Eng. 47, 589–593. doi: 10.1109/10.841330

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yang, Z., LaConte, S., Weng, X., and Hu, X. (2008). Ranking and averaging independent component analysis by reproducibility (RAICAR). Hum. Brain Mapp. 29, 711–725. doi: 10.1002/hbm.20432

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text