Detection of Artifacts and Brain Responses Using Instantaneous Phase Statistics in Independent Components

1. Introduction

Multichannel recordings of magnetoencephalography (MEG) are usually comprised of repetitive events (e.g. external stimuli) in order to evoke the relatively weak magnetic fields of brain responses to a specific task. The analysis of the underlying electrophysiological brain activity of such unaveraged signals is notoriously challenging. During MEG experiments environmental and other external noise sources derogate the signal of interest. Furthermore, brain activity which is not involved in the task processing and therefore not of prime interest (often termed as “brain noise”) also interfere with the weak brain responses.

In order to increase the signal-to-noise ratio (SNR) of the recorded data a widely common strategy is to perform signal averages. Assuming white

In analogy to white light (visible light), where all types of visible wavelengths are represented, the term „white“ is used for types of signal (or noise), where the power spectral density of the signal is uniform. This means the power spectral density of the signal is the same at all frequencies.

noise that is temporally uncorrelated across trials, the SNR improvement gained by the averaging process over N trials scales theoretically with √N. In practice the noise reduction is a little less than 1/√N since, the evoked activity usually varies in its signal strength over time. An important aspect to bear in mind when performing averages is, that it will reveal the most prominent neuromagnetic correlates of brain responses, only. In this way information about weaker brain activity is largely suppressed, especially when multiple strong and weaker sources acting in a coordinated manner. Additionally, the temporal dynamics of each individual response will not be preserved. In contrast, single-trial analysis retains the temporal dynamics of the neuromagnetic responses, but suffers from poor signal-to-noise ratio (SNR), and is therefore rarely applied.

Multichannel MEG recordings are usually comprised of a mixture of the underlying brain activity and field components originated from noise and artifact sources. In MEG, as well as in electroencephalography (EEG), the most prominent biological artifacts originate from eye blinks/movements (ocular artifacts, OA), heart beats (cardiac artifacts, CA) and muscle activity (MA). The signal strength of such biological artifacts may be several orders of magnitude higher than the signal of interest. Therefore, the analysis of MEG/EEG signals requires the identification and elimination of the artificial signals prior to analysis. Independent component analysis (ICA), for example, is widely used to separate brain signals from noise and artifact components. Effective removal of noise and biological artifacts was reported in various publications (Dammers et al., 2010, Escudero et al., 2011, Mantini et al., 2008, Nikulin et al., 2011, Ting et al., 2006).

Most of the proposed ICA based source separation methods are performed to identify and exclude artifacts from the recorded signal, while only a few studies have been performed to elicit brain responses from the decomposed signals (Dammers et al., 2010, Hild & Nagarajan, 2009, Nikulin et al., 2011). Reliable classification of noise, artifacts and brain activity in decomposed MEG signals is not an easy task to do. However, separating brain responses due to an external stimulus from unwanted signal contributions (i.e., the remaining signal components which are not related to the signal of interest, such as brain noise) in unaveraged multichannel MEG recordings, is even more challenging due to the weak signal strength. Artifact signals, such as eye blinks, eye movements, muscle artifacts or the field contribution from the heart are types of signals that are usually much stronger in amplitude than brain signals are. This is one of the reasons why most of the applied artifact identification routines are based on amplitude statistics of the decomposed MEG signal. In previous reports, it has been shown that in the case of cardiac artifacts the ICA decomposition may reveal multiple independent components (ICs), where peak amplitudes of the cardiac artifact components vary quite substantially across the components, and is therefore, difficult to identify by amplitude based methods (Dammers et al., 2008, Sander et al., 2002).

When trying to identify brain responses from independent components different strategies are indispensable. In the last decade, the number of channels in whole head MEG systems has increased quite substantially. Meanwhile modern MEG systems are comprised with about 250 or more channels to cover most of the human head. For source reconstruction a dense inter-channel spacing has a favorable effect on the spatial resolution of the localized activity. On the other hand, higher channel densities require smaller pick-up loops (i.e., the size of the field measuring coil is smaller), which in turn has the effect of measuring a smaller amount of underlying brain activity. If we assume the same amount of environmental noise at each channel, the measured signal will have a smaller SNR when the size of the pick-up loop decreases.

In particular, when

Magnetometers pick-up coils are in principle much more sensitive to magnetic field changes than gradiometers are, and therefore, they are also more sensitive to noise.

magnetometers are used to detect the neuromagnetic activity, such brain responses must be really strong to be identified by visual inspection. Therefore, the high dimensionality of the data and the poor signal-to-noise ratio (SNR) are the two most inviting challenges in MEG single-trial analysis. In order to overcome these problems different strategies have been proposed. Guimaraes and colleagues, for example, used well-known classification methods to identify brain responses in MEG raw signals or signal decompositions (Guimaraes et al., 2007). In this study single-trial classification of MEG recordings have been applied utilizing two different classifier, the linear discriminant classification (LDC) and support vector machine (SVM). Both strategies were applied to classify single-trial brain responses to auditory and visual presentations of words. Another and novel approach utilizing the frequency domain for the identification of neural oscillations from MEG measurements based on spatio-spectral decompositions was recently introduced by Nikulin and colleagues (Nikulin et al., 2011). The method optimizes the spectral SNR by maximizing the signal power at certain peak frequencies, while simultaneously minimizing the power at the neighboring frequency bins. In this way, the method allows to identify components with a “peaky” spectral profile, which is typical for oscillatory processes (Nikulin et al., 2011).

A different way of detecting strong and weak oscillatory neuromagnetic activity is to incorporating the phase dynamics of the data into the analysis (Gross et al., 2006, Matani et al., 2010, Schyns et al., 2011; Tass et al., 1998, Vairavan et al., 2010, Yoshino et al., 2002). Numerical studies have shown that the analysis of phase dynamics reveals synchronization patterns even in cases of weak or absent amplitude correlation (Rosenblum et al., 1996). In 2008, Dammers and colleagues introduced a fully automated approach using cross trial phase statistics (CTPS) for the identification of cardiac artifacts in phase space from decomposed MEG signals (Dammers et al., 2008). In this study, the authors demonstrated that CTPS is a highly sensitive method, which is capable of identifying not only strong but also weak cardiac artifacts that hardly can be identified by visual inspection. Recently, the same method has been applied to extract brain responses that are highly correlated to an external stimulus from a set of decomposed MEG signals in a user-independent fashion (Dammers et al., 2010). Using a single statistical measure the method automatically extracts the components of interest, with respect to repetitive events. When focusing on brain responses, CTPS leads to an enhanced signal-to-noise ratio and therefore enables single-trial data analysis including source localization with improved spatial resolution.

2. Artifact or signal of interest?

According to the experimental design and the neuroscientific question behind the MEG investigation the number and types of stimuli may vary for each experiment. Quite often, within one MEG experiment several conditions (i.e., different stimuli are presented to invoke different brain areas) are employed, where the subject has to cope with different tasks during the experiment. Depending on the duration of the experiment and the number of different conditions used within one experiment, the number of stimuli for each condition may be in the range of a few tens only. In such cases, averaging is less effective to reduce the noise. In addition to noise, the measured activity in all MEG experiments represents a signal mixture comprised of field contributions originating from different brain areas and artifact sources, such as, ocular, muscle and cardiac activity. In addition, the signal strength of such biological artifacts may be several orders of magnitude higher than the signal of interest. Although the signal strength of ocular artifacts (e.g., eye blinks and eye movements) is less compared to the signal of cardiac activity, ocular artifacts are one of the strongest biological artifact signals in MEG due to the close eye-to-sensor distance. Contrary to ocular and muscle artifacts (e.g. swallow and mimic), the cardiac signal is more frequent and cannot be avoided or prevented. In unfavorable cases, the artifact signal (e.g. ocular or cardiac activity) may occur just around the time of the stimulus. It is therefore advisable, and whenever possible, to avoid stimuli intervals, which do have the same frequency rate as the heart beat. However, even when the experimental design includes variable stimuli intervals (often called as stimulus onset asynchrony, SOA) in many MEG experiments it is observed that subjects do perform eye blinks straight after a visual or auditory stimulation. In terms of source localization, this is a severe problem. The artifact signal will be highly correlated with the signal of interest so that even averaging of hundreds of trials may not prevent the source location being shifted towards to source of the artifact. In Fig. 1 the upper row shows strong magnetic field contributions, which are highly correlated with the onsets of auditory stimuli as indicated by the blue vertical lines. Superposition of the EOG signal reveals that the recorded signal is related to ocular artifacts (eye blinks). The bottom row of the figure shows an example of strong but stimulus uncorrelated field contributions from cardiac activity. In either case, the analysis of MEG signals requires the identification and elimination of the artificial signals prior to analysis.

Figure 1.

Signal of interest or artifact?In a) strong magnetic field contributions appear to be highly correlated with the onsets of the auditory stimuli (blue vertical lines). A superposition of the EOG signal (in orange) reveals that the recorded signal is related to ocular artifacts. In b) an example of strong but stimulus uncorrelated field contributions from cardiac activity (in red) is shown.

3. Methodological concepts of signal separation and classification

In the following sections we will give a brief introduction to the basic principles of signal decomposition utilizing independent component analysis (ICA), followed by an introduction to the basic concepts of cross trial phase statistics (CTPS).

3.1, A short introduction to independent component analysis (ICA)

As described in the previous section the recorded signal at each detector consists of field components originating from the acting brain, noise and artifact sources. A variety of methods for the identification and separation of noise and artifacts from the signal of interest have been proposed over the last ten years to overcome this problem. The rejection of corrupted epochs or trials by visual inspection is however still widely applied. The major disadvantage of this method is that it results in loss of data. Apart from the fact that visual inspection cannot be objective, if the number of trials for averaging is low or the analysis must be performed on single-trials, a rejection of corrupted epochs is not an option.

Independent component analysis (ICA) is widely used to separate brain signals from artifact components using utilizing both semi- and fully automated procedures (Dammers et al., 2008, Escudero, et al., 2010, 2011; Hironaga & Ioannides, 2007, Li, et al., 2006, Mantini et al., 2008, Ossadtchi et al., 2004, Ting et al., 2006). Independent component analysis was developed for solving the blind source separation (BSS) problem with the basic assumption, that the recorded data in a sensor array are linear sums of temporally independent components originating from spatially fixed sources (Bell & Sejnowski, 1995, Comon, 1994, Herault & Jutten, 1986, Hyvärinen, 1999). The general idea of ICA based source analysis methods is to separate the signal of interest from unwanted noise and artifact signals, before source localization is applied to a few number of selected components. The challenge here is to recover N independent source signals s(t) = [s ₁(t), s ₂(t),..., s _N (t)] ^T from N linearly mixed observed signals x(t) = [x₁(t), x ₂(t),..., x _N (t)] ^T recorded by N detectors. Let A denote the mixing matrix the system to be solved can be formulated as

$$x=As$$ E1

If M is the number of data points in x (t) at time t using N detectors the mixing matrix A will be of dimension NxN. The key problem in ICA is to find an unmixing matrix W (similar to a pseudo-inverse A^-1, with W≈A^-1) while imposing that the sources in s are statistically independent

$$c=Wx$$ E2

ICA transforms an N-sensor data array into an N-dimensional component space, where each of the components in c carries a minimum amount of mutual information with respect to the temporal dynamic and thus is maximally independent. The independent components are stored in rows of c. Assuming the model is valid each row in c reflects the time courses of the true source, except that scale, sign and order will not be preserved. For solving equation (2), the method makes no assumption about the mixing process except that it is linear including the restriction that the number of sources is less or equal N. The estimated full rank square matrix W in (2) can be treated as a spatial filter matrix which linearly inverts the mixing process. It is important to note that independence is not equivalent to uncorrelatedness, as it would be for Gaussian distributions. In fact estimating independent components is based on maximization of the non-gaussianity of Wx (Hyvärinen & Oja, 2000). However, if A is sufficiently determined (i.e., by solving W^-1) one can remove the mixing effect by applying W to x.

After signal decomposition the problem now is to objectively identify components of interest in such high dimensional data sets, as it is the case for modern MEG systems. Different strategies have been proposed to exactly tackle the problem of unsupervised component extraction to find either artifacts (Dammers et al., 2008, Escudero et al., 2011; Li et al., 2006, Mantini et al., 2008) or signal of interest (Dammers et al., 2010, Hild & Nagarajan, 2009, Nikulin et al., 2011, Vairavan et al., 2010).

Apart from time series analysis applied to the independent components for the identification of the source of interest, it might be helpful to look at the so-called place map of each temporally independent component. Such maps can be constructed from the columns of W^-1. Since each column of matrix W^-1 has N entries a projection of the signal onto the sensor arrangement gives rise to the location of the detectors showing strong signal contributions for the associated component.

Once all components of interest are identified, the separation from unwanted components is done by zeroing all columns of W^-1 except those which correspond to the signal of interest. The zeroing results in a new matrix Ŵ^-1 which filters out the unwanted source contributions, while leaving the signal of interest. The ICA filtering (i.e., artifact removal or extraction of components of interest) of the measured data is then performed by back transformation of c using the filter matrix Ŵ^-1

$$x\text{'}={\widehat{W}}^{-\text{1}}c$$ E3

3.2. A short introduction to cross trial phase statistics

Phase space methods are well established in non linear time series analysis; a good overview can be found in (Kantz & Schreiber, 2003). These methods are based on the famous Takens' Theorem stating that the dynamics of a system can be reconstructed from a single signal measured with a constant sampling rate (Takens, 1981). The concept based on these methods to concentrate on phase dynamics while neglecting amplitude dynamics has proven to be very powerful in analyzing noisy data of nonlinear or even chaotic systems (M. Rosenblum et al., 1996). As physiological systems very often show these features, phase based methods have been widely applied in this field. Especially the research on physiological oscillators and their mutual synchronization has a long tradition (Holst, 1937) and phase based methods for synchronization analysis have contributed a lot to this field of research (Glass, 2001, Mosekilde et al., 2002, Schiek et al., 1998, P. Tass, 1999), while the improvement of these methods still continues (Ocon et al., 2011, Wacker & Witte, 2011, Wagner et al., 2010).

To define instantaneous phases of a signal of interest one starts to transform the signal s using a rotational process by introducing an additional but related signal ŝ as a second coordinate in the so called embedding plane. The most common method to obtain these related signals is the Hilbert transform (Bracewell, 1999). Alternative, if one knows the characteristic period length τ of the process of interest the following definition is also an appropriate choice (Schiek et al., 1998)

$$\hat{s}\left(t\right)=s\left(t-\raisebox{1ex}{$\tau $}\!\left/ \!\raisebox{-1ex}{$4$}\right.\right)$$ E4

A precondition for both methods is to band pass filter the original signal, where the filter limits being defined by the frequency of interest. In the following we will use the discrete approximation of the Hilbert transform:

$$\hat{s}\left(t\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\pi $}\right.{\displaystyle \sum }_{k=1}^{\infty }\frac{s(t-k)-s(t+k)}{k}$$ E5

with t denoting the time step.

This leads to the following phase definition as the rotational angle around the origin within the two-dimensional embedding space:

$$\phi \left(t\right)=\arctan \frac{\hat{s}\left(t\right)}{s\left(t\right)}$$ E6

The idea of the cross trial phase statistics (CTPS) is to test for changes within the phase dynamics induced by the recurring events. Therefore the time series of the instantaneous phase is split into windows, where the center is defined by the event itself. In CTPS the statistics of these cross trial phases are analyzed. Stereotypic responses of a signal to stimuli for example, will reveal itself as a single cluster of phases short after the occurrence of the stimulus. Non stereotypic responses may also be detected by the method, as long as there are only a limited number of possible responses, and furthermore, this number is small compared to the number of stimuli. Two different responses for example would lead to two clusters within the cross trial phase plot. In systems without any stimuli related responses no phase clustering will be observed. In such a case, the phase values will form a more or less uniform distribution (Fig. 2D). Three examples for phase clustering (Fig. 2A-C) and one for a uniform phase distribution are shown in Fig. 2. The differences between the large or weak coupled signals can be detected by visual inspection very easily. With no difficulty one is able to differentiate between rather strict responses (Fig. 2A) or more weak responses (Fig. 2 B-C), where the phase clustering is more broad. Such a moderate phase clustering may be caused by measurement noise and/or by weak or irregular coupling.

In order to quantify these findings (i.e. to quantify the strength of the response to the stimuli) the so called Kuiper test is used. This statistical test is a refinement of the more common Kolmogorow-Smirnow test (KS test). The KS test is widely accepted to quantify the probability that two data sets are samples of the same distribution. As a preparatory step, the two data sets have to be transformed to their normalized cumulative distribution. A good introduction to both tests can be found in (Press et al., 2001). The application of these tests in phase resetting analysis is described in (Tass, 2004). To check for stimuli related phase responses we compare the cross trial phase distribution g(Ф) with the uniform distribution u(Ф).

Let s(t) be the original signal, Ф(t) be the instantaneous phase, normalized by 2π (i.e., the phase values range from 0 to 1), g(Ф) being the normalized cumulative cross trial phase distribution (g(Ф) ~ relative number of cross trial phases smaller than Ф, and g ranging from 0 to 1), and u(Ф) = Ф being the uniform density distribution. Using the absolute maximal difference between the two cumulative distribution

$$D=\underset{0\le \Phi \le 1}{\max }\mid g(\Phi )-u(\Phi )\mid$$ E7

the KS-test calculates the probability of error for rejecting the null hypothesis that g and u are samples of the same distribution with

$$P_{KS}(\lambda )=2{\displaystyle \sum }_{k=1}^{\infty }{(-1)}^{k-1}{e}^{-2{\text{k}}^2{\lambda }^2}$$ E8

with λ = V(√N + 0.155 + 0.24/√N) and N being the number of data points.

The disadvantage of the KS-test is that the accuracy decreases in case where the largest difference D is located close to the tails, i.e. Ф=0 or Ф=1. This however is considered in the Kuiper test by defining D in the following way

$$D=\underset{0\le \Phi \le 1}{\max }\left[g(\Phi )-u(\Phi )\right]+\underset{0\le \Phi \le 1}{\max }\left[u(\Phi )-g(\Phi )\right]$$ E9

and by using a modified probability statistic that reads

$$P_K(\lambda )=2{\displaystyle \sum }_{k=1}^{\infty }(4{\text{k}}^2{\lambda }^2-1)e^{-2{\text{k}}^2{\lambda }^2}$$ E10

Figure 2.

Cross trial phase statistics of strong (A), moderate (B-C) and non-coupled (D) signal responses.

In the case of a disproof of the null hypothesis (i.e. g(Ф) has a non-uniform distribution), P _K becomes 0 in the limit of D→1. Since P _K becomes very small for D→1, we introduced the negative logarithmic value

$${\text{p}}_{\text{K}}=-{\text{log}}_{\text{1}0}\left({\text{P}}_K\right)/N_T$$ E12

with N _T being the number of trials used. This is done to quantify the strength of the response to the stimuli independently from the number of trials used. It is a consequence of the applied statistical method that the test statistic -log₁₀(P _K ) increases the more events (i.e. stimuli) are taken into account for the analysis. In order to have comparable results the test statistic is normalized to the number of trials used, so that the p_K value is ranging between 0 and 1.

4. CTPS based cardiac and ocular artifact rejection

For automatic identification of cardiac and ocular artifact using CTPS data segments (trials) from the decomposed MEG signals needs to be defined around the event of interest (e.g., the R-peak of the QRS-complex or the peak of an eye blink). The window length should be large enough to cover the artifact signal and may be set to about 0.5s and 1s for the cardiac and eye blink artifact, respectively. Each window (trial) should be centered around the (main) peak of the artifact of interest The peak latencies for the two types of artifacts can easily be extracted from auxiliary ECG and EOG channels, while for R-peak detection an ECG channel may not be necessary, as shown by Escudero and colleagues (Escudero, et al., 2011).

For a reasonable definition of the instantaneous phases a bandpass filter should be applied to the ECG and EOG signals, which covers the main energy of the QRS power spectrum. The normal duration of the QRS-complex is typically in the range of about 0.06 – 0.10 s, while for normal subjects the QRS duration should be ≤0.12s (Kashani & Barold, 2005). The total eye blink duration (closing and opening of the eye lid) is even more variable and depends on the subject's condition and the task demands. Under normal condition the full blink duration d_blink is often reported to be in the range of 0.1s<d_blink<1s (Benedetto et al., 2011). Therefore, reasonable frequency ranges for the CTPS analysis to reveal field contributions from the QRS-complex and eye blinks in independent components are 10-20 Hz and 1-5 Hz, respectively. Slight (!) variations in the selection of the frequency range will not result in a complete failure of the method, but may reduce the significance level of the test statistic.

In Fig. 3 typical cardiac and eye blink artifacts are shown for one subject after signal decomposition using the Infomax algorithm (Bell & Sejnowski, 1995). In the first three rows components related to the cardiac artifact are shown, while the last component highly correlates (by means of Pearson's linear correlation coefficient c, with c=0.98) with the subject's eye blink activity. All components in red shown in Fig. 3 were automatically identified using different second and higher order statistics as described in (Dammers et al., 2008). Note, the third cardiac component (in grey) has much weaker amplitudes compared to the first two components and was identified by CTPS analysis only. Previous reports showed that the cardiac activity often splits into multiple components after ICA decomposition, where the signal strength of each IC can vary quite substantially (Dammers et al., 2008, Sander et al., 2002) However, it has been shown that CTPS is capable of identifying weak artifact components, where second and higher order statistics (including correlation analysis) failed (Dammers et al., 2008).

Figure 3.

Identification of ocular and cardiac artifacts.Typical examples of independent components (ICs) showing cardiac (a-c) and ocular (d) artifacts in one subject, during the first minute of an MEG experiment. ICs shown in red were identified using second or higher order statistics, or correlation analysis. In contrast to the amplitude based methods, CTPS analysis was able to identify all four components (a-d). Note, the poor correlation of IC#1 (a) and IC#15 (c) to the ECG signal, where correlation coefficient was 0.02 and 0.15, respectively. In contrast, both ICs were perfectly identified using CTPS analysis (significance level ≥ 0.5).

The number of components related to cardiac activity depends on the subject's individual electrical heart axis and the position of the MEG sensor array (e.g., seated or supine position). Therefore, during analysis one may observe that the latencies across the cardiac components are shifted with respect to the R-peak latencies of the ECG. In such cases, it is much more difficult to identify cardiac components by correlation analysis, even if the amplitude of the signal is large.

Such an example is shown in Fig. 3, where the first cardiac component (IC#1, Fig. 3a) has an absolute correlation coefficient of c=0.02 only, even though this component has the largest amplitude. Further analysis showed that the peak latencies of component IC#1 (Fig. 3a) and IC#15 (Fig. 3c) are correlated to the Q- and S-peak of the ECG signal, respectively (Fig. 4). Therefore, for both components the

The correlation of both signals can be increased by allowing a time dependent lag, but with the expense of being much more sensitive to non-cardiac or noise components.

correlation analysis failed, while the two ICs were identified with no problem as cardiac artifacts using CTPS. The reason here is that CTPS analysis is applied to a predefined time window (e.g., around the R-peak) to search for non-uniform phase distributions. As long as at one latency within the time window, the cross trial phase distribution differs significantly from a uniform distribution, CTPS will be able to identify such a component that is synchronized with the event of interest.

Figure 4.

Additional cardiac component identified by CTPS.Magnetic fields constructed (black) from IC#15 (cf. Fig. 3) are averaged to the onset of the R-peak. For comparison the ECG signal average is superimposed in red. Note the time shift between the R-peak onset (at time t=0s) and the peak latency of the magnetic signal reconstructed from component #15. The maximum field strength of this IC is even more than 180 ft at the MEG sensor level, while the maximum signal strength of IC#15 (7.8 a.u.) is about one order less compared to IC#1 (78.6 a.u.). This artifact component was solely identified by CTPS, while the applied amplitude based methods failed.

An estimation of how well the CTPS based (or any other) artifact rejection method performs, the artifact signal removal can easily be visualized (Fig. 5). For the example described above, the MEG signals are aligned to the onset of the R-peak and the peak latency of the eye blink. Cross trial averages of the MEG signal will enhance the content of the artifact before and after artifact rejection (Fig. 5). For quantification of the artifact removal one can use the rejection performance measure as introduced in (Dammers et al., 2008).

Figure 5.

Artifact rejection performance.MEG signals averaged with respect to the R-peak of the ECG (left) and the eye blink signal (right) before (top row) and after (bottom row) CTPS based artifact removal. The respective reference signal is superimposed in red and blue (in arbitrary scaling) for the cardiac and ocular artifact, respectively. The RMS based rejection performance (bottom row) indicates the quality of the artifact cleaning process, using the performance measure described in (Dammers et al., 2008). The MEG data were recorded using the MAGNES WH3600 from 4D-neuroimaging.

5. CTPS based identification of strong and weak brain responses

In most cases unaveraged neuromagnetic brain responses are too weak (or too noisy) to be reliably identified by its amplitude statistic only. Averaging on the other hand will enhance the most prominent brain activity by suppressing uncorrelated noise, with the expense of not preserving the temporal dynamics of each individual brain response.

In the previous sections we showed that the analysis of the cross trial phase dynamics in decomposed MEG signals discloses artifact components that are highly correlated to the phase dynamics from trials of a reference signal. Such a reference signal may be the QRS-complex or the eye blink signals recorded by the ECG and EOG channel, respectively. In principal, CTPS analysis can be applied to any other event of interest in order to investigate the phase dynamics of signals and its dependency to predefined points in time. For example, burst onsets in the electromyogram (EMG) or stimulus onsets in the trigger channel may be used to define trials of interest for CTPS analysis. In this way, even weak but stimulus dependent brain responses can be disclosed by its phase clustering around the event of interest, while at the same time, the temporal dynamics of the underlying sources will be remained. The CTPS based identification of stereotypic brain responses in MEG data is best performed if the effect of signal mixtures is removed. Therefore, the decomposition of the MEG data using for example independent component analysis (ICA) is a mandatory requisite.

5.1. CTPS analysis applied to auditory cued MEG signals

For evaluation the method is applied to decomposed MEG data (MAGNES 2500WH, 4D-Neuroimaging) from a simple auditory experiment. About 300 click tones (1 kHz, 50 ms duration, SOA of 1.5 s ± 0.5 s) were presented to both ears of a male volunteer. ICA (Infomax) was applied separate noise and artifacts from brain activity. The trigger onset latencies were used to define 300 trials using a time window ranging from stimulus onset to 150 ms post stimulus. From the MEG raw data (frequency range 1-200 Hz) only, one can be rather vague about stimulus dependent activity (Fig. 6a). After source separation (i.e., by means of ICA) a stimulus related component is identified by CTPS analysis (Fig. 6a). Since analysis of phase dynamics requires a frequency band of interest, either prior knowledge of the frequency components or multiple frequency ranges should be used for CTPS analysis. Therefore, within the frequency ranges of 2‒4, 4‒8, 8‒12, 12‒16, and 16‒20 Hz we scanned for the largest p_K value in all ICs. The component with the highest significance level of the

Figure 6.

CTPS based source identification.In a) unaveraged MEG raw signals (gray) are shown for the first 15 s of an auditory experiment. All 148 MEG signals are superimposed, while the blue lines indicate the stimulus onset times. The independent component (IC#20 in red) with the maximum significance level after CTPS analysis is plotted below the MEG signal. The peaked activity of this component reflects brain activity that is highly synchronized with the stimulus. In b) the MEG raw signals (gray) are shown together with MEG signals after back-transformation of the identified component (red). The amplitude varies, but stimulus dependent brain activity is now obvious.

test statistic (Eq. 11) was found in IC#20 at a frequency range of 4‒8 Hz (p_K,max=0.45), while for the frequency ranges 8‒12 and 12‒16 Hz the maximum p_K value was still about 0.44.

Typically, p_K values of ICs related to brain responses are smaller compared to ICs which are related to cardiac activity (where p_K values typically are above 0.5). The reason is that we do not precisely know the latencies of the trial-by-trial brain activity. For example, for auditory cued brain responses we expect strong evoked activity around 100 ms after stimulus onset. The latencies of the true brain responses however jitter a few milliseconds from trial to trial, with the effect, that the phase clustering around the average response latency is not as narrow as it is e.g. for cardiac activity. Nevertheless, the maximum p_k value of IC#20 was found to be p_K,max=0.452, while the mean value across all ICs was p_K,mean=0.016 ±0.038.

Back-transformation of the identified signal component to the MEG sensor space revealed stimulus dependent evoked activity with a much better signal-to-noise ratio (Fig. 6b). Subsequent source analysis of the identified component revealed that the corresponding signal was localized in the primary auditory cortices. In Fig. 7 five trials (trial #6-10) are shown before and after CTPS based source extraction. Compared to trials of the MEG raw signal the signal-to-noise ratio after back-transformation of the identified source component (IC#20) is much better. With such an improvement in SNR CTPS analysis enables source localization on the single-trial level.

Figure 7.

Signal enhancement in MEG single-trials.Selected single-trials from an auditory MEG experiment. In a) trials of unprocessed MEG raw signals (gray) are shown. In b) single-trials of MEG signals after back-transformation (in red) of one source component (IC#20 in Fig. 6) as identified by CTPS analysis. The blue lines indicate the stimulus onset times of the auditory click tones.

5.2. CTPS analysis on MEG signals evoked by voluntary finger movements

In the following experiment CTPS analysis is applied to search for stereotypic brain responses with respect to voluntary finger movements. The MEG data of one subject were recorded using a 248 magnetometer whole head MEG system from 4D-Neuroimaging (MAGNES 3600WH). The right handed subject voluntarily lifted up and down the right index finger 117 times. Between two subsequent finger lifts the subject was asked to perform a short rest of about 1‒2 s.

After ICA decomposition cross trial phase distributions were calculated for each IC within a time window of 600 ms centered around the movement onset detected by a photoelectric barrier. We have analyzed the phase dynamics of all ICs within different frequency bands (cf. section 5.1). Results from CTPS analysis revealed that the phase dynamics for different frequency bands qualitatively resemble, therefore the results presented here are limited to the frequency band corresponding to the largest p_K value (Fig. 8).

Figure 8.

Significance level of components related to voluntary finger movement.CTPS analysis applied to a voluntary finger movement experiment. ICs with p_k values greater or equal the 98^th percentile of all p_k values (indicated by the black horizontal line) are marked in red. ICs that have been identified as artifacts are shown in blue.

For subsequent single-trial data analysis, only ICs related to significant brain responses (as defined by large p_k values in Fig. 8) were back-transformed into the MEG signal space. The separation of these components from other contributions results in a much better SNR in the single-trial data as illustrated in Fig. 9. On average the strongest peak in the MEG signal is found at the time t=0 ms, when the finger is just lifted (Fig. 9a), while the peak latencies of the individual brain responses may vary as indicated by the red arrow in Fig. 9.

For the same experiment, current source localization was then performed to single-trial MEG data, as well as to signal averages aligned to the onsets of voluntary finger movements. The difficulty here is that two nearby sources, i.e. the primary motor (MI) and the somatosensory cortex (SI), are activated at the same time, where a separation of the two focal sources is challenging for any source analysis. Here we used magnetic field tomography (MFT, Ioannides et al., 1990) to reveal source localizations of the underlying neuromagnetic activity from both averaged and single-trial data (Fig. 10). Source reconstruction based on standard averaging revealed distributed activity over both motor and somatosensory areas (Fig. 10a) for the time indicated by the red arrow as depicted in Fig. 9.

In Fig. 10 b-c localization results of a single epoch (trial #5, cf. Fig. 9b) are illustrated, while in Fig. 10c the reconstruction was applied to MEG signals derived from the IC with the maximum p_K value (IC#12, cf. Fig. 8). Due to the poor signal-to-noise ratio for this single-trial the reconstructed activity shown in Fig. 10b is widely distributed over the brain, while its maximum was found to be close to the cerebellum (not shown). In contrast, the CTPS based source localization of the components with the strongest p_K value (IC#12) reveals a more focal source pattern for this single-trial, while its maximum activity (in red) is located in MI (Fig. 10c). For comparison, localization results based on signal averages are shown in (Fig. 10a). Although uncorrelated noise is fairly suppressed after signal averaging, the maximum activity here is distributed over both areas, the motor and somatosensory cortex, since both areas are activated at the same time.

Figure 9.

MEG signal averages vs. single-trial analysis.In a) MEG signal averages aligned to the onset of voluntarily finger movements (indicated by the blue line). In b). unaveraged MEG signals are shown for one selected trial (trial #5). The signal-to-noise ratio does not allow to make any statements about the stimulus related activity. For the same trial the MEG signal derived from the IC with the maximum p_K value is shown in c) (cf. Fig. 8). Note, the shift in peak latency (indicated by the red arrow) of the individual response in trial#5 (b-c) compared to the averaged peak latency in a). The improvement in the SNR after CTPS analysis as shown in c) is evident compared to b).

Localization of the 2^nd to 5^th strongest components (with respect to their p_K value, cf. Fig. 8) revealed “independent“ somatosensory activity (by means of statistical independence as defined by ICA), where the origin of the maximum activity was found to be within the somatosensory cortex, as indicated by the location of the central sulcus (Fig. 11). By adding the current densities of each localized component, the resulting activation pattern closely resemble the activation one gets, when localizing the fields constructed from all four components. However, the major advantages in CTPS based source analysis is i) the method is capable of identifying components that are highly synchronized to an event of interest and ii) nearby sources which are activated at the same time can be separated, where standard averaging based methods may fail.

Figure 10.

Source reconstruction of signal averages vs. single-trial data.All localization results are shown for the latency indicated by the red arrow in Fig. 9. Standard averaging source analysis (a) shows that the activity is distributed in both motor and somatosensory areas. The white arrow indicates the central sulcus. The reconstructed current density distribution of the unaveraged data in b) is fairly widespread over the entire brain, while the strongest activity was found to be in the cerebellum (not shown). In c) the CTPS based source localization of the signal shown in Fig. 9c is much more focal and its maximum activity (in red) is located in the primary motor cortex.

Figure 11.

CTPS based source reconstruction.Source localization of magnetic fields constructed from ICs with p_K values ranging from 0.12 ≤ p_K< p_Kmax as identified by CTPS analysis (cf. Fig. 8). The localization of the corresponding single-trial MEG signal shows that the activity can be attributed to a weak stimulus related somatosensory response, which could not be isolated using averaging based localization (cf. Fig. 10). The white arrow indicates the central sulcus.

6. Conclusion

Cross trial phase statistics (CTPS) applied to decomposed MEG signals allows for both, user independent identification of artifacts, as well as the extraction of strong and weak event related brain responses. With respect to cardiac artifact (CA) rejection, the ICA decomposition of the recorded MEG signal often shows multiple independent components (ICs) that can be attributed to CA. Since, the amplitudes of CA components can vary quite substantially; it is therefore difficult for any automatic procedure that exclusively relies on amplitude information alone. In contrast, CTPS based identification of CA related components is very sensitive. Even in cases of weak signal strength, the method identifies multiple components that can be attributed to cardiac activity, where amplitude based methods fail (Figs. 3-4). For ocular artifacts, CTPS as well as second and higher order statistics are able to identify eye blink components due to its larger amplitudes (Fig. 3). In any case, the performance of the artifact rejection can be visualized (Fig. 5) and estimated using the rejection performance metric described in (Dammers et al., 2008).

Within CTPS analysis the information about the occurrence of an event is of importance. Phase clustering of an event related signal is enhanced if trials of the phase dynamics are perfectly aligned to the time of the event. In the case of cardiac artifact rejection the CTPS based identification is very accurate, since all trials can precisely be aligned by extracting the R-peak latencies of the QRS-complex signal from the ECG signal.

Since CTPS is a highly sensitive approach for extracting weak event-related fields, the same strategy can also be used for the identification of stereotypic brain responses. For this, time information about the occurrence of each brain response must be available. Therefore, the latencies of the brain responses which do correlate to an external event are estimated by make use an external reference signal, such as the trigger or response channel. These channels provide time information about a stimulus presentation or a subject's response to a stimulus (e.g. a button press). However, the delay between for example a stimulus and the event related brain activity is not strictly fix. This results in a drop of p_K values (p_K describes the significance of a component being synchronized to an event of interest) when compared to values of CTPS based cardiac artifact rejection. Nevertheless, as long as the brain responses jitter a few milliseconds only, CTPS is able to identify event related brain activity in unaveraged decomposed MEG data (Figs. 6,10,11).

In this way the method enables single-trial data analysis by providing reconstructed MEG signals with sufficient SNR (cf. Fig. 7) for single-trial source reconstruction. The CTPS based single-trial source localization provides equal or even improved spatial accuracy compared to standard averaging techniques (Fig. 10). Reconstruction of identified single (or a group of) components have the advantage of being able to separate nearby sources although they are activated at the same time (Figs. 10c, 11). Therefore, if single-trial brain dynamics or the source separation of nearby sources is of interest, CTPS based source reconstruction seems to be a promising tool for the identification of event related brain activity in unaveraged decomposed MEG data.